Alexandre B Tsybakov
email: tsybakov@ccr.jussieu.fr

Anatoli B Juditsky
email: anatoli.juditsky@imag.fr

Oleg V Lepski
email: lepski@cmi.univ-mrs.fr

Nonparametric estimation of composite functions

Keywords: multidimensional nonparametric estimation, minimax estimation, adaptive estimation, composite functions, single index model

Nonparametric estimation of composite functions

Introduction

In this paper we study the problem of nonparametric estimation of an unknown function g : R d → R in the multidimensional gaussian white noise model described by the stochastic differential equation

X ε (dt) = g(t)dt + εW (dt), t = (t 1 , . . . , t d ) ∈ D ( 1
)
where D is an open interval in R d containing [-1, 1] d , W is the standard Brownian sheet in R d and 0 < ε < 1 is a known noise level. Our goal is to estimate the function g on the set [-1, 1] d from the observation {X ε (t), t ∈ D}. For d = 2 this corresponds to the problem of image reconstruction from observations corrupted by additive noise. We consider observation set D which is larger than [-1, 1] d in order to avoid the discussion of boundary effects.

To measure the performance of estimators, we will use the risk function determined by the sup-norm • ∞ on [-1, 1] d : for g : R d → R, 0 < ε < 1, p > 0, and for an arbitrary 1 estimator gε based on the observation {X ε (t), t ∈ D} we consider the risk R ε gε , g = E g gεg p ∞ .

Here and in what follows E g denotes the expectation w.r.t. the distribution P g of the observation {X ε (t), t ∈ D} satisfying (1). We will suppose the g ∈ G s , where {G s , s ∈ S} is a collection of functional classes indexed by s ∈ S. The functional classes G s that we will consider consist of smooth composite functions and below we discuss in detail this choice.

For a given class G s we define the maximal risk

R ε gε , G s = sup g∈Gs R ε gε , g . (2)
Our first aim is to study the asymptotics, as the noise level ε tends to 0, of the minimax risk inf gε R ε gε , G s where inf gε denotes the infimum over all estimators of g. We suppose that parameter s is known, and therefore the functional class G s is fixed. In other words, we are interested in minimax estimation of g. We find the minimax rate of convergence φ ε (s) on G s , i.e., the rate which satisfies φ p ε (s) ≍ inf gε R ε gε , G s and we construct an estimator attaining this rate, called rate optimal estimator in asymptotic minimax sense. The estimator depends on parameter s which restricts its application in practice. We discuss approaches to treat adaptation to s and state some conjectures on this issue. We suggest a possible construction of such an adaptive procedure. Further details will be given in a forthcoming paper.

Motivation

It is well known that the main difficulty in estimation of multivariate functions is the curse of dimensionality: the best attainable rate of convergence of the estimators decreases very fast as the dimension grows. To illustrate this effect, suppose, for example, that the underlying function g belongs to G s = H d (α, L), s = (α, L), α > 0, L > 0, where H d (α, L) is an isotropic Hölder class of functions. We give the exact definition of this functional class later. Here we only mention that H d (α, L) consists of functions g with bounded partial derivatives of order ≤ ⌊α⌋ and such that, for all x, y ∈ D, |g(y) -P g (x, yx)| ≤ L xy α , where P g (x, yx) is the Taylor polynomial of order ≤ ⌊α⌋ obtained by expansion of g around the point x, and • is the Euclidean norm in R d . Parameter α characterizes the isotropic (i.e., the same in each direction) smoothness of function g.

If we use the risk [START_REF] Bertin | Asymptotically exact minimax estimation in sup-norm for anisotropic Hölder classes[END_REF], uniformly on H d (α, L) the rate of convergence of estimators cannot be asymptotically better than ψ ε,d (α) = ε ln (1/ε) 2α/(2α+d) (cf. (9; 10; 23; 20; 5)). This is also the minimax rate on H d (α, L): it is achieved, for example, by a kernel estimator with properly chosen bandwidth and kernel. More results on asymptotics of the minimax risks in estimation of multivariate functions can be found in (11; 12; 2; 3). It is clear that if α is fixed and d is large enough this asymptotics is too pessimistic to be used for real data: the value ψ ε,d (α) is small only if the noise level ε is unreasonably small. On the other hand, if the noise level ε is realistically small and α is small the above asymptotics can be of no use already in dimensions 2 or 3.

At the origin of this phenomenon is the fact that the d-dimensional isotropic Hölder class H d (α, L) is too massive in terms of its metric entropy. To "overcome" the curse of dimensionality one usually considers models with slimmer functional classes (i.e., classes with smaller metric entropy). There are several ways to do it.

• A first way is to impose a restriction on the smoothness parameter of the functional class. For the class H d (α, L), a convenient restriction is to assume that the smoothness α increases with the dimension, and thus the class becomes smaller (its metric entropy decreases). For instance, we can suppose that α = κd with some fixed κ > 0.

Then the dimension disappears from the expression for ψ ε,d (α), which means that we escape from the curse of dimensionality. However, the condition α = κd or other similar restrictions linking smoothness and dimension are usually difficult to motivate. One interesting related example seems to be the class of functions with absolutely integrable multivariate Fourier transform [START_REF] Barron | Universal approximation bounds for superpositions of a sigmoidal function[END_REF].

• Another way is to impose a structural assumption on the function g to be estimated.

Two classical examples are provided by the single index and additive structures (cf., e.g., (24; 6; 8)).

The single index structure is defined by the following assumption on g: there exist a function F 0 : R → R and a vector ϑ ∈ R d with ϑ = 1 such that g(x) = F 0 (ϑ T x).

The additive structure is defined by the following assumption: there exist functions F i : R → R, i = 1, . . . , d, such that g(x) = F 1 (x 1 ) + • • • + F d (x d ), where x j is the j-th component of x ∈ R d .

If we suppose that F i ∈ H 1 (α, L), i = 0, . . . , d, then in both cases function g can be estimated with the rate ε ln (1/ε) 2α/(2α+1) which does not depend on the dimension and coincides with the minimax rate ψ ε,1 (α) of estimation of functions on R.

In general, under structural assumptions the rate of convergence of estimators improves, as compared to the slow d-dimensional rate ψ ε,d (α). For the above examples the rate does not depend on the dimension.

However, it is often quite restrictive to assume that g has some simple structure, such as the single index or additive one, on the whole domain of its definition. In what follows we refer to this assumption as global structure.

A more flexible way of modeling is to suppose that g has a local structure. For instance, we can assume that g is well approximated by some single index or additive structure (or by a combination both) in a small neighborhood of a given point x. Local structure depends on x and remains unchangeable within the neighborhood. Such an approach can be used to model much more complex objects than the global one. However, the form of the ddimensional neighborhood and the local structure should be chosen by the statistician in advance, which makes the local approach rather subjective.

In the present paper we try to find a compromise between the global and local modeling. Our idea is to consider a sufficiently general global model which would generate suitable local structures, and thus would allow us to construct estimators with nice statistical properties. Such a global model should satisfy the following requirements.

(i) The optimal performance of estimators under this model is better than the performance of methods based only on the smoothness properties of the underlying function.

(ii) The model contains "global" parameters such that their choice automatically generates interpretable local structures and associated domains of localization (neighborhoods). Adaptation to such global parameters, if it is feasible, means adaptation to different local models simultaneously.

(iii) The model contains a tuning parameter such that its choice allows one to reduce the influence of the dimension.

We argue that this program can be realized for global models where the underlying function g is a composition of two smooth functions.

Smooth composite functions

We now define our global structural model. We will assume that g is a composite function, i.e., that g(t) = f (G(t)) for all t ∈ R d where f : R → R and G : R d → R with d ≥ 2.

We will further suppose that f and G are smooth functions such that f ∈ H 1 (γ, L 1 ) and

G ∈ H d (β, L 2 )
= (γ, β) ∈ R 2 + and L = (L 1 , L 2 ) ∈ R 2 + .
This model is a generalization of the single index model: instead of the linear function we have here a general G(•).

The performance of an estimation procedure will be measured by the sup-norm risk (2) where we set s = (A, L) and G s = H(A, L). The global parameter of the model is s = (A, L), and we will show that the choice of A leads to different local structures. Note also that the value of s determines the quality of estimation associated to our model, i.e., the rate of convergence of the minimax risk.

We start with the following definitions:

Definition 1. Fix α > 0 and L > 0. Let ⌊α⌋ be the largest integer which is strictly less than α, and for k

= (k 1 , . . . , k d ) ∈ N d set | k| = k 1 + • • • + k d . The isotropic Hölder class H d (α, L)
is the set of all functions G : R d → R having on R d all partial derivatives of order ⌊α⌋ and such that

0≤| k|≤⌊α⌋ sup x∈R d ∂ | k| G(x) ∂x k 1 1 • • • ∂x k d d ≤ L, G(y) - 0≤| k|≤⌊α⌋ ∂ | k| G(x) ∂x k 1 1 • • • ∂x k d d d j=1 (y j -x j ) k j k j ! ≤ L y -x α , ∀ x, y ∈ R d . ( 3
)
where x j and y j are the jth components of x and y.

Definition 2. Let G be some class of functions on R d . We say that G is a class of effective smoothness α > 0 if G ⊆ H d (α, L) for some L > 0, and

G ⊆ H d (α ′ , L ′ ) for all α ′ > α, L ′ > 0.
Observe that for any α > 0, L > 0, the effective smoothness of the single index and additive classes is equal to α. Next, it is easy to see that for any

A = (γ, β) ∈ R 2 + and L = (L 1 , L 2 ) ∈ R 2 + the effective smoothness of H(A, L) is α γ,β γβ if 0 < γ, β ≤ 1, min(γ, β) otherwise. ( 4
)
We say that there is an improvement in statistical performance of estimators due to the structure if the minimax rate of convergence associated to the class G is o(ψ ε,d (α)), as ε → 0, where α is the effective smoothness of G. In other words, the knowledge of a structure allows certain improvement if the best estimator, based only on the smoothness properties, converges slower than the best estimator which takes into account the whole structure. For example, for the classes of functions with single-index or additive structure there is always an improvement due to the structure, because the corresponding minimax rate is

ε ln (1/ε) 2α/(2α+1) = o(ψ ε,d (α)).
For the class H(A, L) of composite functions this property is not always true. For certain values of A = (γ, β) no improvement due to the structure can be expected. This happens if our structural assumption is essentially equivalent to the fact that g belongs to some isotropic Hölder class of functions of full dimension d, and the knowledge of the composition structure does not help to improve the statistical analysis. Such an effect appears in the following two zones of (γ, β).

1 • . Zone of slow rate: 0 < γ, β ≤ 1. Clearly, in this zone H(A, L) ⊂ H d (γβ, L 3 )
, where L 3 is a positive constant depending only on γ, β and L. Due to this inclusion a standard kernel estimator with properly chosen bandwidth and the boxcar kernel converges with the rate

ψ ε,d (γβ) = ε ln (1/ε) 2γβ/(2γβ+d) .
It is not hard to see (cf. Section 8) that this rate is optimal, i.e., that a lower bound on the minimax risk holds with the same "slow" rate ψ ε,d (γβ) (note that γβ ≤ 1). As the effective smoothness of H(A, L) for 0 < γ, β ≤ 1 equals to γβ, there is no improvement due to the structure. 2 • . Zone of inactive structure: γ ≥ β, γ ≥ 1.

In this zone we easily get the inclusions

H d (β, L 4 ) ⊂ H(A, L) ⊂ H d (β, L 5 )
, where L 4 and L 5 are positive constants depending only on β and L. To show the left inclusion it suffices to fix a linear function f . Therefore, the asymptotics of the minimax risk on H(A, L) is the same as for any isotropic Hölder class H d (β, •). In particular, a standard kernel estimator converges with the rate ψ ε,d (β). Note that here we estimate as if there were no structure, and the asymptotics of the minimax risk does not depend on γ. This explains why we refer to this zone as that of inactive structure.

We finally remark that if β ≤ 1 the composite function g is rather nonsmooth. The effective smoothness equals to (1∧γ)β, and in view of the above discussion, the minimax rate of convergence of estimators on H(A, L) is the same as on the Hölder class H d ((1 ∧ γ)β, •). This is a very slow rate ψ ε,d ((1 ∧ γ)β). Therefore, only for β > 1 one can expect to find estimators with interesting statistical properties.

Main results

In this section we state the main results and outline the estimation method. The formal description of the estimation procedure and the proofs are deferred to Sections 6 and 8 -9 respectively.

Lower bound for the risks of arbitrary estimators

For any A = (γ, β) ∈ R 2 + define π(A) = 2γ 2γ + 1 + (d -1)/β ∧ 2 2 + d/β ∧ 2 2 + d/(γβ) , (5)
and

φ ε (γ, β) = ε ln(1/ε) π(A) . (6)
In an expanded form, we may write

φ ε (γ, β) =        ε ln (1/ε) 2γ 2γ+1+(d-1)/β if β > 1, β ≥ d(γ -1) + 1, ε (ln 1/ε) 2 2+d/β if γ > 1, β < d(γ -1) + 1, ε ln (1/ε) 2 2+d/(γβ) if (γ, β) ∈ (0, 1] 2 . ( 7
)
The boundaries between the zones of these three different rates in R 2 + are presented by the dashed lines in Figure 1.

An asymptotic lower bound for the minimax risk on H(A, L) is given by the following theorem.

Theorem 1. For any A = (γ, β) ∈ R 2 + and any p > 0 we have

lim inf ε→0 inf gε sup g∈H(A,L) E g φ -1 ε (γ, β) gε -g ∞ p > 0,
where inf gε denotes the infimum over all estimators of g.

This theorem shows that the rate of convergence φ ε (γ, β) cannot be improved on any estimators. We will next claim that for 0 < γ, β ≤ 2 there exist estimators attaining this rate. Before stating the corresponding result, we make some remarks on the properties of the rate φ ε (γ, β).

Remark 1. There is an improvement due to the structure everywhere except for the trivial cases 1 • and 2 • discussed in Section 3. This corresponds to the zone {A = (γ, β) : β > γ, β ≥ 1} which we refer to as zone of improved rate (cf. Figure 1). Indeed, when A belongs to this zone the effective smoothness is α γ,β = γ (cf. (4)), and hence φ ε (γ, β) = o(ψ ε,d (α γ,β )), as ε → 0. Thus, the requirement (i) of Section 2 is met.

