Christophe Duhamel

Philippe Lacomme

Hélène Toussaint

A GRASP  ELS for the vehicle routing problem with three-dimensional loading constraints

Keywords: Vehicle Routing, GRASP, Evolutionary local search, 3L-CVRP, 3D orthogonal packing

This paper addresses an extension of the Capacitated Vehicle Routing Problem where the client demand consists of three-dimensional weighted items (3L-CVRP). The objective is to design a set of trips for a homogenous fleet of vehicles based on a depot node which minimizes the total transportation cost. Items in each vehicle trip must satisfy the three-dimensional orthogonal packing constraints. A GRASP ELS algorithm is proposed to compute the best possible solution. We propose a new method to address the 3D packing which allows items to be rotated or not. It is based on a relaxation of the 3D problem in which items coordinates are first computed before getting compatible z-coordinates. Additional techniques are used to reduce as much as possible the time to check the 3D packing feasibility of trips. The effectiveness of our approach is evidenced through computational experiments on 3L-CVRP instances from the literature. New realistic instances are also proposed. These instances are based on the 96 French districts and encompass both small scale instances and large scale instances with up to 200 nodes

Introduction

Capacitated Vehicle Routing Problem and extensions with packing constraints

The Capacitated Vehicle Routing Problem (CVRP) is a classical NP-hard node routing problem which received a considerable amount of attention for decades [START_REF] Baldacci | Routing a Heterogeneous Fleet of Vehicles[END_REF] [2] [START_REF] Prins | Two memetic algorithms for heterogeneous fleet vehicle routing problems[END_REF]: it consists in optimally organizing vehicles trips in order to deliver goods required by a set of clients. It can be fully defined by considering a depot and a set of clients. Each one corresponds to a node of a complete graph where V is a set of n+1 nodes, 0 being the depot and nodes 1...n being the clients. Each edge has a finite cost and each node is given a demand . A fleet of homogeneous vehicles of limited capacity is located at the depot. The objective is to design a set of trips of minimal total cost to service all clients. A trip is a cycle performed by one vehicle. It starts at the depot, visits a subset of nodes, before returning to the depot. The trip total load is upper bounded by the vehicle capacity. Since split deliveries are not allowed, each client is serviced by exactly one vehicle. As stressed in [START_REF] Cordeau | New heuristics for the vehicle routing problem[END_REF], exact methods can only solve small to medium instance. Thus, medium and large CVRP instances are typically addressed by metaheuristics.

The 2L-CVRP is an extension of the CVRP which includes two-dimensional orthogonal rectangle loading constraints (the 2L constraints). This problem is essentially addressed in [START_REF] Gendreau | A Tabu Search Heuristic for the Vehicle Routing Problem with Two-Dimensional Loading Constraints[END_REF][6] [START_REF] Fuellerer | Ant colony optimization for the two-dimensional loading vehicle routing problem[END_REF]. It can be reduced to the CVRP when the size of the items is not considered or when items are 1 1 squares, thus dealing only with their weight. The 2L-CVRP resolution has been first addressed by Iori et al. [START_REF] Iori | An exact approach for capacitated vehicle routing problems with two-dimensional loading constraints[END_REF] using a branch and cut approach limited to small scale instances (less than 25 clients). Then Gendreau et al. [START_REF] Gendreau | A Tabu Search Heuristic for the Vehicle Routing Problem with Two-Dimensional Loading Constraints[END_REF] introduced a tabu search algorithm. Zachariadis et al. [START_REF] Zachariadis | A guided Tabu Search for the Vehicle Routing Problem with two dimensional loading constraints[END_REF] developed a guided tabu search. Fuellerer et al. [START_REF] Fuellerer | Ant colony optimization for the two-dimensional loading vehicle routing problem[END_REF] proposed an efficient version of the Ant Colony scheme to solve the 2L-CVRP. Recently, Duhamel et al. [START_REF] Duhamel | A multi-start evolutionary local search for the two-dimensional loading capacitated vehicle routing problem[END_REF] introduced a multi-start evolutionary local search scheme which outperforms all previous published methods. The approach is original as it does not address the 2L-CVRP during the main optimization process but rather a relaxation into the so-called RCPSP-CVRP.

In the RCPSP-CVRP, the two-dimensional packing problem is relaxed into a RCPSP: at each point of the vehicle length the total width used must not exceed the vehicle width. Thus the vehicle width is related to the RCPSP resource availability. At the end of the main optimization process, the RSPCP-CVRP solution is transformed into a 2L-CVRP solution by a dedicated procedure. The authors showed in their experiments that most of the RCPSP-feasible solutions can be efficiently transformed into 2L-CVRP feasible solutions by only considering packing solutions which satisfy the previously computed x-abscissa.

The three-dimensional loading CVRP (3L-CVRP) is an extension of the 2L-CVRP where the height is also considered. More formally, each vehicle of the homogenous fleet is now defined by a weight capacity D and by a volume where is the vehicle length, is the vehicle width and is the vehicle height (related to (x, y, z) coordinates). The demand of each client consists of a set of items of total weight . Each item is a three-dimensional cuboid of length l ik , width w ik and height h ik . Each client must be serviced by exactly one vehicle, which is assigned to a single trip. A trip is a sequence of clients where corresponds to the depot. Each trip must be both "weight-feasible" and "packing-feasible". A trip is "weight-feasible" if the total weight of carried items does not exceed the vehicle capacity, i.e.

. It is "packing-feasible" if the client items can be loaded into the vehicle without overlapping and if it satisfies the classical orthogonal three-dimensional packing constraints. A set of "weight-feasible" and "packing-feasible" trips which involves all the clients defines a solution of the 3L-CVRP.

The 3L-CVRP has been addressed by Gendreau et al. [START_REF] Gendreau | A Tabu Search Algorithm for a Routing and Container Loading Problem[END_REF] and more recently by Fuellerer et al. [START_REF] Fuellerer | Metaheuristics for vehicle routing problems with threedimensional loading constraints[END_REF]. Only medium instances have been considered since three-dimensional packing problems are much harder to solve than their two-dimensional counterparts. The seminal publication of [START_REF] Gendreau | A Tabu Search Algorithm for a Routing and Container Loading Problem[END_REF] introduces a tabu search algorithm that iteratively invokes a tabu search procedure for solving the inner loading sub-problem. Fuellerer et al. [START_REF] Fuellerer | Metaheuristics for vehicle routing problems with threedimensional loading constraints[END_REF] introduce a highly efficient ant colony optimization algorithm which takes advantage of both fast packing heuristics for the loading sub-problem and of effective heuristics for the routing problem. These two publications also consider additional constraints about item fragility, LIFO unloading and support. Note that both instances from the literature and real-world instances were used by [START_REF] Gendreau | A Tabu Search Algorithm for a Routing and Container Loading Problem[END_REF] to evaluate the performance of their method.

Cutting and Packing problems

General Cutting and Packing problems

Packing problems belong to the well-known family of cutting and packing problems. Many packing problems deal with the insertion of rectangular items in a rectangular bin in both two and three dimensions. They mostly differ on the objective function to optimize.

 The Three-Dimensional Bin Packing Problem (3BPP) consists in packing a set of rectangular boxes into a minimal number of identical rectangular boxes [START_REF] Hifi | A Linear Programming Approach for the Three-Dimensional Bin-Packing Problem[END_REF] [13];  The Three-Dimensional Strip Packing Problem (3SPP) consists in packing a set of rectangular boxes into a strip of known width and infinite height so as to minimize the overall height of the packing [START_REF] Bortfeldt | A heuristic for the three-dimensional strip packing problem[END_REF] [START_REF] Allen | A hybrid placement strategy for the three-dimensional strip packing problem[END_REF];  The Three-Dimensional Packing Problem (3PP) consists in checking if a set of rectangular boxes can be packed into one bin (rectangle box) of fixed size, see [START_REF] Almeida | A particular approach for the Three-dimensional Packing Problem with additional constraints[END_REF] for instance.

Several extensions have also been addressed over time, including but not limited to, rotation of items, limitations on the total weight and/or item costs.

The 3D packing sub-problem in the 3L-CVRP

