Nonlinear Approximation via Compositions
Zuowei Shen Haizhao Yang Shijun Zhang


An almost-solvable model of complex network dynamics
    Qi Guo^1, Artur Sowa^1
    December 30, 2023
====================================================




Given a function dictionary 𝒟 and an approximation budget N∈ℕ^+,  nonlinear approximation seeks the linear combination of the best N terms {T_n}_1≤ n≤ N⊆𝒟 to approximate a given function f with the minimum approximation error 
    ϵ_L,f:=min_{g_n}⊆ℝ,{T_n}⊆𝒟f(x)-∑_n=1^N g_n T_n(x).
 Motivated by recent success of deep learning, we propose dictionaries with functions in a form of compositions, i.e., 
    T(x)=T^(L)∘ T^(L-1)∘⋯∘ T^(1)(x)
 for all T∈𝒟, and implement T using ReLU feed-forward neural networks (FNNs) with L hidden layers. We further quantify the improvement of the best N-term approximation rate in terms of N when L is increased from 1 to 2 or 3 to show the power of compositions. In the case when L>3, our analysis shows that increasing L cannot improve the approximation rate in terms of N.

In particular, for any function f on [0,1], regardless of its smoothness and even the continuity, if f can be approximated using a dictionary when L=1 with the best N-term approximation rate ϵ_L,f=(N^-η), we show that dictionaries with L=2 can improve the best N-term approximation rate to ϵ_L,f=(N^-2η). We also show that for Hölder continuous functions of order α on [0,1]^d, the application of a dictionary with L=3 in nonlinear approximation can achieve an essentially tight best N-term approximation rate ϵ_L,f=( N^-2α/d). Finally, we show that dictionaries consisting of wide FNNs with a few hidden layers are more attractive in terms of computational efficiency than dictionaries with narrow and very deep FNNs for approximating Hölder continuous functions if the number of computer cores is larger than N in parallel computing. 



  Deep Neural Networks, ReLU Activation Function, Nonlinear Approximation, Function Composition, Hölder Continuity, Parallel Computing.










§ INTRODUCTION


For non-smooth and high-dimensional function approximation, a favorable technique popularized in recent decades is the nonlinear approximation <cit.> that does not limit the approximants to come from linear spaces, obtaining sparser representation, cheaper computation, and more robust estimation, and therein emerged the bloom of many breakthroughs in applied mathematics and computer science (e.g., wavelet analysis <cit.>, dictionary learning <cit.>, data compression and denoising <cit.>, adaptive pursuit <cit.>, compressed sensing <cit.>).  

Typically, nonlinear approximation is a two-stage algorithm that designs a good redundant nonlinear dictionary, 𝒟, in its first stage, and identifies the optimal approximant as a linear combination of N elements of 𝒟 in the second stage:

    f(x) ≈ g∘ T(x) :=∑_n=1^N g_n T_n(x),

where f(x) is the target function in a Hilbert space ℋ associated with a norm denoted as ·_*, {T_n}⊆𝒟⊆ℋ, T is a nonlinear map from ℝ^d to ℝ^N with the n-th coordinate being T_n, and g is a linear map from ℝ^N to ℝ with the n-th coordinate being g_n∈ℝ. The nonlinear approximation seeks g and T such that

    {{T_n},{g_n}}=_{g_n}⊆ℝ,{T_n}⊆𝒟f(x) - ∑_n=1^N g_n T_n(x)_*,

which is also called the best N-term approximation. 

One remarkable approach of nonlinear approximation is based on one-hidden-layer neural networks that give simple and elegant bases of the form T(x)=σ(Wx+b), where Wx+b is a linear transform in x with the transformation matrix W (named as the weight matrix) and a shifting vector b (called bias), and σ is a nonlinear function (called the activation function). The approximation

    f(x)≈∑_n=1^N g_n T_n(x)= ∑_n=1^N g_nσ(W_n x + b_n)

includes wavelets pursuits <cit.>, adaptive splines <cit.>, radial basis functions <cit.>,  sigmoidal neural networks <cit.>, etc. For functions in Besov spaces with smoothness s, <cit.> constructed an (N^-s/d)[In this paper, we use the big (·) notation when we only care about the scaling in terms of the variables inside (·) and the prefactor outside (·) is independent of these variables.] approximation that is almost optimal <cit.> and the smoothness cannot be reduced generally <cit.>. For Hölder continuous functions of order 1 on [0,1]^d, <cit.> essentially constructed an (N^-1/2d) approximation, which is far from the lower bound (N^-2/d) as we shall prove in this paper. Achieving the optimal approximation rate of general continuous functions in constructive approximation, especially in high dimensional spaces, remains an unsolved challenging problem.



 §.§ Problem Statement





ReLU FNNs have been proved to be a powerful tool in many fields from various points of view <cit.>, which motivates us to tackle the open problem above via function compositions in the nonlinear approximation using deep ReLU FNNs, i.e., 

    f(x) ≈ g∘ T^(L)∘ T^(L-1)∘⋯∘ T^(1)(x),

where T^(i)(x)=σ(W_i x + b_i) with W_i∈^N_i× N_i-1, b_i∈^N_i for i=1,…,L, σ is the ReLU activation function, and f is a Hölder continuous function. For the convenience of analysis, we consider N_i=N for i=1,…,L. Let 𝒟_L be the dictionary consisting of ReLU FNNs g∘ T^(L)∘ T^(L-1)∘⋯∘ T^(1)(x) with width N and depth L. 
To identify the optimal FNN to approximate f(x), it is sufficient to solve the following optimization problem

    ϕ^*=_ϕ∈𝒟_Lf - ϕ_*.
 

The fundamental limit of nonlinear approximation via the proposed dictionary is essentially determined by the approximation power of function compositions in (<ref>), which gives a performance guarantee of the minimizer in (<ref>). Since function compositions are implemented via ReLU FNNs, the remaining problem is to quantify the approximation capacity of deep ReLU FNNs, especially their ability to improve the best N-term approximation rate in N for any fixed L defined as

    ϵ_L,f(N)=min_ϕ∈𝒟_Lf - ϕ_*.


Function compositions can significantly enrich the dictionary of nonlinear approximation and this idea was not considered in the literature previously due to the expensive computation of function compositions in solving the minimization problem in (<ref>). Fortunately, recent development of efficient algorithms for optimization with compositions (e.g., backpropagation techniques  <cit.> and parallel computing techniques <cit.>) makes it possible to explore the proposed dictionary in this paper. Furthermore, with advanced optimization algorithms <cit.>, good local minima of (<ref>) can be identified efficiently <cit.>. 















 §.§ Related work and contribution

The main goal in the remaining article is to quantify the best N-term approximation rate ϵ_L,f(N) defined in (<ref>) for ReLU FNNs in the dictionary 𝒟_L with a fixed depth L when f is a Hölder continuous function. This topic is related to several existing approximation theories in the literature, but none of these existing works can be applied to answer the problem addressed in this paper.


First of all, this paper identifies explicit formulas for the best N-term approximation rate 

    ϵ_L,f(N)≤
    2 ν N^-2α,    when  L≥ 2 and d=1,
    
    2(2√(d))^αν N^-2α/d,    when L≥ 3 and  d>1,

