Understanding Two-Layer Neural Networks with Smooth Activation Functions

Changcun Huang
cchuang@mail.ustc.edu.cn
Shuitu Institute of Applied Mathematics, Chongqing 400700, P.R.C

This paper aims to understand the training solution, which is obtained by the back-propagation algorithm, of two-layer neural networks whose hidden layer is composed of the units with smooth activation functions, including the usual sigmoid type most commonly used before the advent of ReLUs. The mechanism contains four main principles: construction of Taylor series expansions, strict partial order of knots, smooth-spline implementation and smooth-continuity restriction. The universal approximation for arbitrary input dimensionality is proved and the explanation of training solutions is given. Through the principles proposed, the mystery of “black box” of the solution space is largely revealed. The new proofs employed also enrich approximation theory.

Keywords: two-layer network, training solution, universal approximation, smooth activation function, smooth spline.

1 Introduction↩︎

Smooth activation functions (e.g., logistical sigmoid and tanh cases) were the most common choice in neural networks [1]–[3] before the popularity of rectified linear units (ReLUs) [4] since 2010. In engineering, the mechanism of two-layer neural networks with smooth activation functions is elusive and referred to as a “black box”. We study this problem due to the following reasons. First, it is of the simplest network architecture leading to a deeper one and the associated principle may be the foundation of the latter. Second, the activity of neurons [5] could be complex and the logistic sigmoid function \(\sigma(x)=1/(1+e^{-x})\) had been employed to model their behavior (e.g., [6]); so the principle is possibly useful in understanding neuroscience.

1.1 Problem Description↩︎

Given a set \(D\) of data points \((\boldsymbol{x}_i, y_i)\) with \(\boldsymbol{x}_i \in \mathbb{R}^n\) and \(y_i \in \mathbb{R}\), a two-layer neural network \(\mathfrak{N}\) can interpolate them by an output function \[g(\boldsymbol{x})=\sum_j\lambda_j\sigma(\boldsymbol{w}_j^T\boldsymbol{x}+b_j)\] with \(g(\boldsymbol{x}_i) = y_i\), whose parameters are derived from the back-propagation algorithm [7], where \(\sigma(x)\) is a smooth activation function of hidden-layer units.

The mechanism of \(g(\boldsymbol{x})\) interpolating \(D\) is unknown and is part of the “black box” of neural networks. We want to solve this problem by constructing the solution of universal approximation. The interpolation process can be interpreted as realizing a function that passes through the given data; thus if any function can be constructed in the way associated with training solutions, the black box could be uncovered.

In theoretical physics, the capability of explaining or predicting natural phenomena is a criterion for judging whether a proposed theory is successful. We also use this criterion to test the effectiveness of our theory. More specifically, first, given a training solution, we will explain how it is formulated and what the meaning of the obtained parameters is; second, we will check whether the solution predicted by the theory would occur in experiments. Especially, in the realm of artificial neural networks, the first rule above means that we can manually construct a solution in a deterministic way, which is originally obtained by gradient-descent method; and this is one of the ingredients of black box of neural networks.

Hereafter, we will call \(g(\boldsymbol{x})\) of equation 1.1 derived from the back-propagation algorithm a training solution of two-layer neural networks.

1.2 Main Background↩︎

The universal approximation of two-layer neural networks with “sigmoidal unit” [8], whose activation function is, for example, \(\sigma(x)=1/(1+e^{-x})\), had been proved for more than thirty years in several ways, such as [9], [10] and [11]. The case of more general smooth activation functions was also intensively studied (e.g., [12]–[20]). However, those works are mainly of pure mathematical research and have no connection with engineering applications. Thus, so far little was known about the mechanism of training solutions. We try to bridge this gap in this paper.

1.3 Organization of the Paper↩︎

The paper is mainly classified into three parts. The first that includes sections from 2 to 4 is for univariate-function approximation and the second (sections 5 and 6) is for the multivariate case. The third (section 7) is about experimental verification or explanation. The detailed description is as follows.

Section 2 constructs local approximation through Taylor series expansions and the result is applicable to any smooth activation function. The construction method will be generalized to the multivariate case in later section 5. The result of this section has its own meaning and is also the foundation of global solutions of subsequent sections.

Section 3 turns to global approximation, in the sense that a two-layer neural network \(\mathfrak{N}\) with generalized sigmoidal units (definition 3) can approximate a given function \(f(x)\) through a smooth spline \(s(x)\), not necessarily within a sufficiently short interval as the method of Taylor series expansions. The spline \(s(x)\) is constructed by integrating a piecewise linear approximation to the \(k\)th derivative of \(f(x)\) (theorem 4); \(s(x)\) is then implemented by \(\mathfrak{N}\) (theorem 6), with the approximation error controlled by theorem 5; and this leads to universal approximation for univariate functions (theorem 7). Corollary 1 investigates the tanh-unit case. We give another solution in proposition 2 that could be observed in experiments.

Section 4 reduces the spline implementation to the solution of a system of linear equations, by which a theoretical framework containing all the possible solutions is presented, including the special ones of section 3.

Section 5 constructs the Taylor series expansion of a multivariate function \(f(\boldsymbol{x})\) for \(\boldsymbol{x} \in \mathbb{R}^n\) and \(n \ge 2\) via neural networks and is the generalization of section 2. Section 6 deals with the global approximation to \(f(\boldsymbol{x})\) as section 3. A multivariate spline (theorem 11) approximating \(f(\boldsymbol{x})\) is obtained by the integral of a piecewise linear approximation to a directional-derivative hypersurface (definition 14) of \(f(\boldsymbol{x})\). A basic property of two smoothly connected polynomial pieces is given in theorem 12. The principle of smooth-continuity restriction is proposed in theorem 15. Theorem 17 constructs a desired spline on a certain type of partition through neural networks. Theorem 18 proves the universal approximation for arbitrary higher-dimensional input. Theorem 19 and proposition 8 provide some rules for two-sided solutions ubiquitous in experiments.

Section 7 verifies the theory by some typical experimental examples or uses the preceding results to explain experimental phenomena. Section 8 summarizes the main principles of the black box of two-layer neural networks. Section 9 ends this paper by a discussion.

2 Univariate Local Approximation↩︎

Throughout this paper, the “smooth function” \(f: [0, 1] \to \mathbb{R}\) means that the derivatives of \(f(x)\) up to some order required are continuous, and if its \(k\)th-order derivative satisfies the continuous property, we sometimes say that it is smooth with order \(k\) or \(f(x) \in C^k([0, 1])\), where \(C^k([0, 1])\) denotes the set of functions on \([0, 1]\) whose \(k\)th derivative is continuous; \(C^k(\mathbb{R})\) represents the set of functions defined on \((-\infty, +\infty)\) with smoothness order \(k\). The multivariate case is similar and the rigorous definition will be given in section 5.

The term “local approximation” is in the sense that the approximation is based on Taylor series expansion whose approximation error can only be ensured within a sufficiently small neighborhood of a point, and so is the multivariate case of section 5. Note that the results of this section are applicable to any smooth activation function, including the generalized sigmoidal type to be discussed later.

2.1 Polynomial Construction↩︎

The main idea is to use weighed activation functions of the hidden-layer units of a two-layer neural network to construct a desired Taylor-series expansion approximating a given function; theorem 1 below is a foundation.

In this paper, the approximation error is measured by the \(L^2\)-norm distance, that is, \(\lVert f(\boldsymbol{x})\rVert_2=(\int_{U} |f(\boldsymbol{x})|^2 d\boldsymbol{x})^{1/2}\) for \(\boldsymbol{x}\in \mathbb{R}^n\), including the special case \(n=1\) of this section.

To satisfy equation 2.6, we obtain a system of linear equations, with each being \[c_{1\nu}\lambda_1 + c_{2\nu}\lambda_2 + \dots + c_{k\nu}\lambda_k = a_{\nu}\] for \(\nu = 0, 1, \dots, m\), which can be expressed as \[A\boldsymbol{\lambda} = \boldsymbol{a},\] where \(A\) is from equation 2.4 and \(\boldsymbol{a} = [a_0, a_1, \dots, a_m]^T\). So \[\begin{align} \lVert T(x) - \sum_{i}&\lambda_i\phi_i(x) \rVert_2 = \lVert \sum_{i}\lambda_it_i(x) - \sum_{i}\lambda_i\phi_i(x) \rVert_2 \\ &= \lVert \sum_i\lambda_ir_i(x)\rVert_2\le \sum_i|\lambda_i|\varepsilon_i, \end{align}\] which is equation 2.7.

2.2 Generalized Wronskian Matrix↩︎

Wronskian determinant was studied in the theory of total positivity [21] and the associated term “Wronskian matrix” in [22] is a square matrix. We borrow that terminology with slight modifications.

By theorem 1, we have \[\begin{align} \lVert f(x) - \sum_{i}&\lambda_i\phi_i(x) \rVert_2 = \lVert f(x) - T(x) + T(x) - \sum_{i }\lambda_i\phi_i(x) \rVert_2 \\ &\le \lVert f(x) - T(x)\rVert_2 + \lVert T(x) - \sum_{i}\lambda_i\phi_i(x) \rVert_2 \\ &= \lVert R(x)\rVert_2 + \epsilon(m) = \varepsilon. \end{align}\] Since \(\phi_i(x)\)’s are smooth, \(\phi_i^{(m+1)}(\xi_i)\) of equation 2.11 is bounded. If \(m\) and \(\lambda_i\)’s are fixed, when \(|x-x_0|\) is small enough, the weighted error \(\epsilon(k)\) could approach zero with arbitrary precision, because \(\varepsilon_i\)’s of equation 2.7 tend to be zero; and so is \(\lVert R(x)\rVert_2\).

2.3 Construction of Local Approximation↩︎

We know that \[\psi^{(k)}(wx) = w^k\frac{d^k\psi(y)}{dx^k} = w^k\psi^{(k)}(y),\] where \(y = wx\). When \(w\) is sufficiently small, we have \(y \in [0, 1]\). Then \(\psi^{(k)}(y)\) is bounded, that is, \(|\psi^{(k)}(y)| \le M\), owing to the smoothness property of \(\psi(x)\); the conclusion is obvious.

