Title: ADAPTIVITY OF DEEP RELU NETWORK FOR LEARN-ING IN BESOV AND MIXED SMOOTH BESOV SPACES: OPTIMAL RATE AND CURSE OF DIMENSIONALITY

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. Function class Hölder Barron class m-Sobolev (0 < β ≤ 2) m-Besov (0 < β) Approximation

1. INTRODUCTION:

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 L 1 (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 Hö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 Table 1 : Comparison between the performances of deep learning and linear methods. N is the number of parameters to approximate a function in a Besov space (B s p,q ([0, 1] d )), and n is the sample size. The approximation error is measured by L r -norm. The Õ symbol hides the poly-log order.

Model:

Deep learning Linear method Approximation error rate Õ(N

-s d ) Õ N -s d +( 1 p -1 r )+

Estimation error rate Õ(n -2s 2s+d ) Ω n -2s-(2/(p∧1)-1)

2s+1-(2/(p∧1)-1)

functions in the Hö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 Hö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 Hö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 (B s p,q ) for s > 0 and 0 < p, q ≤ ∞ with s > d(1/p -1/r) + where L r -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.

2. SET UP OF FUNCTION SPACES:

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 Ω ⊂ R d be a domain of the functions. Throughout this paper, we employ -

Ω = [0, 1] d . For a function f : Ω → R, let f p := f L p (Ω) := ( Ω |f | p dx) 1/p for 0 < p < ∞. For p =

2β 2β+1 × log(n) 2(d-1)(u+β) 1+2β ) ∞, we define f ∞ := f L ∞ (Ω) := sup x∈Ω |f (x)|. For α ∈ R d , let |α| = d j=1 |α j |. Let C 0 (Ω) be the set of continuous functions equipped with L ∞ -norm: C 0 (Ω) := {f : Ω → R | f is continuous and f ∞ < ∞} 2 . For α ∈ Z d + , we denote by D α f (x) = ∂ |α| f ∂ α 1 x1...∂ α d x d (x) 3 . Definition 1 (Hölder space (C β (Ω))). Let β > 0 with β ∈ N be the smoothness parameter. For an m times differentiable function f : R d → R, let the norm of the Hölder space C β (Ω) be f C β := max |α|≤m D α f ∞ + max |α|=m sup x,y∈Ω |D α f (x)-D α f (y)| |x-y| β-m

, where m = β (the largest integer less than β). Then, (β-)Hölder space C β (Ω) is defined as

C β (Ω) = {f | f C β < ∞}.

The parameter β > 0 controls the "smoothness" of the function. Along with the Hölder space, the Sobolev space is also important. Definition 2 (Sobolev space (W k p (Ω))). Sobolev space (W k p (Ω)) with a regularity parameter k ∈ N and a parameter 1 ≤ p ≤ ∞ is a set of functions such that the Sobolev norm

f W k p := ( |α|≤k D α f p p )

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 ∈ L p (Ω) for some p ∈ (0, ∞], the r-th modulus of smoothness of f is defined by w r,p (f, t) = sup

h∈R d : h 2 ≤t ∆ r h (f ) p , where ∆ r h (f )(x) = r j=0 r j (-1) r-j f (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 :

= α + 1, let the seminorm | • | B α p,q be |f | B α p,q := ∞ 0 (t -α w r,p (f, t)) q dt t 1 q (q < ∞), sup t>0 t -α w r,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 ∈ L p (Ω) | f B α p,q < ∞}.

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; Giné & Nickl, 2015) and approximation theory (Temlyakov, 1993a) . These spaces are closely related to each other as follows (Triebel, 1983) :

EQUATION

Ω). In particular, under the same condition, from the definition of • B s p,q , it holds that B s p,q (Ω) → L r (Ω).

(1)