The packing problem within the 3L-CVRP falls into the last category (3PP) since each trip has to be "packing-feasible". A 3PP instance consists of a set of items which have to be packed into a bin of length , of width and of height . An item i has a length l i , a width w i and a height h i ( .

A 3PP solution can be fully defined by the position of each item i, denoted (x i , y i , z i ), into the bin. This position corresponds to the coordinates of its bottom-left corner. Item rotation is only allowed in the (x, y) plane as rotations in other planes may be prohibited in the corresponding real-life application (items often have a "top" side for instance). Moreover the packing must be orthogonal, i.e. the items must be placed with their edges parallel to the sides of the bin.

Some authors have added extra constraints:

-fragility: the items tagged as "fragile" cannot be put under another item;

-support: each item must have a minimum "supporting area", i.e. a given percent of its basis must be defined by the top of other items (or by the floor of the bin);

-LIFO: the items of any client in the trip can be unloaded by only using straight movements, i.e. the items of a client i are not blocked by items of yet unvisited clients.

Such constraints correspond to realistic considerations in the industrial context of transportations and logistics. They are mandatory in many situations as CVRP solutions involving fully-loaded or nearly fully-loaded vehicles may not be 3L-CVRP feasible in practice, thus greatly reducing the interest of many CVRP commercial solvers.

GRASP ELS framework for the 3L-CVRP

GRASP ELS Principle

The GRASP ELS [START_REF] Prins | A GRASP×Evolutionary Local Search Hybrid for the Vehicle Routing Problem[END_REF] is a hybridization of the GRASP metaheuristic and of the ELS metaheuristic combining the positive features of both methods. The GRASP (Greedy Randomized Adaptive Search Procedure) [START_REF] Feo | Greedy randomized adaptive search procedures[END_REF] is a multi-start Local Search metaheuristic. At each iteration, an initial solution is constructed by using a greedy randomized heuristic. It is then improved by a local search and the best solution obtained at the end of each GRASP iteration is kept. The ELS (Evolutionary Local Search) [START_REF] Prins | A simple and effective evolutionary algorithm for the vehicle routing problem[END_REF] is an extension of the ILS (Iterated Local Search, [START_REF] Lourenço | Iterated Local Search[END_REF]). At each iteration of the ELS, several copies of the current solution are done. Each copy is modified (mutation) before being improved by a Local Search. The best resulting solution is kept as the new current solution. The purpose of the ELS is to better investigate the neighbourhood of the current local optimum before leaving it, while the GRASP aims at managing the diversity during the solution space exploration. The framework we promote is a multi-start ELS in which the ELS is applied to the initial solutions generated by greedy randomized heuristics. Such an approach can also be viewed as a GRASP ELS in which the ELS is used as Local Search. Besides combining GRASP with ELS, another important feature of our approach is the alternation between two solution spaces: the giant tour space and the 3L-CVRP solution space. By defining genuine exploration on those two search spaces and by defining projections from one search space into the other one, one can more easily avoid being trapped in local optima. The high quality solutions obtained by Prins [START_REF] Prins | A simple and effective evolutionary algorithm for the vehicle routing problem[END_REF] for the VRP, alternating between two search spaces (giant tour and VRP solutions) is a clear illustration of approaches which manage alternation between a set of giant tours and a set of solutions.

Two solution representations are used: solutions encoded as giant tours (TSP tours on the n clients) and 3L-CVRP solutions encoded as the set of trips (see Figure 1). where is the i th trip and where is the number of trips in T. It first generates a new giant trip by modifying the concatenation order. Then some clients are exchanged to get the new giant trip

Proposal for a new vehicle loading resolution approach

The approach we propose shares some similarities with the method we developed for the 2D packing problem in the 2L-CVRP [START_REF] Duhamel | A multi-start evolutionary local search for the two-dimensional loading capacitated vehicle routing problem[END_REF]. For the 2L-CVRP, the original 2PP is first relaxed into a RCPSP with one resource, leading to the RCPSP-CVRP. A solution to the RCPSP-CVRP is then computed before being transformed back into a 2L-CVRP by using an efficient procedure. In most of the cases, the resulting 2L-CVRP solution is packing-feasible which means no other subsequent RCPSP-CVRP solution has to be investigated.

Unfortunately similar idea cannot be successfully applied to the 3L-CVRP. One should think that relaxing the 3PP sub-problem into a RCPSP with two resources (for example the width and the height) would also lead to the RCPSP-CVRP and most of the previous work could be re-used as well. However, the transformation of a RCPSP solution into a 3PP solution is often not possible as all the items are likely to be packed at the same location. Thus we propose a variation based on a 2-step procedure.

General process to solve the 3PP

Let be a set of items. The following two steps are performed to compute a solution to the 3PP:

-Step 1: (x i , y i ) positions are computed for each item i. The 3D geometry of the items is relaxed and the height of the item is considered as a cost c i = h i . Thus the following sub-problem has to be solved: "Let I be a set of rectangular items i defined by their length l i , their width w i and their cost c i , and let a rectangular bin be defined by its length L, its width W and its capacity C. Find a position (x i , y i ) for each item i of I in the bin such that (i) the packing is orthogonal, (ii) the sum of the overlapping items costs does not exceed C". This step is addressed in part 2.2.2.

-Step 2: given the (x i , y i ) positions obtained in Step 1, the z i coordinates are computed such that (x i , y i , z i ) positions lead to a 3PP solution for the set of items I. Thus a 3PP has to be solved in this step, except that the solution is already partially defined. The resolution is fully detailed in part 2.2.3.

To the best of our knowledge, this kind of approach is original. However Gilmore and Gomory proposed in 1965 a stack building approach [START_REF] Gilmore | Multistage cutting stock problems of two and more dimensions[END_REF]. It consists in packing items stack after stack by solving a two-dimensional packing problem for each stack. The method we introduce is quite different since it does not solve as many two-dimensional packing problems. In fact, only one problem need to be solved in step 1 (which can be seen as a 3PP relaxation and not as a 2PP) and the solution is then transformed into a 3PP solution in Step 2.

Step 1: solving the relaxed 3PP

As stressed in section 2.2.1, the arrangement problem introduced in Step 1 is considered. It is defined as follows: let I be a set of rectangular items i defined by their length l i , their width w i and their cost c i . Let a rectangular bin be defined by its length L, its width W and its capacity C. The problem consists in finding a (x i , y i ) position for each item i of I in the bin such that (i) the packing is orthogonal, (ii) the sum of the overlapping items costs does not exceed C.

The arrangement problem has to be solved each time the packing feasibility is checked. Since the check has to be done each time a solution is modified, its time efficiency is crucial. Thus, for time efficiency, we propose a greedy (heuristic) approach where items are scanned in an ordered list O. The items in O are considered and tentatively placed into the bin while satisfying constraints (i) and (ii). This process is done by the Solve_x_y_coordinate procedure (see Algorithm 1).

The Solve_x_y_coordinate main loop uses a current position in the bin denoted by (posx, posy). It tries to pack as many items from O as possible at this position. Any successfully packed item from O is removed from O. The (posx, posy) position is first initialized at the origin (0, 0). It is then updated according to an increasing order of x-coordinates and y-coordinates. The way (x, y) coordinates are scanned allow us to state that an item i can be packed at the position (x, y) if:

where is the sum of the items costs which are overlapping at the position (x, k).

The way the positions are scanned in the arrangement is crucial. One must look for empty spaces reduction above the items while limiting the items stow in order to be able to successfully solve the 3 dimensional packing in the following Step 2.

Algorithm 1: packing items in step 1

The main drawback of this approach is its greediness (heuristic). This means the local choices may lead to a packing failure although packing could be done. To prevent such wrong answers, one could consider a backtracking mechanism (like a tree search). However this would be computationally too expensive since Solve_x_y_coordinate is called a lot of times during the GRASP process. A partial workaround based on a look-ahead mechanism has been added. It consists in adding an extra condition when trying to pack one of the last three items from O: the candidate item i can be packed at the position (posx, posy) only if the remaining items from O can be packed afterwards. Setting a limit of three remaining items has experimentally shown to be a good compromise between efficiency and time consumption.

A post processing step consists in spreading items over the bin. Indeed the way x and y coordinates are scanned leads to the items being packed as long as at the bottom-left side of the bin. As a consequence, the opposite area (top-right part of the bin) is not exploited the best possible way. Thus packed items are scanned in the decreasing order of their right edge position. Each item is then shifted as much as possible to its right (without introducing new overlaps). The same process is applied on y coordinates. This step reduces the number of overlapping items and makes the problem at step 2 easier to be solved.

Step 2: solving the 3PP using the partial solution computed at step 1

This step aims at computing a solution to the 3PP by using the solution found at Step 1. It consists in computing the position of the items. The x and y positions have already been computed in Step 1. The idea is to scan the coordinates, starting from 0. For each value, as much items as possible are packed respecting their position. This process ends when all items are packed or when the top of the bin is reached. The Solve_z_coordinate procedure is fully described in Algorithm 2.

Whole packing feasibility check

As previously mentioned, a trip is feasible if (i) the total weight of the clients items does not exceed the vehicle capacity and if (ii) the items can be packed into the vehicle with respect to the 3PP constraints. Checking the first constraint is trivial. Checking the second constraint is trickier and we use the method described above. The global check is done by the 3D_Check_trip procedure (see Algorithm 3). The procedure iteratively generates an ordered list O before checking it. It stops as soon as a packing has been found or when the maximal number of attempts has been reached. The procedure Solve_x_y_coordinate tries to identify a packing which relies on the ordered set O. Upon success, Solve_z_coordinate is called. Otherwise, the Random_Neighboord_Generation generates a new list O' by randomly exchanging some items in O. Rotations are addressed by a random selection of item in O and by swapping their length and width.

Preliminary computation and storage

A lot of trips are evaluated during the optimization process. Moreover, same trip can be evaluated several times. Thus, a way to save time consists in avoiding unprofitable calls to 3D_Check_trip (several runs with identical parameters) by saving the result (true or false) of each trip feasibility check. A dedicated data structure is used and it is updated along the GRASP ELS process.

A combination of data structures can be introduced: three matrices are dedicated to trips which deliver from 2 to 4 costumers. Trips with a single client are trivially feasible, unless the instance is unfeasible. Note that items for one customer can be packable or not depending if items rotations are allowed or not. These matrices provide a O(1) check if the trip has already been checked, either being packingfeasible or not. Otherwise the feasibility check is performed and the result is stored into the corresponding matrix. The major drawback is the huge memory footprint, especially for the last 4dimensional matrix. Another data structure is used for trips involving more than 4 clients. It is a redblack tree (self-balancing binary search tree), see the seminal contributions [START_REF] Guibas | A Dichromatic Framework for Balanced Trees[END_REF] [START_REF] Bayer | Organization and Maintenance of Large Ordered Indices[END_REF]. In associative data structures, each element is associated to a key which is used to find it back. Here the key corresponds to the set of clients of the trip without any relative order consideration. In order for the storage to be efficient, the relation between the keys and the trips should be as close as possible to a 1-1 correspondence. We propose the following key computation: given a trip , its key is generated by first computing the number of clients n(t) in the trip. Then the client identification numbers are concatenated in the increasing order, leading to a value . For example, if , then and . Such an order is total since it is always possible to compare two different trips and :

The search in a red-black tree is done in O(log(n)) where n is the size of the tree.

The Load_Resolution procedure (see Algorithm 4) is in charge of evaluating a trip. This happens if the trip has never been evaluated or if it has been submit to less than p unproductive attempts there have been less than p failed evaluation (packing) attempts. For convenience, Store(t) denotes the storing evaluation of the trip . It is independent of the structure used to store the trip. Store(t)has the following meaning:

3D packing resolution example

Let us consider the instance E023-05s.DAT from [START_REF] Gendreau | A Tabu Search Algorithm for a Routing and Container Loading Problem[END_REF]: 5 clients have to be serviced for a total of 12 items, detailed in Table 1 Let us consider the ordered list which leads to an arrangement solution. The next figures (from Figure 2 to Figure 7) illustrate the evolution of the arrangement process at different steps. The large rectangular area ( corresponds to the bin while the small rectangles inside it are the items already packed. The number in the small rectangles is the total cost for the associated area of the bin. Let us remind that the item cost corresponds to its height. The limit on the cost (the bin height) is set to 30. For each figure, the last packed item is filled with dotted lines. No items can be packed in (0, 8). Thus the next position investigated is (0, 12) and all the remaining items in the list are scanned: the first packable item is and the second one is . Then the position (0, 18) is eligible for packing , which leads to the arrangement in Figure 3. The method skips to abscissa 14, considering positions (14;0), (14;8) and (14;12). The item can be placed at (14;12) leading to the packing solution of Figure 4. No item can be put at the next positions investigated. The first interesting position is [START_REF] Feo | Greedy randomized adaptive search procedures[END_REF][START_REF] Hifi | A Linear Programming Approach for the Three-Dimensional Bin-Packing Problem[END_REF] where item can be placed. The Solve_x_y_coordinate procedure has produced a compact arrangement and the computed position for each item are given in Table 2

Items shift

Shifting the items is done iteratively along the x-axis and then along the y-axis until no further shift can be done. This process leads to the new items coordinates in Table 3. The sequence from Figure 8 to Figure 13 illustrates the way the packing is built by Solve_z_coordinate. The ordered set of items is . First coordinate is investigated and as many items as possible are packed at this current z according to their (x, y) position and according to the O. Thus and are packed at (see Figure 8). The current is updated to the smallest available height, i.e.

Items B1 B2 B3 B4 B5 B6 B7 B8 B9 B10 B11 B12 x-coordinate 24
. Items and can be packed leading to the partial vehicle load shown in Figure 9. This process goes on with z = 9 where only can be packed (Figure 10), where only is packed (see Figure 11) and finally where and are packed. This leads to the final items packing in Figure 13.

Packing the items according to increasing values on the z-coordinate strategy usually produces dense layers with as much items as possible packed at the same time.

The final 3D-loading solution is shown in Table 4 Split first builds an auxiliary digraph H T = (X, Y, Z) where X is a set of n+1 nodes indexed from 0 to n. Node 0 is a dummy node, while the nodes 1…n correspond to the client sequence of the giant tour . An arc (i,j) belongs to Y if a trip servicing clients v i+1 to v j (included) is both weightfeasible and 3D-feasible. The weight of the arc corresponds to the trip cost . Optimally splitting T can be done by computing a min-cost path from node 0 to node n in H. An initial label is set at node 0. The labels are then propagated from node to node in H using the arcs. The best label at node n is kept as the optimal split. Let be the p th label assigned to node i. It corresponds to a feasible split of the initial clients t 1 ...t i into trips.

is the number of vehicles still available, is the cost of the trips previously built and is the reference to its father label, e.g. , the k th label at node j. The initial label at node 0 is defined as . It corresponds to the empty solution where all the vehicles are available. Propagating the label along the arc produces the label the following way:   Since a lot of labels are generated and stored at each node, the computational time can quickly grow. Thus dominance rules must be defined in order to keep a good time efficiency. A label is said to dominate the label if one of the following conditions holds:

The critical path leading to the best final label defines the trips of the 3L-CVRP solution. The procedure Split is detailed in Algorithm 5. For each node i, NB [i] gives the number of associated labels. The procedure Check_Domination_On_Node checks if the new label L is dominated by another label at node j. The procedure Insert inserts this label into the set of labels from node j and removes the dominated labels. The number of labels is updated accordingly.

1.

procedure Split 2.

input parameters 3.

T: giant tour 4.

output parameters 5.

S: 3L-CVRP solution 6.

global parameter 7.

D : maximal vehicle weight capacity 8.

V : vehicle volume 9.

d i : total items weight of client i 10.

v i : total volume of items located at client i 11.

c ij : cost from client i to j 12.

n : number of clients 13. begin 14.

, S :=  15. pos_last := 0 16. for i := 1 to n do L i :=  endfor 17. for i := 0 to n -1 do 18.

j := i + 1 Table 5: parameters setting for the classical instances

Implementation and classical benchmarks used

We report results on the set of instances used in [START_REF] Gendreau | A Tabu Search Algorithm for a Routing and Container Loading Problem[END_REF] and then in [START_REF] Fuellerer | Metaheuristics for vehicle routing problems with threedimensional loading constraints[END_REF]. The number of clients varies from 15 to 100 and the total number of boxes varies from 32 to 198. The number of vehicles varies from 5 for the small instances to 28 for the largest ones. These instances can be downloaded at http://www.or.deis.unibo.it/research_pages/ORinstances/. [START_REF] Fuellerer | Metaheuristics for vehicle routing problems with threedimensional loading constraints[END_REF] and with the Tabu Search of [START_REF] Gendreau | A Tabu Search Algorithm for a Routing and Container Loading Problem[END_REF].

The GRASP  ELS is a random search algorithm. To provide a fair comparative study with Fuellerer's proposal, each instance has been solved ten times, the same way they did in their experiments. We report the average cost as well as the average CPU time to get the best solution over the 10 replications. Note that the best found solution over the 10 runs is also kept with the corresponding CPU time to reach it. The computational time of each method has been scaled by the speed factor presented in Table 7. This coefficient takes into account the MIPS performance of each processor.

Gendreau et al. [START_REF] Gendreau | A Tabu Search Algorithm for a Routing and Container Loading Problem[END_REF] Fuellerer et al. [ All previously published methods were benchmarked over 1 hour of computational time, i.e. 1 hour of computation is assigned for one run of the methods. Since the reference results [START_REF] Fuellerer | Metaheuristics for vehicle routing problems with threedimensional loading constraints[END_REF] have been obtained on a computer which is 1.5 times faster than ours, the GRASP  ELS time limit is set to 1h30.

Average results for 3L-CVRP instances

A summary of the results is presented in Table 8 for the three methods. For each method, the number of time the method gives the best published method (line 2), the number of time the method compete with the best one (line 3) and the number of time the method is worst (line 4) are reported. The results show that the GRASP  ELS find the best solution for 16 out of 27 instances and outperforms both the Tabu search from Gendreau et al. [START_REF] Gendreau | A Tabu Search Algorithm for a Routing and Container Loading Problem[END_REF] and the Ant colony Scheme from Fuellerer et al. [START_REF] Fuellerer | Metaheuristics for vehicle routing problems with threedimensional loading constraints[END_REF]. The average value 846.1 is also the best. Two versions of GRASP ELS are provided in order to evaluate the consequence of allowing items rotation or not. Quite surprisingly, forbidding items rotation does not deteriorate that much the solution. On average, our method produces solutions that are 3 units higher. Thus rotations do not seem to play a significant role for this set of instances. When carefully checking the results (see the Appendix), one can note that the solution is the same, with and without items rotation, for half the instances.

Hash function performances on results

The hash function kept in memory the 3D packing results leading to a learning algorithm in order to save time during process. The saving time increased over replications since the hash tables are not erased between replications.

Impact of hash function can be easily evidenced in numerous instances including the instance 07 shown in Table 9. For this instance the limit is fixed to 63000 iterations. The total time to perform those iterations is about 3079.4 s in the first replication. It quickly drops in the second iterations (1475.4 s). The time difference corresponds to the packing results kept in memory in the first replication since both replications have exactly the same parameters. We can notice that the total time decreases over the 10 replications, dropping from 3000 s to 700 s. Table 9: GRASP ELS performances over the 10 iterations

Example of a 3L-CVRP solution

Let us consider the instance 08 with 22 clients to service, 43 items to load and 8 vehicles available.

The GRASP  ELS provides a solution of value 730 which is better than the solution given by the Ant Colony Scheme [START_REF] Fuellerer | Metaheuristics for vehicle routing problems with threedimensional loading constraints[END_REF]. This solution is made of 5 trips: Figure 14 provides a graphical representation of the trips. For each trip, table 10 reports the total items weight, the total item volume and the trip cost. Let us note that a 3D visualization tool can be obtained at http://www.isima.fr/~toussain/.

New benchmarks

Using the GIS system developed by Bajart and Charles [START_REF] Bajart | Systèmes d'Information Géographique. 3 rd year project report[END_REF], shortest paths are computed between cities with more than 100 or 500 citizens for the 96 French counties. The shortest paths are computed using the Google web service and they correspond to the roadmap distance in kilometers between cities. Thus, 96 realistic instances are provided in terms of distances, with size varying from 60 to 255 nodes. To the best of our knowledge, those are the first available instances based on real counties. They can be divided into 4 subsets: -DLT_3LCVRP_1: 13 small instances with less than 100 nodes; -DLT_3LCVRP_2: 40 medium instances with 100 to 150 nodes; -DLT_3LCVRP_3: 33 large instances with 150 to 200 nodes; -DLT_3LCVRP_4: 11 very large instances with more than 200 nodes.

The results for the 96 instances are available in Appendix 2.

Table 12: GRASPELS performance for the new instances

For the classical instances, allowing items rotation slightly improves the results. For this new set of instances, eight instances cannot be solved if rotations are forbidden since the items of some clients cannot be packed with the heuristic we introduced (see

Concluding remarks

This article considers an extension of the well-known CVRP in which three dimensional packing constraints must be addressed in each trip servicing clients. This problem deals with two combinatorial optimization problems: vehicle routing and three-dimensional packing. The method we propose compete with the best published methods but the method is currently dedicated to the 3L-CVRP with no extra constraints. It is based on an original resolution of the 3PP based on a dedicated heuristic for the vehicle loading resolution. We are currently investigating the 3L-CVRP with additional constraints, trying to extend the original 3D-packing scheme we introduce.

  y i : y-position of item i z i : z-position of item i ok : boolean (true upon success) begin O := items from cli k := 1, l := 1, j := 1 //number of iterations ok := false while (k < nm) and (ok = false) //main loop while (l < nm1) && (ok = false) //search for x and y coordinates O := Random_Neighboord_Generation(O) (ok, x, y) := Solve_x_y_coordinate (O, V) l := l+1 endwhile if (ok = true) then //search for z coordinate ok = false while (j < nm2) and (ok = false) O := Random_Neighboord_Generation(O) (ok, z) = Solve_z_coordinate(O, x, y, V) j

Figure 2 :

 2 Figure 2: putting the first three items

  Figure 3: adding , and

  Figure 4: adding

  Figure 5: adding

Figure 6 Figure 7 :

 67 Figure 6: adding

Figure 8 :Figure 9 : 8 Figure 10 : 9 Figure 11 : 17 Figure 12

 898109111712 Figure 8: Packing items at z = 0 Figure 9: Packing items at z = 8

Algorithm 5 :

 5 check := (trip_load  D) and (trip_volume 9 (j  pos_last) and (size > 1) then 40.res := Solve_3D(set_boxes) 41.else res := true 42. endif 43.if (res = true) then // 3D packing successfully solved 44.for p := 1 to NB i do Split for the 3L-CVRP 3 Computational experiments All procedures have been implemented in C++ and compiled using g++. Numerical experiments have been carried out on a 2.1 GHz Opteron computer running Linux operating system. The CPU power has been evaluated at around 4140 Mflops/s. The numerical experiments are based on two instance sets:  a set of instances previously introduced in [10];  a new set of instances based on the 96 French counties. To the best of our knowledge, this is the first step towards the definition of realistic and available instances for the 3L-CVRP. They are available for further experiments at http://www.isima.fr/~toussain.

Trip 1 :

 1 Depot, 14, 17, 22, 20, 19, Depot Trip 2: Depot, 11, 13, 9, 5, 4, 7, Depot Trip 3: Depot, 16, 15, 3, 2, 1, 6, 12, Depot Trip 4: Depot, 21, 8, 10, Depot Trip 5: Depot, 18, Depot

Figure 14 :

 14 Figure 14: Solution for the instance 08

Figure 1: GRASP ELS with alternation between the two search spaces

	problems. As a giant trip is not a direct representation of a 3L-CVRP solution, we have chosen the
	inner ELS to work on 3L-CVRP solutions while GRASP focuses on giant tours.
	np GRASP iterations			Start
	Search space	Search space
	of giant trips		of solutions
					Randomized
					Heuristic
	ne ELS iterations	T (giant trip)	Concat	S (solution : set of trips)
			Mutation
			T' (giant trip)	Split	S' (solution : set of trips)
	nd	neighbourhood	iterations		Local Search
				Concat	S'' (solution : set of trips)
			T'' (giant trip)
			Update the set of nd giant trips (T'')
	Best giant trip
				end
	A random heuristic is required to generate an initial solution S (set of trips) at each iteration of
	GRASP. It is then transformed into a giant trip T before being perturbed in a way similar to the
	mutation operator in Genetic Algorithms. The resulting giant tour is split into 3L-CVRP trips which
	provides a solution S'. Then S' is improved using a Local Search operating on 3L-CVRP trips. The
	new solution S'' is associated to the giant trip T'' by trips concatenation and it becomes the incumbent
	solution (S,T). During ELS, nd "children" are generated out of S, each one being mutated and
	improved by the local search. The best child replaces S. The process is iterated until ne iterations are Converting a 3L-CVRP solution into a giant tour is done by the Concat procedure. It consists in removing the depot from each trip and done. The incumbent solution is updated before starting a new GRASP iteration.
	then concatenating the resulting trips into a single one. The reverse operation, i.e. converting a giant The Local Search is defined as a first improvement descent method using several classical VRP tour into a 3L-CVRP, requires more work. It is usually done by a dedicated splitting procedure (Split) neighborhoods to improve the initial 3L-CVRP solution: 2-Opt within a trip, 2-Opt between two trips, and it relies on dynamic programming. Such an approach has been successfully applied to numerous routing problems including the Capacitated Arc Routing Problem, the Vehicle Routing Problem, the Swap within a trip and Swap between two trips.
	Location Routing Problem for instance, see [21] for a recent state of the art of Split in routing

The random heuristic is indeed a randomized version of both the Path-Scanning heuristic and the heuristic of Golden et al. Thus, each call is likely to produce a different solution. The mutation operator is defined on the giant tour ,

Algorithm 2 : computing z coordinate (step 2)

	1.	procedure Solve_z_coordinate
	2.	input parameters
	3.	O : ordered set of items
	4.	x : set of positions in x (x i = x-position of item i)
	5.	y : set of positions in y (y i = y-position of item i)
	6.	B : bin
	7.	output parameters
	8.	ok : boolean (true upon 3BPP success)
	9.	z : set of positions in z (z i = z-position of item i)
	10.	local parameters
	11.	h : array [1…L][1…W]	//h[x][y] = height already reached at (x,y)
	12.	begin
	13.	z := 0
	14.	ok := true
	15.	while ( ok = true ) do
	16.	for ( k := 1 to Card(O) ) do
	17.	item := O[k]
	18.	if ( item can be packed in position (item.x, item.y, z) ) then
	19.	update h
	20.	z item := z
	21.	remove item from O
	22.	endif
	23.	if (z + item.height > B.height) then
	24.	ok := false
	25.	endif
	26.	endfor
	27.	endwhile
	28.	end
	29.

endif endcase endswitch end Algorithm 4: Vehicle Load Resolution

	1.	procedure Load_Resolution
	2.	input parameters
	3.	t : trip
	4.	p : number of 3D trip evaluation attempts
	5.	nm, nm1, nm2 : maximal number of attempts for 3D_Check_trip procedure
	6.	V : vehicle (bin)
	7.	output parameters
	8.	ok : boolean (true upon success)
	9.	local parameters
	10.	cli : set of clients in trip t
	11.	begin
	12.	ok := false
	13.	switch case:
	14.	case Store(t) = 1
	15.	ok := true
	16.	endcase
	17.	case Store(t) = -p
	18.	ok := false
	19.	endcase
	20.	case (Store(t) ≠ -p) and (Store(t) ≠ 1)
	21.	(x, y, z, ok) = 3D_Check_trip (cli, nm, nm1, nm2, V)
	22.	if (ok = true) then
	23.	Store(t) := 1
	24.	else
	25.	Store(t) := Store(t) -1
	26.
	27.
	28.
	29.

Table 1 : set of items to pack 2.3.1 Solving the arrangement problem (Solve_x_y_coordinate)

 1 .

		1	2	3	4	5	6	7	8	9	10	11	12
	Client	C20	C20	C20	C1	C13	C13	C13	C7	C7	C22	C22	C22
	Items	B1	B2	B3	B4	B5	B6	B7	B8	B9	B10	B11	B12
	Length	36	29	29	34	14	24	15	15	22	18	22	19
	Width	10	10	8	10	9	7	10	11	6	13	8	12
	Height	10	8	9	13	11	7	8	17	12	11	11	17

Table 2 : Items position after resolution of the arrangement problem

 2 .

	Items	B1	B2	B3	B4	B5	B6	B7	B8	B9	B10	B11	B12
	x-coordinate	18	0	0	19	0	0	19	34	36	0	14	0
	x-coordinate	12	8	0	0	12	18	0	0	10	12	12	0

Table 3 : Items position after items shift 2.3.3 Items packing in z (Solve_z_coordinate procedure)

 3

			1	1	26	2	0	30	45	38	6	16	7
	x-coordinate	15	8	0	5	16	18	5	4	19	12	17	0

Table 4 : 3D-packing solution 2.4

 4 . Split procedure As previously mentioned, Split is a key-procedure which converts a giant tour into a 3L-CVRP solution (with respect to the sequence). It is based on the classical Split procedure[START_REF] Prins | A simple and effective evolutionary algorithm for the vehicle routing problem[END_REF][25][START_REF] Beasley | Route-first cluster-second methods for vehicle routing[END_REF], tuned to address the specific 3L-CVRP constraints.

	Items	B1	B2	B3	B4	B5	B6	B7	B8	B9	B10	B11	B12
	x-coordinate	24	1	1	26	2	0	30	45	38	6	16	7
	y-coordinate	15	8	0	5	16	18	5	4	19	12	17	0
	z-coordinate	19	0	0	17	8	0	0	0	0	19	8	9

Table 5

 5 gives the set of parameters used for the two set of instances.

		Parameters definition	Parameters value
	np	number of GRASP iterations	60
	ne	number of ELS iterations	15 + min(6, nbVehicule)
	nd	number of neighborhoods	10
	p	maximal number of 3D trip evaluation	5

Table 6 : instances characteristics

 6 The details of the GRASP  ELS solutions are available at http://www.isima.fr/~lacomme, http://www.isima.fr/~toussain and http://www.isima.fr/~duhamel. The GRASP  ELS is compared with the Ant Colony Scheme of

Table 7 : comparative performance of processors

 7

Table 8 : average GRASPELS performance, with and without rotations

 8

		Gendreau et al. [10]	Fuellerer et al. [11] GRASP  ELS GRASP  ELS
	rotation	yes	yes	yes	no
	nb best	0	2	16	14
	nb of equal	7	7	8	6
	Nb of worst	20	18	3	7
	avg value	876.31	856.7	847.04	848.88
	best value	?	?	841.96	845.48

Table 10 : trips details

 10 Let us consider the trip 4. It consists in servicing clients 21, 8 and 10. Table11gives the list of the boxes for each client, along with their dimensions.

	Client 21	Client 8	Client 10
	box 1: 24 15 8	box 1: 36 11 13	box 1: 18 11 8
	box 2: 13 14 14	box 2: 27 11 17
		box 3: 34 7 16

Table 11 : list of boxes to pack for each client

 11

	Client 21, box 1: (31;0;0)
	Client 21, box 2: (18;11;0)
	Client 8, box 1: (0;0;17)
	Client 8, box 2: (0;0;0)
	Client 8, box 3: (0;11;14)
	Client 10, box 1: (0;11;0)

A feasible 3PP solution considering those boxes is as follows:

  Table 12 gives the whole results with rotations allowed. Results without rotations are available at http://www.isima.fr/~toussain/

	DLT_3LCVRP_1	DLT_3LCVRP_2	DLT_3LCVRP_3 DLT_3LCVRP_4
	1069.24	2522.99	3936.89	5370.56
	3462.26	4949.80	5220.17	5493.75
	1038.20	2457.20	3520.25	4424.02
	3519.08	5148.23	5357.04	5532.92

Table 13 )

 13 . For one instance, GRASPELS found a solution with 16 vehicles while only 15 vehicles are available.

	instances	client
	DLT_3LCVRP_2b	10
	DLT_3LCVRP_12	104
	DLT_3LCVRP_09	205
	DLT_3LCVRP_21	119
	DLT_3LCVRP_30	25
	DLT_3LCVRP_40	116
	DLT_3LCVRP_49	83
	DLT_3LCVRP_50	44

Table 13 : Client packing failure with the heuristic if rotations are not allowed

 13

Appendix 1 instance (Gendreau et al., 2006) (Fuellerer et al., 2010)

							GRASP ELS
	01	297.65	3.40	297.65	1.00	297.65	3.53	297.65	0.0
	02	334.96	0.60	334.96	0.10	335.67	0.06	334.96	0.0
	03	362.27	448.10	362.27	16.20	362.27	13.99	362.27	0.2
	04	430.89	11.10	430.89	0.50	430.88	0.40	430.88	0.0
	05	395.64	0.50	406.50	9.60	379.43	8.16	379.43	0.1
	06	495.85	14.70	495.85	1.20	495.85	0.30	495.85	0.0
	07	742.23	1.80	732.52	18.10	725.43	237.39	725.43	4.9
	08	735.14	104.90	735.14	13.30	735.14	36.63	735.14	1.1
	09	630.13	977.80	630.13	3.70	630.13	2.18	630.13	0.1
	10	717.90	410.70	711.45	92.60	687.57	589.11	687.57	32.1
	11	718.24	208.10	718.25	81.90	718.24	1453.35	718.24	1.8
	12	614.60	1 302.70	612.63	7.50	610.05	19.66	610.00	2.0
	13	2 316.56	2 317.30	2391.77	174.50	2306.04	1242.44	2306.04	86.9
	14	1 276.60	2 121.30	1222.17	425.90	1186.96	2423.64	1184.44	3600.2
	15	1 196.55	2 916.14	1182.86	645.00	1161.20	2144.72	1161.11	689.3
	16	698.61	863.00	698.61	2.80	698.61	2.87	698.61	0.0
	17	906.42	753.20	862.18	3.10	861.80	8.58	861.79	1.2
	18	1 124.33	2198.90	1112.18	1484.60	1084.26	1893.69	1078.41	2030.8
	19	680.29	1 390.30	671.60	414.40	670.44	3322.67	658.34	3429.6
	20	529.00	7 007.50	515.39	1436.70	510.95	2892.97	503.30	1469.7
	21	1 004.40	6 262.50	951.87	2105.70	943.05	4173.74	921.25	4697.4
	22	1 068.96	2 078.70	1030.12	1218.40	1029.87	3561.80	1009.45	3348.3
	23	1 012.51	4 314.10	971.05	1231.70	987.06	3120.66	976.46	1889.1
	24	1 063.61	1 052.50	1057.39	184.70	1056.33	2610.20	1047.75	682.8
	25	1 371.32	500.90	1207.97	3986.10	1232.73	4489.01	1219.77	4658.4
	26	1 557.12	1 075.00	1453.39	2843.60	1415.15	3484.63	1393.76	3066.6
	27	1 378.52	3 983.20	1333.16	2208.30	1317.38	3372.87	1304.82	2422.3
	Avg. Cost	876.31		856.67		847.04		841.96
	Avg. Time		1567.4		689.3		1522.56		1189.44
	Avg. Norm.		1504.1		689.3		1004.89		785.03
	Time

Table A1 : Solution values (rotations allowed)

 A1


Dr Pierre Lefebvre
email: pierre.lefebvre@univ-montp2.fr

S Albert

J Ristić

S Fernández-Garrido

J Grandal

M.-A Sánchez- García

E Calleja

Oxygen photo-adsorption related quenching of photoluminescence in group-III nitride nanocolumns

Keywords: nanocolumns, nanowires, group III nitrides, GaN, InGaN, photoluminescence, surface effects

GaN and InGaN nanocolumns of various compositions are studied by roomtemperature photoluminescence (PL) under different ambient conditions. GaN nanocolumns exhibit a reversible quenching upon exposure to air under constant UV excitation, following a t -1/2 time dependence and resulting in a total reduction of intensity by 85-90%, as compared to PL measured in vacuum, with no spectral change. This effect is not observed when exposing the samples to pure nitrogen. We attribute this effect to photoabsorption and photodesorption of oxygen that modifies the surface potential bending. InGaN nanocolumns, under the same experimental conditions do not show the same quenching features: The high-energy part of the broad PL line is not modified by exposure to air, whereas a lower-energy part, which does quench by 80-90%, can now be distinguished. We discuss the different behaviors in terms of carrier localization and possible composition or strain gradients in the InGaN nanocolumns.

Introduction

In the past decade, semiconductor nanocolumns (NCs) -also called "nanorods" -or nanowires (NWs) based on wide band-gap materials such as ZnO [1], diamond [2] and group-III nitrides have received increasing interest due to their high aspect ratio and compatibility with silicon-based technology, paving the way to a variety of novel applications [3]. In case of nitride NCs, in addition to superior crystal quality, the main reason for this interest is certainly the potential for improving efficiency and reducing costs of light emitting devices in the visible and ultraviolet ranges of the spectrum [4,5]. Indeed, the nearly-perfect crystal quality of NCs is expected to limit the technological steps that are needed, in planar epilayers, for reducing the densities of nonradiative defects [5][6][7][8].

One of the most developed ways to produce vertically aligned, self-assembled III-N NCs is the plasma-assisted molecular beam epitaxy (PA-MBE) [4][5][6][7][8][9][10][11]. This technique [12][13][14][15] and also vapor phase deposition methods [16][17][18] have been used to produce InGaN NCs [13,14,16,18] and nanocolumnar InGaN/GaN quantum heterostructures [12,15,17], providing efficient room-temperature (RT) light emission, probably resulting from carrier localization, as already well-known for InGaN alloys [19]. Emission over the entire visible spectrum was achieved by control of the In incorporation. Now, the high effective surface of dense NC assemblies suggests that surface should play a much more prominent role than in planar structures. As an example of such surface effects, in undoped GaN planar epilayers, Foussekis et al. [20] recently reported that photo-adsorption of oxygen causes significant quenching of the RT PL. This was explained by the upward bending of the near-surface potential, which is increased by the photoinduced adsorption of species in air ambient. This creates a strong electric field that separates electrons and holes, reducing the effective volume for radiative recombination of free excitons, which is the dominant light emission mechanism in high-quality GaN epilayers at RT. The critical effect of photo-absorption on RT PL intensity has been observed in ZnO NWs [1] as well. It has been showed, also, that surface effects can modify excitonic states at low temperatures [21], e.g. via the build-up of surface-excitons [22]. The existence of the latter was also demonstrated in GaN NWs [23], simultaneously with dramatic ageing effects (radical changes in micro-PL spectra, after long exposure to ambient air) resulting from surface alteration. Surface excitons in GaN NWs and NCs (observed at low temperatures, thus in vacuum) were proved to involve donors randomly distributed at variable distance from the surface [23,24]. For completeness, let us mention that surface modifications of transport properties in devices embedding disordered nitride or ZnO NWs, when exposed to certain gases, led to new concepts of gas sensors [25]. In this context, it is important to learn more about possible surface-related effects affecting optical properties of nitride-based NCs.

In this paper, we report on the strong quenching of RT excitonic PL induced by exposing dense ensembles of GaN and InGaN NCs to ambient air, whereas exposure to pure nitrogen has no such effect. We emphasize the significant differences between the two types of NCs, and we discuss the microscopic mechanisms at work.

Samples and experimental details

Our samples were grown by plasma-assisted molecular beam epitaxy on Si(111) substrates. Effective N-rich conditions allow for the self-organized growth of NCs with typical densities, heights and diameters of a few 10 10 cm -2 , 1.0-1.3 µm and 30-50 nm, respectively [11]. We present below the PL results on InGaN NCs grown on top of pure GaN NCs, both parts having similar nominal heights of ~0.6-0.7 µm. The In content was controlled by the growth temperature: for the chosen N, Ga and In fluxes, the growth temperatures were 800°C for GaN and between 650 an d 750°C for InGaN. The lower the growth temperature, the higher the indium incorporation (up to 28% in our samples) [13].

Detailed studies of the morphology and composition of these NCs, of their growth mechanism and of the temperature-dependent PL spectra will be subject of forthcoming publications.

To introduce this work, let us mention that, at T = 8K, all our GaN NCs (including the ones on which InGaN NCs were grown) show an intense donor-bound exciton line at 3.4705 eV, with FWHM of 2.0-2.5 meV, accompanied by free A and B excitons, easily observable.

When the temperature is increased to 300K, upon ionization of the donor-bound excitons, the free excitons dominate the PL spectrum. As for InGaN NCs, their PL features are comparable to those of InGaN epilayers: the linewidth increases with In content, due to enhanced composition fluctuation effects and the PL intensity drops with increasing temperature, but the larger the In composition, the smaller the drop, which can be attributed to increased carrier localization. However, the PL peak energy does not clearly follow the usual s-shaped dependence versus temperature, in the present experimental conditions. Preliminary cathodoluminescence results [26] suggest that the average In content of different NCs can be quite different, thus making complicated the interpretation of T-dependent macroscopic PL spectroscopy.

In the following, we focus on the changes of RT PL when the samples pass from vacuum (pressure of ~10 -4 Pa) to air ambient, under uninterrupted laser excitation. The samples were placed in a chamber where the conditions could be rapidly changed between vacuum and air (or pure nitrogen) at atmospheric pressure and excited by He-Cd laser (λ = 325 nm) with typical power density of 5 W.cm -2 . These conditions are comparable to those of Refs. 1 and 20.

Results

Figure 1 shows the effect of exposure to ambient air on pure GaN NCs. In fact, we observed the same effect in all our GaN NCs, including those followed by InGaN NCs: the intensity of the band-edge excitonic PL quenches by ~85-90 % with no change in spectral shape. The PL intensity drop is fast at the beginning (80 % of the signal is lost after ~30 s), but it slows down progressively with time: a good fit is given by a t -1/2 decay, with a residual signal of ~10% of the initial one, as measured after 2 hours (see the inset in Fig. 1). This quenching dynamics suggests some self-limited chemisorption process in which photocreated free carriers are involved in the adsorption of species at the surface. The accumulation of those species with time raises a potential barrier [1,20] thus hindering further access of free carriers to the surface, which therefore limits additional adsorption.

To confirm the photochemical nature of the process, we checked experimentally that exposure to air alone is not sufficient to induce the quenching effect: starting from vacuum conditions and maximum PL signal, we blocked the laser excitation and exposed the NCs to ambient air for one hour. When the laser was directed again onto the sample (still in air), the PL intensity was not at 10% of the maximum value, as might be expected from the dynamics shown in Fig. 1. Instead, the signal reached 60% of its maximum within a few seconds and, then only, it started decaying, finally reaching 10% of the maximum value. This result is in favor of a photo-adsorption mechanism [20].

The observed effect is reversible: when the sample was placed in vacuum again, under laser excitation, the PL intensity recovered its initial value. This observation suggests the same photodesorption mechanism as the one proposed in Refs. 1 and 20. It also agrees with the recent observation of the increase of low-temperature PL intensity of GaN NCs upon continuous UV excitation, also explained [27] by oxygen photodesorption effects.

To confirm the role played by oxygen, we submitted the samples, under constant laser excitation, to cycles of (i) exposure to industrial grade purity N 2 (O 2 < 0.01 ppm ; H 2 O < 0.02 ppm ; CO/CO 2 < 0.5 ppm) for at least 20 min, followed by (ii) exposure to vacuum for similar times. None of the cycles produced PL intensity changes larger than 5%, contrary to the drastic changes provoked by air. Therefore, oxygen is clearly involved in the quenching mechanism, of which a modeling was proposed in Refs. 1 and 20. From the viewpoint of band structure effects, this model has similarities with the one proposed by Calarco et al. [28] to explain transport measurements in GaN NWs. The model suggests that the observed PL quenching results from the decrease of the volume in which stable excitons can exist, because they are ionized by a strong near-surface electric field. The critical ionization field for excitons in GaN can be roughly estimated as F C = E B / (e.a B ) ~ 10 5 V/cm, where E B and a B are the exciton binding energy (~26 meV) and Bohr radius (~3 nm). Following this model, our experiments indicate that, after exposure to air, the volume of the surface "exciton ionization" layer would correspond to ~90% of the total NC volume.

The observation by Pfüller et al. [27] of an increase of the recombination rate of donorbound excitons at low temperatures, induced by oxygen desorption, encourages us to favor such a model based on the presence of surface electric fields. Nevertheless, we cannot discard completely the possibility that oxygen adsorption creates nonradiative recombination channels near the surface of the GaN NCs.

Let us now discuss the RT PL results obtained in InGaN-on-GaN NC samples, shown in Fig. 2. For small In content, the PL intensity of InGaN NCs was much weaker than the one from the underlying GaN NCs. For this reason, subtraction of a baseline was necessary to analyze the PL of InGaN NCs with In content of ~10% (insert in Fig. 2(a)). For higher In contents (Figs. 2(b) and (c)), the InGaN related PL spectra are more clearly separated from the GaN one. Upon exposure to air, the PL from InGaN NCs showed much less quenching than in the case of GaN. In fact, the higher energy parts of the InGaN PL spectra were nearly unaffected by this exposure to air. On the other hand, our experiments reveal the existence of a lower-energy contribution that is quenched by ~80-90%: the spectrum of InGaN NCs does change. The insert in Fig. 2(a) shows a weak PL contribution at 2.68 eV, for the sample placed in vacuum: after exposure to air, this contribution loses intensity relatively to the higher energy side of the InGaN PL, at 2.83 eV.

For the samples with higher In contents, it is easier to visualize the low-energy contribution: the difference between PL spectra in air and in vacuum is shown by the dashdotted lines in Figs. 2(b) and(c). Clearly, this difference exhibits a peak situated typically 0.15 eV below the PL peak in the vacuum.

Discussion

We need to discuss the differences between GaN and InGaN NCs along two lines: (1) the overall reduced quenching upon exposure to air for InGaN NCs and (2) the observation of a low-energy contribution in InGaN NC PL that does exhibit significant quenching.

(1) We can tentatively attribute the first observation to the localization of carriers at deep potential fluctuations, which exists in the ternary alloys but not in the binary compound [19]. The linewidths of the present PL spectra (several hundred meV) provide an estimation of the overall range of those potential fluctuations. A localized electron-hole pair in a potential well of typically several tens of meV of depth, in InGaN, will be less easily split apart by an electric field of a few hundred kV/cm than the free electron-hole pair in GaN, only bound by Coulomb interaction. Therefore, if a strong surface electric field is established upon oxygen accumulation, localized excitons are more robust than free ones against ionization. They are more robust, too, against capture by nonradiative recombination centers, if this turns out to be the mechanism at work upon oxygen adsorption. Carrier localization can be also invoked as a limiting factor of the oxygen adsorption itself, since the transfer of charged carriers towards the surface was proposed as a crucial mechanism favoring this adsorption. If the carrier mobility towards the surface is hindered by potential fluctuations, the photochemical process of adsorption can be strongly reduced. Finally, we must mention that InGaN is well known to have a stronger residual n-doping than GaN. This may drastically reduce the thickness of the surface space-charge region [28] and therefore minimize the overall sensitivity of the PL in InGaN NCs to oxygen adsorption.

(2) As for the existence of a lower-energy PL contribution that quenches under oxygen exposure, we have to make a preliminary remark: from the above explanation, we can infer that the deeper an exciton is localized, the less it will be affected by surface field effects. Therefore, one would expect the lower-energy part of the PL spectrum to be less sensitive to surface effects, contrary to what we show in Fig. 2. To explain this apparent contradiction, several mechanisms can be invoked. The first explanation supposes the existence, in each of the NCs, of a core and a surface region, the latter being responsible for the lower-energy PL contribution. The simplest reason for this could be In segregation or any type of growth mechanism that would result in In-rich lateral surfaces with PL at lower energy. We do not discard this possibility, though recently reported results rather indicate the presence of a GaN shell around similar InGaN NCs [START_REF] Bougerol | Communication at the International Workshop on Nitride Semiconductors (IWN 2010)[END_REF]. Another way to obtain lower-energy PL close to the surface is to consider a constant average In content but a gradient of strain, for example via strain relaxation mechanisms during growth. The core of the InGaN NC would be compressively strained by lattice-matching on the underlying GaN NC (which is unstrained, as proved by the excitonic transition energies). This compression would result in higher PL energy, for the core region, whereas strain relaxation near the surface would result in lower PL energies. The second explanation is rather related to different average In compositions in different NCs, as seems to be indicated by preliminary cathodoluminescence results. Though not yet fully established, the dynamics of the growth mechanisms could lead, for example, to a dependence of In incorporation on the diameter of the underlying GaN NC. If the PL sensitivity to air exposure is also diameter-dependent, then our results could also arise from that dependence. Why the higher In composition would then correspond to higher sensitivity to oxygen adsorption, still has to be understood, but this possibility should not be neglected at this stage.

Further structural characterization is needed to ascertain any of these hypotheses: insight into local In content and strain, e.g. by high-resolution transmission electron microscopy, may provide the required information.

Conclusions

In summary, we have observed the influence of the atmospheric environment on the RT PL intensities of GaN and InGaN nanocolumns. We have shown that photoadsorption and photodesorption of oxygen induce strong changes in PL intensities. However, while the PL from GaN NCs placed in air ambient is quenched by 90% from its value in vacuum, the PL from InGaN NCs exhibits two clear-cut contributions, of which only the lower-energy one is quenched. The overall PL intensity of InGaN NCs is, however, much less affected by oxygen adsorption than for GaN NCs. We tentatively assign this observation to carrier localization in the ternary alloy. Concerning the observation of two contributions in the PL spectrum of InGaN NCs, we have discussed the results in terms of composition or strain variation within each single NC and amongst the different NCs.

These results might be of importance in terms of applications to light-emitting devices, since we show that surface adsorption of oxygen affects more RT radiative recombination in GaN NCs than in InGaN NCs.

FIGURE CAPTIONS

FIG. 1: Quenching of the room-temperature PL from GaN NCs, after exposure to air ambient under uninterrupted He-Cd laser excitation. The spectrum shape is unchanged during the slow decay, which is well fitted by a power law (in t -1/2 ). Dash-dotted blue curve shows the difference between the black and the red curves and the dotted green line shows the intensity ratio between them (air/vacuum). In (a), logarithmic scale was used and a baseline was subtracted, as shown in inset, in order to extract the weak InGaN PL from the strong GaN PL.

FIG. 2 :

 2 FIG. 2: PL quenching in superposed GaN and InGaN NC samples, with respective InGaN

Acknowledgments

The authors wish to acknowledge partial financial support by the EU under FP7 contract SMASH CP-IP 228999-2, and by national funding agencies under contracts CAM/ P2009/ESP-1503, MICINN/ CSD2006-19 and MICINN/ PLE2009-0023.


Christophe Chesneau

Maher Kachour

Bertrand Maillot

Nonparametric Estimation for Functional Data by Wavelet Thresholding

Keywords: Functional data, Density estimation, Nonparametric regression, Wavelets, Hard thresholding AMS 2000 Subject Classifications: 62G07, 60B11

This paper deals with density and regression estimation problems for functional data. Using wavelet bases for Hilbert spaces of functions, we develop a new adaptive procedure based on wavelet thresholding. We provide theoretical results on its asymptotic performances.

Introduction

Due to technological progress, in particular the enlarged capacity of computer memory and the increasing efficiency of data collection devices, there is a growing number of applied sciences (biometrics, chemometrics, meteorology, medical sciences. . . ) where collected data are curves which require appropriate statistical tools. Because of this, functional data analysis has known a quite important development in the last fifteen years (see e.g. [START_REF] Ramsay | Functional Data Analysis[END_REF], [START_REF] Ramsay | Applied Functional Data Analysis: Methods and Case Studies[END_REF], [START_REF] Ferraty | Nonparametric functional data analysis[END_REF], Dabo-Niang and Ferraty (2008), [START_REF] Ferraty | Special Issue : Statistical Methods and Problems in Infinite-dimensional Spaces. 1st International Workshop on Functional and Operatorial Statistics (IWFOS'2008)[END_REF], [START_REF] Ferraty | The oxford handbook of functional data analysis[END_REF] and [START_REF] Ferraty | Recent Advances in Functional Data Analysis and Related Topics[END_REF] for monographs and collective books on this specific subject). However, whereas there has been substantial work 1 on the nonparametric estimation of the probability density function for univariate and multivariate random variables since the papers of [START_REF] Parzen | On estimation of a probability density function and mode[END_REF] and [START_REF] Rosenblatt | Remarks on some nonparametric estimates of a density function[END_REF], much less attention has been paid to the infinitedimensional case. The extension of the results from the multivariate framework to the infinite dimensional one is not direct since there is no equivalent of the Lebesgue measure on an infinite dimensional Hilbert space. In fact, the only locally finite and translation invariant measure on an infinite dimensional Hilbert space is the null measure and any locally finite measure µ is even very irregular: denoting by B(x, r) the ball of center x and radius r, we have that, for any point x, any arbitrary large M and any arbitrary small r such that µ(B(x, r)) < ∞, there exist (x 1 , x 2 ) ∈ B(x, r) 2 such that µ(B(x 1 , r/4)) < M × µ(B(x 2 , r/4)). For a coverage of the theme of measures on infinite dimension spaces, we refer to [START_REF] Xia | Measure and integration theory on infinite-dimensional spaces[END_REF], [START_REF] Yamasaki | Measures on infinite-dimensional spaces[END_REF], [START_REF] Dalecky | Measures and differential equations in infinite-dimensional space[END_REF] and [START_REF] Uglanov | Integration on infinite-dimensional surfaces and its applications[END_REF].

The first consistency result for a kernel estimator of the density function for infinite dimensional random variables has been obtained in [START_REF] Dabo-Niang | Estimation de la densité dans un espace de dimension infinie: application aux diffusions[END_REF] where a rate is given in the special case when the kernel is an indicator function and the density is defined with respect to the Wiener measure.

Later, different estimators of the density, based on orthogonal series (see [START_REF] Dabo-Niang | Density estimation by orthogonal series in an infinite dimensional space: application to processes of diffusion type I[END_REF]), delta sequences (see Prakasa Rao (2010b)) or wavelets (see Prakasa Rao (2010a)), have been proposed but none of them is adaptive.

Note that the estimation of the density probability function is nonetheless itself of intrinsic interest but it also has a key role in mode estimation and curve clustering (see [START_REF] Dabo-Niang | Mode estimation for functional random variable and its application for curves classification[END_REF]).

Contrary to the chronology of studies in the multivariate case, in the functional framework, estimators of the regression function have been proposed before those of the density. Ferraty and Vieu introduced the first fully nonparametric estimator of the regression function, at first under the hypothesis that the underlying measure has a fractal dimension in [START_REF] Ferraty | Dimension fractale et estimation de la rgression dans des espaces vectoriels semi-norms[END_REF] and then using only probabilities of small balls in [START_REF] Ferraty | Nonparametric models for functional data, with application in regression, time-series prediction and curve discrimination[END_REF]. However, since these pioneering works, no adaptive estimator has been proposed.

Considering the density estimation problem from functional data, Prakasa Rao (2010a) has recently developed a new procedure based on the multiresolution approach on a separable Hilbert space introduced by [START_REF] Goh | Wavelet bases for Hilbert spaces of functions[END_REF].

This procedure belongs to the family of the linear wavelet estimators. As proved in (Prakasa Rao, 2010a, Theorem 3.1), it enjoys powerful asymptotic properties. However, such a linear wavelet estimator has two drawbacks: it is not adaptive (i.e. its performances are deeply associated to the smoothness of the unknown function) and it is not efficient to estimate functions with complex singularities (the sparsity nature of the wavelet decomposition of the unknown function is not captured). For these reasons, (Prakasa Rao, 2010a, Page 2 lines 14-16) states "it would be interesting to investigate the advantage of these wavelet estimators for functional data by using wavelet thresholding suggested by [START_REF] Donoho | Density estimation by wavelet thresholding[END_REF]". This perspective motivates our study.

Adopting the multiresolution approach on a separable Hilbert space H of [START_REF] Goh | Wavelet bases for Hilbert spaces of functions[END_REF], we construct a new adaptive wavelet procedure using the hard thresholding rule of [START_REF] Donoho | Density estimation by wavelet thresholding[END_REF]. Since H remains an abstract space, we propose to evaluate its asymptotic properties over the intersection of two different kinds of Besov spaces (defined in Section 2). The considered spaces are related to the maxiset approach introduced by [START_REF] Cohen | Maximal spaces with given rate of convergence for thresholding algorithms[END_REF].

They are of interest as they contain a wide variety of unknown functions, complex or not. Finally, we adapt the construction of our wavelet hard thresholding estimator to the problem of regression modeling with functional data. Its asymptotic properties are explored.

The paper is structured as follows. In Section 2, we briefly describe the wavelet bases on H and we define some decomposition spaces. The density estimation problem for functional data via wavelet thresholding is considered in Section 3. The regression one is developed in Section 4. The proofs are gathered in Section 5.

Wavelet Bases on H and Decomposition Spaces

Wavelet Bases on H

Let us briefly describe the construction of wavelet bases on H introduced by [START_REF] Goh | Wavelet bases for Hilbert spaces of functions[END_REF]. Let H be a separable Hilbert space of real-or complex-valued functions defined on a complete separable metric space or a normed vector space S. Since H is separable, it has an orthonormal basis E = {e j ; j ∈ Λ} for some countable index set Λ. As usual, we denote by ., . and ||.|| the inner product and corresponding norm that H is equipped with.

Let {I k ; k ≥ 0} be an increasing sequence of finite subsets of Λ such that k≥0 I k = Λ and, for any k ≥ 0, J k = I k+1 /I k . For any k ≥ 0, we suppose that there exist ζ k,ℓ ∈ S, ℓ ∈ I k and η k,ℓ ∈ S, ℓ ∈ J k , such that the two matrices

A k = (e j (ζ k,ℓ )) (j,ℓ)∈I 2 k , B k = (e j (η k,ℓ )) (j,ℓ)∈J 2 k ,
satisfy one of the two following conditions:

(A1) A * k A k = diag(c k,ℓ ) ℓ∈I k and B * k B k = diag(s k,ℓ ) ℓ∈J k , where c k,ℓ , s k,ℓ ′ , for ℓ ∈ I k and ℓ ′ ∈ J k , are positive constants, (A2) A k A * k = diag(d k,j ) j∈I k and B k B * k = diag(t k,j ) j∈J k , where d k,j , t k,j ′ for j ∈ I k and j ′ ∈ J k , are positive constants.
For any x ∈ S, we set

           φ k (x; ζ k,ℓ ) = j∈I k 1 √ g j,k,ℓ e j (ζ k,ℓ )e j (x), ψ k (x; η k,ℓ ) = j∈J k 1 h j,k,ℓ e j (η k,ℓ )e j (x),
where

g j,k,ℓ =    c k,ℓ if (A1), d k,j if (A2), h j,k,ℓ =    s k,ℓ if (A1), t k,j if (A2).
Then the collection

B = {φ 0 (x; ζ 0,ℓ ), ℓ ∈ I 0 ; ψ k (x; η k,ℓ ), k ≥ 0, ℓ ∈ J k }
is an orthonormal basis for H (see [START_REF] Goh | Wavelet bases for Hilbert spaces of functions[END_REF], Theorem 2 (a))).

Consequently, any f ∈ H can be expressed on B as

f (x) = ℓ∈I 0 α 0,ℓ φ 0 (x; ζ 0,ℓ ) + k≥0 ℓ∈J k β k,ℓ ψ k (x; η k,ℓ ), x ∈ S,
where

α 0,ℓ = f, φ 0 (.; ζ 0,ℓ ) , β k,ℓ = f, ψ k (.; η k,ℓ ) . (2.1)
We formulate the two following assumptions on E:

• there exists a constant C 1 > 0 such that, for any integer k ≥ 0,

j∈I k 1 g j,k,ℓ |e j (ζ k,ℓ )| 2 ≤ C 1 , j∈J k 1 h j,k,ℓ |e j (η k,ℓ )| 2 ≤ C 1 . (2.2)
This assumption is obviously satisfied under (A1) with C 1 = 1. Remark also that the second example in [START_REF] Goh | Wavelet bases for Hilbert spaces of functions[END_REF], Section 4) satisfies both (A2) and (2.2).

• there exists a constant C 2 > 0 such that, for any integer k ≥ 0, sup

x∈S j∈J k |e j (x)| 2 ≤ C 2 |J k |. (2.3)
This assumption is satisfied by the three examples in Goh (2007) (we Rao, 2010a, (3.16)).

have sup x∈S sup j∈J k |e j (x)| ≤ 1). Remark that it contains (Prakasa

Decomposition Spaces

Let s > 0 and r > 0. From the wavelet coefficients (2.1) of a function f ∈ H,

we define the Besov spaces B s ∞ (H) by B s ∞ (H) =    f ∈ H; sup m≥0 |J m | 2s k≥m ℓ∈J k |β k,ℓ | 2 < ∞    (2.4)
and the "weak Besov spaces" W r (H) by

W r (H) =    f ∈ H; sup λ>0 λ r k≥0 ℓ∈J k 1I {|βk,ℓ|≥λ} < ∞    , (2.5)
where 1I A is the indicator function on A.

Such kinds of function spaces are extensively used in approximation theory for the study of non linear procedures such as thresholding and greedy algorithms. See e.g. DeVore (1998) and [START_REF] Temlyakov | The best m-term approximation and greedy algorithms[END_REF]. From a statistical point of view, they are connected to the maxiset approach. See e.g. [START_REF] Cohen | Maximal spaces with given rate of convergence for thresholding algorithms[END_REF], [START_REF] Kerkyacharian | Thresholding algorithms, maxisets and well concentrated bases (with discussion and a rejoinder by the authors)[END_REF] and [START_REF] Autin | Point de vue Maxiset en estimation non paramétrique[END_REF].

3 Density Estimation for Functional Data

Problem statement

Let {Ω, F, P } be a probability space and {X i ; i ≥ 1} be i.i.d. random variables defined on {Ω, F, P } and taking values in a complete separable metric space or a Hilbert space S associated with the corresponding Borel σ-algebra B. Let P X be the probability measure induced by X 1 on (S, B).

Suppose that there exists a σ-finite measure ν on the measurable space (S, B)

such that P X is dominated by ν. The Radon-Nikodym theorem ensures the existence of a nonnegative measurable function f such that

P X (B) = B f (x)ν(dx), B ∈ B.
In this context, we aim to estimate f based on n observed functional data

X 1 , . . . , X n .
We suppose that f ∈ H, where H is a separable Hilbert space of realvalued functions defined on S and square integrable with respect to the σ-finite measure ν.

Moreover, we suppose that there exists a known constant

C f > 0 such that sup x∈S f (x) ≤ C f . (3.1)

Estimator

Adopting the notation of Section 2, we define the wavelet hard thresholding estimator f by

f (x) = ℓ∈I 0 α0,ℓ φ 0 (x; ζ 0,ℓ ) + mn k=0 ℓ∈J k βk,ℓ 1I | βk,ℓ |≥κ ln n n ψ k (x; η k,ℓ ), (3.2)
x ∈ S, where

αk,ℓ = 1 n n i=1 φ k (X i ; ζ k,ℓ ), βk,ℓ = 1 n n i=1 ψ k (X i ; η k,ℓ ), (3.3)
κ is a large enough constant and m n is the integer satisfying

1 2 n ln n < |J mn | ≤ n ln n .
The construction of f consists in three steps: firstly, we estimate the unknown wavelet coefficients (2.1) of f by (3.3), secondly, we select only the "greatest" βk,ℓ via a hard thresholding (the "universal threshold" κ(ln n/n) 1/2 is considered) and thirdly we reconstruct the selected elements of the initial wavelet basis. Details on the wavelet hard thresholding estimator for

H = L p ([a, b]
) and the standard nonparametric models can be found in [START_REF] Donoho | Density estimation by wavelet thresholding[END_REF], [START_REF] Delyon | On minimax wavelet estimators[END_REF], [START_REF] Härdle | Wavelet, Approximation and Statistical Applications[END_REF] and [START_REF] Vidakovic | Statistical Modeling by Wavelets[END_REF].

Note that our wavelet hard thresholding procedure is adaptive i.e. it does not depend on the knowledge of the smoothness of f .

Results

Theorem 3.1 below evaluates the performance of f assuming that f belongs to the decomposition spaces described in Subsection 2.2.

Theorem 3.1 Consider the density estimation problem described in Subsection 3.1. Suppose that E satisfies (2.2) and (2.3). Let f be given by (3.2). Suppose that f satisfies (3.1) and, for any θ ∈ (0, 1),

f ∈ B θ/2 ∞ (H) ∩ W 2(1-θ) (H), where B θ/2 ∞ (H) is (2.4) with s = θ/2 and W 2(1-θ) (H) (2.5) with r = 2(1 -θ).
Then there exists a constant C > 0 such that

E(|| f -f || 2 ) ≤ C ln n n θ for n large enough.
An immediate consequence is the following upper bound result:

if f ∈ B s/(2s+1) ∞ (H) ∩ W 2/(2s+1) (H) for s > 0, then there exists a constant C > 0 such that E(|| f -f || 2 ) ≤ C ln n n 2s/(2s+1)
.

This rate of convergence corresponds to the near optimal one in the "standard" minimax setting (see e.g. [START_REF] Härdle | Wavelet, Approximation and Statistical Applications[END_REF]).

Moreover, applying (Kerkyacharian and Picard, 2000, Theorem 3.2), one

can prove that B θ/2 ∞ (H) ∩ W 2(1-θ) (H) is the "maxiset" associated to f at the rate of convergence (ln n/n) θ i.e. lim n→∞ n ln n θ E(|| f -f || 2 ) < ∞ ⇔ f ∈ B θ/2 ∞ (H) ∩ W 2(1-θ) (H).

A Note on Regression Estimation for Functional Data

Let {Ω, F, P } be a probability space and {(X i , Y i ); i ≥ 1} be i.i.d. replica of a couple of random variables (X, Y ) defined on {Ω, F, P }, where Y is real valued and X takes values in a complete separable metric space or a Hilbert space S associated with the corresponding Borel σ-algebra B such that

Y = f (X) + ǫ,
f denotes an unknown regresion function and ǫ is a random variable independent of X with ǫ ∼ N (0, 1). We suppose that f ∈ H where H is a separable Hilbert space of real-valued functions defined on S. Let P X be the probability measure induced by X 1 on (S, B). Suppose that there exists a σ-finite measure ν on the measurable space (S, B) such that P X is dominated by ν. As a consequence of the Radon-Nikodym theorem, there exists a nonnegative measurable function g such that

P X (B) = B g(x)ν(dx), B ∈ B.
We suppose that g is known.

In this context, we want to estimate f from (X 1 , Y 1 ), . . . , (X n , Y n ).

Note that the kernel estimator of the regression function for functional data has been proposed by [START_REF] Ferraty | Nonparametric models for functional data, with application in regression, time-series prediction and curve discrimination[END_REF].

Here, we suppose that there exist two known constants C f > 0 and c g > 0 such that

sup x∈S f (x) ≤ C f , inf x∈S g(x) ≥ c g . (4.1)
Theorem 4.1 Consider the regression estimation problem described above.

Suppose that E satisfies (2.2) and (2.3). Let f be as in (3.2) with

αk,ℓ = 1 n n i=1 Y i g(X i ) φ k (X i ; ζ k,ℓ ), βk,ℓ = 1 n n i=1 Y i g(X i ) ψ k (X i ; η k,ℓ ),
κ is a large enough constant and m n is the integer satisfying 1 2

n (ln n) 2 < |J mn | ≤ n (ln n) 2 .
Suppose that f and g satisfy (4.1) and, for any θ ∈ (0, 1),

f ∈ B θ/2 ∞ (H) ∩ W 2(1-θ) (H), where B θ/2 ∞ (H) is (2.4) with s = θ/2 and W 2(1-θ) (H) (2.5) with r = 2(1 -θ). Then there exists a constant C > 0 such that E(|| f -f || 2 ) ≤ C ln n n θ for n large enough. Again, note that, if f ∈ B s/(2s+1) ∞ (H) ∩ W 2/(2s+1) (H) for s > 0, then there exists a constant C > 0 such that E(|| f -f || 2 ) ≤ C ln n n 2s/(2s+1)
.

This rate of convergence corresponds to the near optimal one in the "standard" minimax setting (see e.g. [START_REF] Härdle | Wavelet, Approximation and Statistical Applications[END_REF]).

Proofs

In this section, C denotes any constant that does not depend on j, k and n.

Its value may change from one term to another and may depends on φ or ψ.

Proof of Theorem 3.1. The proof of Theorem 3.1 is a consequence of (Kerkyacharian and Picard, 2000, Theorem 3.1) with c(n) = (ln n/n) 1/2 , σ i = 1, r = 2 and the following proposition.

Proposition 5.1 For any k ∈ {0, . . . , m n } and any ℓ ∈ I k or ℓ ∈ J k , let α k,ℓ and β k,ℓ be given by (2.1), and αk,ℓ and βk,ℓ be given by (3.3). Then (i) there exists a constant C > 0 such that

E(| αk,ℓ -α k,ℓ | 2 ) ≤ C ln n n .
(ii) there exists a constant C > 0 such that

E(| βk,ℓ -β k,ℓ | 4 ) ≤ C ln n n 2 .
(iii) for κ > 0 large enough, there exists a constant C > 0 such that

P | βk,ℓ -β k,ℓ | ≥ κ 2 ln n n ≤ 2 ln n n 2 .
Let us now prove (i), (ii) and (iii) of Proposition 5.1 (which corresponds to (Kerkyacharian and Picard, 2000, (3.1) and (3.2) of Theorem 3.1)).

(i) We have

E(α k,ℓ ) = E(φ k (X 1 ; ζ k,ℓ )) = S f (x)φ k (x; ζ k,ℓ )ν(dx) = α k,ℓ . (5.1) So E(| αk,ℓ -α k,ℓ | 2 ) = V (α k,ℓ ) = 1 n V (φ k (X 1 ; ζ k,ℓ )) ≤ 1 n E |φ k (X 1 ; ζ k,ℓ )| 2 .
It follows from (3.1), the fact that E is an orthonormal basis of H and (2.2) that

E |φ k (X 1 ; ζ k,ℓ )| 2 = S |φ k (x; ζ k,ℓ )| 2 f (x)ν(dx) ≤ C f S |φ k (x; ζ k,ℓ )| 2 ν(dx) = C f S j∈I k 1 √ g j,k,ℓ e j (ζ k,ℓ )e j (x) 2 ν(dx) = C f j∈I k 1 g j,k,ℓ |e j (ζ k,ℓ )| 2 ≤ C f C 1 . (5.2)
Therefore there exists a constant C > 0 such that

E(| αk,ℓ -α k,ℓ | 2 ) ≤ C 1 n ≤ C ln n n .
(ii) Proceeding as in (5.1), we show that E

(ψ k (X i ; η k,ℓ )) = β k,ℓ . Hence E(| βk,ℓ -β k,ℓ | 4 ) = 1 n 4 E   n i=1 U i,k,ℓ 4   , (5.3)
where

U i,k,ℓ = ψ k (X i ; η k,ℓ ) -E(ψ k (X i ; η k,ℓ )).
We will bound this last term via the Rosenthal inequality (recalled in the Appendix).

We have E(U 1,k,ℓ ) = 0.

By the Hölder inequality and (5.2) with

ψ k (X 1 ; η k,ℓ ) instead of φ k (X 1 ; ζ k,ℓ ),
we have

E(|U 1,k,ℓ | 2 ) ≤ CE |ψ k (X 1 ; η k,ℓ )| 2 ≤ C.
(5.4)

Let us now investigate the bound of E(|U 1,k,ℓ | 4 ). Observe that, thanks to the Cauchy-Schwarz inequality, (2.2) and (2.3), we have

sup x∈S |ψ k (x; η k,ℓ )| ≤ sup x∈S j∈J k 1 h j,k,ℓ |e j (η k,ℓ )||e j (x)| ≤   j∈J k 1 h j,k,ℓ |e j (η k,ℓ )| 2   1/2   sup x∈S j∈J k |e j (x)| 2   1/2 ≤ C 1/2 1 C 1/2 2 |J k | ≤ C |J mn | ≤ C n ln n .
(5.5)

The Hölder inequality, (5.5) and (5.4) yield

E(|U 1,k,ℓ | 4 ) ≤ CE |ψ k (X 1 ; η k,ℓ )| 4 ≤ CnE |ψ k (X 1 ; η k,ℓ )| 2 ≤ Cn. (5.6)
It follows from the Rosenthal inequality, (5.4) and (5.6) that

1 n 4 E   n i=1 U i,k,ℓ 4   ≤ C 1 n 4 max nE |U 1,k,ℓ | 4 , nE |U 1,k,ℓ | 2 2 ≤ C 1 n 2 ≤ C ln n n 2 .
(5.7) By (5.3) and (5.7), we prove the existence of a constant C > 0 such that

E(| βk,ℓ -β k,ℓ | 4 ) ≤ C ln n n 2 .
(iii) We adopt the same notation as in (ii). Observe that

P | βk,ℓ -β k,ℓ | ≥ κ 2 ln n n = P n i=1 U i,k,ℓ ≥ n κ 2 ln n n .
(5.8)

We will bound this probability via the Bernstein inequality (recalled in the Appendix).

We have E(U 1,k,ℓ ) = 0.

By (5.5),

|U 1,k,ℓ | ≤ C sup x∈S |ψ k (x; η k,ℓ )| ≤ C n ln n . Applying (5.2) with ψ k (X 1 ; η k,ℓ ) instead of φ k (X 1 ; ζ k,ℓ ), we obtain E(|U 1,k,ℓ | 2 ) ≤ C.
It follows from the Bernstein inequality that κ) , (5.9) where w(κ) =

P n i=1 U i,k,ℓ ≥ n κ 2 ln n n ≤ 2 exp   - Cn 2 κ 2 ln n n n + nκ ln n n n ln n   ≤ 2n -w(
Cκ 2 1 + κ .

Since lim κ→∞ w(κ) = ∞, combining (5.17) and (5.19), and taking κ such that w(κ) = 2, we have

P | βk,ℓ -β k,ℓ | ≥ κ 2 ln n n ≤ C 1 n 2 ≤ C ln n n 2 .
The points (i), (ii) and (iii) of Proposition 5.1 are proved. The proof of Theorem 3.1 is complete.

Proof of Theorem 4.1. As in the proof of Theorem 3.1, we only need to prove (i), (ii) and (iii) of Proposition 5.1.

(i) Since X 1 and ǫ 1 are independent and E(ǫ 1 ) = 0, we have

E(α k,ℓ ) = E Y 1 g(X 1 ) φ k (X 1 ; ζ k,ℓ ) = E f (X 1 ) g(X 1 ) φ k (X 1 ; ζ k,ℓ ) = S f (x) g(x) φ k (x; ζ k,ℓ )g(x)ν(dx) = α k,ℓ . (5.10) So E(| αk,ℓ -α k,ℓ | 2 ) = V (α k,ℓ ) = 1 n V Y 1 g(X 1 ) φ k (X 1 ; ζ k,ℓ ) ≤ 1 n E Y 1 g(X 1 ) φ k (X 1 ; ζ k,ℓ ) 2 . It follows from (4.1), |Y 1 | ≤ C f + |ǫ 1 |, g(X 1 ) ≥ c g , the independence between X 1 and ǫ 1 , E(ǫ 2 1 ) = 1, the fact that E is an orthonormal basis of H and (2.2) that E Y 1 g(X 1 ) φ k (X 1 ; ζ k,ℓ ) 2 ≤ (C 2 f + 1) 1 c g E |φ k (X 1 ; ζ k,ℓ )| 2 1 g(X 1 ) = (C 2 f + 1) 1 c g S |φ k (x; ζ k,ℓ )| 2 1 g(x) g(x)ν(dx) = C S |φ k (x; ζ k,ℓ )| 2 ν(dx) = C S j∈I k 1 √ g j,k,ℓ e j (ζ k,ℓ )e j (x) 2 ν(dx) = C j∈I k 1 g j,k,ℓ |e j (ζ k,ℓ )| 2 ≤ C.
(5.11)

Therefore there exists a constant C > 0 such that

E(| αk,ℓ -α k,ℓ | 2 ) ≤ C 1 n ≤ C ln n n .
(ii) Proceeding as in (5.10), we show that

E (Y i ψ k (X i ; η k,ℓ )/g(X i )) = β k,ℓ . Set U i,k,ℓ = Y i g(X i ) ψ k (X i ; η k,ℓ ) -E Y i g(X i ) ψ k (X i ; η k,ℓ ) .
and observe that

E(| βk,ℓ -β k,ℓ | 4 ) = 1 n 4 E   n i=1 U i,k,ℓ 4   .
(5.12)

We will bound this last term via the Rosenthal inequality (recalled in the Appendix).

We have E(U 1,k,ℓ ) = 0.

By the Hölder inequality and (5.11) with

ψ k (X 1 ; η k,ℓ ) instead of φ k (X 1 ; ζ k,ℓ ),
we obtain

E(|U 1,k,ℓ | 2 ) ≤ CE Y 1 g(X 1 ) ψ k (X 1 ; η k,ℓ ) 2 ≤ C. ( 5
(x; η k,ℓ )| ≤ sup x∈S j∈J k 1 h j,k,ℓ |e j (η k,ℓ )||e j (x)| ≤   j∈J k 1 h j,k,ℓ |e j (η k,ℓ )| 2   1/2   sup x∈S j∈J k |e j (x)| 2   1/2 ≤ C 1/2 1 C 1/2 2 |J k | ≤ C |J mn | ≤ C n (ln n) 2 . (5.14)
The Hölder inequality, (5.14) and (5.13) yield

E(|U 1,k,ℓ | 4 ) ≤ CE |ψ k (X 1 ; η k,ℓ )| 4 ≤ CnE |ψ k (X 1 ; η k,ℓ )| 2 ≤ Cn.
(5.15)

It follows from the Rosenthal inequality, (5.13) and (5.15) that

1 n 4 E   n i=1 U i,k,ℓ 4   ≤ C 1 n 4 max nE |U 1,k,ℓ | 4 , nE |U 1,k,ℓ | 2 2 ≤ C 1 n 2 ≤ C ln n n 2 .
(5.16) By (5.12) and (5.16), we prove the existence of a constant C > 0 such that

E(| βk,ℓ -β k,ℓ | 4 ) ≤ C ln n n 2 . (iii) We adopt the same notation as in (ii). Since E (Y i ψ k (X i ; η k,ℓ )/g(X i )) = β k,ℓ , we can write U i,k,ℓ = V i,k,ℓ + W i,k,ℓ ,
where

V i,k,ℓ = Y i g(X i ) ψ k (X i ; η k,ℓ )1I A i -E Y i g(X i ) ψ k (X i ; η k,ℓ )1I A i , W i,k,ℓ = Y i g(X i ) ψ k (X i ; η k,ℓ )1I A c i -E Y i g(X i ) ψ k (X i ; η k,ℓ )1I A c i , A i = |ǫ i | ≥ c * √ ln n
and c * denotes a constant which will be chosen later.

We have

P | βk,ℓ -β k,ℓ | ≥ κ 2 ln n n = P n i=1 U i,k,ℓ ≥ n κ 2 ln n n ≤ I 1 + I 2 ,
(5.17)

where

I 1 = P n i=1 V i,k,ℓ ≥ κ 4 √ n ln n
and

I 2 = P n i=1 W i,k,ℓ ≥ κ 4 √ n ln n .
Let us now bound I 1 and I 2 .

Upper bound for I 1 . The Markov inequality and the Cauchy-Schwarz inequality yield

I 1 ≤ 4 κ √ n ln n E n i=1 V i,k,ℓ ≤ C n ln n E(|V 1,k,ℓ |) ≤ C n ln n E(|V 1,k,ℓ | 2 ) ≤ C n ln n E Y 1 g(X 1 ) ψ k (X 1 ; η k,ℓ )1I A 1 2 ≤ C n ln n E Y 1 g(X 1 )
ψ k (X 1 ; η k,ℓ ) 4 1/4

(P (A 1 )) 1/2 .

Using (5.15), an elementary Gaussian inequality and taking c * large enough, we obtain

I 1 ≤ C n ln n n 1/4 e -c 2 * ln n/4 ≤ C 1 n 2 .
(5.18)

Upper bound for I 2 . We will bound this probability via the Bernstein inequality (recalled in the Appendix).

We have E(W 1,k,ℓ ) = 0. Lemma 5.2 [START_REF] Petrov | Limit Theorems of Probability Theory: Sequences of Independent Random Variables[END_REF]) Let n be a positive integer and V 1 , . . . , V n be n i.i.d. zero mean independent random variables such that there exists a constant M > 0 satisfying |V 1 | ≤ M < ∞. Then, for any υ > 0,

P n i=1 V i ≥ υ ≤ 2 exp - υ 2 2 nE(V 2
1 ) + υM /3 .

Using ( 4 .

 4 .1) which implies |Y 1 1I A c 1 | ≤ C f +c * √ ln n ≤ C √ ln n and g(X 1 ) ≥ c g ,and (5.14), we obtain|W i,k,ℓ | ≤ C √ ln n sup x∈S |ψ k (x; η k,ℓ )| ≤ C 11) with ψ k (X 1 ; η k,ℓ ) instead of φ k (X 1 ; ζ k,ℓ ), we obtain E(|W 1,k,ℓ | 2 )Since lim κ→∞ w(κ) = ∞, taking κ such that w(κ) = 2, we haveI 2 ≤ 2 1 n 2 . It follows from (5.17), (5.18) and (5.19) that P | βk,ℓβ k,ℓ | Hence the points (i), (ii) and (iii) of Proposition 5.1 are satisfied by our estimators. The proof of Theorem 4.1 is complete.

Appendix

Here we state the two inequalities that have been used for proving the results in earlier section.

Lemma 5.1 [START_REF] Rosenthal | On the subspaces of L p (p ≥ 2) spanned by sequences of independent random variables[END_REF]) Let n be a positive integer, p ≥ 2 and V 1 , . . . , V n be n zero mean i.i.d. random variables such that E(|V 1 | p ) < ∞.

Then there exists a constant C > 0 such that