Remark 2. Parameter β can be viewed as a tuning parameter of the model: its choice can reduce the impact of the dimension d on the accuracy of estimation. Indeed, as the ratio d/β tends to 0, the rate φ ε (γ, β) approaches either the one-dimensional Hölder rate ψ ε,1 (γ) or the "almost parametric" rate ε ln (1/ε). Thus, the requirement (iii) of Section 2 is met. β < d(γ -1) + 1 the rate of convergence φ ε (γ, β) does not depend on γ and coincides with the minimax rate ψ ε,d (β) associated to the d-dimensional Hölder class H d (β, •). This is rather surprising: in this zone the composite function g = f •G can be estimated with the same rate as G, independently of how smooth is f . Such a behavior cannot be explained in terms of the smoothness because in the considered case the effective smoothness α γ,β takes the value γ and not β (cf. (4)). is valid for all positive γ, β. Inspection of the proof shows that this bound is attained at the least favorable functions that, for the particular case of d = 2, are of the form f 0 (ϕ 1 (t 1 ) + ϕ 2 (t 2 )) where f 0 is a function of Hölder smoothness γ and both functions ϕ j are of Hölder smoothness β. So, for d = 2 the lower bound with rate ε ln (1/ε) 2γ 2γ+1+1/β holds for such more restricted class of functions, whatever are γ and β. In particular, when γ = β, this lower rate becomes ε ln (1/ε)

Remark 3. If β ≥ γ > 1,
2β 2 2β 2 +β+1 . Since 2β 2 2β 2 +β+1 < 2β
2β+1 this is always slower than the classical one-dimensional rate ε 2β 2β+1 . On the other hand, a recent result of Horowitz and Mammen [START_REF] Horowitz | Rate-optimal estimation for a general class of non-parametric regression models with unknown link function[END_REF] shows that for γ = β functions of the form f 0 (ϕ 1 (t 1 ) + ϕ 2 (t 2 )) can be estimated at the rate ε 2β 2β+1 in the L 2 -norm. This phenomenon is very surprising because, in contrast to classical nonparametric estimation problems, we observe here a significant (and not only a logarithmic) deterioration of the rate when passing from the L 2 -norm to the L ∞ -norm.

Outline of the estimation method

The exact definition of our estimator is given in Section 6. Here we only outline its construction. We suppose that A = (γ, β) ∈ (0, 2] 2 . The initial building block is a family of kernel estimators. In contrast to the classical kernel construction which involves a unique bandwidth parameter, the kernel K J that we consider is determined by the triplet J = (A, ϑ, λ) where the form parameter A is the couple (γ, β) ∈ (0, 2] 2 , the orientation parameter ϑ is a unit vector in R d and λ is a positive real which we refer to as size parameter. We denote J the set of all such triplets J and consider a family of kernel estimators (ĝ J , J ∈ J) where for any x ∈ [-1, 1] d the estimator ĝJ (x) of g(x) is given by ĝJ (x)

D K J (t -x)X ε (dt).
We will see that, in general, the size parameter λ is not equivalent to the bandwidth of classical kernel estimator. In fact, the value of λ characterizes the bias of the estimator ĝJ when the orientation of the window ϑ is locally "correct". Namely, the kernel K J is chosen in such a way that for each x ∈ [-1, 1] d the bias of ĝJ is of the order O(λ) if ϑ = ϑ x 0 is collinear to the gradient ∇G(x).

The estimation method proceeds in three steps, and the basic device underlying the construction of the optimal estimation method is the notion of the local model. It is an important feature of the composition structure that different local models arise in different subsets of the zone of improved rate.

Step 1. Specifying a collection of local models. The underlying function g with a complicated global structure can have a simple local structure. However, the local structure depends on the function itself. Therefore, g can be only described by a collection of local models. In our case, this collection is indexed by a finite-dimensional parameter which can be considered as a nuisance parameter. Specifically, we pass from the global composition model defined in Section 3 to a family of local models {M J (x), J ∈ J, x ∈ [-1, 1] d } where the type of each local model M J (x), J = (A, ϑ, λ), is determined by A, while ϑ and λ are the local orientation and size parameters. We will show that, depending on the value of A = (γ, β) (cf. Figure 2), our global model induces only two types of local models: a local single index model and a combined local model. The latter combines elements of both single index and additive models. This responds to the requirement (ii) of Section 2.

• . Local single index model: γ ≤ 1, 1 < β ≤ 2.
In this domain of γ, β, using the smoothness properties of functions f and G it is not hard to show that in the ball

B λ,x (A) = t ∈ R d : t -x ≤ λ 1 γβ the composite function g(•) can be approximated with the accuracy O(λ) by the function f G(x) + ϑ T [• -x] . Here ϑ = ϑ x
0 is a unit vector collinear to the gradient ∇G(x). Indeed, since the inner function G belongs to H d (β, L 2 ), for any x, y ∈ D we have

G(t) = G(x) + ∇G(x) T (t -x) + B x (t), with B x (t) ≤ L 2 t -x β . (8)
Next, using the fact that f ∈ H 1 (γ, L 1 ), we conclude that g(t) = f (G(t)) admits the representation

g(t) = Q x (t) + C x (t),
where

Q x (t) = f (G(x) + ∇G(x) T (t -x)) and |C x (t)| ≤ L 1 |B x (t)| γ ≤ L 1 L γ 2 t -x γβ .
In other words, for any kernel K with the support on the ball B λ (A) = t ∈ R d : t ≤ λ 1/γβ and such that K(y)dy = 1,

K(t -x)[g(t) -Q x (t)]dt = O(λ). (9)
We understand the relation [START_REF] Ibragimov | Statistical Estimation. Asymptotic Theory[END_REF] as the definition of the local single index model Q x of g. The choice of the approximation kernel for the function g is naturally suggested by the form of the local model Q x together with the bound (9): the kernel K J can be taken as the indicator function of a hyperrectangle normalized by its volume and oriented in such a way that ∇G(x) is collinear to the first basis vector in R d . The sides of the hyperrectangle are chosen to have the lengths

l 1 = λ 1 γ and l j = λ 1 γβ , j = 1, . . . , d -1. 2 • . Combined local model: 1 < γ ≤ β ≤ 2.
Let M ϑ be an orthogonal matrix with the first column equal to ϑ = ϑ x 0 , and let y = M T ϑ (tx), t ∈ R d . We denote with y j the jth component of y and consider the set

X λ,x (A) = t ∈ R d : |y 1 | ≤ λ 1 β , y ≤ λ 1 γβ , |y 1 | γ-1 y β ≤ λ . ( 10
)
We show that the estimation of the composite function g at x can be reduced to the problem of estimation under the local model

Q x (y) = q x (y 1 ) + P x (y 2 , . . . , y d ),
where

q x ∈ H 1 (γ, L 1 L γ 2 ) and P x ∈ H d-1 (β, 2L 1 L 2 ) on the set X λ,x (A)
. This local model is established in an unknown coordinate system determined by the parameter ϑ = ϑ x 0 . Since y 1 is the coordinate of the projection on ϑ x 0 , the component q x (y 1 ) constitutes an element of single index structure. An element of additive structure comes from the separation of Q x into the sum of two functions depending on non-intersecting sets of coordinates.

The explanation of the local model represented by Q x on the set X λ,x (A) is provided by the following argument. Using the smoothness properties of functions f and G, we obtain due to the inclusions

f ∈ H 1 (γ, L 1 ), G ∈ H d (β, L 2 ): g(t) = f (G(x) + ∇G(x) T (t -x)) + f ′ (G(x) + ∇G(x) T (t -x))B x (t) + C x (t) = f (G(x) + ∇G(x) T (t -x)) + f ′ (G(x))B x (t) + D x (t) + C x (t),
where

|C x (t)| ≤ C(L 1 , L 2 , γ) t -x γβ , |D(x, y)| ≤ C(L 1 , L 2 ) |∇G(x) T (t -x)| ∇G(x) t -x β ,
and the function B x (t), which is defined in [START_REF] Hristache | Direct estimation of the index coefficient in a single-index model[END_REF], belongs to the class H d (β, 2L 2 ). In the transformed coordinates (determined by the orthogonal matrix M ϑ ) we may write

g(t) = g(x + M ϑ y) = q(y 1 ) + Bx (y) + Dx (y) + Cx (y), (11)
where

| Dx (y) + Cx (y)| ≤ C(L 1 , L 2 , γ)(|y 1 | γ-1 y β + y γβ ). ( 12
)
and Bx ∈ H d (β, 2L 2 ). The latter inclusion leads to

Bx (y) -P x (y 2 , . . . , y d ) -y 1 ∂ ∂y 1 Bx (0, y 2 , . . . , y d ) ≤ 2L 2 |y 1 | β , (13)
where P x (y 2 , . . . , y d ) = Bx (0, y 2 , . . . , y d ). Let again K be a kernel such that K(t)dt = 1, supported on X λ,x (A). Then

K(y -x)[g(x + M ϑ y) -Q x (y)]dy = O(λ) ( 14
)
if K is symmetric in y 1 . We understand this property as the definition of the combined local model Q x for the composite function g. We conclude that if A belongs to the zone marked as "combined" in Figure 2, the global structural assumption that the underlying function is a composite one leads automatically to a local structure containing elements of both single index and additive models.

A good kernel K J for the zone of combined local model should be supported on the right window X λ,x (A), possess small bias on both single-index component q x and "regular" component P x and have a small L 2 -norm to ensure small variance of the stochastic term of the estimation error. Construction of such a kernel is a rather involved task (cf. Section 7.2). Using a rectangular kernel, as for the local single-index model, does not give a solution, since it leads to suboptimal estimation rates.

As we see, the definition of local model has two ingredients: the neighborhood (window) and the local structure within the window. For the local single index model the window is just an Euclidean ball, whereas for the combined local model the window is the set X λ,x (A) which has quite a nonstandard form (cf. Figure 3). Step 2.Optimizing the size parameter and specifying candidate estimators.

Once the local model is determined and the corresponding kernel is constructed we can chose the size parameter λ = λ ε (A) in an optimal way. To do it we optimize our sup-norm risk with respect to λ, i.e., we get the value λ which realizes the balance of bias and variance terms of the risk in the ideal case where the orientation ϑ = ϑ x 0 is "correct" for all x. Recall that the kernel K J supported on the window is chosen in such a way that the bias of the kernel estimator ĝJ , for the "correct" orientation ϑ, is of the order O(λ) on every local model. Thus, the bias-variance balance relation for the sup-norm loss can be written in the form

λ ≍ ε ln 1/ε K J 2 . ( 15
)
We will see that K J 2 depends on A and λ but does not depend on ϑ. This will allow us to choose the optimal value λ ε (A) independent of ϑ. For instance, for the local single index model the kernel K J is just a properly scaled and rotated indicator of a hyperrectangle. In this particular case the bias-variance balance (15) can be written in the form

λ ≍ ε ln 1/ε √ volume of hyperrectangle = ε ln 1/ε λ 1 γ + d-1 γβ 1/2 .
Note that in this case λ ε (A) ≍ φ ε (γ, β), where φ ε (γ, β) is defined in [START_REF] Horowitz | Rate-optimal estimation for a general class of non-parametric regression models with unknown link function[END_REF]. On the other hand, to guarantee that the same relation λ ε (A) ≍ φ ε (γ, β) holds in the zone of combined local model need a rather sophisticated construction of the kernel K J (cf. Section 7.2).

With λ ε (A) being chosen, we obtain a family of kernel estimators

ĝJ (x), J = (A, ϑ, λ ε (A)) ∈ J, x ∈ [-1, 1] d . ( 16
)
For a fixed x ∈ [-1, 1] d this family only depends on two parameters, A and ϑ.

Step 3. Aggregating the estimators. We now choose an estimator from the family ( 16) which corresponds to some Ĵ ∈ J selected in a data-dependent way, and define our final estimator as a piecewise-constant approximation of the function x → ĝ Ĵ (x). To choose Ĵ we apply an aggregation procedure which is a special case of the method of aggregation of linear estimators proposed in [START_REF] Lepski | Lectures at the Seminar Paris-Berlin[END_REF].

We introduce a discrete grid on the unit sphere {ϑ ∈ R d : ϑ = 1}, and we divide the domain of definition of x into small blocks. For each block, we consider a finite set of estimators ĝJ (x) extracted from the family [START_REF] Lepski | On problems of adaptive estimation in white Gaussian noise[END_REF], with x which is fixed as the center x 0 of the block and all the ϑ on the grid. We then select a data-dependent θ in the grid applying our aggregation procedure to this finite set. The value of our final estimator g * A,ε on this block is constant and is defined as g * A,ε (x) ≡ ĝ(A, θ,λε(A)) (x 0 ). We thus get a piecewiseconstant estimator g * A,ε on [-1, 1] d which depends only on A and on the observations (the exact definition of g * A,ε is given in Section 6). Remark 5. If A is unknown we need simultaneous adaptation to A and to ϑ, i.e., to the smoothness and to the local structure of the underlying function. Note, however, that parameters A and ϑ are not independent. In particular, A determines the form of the neighborhood where we have an unknown local structure depending on ϑ. This is important because our construction of the family of estimators {ĝ J , J ∈ J} heavily relates on the local representation of the model. For example, if the family {ĝ J , J ∈ J} does not contain an estimator corresponding to the correct local structure, the choice from this family cannot even guarantee consistency. Another difficulty is that different values of A can correspond to different types of local models (cf. Figure 2). Therefore, if A is totally unknown (fully adaptive estimation) then both the type of local structure and the form of the corresponding window are unknown. So, we see that adaptive estimation of composite functions is more difficult than classical adaptation to the unknown smoothness as considered, for example, in [START_REF] Lepski | Asymptotically minimax adaptive estimation I: upper bounds. Optimal adaptive estimates[END_REF][START_REF] Lepski | Asymptotically minimax adaptive estimation II: Statistical models without optimal adaptation[END_REF][START_REF] Lepski | On problems of adaptive estimation in white Gaussian noise[END_REF]. In a forthcoming paper we will show that it is possible to adapt to unknown type of the local structure (including adaptation to the local orientation ϑ) under certain known restrictions on A. We will call this partial adaptation. Partial adaptive procedure can be constructed in a similar way as discussed above. A difference is that we need to introduce a grid not only on the values of ϑ, but also on those of A, and we aggregate estimators corresponding to the product of both grids.

