Synthesis of Feedforward Networks in

Supremum Error Bound

Krzysztof Ciesielski, Jaros�aw P. Sacha, and Krzysztof J. Cios

K. Ciesielski is with the Department of Mathematics, West Virginia University, Morgantown, USA; e-mail: K Cies@math.wvu.edu

J. P. Sacha; e-mail: Jarek@IEEE.org K. J. Cios is with the University of Colorado at Denver, USA; e-mail: KCios@carbon.cudenver.edu

SYNTHESIS OF FEEDFORWARD NETWORKS 2

Abstract The main result of this paper is a constructive proof of a formula for the upper bound of the approximation

error in L_{∞} (supremum norm) of multi-dimensional functions by feedforward networks with one hidden layer of sigmoidal units and a linear output. This result is applied to formulate a new method of neural network synthesis. The result can also be used to estimate complexity of the maximum-error network and/or to initialize that network weights. An example of the network synthesis is given.

Keywords Approximation error, supremum norm, synthesis, feedforward networks, sigmoidal functions, complexity,

Radon transform.

I. Introduction

Universal approximation capabilities for a broad range of neural network topologies have been established by researchers like Cybenko [1], Funahashi [2], Hornik et al. [3], Ito [4], and Leshno et al. [5]. Their work focused on the existence of an approximating network and was fundamental for establishing validity of artificial neural networks.

Problem of complexity and realizability of networks with two hidden layers using step (threshold) units was addressed by Blum and Li [6]. The authors presented a method for construction of such networks and specified upper bound of the approximation error. Korain [7] extended approach of Funahashi [2] and established a method for constructing single hidden layer networks with step function neurons and a linear output; he also estimated complexity and the upper bound of the approximation error. A network using either step, ramp, or sigmoid functions was studied by Bulsari [8]. He gave the upper bound of approximation error for networks with one-dimensional input using step and ramp functions. Realizability of networks based on Kolmogorov’s representation theorem was investigated by Ku̇rková [9] and Sprecher [10]. The latter presented a method for constructing transfer functions for hidden layers neurons. Sprecher [11] discussed numerical aspects of implementing Kolmogorov’s representation approach.

More recently an important result has been established for approximation in the L_{2} norm or mean squared error. Several authors (Jones [12], Barron [13], Girosi and Anzellotti [14],

SYNTHESIS OF FEEDFORWARD NETWORKS 3

Murata [15]) have shown for various classes of sufficiently smooth functions that feedforward networks, with one hidden layer of sigmoidal nodes, achieve integrated square error of order O(1/N), where N is the number of neurons in the hidden layer. This result is significant since it demonstrates that approximation with neural networks in L_{2} norm can be free of the “curse of dimensionality,” i.e., approximation rate is independent of the dimensionality of the space. However, the “curse” was avoided only in terms of the accuracy of approximation but not in terms of computational complexity [13].

Learning algorithms like backpropagation are based on minimization of the error in L_{2} norm (mean square error). A number of applications exist, however, where approximation in the L_{∞} norm (maximum error) is more advantageous then using L_{2} norm. A typical example is manipulation of a robot arm, or other control applications, where one requires that the maximum deviation be within certain bounds, e.g., safety limits. To that end an interesting new learning algorithm that performs minimization of the error in L_{∞} norm have been proposed by Meltser et al. [16]. The algorithm is used for networks with a hidden layer of sigmoidal nodes and a linear output layer.

The main result of this paper is a constructive proof of a formula for the upper bound of the approximation error in L_{∞} (supremum norm) of multi-dimensional functions by feedforward networks with one hidden layer of sigmoidal units and a linear output. This result is used to formulate a new method of neural network synthesis. The method specifies network topology and all its weights and can be used to initialize learning algorithms based on the L_{∞} norm, like the one proposed by Meltser et al. [16], or simply to estimate complexity of a maximum-error network. Note, that the synthesis method presented here does not provide for direct learning from examples, it assumes that full information about the function being approximated is given.

In what follows we shall consider a one-hidden-layer network:

N_{N}(x) = c_{0} + ∑_{i=1}^{N} (c_{i} ϕ (w_{i0} + ∑_{j=1}^{n} w_{ij}x_{j})) , N_{N}: R^{n} ↦ R^{1}.  (1)

SYNTHESIS OF FEEDFORWARD NETWORKS 4

For any bounded subset B of R^{n} we can find the one hidden layer network N for which

sup_{x∈B} |f(x) − N (x)| < ε,

where f: R^{n} ↦ R^{1} is the target function and ε is the desired approximation error. A constructive solution to this problem, in which the complexity of the approximating network (number of neurons in the hidden layer) is O_{ε}(1/ε^{n}), will be presented.

The proposed network synthesis method consist of three stages. First, the problem of multidimensional approximation is replaced by several simpler one-dimensional problems. Second, each of the one-dimensional problems is solved by finding an approximation subnetwork. Third, the subnetworks are combined together to form the final approximation network.

In the first stage function f is approximated by a finite sum of plane waves, or ridge functions b_{l}:

f(x) ≈ ν ∑_{l=1}^{L} b_{l}(x), b_{l}: R^{n} ↦ R, ∀ l,

where ν is a constant. A characteristic feature of a plain wave is that it can be represented by a one-dimensional function and a directional vector [17]. Let us denote this one-dimensional function by b_{l}: R ↦ R and the directional vector by ω^{(l)} ∈ R^{n}. Then

b_{l}(x) = b_{l}(x · ω^{(l)}), ∀ x ∈ R^{n}.

Figure 1 shows example of a plain wave:

x = [ \begin{matrix}x_{1} \\ x_{2} \end{matrix} ] , ω^{(l)} = [ \begin{matrix}2 \\ 1 \end{matrix} ] , b_{l}(x · ω^{(l)}) = sinc(x · ω^{(l)}) = sinc(2x_{1} + x_{2}),

where

sinc(t) = \frac{sin(πt)}{πt}.