The main idea is to control the degree of the Taylor series expansion of each \(\phi_i(x)\) by scaling the parameter \(w_i\) at different levels, as adjusting the scaling parameter of wavelets. To each \(\phi_i(x) = \sigma(w_ix + b_i)\), let \(y = w_ix + b_i\) and \(y_i = w_ix_0 + b_i\). By equation 2.12, the corresponding generalized Wronskian matrix is \[\mathcal{W}(\sigma(y_i)) = \begin{pmat}({}) \sigma(y_1) & \sigma'(y_1)w_1 & \cdots & \sigma^{(m)}(y_1)w_1^{m} \cr \sigma(y_2) & \sigma'(y_2)w_2 & \cdots & \sigma^{(m)}(y_2)w_2^{m} \cr \vdots & \vdots & \ddots & \vdots \cr \sigma(y_{m+1}) & \sigma'(y_{m+1})w_{m+1} & \cdots & \sigma^{(m)}(y_{m+1})w_{m+1}^{m} \cr \end{pmat}.\] Let \(w_1 = \Delta t^{(1+c)}\) for \(c > 0\), \(w_j = \Delta t^{1/(j-1)}\) for \(j = 2, 3, \dots, m\), where \(0 < \Delta t < 1\), and \(w_{m+1} =\Delta t^{-(m-1)}\); then equation 2.24 becomes \[\begin{align} &\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\mathcal{W}(\sigma(y_i), \Delta t) = \\ &\begin{pmat}({}) \sigma(y_1) & \sigma'(y_1)\Delta t^{1+c} & \sigma''(y_1)\Delta t^{2(1+c)} & \cdots & \sigma^{(m)}(y_1)\Delta t^{m(1+c)} \cr \sigma(y_2) & \sigma'(y_2)\Delta t^{1/1} & \sigma''(y_2)\Delta t^{2/1} & \cdots & \sigma^{(m)}(y_2)\Delta t^{m/1} \cr \sigma(y_3) & \sigma'(y_3)\Delta t^{1/2} & \sigma''(y_3)\Delta t^{2/2} & \cdots & \sigma^{(m)}(y_3)\Delta t^{m/2} \cr \vdots & \vdots & \vdots & \ddots & \vdots \cr \sigma(y_{m+1}) & \sigma'(y_{m+1})\Delta t^{1/m-1} & \sigma''(y_{m+1})\Delta t^{2(1/m-1)} & \cdots & \sigma^{(m)}(y_{k})\Delta t^{-(m-1)} \cr \end{pmat}, \end{align}\] which can be expressed as \[\mathcal{W}(\sigma(y_i), \Delta t) = \begin{pmat}({}) \sigma(y_1) & o(\Delta t) & \cdots & o(\Delta t) \cr \sigma(y_2) & \sigma'(y_2)\Delta t & \cdots & o(\Delta t) \cr \vdots & \vdots & \ddots & \vdots \cr \sigma(y_{m+1}) & \sigma'(y_{m+1})\Delta t^{1/m-1} & \cdots & \sigma^{(m)}(y_{m+1})\Delta t^{-(m-1)} \cr \end{pmat},\] where “\(o\)” represents the little-oh notation. Among the entries of \(\mathcal{W}(\sigma(y_i), \Delta t)\), neglecting those infinitesimal ones relative to \(\Delta t\), we have

\[\mathcal{W}'(\Delta t) = \lim_{\Delta t \to 0}\mathcal{W}(\sigma(y_i), \Delta t) = \begin{pmat}({}) \sigma(y_1) & 0 & \cdots & 0 \cr \sigma(y_2) & \sigma'(y_2)\Delta t & \cdots & 0 \cr \vdots & \vdots & \ddots & \vdots \cr \sigma(y_{m+1}) & \sigma'(y_{m+1})\Delta t & \cdots & \sigma^{(m)}(y_{m+1})\Delta t^{-(m-1)} \cr \end{pmat},\] whose determinant \[\det \mathcal{W}'(\Delta t)) = \prod_{i=0}^{m}\sigma^{(i)}(y_{i+1}) \ne 0\] when each bias \(b_i\) in \(y_i = w_ix_0 + b_i\) is properly set to make \(\sigma^{(i)}(y_{i+1}) \ne 0\). Thus, \(\mathcal{W}'(\Delta t)\) is nonsingular.

According to theorems 1 and 2, the key point is to construct a nonsingular Wronskian matrix, for which lemma 2 provided a solution.

2.4 Related Work↩︎

[13] and [15] constructed \(x^k\) through smooth activation functions; [23] also realized Taylor series expansions via a two-layer neural network. Our difference lies in three aspects: (a) The framework of theorem 2 is general and can include their solutions as special cases; (b) Our construction method is novel and can be directly generalized to the multivariate case without requiring additional ridge functions.

3 Univariate Global Approximation↩︎

The approximation via Taylor series expansion is local and the approximation error can only be assured within a sufficiently small neighborhood of a point, for which a new mechanism related to splines is introduced and the associated solution is called a global one. Both theoretical and experimental results will demonstrate the necessity of this introduction.

Section 3.1 provides several definitions, notations and some knowledge of splines. Section 3.2 proposes a spline-construction method for the convenience of network implementation. Section 3.3 paves the way for approximation-error control. Section 3.4 defines local and global units for different solution patterns. Section 3.5 constructs a concrete solution (one-sided case) for universal approximation. Section 3.6 gives another solution (two-sided case) to explain experiments.

3.1 Preliminaries↩︎

We develop the theory on the basis of splines and first introduce some related concepts, notations and elementary results. Write \[\begin{cases} \;x_{+} := \max(0, x) \\ \;\;\;\;x_{+}^m := (x_+)^m \end{cases}\] [24]. The space of the order-\(m\) splines discussed in this paper is \[\begin{align} \mathfrak{S}_1^m(\Delta) = \{s: s(x)=s_i(x) \in \mathcal{P}_m \;&\text{for} \;x \in I_i, D^{\mu}s_{i}(x_{i})=D^{\mu}s_{i+1}(x_{i}) \;\text{for} \;i \ne \zeta, \\ \;&i = 1, 2, \dots, \zeta\}, \end{align}\] where \(\mathcal{P}_m\) is the set of polynomials of degree \(m\), \(I_1 = [x_0, x_1]\) and \(I_j = (x_{j-1}, x_{j}]\) for \(j = 2, 3, \dots, \zeta\) are the intervals formed by \(0 = x_0 < x_1 < x_2 < \cdots < x_{\zeta-1} < x_{\zeta} = 1\), \(D^{\mu}\) is the \(\mu\)th derivative operator for \(\mu = 0, 1, \dots, m-1\), and \[\Delta = \{x_{\nu}: \nu = 1, 2, \dots, \zeta-1\}\] is the set of the knots. Throughout this paper, we call \(I_1\) and \(I_j\)’s the subintervals of \(\Delta\) or we say that they are derived from \(\Delta\). The notation of equation 3.2 is of [22]’s modified version and \(\mathfrak{S}_1^m(\Delta)\) is equivalent to [25]’s \(\mathscr{S}(m, x_0, x_1, \dots, x_{\zeta})\).

Note that each \(s(x) \in \mathfrak{S}_m(\Delta)\) belongs to \(C^{m-1}([0, 1])\). By the property of smooth splines, two adjacent polynomials \(s_{\nu}(x)\) and \(s_{\nu-1}(x)\) of \(s(x)\) for \(\nu = 2, 3, \dots, \zeta\) satisfy the recurrence formula \[s_{\nu}(x) = s_{\nu-1}(x) + \alpha_{\nu-1}(x - x_{\nu-1})^m_+\] for \(x\in I_{\nu}\), which ensures \(D^{\mu}s_{\nu}(x_{\nu-1})=D^{\mu}s_{\nu-1}(x_{\nu-1})\) of equation 3.2 and is the foundation of spline construction in this paper.

3.2 Splines via Local Linear Approximation↩︎

The approximation capability of a spline \(s(x) \in \mathfrak{S}_1^m(\Delta)\) to continuous functions had been proved (e.g., [25]). In order to realize a spline through neural networks, we give a new construction method.

We first give an example. To any point \(x_0 \in (0, 1)\), a linear function \(l(x) = ax + b\) can be found to approximate \(f'(x)\) within a neighborhood \(\delta(x_0)\) of \(x_0\), where parameters \(a\) and \(b\) can be obtained by the first-order Taylor series expansion of \(f(x)\) at \(x_0\). Then we have \[f'(x) = ax + b + o(x-x_0).\] The original function can be represented as \[f(x) = F(x) := \int_{x_0}^xf'(x)dx + f(x_0),\] since \(F(x_0) = f(x_0)\) and \(F'(x) = f'(x)\); so \[f(x) = \frac{ax^2}{2} + bx + c + \int_{x_0}^xr(x)dx,\] where \(c = -a^2x_0^2/2 - bx_0 + f(x_0)\) and \(r(x) = o(x-x_0)\). Because \(r(x)\) is continuous and bounded in \(\delta(x_0)\), it has a maximum \(r(x_{\tau})\) at some point \(x_{\tau} \in \delta(x_0)\), which is also \(o(x-x_0)\). Thus, \[\int_{x_0}^xr(x)dx < r(x_{\tau})(x-x_0) = o\big((x-x_0)^2\big).\] That is, \[f(x) = \frac{ax^2}{2} + bx + c + o\big((x-x_0)^2\big),\]

Divide \([0, 1]\) into \(\zeta\) intervals \(I_i\)’s. If the length of each \(I_i\) is sufficiently small, the linear approximation above can be done for all of \(I_i\)’s, yielding \[\begin{align} f(x) &= \frac{a_ix^2}{2} + b_ix + c_i + r_i(x) \\ &= \frac{a_ix^2}{2} + b_ix + c_i + o\big((x-x_{i0})^2\big) \end{align}\] for \(x, x_{i0} \in I_i\). Let \(p_{i}(x) = a_ix^2/2 + b_ix + c_i\). When \(i = 1\), the parameter \(c_1\) of \(p_1(x)\) is obtained by fulfilling a selected point of \(f(x)\) on \([0, x_1]\). To ensure the continuity between \(p_{j}(x)\) and \(p_{j-1}(x)\) for \(j = 2, 3, \dots, \zeta\), use \(p_{j-1}(x_{j-1}) = p_{j}(x_{j-1})\) to determine \(c_j\). Then all of \(p_{i}(x)\) ’s comprise a smooth spline \(\mathcal{S}(x) \in \mathfrak{S}_1^2(\Delta)\), which can approximate \(f(x)\) with arbitrary precision, provided that the length of each \(I_i\) is sufficiently small; the approximation error is \[\begin{align} E = \lVert \mathcal{S}(x) &- f(x) \rVert_2 = \big(\sum_i\int_{x_{i-1}}^{x_i} (p_i(x) - f(x))^2dx\big)^{1/2} \\ & < \big(\sum_i\int_{x_{i-1}}^{x_i} r_i(x)^2dx\big)^{1/2} \\ & \le \big(\sum_i\int_{x_{i-1}}^{x_i} \varepsilon^2dx\big)^{1/2} = \varepsilon, \end{align}\] where \[\varepsilon = \max_i\{\beta_i = \sup_{x \in I_i}|r_i(x)|: 1\le i\le \zeta\}\] can be arbitrarily small when \(l = \max\{|x_{i+1} - x_{i}|: 0 \le i \le \zeta-1\}\) is small enough.

The above procedure can be easily generalized to higher-degree polynomials. We use a linear function to approximate the \(m-1\)th derivative of \(f(x)\) at each \(I_i\), that is, \[f^{(m-1)}(x) = a_ix + b_i + o(x-x_{i0}),\] for \(x, x_{i0} \in I_i\), thereby \[f(x) = p_{i}(x) + o\big((x-x_{i0})^{m}\big),\] where \(p_{i}(x)\) is a polynomial of degree \(m\) whose coefficients are derived from the repeated integral of \(a_ix + b_i\) of equation 3.13. Polynomials \(p_{i}(x)\)’s for all \(i\) comprise a spline \(s(x) \in \mathfrak{S}_1^m(\Delta)\). The approximation error of \(s(x)\) to \(f(x)\) can be analysed analogously to equation 3.11.

3.3 Zero-Part Error of Units↩︎

To activation function \(\sigma(x)\) of a generalized sigmoidal unit \(\mathcal{U}\), we call the interval \((-\infty, x_0)\) with \(x_0 < 0\) the zero part of \(\mathcal{U}\), if the contribution of \(\mathcal{U}\) on \((-\infty, x_0)\) to function approximation can be ignored due to \(\lim_{x\to-\infty}\sigma(x)=0\). The case of activation function \(\phi(x) = \sigma(wx + b)\) can be similarly defined. The zero part of units should be exclusively designed to control approximation error.

Note that if a function \(f(x) \in C^m([0,1])\), its \(m-1\)th derivative \(f^{(m-1)}(x)\) is smooth, and hence \(f^{(m-1)}(x)\) can be approximated by a piecewise linear function that could lead to a smooth spline of order \(m\) approximating \(f(x)\). The multivariate case is similar. It’s the reason that the functions to be discussed in the remaining part of this paper are often confined to the condition as \(f(x) \in C^m([0,1])\).

Figure 1a gives an example intuitively demonstrating the main procedures of the proof. We first modify \(\phi(x)\) to be \[\phi_t(x) = \sigma(w''x + b'') = \sigma(\rho(wx + b) + \gamma')\] where \(\rho > 1\), satisfying \[\phi_t(x_k) = \phi(x_k) = \varepsilon\] and \[\phi_t(x) < \phi(x) < \varepsilon\] for \(x < x_k\), for which two steps are required. First, scale the variable as \(\phi_s(x) = \sigma(\rho(wx + b))\). As shown in Figure 1a, the black cure is the scaled function \(\phi_s(x)\) and the blue curve is the original \(\phi(x)\). Because \(\sigma(x)\) monotonically increases on \((-\infty, 0]\) and \(\rho(wx_k + b) < wx_k + b < 0\), \(\phi_s(x) < \phi(x)\) for \(x < x_k\). The second step is to translate \(\phi_s(x)\) to \(\phi_t(x)\) by adjusting the bias of \(\phi_s(x)\), subject to \(\phi_t(x_k) = \phi(x_k)\) of equation 3.22; by \[w''x_k + b'' = \rho(wx_k + b) + \gamma' = wx_k + b,\] we obtain the bias modification value \(\gamma' = (wx_k + b)(1-\rho)\) and thus \(w'' = \rho w\) and \(b'' = (1-\rho)wx_k + b\) in equation 3.21.

See the red cure of Figure 1a, illustrating the difference between \(\phi_t(x)\) and \(\phi(x)\). Since \[(w''x + b'') - (wx + b) = (\rho-1)w(x-x_k) < 0\] and \(wx + b < 0\) when \(x < x_k\), we have \[\phi_t(x) = \sigma(w''x + b'') < \phi(x) = \sigma(wx + b)\] for \(x < x_k\), which proves inequality 3.23.

By inequality 3.23, it is obvious that: \[\lim_{\rho(\gamma') \to +\infty}s_{k}(x) = 0\;\text{for}\;x \in (x_{k-1}, x_{k}),\] without including the right endpoint \(x_k\) of \((x_{k-1}, x_{k})\), where \(\rho(\gamma')\) is similar to \(\rho(\gamma)\) of equation 3.18, \[\lim_{\rho(\gamma') \to +\infty}s_{\tau}(x) = 0\;\text{for}\;x \in (x_{\tau-1}, x_{\tau}],\] where \(\tau = 2, 3, \dots, k-1\), and \[\lim_{\rho(\gamma') \to +\infty}s_{1}(x) = 0\;\text{for}\;x \in [0, x_1].\] A difference between equations 3.27 and 3.28 is that the interval of the former doesn’t contain the right endpoint, because \(\phi_t(x_k) = \varepsilon\) (equation 3.22) is fixed by the previous parameter setting; this problem is to be dealt with later. By the recurrence formula of equation 3.17, we have \[s_{\tau}(x) = s_1(x) + \sum_{j = 1}^{\tau-1}c_{j}(x-x_j)^m_+\] with \(\tau\) from equation 3.28. Equations from 3.28 to 3.30 yield the following equations \[s_2(x) = s_1(x) + c_1(x-x_1)^m_+ \approx 0\] for \(x \in (x_1, x_2]\), where “\(\approx\)” means “equal as \(\rho(\gamma') \to +\infty\)” (similarly for the remaining ones in this proof), \[s_3(x) = s_1(x) + c_1(x-x_1)^m_+ + c_2(x-x_2)^m_+ \approx 0\] for \(x \in (x_2, x_3]\), \(\dots\), \[s_{k-1}(x) = s_1(x) + \sum_{j = 1}^{k-2}c_{j}(x-x_j)^m_+ \approx 0\] for \(x \in (x_{k-2}, x_{k-1}]\). In Equation 3.31, because there are infinitely many \(x \in (x_1, x_2]\) fulfilling this formula, we can regard \(s_1(x)\) and \(-c_1(x-x_1)^m_+\) as almost the same polynomial; so \(s_1(x) \approx -c_1(x-x_1)^m\) also holds on \([0, x_1]\); we know \(s_1(x) \approx 0\) on \([0, x_1]\) (equation 3.29) and \(c_1 \approx 0\) follows. Similarly, in equation 3.32, \(s_1(x) + c_1(x-x_1)^m_+ \approx 0\) implies \(c_2 \approx 0\). This process can be repeatedly done until \[c_{\kappa} \approx 0\] for all \(\kappa = 1, 2, \dots, k-2\).

To the case of \(c_{k-1}\), we should solve the aforementioned problem that equation 3.27 is not applicable to the right endpoint \(x_k\) of \((x_{k-1}, x_k)\). First add a knot \(x_k'\) with \[x_{k-1} < x_k' < x_k;\] correspondingly, the new set of the knots is denoted by \(\Delta'\). Then construct a spline \(\mathcal{S}'(x) \in \mathfrak{S}^m_1(\Delta')\) approximating \(\phi_t(x)\). Let \(s_{k}'(x)\) be the polynomial of \(\mathcal{S}'(x)\) on \((x_{k}', x_k]\) and the notations of other polynomials on the rest of the subintervals remain the same as those of \(\mathcal{S}(x)\). When \(\rho(\gamma')\) is sufficiently large, we have \[\lim_{\rho(\gamma') \to +\infty}s_{k}(x) = 0\;\text{for}\;x \in (x_{k-1}, x_{k}']\] as the counterpart of equation 3.27, with the subinterval including its right end point \(x_k'\). Then \[\lim_{\rho(\gamma') \to +\infty}c_{k-1} = 0\] holds by the method of equation 3.34, where \(c_{k-1}\) is from \[s_k(x) = s_{k-1}(x) + c_{k-1}(x-x_{k-1})^m_+\] for \(x \in (x_{k-1}, x_k']\).

To a small enough positive real number \[\phi_t(x_k') = \epsilon < \varepsilon,\] equation 3.36 can be expressed as \[|s_{k}(x)| \le \epsilon.\] Note that to a fixed \(\epsilon\), equation 3.40 holds for arbitrary large \(x_k' \in (x_{k-1}, x_k)\) as \(\rho(\gamma')\) increases; that is, \(\sup{x_k'} = x_k\) but \(x_k' \ne x_k\), or equivalently, \[\lim_{\rho(\gamma') \to +\infty}{\{x_k - x_k'\}} = 0.\]

We also have \[s_k'(x) = s_k(x) + c_{k-1}'(x-x_k')^m_+\] for \(x \in (x_{k}', x_k]\). In equation 3.42, because \(s_k(x)\) approaches zero in terms of a polynomial of order \(m\) (equations 3.36 and 3.38), we can regarded it as constant zero; so equation 3.42 can be written as \[s_k'(x) \approx 0 + c_{k-1}'(x-x_k')^m_+.\]

Given a small positive real number \(\epsilon < \varepsilon\), suppose that \(\rho(\gamma')\) has been big enough to make \(\phi_t(x) \le \epsilon\) on \([0, x_k']\). Then due to the constraint of both \(\phi_t(x_k') = \epsilon\) (equation 3.39) and \(\phi_t(x_k) = \varepsilon\) (equation 3.22), \(s_k'(x)\) on \((x_k', x_{k}]\) should increase from \(\epsilon\) to \(\varepsilon\), that is, \[s_k'(x_k) - s_k'(x_k') = \varepsilon - \epsilon.\] When the length of interval \((x_k', x_{k}]\) becomes smaller as \(\rho(\gamma') \to +\infty\), under the fixed polynomial degree \(m\), \(c_{k-1}'\) of equation 3.43 should become larger to achieve the goal of equation 3.44, resulting in \[\lim_{\rho(\gamma') \to +\infty}c'_{k-1} = +\infty\] and the monotonically increasing of \(c_{k-1}'\) with \(\rho\).

In Figure 1b, the red cures with labels from 1 to 6 are six instantiations of \(\phi_t(x)\) whose parameter \(\rho\) monotonically increases, from which we can see the behavior of \(\phi_t(x)\) with respect to \(\rho(\gamma')\) and obtain an intuitive understanding of equation 3.45.

To the last step, we translate \(\phi_t(x)\) to \(\phi_a(x)\) by adjusting the bias in terms of \[\phi_a(x) = \sigma(w'x + b') = \sigma(w''x + b'' + \gamma''),\] where \(w''\) and \(b''\) is from equation 3.21, to satisfy \[\phi_a(x_k) = \phi_t(x'_k) = \epsilon,\] which implies \[w'x_k + b' = w''x_k' + b''\] or \[\rho(wx_k + b) + \gamma' + \gamma'' = \rho(wx_k' + b) + \gamma',\] leading to \[\gamma'' = \rho w(x_k' - x_k) < 0;\] thus, in equation 3.16, \[\gamma = \gamma'' + \gamma'.\] We next delete the knot \(x_k'\) and reconstruct the polynomials on the subintervals of \((x_{k}, 1]\); adjusting the positions of the knots bigger than \(x_k\) or adding new knots on \((x_{k}, 1]\) if necessary. Then the new spline \(\mathcal{S}(x) \in \mathfrak{S}^m_1(\Delta)\) can be regarded as zero on \([0, x_k]\) and \(s_{k+1}(x) \approx c_{k-1}'(x-x_k')^m_+\); let \(c_k = c_{k-1}'\) and this proves equation 3.20 by equation 3.45. In Figure 1a, the green curve is \(\phi_a(x)\) obtained by slightly translating \(\phi_t(x)\). This completes the proof.

The proof is by lemma 3. Equation 3.52 of this theorem corresponds to equation 3.19 of lemma 3; the output weight \(\lambda\) can be obtained by \[\lambda c_k = d_k,\] through which \(p_{k+1}(x)\) on \((x_k, x_{k+1}]\) is realized by \(\lambda\phi(x)\). To the part \([0, x_k]\), the approximation of \(\lambda\phi(x)\) to zero is via equations 3.18 and 3.20. A key point is that if \(|\lambda|\) increases as \(\rho(\gamma) \to +\infty\), a larger approximation error would occur on \([0, x_k]\); that’s why we developed the second conclusion of lemma 3. According to lemma 3, \(c_k > 0\) monotonically increases with \(\rho(\gamma)\); then in equation 3.53, \(|\lambda|\) decreases as \(c_k\) becomes larger; in combination with the smaller \(\phi(x)\), \(|\lambda\phi(x)|\) is reduced as \(\rho(\gamma)\) increases, resulting in a smaller approximation error on \([0, x_k]\).

3.4 Univariate Local (Global) Unit↩︎

By the proof of theorem 5, on \([0, x_k]\), \[\lim_{\rho(\gamma) \to +\infty}\{\int_{0}^{x_k}(\lambda\phi(x))^2dx\}^{1/2} = 0.\] which means that truncating \(\phi(x)\) on \([0, x_k]\) would not influence the establishment of inequality 3.55 when \(\rho(\gamma)\) is sufficiently large. Due to the proof of lemma 3, \(\lim_{\varepsilon \to 0}z_i(\varepsilon) = x_k \ne 0\) always holds, no matter how small \(\varepsilon\). This completes the proof.

3.5 One-Sided Solutions↩︎

To approximate univariate functions, the term “one sided” means that all \(w_i\)’s of activation functions \(\phi_i(x) = \sigma(w_ix + b_i)\) are greater than (or less than) zero; otherwise, it is the “two-sided” case.

We first see an example of Figure 2. In this case, \(f(x) = s(x)\) is a spline. Suppose that \(\phi_4(x)\) of \(u_4\) is initially absent, and then \(g(x) = \sum_{i=1}^{3}\lambda_i\phi_i(x)\). Because \(\phi_i(x) \in C^{m}([0, 1])\), it can be approximated by a spline \(s_i(x) \in \mathfrak{S}^m_1(\Delta)\) as precisely as possible, with \(p_{k, i}(x)\) the polynomial on the \(k\)th subinterval satisfying \[p_{k, i}(x) = p_{k-1, i}(x) + \beta_{k-1,i}(x - x_{k-1})_+^{m}.\]

The polynomial \(\mathcal{P}_{k-1}(x)\) of inequality 3.64 is the linear combination of \(p_{k-1, i}(x)\)’s for all \(i\), that is, \[\mathcal{P}_{k-1}(x) = \sum_{i}\lambda_ip_{k-1, i}(x).\] Because \(\sum_{i=1}^{3}\lambda_ip_{k, i}(x) \ne \mathcal{P}_k(x)\), inequality 3.66 doesn’t hold, for which we add a new unit to change this situation. In Figure 2, \(u_4\) is the added unit whose activation function is \(\phi_4(x)\). There are two main steps to adjust the parameters of \(\phi_4(x)\) to satisfy \(\sum_{i=1}^{4}\lambda_ip_{k, i}(x) = \mathcal{P}_k(x)\). The first is by lemma 3 to reduce the influence of \(\phi_4(x)\) on \([0, x_{k-1}]\) to a desired accuracy. The second is to set \(\lambda_4\) such that \(\sum_{i=1}^3\lambda_{i}p_{k, i}(x) + \lambda_4p_{k, 4}(x) = \mathcal{P}_k(x)\). Let \(\zeta = 3\). The original polynomial \(\mathcal{P}_k(x)\) before adding \(u_{\zeta + 1}\) can be expressed as \[\begin{align} \sum_i^{\zeta} \lambda_{i} & p_{k,i}(x) = \sum_{i=1}^{\zeta} \lambda_i\big(p_{k-1,i}(x) + \beta_{k-1,i}(x-x_{k-1})_+^{m}\big) \\ &= \mathcal{P}_{k-1}(x) + \gamma(x - x_{k-1})_+^{m} \ne \mathcal{P}_{k}(x), \end{align}\] where \[\gamma = \sum_{i=1}^{\zeta}\lambda_i\beta_{k-1,i}.\]

We have two additional conditions. One is \[\mathcal{P}_{k}(x) = \mathcal{P}_{k-1}(x) + c(x - x_{k-1})_+^{m}.\] The other is that equation 3.67 is applicable to the the new unit \(u_{4}\) as well, that is, \[p_{k,\zeta+1}(x) = p_{k-1, \zeta+1}(x) + \beta_{k-1,\zeta+1}(x - x_{k-1})_+^{m}.\] Then the updated version of equation 3.69 is \[\mathcal{P}'_{k}(x) = \sum_{\mu = 1}^{\zeta + 1} \alpha_{\mu}p_{k{\mu}}(x) = \mathcal{P}'_{k-1}(x) + \gamma'(x - x_{k-1})_+^{m},\] where \[\mathcal{P}'_{k-1}(x) = \mathcal{P}_{k-1}(x) + \lambda_{\zeta+1}p_{k-1,\zeta+1}(x)\] and \[\gamma' = \gamma + \lambda_{\zeta+1}\beta_{k-1,\zeta+1}.\] In equation 3.74, setting \(p_{k-1,\zeta+1}(x)\approx 0\) by lemma 3, where “\(\approx\)” means “approximately equal with arbitrary precision”, we have \[\mathcal{P}'_{k-1}(x) \approx \mathcal{P}_{k-1}(x).\] Since \(\lambda_{\zeta+1}\) of equation 3.75 can be arbitrarily set, we can find a solution of \(\lambda_{\zeta+1}\) to satisfy \(\gamma' = c\) of equation 3.71, which is \[\lambda_{\zeta+1} = (c - \gamma)/(\beta_{k-1,\zeta + 1}).\] Equations from 3.73 to 3.77 realize \(\mathcal{P}'_{k}(x) \approx \mathcal{P}_{k}(x)\). This completes the proof.

According to lemma 2, the first polynomial \(\mathcal{P}_1(x)\) on \([0, x_1]\) of \(s(x)\) can be realized by a set of \(\phi_j(x)\)’s for \(j = -\mu, -\mu+1, \dots, 0\) with \(\mu+1 \ge m+1\). The remaining polynomials are implemented by lemma 4, with each requiring one unit. So the total number of the units is at least \[(m + 1) + (\zeta - 1) = \zeta + m.\] The approximation error is controlled by theorem 5 and lemma 4.

Theorem 4 gave the method of constructing spline \(\hat{f}(x)\) approximating \(f(x)\) and theorem 6 realized it via two-layer neural networks.

Denote the new neural network by \(\mathfrak{N}'\). Under theorem 7, the polynomial \(s_i(x)\) on the \(i\)th subinterval for \(\hat{f}(x)\) can be expressed as \[s_i(x)=\sum_{j=-\mu}^{0}\lambda_j\phi_j(x)+\sum_{j=1}^{i-1}\lambda_j\phi_j(x),\] where \(\phi_j(x)\)’s are the activation function of the hidden-layer units of \(\mathfrak{N}\), among which \(\mu+1\) ones for \(\mu+1\ge m+1\) and \(-\mu\le j\le0\) are of global unit and the remaining ones are of local unit. Letting \(\tau_j(x)\) be the activation function of a generalized tanh unit corresponding to \(\phi_j(x)\), then \[\psi_j(x) = \tau_j(x)-\mathcal{C}\] is of a generalized sigmoidal unit, where \(\mathcal{C}\) is a constant from definition 4. So \(s_i(x)\) of equation 3.81 can also be represented as \[\begin{align} s_i(x)=\sum_{j=-\mu}^{0}\alpha_j\tau_j(x&)+\sum_{j=1}^{i-1}\alpha_j\psi_j(x) = \sum_{j=-\mu}^{0}\alpha_j\tau_j(x)+\sum_{j=1}^{i-1}\alpha_j(\tau_j(x)-\mathcal{C})\\ &=\sum_{j=-\mu}^{0}\alpha_j\tau_j(x) + \sum_{j=1}^{i-1}\alpha_j\tau_j(x) -\sum_{j=1}^{i-1}\alpha_j\mathcal{C}. \end{align}\]

Especially, when \(i=\zeta\), we have \[s_{\zeta}(x)= \sum_{j=-\mu}^{0}\alpha_j\tau_j(x) + \sum_{j=1}^{\zeta-1}\alpha_j\tau_j(x) + \mathfrak{C},\] where \[\mathfrak{C}=-\sum_{j=1}^{\zeta-1}\alpha_j\mathcal{C}\] is a constant when \(\alpha_j\)’s for \(1\le j\le \zeta-1\) are fixed. If \(i < \zeta\), the new added \(-\alpha_{i+1}\mathcal{C}\) would not influence \(s_{i}(x)\) because \(\tau_{i+1}(x)-\mathcal{C}\approx 0\); so \(\mathfrak{C}\) stop change when \(i=\zeta\), and its update when \(i < \zeta\) cannot affect the previously established result. The polynomial \(s_1(x)\) of \(\hat{f}(x)\) on \([0, x_1]\) can be implemented by a set of any type of smooth activation (theorem 1), for which we can replace \(\phi_i(x)\) of equation 3.81 by \(\tau_j(x)\) of equation 3.83 when \(-\mu\le j\le0\).

Therefore, if the constant \(\mathfrak{C}\) of equation 3.85 can be produced by generalized tanh units, the universal approximation of \(\mathfrak{N}'\) could be ensured. We provide two solutions. One is that add a unit whose activation function over \([0, 1]\) is nearly a constant \(\beta \ne 0\) and set its output weight to be \(\mathfrak{C}/\beta\). The other is by adding some global units whose linear combination is \(\mathfrak{C}\). This completes the proof.

3.6 A Two-Sided Solution↩︎

Figure 3 gives an example of this lemma, in which \(\phi'_4(x)\) is the activation function of new added unit \(u'_{4}\) whose zero-error point is \(x_2=0.4\). After \(u'_{4}\) having been introduced, the output weights \(\lambda_i\)’s for \(i = 1, 2, 3\) of \(\phi_i(x)\)’s should be altered to compensate the impact of \(\phi'_4(x)\) on \([0, x_2]\); this is not difficult by lemmas 2 and 4. The parameter resetting could influence the polynomial on \((x_2, x_3]\), but the output weight \(\lambda_4\) of \(\phi_4(x)\) can be adjusted to compensate the disturbance (lemma 4), such that the original polynomial can still be generated. This example contains the general principle of this lemma.

This is the repeated application of lemma 5.

4 General Solution Framework↩︎

The solutions constructed in section 3 may not include all those obtained by the back-propagation algorithm. So we develop a general framework by reducing the spline implementation to a system of linear equations, through which the diversity of solutions can be attributed to the nonsingular property of matrices. We first show two examples derived from section 3 and then give the general form.

4.1 One-Sided Example↩︎

This section describes theorem 6 through matrix operations. Let \(g(x) = \sum_{i = 1}^{\Theta}\lambda_i\phi_i(x)\) be the output of a two-layer network \(\mathfrak{N}\) with generalized sigmoidal units, where \(\phi_i(x) = \sigma_i(w_ix + b_i) \in C^{(m)}(\mathbb{R})\) is the activation function of the \(i\)th unit of the hidden layer. The set \(\Delta = \{x_{j}: j = 1, 2, \dots, \zeta-1\}\) includes all the knots. Denote by \(u_{k}\)’s for \(k = 1, 2, \dots, \mu\) the global units, with \(\mu \ge m+1\), and by \(u_{j}\)’s for \(j = \mu+1, \mu+2, \dots, \mu + \zeta-1\) the local units.

By theorem 4, each activation function \(\phi_i(x)\) has an approximation via a spline \(s_i(x) \in \mathfrak{S}^m_1(\Delta)\) whose polynomial on subinterval \(I_{\nu}\) is denoted by \(p_{i\nu}(x)\) for \(\nu = 1, 2, \dots, \zeta\). Then \[p_{ij}(x) = p_{i, j-1}(x) + \beta_{i,j-1}(x - x_{j-1})_+^{m}\] for \(j = 2, 3, \dots, \zeta\). Let \(\mathcal{S}(x) \in \mathfrak{S}^m_1(\Delta)\) be the spline constructed from \(f(x)\), which is to be implemented by network \(\mathfrak{N}\), and let \(\mathcal{P}_{\nu}(x)\) be the polynomial of \(\mathcal{S}(x)\) on \(I_{\nu}\). We also have \[\mathcal{P}_{j}(x) = \mathcal{P}_{j-1}(x) + c_{j-1}(x - x_{j-1})_+^{m}.\] The goal is to realize \(\mathcal{P}_{j}(x)\)’s by the weighted sum of \(p_{ij}(x)\)’s.

According to section 1, the way of producing \(\mathcal{P}_1(x)\) on \(I_1 = [0, x_1]\) is by \[\mathcal{P}_1(x) = \sum_{k \in G \;\text{or} \; k\le\mu}\lambda_{k}p_{k1}(x),\] where \(G\) is the set of the subscripts of the global units, leading to \[A_1\boldsymbol{\lambda}_1 = \boldsymbol{b}_1\] as equation 2.9, where \(A_1\) is a \((m+1) \times \mu\) matrix, \(\boldsymbol{\alpha}_1 = [\lambda_{1}, \lambda_{2}, \dots, \lambda_{\mu}]^T\) and \(\boldsymbol{b}_1\) is a \((m+1) \times 1\) vector whose entries are the coefficients of the terms of \(\mathcal{P}_1(x)\).

Polynomial \(\mathcal{P}_2(x)\) is obtained by \[\mathcal{P}_{2}(x) = \mathcal{P}_1(x) + \big(\sum_{k \in G}\lambda_{k}\beta_{k1} + \lambda_{\mu+1}\beta_{\mu+1, 1}\big)(x - x_{1})_+^{m}.\] We make \[\sum_{k \in G}\lambda_{k}\beta_{k1} + \lambda_{\mu+1}\beta_{\mu+1, 1} = c_1\] of equation 4.2 by adjusting the free parameter \(\lambda_{\mu+1}\). The remaining \(\mathcal{P}_{j}(x)\)’s are similar. Through \[\mathcal{P}_{j}(x) = \mathcal{P}_{j-1}(x) + \big(\sum_{k < \mu + j-1}\lambda_{k}\beta_{k,j-1} + \lambda_{\mu + j-1}\beta_{\mu + j-1, j-1}\big)(x - x_{j-1})_+^{m},\] we obtain \[\sum_{k < \mu + j-1}\lambda_{k}\beta_{k,j-1} + \lambda_{\mu + j-1}\beta_{\mu + j-1, j-1} = c_{j-1}\] and set \(\lambda_{\mu + j-1}\) to generate \(c_{j-1}\).

Equations 4.4 and 4.8 can be combined to a integrated system of linear equations \[\mathcal{A}\boldsymbol{\lambda} = \boldsymbol{b},\] where

\[\mathcal{A} = \begin{pmat}({ccc|}) &A_1 & & & &\boldsymbol{0} \cr\- \beta_{11} & \beta_{21} & \cdots & \beta_{\mu,1} &\beta_{\mu+1,1} &0 &\cdots &0\cr \beta_{12} & \beta_{22} & \cdots & \beta_{\mu,2} &\beta_{\mu+1, 2} &\beta_{\mu+2, 2} &\cdots &0\cr \vdots & \vdots & \ddots & \vdots &\vdots &\vdots &\ddots &\vdots \cr \beta_{1, \zeta-1} & \beta_{2, \zeta-1} & \cdots & \beta_{\mu, \zeta-1} & \beta_{\mu + 1,\zeta-1} & \beta_{\mu + 2, \zeta-1} &\cdots &\beta_{\mu + \zeta-1,\zeta-1}\cr \end{pmat}\] is a \((\zeta + m) \times (\mu + \zeta - 1)\) matrix and \(\boldsymbol{b}\) is a \((\zeta + m) \times 1\) vector including the entries of \(\boldsymbol{b}_1\) in equation 4.4 as well as \(c_{j-1}\)’s of equation 4.8. Equation 4.10 can be written as \[\mathcal{A} = \begin{pmat}({}) A_1 & & & B_1 \cr C_1 & & & D_1 \cr \end{pmat} = \begin{pmat}({}) A_1 & & & \boldsymbol{0} \cr C_1 & & & D_1 \cr \end{pmat},\] where \(C_1\) and \(D_1\) equals the corresponding blocks of \(\mathcal{A}\), respectively.

Since \(\det{A} = \det{A_1}\det{D_1} = \det{A_1}\prod_{\nu=1}^{\zeta-1}\beta_{\mu+\nu, \nu}\), the conclusion follows.

4.2 A Two-Sided Example↩︎

This section corresponds to the solution of section 3.6. Under the notations of section 3.6, to a unit \(u_{k}\) , add a negative unit \(u_{k}'\) whose zero-error point is knot \(x_k\). Let \(\phi_{k}(x) = \sigma(w'_{k}x + b_{k}')\) be the activation function of \(u_{k}'\), and \(p_{k\tau}'(x)\) for \(\tau = 1, 2, \dots, \zeta\) be the \(\tau\)th polynomial of a spline approximating \(\phi_{k}'(x)\) on subinterval \(I_{\tau}\), satisfying \[p_{k j}'(x) = p_{k,j-1}'(x) + \beta_{k,j-1}'(x - x_{j-1})_+^{m}.\] for \(2 \le j \le k\). After introducing \(u_k'\), matrix \(\mathcal{A}\) of equation 4.11 is changed into \[\mathcal{A}' = \begin{pmat}({}) A_1' & & & B_1' \cr C_1' & & & D_1' \cr \end{pmat},\] whose blocks correspond to those of equation 4.11, respectively, among which \(A_1'=A_1\) and \(C_1'=C_1\). Blocks \(B_1'\) and \(D'_1\) is obtained as follows. Equation 4.3 becomes \[\mathcal{P}_1(x) = \sum_{\nu \in G}\lambda_{\nu}p_{\nu1}(x) + \lambda_{k}'p_{k1}'(x)\] and correspondingly, the submatrix \(B_1\) of equation 4.11 should be modified to \[B_1' = \begin{pmat}({cc||}) 0 &\cdots &0 & \gamma_{k1} &0 &\cdots &0 \cr\; 0 &\cdots &0 &\gamma_{k2} &0 &\cdots &0 \cr\; \vdots &\ddots &\vdots &\vdots &\vdots &\ddots &\vdots \cr\; 0 &\cdots &0 &\gamma_{k,m+1} &0 &\cdots &0 \cr \end{pmat},\] by adding a column between the original \(k-1\)th and \(k\)th columns, whose entries \(\gamma_{ki}\)’s for \(1 \le i \le m+1\) are derived from the coefficients of the polynomial \(p_{k1}'(x)\) of equation 4.14 weighted by \(\lambda_{k}'\). The submatrix \[\begin{align} &\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad D_1' = \\ &\begin{pmat}({ccc||}) \beta_{\mu+1,1} &0 &0 &\cdots &\beta_{\mu+k,1}'&0 &0 &\cdots &0 \cr \beta_{\mu+1,2} & \beta_{\mu+2,2} &0 &\cdots &\beta_{\mu+k,2}'&0 &0 &\cdots &0\cr \vdots &\vdots &\vdots &\ddots &\vdots &\vdots &\ddots &\vdots \cr \beta_{\mu+1,k} & \beta_{\mu+2,k} & \cdots &\beta_{\mu+k-1,k} &\beta'_{\mu+k,k} &0 &0 &\cdots &0\cr \beta_{\mu+1,k+1} & \beta_{\mu+2,k+1} & \cdots &\beta_{\mu+k-1,k+1} &0 &\beta_{\mu+k,k+1} &0 &\cdots &0\cr \vdots &\vdots &\vdots &\ddots &\vdots &\vdots &\ddots &\vdots \cr \beta_{\mu + 1,\zeta-1} & \beta_{\mu + 2,\zeta-1} & \cdots &\beta_{\mu+k-1,\zeta-1} &0 &\beta_{\mu+k,\zeta-1} &\beta_{\mu+k+1,\zeta-1} &\cdots &\beta_{\psi,\zeta-1}\cr \end{pmat}, \end{align}\] where \(\psi = \mu+ \zeta-1\), which is different from \(D_1\) of equation 4.10 in adding a column \[[\beta_{\mu+k,1}', \beta_{\mu+k,2}', \dots, \beta'_{\mu+k,k}, \dots,0,0]^T.\]

The reason is that it has a nonsingular square submatrix \(\mathcal{A}\) of equation 4.11 obtained by deleting the added column in \(B_1'\) and \(D_1'\).

4.3 General Solution Form↩︎

The proof is the generalization of sections 4.1 and 4.2.

5 Multivariate Local Approximation↩︎

In lemma 2 of section 1, Taylor series expansions of a univariate function were realized by a two-layer neural network with smooth activation functions, and this section generalizes that idea to the multivariate case. The obtained result is the basis of multivariate-spline construction of section 6.

5.1 Theoretical Framework↩︎

We require the number of the coefficients of a multivariate polynomial. To an \(n\)-variate homogeneous polynomial of degree \(m\), the number of its coefficients is \(\mathscr{N}(n, m) = \binom{n+m-1}{m}\), whose calculation can be reduced to the combination problem of multisets ([26]’s theorem 3.5.1). The number of the coefficients of an \(n\)-variate polynomial \(p(\boldsymbol{x})\) of degree \(m\) is \[\mathcal{N}(n, m) = \binom{n+m}{m},\] because \(p(\boldsymbol{x})\) comprises \(m+1\) sets of homogeneous polynomials whose degrees are \(m, m-1, \dots, 0\), respectively, such that \[\mathcal{N}(n, m) = \binom{n+m-1}{m} + \binom{n+m-2}{m-1} + \dots + \binom{n}{1} + \binom{n-1}{0} = \binom{n+m}{m},\] by the identity \(\binom{n-1}{0} = \binom{n}{0}\) and by repeatedly applying Pascal’s formula from the right side of the middle term.

Note that in equation 5.4, the term \(\partial x_1^{\alpha_1}\partial x_2^{\alpha_2}\dots\partial x_n^{\alpha_n}\) has an one-to-one correspondence with \(x_1^{\alpha_1}x_2^{\alpha_2}\dots x_n^{\alpha_n}\). So the conclusion is obvious by equation 5.1.

It’s easy to verify:

The proof is similar to that of theorem 1 and the difference lies in multivariate Taylor series expansions. So we have \[\mathcal{W}(u_i(\boldsymbol{x}_0))^T\boldsymbol{\lambda} = \boldsymbol{a},\] corresponding to equation 2.13 of the univariate case.

5.2 Construction of Local Approximation↩︎

The proof is composed of three parts. Part 1 illustrates the main idea by an example. Part 2 deals with the order of the terms of a multivariate polynomial. Part 3 proposes the general construction method.

Part 1. For simplicity, we use the following notations in the next equation 5.11: \(\sigma'_{x_{\nu}}(y):=\partial\sigma(y)/\partial x_{\nu}\) for \(\nu = 1, 2\), \(\sigma''_{x_{\nu}x_{\mu}}(y):=\partial^2\sigma(y)/\partial x_{\nu}^{\alpha_1}\partial x_{\mu}^{\alpha_2}\) for \(\nu, \mu = 1, 2\), with \(\alpha_1 + \alpha_2 = 2\) and \(0 \le \alpha_1, \alpha_2 \le 2\). Analogous to lemma 2, write \[\begin{align} &\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \mathcal{W}(\sigma(y_i)) = \\ &\begin{pmat}({}) \sigma(y_1) &\sigma'_{x_1}(y_1)w_{11} &\sigma'_{x_2}(y_1)w_{12} &\sigma''_{x_1^2}(y_1)w_{11}^2 &\sigma''_{x_1x_2}(y_1)w_{11}w_{12} & \sigma''_{x_2^2}(y_1)w_{12}^2 \cr \sigma(y_2) &\sigma'_{x_1}(y_2)w_{21} &\sigma'_{x_2}(y_2)w_{22} &\sigma''_{x_1^2}(y_2)w_{21}^2 &\sigma''_{x_1x_2}(y_2)w_{21}w_{22} & \sigma''_{x_2^2}(y_2)w_{22}^2 \cr \sigma(y_3) &\sigma'_{x_1}(y_3)w_{31} &\sigma'_{x_2}(y_3)w_{32} &\sigma''_{x_1^2}(y_3)w_{31}^2 &\sigma''_{x_1x_2}(y_3)w_{31}w_{32} & \sigma''_{x_2^2}(y_3)w_{32}^2 \cr \sigma(y_4) &\sigma'_{x_1}(y_4)w_{41} &\sigma'_{x_2}(y_4)w_{42} &\sigma''_{x_1^2}(y_4)w_{41}^2 &\sigma''_{x_1x_2}(y_4)w_{41}w_{42} & \sigma''_{x_2^2}(y_4)w_{42}^2 \cr \sigma(y_5) &\sigma'_{x_1}(y_5)w_{51} &\sigma'_{x_2}(y_5)w_{52} &\sigma''_{x_1^2}(y_5)w_{51}^2 &\sigma''_{x_1x_2}(y_5)w_{51}w_{52} & \sigma''_{x_2^2}(y_5)w_{52}^2 \cr \sigma(y_6) &\sigma'_{x_1}(y_6)w_{61} &\sigma'_{x_2}(y_6)w_{62} &\sigma''_{x_1^2}(y_6)w_{61}^2 &\sigma''_{x_1x_2}(y_6)w_{61}w_{62} & \sigma''_{x_2^2}(y_6)w_{62}^2 \cr \end{pmat}, \end{align}\] where \(y_i = \boldsymbol{w}_i^T\boldsymbol{x}_0 + b_i\) for \(1 \le i \le 6\) and \(w_{ij}\) for \(j = 1, 2\) is the entry of \(\boldsymbol{w}_i = [w_{i1}, w_{i2}]^T\), which is the generalized Wronskian matrix of \(\phi_i(x)\)’s with respect to order \(2\). To the diagonal entries of \(\mathcal{W}(\sigma(y_i))\), the partial derivatives of \(\sigma(y)\) on a fixed point could be nonzero by adjusting the bias parameter; thus we only consider the construction of \(w_{ij}\)’s to make \(\mathcal{W}(\sigma(y_i))\) approximately nonsingular.

A key point is that only the diagonal and upper-triangular part of \(\mathcal{W}(\sigma(y_i))\) should be taken into consideration, while the lower-triangular part can be ignored. The ultimate goal is to set each diagonal (except for the first and last ones) entry to be \(\Delta t\), and to set each upper-triangular entry to be of \(o(\Delta t)\). Another point is that, with \(w_{ij}\)’s as variables, the polynomial degrees of the diagonal entries are the same as those of the entries of each row, respectively, as shown in an example of Table 1.

The two points above play an important role in the proof of this lemma.

For simplicity of descriptions, let \[W = \begin{pmat}({}) 1 & w_{11} & w_{12} & w_{11}^2 & w_{11}w_{12} & w_{12}^2 \cr 1 & w_{21} & w_{22} & w_{21}^2 & w_{21}w_{22} & w_{22}^2 \cr 1 & w_{31} & w_{32} & w_{31}^2 & w_{31}w_{32} & w_{32}^2 \cr 1 & w_{41} & w_{42} & w_{41}^2 & w_{41}w_{42} & w_{42}^2 \cr 1 & w_{51} & w_{52} & w_{51}^2 & w_{51}w_{52} & w_{52}^2 \cr 1 & w_{61} & w_{62} & w_{61}^2 & w_{61}w_{62} & w_{62}^2 \cr \end{pmat},\] regardless of the coefficients of variable \(w_{ij}\)’s in equation 5.11. We deal with \(W\) from the first row; just letting \(w_{11} = w_{12} = \Delta t^{1 + c}\) with \(c > 0\), all the entries except for the first one is \(o(\Delta t)\). To the second row, set the diagonal entry \(w_{21} = \Delta t\); if \(w_{22} = \Delta t^{1 + c}\), \(W(2, j)\)’s for \(j \ge 3\) would be \(o(\Delta t)\). The third row is similar to the second one. In the fourth row, let \(w_{41} = \Delta t^{1 / 2}\) and \(w_{42} = \Delta t^{\beta}\) , where \(\beta\) can be arbitrarily set as long as \(W(4, 5)\) and \(W(4, 6)\) are \(o(\Delta t)\). The above cases are all trivial, because the diagonal entry has only one variable such that other variables can be freely adjusted without affecting it.

The more difficult case is the fifth row when \(w_{52}\) appears in both diagonal and upper-triangle positions. We should set \(w_{52}\) properly to make \(W(5, 5) = \Delta t\) and \(W(5, 6) = o(\Delta t)\) simultaneously. Let \(w_{51} = \Delta t^{c}\) for \(c > 0\) and \(w_{52} = \Delta t^{1 - c}\); then \(w_{51}w_{52} = \Delta t\) and if \(c < 1/2\), \(w_{52}^2 = \Delta t^{2(1 - c)}\) whose exponential \(2(1 - c) > 1\) , that is, \(w_{52}^2 = o(\Delta t)\). The last row is trivial by setting \(w_{62} = 1/\Delta t^4\) or \(w_{62} = \Delta t^{-4}\). So we finally get \[W = \begin{pmat}({}) 1 & o(\Delta t) & o(\Delta t) & o(\Delta t) & o(\Delta t) & o(\Delta t) \cr * & \Delta t & o(\Delta t) & o(\Delta t) & o(\Delta t) & o(\Delta t) \cr * & * & \Delta t & o(\Delta t) & o(\Delta t) & o(\Delta t) \cr * & * & * & \Delta t & o(\Delta t) & o(\Delta t) \cr * & * & * & * & \Delta t & o(\Delta t) \cr * & * & * & * & * & 1/\Delta t^4 \cr \end{pmat},\] where the “\(*\)” entries mean that they cannot influence the determinant of the matrix when \(\Delta t \to 0\). So \(\lim_{\Delta t \to 0}\det W(\Delta t) = 1\). By this result, we have \[\lim_{\Delta t \to 0}\det \mathcal{W}(\sigma(y_i), \Delta t) = \alpha \ne 0,\] where \(\alpha\) is the product of the coefficients of the diagonal entries of \(\mathcal{W}(\sigma(y_i))\).

Part 2. Write a polynomial of degree \(2\) with three variables in this form \[p(\boldsymbol{w}) = p_0 + p_1 + p_2 + p_3,\] where \(\boldsymbol{w} = [w_1, w_2, w_3]^T\), \(p_0 = 1\), \(p_1 = w_1 + w_2 + w_3\), \[p_2 = \big\{w_1^2 + w_1(w_2 + w_3)\big\} + \big\{w_2^2 + w_2w_3\big\} + w_3^2,\] and \[p_3 = \big\{w_1^3 + w_1^2(w_2 + w_3) + w_1(w_2^2 + w_2w_3 + w_3^2)\big\} + \big\{w_2^3 + w_2^2w_3 + w_2w_3^2\big\} + w_3^3.\]

The order of the terms of \(p(\boldsymbol{w})\) in equation 5.15 is important to this lemma. Each \(p_i\) for \(i = 0, 1, 2, 3\) is a homogeneous polynomial of degree \(i\) and \(p_i\)’s are arranged in ascending order of their degrees. The terms of each \(p_i\) are also arranged in a certain order. For example, to \(p_2\) of equation 5.16, the first term \(w_1^2\) is a univariate-\(w_1\) polynomial of degree 2. Then reduce the degree of \(w_1^2\) and simultaneously increase that of other variables with the constraint that the total degree remains invariant, resulting in two terms \(w_1(w_2 + w_3)\). The order of introducing new variables in \(p_{2}\) is by the ascending order of the subscripts of \(w_i\)’s for \(i > 1\). When the degree of \(w_1\) is reduced to zero, variable \(w_2\) plays the role of \(w_1\) and the similar process of creating new terms is repeated as above. A difference is that when introducing new variables for \(w_2\), only those whose subscripts are greater than 2 should be considered. The last case \(w_3^2\) has only one term because there’s no variable left with subscript greater than 3.

A more complex example is \(p_3\) of equation 5.17. First dealt with variable \(w_1\) by repeatedly reducing its degree from 3 to 1 and simultaneously increasing the degree of other variables. Note that in \(w_1(w_2^2 + w_2w_3 + w_3^2)\), the terms of \(w_2^2 + w_2w_3 + w_3^2\) with degree 2 should be arranged by the rule of \(p_2\) discussed above; so this is a recursive procedure obeying a consistent principle. When the degree of \(w_1\) is zero, \(w_2\) follows and would be processed analogous to \(w_1\). And also, a difference from \(w_1\) is that only variable \(w_3\) whose subscript is greater than \(w_2\) should be introduced when decreasing the degree of \(w_2^3\). The last step is for \(w_3\), which is trivial as the case of \(p_2\).

We summarize a rule for the order of the terms of \(p(\boldsymbol{w})\) above, which is in fact the strict partial order “\(\prec\)” of definition 10. So to any \(n\)-variate polynomial \(p(\boldsymbol{w})\) of degree \(m\), we can arrange its terms in the form \(\{q_k: 1 \le k \le \tau\}\) such that \[(q_i = w_1^{k_1}w_2^{k_2}\dots w_n^{k_n}) \prec (q_{i+1} = w_1^{k_1'}w_2^{k_2'}\dots w_n^{k_n'})\] for \(i= 1,2,\dots, \tau-1\), where \(q_{i}\) and \(q_{i+1}\) are arbitrary two terms of \(p(\boldsymbol{w})\) without considering the coefficients.

Part 3. On the basis of the preceding two parts, we give the general construction method. Given a generalized Wronskian matrix \(\mathcal{W}\) with respect to order \(m\), the associated matrix \(W\) is obtained by ignoring the coefficients of the terms, analogous to equation 5.12. To each row of \(W\), the entries are regard as the terms of a polynomial and must be arranged by the order “\(\prec\)” of definition 10 as in equation 5.18.

The main procedure is as follows. We process each row of \(W\) one by one according to the ascending order of \(i = 1, 2, \dots, \tau = \binom{n+m}{m}\). The case of \(i = 1\) is trivial as the example of part 1. In the \(\nu\)th step for \(\nu = 2, \dots, \tau-1\), let \(W(\nu, \nu)=\Delta t\); and the terms bigger than \(W(\nu, \nu)\) are set to be \(o(\Delta t)\). Finally, let \(W(\tau, \tau) = 1/\Delta t^{\tau-2}\). The remaining proof is mainly for the second case.

Let \(p(\boldsymbol{w})\) be the polynomial of the \(\nu\)th row. For simplicity, denote the term \(W(\nu, \nu)\) of \(p(\boldsymbol{w})\) by \(w_1^{k_1}w_2^{k_2}\dots w_n^{k_n}\) or \(\prod_{j=1}^nw_j^{k_j}\), disregarding the variables’ row index (or the first subscript); similarly, each term of \(p(\boldsymbol{w})\) bigger than \(W(\nu, \nu)\) is denoted by \(w_1^{k_1'}w_2^{k_2'}\dots w_n^{k_n'}\).

When \(w_1^{k_1}w_2^{k_2}\dots w_n^{k_n}\) is a univariate polynomial with only one of \(k_j\)’s being nonzero, this case is trivial. For example, suppose that \(w_1^{k_1}w_2^{k_2}\dots w_n^{k_n} = w_{\mu}^k\) for \(1\le \mu \le n\). Just let \(w_{\mu} = \Delta t^{1/k}\). The terms of this row bigger than \(w_{\mu}^k\) have three possibilities: (a) it’s a univariate polynomial with variable \(w_{\mu}\) whose exponential is greater than \(k\) and thus is \(o(\Delta t)\); (b) it’s a multivariate polynomial having other variables; in this case we can set the exponential of any other variable to be \(\Delta t^{\alpha}\), where \(\alpha\) can be arbitrarily large such that the polynomial becomes \(o(\Delta t)\); (c) it’s a polynomial that does not contain \(w_{\mu}\) and can also be freely adjusted without any restriction.

The difficult case is when \[W(\nu, \nu)= w_1^{k_1}w_2^{k_2}\dots w_n^{k_n}\] is a multivariate polynomial. Denoted by \(S\) the set of the terms of \(p(\boldsymbol{w})\) bigger than \(w_1^{k_1}w_2^{k_2}\dots w_n^{k_n}\). Suppose that \[\sum_{i=1}^nk_i = k.\] The set \(S\) can be classified into three categories, including \[S = S_1 \cup S_2 \cup S_3,\] in which \[S_1 = \{w_1^{k_1'}w_2^{k_2'}\dots w_n^{k_n'}: 0<k_1'<k_1, \;\sum_{i=1}^nk_i' = k \},\] \[S_2 = \{w_1^{k_1'}w_2^{k_2'}\dots w_n^{k_n'}: 0<k_1' \le k_1, \;\sum_{i=1}^nk_i' > k \},\] and \[S_3 = S-(S_1 \cup S_2),\] where \(k_1\) is from equation 5.19.

Let \[w_1 = \Delta t^{\frac{1+c_1}{k+1}},\] and \[w_j = \Delta t^{\frac{1+c_j}{k}},\] for \(j = 2, 3, \dots, n\), where \[c_j > c_1 > 0.\] We should adjust the parameters of \(c_i\)’s for \(i=1,2,\dots,n\) such that \[w_1^{k_1}w_2^{k_2}\dots w_n^{k_n} = \Delta t\] and each element of \(S\) becomes \(o(\Delta t)\). To satisfy equation 5.28, we have \[\frac{1+c_1}{k+1}k_1 + \frac{1+c_2}{k}k_2 + \dots + \frac{1+c_n}{k}k_n = 1,\] which can be simplified into \[c_1\frac{k_1}{k+1} + \sum_{j = 2}^nc_j\frac{k_j}{k} = \frac{k_1}{k(k+1)}\] by equation 5.20.

We now prove that if equation 5.30 has a solution of \(c_i\)’s, all the elements of \(S\) would be \(o(\Delta t)\). To the case of \(S_1\) of equation 5.22, when \(k_1\) is reduced, some of \(k_j\)’s would be increased, due to \(\sum_ik_i' = \sum_{i}k_i = k\). Then equation 5.29 is changed into \[\frac{1+c_1}{k+1}k_1' + \frac{1+c_2}{k}k_2' + \dots + \frac{1+c_n}{k}k_n' > 1,\] because inequality 5.27 leads to \((1+c_j)/k > (1+c_1)/(k+1)\); thus, each element of \(S_1\) is \(o(\Delta t)\).

For \(S_2\), let \(\gamma = \min\{w_i: 1\le i\le n\} = \Delta t^{\frac{1+c_1}{k+1}}\). So \[w_1^{k_1'}w_2^{k_2'}\dots w_n^{k_n'} \ge \gamma^{k+1} = \Delta t^{1+c_1},\] which is also \(o(\Delta t)\). The elements of \(S_3\) can be dealt with similarly to \(S_2\) and the difference only lies in the value of \(\gamma\).

Therefore, the key point is the existence of \(c_i\)’s of equation 5.30 under the constraint of inequality 5.27. We can make \(c_1\) as small as possible until all the \(c_j\)’s are greater than it; noting that the extreme case is that \(c_1 = 0\) and \(c_j >0\), before reaching this limit, there’s a large space to tune them to reach a solution. For example, we can first find an initial solution with \(c_i>0\) for all \(i\) but cannot satisfy inequality 5.27 for all \(j\)’s. Choose one of \(c_j\)’s smaller than \(c_1\), say \(c_{\mu}\). Decrease \(c_1\) and simultaneously increase \(c_{\mu}\) to compensate the reduction of \(c_1k_1/(k+1)\) until \(c_{\mu} > c_1\). Each \(c_j\) that doesn’t fulfil inequality 5.27 should be dealt with analogously to \(c_{\mu}\) and a solution of \(c_i\)’s can be finally obtained. This completes the proof.

Lemma 6 and theorem 9 imply this conclusion.

6 Multivariate Global Approximation↩︎

This section turns to global approximation for multivariate functions, also in terms of splines as the univariate case of section 3. Section 6.1 constructs a multivariate spline approximating a given function to be realized by two-layer neural networks. Section 6.2 gives a fundamental result, a recurrence formula of polynomial pieces of a multivariate spline analogous to equation 3.4. Sections 6.3 deals with the zero-error part of a unit for higher-dimensional input to control approximation error.

Section 6.4 presents a property of solutions that is useful in locally determining the output weight of units. Section 6.5 proposes a most distinguished character of two-layer neural networks—smooth-continuity restriction. Section 6.6 discusses spline construction over a single strict partial order, the simplest multivariate solution leading to more complex ones. Section 6.7 proves universal approximation. Section 6.8 investigates the two-sided solutions that are ubiquitous in engineering. Section 6.9 is a short discussion about generalized tanh units.

6.1 Multivariate Splines↩︎

Two main techniques had been developed for multivariate splines. One is by tensor product de?, Boor2001? and the other is by a relationship between adjacent polynomial pieces [27]. In order to implement a spline via neural networks, we propose a new method based on the integral of linear hyperplanes, with an advantage that an approximation to a given function can be simultaneously obtained.

In the rest of the paper, some definitions are not formally presented because they can be found in [28], such as regions, adjacent regions, negative units, one-sided bases and boundaries. But for a self-contained purpose, they would be shortly described when first used.

A region is a connected part of \(\mathbb{R}^n\) divided by a set of \(n-1\)-dimensional hyperplane, including its boundary. A knot \(l\) in \(\mathbb{R}^n\) is an \(n-1\)-dimensional hyperplane \(\boldsymbol{w}^T\boldsymbol{x} + b = 0\) and can be regarded as being from the activation function \(\sigma(\boldsymbol{w}^T\boldsymbol{x} + b )\) of a unit \(\mathcal{U}\). Denote by \(l^+\) the region whose output of \(\mathcal{U}\) is positive, and by \(l^0\) the region \(\mathbb{R}^n-l^+\), both of which should include their boundary via a closure operation of sets.

Two regions of \(R_1\) and \(R_2\) are said to be adjacent if \(R_1 \cap R_2\) is \(n-1\)-dimensional as a part of the \(n-1\)-dimensional hyperplane separating them.

A multivariate smooth function was defined in definition 9 and the following proposition gives an equivalent definition that is more suitable for use in this paper.

Let \[D^{(k, \boldsymbol{d})} f(\boldsymbol{x}_0)= \sum_{\{\alpha_i\} \in A} \frac{\partial^{k} f(\boldsymbol{x}_0)}{\partial x_1^{\alpha_1}\partial x_2^{\alpha_2}\dots\partial x_n^{\alpha_n}}d_1^{\alpha_1}d_2^{\alpha_2} \dots d_n^{\alpha_n} = D^k f(\boldsymbol{x}_0+\boldsymbol{d}t) \Big |_{t=0},\] be the \(k\)th derivative of \(f(t)\) of equation 6.3 with respect to direction \(\boldsymbol{d}\), where \(k = 0, 1, \dots, m\), \(\boldsymbol{d} = [d_1, d_2, \dots, d_n]^T\), \(0\le \alpha_i\le k\) for \(i = 1, 2, \dots, n\) and each element of \(A\) is a set of \(\alpha_i\)’s satisfying \(\alpha_1+\alpha_2+\dots+\alpha_n=k\). By equation 6.4, it is obvious that definition \(9\) implies equation 6.3. Conversely, for example, to dimension \(x_1\), letting \(\boldsymbol{d}=[1,0,0,\dots,0]^T\), then equation 6.3 results in the continuity of \(\partial^{k} f(\boldsymbol{x})/\partial x_1^{\alpha_1}\partial x_2^{\alpha_2}\dots\partial x_n^{\alpha_n}\) with respect to variable \(x_1\); and other variables are similar. Thus, if and only if the condition of definition 9 is fulfilled, equation 6.3 holds:

By equation 6.4, write \[D^k(f(t), \boldsymbol{d}, \boldsymbol{x}_0) = \sum_{\{\alpha_i\} \in A} g_{\nu}(\boldsymbol{x})d_1^{\alpha_1}d_2^{\alpha_2} \dots d_n^{\alpha_n}\Big |_{\boldsymbol{x}=\boldsymbol{x}_0+\boldsymbol{d}t},\] where \[g_{\nu}(\boldsymbol{x}) = \frac{\partial^{k} f(\boldsymbol{x})}{\partial x_1^{\alpha_1}\partial x_2^{\alpha_2}\dots\partial x_n^{\alpha_n}}\] for \(\nu = 1, 2, \dots, |A|\). The hypersurface \(\mathcal{F}^{(k, f, \boldsymbol{d})}(\boldsymbol{x})\) is made up of all the values of \(D^k(f(t), \boldsymbol{d}, \boldsymbol{x}_0)\) for different \(t\) and \(\boldsymbol{x}_0\), or \[\mathcal{F}^{(k, f, \boldsymbol{d})}(\boldsymbol{x}) = \sum_{\{\alpha_i\} \in A} g_{\nu}(\boldsymbol{x})d_1^{\alpha_1}d_2^{\alpha_2} \dots d_n^{\alpha_n}.\] Since \(f(\boldsymbol{x}) \in C^m([0, 1]^n)\), \(g_{\nu}(\boldsymbol{x})\)’s is smooth with orders \(m-k\) and the conclusion follows by equation 6.7.

This is a special case of equation 6.7 when \(f(\boldsymbol{x}) = p(\boldsymbol{x})\) is a multivariate polynomial such that \(g_{\nu}(\boldsymbol{x})\) is a linear function of the entries of \(\boldsymbol{x} = [x_1, x_2, \dots, x_n]^n\).

We use equation 6.7 to express the condition of this lemma as \[\sum_{\{\alpha_i\} \in A} g^{(1)}_{\nu}(\boldsymbol{x})d_1^{\alpha_1}d_2^{\alpha_2} \dots d_n^{\alpha_n} = \sum_{\{\alpha_i\} \in A} g^{(2)}_{\nu}(\boldsymbol{x})d_1^{\alpha_1}d_2^{\alpha_2} \dots d_n^{\alpha_n}\] for \(\boldsymbol{x} \in l\), where \(g^{(i)}_{\nu}(\boldsymbol{x}) = \partial^{k} p_i(\boldsymbol{x})/\partial x_1^{\alpha_1}\partial x_2^{\alpha_2}\dots\partial x_n^{\alpha_n}\) for \(i = 1, 2\) is also a polynomial. Because equation 6.8 holds for infinitely many points on knot \(l\), we have \(g^{(1)}_{\nu}(\boldsymbol{x}) = g^{(2)}_{\nu}(\boldsymbol{x})\), regardless the entries \(d_i\)’s of \(\boldsymbol{d}\), which implies this lemma.

The proof is based on the univariate case of theorem 4. By lemma 7, the directional-derivative hypersurface \(\mathcal{F}^{(m-1, f, \boldsymbol{d})}(\boldsymbol{x})\) of equation 6.7 with respect to any direction \(\boldsymbol{d}\) is smooth. Then a continuous linear spline \(\mathscr{S}(\boldsymbol{x}) \in \mathfrak{S}_n^1(H, \mathcal{R})\) approximating \(\mathcal{F}^{(m-1, f, \boldsymbol{d})}(\boldsymbol{x})\) can be constructed (lemma 6 of [28]) and the approximation error could be arbitrarily small, provided that the volume of each region of \(\mathcal{R}\) is small enough; we use the formula \[\mathscr{S}(\boldsymbol{x}) \approx \mathcal{F}^{(m-1, f, \boldsymbol{d})}(\boldsymbol{x})\] to represent this meaning.

We can find an \(n-1\)-dimensional hyperplane \(l\) with \(U \subset l^+\). To arbitrary point \(\boldsymbol{x}_0 \in l\), the ray \(\boldsymbol{x}_0 + t\boldsymbol{d}\) can run through \(l^+\) as \(\boldsymbol{x}_0\) and \(t > 0\) both arbitrarily change, that is, \[l^+ = \{\boldsymbol{x}_0 + t\boldsymbol{d}: \boldsymbol{x}_0 \in l, t>0\},\] which contains \(U\) as a subset. To a fixed \(\boldsymbol{x}_0\), write \[S_t := \{t: \big(U \cap (\boldsymbol{x}_0 + t\boldsymbol{d})\big) \ne \emptyset\},\] whose minimum and maximum elements are denoted by \(t_i\) and \(t_a\), respectively. Let \(\mathcal{L} = U \cap (\boldsymbol{x}_0 + t\boldsymbol{d})\) , which is a line segment; \(\mathcal{L}\) passes though some regions of \(U\) and a univariate linear spline \(g(t)\) defined on \([t_i, t_a]\) can be derived from \(\mathscr{S}(\boldsymbol{x})\) of equation 6.9, whose values are approximately the \(m-1\)th directional derivative of \(f(\boldsymbol{x})\) with respect to \(\boldsymbol{d}\).

Then we use a method analogous to theorem 4 to integrate \(g(t)\) to obtain a smooth spline \(G(t)\) approximating \(f(\boldsymbol{x})\) on \(\mathcal{L}\); and as \(\boldsymbol{x}_0\) changes in \(l\), \(G(t)\) can form a multivariate smooth spline \(s(\boldsymbol{x}) \in \mathfrak{S}_n^m(H; \mathcal{R})\) approximating \(f(\boldsymbol{x})\). To this purpose, we first need a result that an integral of a directional-derivative hyperplane could be an algebraic hypersurface or a polynomial. We show an example of one region. On a certain region \(r\) of \(\mathcal{R}\), \(g(t)\) could form an \(n\)-dimensional hyperplane \(h(\boldsymbol{x})\) as \(\boldsymbol{x}_0\) varies within the associated domain and is the linear function on \(r\) approximating the directional-derivative of \(f(\boldsymbol{x})\). Let \(\mathscr{L}\) be one of the boundaries of \(r\), on which the function values and the directional-derivative values of \(f(\boldsymbol{x})\) or their approximations via polynomials are given; the two types of values are collectively called “boundary condition” in this proof. We want to obtain the expression of the integral of \(h(\boldsymbol{x})\) along direction \(\boldsymbol{d}\) under a boundary condition.

First see the case of \(m=2\) and \(n=2\) when \(h(\boldsymbol{x})\) is a plane and \(\mathscr{L}\) is a line segment. Let \(\boldsymbol{y}_0 \in \mathscr{L}\) and write \(\boldsymbol{y}_0=\boldsymbol{O}+t_1\boldsymbol{d}_1\), where \(\boldsymbol{O}\) is an arbitrary fixed point of \(\mathscr{L}\) and \(\boldsymbol{d}_1\) is the direction of \(\mathscr{L}\). Given a \(\boldsymbol{y}_0\), we have \[h(\boldsymbol{x})=h(t)=at+b(\boldsymbol{y}_0),\] where \(b(\boldsymbol{y}_0)\) is the function value of \(h(\boldsymbol{x})\) at \(\boldsymbol{y}_0\). We can write \(h(\boldsymbol{x})=at+b(t_1)=h(t,t_1)\), where \(b(t_1)\) is a linear function of \(t_1\). Thus, \(h(t,t_1)=h(\boldsymbol{x})\) is a linear function of \(t\) and \(t_1\) and \(\boldsymbol{t}=[t,t_1]^T\) together with origin point \(\boldsymbol{O}\) can be regarded as a new coordinate system derived from an affine transformation of \(\boldsymbol{x}\), namely \(\boldsymbol{t}=A\boldsymbol{x}+\boldsymbol{O}\). Then we have \[f(\boldsymbol{x}) = f(t, t_1) \approx p(t,t_1) = \int atdt + b(t_1)t + c(t_1),\] in which \[c(t_1)= \int a_1t_1dt_1 + b_1t_1 + c_1,\] where the parameters \(a_1\) and \(b_1\) are from \(h(\boldsymbol{O}+t_1\boldsymbol{d}_1)=a_1t_1+b_1\) and the constant \(c_1=f(\boldsymbol{O})\) or \(c_1 \approx f(\boldsymbol{O})\) when \(c_1\) is of a polynomial approximation to \(f(\boldsymbol{x})\); thus, \(p(t, t_1)\) is a polynomial of variables \(t\) and \(t_1\). The algebraic expression \(p(\boldsymbol{x})\) with respect to variable \(\boldsymbol{x}\) can be obtained by the affine transformation from \(\boldsymbol{x}\) to \(\boldsymbol{t}\) as mentioned above.

The case of arbitrary \(n \ge 3\) when \(m=2\) is only different in the expressions \(b(t_1)\) and \(c(t_1)\) and they would be changed into \(b(\boldsymbol{t}_r)\) and \(c(\boldsymbol{t}_r)\), respectively, where \(\boldsymbol{t}_r=[t_1, t_2, \dots, t_{n-1}]^T\). The term \(b(\boldsymbol{t}_r)\) is a linear function of the variables in \(\boldsymbol{t}_r\), while \(c(\boldsymbol{t}_r)\) is a polynomial of degree 2 that can be obtained similarly to equation 6.10. The polynomial \(p(\boldsymbol{x})\) can be constructed iteratively on the base of \(n=2\) of equation 6.10. For instance, to \(n=3\), we first use equation 6.10 to construct \(c(\boldsymbol{t}_r)\), then integrate each \(g(t)\) to form \(p(\boldsymbol{t})\). Analogously, when \(n=4\), \(c(\boldsymbol{t}_r)\) is obtained by the method of \(n=3\), after which the integral of \(g(t)\) is operated on the remaining dimension. This process can be repeatedly done until the desired dimensionality is reached.

To \(m \ge 3\), the difference lies in the number of integral operations. Based on the discussion above, iterative integrals as equation 6.10 always yield a polynomial, since the integrand is a polynomial of the previous step and the boundary condition newly introduced is also of a polynomial. Note that to a new integral or when \(m\) is increased by one, a derivative value or function value of \(f(\boldsymbol{x})\) should be used to determine the constant term of the associated polynomial, if necessary.

We now go back to the case of \(m=2\) and \(n=2\) to complete the spline construction. To produce a smooth spline \(s(\boldsymbol{x}) \in \mathfrak{S}_2^2(H; \mathcal{R})\) on all the regions of \(U\), use theorem 4 to integrate \(g(t)\)’s to yield smooth \(G(t)\)’s for all \(\boldsymbol{x}_0 \in l\); to ensure the smooth continuity in the whole input space, the constant term of the first polynomial of each \(G(t)\) (as \(c_1\) of equation 6.11) should be determined by a polynomial approximating \(f(\boldsymbol{x})\), except for one case that an initial value must be introduced from \(f(\boldsymbol{x})\).

We can also regard the construction process as being composed of the steps, with each corresponding to a certain \(h(\boldsymbol{x})\) that leads to different \(h(t)\) and \(c(t_1)\). In this view, the expression of each polynomial can be simultaneously obtained. The function value of \(f(\boldsymbol{x})\) is used only once for determining the constant term of a polynomial, and other cases are through the polynomials approximating \(f(\boldsymbol{x})\).

Generally, a spline \(s(\boldsymbol{x})\) for \(n \ge 3\) or \(m \ge 3\) can also be iteratively constructed on the basis of \(n=2\) and \(m=2\), as the polynomial construction above. The constructed \(s(\boldsymbol{x})\) is smooth with degree \(m-1\) with respect to direction \(\boldsymbol{d}\); and by lemma 9, \(s(\boldsymbol{x})\) is also smooth at any other direction and so \(s(\boldsymbol{x}) \in \mathfrak{S}_n^m(H; \mathcal{R})\).

To the approximation error, first note that by the construction method of \(\mathscr{S}(\boldsymbol{x})\), \(\varepsilon_1 = \max_{\boldsymbol{x}\in U}{|\mathscr{S}(\boldsymbol{x})-\mathcal{F}^{(m-1, f, \boldsymbol{d})}(\boldsymbol{x})|}\) can be arbitrarily small, and so is \(\varepsilon_2=\max_{\boldsymbol{x}\in U}|s(\boldsymbol{x})-f(\boldsymbol{x})|\) by the proof of theorem 4. The final approximation error \(\varepsilon=(\int_U(s(\boldsymbol{x})-f(\boldsymbol{x}))^2d\boldsymbol{x})^{1/2} < \varepsilon_2\). This completes the proof.

6.2 Relations between Adjacent Polynomial Pieces↩︎

This conclusion is by lemma 8 and definition 12.

We first prove that equation 6.12 implies the smoothness between \(p_1(\boldsymbol{x})\) and \(p_2(\boldsymbol{x})\). By equation 6.12, each line \(l = \boldsymbol{x}_0 + t\boldsymbol{d}\) penetrating \(\mathcal{L}\) yields two univariate polynomials \(p_1(l)\) and \(p_2(l)\) satisfying \(p_2(l) = p_1(l) + \lambda(\boldsymbol{w}^Tl + b)_+^m\). If \(\boldsymbol{x}_0 = l \cap \mathcal{L}\), the previous formula can be simplified into \(p_2(t) = p_1(t) + \lambda't^m_+\), where \(\lambda' = \lambda \boldsymbol{w}^T\boldsymbol{d}\), which means that \(p_2(t)\) and \(p_2(t)\) are smoothly continuous at \(\boldsymbol{x}_0\) by order \(m-1\), regardless of direction \(\boldsymbol{d}\). Note that different direction \(\boldsymbol{d}\) would result in distinct \(\lambda'\), but the set of all \(\boldsymbol{d}\) shares the same parameter \(\lambda\); this fact would be helpful to understand the converse conclusion of this theorem.

The converse is also true. By [28]’s theorem 2, if \(\mathscr{P}_2(\boldsymbol{x})\) and \(\mathscr{P}_1(\boldsymbol{x})\) are two linear functions and continuous at \(\mathcal{L}\), then \[\mathscr{P}_2(\boldsymbol{x}) = \mathscr{P}_1(\boldsymbol{x}) + \alpha(\boldsymbol{w}^T\boldsymbol{x} + b)_+\] for \(\boldsymbol{x} \in \mathcal{L}^+\), where \(\alpha\) is a constant. To equation 6.12, let \(q_1(\boldsymbol{x})\) and \(q_2(\boldsymbol{x})\) be the \(m-1\)th directional-derivative hyperplanes of \(p_1(\boldsymbol{x})\) and \(p_2(\boldsymbol{x})\) for some direction \(\boldsymbol{d}\), respectively. By lemmas 8 and 10, both \(q_1(\boldsymbol{x})\) and \(q_2(\boldsymbol{x})\) are \(n-1\)-dimensional hyperplanes and are continuous at \(\mathcal{L}\). Then by equation 6.13, we have \[q_2(\boldsymbol{x}) = q_1(\boldsymbol{x}) + \beta(\boldsymbol{w}^T\boldsymbol{x} + b)_+\] for \(\boldsymbol{x} \in \mathcal{L}^+\). When \(\boldsymbol{x} = \boldsymbol{x}_0 + t\boldsymbol{d}\) and \(\boldsymbol{x}_0 \in \mathcal{L}\), equation 6.14 becomes \[q_2(t) = q_1(t) + \gamma t_+,\] where \(\gamma=\beta \boldsymbol{w}^T\boldsymbol{d}\), in which \(q_2(t)\) and \(q_1(t)\) comprise a continuous linear spline.

Then by the method of theorem 4, integrate equation 6.15 repeatedly by \(m-1\) times and obtain \[p_2(\boldsymbol{x})=p_2(t) = p_1(t) + \frac{\gamma}{m!}t^{m}_+ = p_1(\boldsymbol{x}) + \frac{\beta}{m!(\boldsymbol{w}^T\boldsymbol{d})^{m-1}}(\boldsymbol{w}^T\boldsymbol{x}+b)^{m}_+\] for \(\boldsymbol{x} = \boldsymbol{x}_0 + t\boldsymbol{d}\), which is equivalent to equation 6.12 on line \(l=\boldsymbol{x}_0 + t\boldsymbol{d}\), with \[\lambda = \frac{\beta}{m!(\boldsymbol{w}^T\boldsymbol{d})^{m-1}}.\]

Note that equation 6.16 holds for arbitrary \(\boldsymbol{x}_0 \in \mathcal{L}^+\) and thus we have proved equation 6.12 for a certain direction \(\boldsymbol{d}\). Next step is to prove that the obtained parameter \(\lambda\) is invariant with respect to \(\boldsymbol{d}\). The parameter \(\lambda\) is unique, because different \(\lambda\) leads to different \(p_2(\boldsymbol{x})\) and it’s an one-to-one correspondence when \(\boldsymbol{w}^T\boldsymbol{x} + b\) is fixed. This completes the proof.

6.3 Zero-Part Error of Units↩︎

The proof is by reducing the multivariate input to the univariate case of lemma 3. Given a line \(h = \boldsymbol{x}_0 + t\boldsymbol{d}\) with \(\boldsymbol{x}_0 \in l_0\), the activation function \(\phi(\boldsymbol{x})\) is reduced to \[\phi(t) = \sigma(\boldsymbol{w}^T(\boldsymbol{x}_0 + t\boldsymbol{d}) + b) = \sigma(wt),\] where \(w = \boldsymbol{w}^T\boldsymbol{d}\), which can be regarded as a univariate activation function of lemma 3. Correspondingly, \(h \cap \mathcal{L}\) is a point knot \[t_k = (b-b_k)/w\] obtained by \(\boldsymbol{w}^T(\boldsymbol{x}_0+t\boldsymbol{d})+b_k = 0\), as the counterpart \(x_k\) of lemma 3. Similarly, \(t_{k+1} = (b-b_{k+1})/w\) determines the point \(h \cap l_{k+1}\). Then we use the method of lemma 3 to change \(\sigma(wt)\) into \[\phi_{a}(t) = \sigma(\rho wt + \gamma),\] such that \[\lim_{\rho(\gamma) \to +\infty}\phi_{a}(t) = 0\] for \(t \le t_k\), \[\phi_{a}(t) \approx c'(t-t_k)_+^m\] for \(t \in (t_k, t_{k+1}]\), and \[\lim_{\rho(\gamma) \to +\infty}c' = +\infty.\]

When \(\boldsymbol{x} = \boldsymbol{x}_0 + t\boldsymbol{d}\), equations 6.24 and 6.25 can be expressed as equations 6.18 and 6.19, respectively. Since \(\boldsymbol{w}^T(\boldsymbol{x}_0 + t\boldsymbol{d}) + b_k = w(t-t_k)\), equation 6.26 is equivalent to \[\phi_{a}(\boldsymbol{x}) \approx c(\boldsymbol{w}^T\boldsymbol{x} + b_k)_+^m\] for \(\boldsymbol{x} = \boldsymbol{x}_0 + t\boldsymbol{d}\) with \(t \in (t_k, t_{k+1}]\), where \[c = c' / w;\] thus, \[\lim_{\rho(\gamma) \to +\infty}c = +\infty\] by equation 6.27. So we proved this lemma over the line \(\boldsymbol{x} = \boldsymbol{x}_0 + t\boldsymbol{d}\).

Note that equations from 6.24 to 6.27 as well as equation 6.29 are all independent of the selection of \(\boldsymbol{x}_0\). When \(\boldsymbol{x}_0\) runs through \(l_0\) and \(t\) arbitrarily changes, \(p_{k}(\boldsymbol{x}) = c(\boldsymbol{w}^T\boldsymbol{x} + b_k)_+^m\) of equation 6.28 forms an \(n\)-variate polynomial of degree \(m\), which is composed of identical \(q_{k}(t) = c'(t-t_k)^m_+\) of equation 6.26 at different \(\boldsymbol{x}_0\). The formed \(p_{k}(\boldsymbol{x})\) can approximate \(\phi_{a}(\boldsymbol{x})\) as precisely as possible on \(R_2\) according to equation 6.26, and the approximation-error analysis is similar to the proof of theorem 11. To arbitrary another direction \(\boldsymbol{d}'\), the above process can also yield a multivariate polynomial \(p_{k}'(\boldsymbol{x})\) and it would be equal to \(p_{k}(\boldsymbol{x})\), due to the uniqueness of the polynomial approximating \(\phi_{a}(\boldsymbol{x})\) on \(R_2\) under arbitrary precision; so the parameter \(c\) is unique with respect to different \(\boldsymbol{d}\). This completes the proof.

The proof is similar to that of theorem 5.

6.4 Multiple Expressions of Polynomial Pieces↩︎

Hereafter, to simplify the descriptions, given a partition of \(U = [0, 1]^n\) by a set \(H = \{l_i : 1\le i\le\Theta\}\) of \(n-1\)-dimensional hyperplanes, if we say that the units of the hidden layer of a two-layer neural network \(\mathfrak{N}\) with generalized sigmoidal units are derived from \(H\), it means that to each local unit \(u_{\nu}\), \(l_{\nu}\) is the zero-error hyperplane of \(u_{\nu}\); and based on the results of sections 6.3, the activation function of \(u_{\nu}\) on \(l_{\nu}^0\) can be regarded as zero. To global unit \(u_{\mu}\), \(l_{\mu}\) is the original meaning.

The following results not only introduce a parameter-setting technique, but also involve a local property of parameters that is useful in understanding training solutions. In this section, “\(f(\boldsymbol{x})\approx g(\boldsymbol{x})\)” means that \(f(\boldsymbol{x})\) can approximate \(g(\boldsymbol{x})\) with arbitrary precision.

We first see an example of Figure 4. To the polynomials on region \(R_3\), we have \[p_{j3}(\boldsymbol{x}) = p_{j2}(\boldsymbol{x}) + \alpha_{j3}(\boldsymbol{w}_2^T\boldsymbol{x} + b_2)^m_+\] for \(j = 1, 2, 3\). Since \(R_2 \subset l_1^+ \cap l_3^+\), \[p_2(\boldsymbol{x}) = \lambda_1p_{12}(\boldsymbol{x}) + \lambda_3p_{32}(\boldsymbol{x}).\] Similarly, \[p_3(\boldsymbol{x}) = \lambda_1p_{13}(\boldsymbol{x}) + \lambda_3p_{33}(\boldsymbol{x}) + \lambda_2p_{23}(\boldsymbol{x})\] due to \(R_3 \subset l_1^+ \cap l_3^+ \cap l_2^+\). Because \(R_2 \subset l_2^0\) and \(R_3 \subset l_2^+\), in equation 6.36, \[p_{23}(\boldsymbol{x}) = \alpha_{23}(\boldsymbol{w}_2^T\boldsymbol{x}+b_2)^m_+ \approx \phi_2(\boldsymbol{x})\] since \(p_{22}(\boldsymbol{x})\approx0\). Equations from 6.36 to 6.39 yield \[\begin{align} p_3(\boldsymbol{x}) &= p_2(\boldsymbol{x}) + (\lambda_1\alpha_{13}+\lambda_3\alpha_{33}+\lambda_2\alpha_{23})(\boldsymbol{w}_2^T\boldsymbol{x} + b_2)^m_+\\ &\approx p_2(\boldsymbol{x}) + (\lambda_1\alpha_{13}+\lambda_3\alpha_{33})(\boldsymbol{w}_2^T\boldsymbol{x} + b_2)^m_+ + \lambda_2\phi_2(\boldsymbol{x}). \end{align}\]

The case of equation 6.34 is similar. This example contains the general principle of this lemma.

This theorem is a more general description of lemma 12.

6.5 Smooth-Continuity Restriction↩︎

Due to lemma 8 and \(s(\boldsymbol{x}) \in C^{m-1}([0,1]^n)\), the conclusion follows.

Denote by \(\phi_i(\boldsymbol{x}) \in C^{m}([0, 1]^n)\) the activation function of the \(i\)th hidden-layer unit \(u_i\) of \(\mathfrak{N}\). By theorem 11, each \(\phi_i(\boldsymbol{x})\) can be approximated by a spline \(s_i(\boldsymbol{x}) \;\in \mathfrak{S}_n^{m}(H, \mathcal{R})\) to any desired accuracy. Let \(g(\boldsymbol{x}) = \sum_{i = 1}^{\Theta}\lambda_i\phi_i(\boldsymbol{x})\) be the function output by \(\mathfrak{N}\), which can be approximated by \[\mathcal{S}(\boldsymbol{x}) = \sum_{i=1}^{\Theta}\lambda_is_i(\boldsymbol{x})\] as precisely as possible, with \(\mathcal{S}(\boldsymbol{x}) \in \mathfrak{S}_n^{m}(H, \mathcal{R})\) as well. Write \(\hat{s}_i(\boldsymbol{x})\), which belongs to \(\mathfrak{S}_n^1(H, \mathcal{R})\) and is composed of the directional-derivative hyperplanes of \(s_i(\boldsymbol{x})\) with respect to some direction \(\boldsymbol{d}\). Letting \[\hat{\mathcal{S}}(\boldsymbol{x}) = \sum_{i=1}^{\Theta}\lambda_i\hat{s_i}(\boldsymbol{x}),\] we have \(\hat{\mathcal{S}}(\boldsymbol{x}) \in \mathfrak{S}_n^1(H, \mathcal{R})\). Since the integral of the right side of equation 6.47 with respect to \(\boldsymbol{d}\) (as in theorem 12) gives \(\mathcal{S}(\boldsymbol{x})\) of equation 6.46, the directional-derivative hyperplanes of \(\mathcal{S}(\boldsymbol{x})\) (lemma 13) comprise \(\hat{\mathcal{S}}(\boldsymbol{x})\) of equation 6.47.

Let \(q_{i\nu}(\boldsymbol{x})\) be the linear function on \(R_{\nu}\) of \(\hat{s_i}(\boldsymbol{x})\), and \(q_{\nu}(\boldsymbol{x})\) be the one of \(\hat{\mathcal{S}}(\boldsymbol{x})\). According to the “boundary-determination principle” of splines in \(\mathfrak{S}_n^1(H, \mathcal{R})\) (see the proof of [28]’s lemma 6), \(q_{i2}(\boldsymbol{x})\) and \(q_{i4}(\boldsymbol{x})\) for \(\hat{s_i}(\boldsymbol{x})\) can uniquely determine \(q_{i3}(\boldsymbol{x})\) due to the property of continuous piecewise linear functions, and so is the unique determination of \(q_{3}(\boldsymbol{x})\) through \(q_2(\boldsymbol{x})\) and \(q_4(\boldsymbol{x})\) for \(\hat{\mathcal{S}}(\boldsymbol{x})\). Thus, if the linear combination of equation 6.47 via weight parameters \(\lambda_i\)’s are set to produce \(q_2(\boldsymbol{x})\) and \(q_4(\boldsymbol{x})\), \(q_3(\boldsymbol{x})\) is simultaneously implemented. This principle remains invariant under the integral operation of the two-sides of equation 6.47 and this lemma follows.

The boundary of a region \(R\) is a set whose each element is part of a hyperplane forming \(R\). For example, in Figure 4, the boundary of \(R_3\) is \(\{l_1\cap R_3, l_2\cap R_3\}\) with two elements.

This conclusion is due to the proof of lemma 14. To the example of Figure 4, this theorem means that if the polynomials on the boundary of \(R_3\) formed by \(l_1\) and \(l_2\) have been realized by \(\mathfrak{N}\), the one on \(R_3\) is simultaneously implemented.

6.6 Single Strict Partial Order↩︎

Let \(H = \{l_1, l_2, \dots, l_\zeta\}\) be a set of \(n-1\)-dimensional hyperplane satisfying \[l_1 \prec l_2 \prec \dots \prec l_{\zeta},\] where “\(\prec\)” represents a strict partial order [28] of higher-dimensional knots analogous to the “\(<\)” (less than) relation of real numbers. The region \(R_i\)’s for \(i = 1, 2, \dots, \zeta\) are the corresponding ordered regions and form a set \[\mathcal{R}:=\{R_i: 1\le i\le \zeta\}.\] To each \(s(\boldsymbol{x}) \in \mathfrak{S}^m_n(H, \mathcal{R})\) over a single strict partial order, according to theorem 12. we get the recurrence formula \[s_{\nu}(\boldsymbol{x}) = s_{\nu-1}(\boldsymbol{x}) + \alpha_{\nu-1}(\boldsymbol{w}_{\nu-1}^T\boldsymbol{x}+b_{\nu-1})^m_+\] for \(\nu = 2, 3, \dots, \zeta-1\) analogously to the univariate case of equation 3.4, where \(s_{\nu}(\boldsymbol{x})\) is the polynomial on \(R_{\nu}\). Equation 6.50 is a foundation of solution construction for higher-dimensional input.

The definition of one-sided bases of \(\mathfrak{S}^m_n(H, \mathcal{R})\) is similar to that of \(\mathfrak{S}^1_n(H, \mathcal{R})\) for linear splines ([28]’s definition 10). Let \[\mathfrak{B}_1 = \{\rho_{j}(\boldsymbol{x}) = \sigma(\boldsymbol{w}_{j}^T\boldsymbol{x} + b_{j}): -\xi \le j \le \zeta, \;\binom{n+m}{m} \le \xi+1\}\] be a set of one-sided bases of \(\mathfrak{S}_n^m(H, \mathcal{R})\), in which \(\rho_{\nu}(\boldsymbol{x})\)’s for \(\nu = -\xi, -\xi+1, \dots, 0\) are global bases and the remaining ones are local bases. Correspondingly, the associated units \(u_j\)’s are also classified into global and local units that are consistent with definition 15.

On the basis of the strict partial order “\(\prec\)” of equation 6.48 that is equivalent to “\(<\)”, theorems 10, 12 and 13 convert this problem to the univariate case of theorem 6.

The solution of a spline \(\mathcal{S}(\boldsymbol{x}) \in \mathfrak{S}_n^m(H, \mathcal{R})\) over a single strict partial order, including the one of theorem 16 as well as the two-sided case, also has a general framework analogous to the univariate case of theorem 8. To a generalized sigmoidal unit \(\mathfrak{U}\) with activation function \(\phi(\boldsymbol{x}) = \sigma(\boldsymbol{w}^T\boldsymbol{x}+b) \in C^{m}([0, 1]^n)\), the indicator function \(\mathbb{I}: N \to \{0, 1\}\) for \(\mathfrak{U}\), where \(N = \{1, 2, \dots, \zeta\}\), is defined as \[\mathbb{I}(i) = \begin{cases} \;\;0, \;\text{if} \;\big(\int_{R_i}\phi(\boldsymbol{x})^2d\boldsymbol{x}\big)^{1/2} < \varepsilon \\ \;\;1, \;\;\;\;\;\;\;\;\;\;\;\text{otherwise} \end{cases}\] for \(i = 1, 2, \dots, \zeta\). To each local unit \(u_j\)’s of equation 6.51, let \(\mathbb{I}_{j}\) be its indicator function. As the counterpart of equation 4.20, write \[\mathscr{A} = \begin{pmat}({}) \mathcal{A} & & & \mathcal{B} \cr \mathcal{C} & & & \mathcal{D} \cr \end{pmat},\] in which \(\mathscr{A} = \mathcal{W}(u_i(\boldsymbol{x}_0))^T\) of equation 5.10, while \(\mathcal{B}\), \(\mathcal{C}\) and \(\mathcal{D}\) can be obtained similarly to those of equation 4.20, respectively. A two-layer neural network with generalized sigmoidal units have a solution of \(\mathcal{S}(\boldsymbol{x})\) if and only if the rank of \[\mathscr{A}\boldsymbol{\lambda} = \boldsymbol{b}\] is \[\gamma = \zeta-1+\binom{n+m}{m}.\]

6.7 Multivariate Universal Approximation↩︎

The concepts of “standard partition” and “universal global hyperplane” were defined in [28]. Intuitively, when the input is two-dimensional, a standard partition is a set of regions that appear to be formed by horizontal and vertical lines as in Figure 5; a universal global hyperplane \(l\) of \(U=[0,1]^n\) means that \(U \subset l^+\).

In combination with smooth-continuity restriction of theorem 15, the proof is analogous to [28]’s lemma 10.

The boundary of a partition can be found in [28]’s definition 23. As an example in Figure 5, the set of the red rectangles is the boundary of the partition of \(U\).

This corollary is a direct consequence of the proof of theorem 17.

The proof is formally the same as that of [28]’s theorem 10. The preceding results of section 6 provided all the prerequisites to prove this theorem.

6.8 Two-Sided Bases↩︎

To two-sided solutions over a single strict partial order, equation 6.55 of section 6.6 has provided a general principe; this section investigates the more complicated case for multiple strict partial orders. Given a unit \(\mathcal{U}\) with activation function \(\phi(\boldsymbol{x})=\sigma(\boldsymbol{w}^T\boldsymbol{x}+b)\), if \(\phi(\boldsymbol{x})\) is changed into \(\phi(\boldsymbol{x})=\sigma(-\boldsymbol{w}^T\boldsymbol{x}-b)\), we call the associated unit a negative form of \(\mathcal{U}\), denoted by \(-\mathcal{U}\). Correspondingly, the hyperplane \(-l\) with equation \(-\boldsymbol{w}^T\boldsymbol{x}-b=0\) is the negative form of \(l\) whose equation is \(\boldsymbol{w}^T\boldsymbol{x}+b=0\).

Before modification, \(\mathcal{S}(\boldsymbol{x})\) is constructed by theorem 17. The principle of this theorem can be illustrated by the example of Figure 5. In Figure 5, there are two strict partial orders \(\mathscr{P}_1 = l_{11} \prec l_{12}\) and \(\mathscr{P}_2 = l_{21} \prec l_{22}\); lines \(l_{ij}\) for \(i, j = 1, 2\) form a standard partition of \(U = [0, 1]^2\).

A typical feature of \(\mathscr{P}_1\) and \(\mathscr{P}_2\) is that any line of \(\mathscr{P}_1\) cannot influence the ordered regions of \(\mathscr{P}_2\) in terms of positive outputs, and vice verse; on the other hand, this feature implies that any negative line of \(\mathscr{P}_1\) (or \(\mathscr{P}_2\)) would contribute to a global influence on \(\mathscr{P}_2\) (or \(\mathscr{P}_1\)), such that \(\mathscr{P}_2\) (or \(\mathscr{P}_1\)) can reformulate its parameters to produce the original function without disturbing \(\mathscr{P}_1\) (or \(\mathscr{P}_2\)).

Thus, if we change one of the lines, \(l_{21}\), say, of \(\mathscr{P}_2\) into its negative form, the parameter readjusting of \(\mathscr{P}_1\) can compensate this disturbance and restore the original output, without influencing \(\mathscr{P}_2\) simultaneously. So if the altered units of \(\mathscr{P}_2\) forms a set of two-sided bases, the output of \(\mathscr{P}_2\) can also remain the same. By smooth-continuity restriction, a solution of \(\mathcal{S}(\boldsymbol{x})\) still exists. The principle underlying the above example is related to an order of multiple strict partial orders and can be founded in [28]’s theorem 7.

The influence (if any) of other hyperplanes on \(\mathcal{R}'\) can be regarded as that of global ones for \(\mathcal{R}'\); and thus when only considering \(\mathscr{P}\), the principle of two-sided bases must be satisfied.

6.9 Generalized Tanh-Unit Case↩︎

Under the principle similar to corollary 1, when the input-dimensionality \(n\ge2\), the preceding results of section 6 for two-layer neural network with generalized sigmoid units are all applicable to the generalized tanh-unit case, except that the number of the units required should be increased by at least 1.

7 Explanation of Training Solutions↩︎

This section uses experiments to demonstrate that our theory grasps the main principle of solutions of two-layer neural networks obtained by the back-propagation algorithm. The examples to be given can also be regarded as the phenomena predicted by the theory.

7.1 Preliminaries↩︎

Since the training process deals with discrete data points, we use the following method to check whether a unit is local or global, which can be regarded as a discrete version of definition 5. Suppose that \[\big\{\sum_{i}\big(f(x_i) - \sum_{j}\lambda_j\phi_j(x_i)\big)^2\big\}^{1/2} = \varepsilon.\] To each unit \(u_{\nu}\) whose activation function is \(\phi_{\nu}(x)=\sigma(w_{\nu}x+b_{\nu})\), when \(w_{\nu} > 0\), by truncation operation at \(0< x_k<1\), equation 7.1 becomes \[\big\{\sum_{i:x_i\le x_k}\big(f(x_i) - \sum_{j \ne \nu}\lambda_j\phi_j(x_i)\big)^2 + \sum_{i:x_i>x_k}\big(f(x_i) - \sum_{j}\lambda_j\phi_j(x_i)\big)^2\big\}^{1/2} = \varepsilon'.\] The zero-error point of \(u_{\nu}\) is determined by \[z_{\nu0} = \sup\{x_k: |\varepsilon' - \varepsilon| < \gamma_1\varepsilon\},\] where \(0 < \gamma_1 < 1\) is a threshold.

Corresponding to equation 3.57, we check if \[|\alpha_{\nu}\phi_{\nu}(0)| > \gamma_2\varepsilon,\] where \(0 < \gamma_2 < 1\), to exclude the case that as \(\varepsilon \to 0\), a unit reduces its truncation error by decreasing \(x_k\) until \([0, x_k]\) shrinks to a point. If inequality 7.4 is satisfied, it is not a local unit. Thus, if \(0 < z_{\nu0} < 1\) and \(|\alpha_{\nu}\phi_{\nu}(0)| \le \gamma_2\varepsilon\), \(u_{\nu}\) is a local unit; if \(z_{\nu0}\le 0\) or \(z_{\nu0}\ge 1\) or \(|\alpha_{\nu}\phi_{\nu}(0)|>\gamma_2\varepsilon\), it is a global unit.

If \(w_{\nu} < 0\), equations from 7.2 to 7.4 are modified to \[\big\{\sum_{i:x_i\ge x_k}\big(f(x_i) - \sum_{j \ne \nu}\lambda_j\phi_j(x_i)\big)^2 + \sum_{i:x_i< x_k}\big(f(x_i) - \sum_{j}\lambda_j\phi_j(x_i)\big)^2\big\}^{1/2} = \varepsilon',\] \[z_{\nu0} = \inf\{x_k: |\varepsilon' - \varepsilon| < \gamma_1\varepsilon\},\] and \[|\alpha_{\nu}\phi_{\nu}(1)| > \gamma_2\varepsilon,\] respectively. Equation 7.1 is a local approximation provided that all \(u_{j}\)’s are global unit; otherwise, it is a global approximation with local units.

We give a criterion to determine whether a unit is inactivated. Let \[\big\{\big(\sum_{i}\big(f(x_i)-\sum_{j \ne \nu}\lambda_j\phi_j(x_i)\big)^2\big\}^{1/2} = \varepsilon''.\] If \[|\varepsilon'' - \varepsilon| < \gamma_3\varepsilon,\] where \(0 < \gamma_3 < 1\), \(u_{\nu}\) is inactivated.

When the input-dimensionality \(n\ge2\), equation 7.1 becomes \[\big\{\sum_{i}\big(f(\boldsymbol{x}_i)-\sum_{j}\lambda_j\phi_j(\boldsymbol{x}_i)\big)^2\big\}^{1/2} = \varepsilon.\] Each unit \(u_j\) leads to an \(n-1\)-dimensional hyperplane \(\boldsymbol{w}_j^Tx + b_j = 0\). To the truncation operation of \(u_{\nu}\) on \(L_{\nu}^0\) whose equation is \(\boldsymbol{w}_{\nu}^Tx + b_{\nu}' = 0\), where \(b_{\nu}'<b_{\nu}\), the higher-dimensional counterpart of equation 7.2 is \[\big\{\sum_{i:\boldsymbol{x}_i \in L_{\nu}^0}\big(f(\boldsymbol{x}_i)-\sum_{j \ne \nu}\lambda_j\phi_j(\boldsymbol{x}_i)\big)^2 + \sum_{i:\boldsymbol{x}_i \in L_{\nu}^+}\big(f(\boldsymbol{x}_i)-\sum_{j}\lambda_j\phi_j(\boldsymbol{x}_i)\big)^2\big\}^{1/2} = \varepsilon'.\] The parameter \(b_{\nu0}\) of the zero-error hyperplane is determined by \[b_{0\nu} = \sup\{b_{\nu}': |\varepsilon' - \varepsilon| < \gamma_1\varepsilon\}.\]

To the condition of equation 6.32, we make an analogy with inequality 7.4 as follows. Let \[\mathcal{B}_{\nu} = \inf\{b_{\nu}': U \cap L_{\nu}^0 \ne \emptyset\}\] and \(\mathscr{L}_{\nu}\) be the hyperplane of equation \(\boldsymbol{w}_{\nu}^Tx + \mathcal{B}_{\nu} = 0\). Then use \[\sum_{x_i \in S}|\alpha_{\nu}\phi_{\nu}(x_i)| > \gamma_2\varepsilon,\] where \(S = U \cap \mathscr{L}_{\nu}\), to play the role as inequality 7.4.

The difference between \[\big\{\sum_{i}\big(f(\boldsymbol{x}_i) - \sum_{j \ne \nu}\alpha_j\phi_j(\boldsymbol{x}_i)\big)^2\big\}^{1/2} = \varepsilon''.\] and \(\varepsilon\) of equation 7.10 can determine whether a unit is inactivated as inequality 7.9.

7.2 One-Dimensional Input↩︎

In practice, equation 7.3 is used in terms of \(z_{\nu0} = \max\{x_k: |\varepsilon'-\varepsilon| < \gamma_1\varepsilon\}\) by discretizing \(x_k \in [0, 1]\) in step length 0.01, and equation 7.6 is similar. The activation function of hidden-layer units of a two-layer neural network \(\mathfrak{N}\) to be trained is the logistic sigmoid type \(\sigma(x) = 1 / (1 + e^{-x})\); \(\mathfrak{N}\) is sometimes called a logistical sigmoidal neural network in the rest of section 7.

In Figure 6a, the circle points are from the discretization of a continuous function \(y = x^3 + 3\) on \([0, 1]\) with step length \(0.01\), and are fitted by \(\mathfrak{N}\) with \(\Theta = 10\) units in the hidden layer. In order to show \(y(x)\) and the activation functions on the same graph, the former is normalized according to the latter by scaling its function values to have maximum 1. The green asterisks are produced by \(\mathfrak{N}\). The weights and biases of \(\mathfrak{N}\) are randomly initialized by the uniform distribution \(U(-1, 1)\) and are optimized by the back propagation algorithm with the learning rate \(c = 0.05\) and the number of training iterations \(N = 5000\). Set \(\gamma_1 = \gamma_2 = \gamma_3 = 0.01\).

If a unit is a local one, a black line segment will be placed at its zero-error point (such as the examples of Figure 6b) and the corresponding activation-function curve is depicted by red colors; otherwise it is a global unit represented by blue curves. The number nearby a line or curve is the index of the unit. In Figure 6d, the nearly identical units \(u_2\) and \(u_5\) are too close to \(u_3\) such that only \(u_3\) is labeled; the zero-error points of \(u_2\) and \(u_5\) are also almost the same. Except for the fitted functions and possibly different parameter settings, the above descriptions are applicable to the remaining subfigures.

Note that there exist almost identical activation functions but with different zero-error points; the reason is that their output weights are different, while the zero-error point is determined by both the activation function and its output weight. The following items are the explanation of experimental results.

• Local approximation. Since there’s no local unit in Figure 6a, it is a local approximation, and similarly for Figure 6c for function \(y = e^{2x}\).

• Global approximation. In Figure 6b the function is \(y = 32x^3 + 3\) and the parameter settings are the same as those of Figure 6a except for \(\gamma_3 = 0.05\). There are 6 global units and 3 local units, while \(u_9\) is inactivated. The function of Figure 6d is \(y = e^{8x}\) with parameters \(c = 0.0001\),\(N = 5000\), \(\gamma_1 = \gamma_2 = 0.0001\) and \(\gamma_3 = 0.01\). There’s no global unit and the function values on the leftmost part of \([0, 1]\) are approximately zero compared to the larger ones on the right. Figure 6d clearly demonstrates the pattern of global approximations.

• Two-sided bases. Figure 6e is of a two-sided solution that can be explained by proposition 2, in which \(y = 30(\sin15x + 1)\), \(\Theta = 10\), \(c = 0.001\), \(N = 10000\), \(\gamma_1 = 0.01, \gamma_2 = 0.001\), and \(\gamma_3 = 0.05\). There exist 7 local units, 2 global units and 1 inactivated unit. Compared with Figure 3 of proposition 2, Figure 6e has the similar pattern that positive unit \(u_1\) and negative unit \(u_4\) share the same zero-error point \(x = 0.32\). In Figure 6f, \(y = 30\sin(6x + 3) + 3\), \(c = 0.01\), \(N = 5000\) and other parameters are set as Figure 6e. We didn’t find a theoretical solution for the pattern of Figure 6f, but it can be included in the framework of theorem 8 and would have a solution provided the condition of theorem 8 is satisfied.

7.3 Two-Dimensional Input↩︎

Similarly to the one-dimensional case, by discretizing \(b_{\nu}'\), equation 7.12 becomes \(b_{0\nu} = \max\{b_{\nu}': |\varepsilon' - \varepsilon| < \gamma_1\varepsilon\}\). The circle points of Figure 7a are from the discretization of \(z = 16(x^3 + y^3)+3\) on \(U = [0, 1]^2\) with step length 0.1 in both \(x\) and \(y\) dimensions. The asterisks are the outputs of a two-layer logistical sigmoidal neural network \(\mathfrak{N}_2\) with \(\Theta = 20\) units in the hidden layer. The parameters of \(\mathfrak{N}_2\) are initialized by uniform distribution \(U(-1, 1)\) and are trained by the back propagation algorithm with learning rate \(c = 0.01\) and learning step \(N = 5000\). Set \(\gamma_1 = \gamma_2 = \gamma_3 = 0.01\).

Figure 7b is the solution of Figure 7a expressed by the zero-error lines (red ones in the figure) of local units; global units whose corresponding line is out of \([0,1]^2\) as well as inactivated units are not depicted in the figure. There are 6 local units, 10 global units and 4 inactivated units. The number of the global units is \(\binom{5}{3} = 10\), equal to the number of the coefficients of a 2-variate polynomial of degree 3; thus any initial polynomial piece of this type could be produced, if the rank of the associated generalized Wronskian matrix satisfies the condition of theorem 9.

We see that the partition of Figure 7b resembles a standard partition of theorem 17 (e.g., Figure 5); for simplicity, we also call Figure 7b a standard partition in this section. Notice that both sigmoidal and ReLU (Figure 11b of [28]) two-layer neural networks have a standard-partition solution, providing an experimental evidence that the two types of neural networks obey a similar rule related to continuity restriction, as indicated in this paper’s theorem 15 and [28]’s theorem 9.

To the two-sided bases, a strict partial order \(-l_{20} \prec l_{15} \prec l_{9} \prec l_{7}\) can be obtained by changing \(l_{20}\) into its negative form \(-l_{20}\); so the original four lines yield two-sided bases if the associated matrix \(\mathscr{A}\) of 6.55 satisfy the condition of equation 6.56. The change from \(-l_{20}\) to \(l_{20}\) does not influence the existence of solutions over the other strict partial order \(l_2 \prec l_{11}\), as explained in proposition 8.

In Figure 7c, the function is \(z = \sin{3(x + y + 1)} + 3\) and the parameter setting is the same as that of Figure 7a. Figure 7d is the solution of Figure 7c and can be nearly regarded as a standard partition. The numbers of local, global and inactivated units are 6, 7 and 7, respectively. We depict one global unit \(u_{19}\) (blue line) because its line lies in \(U=[0,1]^2\). The seven global units can form any linear function (the least number required is 4), provided the condition of theorem 9 is fulfilled; and they can also produce some polynomials of degree \(m\ge2\).

We use the example of Figures 7e and 7f to show how the principle of smooth continuity restriction of theorem 15 applies to non-standard partition case. The data points of Figure 7e are from Figure 7a and the parameters are also identical except for the number of units \(\Theta = 10\). The solution of Figure 7f includes 5 local units, 4 global units and 1 inactivated unit and is a combination of a standard partition \(\mathcal{P}\) on \(U \cap l_3^+\) with five regions in \(U \cap l_3^0\). The polynomials on \(\mathcal{P}\) can be obtained by the method of Figure 7b or 7d; after that, suppose that the output weight of \(u_3\) is set to yield the polynomial on \(l_9^+l_3^0l_4^0 \cap U\). Then the polynomial on \(l_4^+l_3^0l_2^+ \cap U\) is determined by smooth-continuity restriction, because all the polynomials on the adjacent regions have already been previously set; the remaining three regions can be similarly dealt with.

7.4 Tanh-Unit Case↩︎

To the solution of two-layer neural networks composed of generalized tanh units, corollary 1 indicates that there exist units exclusively for producing constant term \(\mathfrak{C}\) of equation 3.84. We now use experiments to verify that.

In Figure 8a, similarly to the figures of section 7.2, the blue circles are from function \(y = \sin{20x}\); the function values shown in the figure are normalized into \([0, 1]\) for visualization convenience. The hidden-layer units of network \(\mathfrak{N}\) are of the tanh type whose activation function is \[\sigma(x)=(e^{x}-e^{-x})/(e^{x}+e^{-x});\] the parameter settings of \(\mathfrak{N}\) include unit number \(\Theta = 20\), learning rate \(c = 0.05\) and learning steps \(N = 5000\); the green asterisks are produced by \(\mathfrak{N}\). Use the method of section 7.1 to find out the inactivated units and do not depict them.

To judge whether a unit \(u_i\) is for constant-term production, the criterion is \[|(\phi_i(B) - M_i)/M_i)|<\gamma_4,\] where \(\phi_i(x)=\sigma(w_ix+b_i)\) is the activation function of \(u_i\), \(M_i=\max_j{\phi_i(x_j)}\), \(\gamma_4\) is the threshold and \(B=0, 1\) are the two endpoints of \([0, 1]\). If equation 7.17 holds and if the unit is activated, \(u_i\) is for producing a constant term. We set \(\gamma_4 = 0.05\). Under this criterion, in Figure 8a, \(u_{13}\) is the unit predicted by corollary 1, whose activation function is a bold red line. The blue curves correspond to other units.

Note that from equation 7.16, \(\lim_{x\to +\infty}\sigma(x)=1\) is also a constant. Thus, a tanh unit can generate a constant in two ways including \(\sigma(-\infty)=-1\) and \(\sigma(+\infty)=1\); the associated experimental phenomenon is shown in Figure 8b. The function for Figure 8b is \(y=e^{8x}\) whose discrete values are also normalized when depicted. The parameters of network \(\mathfrak{N}\) are \(\Theta=10\), \(c=0.00001\) and \(N=5000\). The threshold \(\gamma_4=0.05\). As can be seen in Figure 8b, units \(u_3\) and \(u_5\) are of nearly constant by \(\sigma(+\infty)=1\), while \(u_9\) and \(u_{10}\) by \(\sigma(-\infty)=-1\). Although one unit is enough for a constant as in Figure 8a, redundant ones can also do that job.

8 Summary for Black Box↩︎

We summarize the main thoughts of the theory as follows. Let \(\mathfrak{N}\) be a two-layer neural network with generalized sigmoidal or tanh units and \(f(\boldsymbol{x}) \in C^m([0,1]^n)\) for \(n\ge1\) be a smooth function. The hidden-layer units of \(\mathfrak{N}\) are denoted by \(u_i\) for \(i=1,2,\dots,\Theta\), whose activation function is \(\phi_i(\boldsymbol{x})=\sigma(\boldsymbol{w}_i^T\boldsymbol{x}+b_i)\).

• Local approximation. Network \(\mathfrak{N}\) can approximate \(f(\boldsymbol{x})\) by realizing its Taylor series expansion. This type of solution not only appears in training solutions but also is a necessary part of global solutions associated with splines.

• Global approximation Over a single strict partial order, network \(\mathfrak{N}\) can approximate \(f(\boldsymbol{x})\) via smooth splines constructed from a piecewise linear approximation to \(f(\boldsymbol{x})\)’s directional-derivative hypersurface. The hidden-layer units of \(\mathfrak{N}\) can be classified into two categories. One is of global units that approximate \(f(\boldsymbol{x})\) through Taylor series expansions; the other is of local units, implementing polynomial pieces by the principle of smooth splines.

• Meaning of the parameters. The weight vector \(\boldsymbol{w}_i\) and bias \(b_i\) of a local unit \(u_i\) are used to generate a knot of splines, while the output weight \(\lambda_i\) produces the polynomial piece corresponding to that knot. The geometric meaning of \(\lambda_i\) can be reduced to the case of two-layer ReLU networks by directional-derivative operations, which is related to the two angles derived from adjacent linear pieces.

• Smooth-continuity restriction. This principle mostly distinguishes network \(\mathfrak{N}\) from other types of multivariate-function approximation. It provides a new representation of smooth functions that is globally correlated and contributes to a fundamental property called “boundary-determination principle” analogous to the boundary-value problem of differential equations.

• Mechanism of realizing splines. Over a single strict partial order, each activation function \(\phi_i(\boldsymbol{x})\) is approximated by a smooth spline \(s_i(\boldsymbol{x})\) and the linear combination of \(s_i(\boldsymbol{x})\)’s leads to a desired spline \(\mathcal{S}(\boldsymbol{x})\) approximating \(f(\boldsymbol{x})\). To obtain the weights of \(s_i(\boldsymbol{x})\)’s, the recurrence relation of the polynomial pieces of a spline and the zero-error part of \(\phi_i(\boldsymbol{x})\) play a central role.

9 Discussion↩︎

To a two-layer neural network \(\mathfrak{N}\) composed of units with smooth activation functions, the conventional method for universal approximation usually resorts to Taylor series expansions, Fourier analysis or step functions. Our investigation of training solutions contributed to the spline solution. Furthermore, for higher-dimensional input, the underlying principle is more than a spline, but also includes smooth-continuity restriction, a significant feature distinguishing \(\mathfrak{N}\) from other types of multivariate-function approximation.

Network \(\mathfrak{N}\) shares some identical or similar basic principles with a two-layer ReLU network, such as multiple strict partial orders, (smooth) continuity restriction, zero-error hyperplanes, polynomial construction and spline implementation. Thus, the future-research problems proposed for ReLU networks in [28] are also applicable to network \(\mathfrak{N}\) with smooth activation functions.

Besides the ones constructed in this paper, more solutions to be discovered are contained in the matrix-form models, including the generalized Wronskian matrix (equations 2.13 and 5.10) and spline matrix (equations 4.20 and 6.54). The study of these matrices would yield more concrete solutions that may explain engineering applications.

Our logically established theory successfully explained the experimental solution obtained by the back-propagation algorithm; and the examples of section 7 showed that the theory can manually construct training solutions through a deterministic way, instead of the usual non-deterministic gradient-descent method. Those achievements demonstrate that the research methodology of theoretical physics can be no doubt applied in the realm of artificial neural networks.

References↩︎