Upper bounds on the risk of the estimators

Define the following three domains of values of A = (γ, β) contained in (0, 2] 2 (cf. Figure 4).

P 1 = A : γ ≤ 1, 1 < β ≤ 2 , P 2 = A : 1 < γ ≤ β ≤ 2, β ≥ d(γ -1) + 1 , P 3 = A : 1 < γ ≤ β ≤ 2, β < d(γ -1) + 1 . ( 17
)
In view of the above discussion, these are exactly the zones where improved rates occur and where the local structure is active. For the sake of completeness, we consider also the remainder zone (zone of no local structure):

P 4 = (0, 1] 2 ∪ {(γ, β) : 1 ≤ β < γ ≤ 2}.
As we will see it later, the optimal kernels K J are defined differently for each of these zones. Theorem 2. Let φ ε (γ, β) be as in [START_REF] Horowitz | Rate-optimal estimation for a general class of non-parametric regression models with unknown link function[END_REF]. For any A = (γ, β) ∈ (0, 2] 2 \ P 2 and any p > 0 the estimator g * A,ε satisfies

lim sup ε→0 sup g∈H(A,L) E g φ -1 ε (γ, β) g * A,ε -g ∞ p < ∞.
For any A = (γ, β) ∈ P 2 and any p > 0 the estimator g * A,ε satisfies

lim sup ε→0 sup g∈H(A,L) E g ln ln (1/ε) -1 φ -1 ε (γ, β) g * A,ε -g ∞ p < ∞.
Combining Theorems 1 and 2 we conclude that φ ε (γ, β) is the minimax rate of convergence for the class H(A, L) if A = (γ, β) ∈ (0, 2] 2 \ P 2 , and that it is near minimax (up to the ln ln(1/ε) factor) if A = (γ, β) ∈ P 2 . Therefore, our estimator g * A,ε is respectively rate optimal or near rate optimal on H(A, L).

Theorem 2 can be viewed as a result on adaptation to the unknown local structure of the function to be estimated: the estimator g * A,ε locally adapts to the "correct" orientation ϑ which is a vector collinear to the gradient ∇G(x) in a neighborhood of x.

Extensions

Related statistical models

1. We consider here the Gaussian white noise model because its analysis requires a minimum of technicalities. Composition structures can be studied for more realistic models, such as nonparametric regression with random design, nonparametric density estimation and classification. Note that our theorems can be directly transposed to gaussian nonparametric regression model with fixed equidistant design using the equivalence of experiments argument (cf. (4; 21)). Note also that results similar to ours have been recently obtained for the problem of testing hypotheses about composite functions in the Gaussian white noise model [START_REF] Lepski | Hypothesis testing under composite functions alternative[END_REF].

2. We restrict our study to the sup-norm loss and to the Hölder smoothness classes. A natural extension would be to consider models where the risk is described by other norms and other smoothness classes. Typical candidates here are Sobolev and Besov classes, the classes of monotone or convex functions. The case of functional classes with anisotropic smoothness is of interest as well. Estimation in other norms may lead to unexpected effects (cf. Remark 4).

3. We consider only the simplest composition f (G), where f : R → R and G : R d → R.

A more general description could be f (G 1 , . . . , G k ), where f : R k → R and G s : R ds → R, s = 1, . . . , k, and

d 1 + • • • + d k = d.

4.

A related more complex modeling can be based on Kolmogorov's theorem of representation of a continuous function of several variables by compositions and sums of functions of one variable addition (13; 22).

Possible refinements

1. In this paper we treat only the case A ∈ (0, 2] 2 . Extension to A / ∈ (0, 2] 2 remains an open problem. However, our lower bound (Theorem 1) is valid for all A ∈ R 2 + . We believe that it cannot be improved. This conjecture is supported by recent results on a hypothesis testing problem with composite functions [START_REF] Lepski | Hypothesis testing under composite functions alternative[END_REF] which is closely related to our estimation problem. The upper bound proved in [START_REF] Lepski | Hypothesis testing under composite functions alternative[END_REF] for all A ∈ R 2 + in the problem of testing hypotheses coincides with our lower bound.

2. The rate of convergence of the minimax procedure (cf. Theorem 2) in the zone P 2 contains an additional ln ln(1/ε) factor, as compared to the lower bound of Theorem 1. This deterioration of the rate is due to the method of aggregation of linear estimators that we use and does not seem to be unavoidable.

3. In Remark 5 we discussed a possible construction for partial adaptation. The ultimate goal of adaptation is, however, to find an estimator which is totally parameter free. Such an estimator should achieve the minimax rate (6) simultaneously for all A ∈ R 2 + .

Definition of the estimator

We first introduce some notation. For a bounded function K ∈ L 1 (R d ) and p ≥ 1 we denote by K p its L p -norm and by K * g its convolution with a bounded function g:

K p = |K(t)| p dt 1/p , [K * g](x) = K(t -x)g(t)dt, x ∈ R d
(here and in the sequel = R d ). We denote J (A, ϑ, λ) where A = (γ, β) ∈ (0, 2] 2 , ϑ is a unit vector in R d and λ > 0. The class of all such triplets J is denoted by J.

Given a unit vector ϑ, let M ϑ ∈ R d×d stand for an orthogonal matrix with the first column equal to ϑ. The collection of the kernels we consider in the sequel is defined as

K J (x) = K (A,λ) (M T ϑ x)
where K (A,λ) : R d → R is a kernel that will be defined in Section 7. Next, for any J ′ , J , ∈ J and all t ∈ R d we define the convoluted kernel

K J ′ * J (t) = K J ′ (t -y)K J (y)dy and the difference ∆ J ′ K J ′ * J = K J ′ * J -K J ′ .
Note that, by definition, the kernel K J is symmetric, i.e., K J (t) = K J (-t), and

K J ′ * J = K J * J ′ . (18)
For all J ∈ J and all

x ∈ [-1, 1] d set ĝJ (x) = D K J (t -x)X ε (dt),
and for all J ′ , J ∈ J define the convoluted estimator

ĝJ ′ * J (x) = D K J ′ * J (t -x)X ε (dt).
In what follows we assume that ε is small enough so that ln ln(1/ε) > 0 and that in the above expressions and in all the subsequent expressions containing convolutions with the kernels we can replace D by R d (this is possible for small ε since all the kernels that we consider are compactly supported). We also define

∆ J ′ ĝJ ′ * J (x) = ĝJ ′ * J (x) -ĝJ ′ (x)
and set

TH ε (J ′ , J ) = C(p, d) ( K J ′ 1 + K J 1 ) K J ′ 2 ε ln (1/ε),
where C(p, d) = 2 + √ 4p + 8d. To define the estimator we first introduce a discrete grid on the set of indices J. We discretize only the ϑ-coordinate of J . Recall that ϑ takes values on the Euclidean unit sphere S in R d .

Discretization Let S ε ⊂ S be an ε 2 -net on S, i.e., a finite set such that

∀ϑ ∈ S ∃ϑ ′ ∈ S ε : ϑ -ϑ ′ ≤ ε 2 ,
and card(S ε ) ≤ ( √ dε -2 ) d . W.l.o.g. we will assume that (1, 0, . . . , 0) ∈ S ε . Fix A ∈ (0, 2] 2 and define λ ε (A) as a solution in λ of the bias-variance balance equation

c 11 λ = ε ln (1/ε) K (A,λ) 2 (19)
where c 11 is a constant in Proposition 1 below, depending only on A, L and d. Finally we define the following grid on the values of J :

J grid J = (A, ϑ, λ ε (A)) : ϑ ∈ S ε ⊂ J.
Acceptability For a given x ∈ [-1, 1] d we define a subset Tx of J grid as follows:

J ∈ Tx ⇐⇒ ∆ J ′ ĝJ ′ * J (x) ≤ TH ε (J ′ , J ), ∀J ′ ∈ J grid .
Any value J belonging to Tx is called acceptable. Note that the threshold TH ε (J ′ , J ) can be bounded from above and replaced in all the definitions by a value that does not depend on J , J ′ ∈ J grid . In fact, either

TH ε (J ′ , J ) ≍ λ ε (A) if A ∈ P 1 ∪ P 3 or TH ε (J ′ , J ) ≍ ln ln (1/ε)λ ε (A) if A ∈ P 2 .
Estimation at a fixed point For any x ∈ [-1, 1] d such that Tx = ∅ we select an arbitrary Ĵx from the set Tx . Note that the set Tx is finite, so a measurable choice of Ĵx is always possible; we assume that such a choice is effectively done. We then define the estimator g * * (x) as follows:

g * * (x) ĝ Ĵx (x) if Tx = ∅, 0 if Tx = ∅. ( 20
)
Global estimator The estimator g * * is defined for all x ∈ [-1, 1] d and we could consider x → g * * (x), x ∈ [-1, 1] d , as an estimator of the function g. However, the measurability of this mapping is not a straightforward issue. To skip the analysis of measurability, we use again a discretization. Introduce the following cubes in R d :

Π ε (z) = d k=1 ε 2 (z k -1), ε 2 z k , z = (z 1 , . . . , z d ) ∈ Z d .
For any x ∈ [-1, 1] d we consider z(x) ∈ Z d such that x belongs to the cube Π ε (z(x)), and a piecewise constant estimator g * * (z(x)). Our final estimator is a truncated version of g * * (z(x)):

g * A,ε (x) g * * (z(x)) if |g * * (z(x))| ≤ ln ln(1/ε), ln ln(1/ε) sign(g * * (z(x))) if |g * * (z(x))| > ln ln(1/ε). (21)
Thus, the resulting procedure g * A,ε is piecewise constant on the cubes Π

ε (z) ⊂ [-1, 1] d , z ∈ Z d .

Construction of the kernel

In this section, as well as in the Appendix, we will distinguish between the couple of "true" parameters A 0 = (γ, β) and a variable couple of parameters A = (a, b) ∈ (0, 2] 2 . This is done to state the lemmas in a form convenient to be applied in the context of adaptation to unknown (γ, β) which will be treated in our forthcoming work.

Depending on the value of A we use different constructions of K (A,λ) . Our objective is to obtain K J with suitable approximation properties for each J ∈ J. Let us summarize here the main requirements on the kernel:

1. Convolution of the kernel K (A,λ) with the "local model" of g corresponding to A should approximate g with the accuracy O(λ). Furthermore, the kernel should be localized, i.e., it should vanish outside of the window where the local structure is valid.

2. A basic characteristic of the kernel is its L 2 -norm which determines the variance of the kernel estimator. Our objective is to achieve its minimal value.

3. As we will see it later, the L 1 -norm of the kernel is also an important parameter of the proposed estimation procedure. Our objective will be to keep the L 1 -norm as small as possible.

We use different kernels K (A,λ) for A belonging to different zones P i (cf. Figure 4). The construction of K (A,λ) is trivial when A is in the zone P 4 of no local structure. In this case a basic boxcar kernel tuned to the effective smoothness of the composite function can be used. Observe that when A ∈ (0, 1] 2 the effective smoothness of the composite function equals to ab, and when A = (a, b) satisfies 1 < b ≤ a ≤ 2 the effective smoothness is b. So, we define the kernel K (A,λ) for the zone P 4 as follows:

K (A,λ) (y) =    2λ 1 ab -d I -λ 1 ab ,λ 1 ab d (y) if A = (a, b) ∈ (0, 1] 2 , 2λ 1/b -d I [-λ 1/b ,λ 1/b ] d (y) if 1 < b < a ≤ 2.
Here I A (•) stands for the indicator function of a set A. The following lemma is straightforward.

Lemma 1. For any

A 0 = (a, b) ∈ P 4 , λ > 0 and x ∈ [-1, 1] d , we have sup g∈H(A,L) |[K (A,λ) * g](x) -g(x)| ≤ c 0 λ,
where the constant c 0 depends only on L and d. Furthermore,

K (A,λ) 1 = 1 and K (A,λ) 2 =      2λ 1 ab -d/2 , (a, b) ∈ (0, 1] 2 2λ 1 b -d/2 , 1 < b < a ≤ 2.
We turn now to the analysis of cases with active local structure. We start with the zone P 1 of local single model.

Kernel for the local single index model

The zone of local single-index model is