for any N∈ℕ^+ and a Hölder continuous function f of order α with a constant ν, while existing theories <cit.> can only provide implicit formulas in the sense that the approximation error contains an unknown prefactor and work only for sufficiently large N or L larger than some unknown numbers. For example, the approximation rate in <cit.> via a narrow and deep ReLU FNN is c(d) L^-2α/d with c(d) unknown and for L larger than a sufficiently large unknown number ℒ; the approximation rate in <cit.> via a wide and shallow (c_1(d)-layer) ReLU FNN is c_2(d) N^-α/d with c_1(d) and c_2(d) unknown and for N larger than a sufficiently large unknown number 𝒩. For another example, given an approximation error ϵ, <cit.> proved the existence of a ReLU FNN with a constant but still unknown number of layers approximating a C^β function within the target error. Similarly, given the ϵ error, <cit.> estimate the scaling of the network size in ϵ and the scaling contains unknown prefactors. Given an arbitrary L and N, no existing work can provide an explicit formula for the approximation error to guide practical network design, e.g., to guarantee whether the network is large enough to meet the accuracy requirement. This paper provides such formulas for the first time and in fact the bound in these formulas is asymptotically tight as we shall prove later. 

Second, our target functions are Hölder continuous, while most of existing works aim for a smaller function space with certain smoothness, e.g. functions in C^α([0,1]^d) with α≥ 1 <cit.>, band-limited functions <cit.>, Korobov spaces <cit.>, or Besev spaces <cit.>. To the best of our knowledge, there is only one existing article <cit.> concerning the approximation power of deep ReLU FNNs for C([0,1]^d). However, the conclusion of <cit.> only works for ReLU FNNs with a fixed width 2d+10 and a sufficiently large L, instead of a fixed L and an arbitrary N as required in the nonlinear approximation (see Figure <ref> for the comparison of the conclusion of <cit.> and this paper).

As we can see in Figure <ref>, the improvement of the best N-term approximation rate in terms of N when L is increased from 1 to 2 or 3 is significant, which shows the power of depth in ReLU FNNs. However, in the case when L>3, our analysis shows that increasing L cannot improve the approximation rate in terms of N. As an interesting corollary of our analysis, for any function f on [0,1], regardless of its smoothness and even the continuity, if f can be approximated using using functions in 𝒟_1 with the best N-term approximation rate ϵ_L,f=(N^-η), we show that functions in 𝒟_2 can improve the best N-term approximation rate to ϵ_L,f=(N^-2η). Extending this conclusion for a general d dimensional function is challenging and we leave it as future work.

 From the point of view of analysis techniques, this paper introduce new analysis methods merely based on the structure of FNNs, while existing works <cit.> rely on constructing FNNs to approximate traditional basis in approximation theory, e.g., polynomials, splines, and sparse grids, which are used to approximate smooth functions.  




Finally, we analyze the approximation efficiency of neural networks in parallel computing, a very important point of view that was not paid attention to in the literature. In most applications, the efficiency of deep learning computation highly relies on parallel computation. We show that a narrow and very deep neural network is inefficient if its approximation rate is not exponentially better than wide and shallower networks.  Hence, neural networks with (1) layers are more attractive in modern computational platforms, considering the computational efficiency per training iteration in parallel computing platforms. Our conclusion does not conflict with the current state-of-the-art deep learning research since most of these successful deep neural networks have a depth that is asymptotically (1) relative to the width.


















 §.§ Organization


The rest of the paper is organized as follows. Section <ref> summarizes the notations throughout this paper. Section <ref> presents the main theorems while Section <ref> shows numerical tests in parallel computing to support the claims in this paper. Finally, Section <ref> concludes this paper with a short discussion.



§ PRELIMINARIES



