OPTIMAL RATE AND CURSE OF DIMENSIONALITY
Year: 2019

Abstract:
Deep learning has shown high performances in various types of tasks from visual recognition to natural language processing, which indicates superior flexibility and adaptivity of deep learning. To understand this phenomenon theoretically, we develop a new approximation and estimation error analysis of deep learning with the ReLU activation for functions in a Besov space and its variant with mixed smoothness. The Besov space is a considerably general function space including the Hölder space and Sobolev space, and especially can capture spatial inhomogeneity of smoothness. Throughout the analysis in the Besov space, it is shown that deep learning can achieve the minimax optimal rate and outperform any nonadaptive (linear) estimator such as kernel ridge regression, which shows that deep learning has higher adaptivity to the spatial inhomogeneity of the target function than other estimators such as linear ones. In addition to this, it is shown that deep learning can avoid the curse of dimensionality if the target function is in a mixed smooth Besov space. We also show that the dependency of the convergence rate on the dimensionality is tight due to its minimax optimality. These results support high adaptivity of deep learning and its superior ability as a feature extractor.

heading: 1 INTRODUCTION
text: Deep learning has shown great success in several applications such as computer vision and natural language processing. As its application range is getting wider, theoretical analysis to reveal the reason why deep learning works so well is also gathering much attention. To understand deep learning theoretically, several studies have been developed from several aspects such as approximation theory and statistical learning theory. A remarkable property of neural network is that it has universal approximation capability even if there is only one hidden layer (Cybenko, 1989; Hornik, 1991; Sonoda & Murata, 2017). Thanks to this property, deep and shallow neural networks can approximate any function with any precision (of course, the meaning of the terminology “any” must be rigorously defined like “any function in L1(R)”). A natural question coming next is its expressive power. It is shown that the expressive power of deep neural network grows exponentially against the number of layers (Montufar et al., 2014; Bianchini & Scarselli, 2014; Cohen et al., 2016; Cohen & Shashua, 2016; Poole et al., 2016) where the “expressive power” is defined by several ways.
The expressive power of neural network can be analyzed more precisely by specifying the target function’s property such as smoothness. Barron (1993; 1994) developed an approximation theory for functions having limited “capacity” that is measured by integrability of their Fourier transform. An interesting point of the analysis is that the approximation error is not affected by the dimensionality of the input. This observation matches the experimental observations that deep learning is quite effective also in high dimensional situations. Another typical approach is to analyze function spaces with smoothness conditions such as the Hölder space. In particular, deep neural network with the ReLU activation (Nair & Hinton, 2010; Glorot et al., 2011) has been extensively studied recently from the view point of its expressive power and its generalization error. For example, Yarotsky (2017) derived the approximation error of the deep network with the ReLU activation for
functions in the Hölder space. Schmidt-Hieber (2018) evaluated the estimation error of regularized least squares estimator performed by deep ReLU network based on this approximation error analysis in a nonparametric regression setting. Petersen & Voigtlaender (2017) generalized the analysis of Yarotsky (2017) to the class of piece-wise smooth functions. Imaizumi & Fukumizu (2018) utilized this analysis to derive the estimation error to estimate the piece-wise smooth function and concluded that deep leaning can outperform linear estimators in that setting. Although these error analyses are standard from a nonparametric statistics view point and the derived rates are known to be (near) minimax optimal, the analysis is rather limited because they are given mainly based on the Hölder space. However, there are several other function spaces such as the Sobolev space and the space of finite total variations. A comprehensive analysis to deal with such function classes from a unified view point is required.
In this paper, we give generalization error bounds of deep ReLU networks for a Besov space and its variant with mixed smoothness, which includes the Hölder space, the Sobolev space, and the function class with total variation as special cases. By doing so, (i) we show that deep learning1 achieves the minimax optimal rate on the Besov space and notably it outperforms any linear estimator such as the kernel ridge regression, and (ii) we show that deep learning can avoid the curse of dimensionality on the mixed smooth Besov space and achieves the minimax optimal rate. As related work, Mhaskar & Micchelli (1992); Mhaskar (1993); Chui et al. (1994); Mhaskar (1996); Pinkus (1999) also developed an approximation error analysis which essentially leads to analyses for Besov spaces. However, the ReLU activation is basically excluded and comprehensive analyses for the Besov space have not been given. As a summary, the contribution of this paper is listed as follows:
(i) To investigate adaptivity of deep learning, we give an explicit form of approximation and estimation error bounds for deep learning with the ReLU activation where the target functions are in the Besov spaces (Bsp,q) for s > 0 and 0 < p, q ≤ ∞ with s > d(1/p− 1/r)+ where Lr-norm is used for error evaluation. In particular, deep learning outperforms any linear estimator such as kernel ridge regression if the target function has highly spatial inhomogeneity of its smoothness. See Table 1 for the overview.
(ii) To investigate the effect of dimensionality, we analyze approximation and estimation problems in so-called the mixed smooth Besov space by ReLU neural network. It is shown that deep learning with the ReLU activation can ease the curse of dimensionality and achieve the near minimax optimal rate. The theory is developed on the basis of the sparse grid technique (Smolyak, 1963). See Table 2 for the overview.