P 1 = {A = (a, b) : a ≤ 1, 1 < b ≤ 2 . For any A ∈ P 1 and λ > 0 consider the hyperrectangle Π λ (A) = -λ 1/a , λ 1/a × -λ 1 ab , λ 1 ab d-1
and define the kernel K (A,λ) as follows:

K (A,λ) = 2 d λ 1 a + d-1 ab -1 I Π λ (A) (y), y ∈ R d . ( 22
)
Approximation property of the kernel K (A,λ) Let q : R → R and B : R d → R be functions such that, for given a ∈ (0, 1],

|q(x) -q(y)| ≤ L|x -y| a , ∀ x, y ∈ R d , sup x∈R d |B(x)| ≤ c 1
where c 1 > 0, L > 0 are constants. We denote by A(a) the set of all pairs of functions (q, B) satisfying these restrictions. Define Q(y) = q(y 1 ) + B(y) y ab , ∀ y ∈ R d .

We have the following evident result:

Lemma 2. For any A = (a, b) ∈ P 1 and λ > 0 we have

(i) sup (q,B)∈A(a) K (A,λ) * Q (0) -q(0) ≤ c 2 λ
where c 2 is a constant depending only on L, c 1 and d. Moreover,

(ii) K (A,λ) 1 = 1 and K (A,λ) 2 = 2 d λ 1 a + d-1 ab -1/2 .

Kernels for the combined local model

The zone of combined local model is

P 2 ∪ P 3 = {A = (a, b) : 1 < a ≤ b ≤ 2}.
The definition of the kernel in this case is more involved. Indeed, taking K (A,λ) as a simple product of boxcar kernels [START_REF] Sprecher | An improvement in the superposition theorem of Kolmogorov[END_REF] results for A ∈ P 2 ∪ P 3 in too large approximation error.

Our aim is to construct a smoothing kernel K (A,λ) : R d → R with the following properties:

for some c > 0, it should vanish outside the set (cf. [START_REF] Ibragimov | Bounds for the risks of nonparametric regression estimates[END_REF])

y ∈ R d : |y 1 | ≤ cλ 1 b , y ≤ cλ 1 ab , |y 1 | a-1 y b ≤ cλ .
for a function q(y 1 ) of the first component y 1 of y ∈ R d , the "characteristic size" of K (A,λ) should be λ . Namely, we want to ensure the relations K (A,λ) (y)q(y 1 )dy = (2λ

1 a ) -1 λ 1 a -λ 1 a
q(y 1 )dy 1 , and

K (A,λ) (y)Q(y 2 , ..., y d )dy = (2λ 1 b ) -(d-1) λ 1 b -λ 1 b ... λ 1 b -λ 1 b Q(y 2 , ..., y d )dy 2 ...dy d .
These properties are crucial to guarantee that the bias of kernel approximation is of the order O(λ) (cf. Lemma 3 below). Note that the simple rectangular kernel [START_REF] Sprecher | An improvement in the superposition theorem of Kolmogorov[END_REF] used for the local single index model can attain such a bias, but only at the price of too large L 2 -norm (which characterizes the variance). We now give an example showing how a kernel with the required properties can be constructed in a particular case.

The two-step kernel Set

u 1 = λ 1 a , u 2 = λ 1 b , v 1 = λ b-a+1 b 2 , v 2 = 1 2 λ 1 b , (23)
Π 1,1 = [0, u 1 ] × [v 2 , v 1 ] d-1 , µ 1,1 = u 1 (v 1 -v 2 ) d-1 ; Π 2,2 = [u 1 , u 2 ] × [0, v 2 ] d-1 , µ 2,2 = (u 2 -u 1 )v d-1 2 ; Π 2,1 = [u 1 , u 2 ] × [v 2 , v 1 ] d-1 , µ 2,1 = (u 2 -u 1 )(v 1 -v 2 ) d-1 . Next, we define, for y ∈ R d + , Λ(y) = µ -1 1,1 I Π 1,1 (y) -µ -1 2,1 I Π 2,1 (y) + µ -1 2,2 I Π 2,2 (y). ( 24
) - + + v 1 v 2 u 1 u 2 Π 1,1 Π 2,2 Π 2,1
K (A,λ) (y) = 2 -d Λ(|y|). ( 25
)
We will call this kernel the two-step kernel (cf. Figure 5). Its key property is as follows.

First, for any integrable function q(y 1 ) of the first coordinate y 1 we have

K (A,λ) (y)q(y 1 )dy = 1 2u 1 u 1 -u 1 q(y 1 )dy 1 ,
since the integral of q over Π 2,1 is exactly the same as that over Π 2,2 . Further, for any integrable function Q(y 2 , ..., y d ) of y 2 , ..., y d ,

K (A,λ) (y)Q(y 2 , ..., y d )dy = (2v 2 ) -(d-1) v 2 -v 2 ... v 2 -v 2 Q(y 2 , ..., y d )dy 2 ...dy d ,
since the integral of Q over Π 2,1 is exactly the same as that over Π 1,1 . In words, the negative term -µ -1 2,1 I Π 2,1 (y) in [START_REF] Stone | Additive regression and other nonparametric models[END_REF] allows us to compensate the excess of the bias introduced by the two other terms, so that the resulting bias remains of the order O(λ) (cf. Lemma 3 below).

For the two-step kernel [START_REF] Tsybakov | Introduction à l'estimation non-paramétrique[END_REF] we have

K (A,λ) (y)dy = 1, K (A,λ) 1 = 3, K (A,λ) 2 2 = µ -1 1,1 + µ -1 2,2 + µ -1 2,1 .
We now define ρ = (d -1)(a -1) b and consider the subset {A = (a, b) : ρ ≥ (ba)/a} of P 3 . It is easy to see that for ρ ≥ (ba)/a we have

K (A,λ) 2 2 = O(λ -d b ).
Since a ≤ b for A ∈ P 3 , this result is better than part (ii) of Lemma 2 where K (A,λ) is a rectangular kernel. But we need the condition ρ ≥ (ba)/a. It is clearly satisfied when ρ ≥ 1 (recall that a > 1, b ≤ 2). For smaller values of ρ we need to add extra "steps" in the construction, i.e., to introduce piecewise constant kernels with more and more pieces of the pavement, in order to get the bias compensation property as discussed above. For instance,

if ρ + ρ 2 ≥ b-a a (since (b -a)/a < 1, this is certainly the case when ρ ≥ √ 5- 1
2 ) we need a pavement of five sets Π i,j in order to obtain a piecewise constant kernel with the required statistical properties, and so on. We come to the following construction of the kernel.

Generic construction Define a piecewise constant kernel K (A,λ) as follows. Fix an integer r that we will further call number of steps (of kernel construction). Let (u j ) j=1,...,r and (v j ) j=1,...,r+1 be, respectively, a monotone increasing and a monotone decreasing sequence of positive numbers with

u 1 = λ 1 a , v r = λ 1 b /2 and v r+1 = 0. We set Π 1,1 = [0, u 1 ] × [v 2 , v 1 ] d-1 , µ 1,1 = u 1 (v 1 -v 2 ) d-1 .
For i = 2, . . . , r and j = i -1, i we define

Π i,j = [u i-1 , u i ] × [v j+1 , v j ] d-1 , µ i,j = (u i -u i-1 )(v j -v j+1 ) d-1 . For y ∈ R d + consider Λ 1 (y) = 1 µ 1,1 I Π 1,1 (y); Λ i (y) = 1 µ i,i I Π i,i (y) - 1 µ i,i-1 I Π i,i-1 (y), i = 2, ..., r.
The kernel K (A,λ) is defined for y = (y 1 , . . . , y d ) ∈ R d as follows:

K (A,λ) (y) = 2 -d r i=1 Λ i (|y|) ( 26
)
where |y| = (|y 1 |, . . . , |y d |). Clearly,

K (A,λ) (y)dy = 1, K (A,λ) 1 = 2r -1.

Construction of the kernel for

A ∈ P 3 = {A : 1 < a ≤ b ≤ 2, b < d(a -1) + 1} If ρ ≥ b-a
a we define K (A,λ) as a two-step kernel, i.e., we set r = 2 and take (u j ) and (v j ) as in [START_REF] Stone | Optimal global rates of convergence for nonparametric regression[END_REF].

If ρ < b-a a we use another definition. We introduce the sequence (α k ) k≥0 as follows:

α 0 = b -1 , α k+1 = α k ρ + b -1 = b -1 k+1 i=0 ρ i , k = 1, 2, . . . . (27)
The sequence (α k ) is monotone increasing and, since b < d(a -1) + 1, we have

lim k→∞ α k = ∞ if ρ ≥ 1, lim k→∞ α k = (b -(a -1)(d -1)) -1 > 1 a if ρ < 1. ( 28
)
Thus we can define an integer r ≥ 2 such that

α r-1 ≥ 1 a > α r-2 . ( 29
)
Note that r depends only on A = (a, b) and d. Now we set 26), with the sequences (u j ) and (v j ) as in (30).