For the purpose of convenience, we present notations and elementary lemmas used throughout this paper as follows.



 §.§ Notations


 
     
     
  * Matrices are denoted by bold uppercase letters, e.g., A∈ℝ^m× n is a real matrix of size m× n, and A^T denotes the transpose of A. Correspondingly, A(i,j) is the (i,j)-th entry of A; A(:,j) is the j-th column of A; A(i,:) is the i-th row of A.
     
  * Vectors are denoted as bold lowercase letters, e.g., v∈^n is a column vector of size n and v(i) is the i-th element of v. v=[v_1,⋯,v_n]^T=[[ v_1;   ⋮; v_n ]] are vectors consisting of numbers {v_i} with v(i)=v_i.
     
     
  * The Lebesgue measure is denoted as μ(·). 
     
     
  * The set difference of A and B is denoted by A\ B. A^c denotes [0,1]^d\ A for any A⊆ [0,1]^d. 
     
     
  * For a set of numbers A, and a number x, A-x:={y-x:y∈ A}.
     
     
  * For any ξ∈, let ⌊ξ⌋:=max{i: i≤ξ, i∈ℤ} and ⌈ξ⌉:=min{i: i≥ξ, i∈ℤ}.
     
     
     
  * Assume n∈^n, then f(n)=𝒪(g(n)) means that there exists positive C independent of n, f, and g such that f(n)≤ Cg(n) when n(i) goes to +∞ for all i. 
     
     
  * Define Lip(ν,α,d) as the class of functions defined on [0,1]^d satisfying the uniformly Lipchitz property of order α with a Lipchitz constant ν>0.
     That is, any f∈Lip(ν,α,d) satisfies
     
    |f(x)-f(y)|≤νx-y_2^α,  for any x,y∈ [0,1]^d.

     
     
  * Let CPL(N) be the set of continuous piecewise linear functions with N-1 pieces mapping [0,1] to . The endpoints of each linear piece are called “break points" in this paper.
     
     
  * Let σ:→ denote the rectified linear unit (ReLU), i.e. σ(x)=max{0,x}.  With the abuse of notations, we define σ:^d→^d as σ()=[[ max{0,x_1};          ⋮; max{0,x_d} ]] for any =[x_1,⋯,x_d]^T∈^d.
     
     
  * We will use NN as a ReLU neural network for short and use Python-type notations to specify a class of NN's, e.g., NN(c_1; c_2; ⋯; c_m) is a set of ReLU FNN's satisfying m conditions given by {c_i}_1≤ i≤ m, each of which may specify the number of inputs (), the total number of nodes in all hidden layers (#node), the number of hidden layers (#layer), the total number of parameters (#parameter), and the width in each hidden layer (widthvec), the maximum width of all hidden layers (maxwidth) etc. For example, NN(=2;=[100,100]) is a set of NN's ϕ satisfying:
     
         
  * ϕ maps from ^2 to .
         
  * ϕ has two hidden layers and the number of nodes in each hidden layer is 100.
     
     
  * [n]^L is short for [n,n,⋯,n]∈^L. For example, 
    NN(=d;=[100,100,100])=NN(=d;=[100]^3).

     
  * For ϕ∈NN(=d;=[N_1,N_2,⋯,N_L]), if we define N_0=d and N_L+1=1, then the architecture of ϕ can be briefly described as follows:
    
    x=h̃_0 ^W_1, b_1h_1⟶^σh̃_1  ⋯^W_L, b_Lh_L⟶^σh̃_L ^W_L+1, b_L+1ϕ(x)=h_L+1,

    where W_i∈^N_i× N_i-1 and b_i∈^N_i are the weight matrix and the bias vector in the i-th linear transform in ϕ, respectively, i.e., h_i :=W_i h̃_i-1 + b_i for i=1,2,⋯,L+1 and h̃_i=σ(h_i) for i=1,2,⋯,L.

     


 §.§ Lemmas




Let us study the properties of ReLU FNNs with only one hidden layer to warm up in Lemma <ref> below. It indicates that CPL(N+1)=NN(=1;=[N]) for any N∈^+. 


    
    Suppose ϕ∈NN(=1;=[N]) with an architecture:
    
    x ^W_1, b_1h^σh̃^W_2, b_2ϕ(x).

Then ϕ is a continuous piecewise linear function. Let W_1=[1,1,⋯,1]^T∈^  N× 1, then we have:
    
        
  * Given a sequence of strictly increasing numbers x_0, x_1,⋯, x_N, there exists b_1∈^  N independent of W_2 and b_2 such that the break points of ϕ are exactly x_0, ⋯, x_N on the interval [x_0,x_N][We only consider the interval [x_0,x_N] and hence x_0 and x_N are treated as break points. ϕ(x) might not have a real break point in a small open neighborhood of x_0 or x_N.].
        
  * Suppose {x_i}_i∈{0,1,⋯,N} and b_1 are given in (1). 
        Given any sequence {y_i}_i∈{0,1,⋯,N}, there exist W_2 and b_2 such that ϕ(x_i)=y_i for i=0,1,⋯,N and ϕ is linear on [x_i,x_i+1] for i=0,1,⋯,N-1. 
    


Part (1) in Lemma <ref> follows by setting b_1=[-x_0,-x_1,⋯,-x_N-1]^T. The existence in Part (2) is equivalent to the existence of a solution of linear equations, which is left for the reader. Next, we study the properties of ReLU FNNs with two hidden layers. In fact, we can show that the closure of  NN(=1;=[2m,2n+1]) contains CPL(mn+1) for any m,n∈^+, where the closure is in the sense of L^p-norm for any p∈ [1,∞). The proof of this property relies on the following lemma.

    
    For any m,n∈^+, given any m(n+1)+1 samples (x_i,y_i)∈^2 with x_0<x_1<x_2<⋯<x_m(n+1) and y_i≥ 0 for i=0,1,⋯,m(n+1), there exists ϕ∈NN(=1;=[2m,2n+1]) satisfying the following conditions:
    
        
  * ϕ(x_i)=y_i for i=0,1,⋯,m(n+1);
        
  * ϕ is linear on each interval [x_i-1,x_i] for i∉{(n+1)j:j=1,2,⋯,m};
        
  * sup_x∈[x_0, x_m(n+1)] |ϕ(x)| ≤ 3max_i∈{0,1,⋯,m(n+1)}y_i ∏_k=1^n(1+max{x_j(n+1)+n-x_j(n+1)+k-1:j=0,1,⋯,m-1}min{x_j(n+1)+k-x_j(n+1)+k-1:j=0,1,⋯,m-1}).
    


    For simplicity of notation, we define I_0(m,n){0,1,⋯,m(n+1)}, I_1(m,n){j(n+1):j=1,2,⋯,m}, and I_2(m,n) I_0(m,n)\ I_1(m,n) for any m,n∈^+.
    Since ϕ∈NN(=1;=[2m,2n+1]), the  architecture of ϕ is

    x ^W_1, b_1h^σh̃^W_2, b_2g^σg̃^W_3, b_3ϕ(x).

Note that g maps x∈ to g(x)∈^2n+1 and hence each entry of g(x) itself is a sub-network with one hidden layer. Denote g=[g_0,g_1^+,g_1^-,⋯,g_n^+,g_n^-]^T, then {g_0,g_1^+,g_1^-,⋯,g_n^+ ,g_n^-}⊆NN(=1;=[2m]). Our proof of Lemma <ref> is mainly based on the repeated applications of Lemma <ref> to determine parameters of ϕ(x) such that Conditions (1) to (3) hold.

Step 1: Determine W_1 and b_1. 


By Lemma <ref>, ∃ W_1=[1,1,⋯,1]^T∈^ 2m× 1 and b_1∈^ 2m× 1 such that sub-networks in {g_0,g_1^+,g_1^-,⋯,g_n^+ ,g_n^-} have the same set of break points: {x_i: i∈ I_1(m,n)∪(I_1(m,n)-1) ∪{0}}, no matter what W_2 and b_2 are. 

Step 2:  Determine W_2 and b_2. 


This is the key step of the proof. Our ultimate goal is to set up g=[g_0,g_1^+,g_1^-,⋯,g_n^+,g_n^-]^T such that, after a nonlinear activation function, there exists a linear combination in the last step of our network (specified by W_3 and b_3 as shown in (<ref>)) that can generate a desired ϕ(x) matching the sample points {(x_i,y_i)}_0≤ i≤ m(n+1). In the previous step, we have determined the break points of {g_0,g_1^+,g_1^-,⋯,g_n^+ ,g_n^-} by setting up W_1 and b_1; in this step, we will identify W_2∈^(2n+1)× 2m and b_2∈^2n+1 to fully determine {g_0,g_1^+,g_1^-,⋯,g_n^+ ,g_n^-}. This will be conducted in two sub-steps.

Step 2.1: Set up.


Suppose f_0(x) is a continuous piecewise linear function defined on [0,1] fitting the given samples f_0(x_i)=y_i for i∈ I_0(m,n), and f_0 is linear between any two adjacent points of {x_i:i∈ I_0(m,n)}. We are able to choose W_2(1,:) and b_2(1) such that g_0(x_i)=f_0(x_i) for i∈ I_1(m,n)∪ (I_1(m,n)-1)∪{0} by Lemma <ref>, since there are 2m+1 points in I_1(m,n)∪ (I_1(m,n)-1)∪{0}. Define f_1:=f_0-g̃_0, where g̃_0=σ (g_0)=g_0 as shown in Equation (<ref>), since g_0 is positive by the construction of Lemma <ref>. Then we have f_1 (x_i)=f_0 (x_i)-g̃_0 (x_i)=0 for i∈  (I_1 (m,n)-n-1) ∪{m (n+1)}. See Figure <ref> (a) for an illustration of f_0, f_1, and g_0.

Step 2.2: Mathematical induction.


For each k∈{1,2,⋯,n}, given f_k, we determine W_2 (2k,:), b_2 (2k), W_2 (2k+1,:), and b_2 (2k+1), to completely specify 
g^+_k and g^-_k, which in turn can determine f_k+1. Hence, it is only enough to show how to proceed with an arbitrary k, since the initialization of the induction, f_1, has been constructed in Step 2.1. See Figure <ref> (b)-(d) for the illustration of the first two induction steps. We recursively rely on the fact of f_k that

    
  * f_k (x_i)=0 for i∈∪_ℓ=0^k-1 (I_1 (m,n)-n-1+ℓ)∪{m (n+1)},
    
  * f_k is linear on each interval [x_i-1,x_i] for i∈ I_2 (m,n)\{0},



to construct f_k+1 satisfying similar conditions as follows: 

    
  * f_k+1 (x_i)=0 for i∈∪_ℓ=0^k (I_1 (m,n)-n-1+ℓ)∪{m (n+1)},
    
  * f_k+1 is linear on each interval [x_i-1,x_i] for i∈ I_2 (m,n)\{0}.



The induction process for W_2 (2k,:), b_2 (2k), W_2 (2k+1,:), b_2 (2k+1), and f_k+1 can be divided into four parts.

Step 2.2.1: Define index sets.


Let Λ^+_k={j: f_k (x_j (n+1)+k)≥ 0,  0≤ j<m} and Λ^-_k={j: f_k (x_j (n+1)+k)< 0, 0≤ j<m}. The cardinality of Λ^+_k∪Λ^-_k is m. We will use Λ^+_k and Λ^-_k to generate 2m+1 samples to determine CPL functions g^+_k(x) and g^-_k(x) in the next step.


Step 2.2.2:  Determine W_2 (2k,:) and b_2 (2k). 


By Lemma <ref>, we can choose W_2 (2k,:) and b_2 (2k) to fully determine g^+_k(x) such that each g^+_k (x_i) matches a specific value for i∈  (I_1 (m,n)-n-1)∪  (I_1 (m,n)-1)∪{m (n+1)}. The values of {g^+_k (x_i):i∈(I_1 (m,n)-n-1)∪  (I_1 (m,n)-1)∪{m (n+1)}} are specified as:   

    
  * If j∈Λ^+_k, specify the values of g^+_k (x_j (n+1)) and g^+_k (x_j (n+1)+n) such that
    g^+_k(x_j (n+1)+k-1)=0 and g^+_k(x_j (n+1)+k)=f_k (x_j (n+1)+k). The existence of these values fulfilling the requirements above comes from the fact that g^+_k (x) is linear on the interval [x_j (n+1),x_j (n+1)+n] and g^+_k (x) only depends on the values of g^+_k(x_j (n+1)+k-1) and g^+_k(x_j (n+1)+k) on [x_j (n+1),x_j (n+1)+n]. Now it is easy to verify that g̃^+_k(x):=σ (g^+_k(x)) satisfies g̃^+_k (x_j (n+1)+k)= f_k (x_j (n+1)+k)≥ 0 and g̃^+_k (x_j (n+1)+ℓ) = 0 for ℓ=0,1,⋯,k-1, and g̃^+_k is linear on each interval [x_j (n+1)+ℓ,x_j (n+1)+ℓ+1] for ℓ=0,1,⋯,n-1.
    
  * If j∈Λ^-_k, let g^+_k (x_j (n+1))=g^+_k (x_j (n+1)+n)=0. Then g̃^+_k(x)=0 on the interval [x_j (n+1),x_j (n+1)+n].
    
  * Finally, specify the value of g^+_k(x) at x=x_m(n+1) as 0.




Step 2.2.3: Determine W_2 (2k+1,:) and b_2 (2k+1). 


Similarly, we choose W_2 (2k+1,:) and b_2 (2k+1) such that g^-_k(x) matches specific values as follows:

    
  * If j∈Λ^-_k, specify the values of g^-_k (x_j (n+1)) and g^-_k (x_j (n+1)+n) such that 
    g^-_k(x_j (n+1)+k-1)=0  and g^-_k(x_j (n+1)+k)=-f_k (x_j (n+1)+k). Then g̃^-_k(x):=σ (g^-_k(x)) satisfies g̃^-_k (x_j (n+1)+k)=-f_k (x_j (n+1)+k)> 0 and g̃^-_k (x_j (n+1)+ℓ)=0 for ℓ=0,1,⋯,k-1, and g̃^-_k(x) is linear on each interval [x_j (n+1)+ℓ,x_j (n+1)+ℓ+1] for ℓ=0,1,⋯,n-1.
    
  * If j∈Λ^+_k, let g^-_k (x_j (n+1))=g^-_k (x_j (n+)+n)=0. Then g̃^-_k(x)=0 on the interval [x_j (n+1),x_j (n+1)+n].
    
  * Finally, specify the value of g^-_k(x) at x=x_m(n+1) as 0.
 

Step 2.2.4: Construct f_k+1 from g^+_k and g^-_k.


For the sake of clarity,  the properties of g^+_k and g^-_k constructed in Step 2.2.3 are summarized below:

    
  * f_k (x_i)=g̃_k^+ (x_i)=g̃_k^- (x_i)=0 for  i∈∪_ℓ=0^k-1 (I_1 (m,n)-n-1+ℓ)∪{m (n+1)};
    
    
  * If j∈Λ^+_k,  g̃_k^+ (x_j (n+1)+k)=f_k (x_j (n+1)+k)≥ 0  and g̃_k^- (x_j (n+1)+k)=0;

    
  * If j∈Λ^-_k,  g̃_k^- (x_j (n+1)+k)=-f_k (x_j (n+1)+k)>0  and g̃_k^+ (x_j (n+1)+k)=0;
    
    
  * g̃_k^+ and g̃_k^- are linear on each interval [x_j (n+1)+ℓ,x_j (n+1)+ℓ+1] for ℓ=0,1,⋯,n-1, j∈Λ^+_k∪Λ^-_k={0,1,⋯,m-1}. In other words, g̃_k^+ and g̃_k^- are linear on each interval [x_i-1,x_i] for i∈ I_2 (m,n)\{0}. 
 
See Figure <ref> (a)-(c) for the illustration of g_0, g^+_1, g^-_1, g^+_2, and g^-_2, and to verify their properties as listed above.

Note that Λ^+_k∪Λ^-_k={0,1,⋯,m-1}, so f_k (x_i)-g̃_k^+ (x_i)+g̃_k^- (x_i)=0 for i∈∪_ℓ=0^k (I_1 (m,n)-n-1+ℓ)∪{m (n+1)}.
Now we define f_k+1:=f_k-g̃^+_k+g̃^-_k, then

    
  * f_k+1 (x_i)=0 for i∈∪_ℓ=0^k (I_1 (m,n)-n-1+ℓ)∪{m (n+1)};
    
  * f_k+1 is linear on each interval [x_i-1,x_i] for i∈ I_2 (m,n)\{0}.

See Figure <ref> (b)-(d) for the illustration of f_1, f_2, and f_3, and to verify their properties as listed just above. This finishes the mathematical induction process. As we can imagine based on Figure <ref>, when k increases, the support of f_k shrinks to the “don't-care" region. 

Step 3: Determine W_3 and b_3.


With the special vector function g=[g_0,g_1^+,g_1^-,⋯,g_n^+,g_n^-]^T constructed in Step 2, we are able to specify W_3 and b_3 to generate a desired ϕ(x) with a well-controlled L^∞-norm matching the samples {(x_i,y_i)}_0≤ i≤ m(n+1).

In fact, we can simply set W_3=[1,1,-1,1,-1,⋯,1,-1]∈^1× (2n+1) and b_3=0, which finishes the construction of ϕ(x). The rest of the proof is to verify the properties of ϕ(x). Note that ϕ=g̃_0+∑_ℓ=1^ng̃^+_ℓ-∑_ℓ=1^ng̃^-_ℓ. By the mathematical induction, we have:

    
  * f_n+1=f_0-g̃_0-∑_ℓ=1^ng̃^+_ℓ+∑_ℓ=1^ng̃^-_ℓ;
    
    
  * f_n+1 (x_i)=0 for i∈∪_ℓ=0^n (I_1 (m,n)-n-1+ℓ)∪{m (n+1)}=I_0 (m,n);
    
    
  * f_n+1 is linear on each interval [x_i-1,x_i] for i∈ I_2 (m,n)\{0}.



Hence, ϕ=g̃_0+∑_ℓ=1^ng̃^+_ℓ-∑_ℓ=1^ng̃^-_ℓ=f_0-f_n+1. Then ϕ satisfies Conditions (1) and (2) of this lemma. It remains to check that ϕ satisfies Condition (3). 

By the definition of f_1, we have 

    sup_x∈[x_0, x_m(n+1)]|f_1 (x)| ≤ 2max{y_i:i∈ I_0 (m,n)}.

By the induction process in Step 2, for k∈{1,2,⋯,n}, it holds that

    sup_x∈[x_0, x_m(n+1)]|g̃^+_k(x)|≤max{x_j(n+1)+n-x_j(n+1)+k-1:j=0,1,⋯,m-1}min{x_j(n+1)+k-x_j(n+1)+k-1:j=0,1,⋯,m-1}sup_x∈[x_0, x_m(n+1)]|f_k(x)|

and 

    sup_x∈[x_0, x_m(n+1)]|g̃^-_k(x)|≤max{x_j(n+1)+n-x_j(n+1)+k-1:j=0,1,⋯,m-1}min{x_j(n+1)+k-x_j(n+1)+k-1:j=0,1,⋯,m-1}sup_x∈[x_0, x_m(n+1)]|f_k(x)|.

Since either g̃^+_k(x) or g̃^-_k(x) is equal to 0 on [0,1], we have 

    sup_x∈[x_0, x_m(n+1)]|g̃^+_k(x)-g̃^-_k(x)|≤max{x_j(n+1)+n-x_j(n+1)+k-1:j=0,1,⋯,m-1}min{x_j(n+1)+k-x_j(n+1)+k-1:j=0,1,⋯,m-1}sup_x∈[x_0, x_m(n+1)]|f_k(x)|.

Note that f_k+1=f_k-g̃^+_k+g̃^-_k, which means 

    sup_x∈[x_0, x_m(n+1)]|f_k+1(x)|   ≤sup_x∈[x_0, x_m(n+1)]|g̃^+_k(x)-g̃^-_k(x)|+sup_x∈[x_0, x_m(n+1)]|f_k(x)|
       ≤(max{x_j(n+1)+n-x_j(n+1)+k-1:j=0,1,⋯,m-1}min{x_j(n+1)+k-x_j(n+1)+k-1:j=0,1,⋯,m-1}+1) sup_x∈[x_0, x_m(n+1)]|f_k(x)|

for k∈{1,2,⋯,n}. Then we have

    sup_x∈[x_0, x_m(n+1)]|f_n+1(x)|
    ≤ 2max_i∈ I_0 (m,n) y_i∏_k=1^n(max{x_j(n+1)+n-x_j(n+1)+k-1:j=0,1,⋯,m-1}min{x_j(n+1)+k-x_j(n+1)+k-1:j=0,1,⋯,m-1}+1) .

Hence 

    sup_x∈[x_0, x_m(n+1)]|ϕ(x)|
       ≤sup_x∈[x_0, x_m(n+1)]|f_0(x)|+sup_x∈[x_0, x_m(n+1)]|f_n+1(x)|
       ≤ 3max_i∈ I_0 (m,n) y_i∏_k=1^n(max{x_j(n+1)+n-x_j(n+1)+k-1:j=0,1,⋯,m-1}min{x_j(n+1)+k-x_j(n+1)+k-1:j=0,1,⋯,m-1}+1) .

So, we finish the proof.




§ MAIN RESULTS



We present our main results in this section. First, we quantitatively prove an achievable approximation rate in the N-term nonlinear approximation by construction, i.e., the lower bound of the approximation rate. Second, we show a lower bound of the approximation rate asymptotically, i.e., no approximant exists asymptotically following the approximation rate. Finally, we discuss the efficiency of the nonlinear approximation considering the approximation rate and parallel computing in FNNs together. 




 §.§ Quantitative Achievable Approximation Rate


    

    
    For any N∈^+ and f∈Lip(ν,α,d) with α∈(0,1], we have[It is easy to generalize the results in Theorem <ref> and Corollary <ref> from L^1 to L^p-norm for p∈[1,∞) since μ([0,1]^d)=1.]:
    
    
  * If d=1, ∃ ϕ∈NN(=1;=[2N,2N+1]) such that
    
    ϕ-f_L^1([0,1])≤ 2 ν N^-2α,   for any  N∈^+;

    
  * If d>1, ∃ ϕ∈NN(=d;=[2d⌊ N^2/d⌋,2N+2,2N+3]) such that
    
    ϕ-f_L^1([0,1]^d)≤ 2(2√(d))^αν N^-2α/d,   for any  N∈^+.

    



Without loss of generality, we assume f(0)=0 and ν=1. 

Step 1: The case d=1.


Given any f∈Lip (ν,α,d) and N∈^+, we know |f (x)|≤ 1 for any x∈ [0,1] since f(0)=0 and ν=1. Set f̅=f+1≥ 0, then 0≤f̅(x)≤ 2 for any x∈ [0,1]. Let X={i/N^2:i=0,1,⋯,N^2}∪{i/N-δ:i=1,2,⋯,N}, where δ is a sufficiently small positive number depending on N, and satisfying (<ref>). Let us order X as x_0< x_1<⋯<x_N (N+1). By Lemma <ref>, given the set of samples {(x_i,f̅ (x_i)):i∈{0,1,⋯,N (N+1)}}, there exists ϕ∈NN (=1;=[2N,2N+1]) such that

    
  * ϕ (x_i)=f̅ (x_i) for i=0,1,⋯,N(N+1);
    
  * ϕ is linear on each interval [x_i-1,x_i] for i∉{(N+1)j:j=1,2,⋯,N};
    
  * ϕ has an upper bound estimation: sup{ϕ (x):x∈[0,1]}≤ 6(N+1)!.
 
It follows that

    |f̅(x)-ϕ(x)|≤ (x_i-x_i-1)^α≤ N^-2α,   if  x∈ [x_i-1,x_i],  for 
        i∉{(N+1)j:j=1,2,⋯,N}.

Define H_0=∪_i∈{(N+1)j:j=1,2,⋯,N} [x_i-1,x_i],  then

    |f̅ (x)-ϕ (x)|≤ N^-2α,     for any  x∈ [0,1]\ H_0,
    
by the fact that f̅∈Lip(1,α,d) and points in X are equispaced. By  μ(H_0)≤ Nδ, it follows that 


    f̅-ϕ_L^1([0,1])   =∫_H_0|f̅ (x)-ϕ (x)|dx+∫_[0,1]\ H_0|f̅ (x)-ϕ (x)|dx
       ≤ Nδ  (2+6(N+1)!)+N^2  (N^-2α) N^-2≤ 2N^-2α,

where the last inequality comes from the fact δ is small enough satisfying

    Nδ  (2+6(N+1)!)≤ N^-2α.

Note that f- (ϕ-1)=f+1-ϕ=f̅-ϕ. Hence, ϕ-1∈NN (=1;=[2N,2N+1]) and f- (ϕ-1)≤ 2N^-2α. So, we finish the proof for the case d=1.

 Step 2: The case d>1.    


The main idea is to project the d-dimensional problem into a one-dimensional one and use the results proved above. For any N∈^+, let n=⌊ N^2/d⌋ and  δ be a sufficiently small positive number depending on N and d, and satisfying (<ref>). We will divide the d-dimensional cube into n^d small non-overlapping sub-cubes (see Figure <ref> for an illustration when d=3 and n=3), each of which is associated with a representative point, e.g., a vertex of the sub-cube. Due to the continuity, the target function f can be represented by their values at the representative points. We project these representatives to one-dimensional samples via a ReLU FNN ψ and construct a ReLU FNN ϕ̅ to fit them. Finally, the ReLU FNN ϕ on the d-dimensional space approximating f can be constructed by ϕ=ϕ̅∘ψ. The precise construction can be found below.

By Lemma <ref>, there exists ψ_0∈NN (=1;=[2n]) such that

    
  * ψ_0 (1)=n-1, and ψ_0 (in)=ψ_0 (i+1n-δ)=i for i=0,1,⋯,n-1;
    
  * ψ_0 is linear between any two adjacent points of {in: i=0,1,⋯,n}∪{in-δ:i=1,2,⋯,n}.



Define the projection map[The idea constructing such ψ comes from the binary representation.] ψ by

    ψ (x)=∑_i=1^d1n^iψ_0 (x_i),  for x= [x_1,x_2,⋯,x_d]^T∈ [0,1]^d.

Note that ψ∈NN (=1;=[2dn]). Given f∈Lip (ν,α,d), then |f (x)|≤√(d) for x∈ [0,1]^d since f(0)=0, ν=1, and α∈(0, 1]. Define f̅=f+√(d), then 0≤f̅ (x) ≤ 2√(d) for x∈ [0,1]^d. Hence, we have

    {(∑_i=1^d θ_in^i, f̅ (θn)):θ= [θ_1,θ_2,⋯,θ_d]^T∈{0,1,⋯,n-1}^d }∪{ (1,0)}

as a set of n^d+1 samples of a one-dimensional function.  By Lemma <ref>, ∃ ϕ̅∈NN (=1;=[2⌈ n^d/2⌉, 2⌈ n^d/2⌉+1]) such that 

    ϕ̅(∑_i=1^d θ_in^i)=f̅ (θn),  for θ= [θ_1,θ_2,⋯,θ_d]^T∈{0,1,⋯,n-1}^d,

and 

    sup_t∈ [0,1]  |ϕ̅ (t)|≤ 6√(d)(⌈ n^d/2⌉+1)!.

Since the range of ψ on [0,1]^d is a subset of [0,1], ∃ ϕ∈NN (=d;=[2nd,2⌈ n^d/2⌉,2⌈ n^d/2⌉+1])
defined via ϕ(x)=ϕ̅ψ(x)  for x∈ [0,1]^d
such that

    sup_x∈ [0,1]^d  |ϕ (x)|≤ 6√(d)(⌈ n^d/2⌉+1)!.

Define H_1=∪_j=1^d{x= [x_1,x_2,⋯,x_d]^T∈ [0,1]^d:x_j∈∪_i=1^n[in-δ,in]}, which separates the d-dimensional cube into n^d important sub-cubes as illustrated in Figure <ref>. To index these d-dimensional smaller sub-cubes, define Q_θ={x= [x_1,x_2,⋯,x_d]^T∈ [0,1]^d:x_i∈[θ_in,θ_i+1n-δ],  i=1,2,⋯,d}
for each d-dimensional index θ= [θ_1,θ_2,⋯,θ_d]^T∈{0,1,⋯,n-1}^d.  
 By (<ref>), (<ref>), and the definition of ψ_0, for any x= [x_1,x_2,⋯,x_d]^T∈ Q_θ, we have ϕ (x)=ϕ̅ (ψ (x))
=ϕ̅ (∑_i=1^d 1n^iψ_0 (x_i))
=ϕ̅ (∑_i=1^d 1n^iθ_i)
=f̅ (θn). Then 

    |f̅ (x)-ϕ (x)|=|f̅ (x)-f̅ (θn)|≤  (√(d)/n)^α,   for any x∈ Q_θ.




Because μ (H_1)≤ dnδ, [0,1]^d=∪_θ∈{0,1,⋯,n-1}^dQ_θ∪ H_1, (<ref>), and (<ref>), we have

    f̅-ϕ_L^1([0,1]^d)   =∫ _H_1 |f̅-ϕ| dx+∫ _[0,1]^d\ H_1 |f̅-ϕ| dx
       ≤μ (H_1) (2√(d)+6√(d)(⌈ n^d/2⌉+1)!)
    +∑_θ∈{0,1,⋯,n-1}^d∫_Q_θ|f̅-ϕ| dx
       ≤ 2nδ d√(d)(1+3(⌈ n^d/2⌉+1)!)+∑_θ∈{0,1,⋯,n-1}^d (√(d)/n)^αμ (Q_θ)
       ≤  2d^α/2n^-α,

where the last inequality comes from the fact that δ is small enough such that




    2nδ d√(d)(1+3(⌈ n^d/2⌉+1)!)≤ d^α/2n^-α.


Note that f- (ϕ-√(d))=f̅-ϕ. Hence, ϕ-√(d)∈NN (=d;=[2nd,2⌈ n^d/2⌉,2⌈ n^d/2⌉+1]) and f- (ϕ-√(d))_L^1([0,1]^d)≤ 2d^α/2n^-α. Since n=⌊ N^2/d⌋, we have ⌈ n^d/2⌉≤ N+1. Therefore, 

    ϕ-√(d)   ∈NN (=d;=[2d⌊ N^2/d⌋,2N+2,2N+3])
    
and 

    f- (ϕ-√(d))_L^1([0,1]^d)≤ 2d^α/2n^-α= 2d^α/2⌊ N^2/d⌋^-α≤ 2d^α/2 (N^2/d/2)^-α=2 (2√(d))^α N^-2α/d,
 where the second inequality comes from the fact ⌊ x⌋≥x2 for x∈ [1,∞). So, we finish the proof when d>1.




Theorem <ref> shows that for f∈Lip(ν,α,d) ReLU FNNs with two or three function compositions can achieve the approximation rate (ν N^-2α/d). 
Following the same proof as in Theorem <ref>, we can show that:



∀ m,n∈^+, the closure of NN(=1;=[2m,2n+1]) contains CPL(mn+1) in the sense of L^1-norm.


An immediate implication of Corollary <ref> is that, for any function f on [0,1], if f can be approximated via one-hidden-layer ReLU FNNs with an approximation rate (N^-η) for any η>0, the rate can be improved to (N^-2η) via one more  composition. 




 §.§ Asymptotic Unachievable  Approximation Rate


In Section <ref>, we have analyzed the approximation capacity of ReLU FNNs in the nonlinear approximation for general continuous functions by construction. In this section, we will show that the construction in Section <ref> is essentially and asymptotically tight via showing the approximation lower bound in Theorem <ref> below.




    
    For any L∈^+, ρ>0, and C>0, ∃ f∈Lip(ν,α,d)  with α∈(0,1],  for all N_0>0, there exists N≥ N_0 such that 
    
    inf_ϕ∈NN(=d;≤ N;≤ L)ϕ-f_L^∞([0,1]^d)≥ Cν N^-(2α/d+ρ).




    
        
        ∀ L∈^+, ρ>0, and C>0, there exists N_0=N_0(L,ρ,C)>0 and f∈Lip(ν,α,d)  with α∈(0,1] such that for any integer N≥ N_0
        
    inf_ϕ∈NN(=d;≤ N;≤ L)ϕ(x)-f(x)_L^∞([0,1]^d)≥ Cν N^-(2α/d+ρ).

    


The proof of Theorem <ref> relies on the  nearly-tight VC-dimension bounds of ReLU FNNs given in <cit.> and is similar to Theorem 4 of <cit.>. Hence, we only sketch out its proof and a complete proof can be found in <cit.>. 


We will prove this theorem by contradiction.
Assuming that Theorem <ref> is not true, we can show the following claim, which will lead to a contradiction in the end.

    
    There exist
     L∈^+,  ρ>0, and C>0, ∀ f∈Lip (ν,α,d) with α∈ (0,1],  then ∃ N_0>0, for all N≥ N_0,  there exists
    ϕ∈NN(=d;≤ N;≤ L) such that
    
    f-ϕ_L^∞([0,1]^d)≤ Cν N^-(2α/d+ρ).


If this claim is true, then we have a better approximation rate. So we need to disproof this claim in order to prove Theorem <ref>.

Without loss of generality, we assume ν=1; in the case of ν≠ 1, the proof is similar by rescaling f∈Lip(ν,α,d) and FNNs with ν. Let us denote the VC dimension of a function set ℱ by VCDim (ℱ). By <cit.>, there exists C_1>0 such that

    VCDim(NN (=d;≤ N; ≤ L))
       ≤ C_1 ( (LN+d+2)(N+1))Lln ( (LN+d+2)(N+1)) b_u,

which comes from the fact the number of parameter of a ReLU FNN in (=d;≤ N;≤ L) is less than (LN+d+2)(N+1).


    
        
        There exist L∈^+, ρ>0, and C_2>0 such that ∀ N_0>0 and f∈Lip (1,α,d),  
        ∃  N≥ N_0 and ϕ∈NN (=d;≤ N;≤ L) such that
        
    f-ϕ_L^∞([0,1]^d)≤ C_2 N^-(2α/d+ρ).

    


One can estimate a lower bound of 

    VCDim(NN (=d;≤ N;≤ L))

using Claim <ref>,

More precisely, we construct {f_β:β∈ℬ}⊆Lip(1,α,d) shattering ⌊ N^2/d+ρ/(2α)⌋ ^d points, ℬ is an set, then we get corresponding {ϕ_β:β∈ℬ} by Claim <ref>, which can shatter ⌊ N^2/d+ρ/(2α)⌋ ^d points. 
and this lower bound can be b_ℓ⌊ N^2/d+ρ/(2α)⌋ ^d, which is asymptotically larger than

    b_u:=C_1 ( (LN+d+2)(N+1))Lln( (LN+d+2)(N+1))=(N^2ln N),

 leading to a contradiction that disproves the assumption that “Theorem <ref> is not true". 


Theorem <ref> shows that the N-term approximation rate via two or three-hidden-layer ReLU FNNs can achieve (N^-2α/d), while Theorem <ref> shows that the rate cannot be improved to (N^-(2α/d+ρ)) for any ρ>0. 
It was conjectured in the literature that function compositions can improve the approximation capacity exponentially. For general continuous functions, Theorem <ref> shows that this conjecture is not true, i.e., if the depth of the composition is L=(1), the approximation rate cannot be better than (N^-2α/d), not to mention (N^-Lα /d), which implies that adding one more layer cannot improve the approximation rate when N is large and L>2.

Following the same proof as in Theorem <ref>, we have the following corollary, which shows that the result in Corollary <ref> cannot be improved. 



∀ ρ>0, C>0 and L∈^+, ∃  N_0(ρ,C,L) such that for any integer N≥ N_0, CPL(CN^2+ρ) is not contained in the closure of NN(=1;≤ N;≤ L) in the sense of L^∞-norm.





 §.§ Approximation and Computation Efficiency in Parallel Computing



In this section, we will discuss the efficiency of the N-term approximation via ReLU FNNs in parallel computing. This is of more practical interest than the optimal approximation rate purely based on the number of parameters of the nonlinear approximation since it is impractical to use FNNs without parallel computing in real applications. Without loss of generality, we assume ν=1, N≫ 1, and d≫ 1. 

Let us summarize standard statistics of the time and memory complexity in parallel computing <cit.> in one training iteration of ReLU FNNs with (N) width and (L) depth using m computing cores and (1) training data samples per iteration. Let T_s(N,L,m) and T_d(N,L,m) denote the time complexity in shared and distributed memory parallel computing, respectively. Similarly, M_s(N,L,m) and M_d(N,L,m) are the memory complexity for shared and distributed memory, respectively. M_s(N,L,m) is the total memory requirement; while M_d(N,L,m) is the memory requirement per computing core. Then


    T_s(N,L,m)={[ 𝒪(L(N^2/m+lnmN)),       m∈ [1,N^2],;        𝒪(Lln N) ,       m∈ (N^2,∞);;                   ].


    T_d(N,L,m)={[ 𝒪(L(N^2/m+t_sln m+t_w N√(m)ln m)),                        m∈ [1,N^2],;                         𝒪(Lln N) ,                        m∈ (N^2,∞);;                                    ].


    M_s(N,L,m)=𝒪(LN^2),  for all  m∈^+ ;

and

    M_d(N,L,m)=
    𝒪(LN^2/m+1),  for all  m∈^+ ,

where t_s and t_w are the “start-up time” and “per-word transfer time” in the data communication between different computing cores, respectively (see <cit.> for a detailed introduction). 
 
In real applications, a most frequently asked question would be: given a target function f∈Lip(ν,α,d), a target approximation accuracy ϵ, and a certain amount of computational resources, e.g., m computer processors, assuming the computer memory is enough, what is a good choice of FNN architecture we should use to reduce the running time of our computers? Certainly, the answer depends on the number of processors m and ideally we hope to increase m by a factor of r to reduce the time (and memory in the distributed environment) complexity by the same factor r, which is the scalability of parallel computing. 
 
 We answer the question raised just above using FNN architectures that almost have a uniform width since the optimal approximation theory of very deep FNNs <cit.> and this manuscript both utilize a nearly uniform width. Combining the theory in <cit.> and ours, we summarize several statistics of ReLU FNNs in parallel computing in Table <ref> and <ref> when FNNs nearly have the same approximation accuracy. For shared memory parallel computing, from Table <ref> we see that: if computing resources are enough, shallower FNNs with (1) hidden layers require less and even exponentially less running time than very deep FNNs; if computing resources are limited, shallower FNNs might not be applicable or are slower, and hence very deep FNNs are a good choice. For distributed memory parallel computing, the conclusion is almost the same by Table <ref>, except that the memory limit is not an issue for shallower FNNs if the number of processors is large enough. In sum, if the approximation rate of very deep FNNs is not exponentially better than shallower FNNs, very deep FNNs are less efficient than shallower FNNs theoretically if computing resources are enough. 
 
  






§ NUMERICAL EXPERIMENTS



In this section, we provide two sets of numerical experiments to compare different ReLU FNNs using shared memory GPU parallel computing. The numerical results for distributed memory parallel computing would be similar. All numerical tests were conducted using Tensorflow and an NVIDIA P6000 GPU with 3840 CUDA parallel processing cores.

Since it is difficult to generate target functions f∈Lip(ν,α,d) with fixed ν and α, we cannot directly verify the nonlinear approximation rate, but we are able to observe numerical evidence close to our theoretical conclusions. Furthermore, we are able to verify the running time estimates in Section <ref> and show that, to achieve the same theoretical approximation rate, shallow FNNs are more efficient than very deep FNNs.

In our numerical tests, we generate 50 random smooth functions as our target functions using the algorithm in <cit.> with a wavelength parameter λ=0.1 and an amplitude parameter √((2/λ)) therein. These target functions are uniformly sampled with 20000 ordered points {x_i} in [0,1] to form a data set. 
The training data set consists of samples with odd indices i's, while the test data set consists of samples with even indices i's. The loss function is defined as the mean square error between the target function and the FNN approximant evaluated on training sample points. The ADAM algorithm <cit.> with a decreasing learning rate from 0.005 to 0.0005, a batch size 10000, and a maximum number of epochs 20000, is applied to minimize the mean square error. The minimization is randomly initialized by the “normal initialization method"[See <https://medium.com/prateekvishnu/xavier-and-he-normal-he-et-al-initialization-8e3d7a087528>.].
The test error is defined as the mean square error evaluated on test sample points. The training and test data sets are essentially the same in our numerical test since we aim at studying the approximation power of FNNs instead of the generalization capacity of FNNs. Note that due to the high non-convexity of the optimization problem, there might be chances such that the minimizers we found are bad local minimizers. Hence, we compute the average test error of the best 40 tests among the total 50 tests of each architecture. 

To observe numerical phenomena in terms of N-term nonlinear approximation,  in the first set of numerical experiments, we use two types of FNNs to obtain approximants: the first type has L=(1) layers with different sizes of width N; the second type has a fixed width N=12 with different numbers of layers L. Numerical results are summarized in Table <ref>. To observe numerical phenomena in terms of the number of parameters in FNNs, in the second set of numerical experiments, we use FNNs with the same number of parameters but different sizes of width N and different numbers of layers L. Numerical results are summarized in Table <ref>. 



By the last columns of Table <ref>, we verified that as long as the number of computing cores m is larger than or equal to N^2, the running time per iteration of FNNs with (1) layers is (ln N), while the running time per iteration of FNNs with (N) layers and (1) width is (N). By the last columns of Table <ref>, we see that when the number of parameters is the same, very deep FNNs requires much more running time per iteration than shallower FNNs and the difference becomes more significant when the number of parameters increases. Hence, very deep FNNs are much less efficient than shallower FNNs in parallel computing. 

Besides, by Table <ref> and <ref>, the test error of very deep FNNs cannot be improved if the depth is increased and the error even becomes larger when depth is larger. However, when the number of layers is fixed, increasing width can reduce the test error. More quantitatively, we define the improvement ratio of an FNN with width N and depth L in Table <ref> as the ratio of the test error of an FNN in NN(=1;=[N/2]^L) (or NN(=1;=[N]^L/2)) over the test error of the current FNN in NN(=1;=[N]^L). Similarly, the improvement ratio of an FNN with a number of parameters W in Table <ref> is defined as the ratio of the test error of an FNN with the same type of architecture and a number of parameters W/2 over the test error of the current FNN. According to the improvement ratio in Table <ref> and <ref>, when L=(1), the numerical approximation rate in terms of N is in a range between 2 to 4. Due to the high non-convexity of the deep learning optimization and the difficulty to generate target functions of the same class with a fixed order α and constant ν, we cannot accurately verify the approximation rate. But the statistics of the improvement ratio can roughly reflect the approximation rate and the numerical results stand in line with our theoretical analysis. 







§ CONCLUSIONS


 We studied the approximation and computation efficiency of nonlinear approximation via compositions, especially when the dictionary 𝒟_L consists of deep ReLU FNNs with width N and depth L. Our main goal is to quantify the best N-term approximation rate ϵ_L,f(N) for the dictionary 𝒟_L when f is a Hölder continuous function. This topic is related to several existing approximation theories in the literature, but they cannot be applied to answer our problem. By introducing new analysis techniques that are merely based on the FNN structure instead of traditional approximation basis functions used in existing work, we have established a new theory to address our problem. 
 
 In particular, for any function f on [0,1], regardless of its smoothness and even the continuity, if f can be approximated using a dictionary when L=1 with the best N-term approximation rate ϵ_L,f=(N^-η), we showed that dictionaries with L=2 can improve the best N-term approximation rate to ϵ_L,f=(N^-2η). We also showed that for Hölder continuous functions of order α on [0,1]^d, the application of a dictionary with L=3 in nonlinear approximation can achieve an essentially tight best N-term approximation rate ϵ_L,f=( N^-2α/d), and increasing L further cannot improve the approximation rate in N. Finally, we showed that dictionaries consisting of wide FNNs with a few hidden layers are more attractive in terms of computational efficiency than dictionaries with narrow and very deep FNNs for approximating Hölder continuous functions if the number of computer cores is larger than N in parallel computing. 

Our results were based on the L^1-norm in the analysis of constructive approximations in Section <ref>; while we used L^∞-norm in the unachievable approximation rate in Section <ref>. In fact, we can define a new norm that is weaker than the L^∞-norm and stronger than the L^1-norm such that all theorems in this paper hold in the same norm. The analysis based on the new norm can be found in <cit.> and our future work, and it shows that the approximation rate in this paper is tight in the same norm. Finally, although the current result is only valid for Hölder continuous functions, it is easy to generalize it to general continuous functions by introducing the moduli of continuity, which is also left as future work.



 
 
















§ ACKNOWLEDGMENTS

H. Yang thanks the Department of Mathematics at the National University of Singapore for the startup grant, the Ministry of Education in Singapore for the grant MOE2018-T2-2-147, NVIDIA Corporation for the donation of the Titan Xp GPU, and NSCC <cit.> for extra GPU resource. We would like to acknowledge the editor and reviewers for their help in improving this manuscript.

apacite