In the second stage of the method, one-dimensional functions b_{l} associated with plane waves are approximated be single-input single-output subnetworks ℵ^{(l)} : R ↦ R. First, functions b_{l} are represented as a difference of two monotonic functions. Second, neural networks are synthesized to approximate each of these monotonic functions. The most significant implication of the

Fig. 1. A plain wave b_{l}(x_{1}, x_{2}) = sinc(2x_{1} + x_{2}).

second part of the method is constructive specification of the upper bound of an approximation error for networks with sigmoidal functions.

The third stage combines subnetworks ℵ^{(l)} into the final approximation network N:

N (x) = ν ∑_{l=1}^{L} ℵ^{(l)} (x · ω^{(l)}) .

The reduction of dimensionality (the first approximation stage) is done by means of the inverse Radon transform of function f. It is based on the proposed here revision of Carroll and Dickinson [18] method. Some of the constrains posed by these authors on the function being approximated will be eliminated and, at the same time, the upper bound of the approximation error, by a finite number of plain waves, will be lowered.

Carroll and Dickinson [18] used inverse Radon transform approach to provide an alternative proof of universal approximation properties of neural networks with mean square error norm. They first decompose a mutltidimensional function into scalar functions using inverse Radon transform. Then they approximate scalar functions with single layer neural networks using steplike approximation — weights of sigmoidal neurons are made to tend to infinity so that the neural transfer functions approach step functions. However, they do not provide actual values of network parameters needed to achieve the desired accuracy. (This part of their proof is existential.)

SYNTHESIS OF FEEDFORWARD NETWORKS 6

The method presented in this paper was inspired by the work of Carroll and Dickinson [18], but is similar to it only in the way it uses inverse Radon transform to decomposed multidimensional functions into scalar functions. In contrast to Carroll and Dickinson [18] we are interested in maximum approximation error rather than mean square error. We also perform approximation of scalar functions differently. We use finite neural weights and provide the actual values needed to achieve the desired error rate. The weights are assigned so that derivative of the neural transfer function near the center of the scaled-shifter sigmoid is finite and close to the derivative of the approximated function for the similar value of the argument.

Application of the inverse Radon transform for construction of representation networks was also studied by Ito [19]. He concentrated on delivering the results for uniform approximation on the entire space R^{n} (not on a compact subset) and use of a step and sigmoidal functions without scaling. However, his mainly constructive proof is existential in the part specifying partitioning for the approximation of the integral representation (not only of the unit sphere, as in our case, but of the entire space R^{n}), and no estimation of the approximation error was given.

The method of network synthesis presented in this paper is non-iterative. The decision on the network size is made before the weights of the network are calculated. This decision is usually based on the required upper bound of the approximation error. The method can also be used to estimate number of neurons in the hidden layer given the approximation error limit and the target function, without the need of calculating network weights.

We start by presenting our initial result: the constructive proof that any one (real) variable monotone function can be approximated, with an arbitrary accuracy, by a neural network with sigmoidal transfer functions. Next, we show how this result can be used to approximate multidimension functions. We prove the upper bound of the maximum error of the multi-dimensional approximation and show an illustrative example of applying our method for network synthesis. Section V gives an example of applying the proposed method. Some remarks on practical use of the method are presented in Section VI. Proofs of the theorems and corollaries are given in the Appendix.

SYNTHESIS OF FEEDFORWARD NETWORKS 7

II. Approximation of One-dimensional Monotone Functions We shall start with discussion of approximation with ramp transfer functions. It is included

here because it illustrates well the intuition behind the more complex approach — approximation with sigmoidal transfer functions.

Let ψ: R ↦ R be continuous strictly increasing function defined on an interval K = [a, b]. Function ψ maps K onto interval Q = [ψ(a), ψ(b)]. Let us assume that the network used for approximation of ψ has N neurons in the hidden layer. Next, denote:

\begin{aligned}Δ & = ψ(b) − ψ(a), \\ δ_{i} & = ψ(a) + i \frac{Δ}{N} for i = 0, 1, …, N, \end{aligned}

and divide interval Q into N equal subintervals Q_{i} = [δ_{i−1}, δ_{i}], i = 1, 2, …, N. Define:

γ_{i} = ψ^{−1}(δ_{i}) for i = 0, 1, …, N.

A. Ramp Transfer Functions

Results presented in this section are similar to well-known results using polynomial and spline approximation. They are presented here in the context of neural networks. The following definition of the ramp function is used:

ϕ(x) = \{ \begin{matrix}−1 & x < −1, \\ x & −1 ≤ x ≤ 1, \\ 1 & x > 1. \end{matrix}  (2)

Let each of the neurons in the hidden layer correspond to one of the intervals Q_{i} defined above; they “approximate” function ψ in these intervals. Assume that the mapping of the network ℵ and the function ψ take the same values at each γ_{i}:

ℵ(γ_{i}) = ψ(γ_{i}) = δ_{i} for i = 0, 1, …, N.

Figure 2 illustrates this situation (the approximation network has N = 5 neurons in the hidden layer). Circles represent points (γ_{i}, δ_{i}) in which values of ψ and ℵ are assumed to be the same.

Fig. 2. Approximation of an increasing function with ramp transfer functions.

Shaded areas correspond to regions where the target function and its approximation are allowed to differ. Since both ψ and ℵ are monotonically increasing the difference is less then δ_{i} − δ_{i−1} = \frac{Δ}{N}. This value constitutes the upper error bound of approximating a monotonic function with ramp functions. It is summarized in the following proposition.

Proposition 1: Let function ψ: R^{1} ↦ R^{1} be continuous strictly increasing on [a, b], a neural network given by (1), the ramp transfer function ϕ given by (2), and the network weights be assigned as follows:

c_{0} = \frac{ψ(a) + ψ(b)}{2},
c_{i} = \frac{ψ(b) − ψ(a)}{2N},
γ_{i} = ψ^{−1} [ψ(a) + i \frac{ψ(b) − ψ(a)}{N}] ,
w_{i0} = −\frac{γ_{i} + γ_{i−1}}{γ_{i} − γ_{i−1}} ,
w_{i1} = \frac{2}{γ_{i} − γ_{i−1}} .

Then, the upper bound of the approximation error of function ψ by network ℵ on interval [a, b] is:

sup_{x∈[a,b]} |ℵ(x) − ψ(x)| ≤ \frac{ψ(b) − ψ(a)}{N}.

□

Fig. 3. Approximation of an increasing function with sigmoidal transfer functions.

B. Sigmoidal Transfer Functions Let us assume that a sigmoidal transfer function is the hyperbolic tangent,

ϕ(x) = tanh(x).

The results presented in this section work for any other increasing continuous function ϕ with bounded range, but parameters have to be slightly changed.^{1} We derive the result for the hyperbolic tangent since it is commonly used in neural networks.

Since ramp functions vary only on a finite interval, it was reasonable to require that ψ and ℵ take the same values at points (γ_{i}, δ_{i}), see Figure 2. However, this requirement is not practical for a transfer function like hyperbolic tangent or other sigmoids which change their values on an infinite interval. In order to find the solution and the error’s upper bound for sigmoids we shall relax this matching requirement at a cost of a slightly higher value of the approximation error, but without compromising the approximation rate. The following theorem specifies the network weights and the the upper bound of the approximation error.

Theorem 2: Let ψ: [a, b] → R be a strictly increasing continuous function, and for every

1The parameters remain the same for any such ϕ for which ϕ(−x) = −ϕ(x) and lim_{x→∞} ψ(x) = 1.

natural number N let ℵ_{N}: [a, b] → R be defined by the formula (1):

ℵ_{N}(x) = c_{0} + ∑_{i=1}^{N} c_{i}ϕ(w_{i1}x + w_{i0})

where ϕ(x) = tanh x and parameters c_{i}’s and w_{jk}’s are defined as below. Then

sup_{x∈[a,b]} |ℵ_{N}(x) − ψ(x)| < \frac{3Δ}{N},

where

Δ = ψ(b) − ψ(a).

In particular ℵ_{N} converges to ψ (as N → ∞) in sup-norm. Here we use, for i = 1, …, N,

d_{N} = ϕ^{−1} (1 − \frac{1}{2N} ) ,
γ_{0} = a,
γ_{i} = ψ^{−1} (ψ(a) + i\frac{Δ}{N} ) ,
w_{i0} = d_{N} − \frac{2d_{N}}{γ_{i} − γ_{i−1}} γ_{i} = −d_{N} \frac{γ_{i} + γ_{i−1}}{γ_{i} − γ_{i−1}} ,  (3)
w_{i1} = \frac{2d_{N}}{γ_{i} − γ_{i−1}} ,
c_{i} = \frac{Δ}{2N} ,
c_{0} = ψ(a) − ∑_{i=1}^{N} c_{i}ϕ(w_{i1}a + w_{i0}).
□

III. Approximation of One-dimensional Functions

Any function of bounded variation can be represented as a difference of two monotonic increasing functions [20]. Thus we can approximate any continuous function of bounded variation with a neural network using the result presented in the previous section. Let g be a continuous function of bounded variation. Function g can be represented as difference of two monotone increasing functions ψ^{+} and ψ^{−}

g(x) = ψ^{+}(x) − ψ^{−}(x)  (4)

using Jordan decomposition^{2} [20]:

\begin{aligned}ψ^{+}(x) & = \frac{1}{2} (T_{g}[a, x] + g(x)) \\ ψ^{−}(x) & = \frac{1}{2} (T_{g}[a, x] − g(x)) \end{aligned}  (5)

where T_{g}[a, b] is the total variation of function g on interval [a, b]:

T_{g}[a, b] = sup ∑_{i=0}^{k−1} |f(α_{i+1}) − f(α_{i})|, a = α_{0} < α_{1} < ⋯ < α_{k} = b.

To ensure that ψ^{+} and ψ^{−} are strictly increasing we can modify the above by the following assignment:

\begin{aligned}ψ^{+}(x) & = \frac{1}{2} (T_{g}[a, x] + g(x) + x), \\ ψ^{−}(x) & = \frac{1}{2} (T_{g}[a, x] − g(x) + x). \end{aligned}  (6)

Strictly increasing functions ψ^{+} and ψ^{−} can be approximated by networks ℵ_{N^{+}}^{+} and ℵ_{N^{−}}^{−}, respectively. The network approximating function g is given by

ℵ_{N}(x) = ℵ_{N^{+}}^{+}(x) − ℵ_{N^{−}}^{−}(x),  (7)

where N = N^{+} + N^{−}.

Lemma 3: Let g: [a, b] → R be of bounded variation. For every natural number N let ℵ_{N}: [a, b] → R be defined by the formula (7). Then ℵ_{N} is in the form of formula (1):

ℵ_{N}(x) = c_{0} + ∑_{i=1}^{N} c_{i}ϕ(w_{i1}x + w_{i0})

where ϕ(x) = tanh x and coefficients are chosen from those for ℵ_{N^{+}}^{+} and ℵ_{N^{−}}^{−}. Moreover

sup_{x∈[a,b]} |ℵ_{N}(x) − g(x)| < \frac{6}{N} T_{g}[a, b].  (8)

□

2If ψ is continuously differentiable then the decomposition can be easier described by the formulas: ψ^{+}(x) = g(a) + ∫_{a}^{x} max{0, g′(t)} dt and ψ^{−}(x) = ∫_{a}^{x} min{0, −g′(t)} dt.

SYNTHESIS OF FEEDFORWARD NETWORKS 12

IV. Approximation of Multi-dimensional Functions

Multi-dimensional functions are first decomposed into one-dimensional functions using plain waves. One-dimensional functions are approximated using results from the preceding sections.

A. Plain Wave Decomposition

Let function f: R^{n} ↦ R be infinitely differentiable rapidly decreasing.^{3} Let f̌ denote the Radon transform of function f [21]:

f̌(ω, t) = ∫_{ω·x = t} f(x) ξ  (9)

where ω is a unit vector in R^{n}, ω · x is a scalar product, and ξ is a differential form used to integrate f over all hyperplanes ω · x = t.

The Radon transform inversion formula

f(x) = ∫_{Φ} b(ω · x, ω) ξ(ω)  (10)

is an integration of plain waves, b(ω · x, ω), over any closed surface Φ enclosing the origin of space R^{n}. A plain wave can be computed as follows [21]:

b(t, ω) = \{ \begin{matrix}\frac{(−1)^{\frac{1}{2} (n−1)}}{2(2π)^{n−1}} (\frac{∂}{∂t})^{n−1}f̌(ω, t) & for odd n, \\ \frac{(−1)^{\frac{1}{2} n}(n−1)!}{(2π)^{n}} ∫_{Φ} ∫_{−∞}^{+∞} [f̌(ω, t)\{t − ω · x\}^{−n}dt] ξ(ω) & for even n, \end{matrix}  (11)

where the integral over t is understood in terms of its regularization, and

ξ(ω) = ∑_{k=1}^{n} (−1)^{k−1}ω_{k}dω_{1} ⋯ dω_{k−1}dω_{k+1} ⋯ dω_{n}.

In the context of the inverse Radon transform the plain waves b(ω · x, ω) are traditionally called the back projection data. We use integral representation (10) as a basis for reduction of dimensionality, utilizing the fact that b(ω · x, ω) is a one-dimensional function for any fixed ω.

Assume some suitable partitioning of the integration surface

Φ = {Φ_{l}}_{l=1}^{L}.

3Function f defined on R^{n} is rapidly decreasing when lim_{x→∞} |x_{1}^{k_{1}} ⋯ x_{n}^{k_{n}} | f(x) = 0 for any nonnegative k_{i}’s.

Now, approximate integral (10) by the sum:

f(x) ≈ f̄(x) = ∑_{l=1}^{L} μ(Φ_{l}) b(ω^{(l)} · x, ω^{(l)}).  (12)

where

μ(Φ_{l}) = ∫_{Φ_{l}} ξ(ω).

The upper bound of the approximation error of (12) is given by the following theorem. This theorem is a correction of the result presented in [18].

Theorem 4: Let f: R^{n} ↦ R be infinitely differentiable rapidly decreasing. Let f̌ be the Radon transform of f and b(t, ω) be the corresponding back-projection data. Set the number of partitions to L = k^{n−1} where k is a positive integer. Then the error of approximation (12) is bounded by:

|f(x) − f̄(x)| ≤ \frac{α}{k},  (13)

α = \frac{2π^{\frac{n+2}{2}} \sqrt{n}}{Γ (\frac{n}{2})} max_{t} max_{ω} |∇_{ω}b(t, ω)| ,

where the gamma function Γ is given by:

Γ(x) = ∫_{0}^{∞} z^{x−1}e^{−z}dz.
□

Notice, that α is constant for fixed f, hence the error can be arbitrarily small as the number of partitions increases.

B. Multi-dimensional Representation by Neural Networks

For a fixed ω^{(l)} = (ω_{1}^{(l)}, …, ω_{n}^{(l)}), l = 1, …, k^{n−1}, the function b(ω^{(l)} · x, ω^{(l)}) is onedimensional. To simplify notation denote t = ω^{(l)} · x, and b_{l}(t) = b(ω^{(l)} · x, ω^{(l)}). Each function b_{l} can be approximated using decomposition discussed in Section II and Theorem 2 by appropriate ℵ_{η}^{(l)} : R ↦ R. We shall assume that each of the subnetworks ℵ_{η}^{(l)} has the same number of neurons in the hidden layer, η. To represent n-dimensional back-projection data mappings b(ω^{(l)} · x, ω^{(l)}) there is a need to implement the inner product t = ω^{(l)} · x. This can be done by

SYNTHESIS OF FEEDFORWARD NETWORKS 14

extending the number of inputs of the subnetwork from 1 to n, and identifying the weights w_{ij}, i = 1, …, ηk^{n−1}, with w_{(i,l)j}, i = 1, …, η and l = 1, …, k^{n−1}, defined as follows

w_{(i,l)j} ≐ ω_{j}^{(l)} w_{i1}^{(l)}, j = 1, …, n.  (14)

Weights c_{i} and w_{i0}, for i = 1, …, ηk^{n−1}, are identified with appropriate c_{i}^{(l)} and w_{i0}^{(l)}, i = 1, …, η and l = 1, …, k^{n−1}, and remain unchanged. The mapping of network N_{N}: R^{n} ↦ R, the final approximated representation of function f: R^{n} ↦ R, is the superposition of the mappings ℵ_{η}^{(l)}

N_{N}(x) = ∑_{l=1}^{k^{n−1}} μ(Φ_{l}) ℵ_{η}^{(l)}(x) = c_{0} + ∑_{i=1}^{N} (c_{i} ϕ (w_{i0} + ∑_{j=1}^{n} w_{ij}x_{j})) ,  (15)

where c_{0} = ∑_{l=1}^{k^{n−1}} μ(Φ_{l}) c_{0}^{(l)}.

Theorem 5: Let f: R^{n} ↦ R be infinitely differentiable rapidly decreasing. Let N_{N} be a single hidden layer network, given by (1), with N = ηk^{n−1} neurons in the hidden layer, ψ(x) = tanh(x). Then there exists an approximated representation of function f by network N_{N} such that

|f(x) − N_{N}(x)| ≤ \frac{α}{k} + \frac{β}{η}

where

α = \frac{2π^{\frac{n+2}{2}} \sqrt{n}}{Γ(\frac{n}{2})} max_{t} max_{ω} |∇_{ω}b(t, ω)| ,
β = \frac{12π^{\frac{n}{2}}}{Γ (\frac{n}{2})} max_{l} T_{b_{l}}[K_{l}],

and K_{l} is the respective domain of function b_{l}(t). □

Observe that for a given function f parameters α and β are constant. Theorem 5 specifies the upper bound of the multi-dimensional approximation error using single hidden-layer network. On the other hand, it can be used to find the number of hidden layer neurons needed to guarantee the required approximation error.

Corollary 6: The upper bound of the approximation error of infinitely differentiable rapidly decreasing function f: R^{n} ↦ R by a single hidden layer network using sigmoidal (hyperbolic

SYNTHESIS OF FEEDFORWARD NETWORKS 15

tangent) function neurons is proportional to the inverse of the n^{th}–root of the number of neurons in the hidden layer:

ε ∼ \frac{1}{\sqrt[n]{N}}

In other words function f can be realized with accuracy ε by a single hidden-layer network using N = O_{ε}(1/ε^{n}) neurons in the hidden layer. □

V. Example

This section shows an example of the network synthesis method proposed in the paper. The method performs three steps: Step 1: Find approximation of function f by a finite number of plain waves b(t, ω^{(l)}) using

formulas (9), (11), and (12).

Step 2: Decompose one-dimensional functions b_{l}(t), associated with plain waves b(t, ω^{(l)}), into monotonic functions b_{l}^{+}(t) and b_{l}^{−}(t) using formula (6). Then find subnetworks approximating each of these functions using Theorem 2.

Step 3: Find final approximation network by combining subnetworks found in Step 2 using

formulas (14) and (15). Let us consider the following two-dimensional function f: R^{2} ↦ R

f(x_{1}, x_{2}) = 3 (1 − x_{1})^{2} exp[−x_{1}^{2} − (x_{2} + 1)^{2}]
− 10 (\frac{1}{5}x_{1} − x_{1}^{3} − x_{2}^{5}) exp[−x_{1}^{2} − x_{2}^{2}]  (16)

− \frac{1}{3} exp[−(x_{1} + 1)^{2} − x_{1}^{2}],

A plot of the function f is shown in Figure 4. A two-dimensional example is chosen so that all of the steps of the method can be graphically illustrated. The method was numerically implemented using MATLAB.

Step 1. Function f is approximated by a finite number of plain waves, using formula (12). Let as assume initially that the number of plain waves used is L = 6. Figures 5 through 7 show some of the plain waves b(t, ω)|_{t=x·ω} obtained from f. Note, that in R^{2} a unit vector can be

Fig. 4. Function f(x_{1}, x_{2}).  Fig. 5. Plain wave b(x_{1}, x_{2})|_{θ_{1}=0^{o}}.

Fig. 6. Plain wave b(x_{1}, x_{2})|_{θ_{1}=−30^{o}}.  Fig. 7. Plain wave b(x_{1}, x_{2})|_{θ_{1}=−60^{o}}.

parameterized by a single parameter θ, which is the angle of rotation of the unit vector. We have:

ω = [ \begin{matrix}cos(θ) \\ sin(θ) \end{matrix} ] .

Each of the plain waves is represented by a one-dimensional function b_{l}(t) and a unit vector ω^{(l)} or, equivalently, by a rotation angle θ. For L = 6

θ ∈ {0^{o}, 30^{o}, 60^{o}, 90^{o}, 120^{o}, 150^{o}}.

Figures 8 and 9 show approximation of function f(x_{1}, x_{2}) using 6 and 12 plain waves, respectively.

Fig. 8. Representation with 6 plain waves.  Fig. 9. Representation with 12 plain waves.

Fig. 10. b_{l}(t)|_{θ_{1}=−30^{o}}.  Fig. 11. b_{l}_{+}(t)|_{θ_{1}=−30^{o}}.  Fig. 12. b_{l}_{−}(t)|_{θ_{1}=−30^{o}}.

Step 2. Each of the one-dimensional functions b_{l}(t) is decomposed into two monotonic increasing functions, b_{l}^{+}(t) and b_{l}^{−}(t), using formula (5). Then, these monotonic functions are approximated by single hidden-layer networks, as is illustrated in Figures 10 through 12. Solid lines represent functions b_{l}(t), b_{l}^{+}(t) and b_{l}^{−}(t), respectively. Dotted lines represent their approximation with neural networks based on Theorem 2. In this example nine neurons in the hidden layer are used. Figure 13 compares the upper bound of the approximation error obtained from formula (8) with the actual approximation error for different numbers of neurons in the hidden layer of the approximating network. Note, that the scale on the vertical axis is logarithmic.

Step 3. The networks approximating functions b_{l}(t) are combined together, using formulas (14) and (15), to form the final approximation network. Examples of approximation of the function (16) with L = 12, N = 9 and L = 12, N = 18 are shown in Figures 14 and 15, respectively.

Fig. 13. Approximation error and its upper bound.

Fig. 14. Approximation with 12 plain waves and 9 neurons.

Fig. 15. Approximation with 12 plain waves and 18

neurons.

VI. Remarks

Although the emphasis of the paper is on giving analytical solution of the problem of neural network synthesis assuming a perfect information about the function being approximation, we would like to comment on what happens when this is not the case. We addressed assumptions on the function, computational error, and complexity of the method. It should be stress that we just touch here the tip of the iceberg and that more precise treatment of these issues calls for a separate study.

SYNTHESIS OF FEEDFORWARD NETWORKS 19

A. Information about a function and computational error

In general, to apply our method for estimating function f: R^{n} ↦ R^{1} on a bounded subset B of R^{n} we need to assume that our function is rapidly decreasing^{4} infinitely many times differentiable on some open subset U containing the closure of B. This requirement comes from the use of Radon transform for reduction of the multi-dimensional case to one dimensional problem. This assumption serves well the current discussion, however, it can be relaxed in more general treatments of the Radon transform [21], [17], [22].

Consider case when our original function f is one dimensional. It can be equal to a function g: R^{1} ↦ R^{1}. This reduces our requirements on f essentially, since from g we require only that it is continuous and of bounded variation (see Section III). The situation becomes especially easy when the number M = sup_{x,y∈B, x≠y} |\frac{g(y)−g(x)}{y−x} | is finite, for example, when g has a continuous derivative. In this case for any M < m < ∞ the representation g(x) = ψ^{+}(x) − ψ^{−}(x) of g as a difference of two strictly monotone increasing functions can be achieved just by putting ψ^{+}(x) = \frac{1}{2}(mx + g(x)) and ψ^{−}(x) = \frac{1}{2}(mx − g(x)). Note, that an error ε of an estimation of M has no impact on this representation as long as M + ε < m. It is also interesting to realize that our estimation of a strictly increasing function ψ: [a, b] ↦ R obtained this way (i.e., ψ equals ψ^{+} or ψ^{−}) has very clear computational complexity and does not depend much on error of calculating function ψ. More precisely, in order to find reasonable estimations of γ_{i}’s (see Section IIB) used in calculating coefficients in our representation we need to calculate the values of ψ only at points a+jδ for j = 0, …, k, where k is the smallest number for which \frac{m(b−a)}{k} < Δ, and δ = \frac{b−a}{k}.^{5} Thus the number of computations needed for finding an approximation of ψ is of order: constant (≤ 10) times \frac{m(b−a)}{Δ}. Finally note that if the values of ψ are calculated within an \frac{Δ}{N} error, this can increase the error of the ℵ estimation of ψ just by \frac{Δ}{N}. (This is the case since our

4In practice, we do not need assume the “rapidly decreasing” part, since we estimate our function only on the bounded set B. More precisely, instead of estimating f we can always apply our method to the function f_{0} = f · g, where g: R^{n} ↦ R^{1} is infinitely many times differentiable, is equal 1 on B, and is equal 0 on the complement of U. Then f_{0} is equal to f on B, while it satisfies the “rapidly decreasing” requirement.

5Estimate γ_{i} by the first a + jδ for which |ψ(a + jδ) − (ψ(a) + i \frac{Δ}{N} )| < \frac{Δ}{N}. This will increase the error estimate in Theorem 2 from \frac{3Δ}{N} to \frac{4Δ}{N}.

SYNTHESIS OF FEEDFORWARD NETWORKS 20

method for finding the parameters is non-iterative; and in calculating the values of ℵ_{N}(x) the method essentially relies only to on the information to which interval [γ_{i−1}, γ_{i}] point x belongs.)

B. Complexity of the method

Let assume that computation of a plain wave b_{l} from a function f has complexity O_{b_{l}}, and the complexity evaluation of function b_{l} for a given x is negligible. Then complexity of onedimensional approximation, O_{1D} is proportional to the number of neurons in the hidden layer η

O_{1D} = O(η).

Complexity of multi-dimensional approximation is proportional to complexity O_{b_{l}}, number of one dimensional approximations k^{n−1} and complexity of each of these approximations O(η)

O_{nD} = η · k^{n−1} · O_{b_{l}} = N · O_{b_{l}},

where N is the total number of neurons in a hidden layer of a neural network. Thus complexity of approximating multi-dimensional function f is directly proportional to the complexity of computing a plain wave b_{l}. Applications of Radon transform are gaining in popularity in recent years. Initially there were restricted to subjects related to computational medical imaging. Recently the transform is used in many new fields like optics and holographic interferometry, geophysics, radio astronomy, and pure mathematics. A good discussion of relatively recent advances in Radon transform and its applications can be found in [23], [22]. In general, computation of plain waves and inverse Radon transform is a challenging problem. However, there have been a significant progress in addressing practical issues of inverse Radon transform computations, for instance wavelet transform based methods introduced in [24], [25].

VII. Conclusions

The goal of this paper was twofold. First, it was to determine the upper bound of the approximation error for single hidden-layer feedforward networks using sigmoidal transfer function neurons. Second, it was to develop the method for synthesis of networks to approximate func-

tions of several variables.

A constructive proof have been presented that the upper bound of the maximum approxi-

mation error using sigmoidal (hyperbolic tangent) function neurons is proportional to 1/\sqrt[n]{N},

where n is dimension of the space and N is the number of neurons in the hidden layer. The

proof of the error’s upper bound was used to introduce a new method of network synthesis.

The presented method for network synthesis consists of three main parts. In part one, the

dimensionality of the approximation problem is reduced by applying the inverse Radon trans-

form. In the second part, one-dimensional problems are approximately represented by single

hidden-layer subnetworks. Then these subnetworks are combined into the final approximation

network. The method fully determines the network topology, including the neural weights.

References

[1] G. Cybenko, “Approximation by superpositions of a sigmoidal function,” Mathematics of Control, Signals and Systems, vol. 2, pp. 303–314, 1989.

[2] K.-I. Funahashi, “On the approximate realization of continuous mappings by neural networks,” Neural Networks, vol. 2, pp. 183–192, 1989.

[3] K. Hornik, M. Stinchcombe, and H. White, “Universal approximation of unknown mapping and its derivatives using multilayer feedforward networks,” Neural Networks, vol. 3, pp. 551–560, 1990.

[4] Y. Ito, “Approximation of continuous functions on R^{d} by linear combination of of shifted rotations of sigmoid function with and without scaling,” Neural Networks, vol. 5, pp. 105–115, 1992.

[5] M. Leshno, V. Y. Lin, A. Pinkus, and S. Schocken, “Multilayer feedforward networks with a nonpolynomial activation function can approximate any function,” Neural Networks, vol. 6, pp. 861–867, 1993.

[6] E. K. Blum and L. K. Li, “Approximation theory and feedforward networks,” Neural Networks, vol. 4, pp. 511–515, 1991.

[7] P. Korain, “On the complexity of approximating mappings using neural networks,” Neural Networks, vol. 6, pp. 649–653, 1993.

[8] A. Bulsari, “Some analytical solutions to the general approximation problem for feedforward neural networks,” Neural Networks, vol. 6, pp. 991–996, 1993.

[9] V. Ku̇rková, “Kolmogorov’s theorem and multilayer neural networks,” Neural Networks, vol. 5, pp. 501–506, 1992.

[10] D. A. Sprecher, “A universal mapping for Kolmogorov’s superposition theorem,” Neural networks, vol. 6, pp. 1089–1094, 1993.

[11] D. A. Sprecher, “A numerical implementation of Kolmogorov’s superpositions,” Neural Networks, vol. 9, no. 5, pp. 765–772, 1996.

[12] L. K. Jones, “A simple lemma on greedy approximation in hilbert space and convergence rates for projection pursuit regression and neural network training.,” Annals of Statistics, vol. 20, no. 1, pp. 608–613, 1992.

[13] A. R. Barron, “Universal approximation bounds for superpositions of a sigmoidal function,” IEEE Transactions on Information Theory, vol. 39, no. 3, pp. 930–945, May 1993.

[14] F. Girosi and G. Anzellotti, “Convergence rates of approximation by translates,” Tech. Rep. A.I. Memo No. 1288, MIT, Artificial Intelligence Laboratory, April 1995.

[15] N. Murata, “An integral representation of functions using three-layered networks and their approximation bounds,” Neural Networks, vol. 9, no. 6, pp. 947–956, 1996.

[16] M. Meltser, M. Shoham, and L. M. Manevitz, “Approximating functions by neural networks: A constructive solution in the uniform norm,” Neural Networks, vol. 9, no. 6, pp. 965–978, 1996.

[17] S. Helgason, The Radon Transform, Brikhäuser, 1980.

[18] S. M. Carroll and B. W. Dickinson, “Construction of neural nets using the Radon transform,” in ’89 IJCNN Proceedings, 1989, vol. I, pp. 607–611.

[19] Y. Ito, “Representation of functions by superposition of a step or sigmoidal functions and their applications to neural network theory,” Neural Networks, vol. 4, pp. 385–394, 1991.

[20] E. J. McShane, Integration, Princeton University Press, Princeton, 1944.

[21] I. M. Gel’fand, M. I. Graev, and N. Y. Vilenkin, Generalized Functions, vol. 5, Academic Press, 1966.

[22] S. R. Deans, “Radon and Abel transforms,” in The Transforms and Applications Handbook, A. D. Poularikas, Ed., pp. 631–717. CRC Press, 1996.

[23] S. R. Deans, The Radon Transform and Some of Its Applications, Krieger Publishing Company, Malabar, Florida, 1993.

[24] G. Kaiser and R. F. Streater, “Windowed Radon transform, analytic signals, and thw wave function,” in Wavelets: A Tutorial in Theory and Applications, C. K. Chui, Ed., pp. 399–441. Academic Press, San Diego, CA, 1992.

[25] M. Holschneider, “Inverse Radon transform through inverse wevelet transform,” Inverse Problems, vol. 7, pp. 853–861, 1991.

Appendix

Proof of Proposition 1: The proof is identical to that of Theorem 2 if you define d_{N} = 1. This

also makes respective weights the same. □

Proof of Theorem 2: Note that the choice of c_{0} guarantees that ℵ_{N}(a) = ψ(a) and the choice

of w_{jk}’s is made to ensure that L_{i}(x) = w_{i1}x + w_{i0} maps [γ_{i−1}, γ_{i}] onto [−d_{N}, d_{N}], that is that

F_{i}(x) = ϕ(L_{i}(x)) = ϕ(w_{i1}x + w_{i0}) maps [γ_{i−1}, γ_{i}] onto [−1 + \frac{1}{2N} , 1 − \frac{1}{2N} ].

SYNTHESIS OF FEEDFORWARD NETWORKS 23

To see the estimation fix N, 0 < j ≤ N, and x ∈ [γ_{j−1}, γ_{j}] and note that

ℵ_{N}(x) = ψ(a) + ∑_{i=1}^{N} \frac{Δ}{2N} [ϕ(w_{i1}x + w_{i0}) − ϕ(w_{i1}a + w_{i0})]

= ψ(a) + \frac{Δ}{2N} ∑_{i=1}^{N}[F_{i}(x) − F_{i}(a)]

= ψ(a) + \frac{Δ}{2N} (∑_{i=1}^{j−1}[F_{i}(x) − F_{i}(a)] + [F_{j}(x) − F_{j}(a)] + ∑_{i=j+1}^{N}[F_{i}(x) − F_{i}(a)])

Now, for i > j we have x ≤ γ_{j} ≤ γ_{i−1} and so

0 < F_{i}(x) − F_{i}(a) < F_{i}(γ_{i−1}) − (−1) = (−1 + \frac{1}{2N} ) + 1 = \frac{1}{2N} .

Similarly for i < j we have x ≥ γ_{j−1} ≥ γ_{i} and so

2 > F_{i}(x) − F_{i}(a) > F_{i}(γ_{i}) − F_{i}(γ_{i−1}) = 2 − \frac{1}{N} .

Thus,

|ℵ_{N}(x) − ψ(a) − Δ\frac{2(j − 1)}{2N}| ≤ \frac{Δ}{2N} (∑_{i=1}^{j−1} |[F_{i}(x) − F_{i}(0)] − 2| + |F_{j}(x) − F_{j}(0)|

+ ∑_{i=j+1}^{N} |F_{i}(x) − F_{i}(0)|)

≤ \frac{Δ}{2N} (\frac{j − 1}{N} + 2 + \frac{N − (j + 1)}{N})
≤ \frac{3Δ}{2N} ≤ \frac{2Δ}{N}

Since we also have

|ψ(x) − ψ(a) − Δ\frac{2(j − 1)}{2N}| = |ψ(x) − ψ(γ_{j−1})| ≤ \frac{Δ}{N}

we conclude that indeed |ℵ_{N}(x) − ψ(x)| ≤ \frac{Δ}{N} + \frac{2Δ}{N} = \frac{3Δ}{N} finishing the argument. □ Proof of Lemma 3: Using Eq. (4), (6), Theorem 2, and assuming N^{+} = N^{−} (N^{+} and N^{−} are

SYNTHESIS OF FEEDFORWARD NETWORKS 24

number of neurons in approximating networks, N^{+} + N^{−} = N) we have:

sup_{x∈[a,b]} |ℵ_{N}(x) − g(x)| ≤ sup_{x∈[a,b]} |ℵ_{N}^{+}(x) − ψ^{+}(x)| + sup_{x∈[a,b]} |ℵ_{N}^{−}(x) − ψ^{−}(x)|
< 3 \frac{ψ^{+}(b) − ψ^{+}(a)}{N^{+}} + 3 \frac{ψ^{−}(b) − ψ^{−}(a)}{N^{−}}
= \frac{6}{N} \frac{1}{2} (T_{g}[a, b] + g(b) − T_{g}[a, a] − g(a)
+ T_{g}[a, b] − g(b) − T_{g}[a, a] + g(a)])
= \frac{6}{N} T_{g}[a, b].
□

Proof of Theorem 4: Assume fixed x. Let Φ be the unit sphere and assume uniform partitioning

|f(x) − f̄(x)| = |∑_{l=1}^{L} ∫_{Φ_{l}} b(ω · x, ω) ξ(ω) − ∑_{l=1}^{L} μ(Φ_{l})b(ω^{(l)} · x, ω^{(l)})|

≤ ∑_{l=1}^{L} ∫_{Φ_{l}} |b(ω · x, ω) − b(ω^{(l)} · x, ω^{(l)})| ξ(ω)

≤ ∑_{l=1}^{L} μ(Φ_{l}) max_{ω∈Φ_{l}} d_{μ}(ω, ω^{(l)}) max_{ω∈Φ_{l}} |∇_{ω}b(ω · x, ω)|

≤ μ(Φ) max_{l=1...L} max_{ω∈Φ_{l}} d_{μ}(ω, ω^{(l)}) max_{ω} |∇_{ω}b(ω · x, ω)| ,

where d_{μ}(ω, ω^{(l)}) is a distance between points ω and ω^{(l)} on the unit sphere. Since the partitioning is uniform we have

max_{l=1...L} max_{ω∈Φ_{l}} d_{μ}(ω, ω^{(l)}) = \frac{π\sqrt{n}}{k}.

The surface of the unit sphere is

Ω_{n} = μ(Φ) = \frac{2π^{\frac{n}{2}}}{Γ(\frac{n}{2})}.

Hence

|f(x) − f̄(x)| ≤ \frac{2}{k} \frac{π^{\frac{n+2}{2}} \sqrt{n}}{Γ(\frac{n}{2})} max_{ω} |∇_{ω}b(ω · x, ω)| .

SYNTHESIS OF FEEDFORWARD NETWORKS 25

Maximizing over x we obtain the upper bound of approximation error:

|f(x) − f̄(x)| ≤ \frac{2}{k} \frac{π^{\frac{n+2}{2}} \sqrt{n}}{Γ(\frac{n}{2})} max_{x} max_{ω} |∇_{ω}b(ω · x, ω)|

= \frac{2}{k} \frac{π^{\frac{n+2}{2}} \sqrt{n}}{Γ(\frac{n}{2})} max_{t} max_{ω} |∇_{ω}b(t, ω)| .

□

Proof of Theorem 5: Form Lemma 3 we have:

|b_{l}(t) − ℵ^{(l)}(t)| < \frac{6}{N_{l}} T_{b_{l}}[K_{l}].

The error of approximation of function f̄by network N, formula (12), is the sum of approximation errors of functions b_{l} by networks ℵ^{(l)}. Assume that each network ℵ^{(l)} has the same number of neurons N_{l} = η, hence

|f̄(x) − N_{N}(x)| ≤ k^{n−1} μ(Φ_{l}) \frac{6}{η} max_{l} T_{b_{l}}[K_{l}] = \frac{12}{η} \frac{π^{\frac{n}{2}}}{Γ (\frac{n}{2})} max_{l} T_{b_{l}}[K_{l}].  (17)

Thus from inequalities (13) and (17)

|f(x) − N_{N}(x)| ≤ \frac{2}{k} \frac{π^{\frac{n+2}{2}} \sqrt{n}}{Γ(\frac{n}{2})} max_{t} max_{ω} |∇_{ω}b(t, ω)| + \frac{12}{η} \frac{π^{\frac{n}{2}}}{Γ (\frac{n}{2})} max_{l} T_{b_{l}}[K_{l}].

□

Proof of Corollary 6: From Theorem 5 we have that approximation’s error upper bound is

ε = α \frac{1}{k} + β \frac{1}{η},  (18)

Let us relax conditions on k and η and assume that they can take any positive value, not only integer. Further assume that number of units in the hidden layer, N = k^{n−1}η, is constant and then find minimum value of ε.

Substitute η = \frac{N}{k^{n−1}} into (18)

ε = α \frac{1}{k} + β \frac{k^{n−1}}{N} ,  (19)

SYNTHESIS OF FEEDFORWARD NETWORKS 26

Differentiate ε in respect to k

\frac{dε}{dk} = −α \frac{1}{k^{2}} + β \frac{(n − 1) k^{n−2}}{N}

\frac{dε}{dk} = 0 ⟹ k = \sqrt[n]{\frac{αN}{β(n − 1)}}.  (20)

Check if it is a minium of function ε(k)

\frac{d^{2}ε}{dk^{2}} = 2 α \frac{1}{k^{3}} + β \frac{(n − 1)(n − 2) k^{n−3}}{N} .  (21)

Formula (21) is positive for any positive value of k, hence we have the minimum. Now, substitute (20) into (19)

ε = [α^{\frac{n−1}{n}} (β(n − 1))^{\frac{1}{n}} + β(n − 1) (\frac{α}{β(n − 1)})^{\frac{n−1}{n}} ] \frac{1}{\sqrt[n]{N}} .  (22)

Thus, we have proportional relation

ε ∼ \frac{1}{\sqrt[n]{N}}

since other components of formula (22) are constant for given function f. And

ε ∼ \frac{1}{\sqrt[n]{N}} ⟹ N ∼ \frac{1}{ε^{n}} .

□