u 1 = λ 1 a , u i = λ α r-i , i = 2, . . . , r; v i = λ 1 b u -a-1 b i+1 , i = 1, . . . , r -1. (30) Recall that v r = 1 2 λ 1 b and v r+1 = 0. If ρ < b-a a define the kernel K (A,λ) by (
Properties of the kernel K (A,λ) Let q : R → R and p : R d → R, B : R d → R be functions such that p is continuously differentiable and, for given A = (a, b) ∈ P 2 ∪ P 3 and λ > 0,

q(0) - 1 2λ 1/a λ 1/a -λ 1/a q(z)dz ≤ c 3 λ, (31)
p(z ′ ) -p(z) -[∇p(z)] T (z ′ -z) ≤ L z ′ -z b , ∀ z, z ′ ∈ R d , (32)
sup

x∈R d |B(x)| ≤ c 4 ( 33
)
where c 3 , c 4 and L are positive constants. Let B(A, λ) denote the set of triplets (q, p, B) satisfying ( 31) -(33). Define

Q(y) = q(y 1 ) + p(y) + B(y)|y 1 | a-1 y b , ∀ y ∈ R d .
Lemma 3. Let A = (a, b) ∈ P 3 . Let the kernel K (A,λ) be defined by ( 26), with the sequences (u i ) and (v i ) as in (30) if ρ < b-a a , and with r = 2, (u i ) and (v i ) as in [START_REF] Stone | Optimal global rates of convergence for nonparametric regression[END_REF] if ρ ≥ b-a a . Then, for any λ > 0 small enough, sup (q,p,B)∈B(A,λ)

K (A,λ) * Q (0) -Q(0) ≤ cλ, (34)
K (A,λ) (y)| y m du ≤ c ′ λ m ab , ∀m ∈ R, ( 35
)
where the constant c depends only on c 3 , c 4 , L, d and A, and c ′ depends only on m, d and A. Furthermore,

K (A,λ) 1 ≤ c ′′ and K (A,λ) 2 ≤ c (3) λ -d 2b ( 36
)
where the constants c ′′ and c (3) only depend on A and d.

Note that for ρ ≥ b-a a the kernel K (A,λ) in this lemma is just the two-step kernel. The corresponding pavement {Π i,j } only contains three sets (cf. Figure 5).

The kernel K (A,λ) depends on A = (a, b) in such a way that the constants in the bounds (34) -(36) diverge when A approaches the boundary d(a -1) + 1 = b of the zone P 3 . So, Lemma 3 cannot be extended to A ∈ P 2 .

Construction of the kernel for A ∈ P 2

We consider now another choice of the sequences (u i ) and (v i ) which provides the kernel K (A,λ) with the properties similar to those of Lemma 3 but satisfied for all A ∈ P 2 ∪ P 3 and, what is more, uniformly over this set. The price to pay for the uniformity is an extra log log(1/λ) factor in the bound for the L 1 -norm of K (A,λ) .

If (ba)/a ≤ (1 + ρ)ρ we define the kernel as in Lemma 3. If (ba)/a > (1 + ρ)ρ we use another definition of sequences (u i ) and (v i ). For any 0 < λ < 1 we define

V (λ) = ln (a -1)(b -a) ab 2 ln (1/λ) . ( 37
)
If V (λ) ≤ 0 we define K (A,λ) as a two-step kernel, i.e., we set r = 2 and take (u j ) and (v j ) as in [START_REF] Stone | Optimal global rates of convergence for nonparametric regression[END_REF]. If V (λ) > 0 we define r = r(λ) > 1 by

r = min s ∈ N : s > 1, V (λ) s -1 < 1 2 ln √ 5 + 1 2 . Next, set α = V (λ) r-1 , ν = √ 5+1 2 1/2
and define the sequences (u i ) and (v i ) as follows

u i = λ 1 a exp b a-1 exp(α(i -1)) , i = 1, ..., r, v i = λ 1 ab exp {-ν exp(αi)} , i = 1, . . . , r -1, v r = 1 2 λ 1 b . (38)
Note that u r = λ

1 b . Lemma 4. Let A = (a, b) ∈ P 2 ∪ P 3 .
Let the kernel K (A,λ) be defined by ( 26), with the sequences (u i ), (v i ) as in Lemma 3 if (1+ ρ)ρ ≥ b-a a , with the sequences (u i ), (v i ) as in (38) if (1 + ρ)ρ < b-a a and V (λ) > 0, and with r = 2, (u i ) and (v i ) as in [START_REF] Stone | Optimal global rates of convergence for nonparametric regression[END_REF] if (1 + ρ)ρ < b-a a and V (λ) ≤ 0. Then, for any λ > 0 small enough,

sup (q,p,B)∈B(A,λ) K (A,λ) * Q (0) -Q(0) ≤ c 5 λ, (39)
K (A,λ) (y) y m du ≤ c 6 λ m ab , ∀m ∈ R, (40)
where the constant c 5 depends only on c 3 , c 4 , L and d, and c 6 > 0 depends only on m and d (both constants are explicit in the proof of the lemma). Furthermore,

K (A,λ) 1 ≤ c 7 ln ln λ -1 and K (A,λ) 2 ≤ c 8 λ -b+d-1 2ab (41)
where the constants c 7 and c 8 only depend on d.

Some remarks are in order here.

1. The number of steps r in the construction of the kernel is typically small. In particular, r = 2 if ρ ≥ b-a a , and

r = 3 if (1 + ≥ b-a a > ρ (cf. ( 29
)). Moreover, for 1 < a ≤ b ≤ 2 we have (a -1)(b -1) ab 2 ≤ (b -1) 2 b 3 ≤ 1 8 . Hence, V (λ) ≤ ln √ 5+1 2
for all λ > 3 • 10 -6 which means that, for (1 + ρ)ρ < b-a a , no more than 3 steps of the construction are needed if λ > 3 • 10 -6 . In other words, unless we are not "extremely far" in the asymptotics, the number of steps r does not exceed 3 and thus the L 1 -norm of the resulting kernel K (A,λ) is bounded by 5.

2. In the asymptotics when λ → 0 the number of steps r = r(λ) in the construction and thus the L 1 -norm of the kernel K (A,λ) is at most O(ln ln λ -1 ). As discussed in the previous remark, this behavior starts "extremely far" in the asymptotics, so it has essentially a theoretical interest. In the theory, it results in an extra ln ln ε -1 factor in the upper bound for the adaptive estimation procedure, as compared to the lower bound in [START_REF] Horowitz | Rate-optimal estimation for a general class of non-parametric regression models with unknown link function[END_REF]. It can be shown that for A ∈ P 2 a kernel with the required approximation properties cannot have the L 1 -norm growing slower than ln ln λ -1 , as λ → 0. On the other hand, as we have seen in Lemma 3, for A ∈ P 3 solely, there is a choice of sequences (u j ) and (v j ) such that the L 1 -norm of the kernel is bounded by a constant independent of λ. This constant, however, depends on A = (a, b) and explodes as A approaches the boundary of P 3 .

Basic approximation results

We can now describe the approximation properties of the kernel K J which serve as a main tool in the proof of the properties of the estimator

g * A,ε (x). Let x ∈ [-1, 1] d and A 0 = (γ, β) ∈ (0, 2] 2 be fixed and let g = f • G ∈ H(A 0 , L). We define ϑ x 0    (1, 0, . . . , 0) if β > 1 and ∇G(x) = 0, or β ≤ 1, ∇G(x)/ ∇G(x) if β > 1, ∇G(x) = 0. (42)
The following statement is an immediate consequence of Lemmas 1 -4.

Corollary 1. For all A 0 = (γ, β) ∈ (0, 2] 2 , and all λ > 0 we have

sup x∈[-1,1] d sup g∈H(A 0 ,L) K J x 0 * g (x) -g(x) ≤ c 10 λ
where J x 0 = A 0 , ϑ x 0 , λ and c 10 is a constant depending only on A 0 , L and d.

In other words, the collection {K J , J ∈ J} of the kernels contains an element K J x 0 such that the quality of approximation of g(x) by the "ideal" smoother K J x 0 * g (x) is of the order O(λ). Here we use the term "ideal" because J x 0 = (A 0 , ϑ x 0 , λ) depends on the gradient ∇G(x), and thus on the unknown function g.

In what follows we also need another property of kernels K J .

Proposition 1. For all

A 0 = (γ, β) ∈ (0, 2] 2 , x ∈ [-1, 1] d , λ 0 > 0 and all J = A, ϑ, λ ∈ J such that λ 1 ab ≤ 2λ 1 γβ 0 ∧ 1 we have sup A∈(0,2] 2 sup g∈H(A 0 ,L) ∆ J K J * J x 0 * g (x) ≤ c 11 K J 1 + K J x 0 1 λ 0 (43) + K J 1 K J x 0 1 ε 2 ,
where J x 0 = A 0 , ϑ x , λ 0 , ϑ x is any element of the unit sphere S such that ϑ xϑ x 0 ≤ ε 2 and c 11 is a constant depending only on A 0 , L and d. Furthermore, for any J , J ′ ∈ J we have

∆ J ′ K J ′ * J 2 ≤ K J ′ 1 + K J 1 K J ′ 2 . ( 44
)
8 Proof of Theorem 1

For any β > 0, γ > 0 and any 0 < ε < 1 define the integers

q 1 = ⌈ ε ln(1/ε) - 2 2γβ+β+(d-1) ⌉ .
Consider the regular grid Γ q 1 on [0, 1] d-1 defined by

Γ q 1 2k 1 + 1 2q 1 , . . . , 2k d-1 + 1 2q 1 : k i ∈ {0, . . . , q 1 -1}, i = 1, . . . , d -1 .
Denote by x 1 , . . . , x m , where m = card(Γ q 1 ) = q d-1 1 , the elements of Γ q 1 numbered in an arbitrary order.

Let u : R → R + be an infinitely differentiable function such that u(0

) = 1, u(t) = u(-t) for all t ∈ R, supp u = [-1/2, 1/2] and u(t) is strictly monotone decreasing on [0, 1/2]. Set f 0 (t) = u(t), ∀ t ∈ R, and ϕ 0 (t 2 , . . . , t d ) = 1 2 d j=2 u(t j ), ∀ t ∈ R.
Define the following infinitely differentiable functions of t

= (t 1 , . . . , t d ) ∈ [-1, 1] d : g 0 (t) = L 0 h γ f 0 t 1 h , g k (t) = L 0 h γ f 0 t 1 h + L 0 h β 1 h ϕ 0 t 2 -x k,2 h 1 , . . . , t d -x k,d h 1 , k = 1, . . . , m,
where h = h β 1 , h 1 = 1/q 1 , 0 < L 0 < 1 is a constant to be chosen small enough, and x k,j stands for the jth component of x k . We note that, in view of the above definitions, the sets where the functions g l and g k differ from g 0 are disjoint for l = k, k = 0, l = 0.

It is easy to see that if L 0 is small enough, g k ∈ H(A, L), k = 0, . . . , m. In what follows, we assume that L 0 is chosen in this way. To prove Theorem 1, we follow the scheme of proving lower bounds based on reduction to the problem of distinguishing between m + 1 hypotheses (cf., e.g., [START_REF] Tsybakov | Introduction à l'estimation non-paramétrique[END_REF]). We choose the hypotheses to be determined by g 0 , . . . , g m and we apply Theorem 2.5 of ( 25), where we consider the sup-norm distance

d(g l , g k ) = g l -g k ∞ = sup t∈[-1,1] d |g l (t) -g k (t)|, l, k = 1, . . . , m.
Since the functions g l and g k differ from g 0 on disjoint sets, for any l = k, l, k = 1, . . . , m, we have

d(g l , g k ) = d(g 0 , g k ) ≥ L 0 h γ |f 0 (0) -f 0 (L 0 h β 1 ϕ 0 (0)/h)| = L 0 h γ |f 0 (0) -f 0 (L 0 (1 + o ε (1))/2)|,
where o ε (1) → 0, as ε → 0. Since L 0 > 0 and f 0 is strictly decreasing on [0, ∞), there exists a constant L * > 0 such that, for ε small enough,

d(g l , g k ) ≥ L * h γ ≍ ε ln(1/ε) 2γ 2γ+1+(d-1)/β , l = k, l, k = 0, . . . , m. (45)
Thus, assumption (i) of Theorem 2.5 in ( 25) is satisfied with s = L * h γ /2. It remains to check assumption (ii) of that theorem. The probability measures P g k are gaussian, and the Kullback-Leibler divergence between P g k and P g 0 has the form

K(P g k , P g 0 ) = ε -2 D (g 0 (t) -g k (t)) 2 dt = ε -2 L 2 0 h 2γ D f 0 t 1 h -f 0 t 1 h + w(t 2 , . . . , t d ) 2 dt
where we write for brevity

w(t 2 , . . . , t d ) L 0 h β 1 h ϕ 0 t 2 -x k,2 h 1 , . . . , t d -x k,d h 1 .
Since, for any w ∈ R,

f 0 t 1 h -f 0 t 1 h + w 2 = w 2 1 0 f ′ 0 t 1 h + uw du 2 ≤ w 2 1 0 f ′ 0 t 1 h + uw 2 du
we find

K(P g k , P g 0 ) ≤ ε -2 L 2 0 h 2γ w 2 (t 2 , . . . , t d )dt 2 . . . dt d × 1 0 f ′ 0 t 1 h + uw(t 2 , . . . , t d ) 2 dt 1 du = L 4 0 ε -2 h 2γ+1 h d-1 1 R d-1 ϕ 2 0 (v)dv R |f ′ 0 (v 1 )| 2 dv 1 ≤ c * L 4 0 ln(1/ε)
where c * > 0 is an absolute constant. Next, m = q d-1 1 , so that ln m ≍ ln(1/ε). This and the previous inequality imply that if L 0 is chosen small enough, we have

K(P g k , P g 0 ) ≤ (1/
E g ε ln (1/ε) - 2γ 2γ+1+(d-1)/β gε -g ∞ p > 0, (47)
which is valid for all β > 0, γ > 0 and all p > 0.

We now show that for the trivial cases discussed in Section 2 we can obtain better lower bounds. Consider first the case where 0 < β, γ ≤ 1. Then we use the same technique as above, but we set now q 1 = ⌈ ε ln(1/ε) -2 2γβ+d ⌉. We then introduce a regular grid Γ * q 1 on [0, 1] d defined by

Γ * q 1 2k 1 + 1 2q 1 , . . . , 2k d + 1 2q 1 : k i ∈ {0, . . . , q 1 -1}, i = 1, . . . , d
and denote by x 1 , . . . , x m , where m = card(Γ * q 1 ) = q d 1 , the elements of Γ * q 1 numbered in an arbitrary order. We set

ϕ 0 (t) d j=1 u(t j ), ∀ t ∈ R,
and we choose the functions g k in a different way:

g 0 (t) = |t| γ , g k (t) = t + L 0 h β ϕ 0 t -x k h γ , k = 1, . . . , m,
where h = 1/q 1 . With this choice, clearly,

d(g l , g k ) ≥ L γ 0 h γβ ϕ γ 0 (0) ≍ ε ln(1/ε) 2γβ 2γβ+d , l = k, l, k = 0, . . . , m. (48)
Next,

K(P g k , P g 0 ) = ε -2 D (g 0 (t) -g k (t)) 2 dt ≤ L 2γ 0 ε -2 h 2γβ+d R d ϕ 2γ 0 (v)dv = O ln(1/ε) , as ε → 0. ( 49
)
Using ( 48), (49) and Theorem 2.5 in ( 25), the proof is completed as in the previous case, so that we get the lower bound lim inf

ε→0 inf gε sup g∈H(A,L) E g ε ln (1/ε) -2γβ 2γβ+d gε -g ∞ p > 0, (50)
which is valid for all 0 < β, γ ≤ 1 and all p > 0.

Finally, the second trivial case where (47) can be improved corresponds to γ ≥ β ∨ 1. As observed in Section 2, in this case we have the inclusion H d (β, L 4 ) ⊂ H(A, L) with some constant L 4 > 0, and we can use the standard lower bound for H d (β, L 4 ) (cf. (20; 5; 3)):

lim inf ε→0 inf gε sup g∈H(A,L) E g ε ln (1/ε) -2β 2β+d gε -g ∞ p > 0. ( 51
)
Combining the bounds (47), ( 50) and (51) we obtain the result of Theorem 1.

Proof of Theorem 2

We need the following technical result. Proof is standard (see, e.g., [START_REF] Kerkyacharian | Nonlinear estimation in anisotropic multiindex denoising I[END_REF]).

To prove Theorem 2 we proceed in steps.

1 • .Reduction to the discrete norm. Fix A = (γ, β) ∈ (0, 2] 2 , and suppose that g ∈ H(A, L). Let, for brevity, ḡ * ε = g * A,ε . In view of the construction of the global estimator (cf. ( 21)) we get, for all g ∈ H(A, L),

ḡ * ε -g ∞ ≤ sup z∈Z d max x∈Πε(z)∩[-1,1] d ḡ * ε (x) -g(x) ≤ |ḡ * ε -g| ∞ + Cε 2γ(β∧1) (52)
where

|ḡ * ε -g| ∞ max z∈Zε ḡ * ε (z) -g(z) with Z ε = ε 2 Z d ∩ [-1, 1] d .
Here and in what follows we will use the same notation C for possibly different positive constants depending only on A, L and d. Since ε 2γ(β∧1) = o(φ ε (γ, β)), ε → 0, for all (γ, β) ∈ R 2 + , it is sufficient to prove Theorem 2 with the loss given by the maximum norm | • | ∞ on the finite set Z ε . Thus, w.l.o.g. in what follows we will replace

• ∞ by | • | ∞ .
2 • .Control of large deviations. To any z ∈ Z ε we assign a vector θ z ∈ S ε such that θ zθ z 0 ≤ ε 2 where θ z 0 is defined in (42). Next, we set J z 0 (A, θ z , λ ε (A)). Introduce the random event

F = ∃z ∈ Z ε : J z 0 / ∈ Tz ,
where Tz is the set of acceptable triplets J defined in Section 6. We now show that for all ε > 0 small enough sup g∈H(A,L)

P g (F) ≤ c 12 ε 2p (53)
where the constant c 12 depends only on d. Indeed, in view of the definition of the random set Tz , F ⊆ z∈Zε J ′ ∈J grid ∆ J ′ ĝJ ′ * J z 0 (z) > TH ε J ′ , J z 0 and therefore

P g (F) ≤ z∈Zε J ′ ∈J grid P g ∆ J ′ ĝJ ′ * J z 0 (z) > TH ε J ′ , J z 0 . (54)
Note that

E g ∆ J ′ ĝJ ′ * J z 0 (z) = ∆ J ′ K J ′ * J z 0 * g (z). Applying Proposition 1 with A 0 = A, J z 0 = (A, θ z , λ ε (A)) and λ = λ 0 = λ ε (A) we obtain, sup g∈H(A,L) E g ∆ J ′ ĝJ ′ * J z 0 (z) ≤ c 11 λ ε (A) K J ′ 1 + K J z 0 1 + K J ′ 1 K J z 0 1 ε 2 (55)
Now, due to the construction of the kernel K (A,λ) and the fact that K J 1 = K (A,λε(A)) 1 for all J ∈ J grid , there exists a constant c 13 depending only on A and d such that K * A max J ∈J grid K J 1 satisfies

K * A ≤ c 13 if A ∈ (0, 2] 2 \ P 2 , K * A ≤ c 13 ln ln(1/ε) if A ∈ P 2 .
Since also K J 1 ≥ 1 and λ ε (A)/(ε ln ln(1/ε)) → ∞, as ε → 0, we have, for ε > 0 small enough, sup g∈H(A,L)

E g ∆ J ′ ĝJ ′ * J z 0 (z) ≤ 2c 11 λ ε (A) K J ′ 1 + K J z 0 1 = 2ε ln(1/ε) K (A,λε(A)) 2 K J ′ 1 + K J z 0 1 (56)
where we used that λ ε (A) is a solution of [START_REF] Nemirovskii | On nonparametric estimation of smooth regression functions[END_REF]. Note also that in P g -probability

∆ J ′ ĝJ ′ * J z 0 (z) -E g ∆ J ′ ĝJ ′ * J z 0 (z) ∼ N 0, ε 2 ∆ J ′ K J ′ * J z 0 2 2 . (57)
Using ( 44), ( 54) -( 57) and the definition of the threshold TH ε (•, •) we obtain that, for ε > 0 small enough,

P g (F) ≤ card(Z ε )card(S ε )P |ξ| > (4p + 8d) ln(1/ε) ≤ card(Z ε )card(S ε )ε 2p+4d
where ξ ∼ N (0, 1). This proves (53) since card(Z ε ) ≤ (2ε

E g φ -p ε (γ, β)|ḡ * ε -g| p ∞ I{F} = 0. ( 58
)
We now control the bias of ĝJ z 0 via Corollary 1, its stochastic error via the bounds on K (A,λε(A)) 2 in Lemmas 2 -4 and apply [START_REF] Nemirovskii | On nonparametric estimation of smooth regression functions[END_REF] to get that, for all A = (γ, β)

∈ (0, 2] 2 , lim sup ε→0 sup g∈H(A,L) E g φ -p ε (γ, β)|ĝ J z 0 -g| p ∞ < ∞. (59) 4
• .Final argument. Note that on the event F c the set Tz of acceptable triplets J is non-empty for every z ∈ Z ε , so that Ĵz exists. Thus, on F c we can write, for all z ∈ Z ε ,

ĝ Ĵz (z) -g(z) ≤ ∆ Ĵz ĝ Ĵz * J z 0 (z) + ∆ J z 0 ĝJ z 0 * Ĵz (z) + ĝJ z 0 (z) -g(z) . (60)
Further, on F c the triplet J z 0 is acceptable for all z ∈ Z ε . This and the acceptability (by definition) of Ĵz imply that on F c , for all z ∈ Z ε ,

∆ J z 0 ĝJ z 0 * Ĵz (z) ≤ TH ε J z 0 , Ĵz , ∆ Ĵz ĝ Ĵz * J z 0 (z) ≤ TH ε Ĵz , J z 0 . (61)
This, the definition of the threshold TH ε and the fact that

K J 2 = K (A,λε(A)) 2 for all J ∈ J grid yield that on F c , for all z ∈ Z ε , ĝ Ĵz (z) -g(z) ≤ 4C(p, d)K * A K (A,λε(A)) 2 ε ln(1/ε) + ĝJ z 0 (z) -g(z) = 4C(p, d)c -1 11 K * A λ ε (A) + ĝJ z 0 (z) -g(z) . (62)
We combine (59) and (62) to get, with some constants c 14c 16 independent of ε,

sup g∈H(A,L) E g |ḡ * ε -g| p ∞ I{F c } ≤ c 14 (K * A λ ε (A)) p + c 15 φ p ε (γ, β) ≤ c 16 (K * A φ ε (γ, β)) p . (63)
Theorem 2 follows now from (58) and (63).

A Proofs of auxiliary results

A.1 Proof of Proposition 1

1 • . Preliminary remarks.
For any J ∈ J and any x ∈ [-1, 1] d we may write

∆ J K J * J x 0 * g (x) = K J * J x 0 * g (x) -K J * g (x) (64) = K J (y -x)K J x 0 (t -y)dy g(t)dt -K J * g (x) = K J (y -x) K J x 0 (t -y)g(t)dt dy -K J * g (x) = K J (y -x)g(y)dy -K J * g (x) + K J (y -x) K J x 0 (t -y) g(t) -g(y) dt dy = K J (y -x) K J x 0 (t -y) g(t) -g(y) dt dy = K J (v) K J x 0 (z) g(z + v + x) -g(v + x) dz dv = K (A,λ) (M T ϑ v) K (A 0 ,λ 0 ) (M T ϑ x z)(g(z + v + x) -g(v + x))dzdv. Define G x (•) = G(• + x) and f x (•) = f (• + G(x)). Then g(z + v + x) = f (G x (z + v)) and g(v + x) = f (G x (v)). Note that, for all x ∈ [-1, 1] d , G x ∈ H d (β, L 2 ), f x ∈ H 1 (γ, L 1 ). ( 65
)
If 1 < γ ≤ 2, the second property in (65) implies

f ′ x ∈ H 1 (γ -1, 2L 1 ). (66)
In the case where 1

< β ≤ 2, for all u ∈ R d , x ∈ [-1, 1] d we define Gx (u) = G x (u)-G x (0)- [∇G x (0)] T u. In view of (65), for all x ∈ [-1, 1] d we have ∇ Gx (u) ≤ 2L 2 , ∀ u ∈ R d , (67)
Gx (t) -Gx (u) -[∇ Gx (u)] T (t -u) ≤ L 2 t -u β , ∀ t, u ∈ R d , (68)
⇒ | Gx (u)| ≤ L 2 u β , u ∈ R d .
It follows from the definition of K (A,λ) and Lemmas 1 -4 that

v γβ K (A,λ) (v) dv ≤ c ′ 6 λ γβ ab , ∀ A ∈ (0, 2] 2 , λ > 0, (69)
where c ′ 6 > 0 is a constant depending only on L and d. Furthermore, for any A = (a, b) ∈ (0, 2] 2 and any λ ≤ 1 the support of

K (A,λ) is contained in a ball {u ∈ R d : u ≤ c K λ 1 ab }
where the constant c K > 0 depends only on d. Therefore,

K (A,λ) M T ϑ u = 0, ∀ u, ϑ ∈ R d : u > c K λ 1 ab , ϑ = 1. ( 70
) 2 • . Proof for the zone of combined local model: 1 < γ ≤ β ≤ 2.
Using (65) and the Taylor expansion for G x we obtain, for all

x ∈ [-1, 1] d , z, v ∈ R d , g(z + v + x) = f G x (0) + ∇G x (0) T (z + v) + Gx (z + v) (71) = f x ∇G x (0) T (z + v) + Gx (z + v) .
Note that, by definition, ∇G x (0) = ∇G(x) = ϑ x 0 ∇G(x) . Set ∇G * = ϑ x ∇G(x) and define g

* (z + v + x) = f x ∇G * T (z + v) + Gx (z + v) .
We now approximate g(z + v + x) by g * (z + v + x) in the last line of (64). In view of (70), it suffices to consider there only the values z, v satisfying z , v ≤ c K . For such z, v and all

x ∈ [-1, 1] d , the condition ϑ x 0 -ϑ x ≤ ε 2 and (65) imply |g(z + v + x) -g * (z + v + x)| ≤ 2c K L 1 ∇G(x) ε 2 ≤ 2c K L 1 L 2 ε 2 . ( 72
)
Using ( 65) -( 68), the Taylor expansion for f x and (66), we get that for all x ∈ [-1, 1] d , z, v ∈ R d the following representation holds:

g * (z + v + x) = f x ∇G * T (z + v) + f ′ x ∇G * T (z + v) Gx (z + v) + B x,1 (z, v) z + v γβ = f x ∇G * T (z + v) + f ′ x ∇G * T (z + v) -f ′ x ∇G * T v × Gx (v) + ∇ Gx (v) T z + f ′ x ∇G * T v Gx (z + v) -Gx (v) + f ′ x ∇G * T v Gx (v) + B x,2 (z, v) ∇G * T z γ-1 z β + B x,1 (z, v) z + v γβ , (73)
where, for all

x ∈ [-1, 1] d , z, v ∈ R d , B x,1 (•, •) and B x,2 (•, •) are functions satisfying B x,1 (z, v) ≤ L 1 L γ 2 , B x,2 (z, v) ≤ 2L 1 L 2 . (74)
Putting z = 0 in (73) we obtain

g * (v + x) = f x ∇G * T v + f ′ x ∇G * T v Gx (v) + B x,1 (0, v) v γβ . (75)
From ( 73) and (75) we get, for all

x ∈ [-1, 1] d , z, v ∈ R d , g * (z + v + x) -g * (v + x) = f x ∇G * T (z + v) -f x ∇G * T v + f ′ x ∇G * T (z + v) -f ′ x ∇G * T v Gx (v) + ∇ Gx (v) T z + f ′ x ∇G * T v Gx (z + v) -Gx (v) + B x,2 (z, v) ∇G * T z γ-1 z β + B x,1 (z, v) z + v γβ -B x,1 (0, v) v γβ . ( 76
) Put u = M T ϑ x v, s = M T ϑ x z. We get from (76) that g * (M ϑ x s + M ϑ x u + x) -g * (M ϑ x u + x) = fx (s 1 + u 1 ) -fx (u 1 ) + A u,x (s 1 ) G x (u) + ∇G x (u) T s + f ′ x ∇G(x) u 1 G x (s + u) -G x (u) + Bx,2 (s, u)|s 1 | γ-1 s β + Bx,1 (s, u) s + u γβ -Bx,1 (0, u) u γβ , (77)
where s 1 and u 1 are the first components of s ∈ R d and u ∈ R d respectively, fx (u 1 ) = f x ( ∇G(x) u 1 ), G x (u) = Gx (M ϑ x u), Bx,1 (s, u) = B x,1 (M ϑ x s, M ϑ x u) Bx,2 (s, u) = ∇G(x) γ-1 B x,2 (M ϑ x s, M ϑ x u),

and

A u,x (s 1 ) = f ′ x ∇G(x) (s 1 + u 1 )f ′ x ∇G(x) u 1 . It is easy to see that inequalities (67) and (68) remain valid with G x in place of Gx . Now for all x ∈ [-1, 1] d , s, u ∈ R d we introduce q u,x (s 1 ) = fx (s 1 + u 1 ) -fx (u 1 ) + A u,x (s 1 ) G x (u) + ∇G x (u) T ϑ x s 1 +f ′ x ∇G(x) u 1 ∇G x (u) T ϑ x s 1 , p u,x (s) = f ′ x ∇G(x) u 1 G x (s + u) -G x (u) -∇G x (u) T s , B u,x (s) = Bx,2 (s, u), Q u,x (s) = q u,x (s 1 ) + p u,x (s) + Bx,2 (s, u)|s 1 | γ-1 s β , P u,x (s) = f ′ x ∇G(x) (s 1 + u 1 ) ∇G x (u) T s ⊥

where s ⊥ = ss 1 ϑ x . With this notation (77) can be written as g * (M ϑ x s + M ϑ x u + x)g * (M ϑ x u + x) = Q u,x (s) + P u,x (s) (78) + Bx,1 (s, u) s + u γβ -Bx,1 (0, u) u γβ .

We now prove that, for all x ∈ [-1, 1] d and all u ∈ R d such that u ≤ c K λ 1 ab (cf. (70)), where λ 1 ab ≤ 2λ 1 γβ 0 , the triplet (q u,x , p u,x , B u,x ) belongs to the set B(A 0 , λ 0 ) (cf. definition before Lemma 3), and thus Lemmas 3 or 4 can be applied. We need to check (31) -(33).

Checking (31). In view of (65) we have

| fx (s 1 + u 1 ) -fx (u 1 ) -f ′ x (u 1 )s 1 | ≤ L 1 L 2 |s 1 | γ .
|s 1 | γ ds 1 ≤ L 1 L 2 2 λ 0 . (79)
Next, remark that (66) implies |A u,x (s 1 )| ≤ 2L 1 L γ-1 2 |s 1 | γ-1 . Furthermore, (68) with G x in place of Gx yields |G x (u)| ≤ L 2 u β . Now, q u,x (0) = 0 and using these remarks, (79) and Checking (32) and (33). It suffices to note that, for all x ∈ [-1, 1] d , the first property in (68) with G x in place of Gx and the second property in (65) yield p u,x (s ′ )p u,x (s) -[∇p u,x (s)] T (s ′s) ≤ f ′

x ∇G(x)

u 1 L 2 s ′ -s β ≤ L 1 L 2 s ′ -s β , ∀ s, s ′ ∈ R d .
This proves (32) with b = β and L = L 1 L 2 . Finally, (33) with B = B u,x , c 4 = 2L 1 L γ 2 follows from (74).

We are now in a position to apply Lemmas 3 and 4. We demonstrate this, for example, for Lemma 4. Take there q = q u,x , p = p u,x , B = B u,x for any u ≤ c K λ 1 ab and x ∈ [-1, 1] d . Since Q u,x (0) = 0, the result (39) of Lemma 4 yields K (A 0 ,λ 0 ) (s)Q u,x (s)ds ≤ c 5 λ 0 (81)

where c 5 depends only on L and d. Furthermore, by construction the kernel K (A 0 ,λ 0 ) is symmetric, i.e., K (A 0 ,λ 0 ) (s) = K (A 0 ,λ 0 ) (-s) and hence We finally get (43) from this inequality invoking (69), (64), (72) and the condition λ 1 ab ≤ 2λ 1 γβ 0 and recalling that K (A,λ) 1 = K J 1 for all A ∈ (0, 2] 2 , λ > 0, and K (A 0 ,λ 0 ) 1 = K J x 0 1 . 31 3 • . Proof of (43) for the local single index zone: γ ≤ 1, 1 < β ≤ 2. Using (68) and the second property in (65), for all z, v ∈ R d , x ∈ [-1, 1] d we may write g * (z + v + x) = f x ∇G * T (z + v) +B x,1 (z, v) z + v γβ , where B x,1 satisfies (74). This can be viewed as a simplified version of (73). Following almost the same argument as in 2 • (the main difference is that now we drop all the terms containing f ′ x and B x,2 ) and applying Lemma 2 we obtain (43). 4 • . Proof of (43) for the zone of slow rate: (γ, β) ∈ (0, 1] 2 . Using the Hölder condition on f and G x we obtain, for all z, v ∈ R d , x ∈ [-1, 1] d ,

K (A 0 ,λ 0 ) (s)P u,x (
g(z + v + x) ≡ f G x (z + v) = f G x (0) +B x,1 (z, v) z + v γβ
where B x,1 satisfies (74). Now, (43) easily follows from this relation, (64), (69), the definition of K (A 0 ,λ 0 ) for the zone of slow rate and the condition λ

f G x (z + v) = f G x (v) + f ′ G x (v) G x (z + v) -G x (v) + B x,1 (z, v) z γ = f G x (v) + f ′ G x (v) G x (z + v) -G x (v) -∇G x (v) T z +f ′ G x (v) [∇G x (v) T z + B x,1 (z, v) z γ = f G x (v) + f ′ G x (v) [∇G x (v) T z + B x,2 (z, v) z β + B x,1 (z, v) z γ
where B x,1 satisfies (74) and |B x,2 (•, •)| ≤ L 1 L 2 . Since the kernel K (A 0 ,λ 0 ) is symmetric,

K (A 0 ,λ 0 ) (M T ϑ x z)f ′ G x (v) [∇G x (v)
T z dz = 0. Now, (43) easily follows from these relations, (64), the definition of K (A 0 ,λ 0 ) for the zone of inactive structure and the condition λ ≤ 1.

6 • . Proof of (44). For a function K ∈ L 2 (R d ), let us denote by K its Fourier transform. Using Parceval's identity we obtain, for any J , J ′ ∈ J,

∆ J ′ K J ′ * J 2 = 1 √ 2π ∆ J ′ K J ′ * J 2 = 1 √ 2π K J -1 K J ′ 2 ≤ 1 √ 2π K J ∞ + 1 K J ′ 2 ≤ K J 1 + 1 K J ′ 2 .
Since K J ′ = 1, this proves (44).

A.2 Proof of Lemma 3

First, note that some cases are trivial because the number r of steps of the kernel construction is bounded by 3. In fact, if (ρ + 1)ρ < (ba)/a and V (λ) ≤ ln

  where γ, L 1 , β, L 2 are positive constants. Here and in what follows H 1 (γ, L 1 ) and H d (β, L 2 ) are the Hölder class on R and the isotropic Hölder class on R d respectively (see Definition 1 below). The class of composite functions g with such f and G will be denoted by H(A, L), where A

Remark 4 .

 4 Theorem 1 obviously implies that the lower bound ε ln (1/ε) 2γ 2γ+1+(d-1)/β

Figure 1 .

 1 Figure 1. Zones of improved rate (cyan), of slow rate (grey) and of inactive structure (white). Dashed lines demarcate the zones of three different expressions for the exponent π(A).

Figure 2 .

 2 Figure 2. Types of local structures. 1• . Local single index model: γ ≤ 1, 1 < β ≤ 2.In this domain of γ, β, using the smoothness properties of functions f and G it is not hard to show that in the ball B λ,x (A) = t ∈ R d : tx ≤ λ

Figure 3 .

 3 Figure 3. Window for the combined local model, d = 2.

4 Figure 4 .

 44 Figure 4. Classification of zones within (0, 2] 2 .

1 a

 1 ; for a function Q(y 2 , ..., y d ) of the remaining components y 2 , ..., y d it should be λ 1 b

y 2 y 1 Figure 5 .

 15 Figure 5. Pavement Π i,j for the two-step kernel, d = 2. The kernel vanishes in the blanc zones.

Lemma 5 .

 5 Let ζ = (ζ 1 , . . . , ζ M ) be a gaussian random vector defined on a probability space (Ω, F, P) and such thatEζ m = 0, Eζ 2 m = σ 2 m , m = 1, . . . , M.Let m be a random variable with the values in (1, . . . , M) on the same probability space. Then for all A > 1 and all s > 0 we haveE (|ζ m | s ) ≤ 2A ln (M) s E σ s m + c 12 (A, s)M 1-A max m=1,...,M σ s mwhere c 12 (A, s) > 0 is a constant depending only on A and s.

fx (s 1 + u 1 )L 2 2λ

 112 fx (u 1 ) ds 1 ≤ L 1

( 1 ≤ L 1 λ 0 ≤ c 3 λ 0 where the constant c 3

 1103 67) we get, for u ≤ c K λ |A u,x (s 1 )| |G x (u)| + ∇G x (u) |s 1 | ds depends only on L and d. It can be taken as a maximum of the last expression in square brackets over (γ, β) ∈ [1, 2] 2 .

2 s 2 K

 22 s)ds = 0. (82) Next, using (74) we findBx,1 (s, u) s + u γβ -Bx,1 (0, u) u γβ ≤ 2 γβ L 1 L γ γβ + u γβ .Combining this inequality and (81) -(82) with (78) we get, for allx ∈ [-1, 1] d , u ∈ R d , K (A 0 ,λ 0 ) (s)(g * (M ϑ x s + M ϑ x u + x)g * (M ϑ x u + x))ds ≤ c 5 λ 0 + 2 γβ L 1 L γ (A 0 ,λ 0 ) (s) s γβ ds + K (A 0 ,λ 0 ) 1 u γβ .

5 •

 5 . Proof of (43) for the zone of inactive structure:1 < β ≤ γ ≤ 2. Since f ∈ H 1 (γ, L 1 ) and ∇G x (•) ≤ L 2 , for all z, v ∈ R d , x ∈ [-1, 1] d we may write

√ 5+1 2

 2 we have r ≤ 3 by definition. If (ρ + 1)ρ ≥ (ba)/a we use the kernel as in Lemma 3. But for this kernel the condition (ρ + 1)ρ ≥ (ba)/a implies that, again, r ≤ 3.So, we will treat only the remaining case where (ρ + 1)ρ < (ba)/a and V (λ) > ln√ 5+1 2

  -2 +1) d and card(S ε ) ≤ ( √ dε -2 ) d . 3 • .Two intermediate bounds on the risks. Using that |ḡ * ε | ≤ ln ln(1/ε) and g ∈ H(A, L) is uniformly bounded we deduce from (53) that, for all A = (γ, β) ∈ (0, 2] 2 ,

	lim sup	sup
	ε→0	g∈H(A,L)

Note that, by definition, α < 1 2 ln √ 5+1 2

. Further, for r ≥ 3 we have also the lower bound: α ≥ 1 4 ln

. Thus for r ≥ 3, 0.786 ≤ √ 5 + 1 2

1 • . Proof of (39). From the definition of K (A,λ) we find

Λ i (|y|)q(y 1 )dy = 2 -d Λ 1 (|y|)q(y 1 )dy = 1 u 1 q(y 1 ) + q(-y 1 ) 2 I [0,u 1 ] (y 1 )dy 1

where u 1 = λ 1/a . This and (31) imply K (A,λ) * q (0)q(0) = 2λ 1/a -1

-λ 1/a q(y 1 )dy 1q(0) ≤ c 3 λ.

We now obtain a similar bound for K (A,λ) * p (0)p(0) . Note that, in view of (32), for all z = (z 1 , . . . , z d ) ∈ R d we have

where p(z) = p(0, z 2 , . . . , z d ) and sup z∈R d |B 1 (z)| ≤ L. For the same reason, for all z (d-1) (0, z 2 , . . . , z d ) we have

where as previously |B 2 (•)| ≤ L. Combining (85) and (86) and taking into account that the function K (A,λ) is symmetric, K (A,λ) = 1 and p(0) = p(0) we get

Using this remark and the fact that b a-1 > 1 we find

and therefore u i /u r ≤ e α(i-r) . This and the equality u r = λ 1 b allow us to get

From (87), ( 88) and (90) we get

We now estimate the value

Using (92), we get similarly to (90):

where the last inequality holds for 0 < λ ≤ 1 and we used that b exp(αl) ≥ αl, ν > 1.

Summing up the results of (84), (91), (93) and taking into account (83) we obtain (39). 2 • . Proof of (40). In the same way as above we get, for 0 < λ ≤ 1,

Here and in what follows use the same notation C(d) for possibly different positive constants depending only on d.

By the definition of v r this implies that v r-1v r ≥ λ We deduce that µ r,r-

Note that by (89),

Also, as ν > 1, it is straightforward to check that

Thus, we get

Recall that we are considering the case where ρ(1

2 . This and the choice of parameters α, ν combined with (83) implies

Hence, for i = 2, . . . , r -1 we have

Note that

We deduce from (94) -(97) that

This proves the second inequality in (41). The first inequality becomes obvious if we note that V (λ) ≤ ln ln(1/λ) and so K (A,λ) 1 = 2r -1 ≤ c 7 ln ln(1/λ), for λ small enough, where c 7 is an absolute constant.

A.3 Proof of Lemma 3

Following the same lines as in the proof of (39) in Lemma 4 we obtain the bound (34) of Lemma 3 with c 5 = C(d)(c 3 + Lr + c 4 r). . Using these remarks and acting as in the proof of (40) in Lemma 4 we obtain, for 0 < λ ≤ 1,

2 • . Proof of (36). Observe that α j+1α j > 0 for j = 1, ..., r -1, so that for λ → 0 we have u j /u j-1 → ∞ and v j-1 /v j → ∞. In particular,

for all λ small enough. Next note that, by definition,

Then u 2 ≤ λ b-a abρ and for λ small enough we get by the definition of ρ:


Abdullah A Saed
email: abdullah.saed@mymail.unisa.edu.au

Siu-Wai Ho
email: siuwai.ho@unisa.edu.au

Jean-Marie Gorce
email: jean-marie.gorce@inria.fr3

Chung Shue Chen
email: chung_shue.chen@nokia-bell-labs.com

Minimal Noise Variance Decoder for Uncoordinated Multiple Access in VLC

In a visible light communications system (VLC), light sources are responsible for both illumination, communications and positioning. These light sources inevitably interfere each others at the receiver. To retain the appealing advantage that VLC systems can reuse existing lighting infrastructure, using an extra network to control or synchronize the light sources should be avoided. This paper proposes an uncoordinated multiple access scheme for VLC systems with positioning capability. The proposed scheme does not require a central unit to coordinate the transmission of the transmitters. Transmitters can be asynchronous with one another and with the receiver. Each transmitter is allocated a unique codeword with L chips for a system with up to L-1 2 transmitters where L is prime. Due to the linear growth in complexity with respect to number of transmitters, our proposed scheme is feasible for systems with large numbers of transmitters. Our novel decoder can minimize the effect of additive Gaussian noise at the receiver side. Simulation results show that the proposed decoder outperforms zero-forcing decoder.

I. INTRODUCTION

Visible light communications (VLC) has recently gained a lot attention as a promising future technology to overcome the scarcity of the radio frequency spectrum. Together with having several appealing features, VLC is usable for illumination, communications and positioning. It uses visible light for transmission generated by light emitting diodes (LEDs) as transmitters and detected by photodiodes (PDs) as receivers.

For broadcasting information using VLC, many different schemes are proposed in literature such as expurgated pulse position modulation [START_REF] Noshad | Application of Expurgated PPM to Indoor Visible Light Communications-Part II: Access Networks[END_REF], adaptive modulation based on power control [START_REF] Saed | Adaptive Modulation for Two users in VLC[END_REF] and optical code division multiple access (OCDMA) [START_REF] Shoreh | Design concepts and performance analysis of multicarrier CDMA for indoor visible light communications[END_REF]. All these techniques require a central unit to coordinate the transmission of the LEDs.

Visible light positioning systems can achieve accurate positioning. [START_REF] Hou | Multiple Access Scheme Based on Block Encoding Time Division Multiplexing in an Indoor Positioning System Using Visible Light[END_REF][START_REF] Yasir | Indoor positioning system using visible light and accelerometer[END_REF] present two positioning systems based on the received light intensity, where the receiver needs to distinguish the light intensity received from each LED. To achieve that in [START_REF] Hou | Multiple Access Scheme Based on Block Encoding Time Division Multiplexing in an Indoor Positioning System Using Visible Light[END_REF][START_REF] Yasir | Indoor positioning system using visible light and accelerometer[END_REF], time division multiple access (TDMA) was used, which also requires a central unit to coordinate the transmission of the LEDs.

To reduce system complexity and cost, it is necessary to develop uncoordinated multiple access (UMA) schemes that work without the need for a central unit nor the extra infrastructure used for controlling the LEDs. Such schemes can simply be installed by replacing the lights by LED transmitters. Thus, the installation cost is reduced as well. However, implementing UMA schemes needs to solve the following challenges. 1) Due to the possibility of using the receiver at arbitrarily locations with arbitrarily orientations, the transmitted signal suffers random channel gain that is unknown to the receiver. 2) LEDs may be controlled by different switches so that they may be switched on at different times. As a results, transmitters may be asynchronous with one another and with the receiver. 3) There is no communication link among the transmitters.

Recently, some UMA schemes usable only for positioning are proposed [START_REF] Zhang | Asynchronous indoor positioning system based on visible light communications[END_REF][START_REF] Saed | Coding for Uncoordinated Multiple Access in Visible Light Positioning Systems[END_REF]. In [START_REF] Saed | Coding for Uncoordinated Multiple Access in Visible Light Positioning Systems[END_REF], a unique sinusoidal-based code of length 2N is allocated to each transmitter for a system with N transmitters. The receiver eliminates multiple access interference (MAI) entirely by applying FFT of length 2N . Moreover, its scheme delay, i.e., the minimum received signal duration needed by the receiver to estimate the average received powers from the transmitters, grows linearly with number of transmitters. However, it does not support communications. An UMA scheme usable for communication and positioning is proposed in [START_REF] Yasir | Uncoordinated Multiple Access Scheme for VLC Systems with Positioning Capability[END_REF][START_REF] Ho | Uncoordinated Multiple Access Schemes for Visible Light Communiations and Positioning[END_REF]. In this scheme, each transmitter is assigned a unique code, however the length of the code grows exponentially with the number of transmitters. As a result, the scheme delay grows exponentially with the number of transmitters. In [START_REF] González | Cyclic code-shift extension keying for multi-user optical wireless communications[END_REF], asynchronous OCDMA is proposed for communications. To obtain a probability of error below 10 -3 for a system with 20 transmitters, the code length L ≥ 150 chips is needed so that the system throughput ≤ 0.133 bits/chip. Moreover, if k bits are wanted to be transmitted in each symbol, then 2 k different codes are allocated to each transmitter that increases the system complexity.

In this paper, we modify the sinusoidal-based codes proposed in [START_REF] Saed | Coding for Uncoordinated Multiple Access in Visible Light Positioning Systems[END_REF] and generalize them to further support information broadcast. We also propose a new decoder for our uncoordinated multiple access scheme with sinusoidal-based codes (UMA-SC) and compare it with zero-forcing decoder (ZFD). The proposed decoder not only eliminates MAI totally but also minimizes the noise variance. Our proposed scheme supports communications and positioning and does not need a central unit. Each transmitter is assigned a unique code with transmitters where L is prime. Transmitters can transmit more than one bit per symbol without the need for extra codes, which reduces the system complexity.

This paper is organised as follows. Section II presents the overview of the proposed system. The proposed scheme is detailed in Section III. Simulation results are given in Section IV, and the paper is concluded in Section V.

Notation: S(x, µ) means the cyclic shift of the vector x to the left by µ steps and S(x, µ, m) means the subset of the vector S(x, µ) from the index 0 to m -1. For example, if x = (x 0 , x 1 , x 2 ), then S(x, 1) = (x 1 , x 2 , x 0 ) and S(x, 1, 2) = (x 1 , x 2 ). All the estimates of variables denoted by the symbols of the original variables together with a hat. For example, the estimate of P i is denoted by Pi .

II. SYSTEM MODEL AND PROBLEM FORMULATION A. System Model

In our proposed system, each transmitter has a single LED. We assume that N LEDs are installed in a space. The receiver is a mobile device equipped with a PD, and it may appear at random location. The receiver's position can be estimated by any positioning algorithm based on the received light intensity. Furthermore, we assume that each LED broadcasts different data information, and this information includes the location information of each LED that is required by the positioning algorithm. The message can also depend on the system application. We assume that there is no central unit to coordinate the transmission of the LEDs, and the N LEDs may be switched on and start transmission at different times. All LEDs transmit simultaneously and they interfere the others.

B. Channel Model

Suppose the distance between LED i in a transmitter and a PD in a receiver is d i for 1 ≤ i ≤ N . If φ i is the irradiance angle with respect to the LED i's normal and ψ is the incidence angle with respect to the PD's normal, then the channel gain of a LoS optical wireless channel [START_REF] Komine | Fundamental Analysis for Visible-Light Communication System using LED lights[END_REF] between Transmitter i and a PD is given by

h i = (l + 1)A 2πd 2 i cos l (φ i )T (ψ)g(ψ) cos M (ψ), (1)
where the parameters are explained as follows. The Lambertian parameters of the LED and PD are given by l =

-log 2 log(cos(φ 1/2 )) and M = -log 2 log(cos(ψ 1/2 ))
, where φ 1/2 is the halfpower angle of irradiance of an LED and ψ 1/2 is the halfpower angle of incidence of a PD. The effective area of the PD at the receiver is given by A. The filter gain and concentrator gain are represented by T (ψ) and g(ψ), respectively. When the average transmitting optical flux of LED i is Φ i (in lumens), the average received optical power 1 of the PD is given by

P i = Φ i h i (in lux• m 2
). In our system, we do not use filter and concentrator at the receiver (i.e., T (ψ) = g(ψ) = 1). 1 Here, we use the photometric unit lux•m2 for power. Physical units lux•m 2 and Watts are interchangeable and the constant for conversion depends on the device. In this setting, we assume without loss of generality that the optical power incident on the PD from LED i is given by

P i = Φ i (l + 1)A 2πd 2 i cos l (φ i ) cos M (ψ). (2)
An UMA scheme using visible light is depicted in Fig. 1. Assume that the receiver starts receiving signals from N LEDs at time t = 0. Define b i [k] as the k-th message sent by LED i and x i (t) as the transmitted signal from LED i at time t. Suppose LED i started to transmit at t = -τ i . Since LEDs may begin transmission at different times, τ i may not be the same in general 2 . The superposition of signals from all LEDs is received together with the background light intensity Φ 0 at the receiver. Here, we assume that h i and Φ 0 are constant over a short period of time. This is justified by that it is common to achieve transmission rate over 10 6 symbols per second in VLC [START_REF] Rajagopal | IEEE 802.15.7 visible light communication: modulation schemes and dimming support[END_REF]. If the receiver's displacement and the change in the background light intensity are negligible within 10 -3 seconds, h i and Φ 0 can be seen as invariant for more than 10 3 symbols. Without loss of generality, assume that R p , the responsivity of the PD, is equal to 1. So the received signal is modelled as

y(t) = N i=1 x i (t)h i + Φ 0 + w(t), (3)
where w(t) is the thermal noise which is a real value additive white Gaussian noise (AWGN) with zero mean and variance σ 2 w . At the receiver, the signal y(t) is sampled at a finite rate 1 T .

C. Problem Formulation

Since the LEDs transmit simultaneously, y(t) is the sum of the received powers from different LEDs. Furthermore, τ i is unknown to the receiver as the system is asynchronous.

The receiver also has no information about the channel gains h i because the receiver may appear at random location with arbitrary orientation. Therefore, we need to cleverly design x i (t) such that for all i, the receiver 1) can decode the broadcast information b i [n] and 2) can estimate the average received power P i = Φ i h i for LED i.

Then positioning algorithms such as [START_REF] Yasir | Indoor positioning system using visible light and accelerometer[END_REF] can be applied.

III. PROPOSED MULTIPLE ACCESS SCHEME

In this section, the transmitter and receiver designs in our proposed system will be described. Our proposed transmitter design in Section III-A generalizes the design in [START_REF] Saed | Coding for Uncoordinated Multiple Access in Visible Light Positioning Systems[END_REF] where only positioning but no data communications is supported. A completely different receiver design will be shown in Section III-B.

A. Transmitter Design

Definition 1. Consider a prime number L with L > 2N . For i = 1, 2, . . . , N , a codeword c i with L chips is assigned to LED i where the j-th chip of c i is defined as M +1 so that the expected value of b i [k] is always 1. The instantaneous optical flux generated from Transmitter i is defined as

c i [j] = cos 2πij L + 1, (4)
so that 0 ≤ c i [j] ≤ 2
x i (t) = ∞ k=0 L-1 j=0 Φ i b i [k]c i [j]rect(t -T (j -kL) + τ i ), (5)
where Φ i is a constant for LED i to control its optical flux per chip and

rect(t) = 1 0 ≤ t < T, 0 otherwise. (6)
Therefore, M -PAM is used in each transmitter.

B. Receiver Design

The receiver receives the superposition of signals transmitted from the LEDs within its field of view (FOV). Suppose LED i has transmitted δ i messages before the receiver starts detecting signals. Due to the asynchronous transmission, there exists time lag m i between LED i and the receiver, where m i = τi T mod L. The system parameters are illustrated in Fig. 2. In the following, we present how the receiver first estimates m i and then decodes the received message.

The estimate of m i , denoted by mi , is the µ which yields the maximum correlation between B and G(c i , µ, K), i.e., mi = arg max 0≤µ≤L-1

(B•G(c i , µ, K)), (8)
where • denotes the dot product between two vectors. After the receiver estimates mi for all i, it can start decoding the messages b i [k] as follows.

2) Minimal Noise variance Decoder (MND): Recall that the received signal y(t) is sampled at a rate 1 T . To recover h i v i [n] for 1 ≤ i ≤ N and n ≥ 0, we construct a vector

y[n] = [y(nLT ) y(nLT + T ) • • • y(nLT + (α -1)T )], (9)
where α = (2γ + 1)L + 1 with a positive integer γ to be determined in Section IV-A. Due to (3) and ( 5),

y[n] = ξ ξ ξ[n]C + w, (10)
where

ξ ξ ξ[n] = {h 1 v 1 [n] . . . h 1 v 1 [n + 2γ + 1] . . . h N v N [n] . . . . . . h N v N [n + 2γ + 1] Φ 0 } (11)
is a row vector of length

(2γ + 2)N + 1 with v i [n] = b i [n+δ i ]Φ i so that ξ ξ ξ[n] contains (
2γ+2)N unknown messages together with the undetermined background light intensity Φ 0 .

Here, w is a noise vector of length α (sampling of w(t)) and C is a ((2γ + 2)N + 1) × α matrix

C = [ C 1 C 2 • • • C N 1 α×1 ] , ( 12
)
where M is the transpose of matrix M, 1 α×1 is a column vector of ones with length α and for

1 ≤ i ≤ N Ci =        S(c i , mi , L -mi ) 0 0 • • • 0 0 0 c i 0 • • • 0 0 . . . . . . . . . 0 0 0 c i 0 0 0 0 0 S(c i , 0, mi + 1)        . ( 13
)
For 1 ≤ i ≤ N , the optimal β β β * i is obtained by solving the optimization problem minimize

α-1 j=0 β β β i [j] 2 , ( 14
)
subject to

1 1×(2γ+2) • Ci • β β β i = 1, (15) Cj • β β β i = 0 (2γ+2)×1 , ∀j = i (16) 1 1×α • β β β i = 0, ( 17
) variables: β β β i ,
where 0 (2γ+2)×1 is a column vector of 0's with length 2γ + 2.

Due to (15)-( 17), our decoder gives

y[n]β β β * i = h i n+2γ+1 j=n λ j v i [j] + w • β β β * i , (18)
where λ j equals the multiplication of the j-th row of Ci and β β β * i . Due to (15), n+2γ+1 j=n λ j = 1. Thus the variance contributed by the additive Gaussian noise becomes var(w

• β β β * i ) = σ 2 w α-1 j=0 β β β * i [j] 2 . ( 19
)
So the proposed decoder not only eliminates MAI perfectly and removes the effect from background light, but also minimizes the noise variance in (19) due to ( 16), ( 17) and ( 14), respectively. The receiver builds a set S i = {y[n]β β β * i } for different values of n. When the set S i is sufficiently large, the average of the elements in S i is equal to the expected value of the right side of (18), i.e., h i Φ i . Therefore, h i can be estimated. So positioning can be done by applying algorithm such as [START_REF] Yasir | Indoor positioning system using visible light and accelerometer[END_REF].

3) Comparison with Zero Forcing Decoder: Due to (10), it seems that ZFD can be used to decode the transmitted messages. However, we are going to see that ZFD is worse than our proposed receiver in minimizing (19).

In ZFD, the receiver estimates ξ ξ ξ[n], which has (2γ+2)N +1 messages, as follows

ξ ξ ξ[n] = y[n]C † = ξ ξ ξ[n] + wC † , (20)
where ξ ξ ξ[n] is the estimate of ξ ξ ξ[n] and C † = (C C) -1 C is the pseudo-inverse of the matrix C [START_REF] Tse | Fundamentals of Wireless Communication[END_REF]. Due to (20), ZFD eliminates MAI totally. Now, consider the noise variance if ZFD is used. Define c † i as the i-th column of the matrix C † . The noise in the i-th element in ξ ξ ξ[n] due to the additive Gaussian noise becomes

σ 2 w α-1 j=0 c † i [j] 2 , (21)
which is very similar to (19). The receiver builds a set

S i = {y[n]c † i * } for different values of n, where i * = arg min (i-1)(2γ+2)≤j<i(2γ+2) α-1 j=0 c † i [j] 2 . ( 22
)
So the one with the smallest noise variance belong to LED i is chosen. In contrast, MND finds an optimal way to combine the estimates of h i v i [n] in order to minimize (19). The following example shows how the noise variance changes with i and m i .

IV. SIMULATION RESULTS

A. Simulation Setup

Monte Carlo simulation is used to evaluate the performance of our proposed MND in terms of mean square error (MSE) of the estimated channel gains ĥi and bit error rate (BER) of the decoded data vs system throughput. The performance of ZFD is also shown to compare with MND. We assume that each LED transmits the same average optical flux per chip (i.e., Φ i is a constant for all i). The receiver is equipped with a PD. LEDs start at random times and then continue transmitting. The parameters of the simulation setup and LEDs and receiver positions are listed in Table I.

B. System Performance versus System Throughput

Here, we evaluate the system performance through MSE of the estimated channel gains ĥi and BER of the decoded data vs system throughput. Assume that the chip duration T is equal to the inverse of the modulation frequency of the LEDs so that the throughput measured in bits/chip is equivalent to the throughput measured in bits/sec/Hz. To obtain different system throughput, we fix number of LEDs N = 4 and Codeword length L = 11 and vary the modulation scheme M-PAM, i.e., M = {2, 4, 8, 16, 32, 64}. Let η = log 2 M so that the system throughput = ηN L . 1) MSE of estimated channel gains ĥi : We define the MSE of the estimated channel gains ĥi as follows: Figure 3 shows the MSEs of both decoders versus the system throughput. The proposed MND outperforms ZFD because it minimizes the effect of Gaussian noise in (18) due to (14). Figure 3 also reveals that the MSEs of both decoders improve with reducing the modulation order M because smaller M gives smaller variance in data, i.e., the summation in (18).