EQUATION

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 B 1 1,1 ([0, 1]) is included in the space of bounded total variation and B 1 1,∞ ([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 f o is lower bounded by n -2s/(2s+d) , (Kerkyacharian & Picard, 1992; Donoho et al., 1996; Donoho & Johnstone, 1998; Giné & 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,i h (f )(x) = ∆ r h (f (x 1 , . . . , x i-1 , •, x i+1 , . . . , x d ))(x i ) for f : R d → 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 ∈ R d is defined as ∆ r,e h (f ) = i∈e ∆ r,i hi (f ), ∆ r,∅ h (f ) = f. Then, the mixed modulus of smoothness is defined as

w e r,p (f, t) := sup |hi|≤ti,i∈e ∆ r,e h (f ) p for t ∈ R d + 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 ≤ ∞, α ∈ R d ++ and r i :=

α i + 1, the semi-norm | • | MB α,e

p,q based on the mixed smoothness is defined by

|f | MB α,e p,q :=    Ω [( i∈e t -αi i )w e r,p (f, t)] q dt i∈e ti 1/q (0 < q < ∞), sup t∈Ω ( 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 α,e p,q ,

and thus MB α p,q (Ω) := {f ∈ L p (Ω) | 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 B s p,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 ∂x 1 , ∂f ∂x 2 , ∂ 2 f ∂x 2 1 , ∂ 2 f ∂x 2 2 , ∂ 2 f ∂x 1 ∂x 2 , ∂ 3 f ∂x 1 ∂x 2 2 , ∂ 3 f ∂x 2 1 ∂x 2 , ∂ 4 f ∂x 2 1 ∂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 ∂x 1 , ∂f ∂x 2 , ∂ 2 f ∂x 2 1 , ∂ 2 f ∂x 2 2 , ∂ 2 f ∂x 1 ∂x 2 ,

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 (Dũ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

C d ω | 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 MB s 2,2 (Ω) is equivalent to the tensor product of Sobolev space (Sickel & Ullrich, 2011) which consists of functions f : Ω → R satisfying

C d d i=1 (1 + |ω i | 2 ) s | f (ω)| 2 dω < ∞.

Therefore, our analysis gives a (similar but) different characterization of conditions to avoid curse of dimensionality.

3. APPROXIMATION ERROR ANALYSIS:

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) ∈ R 1×W , b (L) ∈ R, W (1) ∈ R W ×d , b (1) ∈ R W , W ( ) ∈ R W ×W , b ( ) ∈ R W (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

f o ∈ MB α p,q (Ω) by an element f ∈ Φ(L, W, S, B) with precision > 0 measured by L r -norm: min f ∈Φ f -f o r ≤ .

3.1. APPROXIMATION ERROR ANALYSIS FOR BESOV SPACES:

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) ; Bö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. Let N (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 :

N m (x) = (N * N * • • • * N m + 1 times )(x), where f * g(x) := f (x -t)g(t)dt. It is known that N m is a piece-wise polynomial of order m. For k = (k 1 , . . . , k d ) ∈ Z d + and j = (j 1 , . . . , j d ) ∈ Z d + , let M d k,j (x) = d i=1 N m (2 ki x i -j i ). Even for k ∈ Z + , we also use the same notation to express M d k,j (x) = d i=1 N m (2 k x i -j i ).

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 M m k,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

EQUATION

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 W m p 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 ∈

Φ(L 0 , W 0 , S 0 , B 0 ) with L 0 := 3+2 log 2 3 d∨m c (d,m) + 5 log 2 (d ∨ m) , W 0 := 6dm(m+2)+2d, S 0 := L 0 W 2 0 and B 0 := 2(m + 1) m that satisfies M d 0,0 -M L ∞ (R d ) ≤ , and M (x) = 0 for all x ∈ [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; De-Vore et al., 1993; Dũ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 R r (F, H) := sup

f o ∈U (H) inf f ∈F f o -f L r ([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 + 2 log 2 3 d∨m c (d,m) + 5 log 2 (d ∨ m) , W = N W 0 , S = [(L -1)W 2 0 + 1]N, B = O(N (ν -1 +d -1 )(1∨(d/p-s)+) ), then it holds that R r (Φ(L, W, S, B), B s p,q ([0, 1] d )) N -s/d .

Remark 1. By Eq. (1), the condition (4) indicates that f o ∈ B s p,q satisfies f o ∈ L r (Ω). If we set p = q = ∞ and r = ∞, then B s p,q (Ω) = C s (Ω) 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 = r. In particular, the theorem also supports non-continuous regime (s < d/p) in which L ∞ -convergence does no longer hold but instead L r -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 = 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 Bö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 Bö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 f o 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 best N -term approximation error (Kolmorogov width) of the Besov space is lower bounded as inf

S N ⊂B s p,q sup f ∈U (B s p,q ) inf f ∈S N f -f L r (Ω)    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 S N is any N -dimensional subspace of B s p,q (Romanyuk, 2009; Myronyuk, 2016; Vybá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 (Dũ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.

3.2. APPROXIMATION ERROR ANALYSIS FOR M-BESOV SPACE:

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 D k,d := 1 + d-1 k k 1 + k d-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 < ∞

EQUATION

EQUATION

)

The proof is given in Appendix C.3. It holds that N

2 K K (d-1) , which implies 2 -K N -1 log d-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 log s(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 (MB s p,q , L r ) := inf L N sup f ∈U (MB s p,q ) f -L N (f ) r

, where the infimum is taken over all linear oprators L N with rank N from MB s p,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 Dũ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 log d-1 (N )) s      (q ≤ 2 ≤ r ≤ p < ∞), (q ≤ 1, p = r = ∞), (1 < p = r ≤ 2, q ≤ r), (N -1 log d-1 (N )) s (log d-1 (N )) 1/r-1/q (1 < p = r ≤ 2, q > r), (N -1 log d-1 (N )) s (log d-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 log d-1 (N )) s+1/r-1/p (2 ≤ p, 2 ≤ q ≤ r), (N -1 log d-1 (N )) s+1/r-1/p (log d-1 (N )) (1/r-1/q) + (r ≤ 2).

Therefore, the approximation error given in Theorem 1 achieves the optimal linear width ((N -1 log d-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.

4. ESTIMATION ERROR ANALYSIS:

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:

y i = f o (x i ) + ξ i (i = 1, . . . , n), where x i ∼ P X with density 0 ≤ p(x) < R on [0, 1] d , and ξ i ∼ N (0, σ 2 ). The data D n = (x i , y i ) 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 (y i -f (x i )) 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 L 1 -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.

4.1. ESTIMATION ERROR IN BESOV SPACES:

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 f o ∈ B s p,q (Ω) ∩ L ∞ (Ω) and. f o B s p,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 E Dn [ f o -f 2 L 2 (P X ) ] n -2s 2s+d log(n) 3 ,

where E Dn [•] indicates the expectation w.r.t. the training data D n .

The proof is given in Appendix D. The condition

f o ∞ ≤ F is required to connect the empirical L 2 -norm 1 n n i=1 ( f (x i ) -f o (x i )) 2 to the population L 2 -norm f -f o 2 L 2 (P X )

. 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; Giné & 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 (y i ) n i=1 linearly (it can be non-linearly dependent on (x i ) n i=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) = k x,X (k XX + λI) -1 Y ).

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 P X 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 f o ∈U (B s p,q ) E Dn [ f o -f 2 L 2 (P X ) ] n -2s-v 2s+1-v

where v = 2/(p ∧ 2)-1 and f runs over all linear estimators, that is, f depends on (y i ) n i=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 (B s p,q ) is strictly larger than the original set (Donoho & Johnstone, 1998) .

4.2. ESTIMATION ERROR IN MIXED SMOOTH BESOV SPACES:

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 f o ∈ MB s p,q (Ω) ∩ L ∞ (Ω) and f o MB s p,q ≤ 1 and f o ∞ ≤ F

for F ≥ 1, then letting (W, L, S, B) be as in Theorem 1, we obtain

E Dn [ f o -f 2 L 2 (P X ) ] n -2s 2s+1 log(n) 2(d-1)(u+s) 1+2s log(n) 3 .

Under the same assumption, if s > u log 2 (e) is additionally satisfied, we also have

E Dn [ f o -f 2 L 2 (P X ) ] n - 2s-2u log 2 (e)

2s+1+(1-2u) log 2 (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+log 2 (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 Hölder space with s < 2. However, our analysis is applicable to the m-Besov space which is more general than the mixed smooth Hö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 P X 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 f o ∈U (MB s p,q ) E Dn [ f -f o 2 L 2 (P X ) ] ≥ C 1 n -2s 2s+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 (x i , y i ) n i=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.

5. CONCLUSION:

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.

A PROOF OF LEMMA 1

Proof of Lemma 1. First note that N m (x) = 1 m! 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 of N m (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 = log 2 3 D + 5 log 2 (D) , W = 6d, S = LW 2 and B = 1 such that sup x∈[0,1] D φ mult (x 1 , . . . , x D ) - D i=1 x i ≤ ,

and φ mult (0, . . . , 0) = 0 for y ∈ R D such that D j=1 y j = 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, N m (x) = 0 for x ∈ [0, m + 1] gives that N m (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-j m+1 )

x -j m + 1 , . . . , φ (0,1-j m+1 )

x -j m + 1 m-times , we have that f (x) = 0 for all x ≤ 0 and

sup 0≤x≤m+1 |N m (x) -f (x)| ≤ 1 m! m+1 j=0 m + 1 j (m + 1) m ≤ (m + 1) m √ 2πm m+1/2 e -m 2 m+1 ≤ e (2e) m √ m =: ,

where we used m+1 j=0 m+1 j = 2 m+1 and Stirling's approximation m! ≥ √ 2πm m+1/2 e -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-j m+1 ) m + 1 -j m + 1 , . . . , φ (0,1-j m+1 ) m + 1 -j m + 1 =: δ .

It holds that |δ | ≤ . Since it is possible that δ = 0, we modify f (x) so that f (x) = 0 outside the interval [0, m + 1]. Because of this and noting 0 ≤ N m (x) ≤ 1, we see that g(x) := φ (0,1) (f (x) -

δ m+1 φ (0,m+1) (x)) yields sup x∈R |N m (x) -g(x)| ≤ 2 , sup x∈R |g(x)| ≤ 1, g(x) = 0 (∀x ∈ [0, m + 1]).

Hence, by applying φ mult again, we finally obtain that

sup x∈[0,1] d |M d 0,0 (x) -φ mult (g(x 1 ), . . . , g(x d ))| ≤ sup x∈[0,1] d M d 0,0 (x) - d j=1 g(x j ) + sup x∈[0,1] d d j=1

g(x j ) -φ mult (g(x 1 ), . . . , g(x d ))

≤2d + .

We again applying φ (0,1) , we obtain that h(x) = φ (0,1) • φ mult (g(x 1 ), . . . , g( 

x d )) satisfies M d 0,0 - h L ∞ (R d ) ≤ 2d + , h(x) = 0 for all x ∈ [0, m+1] d ,

B PROOF OF PROPOSITION 1:

For the order m ∈ N of the cardinal B-spline bases, let J(k) = {-m, -m + 1, . . . , 2 k -1, 2 k } d and the quasi-norm of the coefficient (α k,j ) k,j for k ∈ Z + and j ∈ J(k) be

(α k,j ) k,j b s p,q . =    k∈Z+   2 k(s-d/p) j∈J(k) |α k,j | p 1/p   q    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 ∈ B s p,q (Ω), there exists

f N that satisfies f -f N L r (Ω) N -s/d f B s p,q (9)

for N 1, and has the following form:

EQUATION

where

(j i ) n k i=1 ⊂ J(k), K = C 1 log(N )/d , K * = log(λN )ν -1 + K + 1, n k = λN 2 -ν(k-K) (k = K + 1, . . . , K * ) for δ = d(1/p -1/r) + and ν = (s -δ)/(2δ)

where the real number constants C 1 > 0 and λ > 0 are chosen to satisfy

K k=1 (2 k + m) d + K * k=K+1 n k ≤ N independently to N . Moreover, we can choose the coefficients (α k,j ) to satisfy (α k,j ) k,j b s p,q f B s p,q .

Proof of Lemma 2. DeVore & Popov (1988) constructed a linear bounded operator P k having the following form:

EQUATION

where α k,j is constructed in a certain way, where for every

f ∈ L p ([0, 1] d ) with 0 < p ≤ ∞, it holds f -P k (f ) L p ≤ Cw r,p (f, 2 -k ). (12) Let p k (f ) := P k (f ) -P k-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 B s p,q if and only if f can be decomposed into

f = ∞ k=0 p k (f ), with the convergence condition (p k (f )) ∞ k=0 b s p (L p ) := ( k∈Z+ (2 sk p k L p ) q ) 1/q < ∞; in particular, f B s p,q (p k (f )) ∞ k=0 b s p (L p ) .

Here, each p k can be expressed as

p k (x) = j∈J(k) α k,j M d k,j (x) for a coefficient (α k,j

) k,j which could be different from (a k,j ) k,j appearing in Eq. ( 11). Hence, f ∈ B s p,q can be decomposed into

EQUATION

with convergence in the sence of L p . Moreover, it is shown that

EQUATION

Based on this decomposition, Dũng (2011b) proposed an optimal adaptive recovery method such that the approximator has the form (10) under the conditions for K, K * , n k given in the statement and satisfies the approximation accuracy ( 9). This can be proven by applying the proof of Theorem 3.1 in Dũng (2011b) to the decomposition (13) instead of Eq. (3.8) of that paper. See also Theorem 5.4 of Dũ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 f N (10). Let the set of indexes (k, j) ∈ Z × Z that consists f N given in Eq. ( 10) be E N ; i.e.,

f N = (k,j)∈E N α k,j M d k,j . Accord- ingly, we set f := (k,j)∈E N α k,j M d k,j . For each x ∈ R d , it holds that |f N (x) -f (x)| ≤ (k,j)∈E N |α k,j ||M d k,j (x) -M d k,j (x)| ≤ (k,j)∈E N |α k,j |1{M d k,j (x) = 0} ≤ (m + 1) d (1 + K * )2 K * (d/p-s)+ f B s p,q log(N )N (ν -1 +d -1 )(d/p-s)+ f B s p,q

, where we used the definition of K * in the last inequality. Therefore, for each

f ∈ U (B s p,q ([0, 1] d )), it holds that f -f L r f -f N L r + f N -f L r log(N )N (ν -1 +d -1 )(d/p-s)+ f B s p,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 | 2 k(d/p-s)+ f B s p,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 d k,j are evaluated. Remembering that M m k,j (x) = M (2 k x 1j 1 , . . . , 2 k x d -j d ), we see that we just need to bound the quantity 2 k (k ≤ K * ). However, this is bounded by

2 k ≤ N ν -1 +d -1

for k ≤ K * . Hence, we obtain the assertion.

C PROOF OF THEOREM 1:

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 M m k,j (x) on a "regular grid" for k = 1, . . . , K and (j 1 , . . . , j d ) ∈ J(k), and take its superposition as f (x) ≈ k=1,...,K j∈J(k) α k,j M m k,j (x), which consists of O(2 Kd ) terms (see Eq. ( 10)). Hence, the number of parameters O(2 Kd ) 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 (Dũng, 1990; 1991; 1992; Temlyakov, 1982; 1993a; b) . Dũ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 Dũ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 ∈ Z d + and j ∈ J d m (k) := {-m, -m + 1, . . . , 2 k1 -1, 2 k1 } × • • • × {-m, -m + 1, . . . , 2 k d -1, 2 k d } as (α k,j ) k,j mb s p,q :=   k∈Z d +   2 (s-1/p) k 1 j∈J d m (k) |α k,j | p 1/p   q   1/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 ∈ MB s p,q (Ω) and K > 0, there exists R K (f ) such that R K (f ) can be represented as

R K (f )(x) = k∈Z d + : k 1 ≤K j∈J d m (k) α k,j M d k,j (x) + k∈Z d + : K< k 1 ≤K * n k i=1 α k,j (k) i M d k,j (k) i (x), where K * = K(1 + 2δ s-δ ) , (j (k) i ) n k i=1 ⊂ J d m (k), and n k = 2 K-s-δ 2δ ( k 1-K)

, and satisfies the following properties:

(i) For p ≥ r, f -R K (f ) r 2 -Ks D (1/ min(r,1)-1/q)+ K,d f MB s p,q . (ii) For p < r, f -R K (f ) r 2 -Ks D (1/r-1/q)+ K,d f MB s p,q (r < ∞), 2 -Ks D (1-1/q)+ K,d f MB s p,q (r = ∞).

Moreover, the coefficients (α k,j ) k,j can be taken to hold (α k,j ) k,j mb s p,q f MB s p,q .

The proof is given in Appendix C.2. The total number of cardinal B-spline bases consisting of R K (f ) can be evaluated as 17)).

2 K+1 K + d -1 d -1 + k:K< k 1≤K * n k 2 K D K,d + 2 K D K * ,d 2 K D K,d (∵ Eq. (

Here, D K,d can be evaluated as

D K,d K d-1 or D K,d d K .

Therefore, the total number of bases can be evaluated as

2 K min{K d-1 , d K }

which is much smaller than 2 Kd 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 R K (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(2 K K d-1 ) while the full grid takes O(2 Kd ) 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 (Dũng et al., 2016) ). An analogous result is already given by Dũ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 Dũng (2011a) . Because of this, our bound is better than the optimal rate of linear methods (Galeev, 1996; Romanyuk, 2001) and non-adaptive methods (Dũng, 1990; 1991; 1992; Temlyakov, 1982; 1993a; b) especially in the regime of p < r (Dũ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 = (k 1 , . . . , k d ) ∈ Z d + , let P (i)

ki f (x) be the function operating P ki defined in (11) to f as a function of x i with other components x j (j = i) fixed, and let

EQUATION

Then, p k can be expressed as

p k (x) = j∈J d m (k) α k,j M d k,j (x).

Let T (i)

ki = I -P (i)

ki and f p,i be the L p -norm of f as a function of x i with other components x j (j = i) fixed (i.e., if p < ∞, f p p,i = |f (x)| p dx i ), then Eq. ( 12) gives

T (i) ki f p,i sup |hi|≤2 -k i ∆ r,i hi (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 -k i ∆ r,i hi (T kj f ) p,i p,j = sup |hi|≤2 -k i T kj ∆ r,i hi (f ) p,j p,i (∵ the definition of ∆ r,i hi and Fubini's theorem) sup |hi|≤2 -k i sup |hj |≤2 -k j ∆ r,i hi (∆ r,j

hj (f )) p,j p,i , for i = j. Thus, applying the same argument recursively, for u ⊂

[d], it holds that i∈u T (i) ki f p w u r,p (f, 2 -k ), for k ∈ Z d + (u). Therefore, since p k = d i=1 (T (i) ki-1 -T (i) ki )f = u⊂[d] (-1) |u| i∈u T (i) ki i ∈u T (i)

ki-1 f , by letting e = {i | k i > 0}, we have that

p k p u⊂[d]   i∈u T (i) ki i ∈u T (i) ki-1   f p u⊂[d] w û r,p (f, 2 -(k u ) û ) e⊂u w u r,p (f, 2 -k u )

where

k u i := k i (i ∈ u) and k u i := k i -1 (i ∈ u), û = {i | k u i ≥ 0}, and (k u ) û is a vector such that (k u ) û,i = k u i for i ∈ û and (k u ) û,i = 0 for i ∈ û. Now let (p k ) k b α q (L p ) = k∈Z d + 2 α k 1 p k L p q 1/q for p k ∈ L p (Ω) (k ∈ Z d + ). Hence, if we set a k = e⊂u w u r,p (f, 2 -k u ) for k ∈ Z d + and e = {i | k i > 0}, we have that (p k ) k b α q (L p ) (a k ) b α q (L p )

f MB s p,q . On the other hand, following the same line of Theorem 2.1 (ii) of Dũng (2011a) , we also obtain the opposite inequality f MB s p,q

(a k ) b α q (L p ) (p k ) k b α

q (L p ) (note that the analogous inequality to Lemma 2.3 of Dũng (2011a) also holds in our setting by replacing q s with p s and ω e r (f, 2 -k ) p by w e r,p ). Therefore, f ∈ MB α p,q if and only if (p k ) k∈Z d + given by Eq. ( 15) satisfies (p k ) k b α q (L p ) < ∞ and f can be decomposed into f = k∈Z d + p k where convergence is in MB α p,q . Moreover, it holds that f MB α p,q (p k ) k b α q (L p ) . This can be shown by Theorem 2.1 of Dũng (2011a) . Moreover, by the quasi-norm equivalence p k p 2 -k 1/p ( j∈J d m (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 Dũ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

R K (f )(x) = k∈Z d + : k 1≤K p k + k∈Z d + :K< k 1≤K * G k (p k ),

where

G k (p k ) is given as G k (p k ) = 1≤i≤n k α k,ji M d k,ji (x)

where (α k,ji )

|J d m (k)| i=1

is the sorted coefficients in decreasing order of their absolute value:

|α k,j1 | ≥ |α k,j2 | ≥ • • • ≥ |α k,j |J d m (k)| |. Then, it holds that p k -G k (p k ) r ≤ p k p 2 δ k 1 n -δ

k , where δ := (1/p -1/r) (see the proof of Theorem 3.1 of Dũng (2011b) and Lemma 5.3 of Dũng (2011a) ). Moreover, we also have

p k r ≤ p k p 2 δ k 1 for k ∈ Z d + with k 1 > K * .

Here, we define N as N = log 2 (K) . Let = (s -δ)/(2δ), and

K * = K(1 + 1/ ) ,

and Then, by Lemma 5.3 of Dũng (2011a) , we have that

n k = 2 K-( k 1-K) for k ∈ Z d + with K + 1 ≤ k 1 ≤ K * .

EQUATION

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 =D k,d ≤ D k,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

EQUATION

(a) Suppose that q ≤ r and r < ∞. Then

f -R K (f ) q L r = f -R K (f ) r q r L r    K< k 1≤K * [2 δ k 1 n -δ k p k L p ] r + K * < k 1 [2 δ k 1 p k L p ] r    q r

(∵ Eq. ( 16))

K< k 1 ≤K * [2 δ k 1 n -δ k p k L p ] q + K * < k 1 [2 δ k 1 p k L p ] q ≤ N -δq 2 -(s-δ)Kq K< k 1 ≤K * [2 -(s-δ-δ )( k 1-K) ≤1 2 s k 1 p k L p ] q + 2 -q(s-δ)K * K * < k 1 [2 s k 1 p k L p ] q (N -δ 2 -(s-δ)K + 2 -(s-δ)K * ) q f q MB s p,q ≤ (N -s ) q f q MB α 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 -R K (f ) r L r K< k 1≤K * [2 δ k 1 n -δ k p k L p ] r + K * < k 1 [2 δ k 1 p k L p ] r (∵ Eq. (16)) ≤ N -δr 2 -(s-δ)Kr K< k 1≤K * [2 -(s-δ-δ )( k 1-K) 2 s k 1 p k L p ] r + K * < k 1 [2 s k 1 p k L p ] r (2 -(s-δ) k 1 ) r ≤ (N -δ 2 -(s-δ)K + 2 -(s-δ)K * ) r K< k 1≤K * [2 -(s-δ-δ )( k 1-K) 2 s k 1 p k L p ] r + K * < k 1 [2 s k 1 p k L p ] r 2 -(s-δ)( k 1 -K * )r ≤ (N -δ 2 -(s-δ)K + 2 -(s-δ)K * ) r    K< k 1≤K * [2 s k 1 p k L p ] rν + K * < k 1 [2 s k 1 p k L p ] 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 r MB s p,q D r(1/r-1/q) K,d

(∵ Eq. ( 17))

(N -s D 1/r-1/q K,d ) r f r MB s p,q .

(c) Suppose that r = ∞. Then, similarly to the analysis in (b), we can evaluate

f -R K (f ) L r N -δ 2 -(s-δ)K K< k 1≤K * [2 -(s-δ-δ )( k 1-K) 2 s k 1 p k L p ] + K * < k 1 [2 s k 1 p k L p ](2 -(s-δ) k 1 ) (N -δ 2 -(s-δ)K + 2 -(s-δ)K * )D (1-1/q)+ K,d f MB s p,q N -s D (1-1/q)+ K,d f MB s p,q . C.3 PROOF OF THEOREM 1 Let Z d + (e) := {k ∈ Z d + | k i = 0, i ∈ e} and for k ∈ Z d + (e), we define 2 -k := (2 -ki 1 , . . . , 2 -ki |e| ) ∈ R |e| + where {i 1 , . . . , i |e| } = e. By defining (g k ) k b s,e q := k∈Z d + (e) (2 s k 1 |g k |) q 1/q

for a sequence (g k ) k∈Z d + (e) , then it holds that |f | MB s,e p,q e⊂{1,...,d}

(w e r,p (f, 2 -k )) k b s,e q .

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 R K be E

K : R K (f ) = (k,j)∈E K α k,j M d k,j (x).

As in the proof of Proposition 1, we approximate R K (f ) by a neural network given as

f (x) = (k,j)∈E K α k,j M d k,j (x). Each coefficient α j,k satisfies |α j,k | 2 k 1(1/p-s)+ f MB s p,q

2 K * (1/p-s)+ . The difference between R K (f ) and f can be evaluated as

|R K (f ) -f (x)| ≤ (k,j)∈E K |α k,j ||M d k,j (x) -M d k,j (x)| ≤ (k,j)∈E K |α k,j |1{M d k,j (x) = 0} (m + 1) d (1 + K * )D K * ,d 2 K * (1/p-s)+ f MB s p,q . Therefore, by taking so that (m + 1) d (1 + K * )D K * ,d 2 K * (1/p-s)+ ≤ 2 -Ks is satisfied, it holds that |R K (f ) -f (x)| 2 -Ks . By the inequality D K * ,d ≤ e K * +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 * n k ≤2 K+1 K + d -1 d -1 + K<κ≤K * 2 K-s-δ 2δ (κ-K) κ + d -1 d -1 ≤2 K+1 D K,d + 2 K (1 -2 -s-δ 2δ ) -1 D K * ,d ≤ (2 + (1 -2 -s-δ 2δ ) -1 )2 K D K * ,d = N.

Since each unit M d k,j requires width W 0 , the whole width becomes W = N W 0 . The number of nonzero parameters to construct M d k,j is bounded by S = (L -1)W 2 0 N + N .

Finally, the magnitudes of the coefficients hidden in M d k,j are evaluated. Remembering that M d k,j (x) = M (2 k1 x 1 -j 1 , . . . , 2 k d x d -j d ), here maximum of 2 kj is bounded by 2 K * N (1+1/ν) . Hence, we obtain the assertion. Similarly, it holds that |α j,k | N (1+1/ν){1∨(1/p-s)+} .

D PROOFS OF THEOREMS 2 AND 3:

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

log N (δ, F, • ∞ ) N log(N )[log(N ) 2 + log(δ -1 )].

From Proposition 1, it holds that

f o -R K (f o ) 2 N -s/d . Note that f -f o 2 L 2 (P X ) f -f o 2 2 . for any f : [0, 1] d → R because p(x) ≤ R.

Therefore, by applying Proposition 4 with δ = 1/n, we have that

EQUATION

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. 

N (δ, F, • ∞ ) 2 K D K,d [K 2 log(2 K D K,d ) + log(δ -1 )] 2 K D K,d (K 3 + log(1/δ)). Letting r = 2, we have that f o -R K (f o ) 2 2 -sK D u K,d

where u = (1 -1/q) + for p ≥ 2 and u = (1/2 -1/q) + for p < 2.

Then, by noting that

f -f o 2 L 2 (P X ) f -f o 2 2 , for any f : [0, 1] d → R,

and by applying Proposition 4 with δ = 1/n, we have that

EQUATION

)

Here, we use the following evaluations for D 19) gives an upper bound of 2 -2sK e 2uK + 2 K e K (K 2 + log(n)) n .

Then, the right hand side is minimized by K = This gives the assertion.

E MINIMAX OPTIMALITY:

Proof of Theorem 4. First note that since P X is the uniform distribution, it holds that • L 2 (P X ) = • L 2 ([0,1] d ) . The -covering number N ( , G, L 2 (P X )) with respect to L 2 (P X ) for a function class G is the minimal number of balls with radius measured by L 2 (P X )-norm needed to cover the set G (van der Vaart & Wellner, 1996) . The δ-packing number M(δ, G, L 2 (P X )) of a function class G with respect to L 2 (P X ) norm is the largest number of functions {f 1 , . . . , f M } ⊆ G such that f i -f j L 2 (P X ) ≥ δ for all i = j. It is easily checked that N (δ/2, G, L 2 (P X )) ≤ M(δ, G, L 2 (P X )) ≤ N (δ, G, L 2 (P X )).

(20)

For a given δ n > 0 and ε n > 0, let Q be the δ n packing number M(δ n , U (MB s p,q ), L 2 (P X )) of U (MB s p,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): . 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 Dũng et al. (2016) ). Hence, it suffices to take

inf f sup f * ∈U (MB s p,q ) E Dn [ f -f * 2 L 2 (P X ) ] ≥ inf f sup f * ∈U (MB s p,q ) δ 2 n 2 P [ f -f * 2 L 2 (P X ) ≥ δ 2 n /

EQUATION

ε n δ n n -2s 2s+1 log(n)

2(d-1)(s+1/2-1/q) + 2s+1 ,

which gives the assertion.

F AUXILIARY LEMMAS:

Let the -covering number with respect to L 2 (P X ) for a function class G be N ( , G, L 2 (P X )) 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 (y i -f (x i )) 2 .

Assume that f o ∞ ≤ F and all f ∈ F satisfies f ∞ ≤ F for some F ≥ 1. If δ > 0 satisfies N (δ, F, • ∞ ) ≥ 3, then it holds that

E Dn [ f -f o 2 L 2 (P X ) ] ≤ C inf f ∈F f -f o 2 L 2 (P X ) + (F 2 + σ 2 ) log N (δ, 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 has a parameters with distance δ: W ( ) -W ( ) ∞ ≤ δ and b ( ) -b ( ) ∞ ≤ δ. Now, not that A k (f ) ∞ ≤ max j W (k-1) j,:

1 A k-1 (f ) ∞ + b (k-1) ∞ ≤ W B A k-1 (f ) ∞ + B ≤ (B ∨ 1)(W + 1) A k-1 (f ) ∞ ≤ (B ∨ 1) k-1 (W + 1) k-1

, and similarly the Lipshitz continuity of B k (f ) with respect to • ∞ -norm is bounded as (BW ) L-k+1 . Then, it holds that

|f (x) -g(x)| = L k=1 B k+1 (g) • (W (k) • +b (k) ) • A k (f )(x) -B k+1 (g) • (W (k) • +b (k) ) • A k (f )(x) ≤ L k=1 (BW ) L-k (W (k) • +b (k) ) • A k (f )(x) -(W (k) • +b (k) ) • A k (f )(x) ∞ ≤ L k=1 (BW ) L-k δ[W (B ∨ 1) k-1 (W + 1) k-1 + 1]

:

In 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.

Since Ω = [0, 1] d in our setting, the boundedness automatically follows from the continuity.

We let N := {1, 2, 3, . . . }, Z+ := {0, 1, 2, 3, . . . }, Z d + := {(z1, . . . , z d ) | zi ∈ Z+}, R+ := {x ≥ 0 | x ∈ R}, and R++ := {x > 0 | x ∈ R}.

The notation → means a continuous embedding; that is, X → Y for two normed spaces X, Y means X is continuously embedded in Y .

We 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.

Figures and Tables:

and h ∞ ≤ 1. Finally, by carefully checking the network construction, it is shown that h ∈ Φ(L, W, S, B) with L = 3 + 2 log 3 d∨m + 5 log 2 (d ∨ m) , 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 .

Proof of Theorem 3. The proof follows the almost same line as the proof of Theorem 2. By noting S = O(2 K D K,d ), L = O(K) and W = O(2 K D K,d ), Lemma 3 gives an upper bound of the covering number as log

K,d : (a) D K,d K d-1 , and (b) D K,d [e(1 + d K )] K . (a) For the evaluation, D K,d K d-1 , we have an upper bound of the right hand side of Eq. (19) as2 -2sK K 2u(d-1) + 2 K K d-1 (K 3 + log(n)) n ,which is minimized by setting K = 1 1+2s log 2 (n) + (2u-1)(d-1) 1+2s log 2 log(n) 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, D K,d [e(1 + d K )] K ≤ e K e d , Eq. (

2u) log 2 (e) log 2 (n) . Then, we have that n 2u) log 2 (e) log(n) 2 .

)5 . The only difference is the assumption of f ∞ ≤ F for f ∈ F and f = f o 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 bylog N (δ, Φ(L, W, S, B), • ∞ ) ≤ S log(δ -1 L(B ∨ 1) L-1 (W + 1) 2L ) ≤ 2SL log((B ∨ 1)(W + 1)) + S log(δ -1 L). 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 A k (f )(x) = η • (W (k-1) η(•) + b (k-1) ) • • • • • (W (1) x + b (1) ),andB k (f )(x) = (W (L) η(•) + b (L) ) • • • • • (W (k) η(x) + b (k) ),for k = 2, . . . , L. Corresponding to the last and first layer, we define B L+1 (f )(x) = x and A 1 (f )(x) = x. Then, it is easy to see that f(x) = B k+1 (f ) • (W (k) • +b (k) ) • A k (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) ),

Summary of relation between related existing work and our work for a mixed smooth Besov space. N is the number of parameters in the deep neural network, n is the sample size. β represents the smoothness parameter, and d represents the dimensionality of the input. The approximation accuracy is measured by L 2 -norm and estimation accuracy is measured by the square of L 2 -norm. See Theorem 3 for the definition of u.

Table content: <table/>

Peer review (iclr2019):
 Summary: ======== The paper presents rates of convergence for estimating nonparametric functions in Besov spaces using deep NNs with ReLu activations. The authors show that deep Relu networks, unlike linear smoothers, can achieve minimax optimality. Moreover, they show that in a restricted class of functions called mixed Besov spaces, there is significantly milder dependence on dimensionality. Even more interestingly, the Relu network is able to adapt to the smoothness of the problem. While I am not too well versed on the background material, my educated guess is that the results are interesting and relevant, and that the analysis is technically sound. Detailed Comments: ================== My main criticism is that the total rate of convergence (estimation error + approximation error) has not been presented in a transparent way. The estimation error takes the form of many similar results in nonparametric statistics, but the approximation error is given in terms of the parameters of the network, which depends opaquely on the dimension and other smoothness parameters. It is not clear which of these terms dominate, and consequently, how the parameters W, L etc. should be chosen so as to balance them. While the mixed Besov spaces enables better bounds, the condition appears quite strong. In fact, the lower bound is better than for traditional Holder/Sobolev classes. Can you please comment on how th m-Besov space compares to Holder/Sobolev classes? Also, can you similiarly define mixed Holder/Sobolev spaces where traditional linear smoothers might achieve minimax optimal results? Minor: - Defn of Holder class: you can make this hold for integral beta if you define m to be the smallest integer less than beta (e.g. beta=7, m=6). Imo, this is standard in most texts I have seen. - The authors claim that the approximation error does not depend on the dimensionality needs clarification, since N clearly depends on the dimension. If I understand correctly, the approximation error is in fact becoming smaller with d for m-Besov spaces (since N is increasing with d), and what the authors meant was that the exponential dependnence on d has now been eliminated. Is this correct? Other - On page 4, what does the curly arrow notation mean? - Given the technical nature of the paper, the authors have done a good job with the presentation. However, in some places the discussion is very equation driven. For e.g. in the 2nd half of page 4, it might help to explain many of the quantities presented in plain words. Confidence: I am reasonably familiar with the nonparametric regression literature, but not very versed on the deep learning theory literature. I did not read the proofs in detail.