heading: 2 SET UP OF FUNCTION SPACES
text: In this section, we define the function classes for which we develop error bounds. In particular, we define the Besov space and its variant with mixed smoothness. The typical settings in statistical learning theory is to estimate a function with a smoothness condition. There are several ways to characterize “smoothness.” Here, we summarize the definitions of representative functional spaces that are appropriate to define the smoothness assumption.
Let Ω ⊂ Rd be a domain of the functions. Throughout this paper, we employ Ω = [0, 1]d. For a function f : Ω → R, let ‖f‖p := ‖f‖Lp(Ω) := ( ∫ Ω |f |pdx)1/p for 0 < p < ∞. For p =
1In this paper, we mainly consider a regularized empirical risk minimization given in Eq. (7) as an estimation procedure, and say “deep learning” to indicate it.
∞, we define ‖f‖∞ := ‖f‖L∞(Ω) := supx∈Ω |f(x)|. For α ∈ Rd, let |α| = ∑d j=1 |αj |. Let C0(Ω) be the set of continuous functions equipped with L∞-norm: C0(Ω) := {f : Ω → R | f is continuous and ‖f‖∞ <∞} 2. For α ∈ Zd+, we denote by Dαf(x) = ∂|α|f ∂α1x1...∂ αdxd (x) 3. Definition 1 (Hölder space (Cβ(Ω))). Let β > 0 with β 6∈ N be the smoothness parameter. For an m times differentiable function f : Rd → R, let the norm of the Hölder space Cβ(Ω) be ‖f‖Cβ := max|α|≤m
∥∥Dαf‖∞ + max|α|=m supx,y∈Ω |Dαf(x)−Dαf(y)||x−y|β−m , where m = bβc (the largest integer less than β). Then, (β-)Hölder space Cβ(Ω) is defined as Cβ(Ω) = {f | ‖f‖Cβ <∞}.
The parameter β > 0 controls the “smoothness” of the function. Along with the Hölder space, the Sobolev space is also important. Definition 2 (Sobolev space (W kp (Ω))). Sobolev space (W kp (Ω)) with a regularity parameter k ∈ N and a parameter 1 ≤ p ≤ ∞ is a set of functions such that the Sobolev norm ‖f‖Wkp := ( ∑ |α|≤k ‖Dαf‖pp) 1 p is finite (where the derivative D is taken in the weak sense).
There are some ways to define a Sobolev space with fractional order, one of which will be defined by using the notion of interpolation space (DeVore, 1998; Adams & Fournier, 2003), but we don’t pursue this direction here. Finally, we introduce Besov space which further generalizes the definition of the Sobolev space. To define the Besov space, we introduce the modulus of smoothness. Definition 3. For a function f ∈ Lp(Ω) for some p ∈ (0,∞], the r-th modulus of smoothness of f is defined by
wr,p(f, t) = sup h∈Rd:‖h‖2≤t
‖∆rh(f)‖p,
where ∆rh(f)(x) =
{∑r j=0 ( r j ) (−1)r−jf(x+ jh) (x ∈ Ω, x+ rh ∈ Ω),
0 (otherwise).
Based on the modulus of smoothness, the Besov space is defined as in the following definition. Definition 4 (Besov space (Bαp,q(Ω))). For 0 < p, q ≤ ∞, α > 0, r := bαc + 1, let the seminorm | · |Bαp,q be
|f |Bαp,q :=
{(∫∞ 0 (t−αwr,p(f, t)) q dt t ) 1 q (q <∞),
supt>0 t −αwr,p(f, t) (q =∞).
The norm of the Besov space Bαp,q(Ω) can be defined by ‖f‖Bαp,q := ‖f‖p+ |f |Bαp,q , and B α p,q(Ω) = {f ∈ Lp(Ω) | ‖f‖Bαp,q <∞}. 2Since Ω = [0, 1]d in our setting, the boundedness automatically follows from the continuity. 3 We let N := {1, 2, 3, . . . }, Z+ := {0, 1, 2, 3, . . . }, Zd+ := {(z1, . . . , zd) | zi ∈ Z+}, R+ := {x ≥ 0 | x ∈ R}, and R++ := {x > 0 | x ∈ R}.
Note that p, q < 1 is also allowed. In that setting, the Besov space is no longer a Banach space but a quasi-Banach space. The Besov space plays an important role in several fields such as nonparametric statistical inference (Kerkyacharian & Picard, 1992; Donoho et al., 1996; Donoho & Johnstone, 1998; Giné & Nickl, 2015) and approximation theory (Temlyakov, 1993a). These spaces are closely related to each other as follows (Triebel, 1983):
• For m ∈ N, Bmp,1(Ω) ↪→Wmp (Ω) ↪→ Bmp,∞(Ω), and Bm2,2(Ω) = Wm2 (Ω) 4. • For 0 < s <∞ and s 6∈ N, Cs(Ω) = Bs∞,∞(Ω).
• For 0 < s, p, q, r ≤ ∞ with s > δ := d(1/p− 1/r)+, it holds that Bsp,q(Ω) ↪→ Bs−δr,q (Ω). In particular, under the same condition, from the definition of ‖ · ‖Bsp,q , it holds that
Bsp,q(Ω) ↪→ Lr(Ω). (1)
• For 0 < s, p, q ≤ ∞, if s > d/p, then
Bsp,q(Ω) ↪→ C0(Ω). (2)
Hence, if the smoothness parameter satisfies s > d/p, then it is continuously embedded in the set of the continuous functions. However, if s < d/p, then the elements in the space are no longer continuous. Moreover, it is known that B11,1([0, 1]) is included in the space of bounded total variation and B11,∞([0, 1]) includes it (Peetre & Dept, 1976). Hence, the Besov space also allows spatially inhomogeneous smoothness with spikes and jumps; which makes difference between linear estimators and deep learning (see Sec. 4.1).
It is known that the minimax rate to estimate fo is lower bounded by n−2s/(2s+d), (Kerkyacharian & Picard, 1992; Donoho et al., 1996; Donoho & Johnstone, 1998; Giné & Nickl, 2015). We see that the curse of dimensionality is unavoidable as long as we consider the Besov space. This is an undesirable property because we easily encounter high dimensional data in several machine learning problems. Hence, we need another condition to derive approximation and estimation error bounds that are not heavily affected by the dimensionality. To do so, we introduce the notion of mixed smoothness. The Besov space with mixed smoothness is defined as follows (Schmeisser, 1987; Sickel & Ullrich, 2009). To define the space, we define the coordinate difference operator as
∆r,ih (f)(x) = ∆ r h(f(x1, . . . , xi−1, ·, xi+1, . . . , xd))(xi)
for f : Rd → R, h ∈ R+, i ∈ [d], and r ≥ 1. By applying this difference operator to each coordinate in a subset e ⊂ {1, . . . , d} of coordinates, the mixed differential operator for a step length h ∈ Rd is defined as
∆r,eh (f) = (∏ i∈e ∆ r,i hi ) (f), ∆r,∅h (f) = f.
Then, the mixed modulus of smoothness is defined as
wer,p(f, t) := sup|hi|≤ti,i∈e ‖∆ r,e h (f)‖p
for t ∈ Rd+ and 0 < p ≤ ∞. This quantity measures how the function f is “rough” in a coordinatewise manner. In contrast to the Besov space (Def. 3), the differentiation is taken over all coordinates in e, which imposes a coordinate-wise smoothness. Letting 0 < p, q ≤ ∞, α ∈ Rd++ and ri := bαic+ 1, the semi-norm | · |MBα,ep,q based on the mixed smoothness is defined by
|f |MBα,ep,q :=  {∫ Ω [( ∏ i∈e t −αi i )w e r,p(f, t)] q dt∏ i∈e ti }1/q (0 < q <∞),
supt∈Ω( ∏ i∈e t −αi i )w e r,p(f, t) (q =∞).
By summing up the semi-norm over the choice of e, the (quasi-)norm of the mixed smooth Besov space (abbreviated to m-Besov space) is defined by
‖f‖MBαp,q := ‖f‖p + ∑
e⊂{1,...,d}
|f |MBα,ep,q ,
4The notation ↪→ means a continuous embedding; that is, X ↪→ Y for two normed spaces X,Y means X is continuously embedded in Y .
and thus MBαp,q(Ω) := {f ∈ Lp(Ω) | ‖f‖MBαp,q < ∞} where 0 < p, q ≤ ∞ and α ∈ R d ++. In this paper, we assume that α1 = · · · = αd. With a slight abuse of notation, we also use the notation MBαp,q for α > 0 to indicate MB (α,...,α) p,q .
It is known that, when p = q, the m-Besov space is characterized as a tensor product space of Bsp,q([0, 1]) (Sickel & Ullrich, 2009). The m-Besov space includes several important models considered in the literature of statistical learning, e.g., the additive model (Meier et al., 2009) and the tensor model (Signoretto et al., 2010). It is known that an appropriate estimator in these models can avoid curse of dimensionality (Meier et al., 2009; Raskutti et al., 2012; Kanagawa et al., 2016; Suzuki et al., 2016). What we will show in this paper supports that this fact is also applied to deep learning from a unifying viewpoint.
The difference between the (normal) Besov space and the m-Besov space can be informally explained as follows. For regularity condition αi ≤ 2 (i = 1, 2), the m-Besov space consists of functions for which the following derivatives are “bounded”:
∂f ∂x1 , ∂f ∂x2 , ∂2f ∂x21 , ∂2f ∂x22 , ∂2f ∂x1∂x2 , ∂3f ∂x1∂x22 , ∂3f ∂x21∂x2 , ∂4f ∂x21∂x 2 2 .
That is, the “max” of the orders of derivatives over coordinates needs to be bounded by 2. On the other hand, the Besov space only ensures the boundedness of the following derivatives:
∂f ∂x1 , ∂f ∂x2 , ∂2f ∂x21 , ∂2f ∂x22 , ∂2f ∂x1∂x2 ,
where the “sum” of the orders needs to be bounded by 2. This difference directly affects the rate of convergence of approximation accuracy. Further details about this space and related topics can be found in a comprehensive survey (Dũng et al., 2016).
Relation to Barron class. Barron (1991; 1993; 1994) showed that, if the Fourier transform of a function f satisfies some integrability condition, then we may avoid curse of dimensionality for estimating neural networks with sigmoidal activation functions. The integrability condition is given by ∫ Cd ‖ω‖|f̂(ω)|dω < ∞, where f̂ is the Fourier transform of a function f . We call the class of functions satisfying this condition Barron class. A similar function class is analyzed by Klusowski & Barron (2016) too. We cannot compare directly the m-Besov space with the Barron class, but they are closely related. Indeed, if p = q = 2 and s = α1 = · · · = αd, then m-Besov space MBs2,2(Ω) is equivalent to the tensor product of Sobolev space (Sickel & Ullrich, 2011) which consists of functions f : Ω → R satisfying ∫ Cd ∏d i=1(1 + |ωi|2)s|f̂(ω)|2dω < ∞. Therefore, our analysis gives a (similar but) different characterization of conditions to avoid curse of dimensionality.

heading: 3 APPROXIMATION ERROR ANALYSIS
text: In this section, we evaluate how well the functions in the Besov and m-Besov spaces can be approximated by neural networks with the ReLU activation. Let us denote the ReLU activation by η(x) = max{x, 0} (x ∈ R), and for a vector x, η(x) is operated in an element-wise manner. Define the neural network with height L, width W , sparsity constraint S and norm constraint B as
Φ(L,W,S,B) := {(W(L)η(·) + b(L)) ◦ · · · ◦ (W(1)x+ b(1)) | W(L) ∈ R1×W , b(L) ∈ R, W(1) ∈ RW×d, b(1) ∈ RW , W(`) ∈ RW×W , b(`) ∈ RW (1 < ` < L),∑L
`=1 (‖W(`)‖0 + ‖b(`)‖0) ≤ S,max ` ‖W(`)‖∞ ∨ ‖b(`)‖∞ ≤ B},
where ‖ · ‖0 is the `0-norm of the matrix (the number of non-zero elements of the matrix) and ‖ · ‖∞ is the `∞-norm of the matrix (maximum of the absolute values of the elements). We want to evaluate how large L,W, S,B should be to approximate fo ∈ MBαp,q(Ω) by an element f ∈ Φ(L,W, S,B) with precision > 0 measured by Lr-norm: minf∈Φ ‖f − fo‖r ≤ .

heading: 3.1 APPROXIMATION ERROR ANALYSIS FOR BESOV SPACES
text: Here, we show how the neural network can approximate a function in the Besov space which is useful to derive the generalization error of deep learning. Although its derivation is rather standard
as considered in Chui et al. (1994); Bölcskei et al. (2017), it should be worth noting that the bound derived here cannot be attained by any non-adaptive method and the generalization error based on the analysis is also unattainable by any linear estimators including the kernel ridge regression. That explains the high adaptivity of deep neural network and how it outperforms usual linear methods such as kernel methods.
To show the approximation accuracy, a key step is to show that the ReLU neural network can approximate the cardinal B-spline with high accuracy. LetN (x) = 1 (x ∈ [0, 1]), 0 (otherwise), then the cardinal B-spline of order m is defined by taking m+ 1-times convolution of N :
Nm(x) = (N ∗N ∗ · · · ∗ N︸ ︷︷ ︸ m+ 1 times )(x),
where f ∗ g(x) := ∫ f(x − t)g(t)dt. It is known that Nm is a piece-wise polynomial of order m.
For k = (k1, . . . , kd) ∈ Zd+ and j = (j1, . . . , jd) ∈ Zd+, letMdk,j(x) = ∏d i=1Nm(2kixi−ji). Even
for k ∈ Z+, we also use the same notation to express Mdk,j(x) = ∏d i=1Nm(2kxi − ji). Here, k controls the spatial “resolution” and j specifies the location on which the basis is put. Basically, we approximate a function f in a Besov space by a super-position of Mmk,j(x), which is closely related to wavelet analysis (Mallat, 1999).
Mhaskar & Micchelli (1992); Chui et al. (1994) have shown the approximation ability of neural network for a function with bounded modulus of smoothness. However, the class of the activation functions in their analysis does not include ReLU but they dealt with activation functions satisfying the following conditions,
lim x→∞ η(x)/xk → 1, lim x→−∞ η(x)/xk = 0, ∃K > 1 s.t. |η(x)| ≤ K(1 + |x|)k (x ∈ R), (3a)
for k = 2 which excludes ReLU. Mhaskar (1993) analyzed deep neural network under the same setting but it restricts the smoothness parameter to s = k + 1. Mhaskar (1996) considered the Sobolev space Wmp with an infinitely many differentiable “bump” function which also excludes ReLU. However, approximating the cardinal B-spline by ReLU can be attained by appropriately using the technique developed by Yarotsky (2017) as in the following lemma. Lemma 1 (Approximation of cardinal B-spline basis by the ReLU activation). There exists a constant c(d,m) depending only on d and m such that, for all > 0, there exists a neural network M̌ ∈ Φ(L0,W0, S0, B0) withL0 := 3+2 ⌈ log2 ( 3d∨m
c(d,m)
) + 5 ⌉ dlog2(d ∨m)e,W0 := 6dm(m+2)+2d,
S0 := L0W 2 0 and B0 := 2(m+ 1) m that satisfies ‖Md0,0 − M̌‖L∞(Rd) ≤ , and M̌(x) = 0 for all x 6∈ [0,m+ 1]d.
The proof is in Appendix A. Based on this lemma, we can translate several B-spline approximation results into those of deep neural network approximation. In particular, combining this lemma and the B-spline interpolant representations of functions in Besov spaces (DeVore & Popov, 1988; DeVore et al., 1993; Dũng, 2011b), we obtain the optimal approximation error bound for deep neural networks. Here, let U(H) be the unit ball of a quasi-Banach space H, and for a set F of functions, define the worst case approximation error as
Rr(F ,H) := sup fo∈U(H) inf f∈F ‖fo − f‖Lr([0,1]d).
Proposition 1 (Approximation ability for Besov space). Suppose that 0 < p, q, r ≤ ∞ and 0 < s <∞ satisfy the following condition:
s > d(1/p− 1/r)+. (4) Assume that m ∈ N satisfies 0 < s < min(m,m − 1 + 1/p). Let δ = d(1/p − 1/r)+ and ν = (s− δ)/(2δ). For sufficiently large N ∈ N and = N−s/d−(ν−1+d−1)(d/p−s)+ log(N)−1, let
L = 3 + 2dlog2 ( 3d∨m
c(d,m)
) + 5edlog2(d ∨m)e, W = NW0,
S = [(L− 1)W 20 + 1]N, B = O(N (ν −1+d−1)(1∨(d/p−s)+)),
then it holds that Rr(Φ(L,W, S,B), B s p,q([0, 1] d)) . N−s/d.
Remark 1. By Eq. (1), the condition (4) indicates that fo ∈ Bsp,q satisfies fo ∈ Lr(Ω). If we set p = q = ∞ and r = ∞, then Bsp,q(Ω) = Cs(Ω) which yields the result by Yarotsky (2017) as a special case.
The proof is in Appendix B. According to the theorem, the approximation error N−s/d can be achieved by the settings L = O(log(N)), W = O(N) and S = O(N log(N)) for N ∈ N. The convergence rate N−s/d is controlled by the smoothness s and the dimensionality d; as the smoothness s goes up, we have better approximation error, and as the dimensionality d goes up, we have worse error. An interesting point is that the statement is valid even for p 6= r. In particular, the theorem also supports non-continuous regime (s < d/p) in which L∞-convergence does no longer hold but instead Lr-convergence is guaranteed under the condition s > d(1/p − 1/r)+. In that sense, the convergence of the approximation error is guaranteed in considerably general settings. Pinkus (1999) gave an explicit form of convergence when 1 ≤ p = r for the activation functions satisfying Eq. (3) which does not cover ReLU and an important setting p 6= r. Petrushev (1998) considered p = r = 2 and activation function with Eq. (3) where s is an integer such that s ≤ k+1+(d−1)/2. Chui et al. (1994) and Bölcskei et al. (2017) dealt with the smooth sigmoidal activation satisfying the condition (3) with k ≥ 2 or a “smoothed version” of the ReLU activation which excludes ReLU; but Bölcskei et al. (2017) presented a general strategy for neural-net approximation by using the notion of best M -term approximation. Mhaskar & Micchelli (1992) gives an approximation bound using the modulus of smoothness, but the smoothness s and the order of sigmoidal function k in (3) is tightly connected and fo is assumed to be continuous which excludes the situation s < d/p. On the other hand, the above proposition does not require such a tight connection and it explicitly gives the approximation bound for Besov spaces. Williamson & Bartlett (1992) derived a spline approximation error bound for an element in a Besov space when d = 1, but the derived bound is only O(N−s+(1/p−1/r)+) which is the rate of non-adaptive methods described below, and approximation by a ReLU activation network is not discussed. We may also use the analysis of Cohen et al. (2001) which is based on compactly supported wavelet bases, but the cardinal B-spline is easy to handle through quasi-interpolant representation as performed in the proof of Proposition 1.
It should be noted that the presented approximation accuracy bound is not trivial because it can not be achieved by a non-adaptive method. Actually, the bestN -term approximation error (Kolmorogov width) of the Besov space is lower bounded as
inf SN⊂Bsp,q sup f∈U(Bsp,q) inf f̌∈SN
‖f − f̌‖Lr(Ω) & N −s/d+(1/p−1/r)+ (1 < p < r ≤ 2, s > d(1/p− 1/r)),
N−s/d+1/p−1/2 (1 < p < 2 < r ≤ ∞, s > d/p), N−s/d (2 ≤ p < r ≤ ∞, s > d/2),
(5) if 1 < p < r ≤ ∞, 1 ≤ q < ∞ and 1 < s, where SN is any N -dimensional subspace of Bsp,q (Romanyuk, 2009; Myronyuk, 2016; Vybáral, 2008). That is, any linear/non-linear approximator with fixed N -bases does not achieve the approximation error N−s/d in some parameter settings such as 1 < p < 2 < r. On the other hand, adaptive methods including deep learning can improve the error rate up to N−s/d which is rate optimal (Dũng, 2011b). The difference is significant when p < r. This implies that deep neural network possesses high adaptivity to find which part of the function should be intensively approximated. In other words, deep neural network can properly extracts the feature of the input (which corresponds to construct an appropriate set of bases) to approximate the target function in the most efficient way.

heading: 3.2 APPROXIMATION ERROR ANALYSIS FOR M-BESOV SPACE
text: Here, we deal with m-Besov spaces instead of the ordinary Besov space. The next theorem gives the approximation error bound to approximate functions in the m-Besov spaces by deep neural network
models. Define Dk,d := ( 1 + d−1k )k ( 1 + kd−1 )d−1 . Then, we have the following theorem.
Theorem 1 (Approximation ability for m-Besov space). Suppose that 0 < p, q, r ≤ ∞ and s <∞ satisfies s > (1/p − 1/r)+. Assume that m ∈ N satisfies 0 < s < min(m,m − 1 + 1/p). Let δ = (1/p − 1/r)+ and ν = (s − δ)/(2δ). For any K ≥ 1, let K∗ = dK(1 + 2δα−δ )e. Then, for N = (2 + (1− 2−ν)−1)2KDK∗,d, if we set
L = 3 + 2 ⌈ log2 ( 3d∨m
c(d,m)
) + 5 + (s+ ( 1p − s)+ + 1)K ∗ + log([e(m+ 1)]d(1 +K∗)) ⌉ dlog2(d ∨m)e,
W = W0N, S = [(L− 1)W 20 + 1]N, B = O(N (ν −1+1)(1∨(1/p−s)+)), then it holds that (i) For p ≥ r, Rr(Φ(L,W,S,B),MBsp,q([0, 1]d)) . 2−KsD (1/min(r,1)−1/q)+ K,d , (6a)
(ii) For p < r, Rr(Φ(L,W,S,B),MBsp,q([0, 1] d)) .
{ 2−KsD
(1/r−1/q)+ K,d (r <∞),
2−KsD (1−1/q)+ K,d (r =∞).
(6b)
The proof is given in Appendix C.3. It holds that N ' 2KK(d−1), which implies 2−K ' N−1 logd−1(N) if N d (see also the discussion right after Theorem 5 in Appendix C.1 for more details of calculation). Therefore, when r q, the approximation error can be evaluated as O(N−s logs(d−1)(N)) for L = O(log(N)), W = O(N) and S = O(N log(N)) for N ∈ N in which the effect of dimensionality d is much milder than that of Proposition 1. This means that the curse of dimensionality is much eased in the mixed smooth space.
The obtained bound is far from obvious. Actually, it is better than any linear approximation methods as follows. Let the linear M -width introduced by Tikhomirov (1960) be λN (MBsp,q, L
r) := infLN supf∈U(MBsp,q) ‖f − LN (f)‖r, where the infimum is taken over all linear oprators LN with rank N from MBsp,q to L
r. The linear N -width of the m-Besov space has been extensively studied as in the following proposition (see Lemma 5.1 of Dũng (2011a), and Romanyuk (2001)). Proposition 2. Let 1 ≤ p, r ≤ ∞, 0 < q ≤ ∞ and s > (1/p− 1/r)+. Then we have the following asymptotic order of the linear width for the asymptotics N d: (a) For p ≥ r,
λN (MB s p,q, L r) '  (N−1 logd−1(N))s
 (q ≤ 2 ≤ r ≤ p <∞), (q ≤ 1, p = r =∞), (1 < p = r ≤ 2, q ≤ r),
(N−1 logd−1(N))s(logd−1(N))1/r−1/q (1 < p = r ≤ 2, q > r), (N−1 logd−1(N))s(logd−1(N))(1/2−1/q)+ (2 ≤ q, 1 < r < 2 ≤ p <∞),
(b) For 1 < p < r <∞,
λN (MB s p,q, L
r) ' {
(N−1 logd−1(N))s+1/r−1/p (2 ≤ p, 2 ≤ q ≤ r), (N−1 logd−1(N))s+1/r−1/p(logd−1(N))(1/r−1/q)+ (r ≤ 2).
Therefore, the approximation error given in Theorem 1 achieves the optimal linear width ((N−1 logd−1(N))s) for several parameter settings of p, q, s. In particular, when p < r, the bound in Theorem 1 is better than that of Proposition 2. This is because to prove Theorem 1, we used an adaptive recovery technique instead of a linear recovery method. This implies that, by constructing a deep neural network accurately, we achieve the same approximation accuracy as the adaptive one which is better than that of linear approximation.

heading: 4 ESTIMATION ERROR ANALYSIS
text: In this section, we connect the approximation theory to generalization error analysis (estimation error analysis). For the statistical analysis, we assume the following nonparametric regression model:
yi = f o(xi) + ξi (i = 1, . . . , n),
where xi ∼ PX with density 0 ≤ p(x) < R on [0, 1]d, and ξi ∼ N(0, σ2). The data Dn = (xi, yi) n i=1 is independently identically distributed. We want to estimate f
o from the data. Here, we consider a regularized learning procedure:
f̂ = argmin f̄ :f∈Φ(L,W,S,B) n∑ i=1 (yi − f̄(xi))2 (7)
where f̄ is the clipping of f defined by f̄ = min{max{f,−F}, F} for F > 0 which is realized by ReLU units. Since the sparsity level is controlled by S and the parameter is bounded by B, this estimator can be regarded as a regularized estimator. In practice, it is hard to exactly compute f̂ . Thus, we approximately solve the problem by applying sparse regularization such as L1-regularization and optimal parameter search through Bayesian optimization. The generalization error that we present here is an “ideal” bound which is valid if the optimal solution f̂ is computable.

heading: 4.1 ESTIMATION ERROR IN BESOV SPACES
text: In this subsection, we provide the estimation error rate of deep learning to estimate functions in Besov spaces. Theorem 2. Suppose that 0 < p, q ≤ ∞ and s > d(1/p− 1/2)+. If fo ∈ Bsp,q(Ω) ∩ L∞(Ω) and. ‖fo‖Bsp,q ≤ 1 and ‖f
o‖∞ ≤ F for F ≥ 1, then letting (W,L, S,B) be as in Proposition 1 with N n d 2s+d , we obtain
EDn [‖fo − f̂‖2L2(PX)] . n − 2s2s+d log(n)3,
where EDn [·] indicates the expectation w.r.t. the training data Dn.
The proof is given in Appendix D. The condition ‖fo‖∞ ≤ F is required to connect the empirical L2-norm 1n ∑n i=1(f̂(xi) − fo(xi))2 to the population L2-norm ‖f̂ − fo‖2L2(PX). It is known that the convergence rate n− 2s
2s+d is mini-max optimal (Kerkyacharian & Picard, 1992; Donoho et al., 1996; Donoho & Johnstone, 1998; Giné & Nickl, 2015). Thus, it cannot be improved by any estimator. Therefore, deep learning can achieve the minimax optimal rate up to log(n)3-order. The term log(n)3 could be improved to log(n)2 by using the construction of Petersen & Voigtlaender (2017). However, we don’t pursue this direction for simplicity.
Here an important remark is that this minimax optimal rate cannot be achieved by any linear estimator. We call an estimator linear when the estimator depends on (yi)ni=1 linearly (it can be non-linearly dependent on (xi)ni=1). Several classical methods such as the kernel ridge regression, the Nadaraya-Watson estimator and the sieve estimator are included in the class of linear estimators (e.g., kernel ridge regression is given as f̂(x) = kx,X(kXX + λI)−1Y ). The following proposition given by Donoho & Johnstone (1998); Zhang et al. (2002) states that the minimax rate of linear estimators is lower bounded by n−{2s−2(1/p−1/2)+}/{2s+1−2(1/p−1/2)+} for d = 1 which is larger than the minimax rate n− 2s 2s+1 if p < 2.
Proposition 3 (Donoho & Johnstone (1998); Zhang et al. (2002)). Suppose that d = 1 and the input distribution PX is the uniform distribution on [0, 1]. Assume that s > 1/p, 1 ≤ p, q ≤ ∞ or s = p = q = 1. Then,
inf f̂ : linear sup fo∈U(Bsp,q)
EDn [‖fo − f̂‖2L2(PX)] & n − 2s−v2s+1−v
where v = 2/(p∧2)−1 and f̂ runs over all linear estimators, that is, f̂ depends on (yi)ni=1 linearly.
When p < 2, the smoothness of the Besov space is somewhat inhomogeneous, that is, a function in the Besov space contains spiky/jump parts and smooth parts (remember that when s = p = q = 1 for d = 1, the Besov space is included in the set of functions with bounded total variation). Here, the setting p < 2 is the regime where there appears difference between non-adaptive methods and deep learning in terms of approximation accuracy (see Eq. (5)). On the other hand, the linear estimator captures only global properties of the function and cannot capture variability of local shapes of the function. Hence, the linear estimator cannot achieve the minimax optimal rate if the function has spatially inhomogeneous smoothness. However, deep learning possesses adaptivity to the spatial inhomogeneity.
Imaizumi & Fukumizu (2018) has pointed out that such a discrepancy appears when the target function is non-smooth. Interestingly, the parameter setting s > 1/p assumed in Proposition 3 ensures smoothness (see Eq. (2)). This means that non-smoothness is not necessarily required to characterize the superiority of deep learning, but non-convexity of the set of target functions is essentially important. In fact, the gap is coming from the property that the quadratic hull of the model U(Bsp,q) is strictly larger than the original set (Donoho & Johnstone, 1998).

heading: 4.2 ESTIMATION ERROR IN MIXED SMOOTH BESOV SPACES
text: Here, we provide the estimation error rate to estimate functions in mixed smooth Besov spaces. Theorem 3. Suppose that 0 < p, q ≤ ∞ and s > (1/p− 1/2)+. Let u = (1− 1/q)+ for p ≥ 2 and u = (1/2 − 1/q)+ for p < 2. If fo ∈ MBsp,q(Ω) ∩ L∞(Ω) and ‖fo‖MBsp,q ≤ 1 and ‖f o‖∞ ≤ F
for F ≥ 1, then letting (W,L, S,B) be as in Theorem 1, we obtain
EDn [‖fo − f̂‖2L2(PX)] . n − 2s2s+1 log(n) 2(d−1)(u+s) 1+2s log(n)3.
Under the same assumption, if s > u log2(e) is additionally satisfied, we also have
EDn [‖fo − f̂‖2L2(PX)] . n − 2s−2u log2(e) 2s+1+(1−2u) log2(e) log(n)3.
The proof is given in Appendix D. The risk bound (Theorem 3) indicates that the curse of dimensionality can be eased by assuming the mixed smoothness compared with the ordinary Besov space (n− 2s 2s+d ). We show that this is almost minimax optimal in Theorem 4 below. In the first bound, the dimensionality d comes in the exponent of poly log(n) term. If u = 0, then the effect of d can be further eased. Actually, in this situation (u = 0), the second bound can be rewritten as
n − 2s 2s+1+log2(e) log(n)3,
where the effect of the dimensionality d completely disappears from the exponent. This explains partially why deep learning performs well for high dimensional data. Here, we again remark the adaptivity of deep learning. Remind that this rate cannot be achieved by the linear estimator for p < 2 when d = 1 by Proposition 3. Montanelli & Du (2017) has analyzed the mixed smooth Hölder space with s < 2. However, our analysis is applicable to the m-Besov space which is more general than the mixed smooth Hölder space and the covered range of s, p, q is much larger.
Minimax optimal rate for estimating a function in the m-Besov space Here, we show the minimax optimality of the obtained bound as follows. Theorem 4. Assume that 0 < p, q ≤ ∞ and s > (1/p− 1/2)+ and PX is the uniform distribution over [0, 1]d. Regarding d as a constant, the minimax learning rate in the asymptotics of n → ∞ is lower bounded as follows: There exists a constant Ĉ1 such that
inf f̂ sup fo∈U(MBsp,q)
EDn [‖f̂ − fo‖2L2(PX)] ≥ Ĉ1n − 2s2s+1 log(n)
2(d−1)(s+1/2−1/q)+ 2s+1 (8)
where “inf” is taken over all measurable functions of the observations (xi, yi)ni=1 and the expectation is taken for the sample distribution.
The proof is given in Appendix E. Because of this theorem, our bound given in Theorem 3 achieves the minimax optimal rate in the regime of p < 2 and 1/2 − 1/q > 0 up to log(n)3 order. Even outside of this parameter setting, the discrepancy between our upper bound and the minimax lower bound is just a poly-log oder. See also Neumann (2000) for some other related spaces and specific examples such as p = q = 2.

heading: 5 CONCLUSION
text: This paper investigated the learning ability of deep ReLU neural network when the target function is in a Besov space or a mixed smooth Besov space. Based on the analysis for the Besov space, it was shown that deep learning using the ReLU activation can achieve the minimax optimal rate and outperform the linear method when p < 2 which indicates the spatial inhomogeneity of the shape of the target function. The analysis for the mixed smooth Besov space showed that deep learning can adaptively avoid the curse of dimensionality. The bound was derived by sparse grid technique. All analyses in the paper adopted the cardinal B-spline expansion and the adaptive non-linear approximation technique, which allowed us to show the minimax optimal rate. The consequences of the analyses partly support the superiority of deep leaning in terms of adaptivity and ability to avoid curse of dimensionality. From more high level view point, these favorable property is reduced to its high feature extraction ability.
Acknowledgment The author is grateful to Satoshi Hayakawa for showing the correct statement of Proposition 4 by modifying the proof of Lemma 10 of Schmidt-Hieber (2018). TS was partially supported by MEXT Kakenhi (25730013, 25120012, 26280009, 15H05707 and 18H03201), Japan Digital Design and JST-CREST.

heading: A PROOF OF LEMMA 1
text: Proof of Lemma 1. First note that Nm(x) = 1m! ∑m+1 j=0 (−1)j ( m+1 j ) (x − j)m+ (see Eq. (4.28) of Mhaskar & Micchelli (1992) for example). Thus, if we can make an approximation of η(x)m, then by taking a summation of those basis, we obtain an approximate ofNm(x). It is shown by Yarotsky (2017); Schmidt-Hieber (2018) that, for D ∈ N and any > 0, there exists a neural network φmult ∈ Φ(L,W, S,B) with L = dlog2 ( 3D ) + 5edlog2(D)e, W = 6d, S = LW 2 and B = 1
such that
sup x∈[0,1]D ∣∣∣∣∣φmult(x1, . . . , xD)− D∏ i=1 xi ∣∣∣∣∣ ≤ , and φmult(0, . . . , 0) = 0 for y ∈ RD such that ∏D j=1 yj = 0. Moreover, for any M > 0, we can realize the function min{M,max{x, 0}} by a single-layer neural network φ(0,M)(x) := η(x) − η(x−M)(= min{M,max{x, 0}}). Thus, for x ∈ R, it holds that
sup x∈[0,M ] ∣∣φmult(φ(0,1)(x/M), . . . , φ(0,1)(x/M))− (φ(0,1)(x/M))m∣∣ ≤ . Now, Nm(x) = 0 for x 6∈ [0,m+ 1] gives that
Nm(x) = 1
m! m+1∑ j=0 (−1)j ( m+ 1 j ) φ(0,m+1−j)(x− j)m
= 1
m! m+1∑ j=0 (−1)j ( m+ 1 j ) (m+ 1)mφ(0,1−j/(m+1))((x− j)/(m+ 1))m.
Therefore, letting
f(x) = 1
m! m+1∑ j=0 (−1)j(m+1)m ( m+ 1 j ) φmult ( φ(0,1− jm+1 ) ( x− j m+ 1 ) , . . . , φ(0,1− jm+1 ) ( x− j m+ 1 ) ︸ ︷︷ ︸
m-times
) ,
we have that f(x) = 0 for all x ≤ 0 and
sup 0≤x≤m+1
|Nm(x)− f(x)| ≤ 1
m! m+1∑ j=0 ( m+ 1 j ) (m+ 1)m ≤ (m+ 1) m √ 2πmm+1/2e−m 2m+1
≤ e (2e) m
√ m =: ′,
where we used ∑m+1 j=0 ( m+1 j ) = 2m+1 and Stirling’s approximation m! ≥ √ 2πmm+1/2e−m in the second inequality. Hence, we also have that, for all x > m+ 1,
f(x) = 1
m! m+1∑ j=0 (−1)j ( m+ 1 j ) (m+ 1)m
× φmult ( φ(0,1− jm+1 ) ( m+ 1− j m+ 1 ) , . . . , φ(0,1− jm+1 ) ( m+ 1− j m+ 1 )) =: δ′.
It holds that |δ′| ≤ ′. Since it is possible that δ′ 6= 0, we modify f(x) so that f(x) = 0 outside the interval [0,m+ 1]. Because of this and noting 0 ≤ Nm(x) ≤ 1, we see that g(x) := φ(0,1)(f(x)− δ′ m+1φ(0,m+1)(x)) yields
sup x∈R |Nm(x)− g(x)| ≤ 2 ′,
sup x∈R |g(x)| ≤ 1, g(x) = 0 (∀x 6∈ [0,m+ 1]).
Hence, by applying φmult again, we finally obtain that
sup x∈[0,1]d
|Md0,0(x)− φmult(g(x1), . . . , g(xd))|
≤ sup x∈[0,1]d ∣∣∣∣∣∣Md0,0(x)− d∏ j=1 g(xj) ∣∣∣∣∣∣+ supx∈[0,1]d ∣∣∣∣∣∣ d∏ j=1 g(xj)− φmult(g(x1), . . . , g(xd)) ∣∣∣∣∣∣ ≤2d ′ + .
We again applying φ(0,1), we obtain that h(x) = φ(0,1) ◦φmult(g(x1), . . . , g(xd)) satisfies ‖Md0,0− h‖L∞(Rd) ≤ 2d ′+ , h(x) = 0 for all x 6∈ [0,m+1]d, and ‖h‖∞ ≤ 1. Finally, by carefully checking the network construction, it is shown that h ∈ Φ(L,W,S,B) with L = 3 + 2dlog ( 3d∨m ) +
5edlog2(d ∨ m)e, W = 6dm(m + 2) + 2d, S = LW 2 and B = 2(m + 1)m. Hence, resetting ← 2d ′ + = (1 + 2de (2e) m
√ m ) , h becomes the desired M̌ .

heading: B PROOF OF PROPOSITION 1
text: For the order m ∈ N of the cardinal B-spline bases, let J(k) = {−m,−m + 1, . . . , 2k − 1, 2k}d and the quasi-norm of the coefficient (αk,j)k,j for k ∈ Z+ and j ∈ J(k) be
‖(αk,j)k,j‖bsp,q . = ∑ k∈Z+ 2k(s−d/p)( ∑ j∈J(k) |αk,j |p )1/pq 1/q .
Lemma 2. Under one of the conditions (4) in Proposition 1 and the condition 0 < s < min(m,m− 1 + 1/p) where m ∈ N is the order of the cardinal B-spline bases, for any f ∈ Bsp,q(Ω), there exists fN that satisfies
‖f − fN‖Lr(Ω) . N−s/d‖f‖Bsp,q (9) for N 1, and has the following form:
fN (x) = K∑ k=0 ∑ j∈J(k) αk,jM d k,j(x) + K∗∑ k=K+1 nk∑ i=1 αk,jiM d k,ji(x), (10)
where (ji) nk i=1 ⊂ J(k), K = dC1 log(N)/de, K∗ = dlog(λN)ν−1e + K + 1, nk = dλN2−ν(k−K)e (k = K+ 1, . . . ,K∗) for δ = d(1/p− 1/r)+ and ν = (s− δ)/(2δ) where the real number constants C1 > 0 and λ > 0 are chosen to satisfy ∑K k=1(2 k + m)d + ∑K∗ k=K+1 nk ≤ N independently to N . Moreover, we can choose the coefficients (αk,j) to satisfy
‖(αk,j)k,j‖bsp,q . ‖f‖Bsp,q .
Proof of Lemma 2. DeVore & Popov (1988) constructed a linear bounded operator Pk having the following form:
Pk(f)(x) = ∑ j∈J(k) ak,jM d k,j(x) (11)
where αk,j is constructed in a certain way, where for every f ∈ Lp([0, 1]d) with 0 < p ≤ ∞, it holds
‖f − Pk(f)‖Lp ≤ Cwr,p(f, 2−k). (12) Let pk(f) := Pk(f)− Pk−1(f), P−1(f) = 0. Then, it is shown that for 0 < p, q ≤ ∞ and 0 < s < min(m,m − 1 + 1/p), f belongs to Bsp,q if and only if f can be decomposed into
f = ∞∑ k=0 pk(f),
with the convergence condition ‖(pk(f))∞k=0‖bsp(Lp) := ( ∑ k∈Z+(2
sk‖pk‖Lp)q)1/q < ∞; in particular, ‖f‖Bsp,q ' ‖(pk(f)) ∞ k=0‖bsp(Lp). Here, each pk can be expressed as pk(x) =∑
j∈J(k) αk,jM d k,j(x) for a coefficient (αk,j)k,j which could be different from (ak,j)k,j appear-
ing in Eq. (11). Hence, f ∈ Bsp,q can be decomposed into
f = ∞∑ k=0 ∑ j∈J(k) αk,jM d k,j(x) (13)
with convergence in the sence of Lp. Moreover, it is shown that ‖pk‖Lp ' (2−kd ∑ j∈J(k) |αk,j |p)1/p and thus
‖f‖Bsp,q ' ‖(αk,j)k,j‖bsp,q . (14)
Based on this decomposition, Dũng (2011b) proposed an optimal adaptive recovery method such that the approximator has the form (10) under the conditions for K,K∗, nk given in the statement and satisfies the approximation accuracy (9). This can be proven by applying the proof of Theorem 3.1 in Dũng (2011b) to the decomposition (13) instead of Eq. (3.8) of that paper. See also Theorem 5.4 of Dũng (2011b). Moreover, the equivalence (14) gives the norm bound of the coefficient (αk,j).
Proof of Proposition 1. Basically, we combine Lemma 1 and Lemma 2. We substitute the approximated cardinal B-spline basis M̌ into the decomposition of fN (10). Let the set of indexes (k, j) ∈ Z×Z that consists fN given in Eq. (10) be EN ; i.e., fN = ∑ (k,j)∈EN αk,jM d k,j . Accord-
ingly, we set f̌ := ∑
(k,j)∈EN αk,jM̌ d k,j . For each x ∈ Rd, it holds that |fN (x)− f̌(x)| ≤ ∑
(k,j)∈EN
|αk,j ||Mdk,j(x)− M̌dk,j(x)|
≤ ∑
(k,j)∈EN
|αk,j |1{Mdk,j(x) 6= 0}
≤ (m+ 1)d(1 +K∗)2K ∗(d/p−s)+‖f‖Bsp,q . log(N)N (ν −1+d−1)(d/p−s)+ ‖f‖Bsp,q ,
where we used the definition of K∗ in the last inequality. Therefore, for each f ∈ U(Bsp,q([0, 1]d)), it holds that
‖f − f̌‖Lr . ‖f − fN‖Lr + ‖fN − f̌‖Lr . log(N)N (ν −1+d−1)(d/p−s)+‖f‖Bsp,q +N −s/d.
By taking to satisfy log(N)N (ν −1+d−1)(d/p−s)+ ≤ N−s/d (i.e., ≤ N−s/d−(ν −1+d−1)(d/p−s)+ log(N)−1), then we obtain the approximation error bound.
Next, we bound the magnitude of the coefficients. Each coefficient αj,k satisfies |αj,k| . 2k(d/p−s)+‖f‖Bsp,q ≤ 2 k(d/p−s)+ . N (ν −1+d−1)(d/p−s)+ for k ≤ K∗. Finally, the magnitudes of the coefficients hidden in M̌dk,j are evaluated. Remembering that M̌ m k,j(x) = M̌(2
kx1 − j1, . . . , 2
kxd − jd), we see that we just need to bound the quantity 2k (k ≤ K∗). However, this is bounded by 2k ≤ Nν−1+d−1 for k ≤ K∗. Hence, we obtain the assertion.

heading: C PROOF OF THEOREM 1
text: C.1 PREPARATION: SPARSE GRID
Here, we give technical details behind the approximation bound. The analysis utilizes the so called sparse grid technique Smolyak (1963) which has been developed in the function approximation theory field.
As we have seen in the above, in a typical B-spline approximation scheme, we put the basis functions Mmk,j(x) on a “regular grid” for k = 1, . . . ,K and (j1, . . . , jd) ∈ J(k), and take its superposition as f(x) ≈ ∑ k=1,...,K ∑ j∈J(k) αk,jM m k,j(x), which consists of O(2
Kd) terms (see Eq. (10)). Hence, the number of parameters O(2Kd) is affected by the dimensionality d in an exponential order. However, to approximate functions with mixed smoothness, we do not need to put the basis on the whole range of the regular grid. Instead, we just need to put them on a sparse grid which is a subset of the regular grid and has much smaller cardinality than the whole set. The approximation algorithm utilizing sparse grid is based on Smolyak’s construction (Smolyak, 1963) and its applications to mixed smooth spaces (Dũng, 1990; 1991; 1992; Temlyakov, 1982; 1993a;b). Dũng (2011a) studied an optimal non-adaptive linear sampling recovery method for the mixed smooth Besov space based on the cardinal B-spline bases. We adopt this method, and combining this with the adaptive technique developed in Dũng (2011b), we give the following approximation bound using a non-linear adaptive method to obtain better convergence for the setting p < r.
Before we state the theorem, we define an quasi-norm of a set of coefficients αk,j ∈ R for k ∈ Zd+ and j ∈ Jdm(k) := {−m,−m+ 1, . . . , 2k1 − 1, 2k1} × · · · × {−m,−m+ 1, . . . , 2kd − 1, 2kd} as
‖(αk,j)k,j‖mbsp,q := ∑ k∈Zd+ 2(s−1/p)‖k‖1( ∑ j∈Jdm(k) |αk,j |p )1/pq1/q .
Theorem 5. Suppose that 0 < p, q, r ≤ ∞ and s > (1/p − 1/r)+. Assume that the order m ∈ N of the cardinal B-spline satisfies 0 < s < min(m,m− 1 + 1/p). Let δ = (1/p− 1/r)+. Then, for any f ∈MBsp,q(Ω) and K > 0, there exists RK(f) such that RK(f) can be represented as
RK(f)(x) = ∑ k∈Zd+: ‖k‖1≤K ∑ j∈Jdm(k) αk,jM d k,j(x) + ∑ k∈Zd+: K<‖k‖1≤K∗ nk∑ i=1 α k,j (k) i Md k,j (k) i (x),
where K∗ = dK(1 + 2δs−δ )e, (j (k) i ) nk i=1 ⊂ Jdm(k), and nk = d2K− s−δ 2δ (‖k‖1−K)e, and satisfies the following properties:
(i) For p ≥ r,
‖f −RK(f)‖r . 2−KsD(1/min(r,1)−1/q)+K,d ‖f‖MBsp,q .
(ii) For p < r,
‖f −RK(f)‖r .
{ 2−KsD
(1/r−1/q)+ K,d ‖f‖MBsp,q (r <∞),
2−KsD (1−1/q)+ K,d ‖f‖MBsp,q (r =∞).
Moreover, the coefficients (αk,j)k,j can be taken to hold ‖(αk,j)k,j‖mbsp,q . ‖f‖MBsp,q .
The proof is given in Appendix C.2. The total number of cardinal B-spline bases consisting of RK(f) can be evaluated as
2K+1 ( K + d− 1 d− 1 ) + ∑ k:K<‖k‖1≤K∗ nk
.2KDK,d + 2 KDK∗,d . 2 KDK,d (∵ Eq. (17)).
Here, DK,d can be evaluated as
DK,d . K d−1 or DK,d . dK .
Therefore, the total number of bases can be evaluated as 2K min{Kd−1, dK}
which is much smaller than 2Kd which is required to approximate functions in the ordinal Besov space (see Lemma 2). In this proposition, K controls the resolution and as M goes to infinity, the approximation error goes to 0 exponentially fast. A remarkable point in the proposition is in the construction of RK(f) in which the superposition is taken over ‖k‖1 ≤ M instead of ‖k‖∞ ≤ K∗ = O(K). Hence, the number of terms appearing in the summation is at most O(2KKd−1) while the full grid takes O(2Kd) terms. This represents how the mixed smoothness is important to ease the curse of dimensionality.
Several aspects of the m-Besov space such as the optimal N -term approximation error and Kolmogorov widths have been extensively studied in the literature (see a comprehensive survey (Dũng et al., 2016)). An analogous result is already given by Dũng (2011a) in which s > 1/p is assumed and a linear interpolation method is investigated. However, our result only requires s > (1/p − 1/q)+. This difference comes from a point that our analysis allows nonlinear adaptive interpolation instead of (linear) non-adaptive sampling considered in Dũng (2011a). Because of this, our bound is better than the optimal rate of linear methods (Galeev, 1996; Romanyuk, 2001) and non-adaptive methods (Dũng, 1990; 1991; 1992; Temlyakov, 1982; 1993a;b) especially in the regime of p < r (Dũng (1992) also deals with adaptive methods but does not cover p < r for adaptive method). See Proposition 2 for comparison.
C.2 PROOF OF THEOREM 5
Proof of Theorem 5. For k = (k1, . . . , kd) ∈ Zd+, let P (i) ki f(x) be the function operating Pki defined in (11) to f as a function of xi with other components xj (j 6= i) fixed, and let
pk := d∏ i=1 (P (i) ki − P (i)ki−1)f. (15)
Then, pk can be expressed as pk(x) = ∑ j∈Jdm(k) αk,jM d k,j(x).
Let T (i)ki = I − P (i) ki and ‖f‖p,i be the Lp-norm of f as a function of xi with other components xj (j 6= i) fixed (i.e., if p <∞, ‖f‖pp,i = ∫ |f(x)|pdxi), then Eq. (12) gives
‖T (i)ki f‖p,i . sup |hi|≤2−ki ‖∆r,ihi (f)‖p,i.
Thus, by applying the same argument again, it also holds
‖‖T (i)ki T (j) kj f‖p,i‖p,j . ‖ sup |hi|≤2−ki ‖∆r,ihi (Tkjf)‖p,i‖p,j
= sup |hi|≤2−ki
‖‖Tkj∆ r,i hi (f)‖p,j‖p,i (∵ the definition of ∆r,ihi and Fubini’s theorem)
. sup |hi|≤2−ki sup |hj |≤2−kj
‖‖∆r,ihi (∆ r,j hj (f))‖p,j‖p,i,
for i 6= j. Thus, applying the same argument recursively, for u ⊂ [d], it holds that∥∥∥∥∥∏ i∈u T (i) ki f ∥∥∥∥∥ p . wur,p(f, 2 −k),
for k ∈ Zd+(u). Therefore, since pk = ∏d i=1(T (i) ki−1 − T (i) ki
)f =∑ u⊂[d](−1)|u| (∏ i∈u T (i) ki ∏ i 6∈u T (i) ki−1 ) f , by letting e = {i | ki > 0}, we have that
‖pk‖p . ∑ u⊂[d] ∥∥∥∥∥∥ ∏ i∈u T (i) ki ∏ i 6∈u T (i) ki−1  f ∥∥∥∥∥∥ p . ∑ u⊂[d] wûr,p(f, 2 −(ku)û) . ∑ e⊂u wur,p(f, 2 −ku)
where kui := ki (i ∈ u) and kui := ki− 1 (i 6∈ u), û = {i | kui ≥ 0}, and (ku)û is a vector such that (ku)û,i = k u i for i ∈ û and (ku)û,i = 0 for i 6∈ û. Now let
‖(pk)k‖bαq (Lp) = ( ∑ k∈Zd+ ( 2α‖k‖1‖pk‖Lp )q)1/q
for pk ∈ Lp(Ω) (k ∈ Zd+). Hence, if we set ak = ∑ e⊂u w u r,p(f, 2
−ku) for k ∈ Zd+ and e = {i | ki > 0}, we have that
‖(pk)k‖bαq (Lp) . ‖(ak)‖bαq (Lp) ' ‖f‖MBsp,q .
On the other hand, following the same line of Theorem 2.1 (ii) of Dũng (2011a), we also obtain the opposite inequality ‖f‖MBsp,q ' ‖(ak)‖bαq (Lp) . ‖(pk)k‖bαq (Lp) (note that the analogous inequality to Lemma 2.3 of Dũng (2011a) also holds in our setting by replacing qs with ps and ωer(f, 2
−k)p by wer,p).
Therefore, f ∈ MBαp,q if and only if (pk)k∈Zd+ given by Eq. (15) satisfies ‖(pk)k‖bαq (Lp) < ∞ and f can be decomposed into f = ∑ k∈Zd+
pk where convergence is in MBαp,q . Moreover, it holds that ‖f‖MBαp,q ' ‖(pk)k‖bαq (Lp). This can be shown by Theorem 2.1 of Dũng (2011a). Moreover, by the quasi-norm equivalence ‖pk‖p ' 2−‖k‖1/p( ∑ j∈Jdm(k)
|αk,j |p)1/p, we also have ‖(αk,j)k,j‖mbαp,q ' ‖f‖MBαp,q .
If p ≥ r, the assertion can be shown in the same manner as Theorem 3.1 of Dũng (2011a). For the setting of p < r, we need to use an adaptive approximation method. In the following, we assume p < r. For a given K, by choosing K∗ appropriately later, we set
RK(f)(x) = ∑
k∈Zd+:‖k‖1≤K
pk + ∑
k∈Zd+:K<‖k‖1≤K∗ Gk(pk),
where Gk(pk) is given as Gk(pk) = ∑ 1≤i≤nk αk,jiM d k,ji(x)
where (αk,ji) |Jdm(k)| i=1 is the sorted coefficients in decreasing order of their absolute value: |αk,j1 | ≥ |αk,j2 | ≥ · · · ≥ |αk,j|Jdm(k)| |. Then, it holds that
‖pk −Gk(pk)‖r ≤ ‖pk‖p2δ‖k‖1n−δk ,
where δ := (1/p − 1/r) (see the proof of Theorem 3.1 of Dũng (2011b) and Lemma 5.3 of Dũng (2011a)). Moreover, we also have
‖pk‖r ≤ ‖pk‖p2δ‖k‖1
for k ∈ Zd+ with ‖k‖1 > K∗. Here, we define N as N = dlog2(K)e. Let = (s− δ)/(2δ), and
K∗ = dK(1 + 1/ )e,
and nk = d2K− (‖k‖1−K)e for k ∈ Zd+ with K + 1 ≤ ‖k‖1 ≤ K∗. Then, by Lemma 5.3 of Dũng (2011a), we have that
‖f −RK(f)‖rLr . ∑
K<‖k‖1≤K∗ ‖pk −Gk(pk)‖rLr + ∑ K∗<‖k‖1 ‖pk‖rLr
. ∑
K<‖k‖1≤K∗ [‖pk‖p2δ‖k‖1n−δk ]
r + ∑
K∗<‖k‖1
[2δ‖k‖1‖pk‖Lp ]r. (16)
In the following, we require an upper bound of ( k+d−1 d−1 ) . Hence, we evaluate this quantity beforehand. This can be upper bounded by the Stering’s formula as( k + d− 1 d− 1 ) ≤ √ 2e 2π ( 1 + d− 1 k )k ( 1 + k d− 1 )d−1 ︸ ︷︷ ︸
=Dk,d
≤ Dk,d.
Let ξ > 0 be a positive real number satisfying 1 + ξ ≥ K∗/K. We can see that ξ can be chosen as ξ = 1/ + o(1). Then, we have that
DK∗,d = DK,d (1 + d−1K∗ )
K∗
(1 + d−1K ) K
(1 + K ∗ d−1 ) d−1 (1 + Kd−1 ) d−1 ≤ DK,d (1 + d−1K ) K∗ (1 + d−1K ) K
( 1
1 + Kd−1 +
K∗
(d− 1)(1 + Kd−1 )
)d−1
≤ DK,d (
1 + d− 1 K
)K∗−K ( d− 1 +K∗
d− 1 +K
)d−1 = DK,d ( 1 +
d− 1 K
)ξK (1 + ξ) d−1
≤ DK,de(d−1)ξ(1 + ξ)d−1 ' DK,d. (17)
(a) Suppose that q ≤ r and r <∞. Then
‖f −RK(f)‖qLr = ‖f −RK(f)‖ r qr Lr
.  ∑ K<‖k‖1≤K∗ [2δ‖k‖1n−δk ‖pk‖Lp ] r + ∑ K∗<‖k‖1 [2δ‖k‖1‖pk‖Lp ]r  q r
(∵ Eq. (16))
. ∑
K<‖k‖1≤K∗ [2δ‖k‖1n−δk ‖pk‖Lp ]
q + ∑
K∗<‖k‖1
[2δ‖k‖1‖pk‖Lp ]q
≤ N−δq2−(s−δ)Kq ∑
K<‖k‖1≤K∗ [2−(s−δ−δ )(‖k‖1−K)︸ ︷︷ ︸ ≤1 2s‖k‖1‖pk‖Lp ]q + 2−q(s−δ)K ∗ ∑ K∗<‖k‖1 [2s‖k‖1‖pk‖Lp ]q
. (N−δ2−(s−δ)K + 2−(s−δ)K ∗ )q‖f‖qMBsp,q ≤ (N−s)q‖f‖qMBαp,q .
(b) Suppose that q > r and r <∞. Then, letting ν = q/r (> 1) and ν′ = 1/(1− 1/ν) = q/(q− r) (note that 1ν + 1 ν′ = 1), we have
‖f −RK(f)‖rLr . ∑
K<‖k‖1≤K∗ [2δ‖k‖1n−δk ‖pk‖Lp ]
r + ∑
K∗<‖k‖1
[2δ‖k‖1‖pk‖Lp ]r (∵ Eq. (16))
≤ N−δr2−(s−δ)Kr ∑
K<‖k‖1≤K∗ [2−(s−δ−δ )(‖k‖1−K)2s‖k‖1‖pk‖Lp ]r + ∑ K∗<‖k‖1 [2s‖k‖1‖pk‖Lp ]r(2−(s−δ)‖k‖1)r
≤ (N−δ2−(s−δ)K + 2−(s−δ)K ∗ )r { ∑ K<‖k‖1≤K∗ [2−(s−δ−δ )(‖k‖1−K)2s‖k‖1‖pk‖Lp ]r
+ ∑
K∗<‖k‖1
[2s‖k‖1‖pk‖Lp ]r2−(s−δ)(‖k‖1−K ∗)r }
≤ (N−δ2−(s−δ)K + 2−(s−δ)K ∗ )r  ∑ K<‖k‖1≤K∗ [2s‖k‖1‖pk‖Lp ]rν + ∑ K∗<‖k‖1 [2s‖k‖1‖pk‖Lp ]rν  1/ν
×  ∑ K<‖k‖1≤K∗ [2−(s−δ−δ )(‖k‖1−K)]rν ′ + ∑ K∗<‖k‖1 [2−(s−δ)(‖k‖1−K ∗)]rν ′  1/ν′
. (N−δ2−(s−δ)K + 2−(s−δ)K ∗ )r‖f‖rMBsp,qD r(1/r−1/q) K,d (∵ Eq. (17)) . (N−sD1/r−1/qK,d ) r‖f‖rMBsp,q .
(c) Suppose that r =∞. Then, similarly to the analysis in (b), we can evaluate ‖f −RK(f)‖Lr
. N−δ2−(s−δ)K ∑
K<‖k‖1≤K∗ [2−(s−δ−δ )(‖k‖1−K)2s‖k‖1‖pk‖Lp ] + ∑ K∗<‖k‖1 [2s‖k‖1‖pk‖Lp ](2−(s−δ)‖k‖1)
. (N−δ2−(s−δ)K + 2−(s−δ)K ∗ )D (1−1/q)+ K,d ‖f‖MBsp,q . N−sD(1−1/q)+K,d ‖f‖MBsp,q .
C.3 PROOF OF THEOREM 1
Let Zd+(e) := {k ∈ Zd+ | ki = 0, i 6∈ e} and for k ∈ Zd+(e), we define 2−k := (2−ki1 , . . . , 2
−ki|e| ) ∈ R|e|+ where {i1, . . . , i|e|} = e. By defining ‖(gk)k‖bs,eq :=(∑ k∈Zd+(e) (2s‖k‖1 |gk|)q )1/q for a sequence (gk)k∈Zd+(e), then it holds that
|f |MBs,ep,q ' ∑
e⊂{1,...,d}
‖(wer,p(f, 2−k))k‖bs,eq .
Then, we can prove Theorem 1 based on Theorem 5 as follows.
Proof of Theorem 1. The result is immediately follows from Theorem 5. Let the set of indexes of (k, j) consisting of RK be EK : RK(f) = ∑ (k,j)∈EK αk,jM d k,j(x). As in the proof of Proposition 1, we approximate RK(f) by a neural network given as
f̌(x) = ∑
(k,j)∈EK
αk,jM̌ d k,j(x).
Each coefficient αj,k satisfies |αj,k| . 2‖k‖1(1/p−s)+‖f‖MBsp,q . 2 K∗(1/p−s)+ . The difference between RK(f) and f̌ can be evaluated as
|RK(f)− f̌(x)| ≤ ∑
(k,j)∈EK
|αk,j ||Mdk,j(x)− M̌dk,j(x)|
≤ ∑
(k,j)∈EK
|αk,j |1{Mdk,j(x) 6= 0}
. (m+ 1)d(1 +K∗)DK∗,d2 K∗(1/p−s)+‖f‖MBsp,q .
Therefore, by taking so that (m + 1)d(1 + K∗)DK∗,d2K ∗(1/p−s)+ ≤ 2−Ks is satisfied, it holds that |RK(f)− f̌(x)| . 2−Ks.
By the inequality DK∗,d ≤ eK ∗+d−1, it suffices to let ≤ e
−K∗(s+(1/p−s)++1) [e(m+1)]d(1+K∗) . The cardinality of
E(K) is bounded as∑ κ=0,...,K 2κ ( κ+ d− 1 d− 1 ) + ∑ k:K<‖k‖1≤K∗ nk
≤2K+1 ( K + d− 1 d− 1 ) + ∑ K<κ≤K∗ 2K− s−δ 2δ (κ−K) ( κ+ d− 1 d− 1 ) ≤2K+1DK,d + 2K(1− 2− s−δ 2δ )−1DK∗,d ≤ (2 + (1− 2− s−δ 2δ )−1)2KDK∗,d = N.
Since each unit M̌dk,j requires width W0, the whole width becomes W = NW0. The number of nonzero parameters to construct M̌dk,j is bounded by S = (L− 1)W 20N +N .
Finally, the magnitudes of the coefficients hidden in M̌dk,j are evaluated. Remembering that M̌dk,j(x) = M̌(2 k1x1− j1, . . . , 2kdxd− jd), here maximum of 2kj is bounded by 2K ∗ . N (1+1/ν). Hence, we obtain the assertion. Similarly, it holds that |αj,k| . N (1+1/ν){1∨(1/p−s)+}.

heading: D PROOFS OF THEOREMS 2 AND 3
text: Proof of Theorem 2. We use Proposition 4. We just need to evaluate the covering number of F̂ = {f̄ | f ∈ Ψ(L,W, S,B)} for (L,W, S,B) given in Theorem 1 where f̄ is the clipped function for a given f . Note that the covering number of F̂ is not larger than that of Ψ(L,W, S,B). Hence, it is sufficient to evaluate that of Ψ(L,W, S,B). From Lemma 3, the covering number is obtained as
logN(δ, F̂ , ‖ · ‖∞) . N log(N)[log(N)2 + log(δ−1)].
From Proposition 1, it holds that
‖fo −RK(fo)‖2 . N−s/d.
Note that ‖f − fo‖2L2(PX) . ‖f − f o‖22.
for any f : [0, 1]d → R because p(x) ≤ R. Therefore, by applying Proposition 4 with δ = 1/n, we have that
EDn [‖f̂ − fo‖2L2(PX)] . N −2s/d +
N log(N)(log(N)2 + log(n))
n +
1 n . (18)
Here, the right hand side is minimized by setting N n d
2s+d up to log(n)3-order, and then have an upper bound of the RHS as
n− 2s 2s+d log(n)3.
This gives the assertion.
Proof of Theorem 3. The proof follows the almost same line as the proof of Theorem 2. By noting S = O(2KDK,d), L = O(K) and W = O(2KDK,d), Lemma 3 gives an upper bound of the covering number as
logN(δ, F̂ , ‖ · ‖∞) . 2KDK,d[K2 log(2KDK,d) + log(δ−1)] . 2KDK,d(K3 + log(1/δ)).
Letting r = 2, we have that ‖fo −RK(fo)‖2 . 2−sKDuK,d
where u = (1− 1/q)+ for p ≥ 2 and u = (1/2− 1/q)+ for p < 2. Then, by noting that
‖f − fo‖2L2(PX) . ‖f − f o‖22,
for any f : [0, 1]d → R, and by applying Proposition 4 with δ = 1/n, we have that
EDn [‖f̂ − fo‖2L2(PX)] . 2 −2sKD2uK,d +
2KDK,d(K 3 + log(δ−1))
n +
1 n . (19)
Here, we use the following evaluations for DK,d: (a) DK,d . Kd−1, and (b) DK,d . [e(1 + dK )] K .
(a) For the evaluation, DK,d . Kd−1, we have an upper bound of the right hand side of Eq. (19) as
2−2sKK2u(d−1) + 2KKd−1(K3 + log(n))
n ,
which is minimized by setting K = d 11+2s log2(n) + (2u−1)(d−1)
1+2s log2 log(n)e up to log log(n)order. In this situation, we have the generalization error bound as
n− 2s 2s+1 log(n) 2(d−1)(u+s) 1+2s log(n)3.
(b) For the evaluation, DK,d . [e(1 + dK )] K ≤ eKed, Eq. (19) gives an upper bound of
2−2sKe2uK + 2KeK(K2 + log(n))
n .
Then, the right hand side is minimized by K = d 11+2s+(1−2u) log2(e) log2(n)e. Then, we have that
n − 2s−2u log2(e) 1+2s+(1−2u) log2(e) log(n)2.
This gives the assertion.

heading: E MINIMAX OPTIMALITY
text: Proof of Theorem 4. First note that since PX is the uniform distribution, it holds that ‖ · ‖L2(PX) = ‖ · ‖L2([0,1]d). The -covering numberN ( ,G, L2(PX)) with respect to L2(PX) for a function class G is the minimal number of balls with radius measured by L2(PX)-norm needed to cover the set G (van der Vaart & Wellner, 1996). The δ-packing number M(δ,G, L2(PX)) of a function class G with respect to L2(PX) norm is the largest number of functions {f1, . . . , fM} ⊆ G such that ‖fi − fj‖L2(PX) ≥ δ for all i 6= j. It is easily checked that
N (δ/2,G, L2(PX)) ≤M(δ,G, L2(PX)) ≤ N (δ,G, L2(PX)). (20)
For a given δn > 0 and εn > 0, let Q be the δn packing number M(δn, U(MBsp,q), L2(PX)) of U(MBsp,q) and N be the εn covering number of that. Raskutti et al. (2012) utilized the techniques developed by Yang & Barron (1999) to show the following inequality in their proof of Theorem 2(b):
inf f̂ sup f∗∈U(MBsp,q) EDn [‖f̂ − f∗‖2L2(PX)] ≥ inf f̂ sup f∗∈U(MBsp,q)
δ2n 2 P [‖f̂ − f∗‖2L2(PX) ≥ δ 2 n/2]
≥ δ 2 n
2
( 1− log(N) + n2σ2 ε 2 n + log(2)
log(Q)
) .
Thus by taking δn and εn to satisfy n
2σ2 ε2n ≤ log(N), (21) 8 log(N) ≤ log(Q), (22) 4 log(2) ≤ log(Q), (23)
the minimax rate is lower bounded by δ 2 n
4 . This can be achieved by properly setting εn ' δn. Now, for given N with respect to δn > 0, M = log(N) satisfies
δn &M −s log(M)(d−1)(s+1/2−1/q)+
(Theorem 6.24 of Dũng et al. (2016)). Hence, it suffices to take
M ' n 1 2s+1 log(n) 2(d−1)(s+1/2−1/q)+ 2s+1 , (24)
εn ' δn ' n− 2s 2s+1 log(n) 2(d−1)(s+1/2−1/q)+ 2s+1 , (25)
which gives the assertion.

heading: F AUXILIARY LEMMAS
text: Let the -covering number with respect to L2(PX) for a function class G be N ( ,G, L2(PX)) as defined in the proof of Theorem 4.
Proposition 4 (Schmidt-Hieber (2018)). Let F be a set of functions. Let f̂ be the least squares estimator in F:
f̂ = argmin f∈F n∑ i=1 (yi − f(xi))2.
Assume that ‖fo‖∞ ≤ F and all f ∈ F satisfies ‖f‖∞ ≤ F for some F ≥ 1. If δ > 0 satisfies

heading: N (δ,F , ‖ · ‖∞) ≥ 3, then it holds that
text: EDn [‖f̂−fo‖2L2(PX)] ≤ C [ inf f∈F ‖f − fo‖2L2(PX) + (F 2 + σ2) logN (δ,F , ‖ · ‖∞) n + δ(F + σ) ] ,
where C is a universal constant.
Proof of Proposition 4. This is almost direct consequence of Lemma 10 of Schmidt-Hieber (2018)5. The only difference is the assumption of ‖f‖∞ ≤ F for f ∈ F and f = fo while Lemma 10 of Schmidt-Hieber (2018) assumed 0 ≤ f(x) ≤ F ′ for F ′ > 1. However, this can be easily fixed by shifting the function value by +F then the range of f is modified to [0, 2F ]. Then, our situation is reduced to that of Lemma 10 of Schmidt-Hieber (2018) by substituting F ′ ← 2F .
Lemma 3 (Covering number evaluation). The covering number of Φ(L,W, S,B) can be bounded by
logN (δ,Φ(L,W, S,B), ‖ · ‖∞) ≤ S log(δ−1L(B ∨ 1)L−1(W + 1)2L) ≤ 2SL log((B ∨ 1)(W + 1)) + S log(δ−1L).
Proof of Lemma 3. Given a network f ∈ Φ(L,W, S,B) expressed as
f(x) = (W(L)η(·) + b(L)) ◦ · · · ◦ (W(1)x+ b(1)),
let Ak(f)(x) = η ◦ (W(k−1)η(·) + b(k−1)) ◦ · · · ◦ (W(1)x+ b(1)),
and Bk(f)(x) = (W(L)η(·) + b(L)) ◦ · · · ◦ (W(k)η(x) + b(k)),
for k = 2, . . . , L. Corresponding to the last and first layer, we define BL+1(f)(x) = x and A1(f)(x) = x. Then, it is easy to see that f(x) = Bk+1(f) ◦ (W(k) · +b(k)) ◦ Ak(f)(x). Now, suppose that a pair of different two networks f, g ∈ Φ(L,W, S,B) given by
f(x) = (W(L)η(·)+b(L))◦· · ·◦(W(1)x+b(1)), g(x) = (W(L) ′ η(·)+b(L) ′ )◦· · ·◦(W(1) ′ x+b(1) ′ ),
has a parameters with distance δ: ‖W(`) −W(`)′‖∞ ≤ δ and ‖b(`) − b(`) ′‖∞ ≤ δ. Now, not that ‖Ak(f)‖∞ ≤ maxj ‖W(k−1)j,: ‖1‖Ak−1(f)‖∞ + ‖b(k−1)‖∞ ≤ WB‖Ak−1(f)‖∞ + B ≤ (B ∨ 1)(W + 1)‖Ak−1(f)‖∞ ≤ (B ∨ 1)k−1(W + 1)k−1, and similarly the Lipshitz continuity of Bk(f) with respect to ‖ · ‖∞-norm is bounded as (BW )L−k+1. Then, it holds that
|f(x)− g(x)|
= ∣∣∣∣∣ L∑ k=1 Bk+1(g) ◦ (W(k) ·+b(k)) ◦ Ak(f)(x)− Bk+1(g) ◦ (W(k) ′ ·+b(k) ′ ) ◦ Ak(f)(x) ∣∣∣∣∣ ≤
L∑ k=1 (BW )L−k‖(W(k) ·+b(k)) ◦ Ak(f)(x)− (W(k) ′ ·+b(k) ′ ) ◦ Ak(f)(x)‖∞
≤ L∑ k=1 (BW )L−kδ[W (B ∨ 1)k−1(W + 1)k−1 + 1]
5We noticed that there exit some technical flaws in the proof of the lemma, e.g., an incorrect application of the uniform bound to derive the risk of an estimator. However, these flaws can be fixed and the statement itself holds with a slight modification.
≤ L∑ k=1 (BW )L−kδ(B ∨ 1)k−1(W + 1)k ≤ δL(B ∨ 1)L−1(W + 1)L.
Thus, for a fixed sparsity pattern (the locations of non-zero parameters), the covering number is bounded by ( δ/[L(B ∨ 1)L−1(W + 1)L] )−S . There are the number of configurations of the spar-
sity pattern is bounded by ( (W+1)L
S
) ≤ (W + 1)LS . Thus, the covering number of the whole space
Φ is bounded as (W + 1)LS { δ/[L(B ∨ 1)L−1(W + 1)L] }−S = [δ−1L(B ∨ 1)L−1(W + 1)2L]S ,
which gives the assertion.

Review summaries:
This paper describes approximation and estimation error bounds for functions in Besov spaces using estimators corresponding to deep ReLU networks. The general idea of connecting network parameters such as depth, width, and sparsity to classical function spaces is interesting and could lead to novel insights into how and why these networks work and under what settings. The authors carefully define Besov spaces and related literature, and overall the paper is clearly written.