MSE = 1 N N i=1 h i -ĥi h i 2 . ( 23
2) BER of the decoded data: Figure 4 illustrates the BERs of both decoders versus the system throughput. It is shown that MND outperforms ZFD because MND minimizes the noise impact in (18) and gives smaller MSE in the estimation of h i . The BER of MND for 0.75 Bits/chip is not shown because even 4.4×10 7 bits had been transmitted in the simulations, all the bits were correctly received. The BERs of both decoders improve with reducing the modulation order M due to the following two reasons. 1) Improving MSE and therefore more accurate thresholds for decoding are obtained. 2) Increasing the Euclidean distances among transmitted symbols.

V. CONCLUSION

This paper has proposed a multiple access scheme for visible light communications with positioning capability. The proposed scheme reduces the system cost and complexity by avoiding using central units and extra infrastructure to coordinate the transmitters. We have compared the performance of our proposed scheme with zero-forcing decoder by both analysis and simulation. We have argued why our scheme always has the smaller effect due to the additive Gaussian noise. Simulation results showed that the proposed decoder outperforms zero-forcing decoder in terms of MSE, BER and system throughput.

Fig. 1 :

 1 Fig. 1: A channel model of an VLC system for broadcasting N data streams by N transmitters which are asynchronous.

Fig. 2 :

 2 Fig. 2: The time tags among the transmitters and the receiver.

and 1 L L- 1 j=0 c i [j] = 1 .

 111 Assume that the chip duration is T . Suppose b i [k], the k-th message sent by Transmitter i, is uniformly distributed in {b , 2b , . . . , M b } where b = 2

1 )

 1 Estimating time lag m i : This is done by performing cross correlation between the received signals and a shifted version of c i as follows. Construct G(c i , µ, K) by concatenating K copies of S(c i , µ). Here µ will be determined below and the value of K will be discussed in Section IV-A. Let B = [y(0) y(T ) y(2T ) . . . y((KL -1)T )].

Example 1 .

 1 Consider a system with N = 2, L = 5, γ = 2, σ 2 w = 1, m 1 = 1 and m 2 = 2. If ZFD is used, the noise variances of LEDs 1 and 2 in (21) for all i are [198.1451 1.5515 0.2202 0.2168 0.2168 0.2169] and [197.5615 0.7271 0.2181 0.2168 0.2168 0.2170], respectively. Then ZFD selects the minimum noise variances

4 )

 4 Decoding Data: To compare the two decoders, suppose { ĥi } is obtained by either MND or ZFD. We use (20) to get an estimate of ξ ξ ξ[n]. Together with the known Φ i and { ĥi }, the data b i [n] can be estimated.

Fig. 4 :

 4 Fig. 4: BER versus system throughput.

TABLE I :

 I Simulation Parameters.belong to LED 1 and LED 2 which are 0.2168 for both LEDs. If MND is used, the noise variances are 0.049 and 0.048 for LEDs 1 and 2, respectively.

	Transmitters configuration	Values
	Number of LEDs N	4
	LEDs positions	(1,1,3), (1,2,3), (2,1,3), (2,2,3)
	Codeword length L	11
	Average transmitting optical	500 lm
	flux per chip Φ i LED half power-angle φ 1/2 Background light intensity Φ 0	π 3 rad 0.2 µA
	τ i	uniform between [0, 1]
	Receiver configuration	Values
	Receiver position	(1.5,1.5,1)
	Number of PD	1
	K	100
	γ	100
	PD's Lambertian parameter M PD's FOV PD's responsivity Rp	1.4 π 2 rad 22 nA/lux
	Receiver's area A	15 mm 2
	σ 2 noise	2.0856 × 10 -6 [V 2 ] [5]

  )

		10 -2				MND
						ZFD
	MSE	10 -3
		0.4 10 -4	0.75	1.1	1.45	1.8	2.2
					Bits/ chip
	Fig. 3: MSE of estimated channel gains ĥi versus system
	throughput.
		0.4 10 -4	0.75	1.1	1.45 Bits/ chip	1.8	2.2

τ i can also be used to capture propagation delay.