Stefania Costache

Nikos Parlavantzas

Christine Morin

Samuel Kortas

Themis: Economy-Based Automatic Resource Scaling for Cloud Systems

High performance computing (HPC) infrastructures are used to execute increasingly complex scientific applications with time-varying resource requirements. A key challenge for infrastructure providers is to distribute resources to such applications so that application performance objectives are met while guaranteeing a high infrastructure utilization. Most of existing solutions provide limited support for meeting application objectives and can lead to inefficient resource usage. In this paper we propose Themis, a system that uses an economic-based approach to automatically allocate resources to the applications that need them the most. Themis relies on a proportional-share auction that ensures fair differentiation between applications while maximizing the resource utilization. To support meeting application performance objectives, Themis provides generic scaling policies based on feedback control loops. These policies adapt the resource demand to the current infrastructure conditions and can be extended for different application types. We have evaluated the performance of the system through simulations, and the results show that Themis can effectively meet application performance objectives while optimizing infrastructure's utilization.

I. INTRODUCTION

High performance computing (HPC) infrastructures are used to execute increasingly complex, dynamic scientific applications, such as scientific workflows and master-worker frameworks. These applications require or can efficiently use time-varying amounts of underlying resources, such as computing nodes, and have to meet various service level objectives (SLOs), such as respecting time constraints. A key challenge for infrastructure providers is to distribute resources to such dynamic applications so that their objectives are met while guaranteeing a high infrastructure utilization. Traditional resource management systems [16] fail to address this challenge. These systems provide a primitive interface for resource control, forcing users to overprovision resources to meet application SLOs, and having a negative impact on utilization.

In contrast to traditional HPC infrastructures, cloud infrastructures support a dynamic provisioning model in which users obtain virtualised resources on-demand and are charged according to a given pricing scheme [1]. This model is attractive for executing dynamic applications for several reasons. First, users can obtain resources "instantly", removing any need for over-provisioning. Second, the pricing scheme gives users the incentive to regulate their resource usage, allowing resources to be allocated to users that need them the most. Third, virtualisation gives providers the flexibility to increase the utilization of hardware resources. However, typical cloud resource management approaches (e.g., Amazon EC2 ondemand instances) suffer from the following limitation. User requests are either fully fulfilled or rejected, forcing the provider to reserve a sufficient amount of idle resources and thus limiting the infrastructure utilization.

An approach that overcomes the previous limitation is using a spot market with a dynamic pricing scheme. In this approach, users place bids for resources, and resource allocations and prices fluctuate depending on the overall demand. When the demand is low, the prices drop and users can request more resources; when the demand is high, prices increase and users are encouraged to request less resources. To apply this approach in managing HPC infrastructures, two issues must be addressed: (i) providing suitable spot market-based mechanisms that distribute virtual resources in a fair and efficient way, and (ii) supporting applications to meet their SLOs while using these mechanisms.

In this paper we present Themis1 , a system that applies a spot market to efficiently manage resources and meet application SLOs. Themis provides: (i) a proportional-share auction based scheduler; (ii) support for building feedback control loops that adapt the application bid and resource demand in order to meet application SLOs. The proportional-share auction ensures a fair use of resources while its fine-grained work-conserving nature maximizes the resource utilization. The feedback control loops are based on generic scaling policies, extended to target specific application types. Policies for two application types are currently provided, namely, for rigid and elastic applications. Simulation experiments show that our system is capable of providing effective support for application SLOs and fairness.

The rest of the paper is organized as follows. Section II introduces the motivations of this work. Section III details the design of our system and Section IV presents the experimental results and discusses the limitations of our solution. Finally, Section V concludes the paper.

II. MOTIVATION AND RELATED WORK

In this section we detail the motivation behind our work and the research related to our goals.

A. Motivation

A general concern in managing HPC infrastructures is how to use all the available capacity to ensure the best application performance, given the workload heterogeneity. Applications executed on these infrastructures can be classified in two general classes: rigid and elastic. A rigid application can use a fixed number of virtual machines during its life-time. An example of a rigid application is an HPC application composed of a fixed set of tasks that need to be executed in parallel. An elastic application can use a dynamic number of virtual machines during its life-time. This number is dependent of the application's goal and the infrastructure resource availability. A master-worker framework can be such an application: a master component submits tasks to worker components that can be added or removed during the application's execution. Both types of applications are highly popular in the scientific community.

Applying the cloud "on-demand" provisioning model on these infrastructures can be beneficial for both users and providers. Elastic applications can scale dynamically and use as many resources as needed at any point. Moreover, as virtualization provides more flexibility to manage resources, fine-grained resource allocation models can be used to maximize the resource utilization. HPC rigid applications can take advantage of these models as they can be executed with fractional amounts of resources, thus avoiding waiting times in queues when there are idle resources. As a result, users can better meet their SLOs while the providers maximize their infrastructure resource utilization. In this context the main challenge that needs to be addressed is what provisioning models to implement and how users can meet their SLOs by using them.

Spot markets can be used as an efficient and fair resource allocation mechanism. With this mechanism resources are "instantly" allocated to users who value them the most while the provider benefits from the mechanism's property to automatically balance the user demand with the infrastructure's resource availability. Spot markets provide dynamic pricing which acts as a generic decentralized allocation mechanism: the resource allocation is not only the task of the resource provider, but it becomes a task distributed between applications/users and the infrastructure. By using this mechanism as a basis for allocating resources, the resource management becomes flexible and generic, allowing more application SLO types to be met.

Our goal is to allocate fine-grained resources in an efficient and fair way by using a spot-market model while providing a generic support for application SLOs.

B. Related Work

We group the related works in two categories: (i) works that focus on providing spot-market mechanisms; (ii) works that focus on providing support for application SLO while acquiring resources from a spot market.

Spot Market Mechanisms: Spot markets are usually implemented through auction mechanisms. In an auction users bid for resources and the provider assigns resources to the users with the highest bids. Recently, Amazon EC2 introduced Spot Instances [7], virtual machines sold through an auction, to sell unused resource capacity and increase the utilization of its infrastructure. To provision spot instances users submit a bid, representing the maximum price they are willing to pay. The spot instances are then available for the user as long as the spot instance price remains below the user's bid. When the price goes above user's bid, the instances are forcibly killed. Many grid resource allocation systems implemented auction protocols. Belagio [START_REF] Auyoung | Resource Allocation in Federated Distributed Computing Infrastructures[END_REF] and Mirage [START_REF] Chun | Mirage: A microeconomic resource allocation system for sensornet testbeds[END_REF] use a combinatorial auction to distribute bundles of distributed resources to users. The same type of auction was proposed in the context of cloud computing to provision classes of virtual machines [START_REF] Zaman | Combinatorial auction-based allocation of virtual machine instances in clouds[END_REF]. These auction types are computationally expensive and to reduce the mechanism complexity users are forced to express their requirements in a coarse-grained manner (e.g., through virtual machine classes).

To overcome this drawback, we use in our work a proportional-share auction to achieve a fine-grained dynamic allocation of resources. This mechanism combines the good properties of the "classic" proportional-share allocation scheme with the efficiency brought by auctions. Through a proportional-share auction, each user submits a bid b i , expressing its payment and receives a share b i /Σb of the resource capacity C. The resource price is given by Σb/C.

Even when the users are price anticipating (i.e., they anticipate the effect of changing their bids on the price), this mechanism achieves a good efficiency [START_REF] Feldman | The proportional-share allocation market for computational resources[END_REF]. Moreover, the mechanism is work-conserving: there are no applications waiting for resources if some fractional amounts are available for use, leading to a maximum utilization of the infrastructure.

Several systems use a proportional-share auction for resource management. In Tycoon [START_REF] Lai | Tycoon: An implementation of a distributed, market-based resource allocation system[END_REF] each host runs a proportional-share auction while agents use a "Best Response" algorithm to select the hosts to which they submit their bids. A dynamic priority scheduler is proposed for Hadoop [START_REF] Sandholm | Dynamic proportional share scheduling in hadoop[END_REF], to allocate map/reduce slots to an user's job. However, these systems do not support SLO management on top of proportionalshare auctions. Libra [START_REF] Yeo | Pricing for utility-driven resource management and allocation in clusters[END_REF] applies a proportional-share scheme on each host to allocate CPU proportional to the application's deadline while a global pricing mechanism is used to balance supply with demand. In contrast with this work, our approach is more flexible as we don't send the application SLO to the resource manager and we allow each agent to adapt itself to the changing resource prices.

SLA-ensurance Mechanisms: Earlier work addressed the problem of providing guidance to the user regarding how much to bid to meet its desired QoS given the resource price volatility [START_REF] Sandholm | A statistical approach to risk mitigation in computational markets[END_REF]. This work was used to provide resource-level QoS guarantees and can be extended with application-level adaptation strategies.

More recent work focused on improving application execution when using spot instances. The main problem that comes from using this model is the lost computation time due to the unpredictable instance termination in case of "out-of-bid" [START_REF] Mattess | Managing peak loads by leasing cloud infrastructure services from a spot market[END_REF]. Our work is complementary. Our spot market is different from spot instances as our resource manager does not forcibly kill virtual machine instances but instead it diminishes the virtual machine's resource share. Our model gives applications more flexibility to react to price changes and we can leverage these approaches to design better policies.

III. DESIGN

Figure 1 gives an overview of our system. Our system relies on the use of a virtual economy. Users receive amounts of currency from a virtual currency manager (i.e., a Virtual Bank), distributed based on predefined policies. To run their applications, users provision virtual clusters for which they purchase amounts of resources from an market-based scheduler that resides on the infrastructure's frontend. To allocate fractional amounts of resources on hosts, the market-based scheduler periodically runs a proportional-share auction. To ensure that each virtual machine receives its paid share and to maximize the resource utilization, the scheduler also performs load balancing. The operations regarding virtual machine management, networking and user management are performed by a Cloud IaaS manager (e.g., OpenNebula [START_REF] Sotomayor | An Open Source Solution for Virtual Infrastructure Management in Private and Hybrid Clouds[END_REF]).

As it is cumbersome for the user to manually scale its virtual cluster whenever its application demands change, each application has an associated application controller. For each application type there is a specific application controller. The application controller monitors the application and uses feedback-based control loops to adapt the application's bid and resource demand to the price variations and user SLO.

A. Virtual Currency Management

To ensure a fair use of resources, users dispose of limited amounts of currency (i.e., credits). When an application needs to be executed, the user creates an account for it and transfers an initial budget of currency. To avoid situations in which this amount is depleted during the application execution, the user can specify a rate at which the application's budget is replenished. However, the accumulated budget cannot exceed the 1: Variables: initial amount given by the user. More complex policies can be applied to prevent hoarding of credits but their investigation is outside the scope of this paper.

B. Resource Allocation and Load-Balancing

We describe next the scheduling algorithm applied by the market-based scheduler. We assume a homogeneous infrastructure and in the following we address only the CPU resource allocation case. The algorithm periodically takes as input the new and existing provisioning requests and outputs the CPU allocation and the physical placement of each virtual machine. The algorithm executes in two phases: (i) it computes the best allocations the virtual machines can get at the user's bid; (ii) it performs load-balancing to ensure the computed allocations within an error bound while considering the virtual machine migration cost.

Resource Allocation: To compute the CPU allocation of each virtual machine we use a proportional-share allocation scheme. However, if we apply the proportional-share policy to allocate to a virtual machine a proportion from the entire cluster, we can reach a situation in which virtual machines cannot be accommodated on the hosts due to their capacity constraints. We illustrate this situation through the following example: we consider 2 hosts each with a capacity of 100 CPU resource units (thus the total cluster capacity is 200 CPU resource units) and 3 virtual machines, each of them having the same bid value b. According to the proportionalshare scheme each virtual machine should receive 66.6 CPU resource units. However, the hosts can only accommodate two virtual machines with this capacity. To accommodate all the virtual machines, their allocations need to be adjusted to fit the host capacity constraints. In our example one host accommodates two virtual machines while the other accommodates one virtual machine. Thus, two virtual machines receive a We obtain this mapping by using a worst-fit decreasing based heuristic. This heuristic sorts the virtual machine requests in descending order by the value of their bid and assigns each request to the host for which the sum of the bids of already assigned requests is minimum. The "ideal" share that each virtual machine should receive is computed by applying the proportional-share scheme on each host.

Load Balancing: As new virtual machine instances need to be created, or the shares for the existing ones change, to ensure that each virtual machine receives its deserved CPU share, some virtual machines need to be migrated from their current hosts to other hosts more fitted for their needs. However, as migration is costly, consuming system resources and affecting the application performance, it is expensive to achieve an optimal mapping of virtual machines to physical hosts (i.e., the mapping that ensures that each virtual machine receives its fair share) at each scheduling interval. There might be cases in which the difference between the actual and ideal allocation for some virtual machines is too small and migrating them to other hosts is useless. Thus, after computing the ideal virtual machine shares, we apply a load-balancing algorithm that tries to reduce the difference between the actual and the ideal allocation (i.e., the allocation error) of each virtual machine while considering the migration cost.

Figure 2 describes the algorithm. The algorithm takes as input the current cluster configuration. Then, at each iteration it tries to accommodate the requirements of the virtual machines that have a negative allocation error (i.e., the ideal allocation is greater than the actual allocation). For this, it first considers the most overloaded hosts. To ensure that the more valuable virtual machines are less migrated, the less valuable virtual machines are the first ones selected for migration. To choose the destination host, a worst-fit heuristic is applied. The solution is improved if the standard deviation over all virtual machine allocation errors is reduced. To memorize the last migration decisions, a tabu list of size 20 is used. The algorithm terminates if the standard deviation becomes smaller than a given threshold or if there is no improvement in the solution for the last 10 iterations. To allow the infrastructure administrator to make a tradeoff between the allocation accuracy and the migration overhead, the algorithm also stops if a maximum number of virtual machine migrations is reached. Then the mapping of virtual machines to physical nodes is sent to the Virtual Infrastructure Manager, indicating which virtual machines should be migrated and created and the CPU share for each of them.

C. Application Control

To ensure that user's SLOs are satisfied, despite resource price dynamicity and application's varying demand, each application is managed by an application controller that scales the resource demand according to application's requirements. To ease the application controller design, we define two classes of scaling policies: (i) vertical and (ii) horizontal. The policies receive as input the performance metric from the application specific sensor and compute the bid and the number of virtual machines for the next time interval. The vertical policies are used to adapt the submitted bid to scale the allocation of a virtual machine. These policies can be used by applications with time-varying resource utilization (e.g., web applications) or rigid applications for which users are willing to make performance-cost tradeoffs (e.g., pay as much as to ensure a specific application progress per time interval). The horizontal policies are used to adapt the number of provisioned virtual machines. These policies can by used by applications that can take advantage of a variable number of virtual machines (e.g., master-worker applications that need to process large amounts of data). For each policy class we define a generic policy that can be extended according to the application's type and SLO.

Figure 3 describes the defined policies and Table I summarizes the information used by them. Both policies can adjust the application bid between a lower (bid min ) and an upper (bid max ) limit. The lower limit is represented by the reserve price of the provider (we assume that there is a minimum price the users are charged). The upper limit is given by the remaining application budget distributed over a time interval decided according to the application type and SLO.

Vertical Scaling: The vertical scaling policy increases the bid if the actual performance metric, v, received from the application is smaller than a reference metric, v ref , and otherwise it decreases it (Lines 1-20, Figure 3). The bid is increased/decreased proportionally to the difference between the actual and reference metric, T. The scaling direction is given by the variable direction. To ensure a faster convergence to the desired value, we combine this policy with a more aggressive increase triggered if the difference between the actual and the reference metric is maintained over a high time period (Lines 16-19, Figure 3). The moment of triggering the aggressive policy is controlled by the variable limit.

Horizontal Scaling: The horizontal scaling policy increases the number of virtual machines if the performance metric of the application is below a given lower limit v ref low and it decreases the number of virtual machines if the performance metric is above an upper limit v ref high (Lines 22-42, Figure 3). As it is difficult to map the difference between the reference and the actual metric to the number of provisioned virtual machines, we scale the virtual machine number based on two thresholds. The use of two thresholds avoids having too many oscilations in the number of virtual machines. Because booting new virtual machines is costly and decisions to reconfigure the application should be taken carefully, the algorithm adjusts the virtual machine number additively. After deciding the desired number of virtual machines, to ensure that the provisioned virtual machines receive a maximum allocation, the application controller adjusts its bid (Line 36, Figure 3). To avoid being charged with a cost higher than the affordable budget, the number of virtual machines is limited at the maximum number the application controller can afford, nvms max .

IV. VALIDATION

To show the efficiency and flexibility of our approach we implemented our framework in CloudSim toolkit [START_REF] Calheiros | Cloudsim: a toolkit for modeling and simulation of cloud computing environments and evaluation of resource provisioning algorithms[END_REF]. We illustrate the use of our framework by considering the case of executing applications under given time constraints (deadlines). We simulate both rigid and elastic applications. Unless otherwise stated, all applications have assigned deadlines. We analyze how our system provides SLO support to applications and how efficient our spot-market model is. To ensure that the policy is usable, we measure the cost in term of virtual machine operations. Finally we show how our approach allows both horizontal and vertical scaling leading to a better resource utilization. We describe next the use of our framework together with the simulation setup and our results.

A. Scaling Examples

We consider the application types defined in Section II, rigid and elastic, and we design two types of application controllers. For the rigid application we design an application controller that estimates the application execution time and by using a vertical scaling policy it keeps the execution time below the user's deadline. The application controller estimates the application execution time by using information

T = | v ref -v v ref | 10: if T < T dif f then 11:
return bid 12:

if step < limit and T < 2 then about the progress the application makes in its computation. This progress can be defined as a number of iterations the application processes per time interval. We assume the application controller has access to application logs to read periodically the remaining number of iterations and to compute the application progress. The user is also required to submit information regarding the total computation amount needed for application to finish its execution and a performance model to associate the progress rate (i.e., progress made per time period) with the CPU allocation. For the elastic application, we design an application controller that estimates the application execution time and uses a horizontal scaling policy to scale the number of workers to keep the execution time below the user's deadline. The application controller communicates with the application master to retrieve information regarding the remaining number of tasks and the task execution time. The control logic behind the policies applied by the application controllers is described in Figure 4.

Scaling a rigid application: To meet its deadline, a rigid application needs to make progress in its computation at a certain rate. For this, the application controller adapts the bids of the virtual machines to keep the progress rate as close as possible to a predefined reference. This reference is computed as the remaining application length (i.e., the total amount of computation) distributed over the remaining time to deadline (Line 3, Figure 4). Based on these metrics, the application controller applies the vertical scaling policy to adapt the virtual machine bids. To avoid spending credits with no benefit, if the application already has a maximum allocation the bid is decreased (Lines 4-5, Figure 4).

When the price is too high to allow continuing the application execution (i.e., the bid the application controller is required to submit for the next time period is bigger than b max ) at the current user's budget, the controller also applies a suspend/resume policy. This policy is as following: if the current and predicted progress become smaller than the reference progress the virtual machine instances are suspended and are resumed when the price drops and enough time has passed.

If the current price is too high when the application is submitted, the controller postpones the application execution for a time interval computed as a part of the time that remains between the deadline and the ideal application execution time. If the execution is postponed too much (the remaining time to deadline becomes smaller than the ideal application execution time) the application controller aborts the application. Otherwise, the application controller starts the application by submitting an initial bid computed to get enough allocation to meet the application reference progress at the estimated resource price.

Scaling a master-worker framework: To meet its deadline, the elastic application needs to have enough workers to finish its tasks in time. To estimate the application execution time, we use the result from [START_REF] Verma | Aria : Automatic resource inference and allocation for mapreduce environments[END_REF] which establishes the performance bounds for the application execution, given the set of workers and the remaining tasks with their average and maximum execution times (Line 11, Figure 4). Based on these metrics, the application controller applies the horizontal scaling policy to adapt the number of virtual machines/workers. If the estimated execution time is greater than the remaining time to deadline the application controller provisions more virtual machines, and otherwise releases them (Line 13, Figure 4). These policies can be easily adapted to other SLOs. For example, if an user wants to execute its application as fast as possible, is enough to replace the progress metric with the allocation that the application receives. Thus, the user ensures that its application gets a maximum allocation possible under its given budget constraint.

B. System Modeling

When simulating our system we pay attention to two issues: (i) simulating the virtual machine operations for the provider; (ii) generating realistic budgets for applications. For the infrastructure provider, we pay attention to the migration operations, as they have an impact on the application performance. We don't model the performance interference or other overheads, as they would appear for other classic algorithms too (e.g., FCFS or other algorithms employed by current cloud resource managers). We model the migration performance overhead as 10% of CPU capacity used by the virtual machines in which the application is running. This is a high estimation, as previous work found that migration incurs an overhead of approximately 8% on the performance of HPC applications [START_REF] Nagarajan | Proactive fault tolerance for hpc with xen virtualization[END_REF]. We estimate the time to migrate the virtual machines as the time to transfer the memory state of the application by using the available network bandwidth.

We generate budgets for application execution requests based on the deadline values. Each application receives a budget computed as B = B base /deadline, where B base is the same initial budget assigned to all applications. This is a realistic way to model the fact that users with stringent deadlines value more their applications and are also more inclined to give a higher budget to them. In all cases the scheduling interval is set to 5 minutes.

C. SLO Satisfaction and Efficiency

To measure the total SLO satisfaction that our mechanisms can provide, we simulate traces of real workloads, to which we assign synthetic budgets and deadlines. We use two usercentric metrics: the percentage of met SLOs (deadlines) and the efficiency of our mechanism. The last metric gives us the level of satisfaction the users achieve from using our spotmarket system, as compared to non-market based systems. To use this metric, we define the value the application brings to the user as the assigned budget if the application finishes before its deadline, and zero otherwise. We consider that users assign to their applications budgets that reflect their true importance: a user with an urgent application assigns a high budget, and thus will gain more value from executing its application, while a user with a best-effort application assigns a low budget. The efficiency of our mechanism is then computed as the ratio between its total provided value and the total value provided by a non-market based mechanism.

Specifically, we compare our mechanism with two wellknown, non-market based policies: First-Come-First-Served and EDF (earliest deadline first). We consider the first model as best-effort provisioning, while the second model is QoSaware. The first model is commonly applied by cloud resource managers. The second model aims to minimize the number of missed deadlines. We use the second model only for comparison reasons. It is not practical for cloud managers to apply this model. If they did, they would support only a single type of application goal (i.e., meeting deadlines).

We simulate a real-world workload of batch applications, taken from an online archive [2] on an infrastructure of 240 hosts. This archive contains public logs of HPC resource management systems. These logs have information regarding the resource demand and execution time of applications submitted to the infrastructure. We used the HPC2N log, as it has detailed information about application memory requirements, needed to estimate the migration overhead. From this trace we selected the first 1000 applications. The applications are assigned synthetic deadlines between 1.2 and 10 times their ideal execution times. The assigned budget is B base = 60 credits/resource unit/scheduling period, renewed at a rate of 20 credits/resource unit/scheduling period. To simulate different utilization levels, we vary the application inter-arrival time (i.e., the time between two consecutive submissions), by multiplying it with a factor (i.e., inter-arrival factor).

Figure 5(a) shows the number of applications that meet their SLOs. In lightly-loaded periods, our mechanism outperforms both EDF and FCFS. When the system becomes overloaded, the mechanism continues to outperform FCFS, but not EDF. Similar results are achieved for the efficiency metric, as illustrated in Figure 5(b). Our mechanism is 2 times more efficient than FCFS in high load periods while it encounters an efficiency loss of up to 25% compared to EDF.

The performance gap between our mechanism and EDF can be explained as follows. The nature of our mechanism is decentralized: there is no centralized entity that decides the allocation of every application to meet a system-wide objective (e.g., minimizing missed deadlines). Rather, each controller acts selfishly and independently from each other to meet its own application SLO. Thus, the solution is less optimal than a centralized solution. However, the main advantage of our solution comes exactly from this decentralization: as each controller is independent, we can support arbitrary application types and SLOs compared to a centralized mechanism that strives to accomplish a fixed global objective. This experiment shows that our mechanism provides good SLO satisfaction to applications and achieves a good efficiency. The performance degradation is the "price" payed by the decentralized nature of our system, that allows each user to behave in a strategic selfish way.

D. Policy Cost

With the same settings from the previous experiment, we measure the cost of our mechanism in term of virtual machine operations. We log the number of migrations and suspend/resume operations for each time interval during the experiment run. Figure 5(c) shows the average number of operations per time interval for each experiment run. Given the large scale of the infrastructure these numbers are reasonable: a maximum of 61 preempted virtual machines and 45 migrations when the system becomes highly loaded. This proves that our system can be used in a real environment, as the virtual machine operations would not have a highly negative impact on the application performance.

E. Flexibility

To prove that our system is not build specifically for one application type, we show the good behavior of our horizontal scaling policies also. For this, we consider an elastic masterworker application that applies the horizontal policy to manage a burst in its workload and finish its execution before the deadline given by the user. The application has 1075 tasks to process in an time interval of 643.25 minutes, with each task having an uniform execution time between 1 and 40 minutes. After the first 85 minutes of execution, an additional workload of 1500 tasks is generated.

Figure 6 shows the number of virtual machines provisioned by our policy in time (expressed, for simplicity, in scheduling intervals). After starting its execution the application controller provisions aproximately 44 virtual machines, as this number is enough for it to finish its execution in a timely manner. When the additional workload is submitted, the application controller starts provisioning more virtual machines to compensate for the increase in demand. The number is increased first incrementally, and then more and more aggressively. In this case, the controller overestimates the number of provisioned resources, so it decreases it afterwards, reaching an equillibrum point at 116 virtual machines. When the application is close to finishing its execution, less tasks remain to be processed and so the controller shuts down the idle virtual machines. This experiment shows that our mechanism allows applications to react reasonably fast to changes in their demand to meet their SLOs. We consider that for HPC workloads these changes don't need to be very fast, so simple policies would suffice.

V. CONCLUSION

In this paper we proposed Themis, a system that shares the resources among applications using a spot market while providing support for application SLOs. Our system is motivated by the current cloud provisioning models, which provide "on-demand" virtualized resources. In contrast with classic resource management systems, clouds provide more flexibility to users in meeting their SLOs. However the provider has to support the burden of managing resources efficiently. To address this problem, Themis brings the following contributions: (i) an implementation of a proportional-share auction that maximizes resource utilization while considering virtual machine migration costs; (ii) a set of feedback-based control policies used to adapt the application bid and resource demand to fluctuations in price.

Our experiments have shown how Themis is capable of providing overall good SLO support. The decentralized nature of markets makes our system flexible to use by different application types. We have shown how applications with dynamic changes in their demand can make use of our system, thus demonstrating that Themis is not designed for only one application type but it provides generic interfaces that can be used by any application.

Currently we tested Themis through simulations. We also plan to evaluate Themis on a real testbed. As next steps, we plan to improve our resource allocation mechanism by support-ing multiple resource types. We also plan to extend our policies for more application types and SLOs. For example, adaptive applications could scale their resource demands to perform computation at a certain progress rate. Another concern that we plan to address is how to allow users to specify placement preferences for their applications. Such preferences are useful for improving application performance or for fault tolerance.

Fig. 1 :

 1 Fig. 1: System overview.

  Fig. 2: VM placement algorithm.

1 : 2 :

 12 Scale-One-VM(v, v ref ) if v < v ref and direction = 1 then if v > v ref and direction = -1 then

Fig. 3 :

 3 Fig. 3: Generic scaling policies.

else 7 : 12 :Fig. 4 :

 7124 Fig. 4: Specific derived policies.

Fig. 5 :

 5 Fig. 5: Application execution statistics.

Fig. 6 :

 6 Fig. 6: Application elastic scaling.

TABLE I :

 I Information used by the generic scaling policies. share of 50 CPU resource units. As we want to keep the proportional-share work-conserving property, the third virtual machine receives 100 CPU resource units.

This work is supported by ANRT through the CIFRE sponsorship No. 0332/2010


Effect of an electric field on a Leidenfrost droplet Franck Celestini * and Geoffroy Kirstetter Laboratoire de Physique de la Matière Condensée, UMR 7366, CNRS, Université de Nice Sophia-Antipolis, Parc Valrose 06108 Nice Cedex 2, France (Dated: May 16, 2012)

We experimentally investigate the effect of an electric field applied between a Leidenfrost droplet and the heated substrate on which it is levitating. We quantify the electro-Leidenfrost effect by imaging the interference fringes between the liquid-vapour and vapour-substrate interfaces. The increase of the voltage induces a decrease of the vapour layer thickness. Above a certain critical voltage the Leidenfrost effect is suppressed and the drop starts boiling. Our study characterizes this way to control and/or to avoid the Leidenfrost effect that is undesirable in many domains such as metallurgy or nuclear reactor safety.

The first observation of the Leidenfrost effect has been reported more than two centuries ago [1] by the German scientist J. G. Leidenfrost. He showed that a drop deposited on a sufficiently hot plate keeps its spheroidal shape and levitates on its own vapour. Because of the weak heat conduction of the vapour as compared to the one of the liquid, the evaporation time is considerably increased and reaches a maximum at T f , the so-called Leidenfrost temperature. For example, a millimetric water drop has a lifetime of several minutes on a substrate kept at T = 250 • C while it evaporates in less than one second when deposited on the same substrate but kept at 150 • C. This subject is still motivating numerous researches. From a fundamental point of view, studies are essentially motivated by the nonwetting situation, considered as the limit of a perfect super-hydrophobicity, and by the high mobility of the droplet due to the weak friction with the substrate on which it is levitating [2][3][4][5]. For many industrial applications the Leidenfrost effect is rather undesirable. As stressed above, the vapour layer in between the liquid and the substrate dramatically reduces the heat conduction. This has important consequences for example in metallurgy for controlling the quenching process of alloys [6] and in the safety of nuclear reactors [7].

In this letter, we study the effect of a tension applied between a Leidenfrost drop and the substrate on which it is levitating. The vapour layer thickness decreases with the applied voltage and a rather low critical voltage (typically 40V for a millimetric drop) makes it possible to suppress the Leidenfrost effect. In spite of numerous studies on the electro-hydrodynamic enhancement of heat transfer in fluid flow [8], only one study [9] has considered this possibility but quantitative studies have not been performed. In this communication, first we present the experimental set-up consisting in the observation of the interferences fringes between the liquid and the substrate. This method allows a precise measurement of the vapour layer thickness. Without any applied voltage we put in evidence an asymmetric underneath profile of the Leidenfrost drop. We then present our results on the relative variation of the vapour layer thickness sustaining the Leidenfrost droplet as a function of the applied voltage. A model is proposed and predicts that the vapour thickness decreases linearly with the square of the applied voltage. This model successfully reproduces the experimental measures performed on water droplets. We finally draw some conclusions on the electo-Leidenfrost phenomenon presented in this letter.

The experimental set-up is presented in Fig. 1. A copper block is kept at a controlled temperature T . A sapphire substrate is inserted in the copper block. To ensure its electrical conductivity, the substrate is coated with a thin (60 nm width) gold film. A water droplet of radius R is generated from a nozzle where a thin tungsten wire (150µm radius) is fixed in order to ensure the electrical conduction with the drop. A low frequency (0.5 Hz) voltage (amplitude ranging between 0 and 30V ) is applied between the drop and the upper coated part of the sapphire. The sapphire is illuminated by the bottom with a white light source. A semi-transparent mirror together with an interferometer filter (λ = 633 nm, width = 30nm) are used to visualized the interference fringes between the liquid-vapour and vapour-sapphire interfaces. A high speed camera is used to record the spatial and time evolution of the interference fringes. The height difference between two white (or black) fringes is of λ/2 = 0.317µm. During the redaction of this letter we have been aware that a similar experimental set-up was used to study the fluctuation of the vapour thickness of a Leidenfrost drop [10].

The interference fringes are present in the nearest contact zone of the droplet with the substrate where the vapour film is thin enough. We assume that the drop has the shape of a sphere flattened on its bottom. The disc between the drop and the film vapour (the so-called Laplacian disc) has a radius l related to the radius R of the drop through [11] :

l ∝ R 2 /κ -1 (1)
where κ -1 = γ/ρ l g is the capillary length, γ the surface tension, g the acceleration of gravity and ρ l the mass liquid density. This relation is obtained assuming R κ -1 and by balancing the gravity force with the Laplace pressure within the drop.

We present in Fig. 2 the interference patterns observed at the bottom of a Leidenfrost droplet of radius R 1mm, contact disc radius l 0.3mm and a temperature sapphire substrate kept at 280 • C. The pictures are taken for V = 0, 15, 25 and 30V . We can observe that the fringes are not centred on the drop center. The height profile of the vapour pocket below the drop is therefore not axis-symmetric. Previous studies [12][13][14] have both theoretically and numerically investigated the equations governing the vapour thickness underneath a Leidenfost droplet. Nevertheless, all studies have supposed an axis-symmetric profile. The asymmetric shape has been systematically observed and is more important when the radius of the droplet decreases. As can be observed on Fig. 2, the electric field tends to reduce the asymmetry and to increase the value of the Laplacian disk radius.

As sketched in the inset of Fig. 3 we define an x diameter axis. At the back of the drop (x = x b ) the height presents a local minimum h = h b , on its the center (x = x c ) the height is maximum at h = h c and on the front of the droplet (x = x c ) the height reaches a minimum at h = h f . The height profile along the x axis is measured and represented in Fig. 3 for a droplet at V = 0. We take h b = h(x b = 0) = 0 as the reference for our relative height interferential measurements. The height difference between the center and the front is δh cf 8µm while the one between the center and the back δh cb 2.5µm. The asymmetry discussed above can be estimated by the the ratio δh cb /δh cf 3.2 and also by (x c -x b )/(x f -x c ) 1.5. Both values are significantly different from the one expected for the axis-symmetric case considered in the previous theoretical studies cited just above.

To evaluate the effect of the applied voltage on the height profile we apply a modulated tension amplitude of the form :

V (t) = V max 2 (1 + cos(2πf (t -t 0 ))). (2)
A low frequency modulation is used f = 0.5 Hz and a maximum voltage V max = 30V . We recorded the evolution of the interferences fringes and present in the supplementary material a movie recorded during several periods of the applied voltage. Using an image analysis, we compute the time evolution of the relative heights h b ,h c and h f . The reference is taken for h c (t = 0) = 0. We represent the relative heights as a function of time in Fig. 4. The dotted line indicates the time t 0 at which the applied voltage is maximum. At t 0 both h f and h b reach a minimum value. The applied voltage induces a capacitive force that attracts the front and the bottom of the droplet toward the substrate. More surprisingly the center of the drop reaches a maximum height when the applied voltage is maximum. We qualitatively explain it considering that the decrease of h b and h f tends to increase the viscous pressure underneath The height reference is taken for hc(t = 0) = 0, the applied voltage has an amplitude of 30V and a frequency of 0.5 Hz.

V = 30 V V = 0 V V = 25 V 0.4 mm l V = 15 V
the drop and therefore to increase height on its center, h c . It is out the scope of this paper to propose a model describing the evolution of the overall droplet profile under the action of the applied voltage. We limit our analysis considering a droplet flattened at its bottom on a disc of radius l and a vapour width h. The relation between R and l without any applied voltage is given by the equation ( 1). For a Leidenfrost droplet of radius R the vapour thickness (without any applied voltage) h 0 can be estimated balancing the gravity force with the viscous pressure of the Poiseuille vapour flow between the drop and the substrate :

F p ∝ ηλ∆T l 4 Lρ v h 4 ,
ρ v being the mass vapour density. It gives an expression for the equilibrium width h 0 satisfying :

ρgR 3 = ηλ∆T l 4 Lρ v h 4 0 (3)
In the previous expression, L is the latent heat, η the vapour viscosity, λ the thermal conductivity and ∆T the temperature difference between the substrate kept at temperature T and the water drop at temperature

T v = 100 • C.
It is important to note that we used the lubrication approximation to obtain the equilibrium thickness h 0 . As can be seen in Fig. 3 and as discussed recently [13] there exist regions with high curvatures below the drop where the lubrication could not be fully justified.

The force, F v due to the applied voltage between the liquid-vapour and vapour-substrate interface can be approximated as the one in between the two armatures of a capacitor. It reads :

F v ∝ -V 2 l 2 /h 2
where is the dielectric constant. We consider the case for which the applied voltage induces a small variation dh of the vapour thickness compared to its equilibrium value h 0 . Using equation ( 3), the balance between the three different forces therefore leads to :

dh ∝ - V 2 4 ρ v L ηλ∆T l 2 h 3 0 (4)
It is worth mentioning that, as for the electro-wetting effect [15], the deviation from the equilibrium situation for the electro-Leidenfrost effect is proportional to the square of the applied voltage.

We use equation ( 2) to eliminate time in the data presented in Fig. 4. We can therefore plot the thickness variation dh as a function of the applied voltage V . This is done in Fig. 4 for the water droplet standing on a substrate kept at 4). Inset : relative variation of the thickness as a function of time for an applied voltage modulated in amplitude with a frequency f = 0.5 Hz. The black line is the best fit to the equation ( 4)

280 • C.
In order to compare the experimental data with our simple model we assume in what follows that the mean thickness variation dh ∝ dh f . The full black line is a best fit to equation ( 4) with a single free parameter h 0 , the equilibrium thickness. The best agreement is found for a value h 0 = 18µm which is close to previous measurements of the vapour thickness for millimetric drops [2]. In the inset of Fig. 4, we plot dh as a function of time and the best fit to our model. Dividing V max by h 0 gives an estimate of the maximum applied electric field E max 1.6 10 6 V m -1 . This value is below the disruptive electric field and therefore justifies that we did not observe any electrical discharges during our experiments. For the data obtained on substrates kept lower temperatures we observed that the model underestimates dh*. At lower temperatures the equilibrium thickness is thinner and the condition dh h 0 is no longer valid for large values of the applied voltage.

Applying a larger voltage suppress the Leidenfrost effect. For the millimetric drops considered in this study, an applied voltage of 40V is sufficient to decrease to 0 the value of h f and to let the drop starts boiling. We illustrate this effect by suddenly applying a voltage of 40V on a millimetric droplet initially standing in the Leidenfrost state on a copper substrate. We represent in Fig. 6 three successive pictures taken at t = 0, 0.1 and 0.4ms. The voltage is applied at t = 0 and we can observe the boiling crisis of the droplet. The movie M2 in the supplementary material displays the evolution of the droplet recorded at 90000 frame per second. The critical voltage for suppressing the Leidenfrost effect could be simply evaluated combining equations ( 3) and (1) with the condition dh/h 0 = 1. Nevertheless this would underestimate the critical voltage. First, equation ( 3) is obtained assuming dh/h 0 1 and, as can be seen in 2 and in the movie M1 of the supplementary material, the applied voltage also tends to increase the value of l. This effect has not be taken into account in our simple model.

In this communication we have investigated the effect of an electric field applied between a Leidenfrost drop and the substrate on which it is standing. Using an interferential technique we have measured the thickness of the vapour layer as a function of the applied tension. The experimental data are in good agreement with the simple model proposed that predicts a decrease of the vapour layer thickness proportional to the square of the applied voltage. For a millimetric drop, applying a voltage of order 40V permits to suppress the Leidenfrost effect and to let the drop starts boiling. This study should find applications in domains where the Leidenfrost effect has to be controlled and/or avoided. The electro-Leidenfrost effect presented in this letter will have to be fully characterised as it has been done for the electro-wetting effect. The electrical conductivity of water as well as the nature of the substrate on which it is standing should have important consequences on its efficiency. Finally, we hope this communication will motivate theoretical and numerical studies concerning the asymmetric underneath shape put in evidence for Leidenfrost droplets without any applied voltage.

We dedicate this paper to the memory of Prof. Richard Kofman who initiated this study.

FIG. 1 :

 1 FIG.1: Experimental set-up.

FIG. 2 :FIG. 3 :

 23 FIG. 2: Interferences patterns at the bottom of a Leidenfrost droplet of radius R 1 mm and Laplacian disk radius l. The interference patterns are represented for V = 0, 15, 25 and 30V .

FIG. 4 :

 4 FIG.4: Time dependence of the relative heights, dhc (red squares), dh b (blue diamonds) and dh f (black circles) for a Leidenfrost droplet deposited on a substrate at T = 280 • C. The dotted line indicates the time t0 at which the applied voltage is maximum. The height reference is taken for hc(t = 0) = 0, the applied voltage has an amplitude of 30V and a frequency of 0.5 Hz.

FIG. 5 :

 5 FIG.5: Thickness variation dh as a function of the applied voltage V for a Leidenfrost droplet standing on a substrate kept at 280 • C. The full line corresponds to the best fit to the equation (4). Inset : relative variation of the thickness as a function of time for an applied voltage modulated in amplitude with a frequency f = 0.5 Hz. The black line is the best fit to the equation(4)

[1] J. G. Leidenfrost, De Aquae Communis Nonnullis Qualitatibus Tractatus (Duisbourg, 1756)


Gérard Biau
email: gerard.biau@upmc.fr

Luc Devroye
email: lucdevroye@gmail.com

Adam Krzyżak
email: krzyzak@cs.concordia.ca

An Affine Invariant k-Nearest Neighbor Regression Estimate

Keywords: Nonparametric estimation, Regression function estimation, Affine invariance, Nearest neighbor methods, Mathematical statistics. 2010 Mathematics Subject Classification: 62G08, 62G05, 62G20

We design a data-dependent metric in R d and use it to define the knearest neighbors of a given point. Our metric is invariant under all affine transformations. We show that, with this metric, the standard k-nearest neighbor regression estimate is asymptotically consistent under the usual conditions on k, and minimal requirements on the input data.

Introduction

The prediction error of standard nonparametric regression methods may be critically affected by a linear transformation of the coordinate axes. It is typically the case for the popular k-nearest neighbor (k-NN) predictor (Fix and Hodges [START_REF] Fix | Discriminatory analysis. Nonparametric discrimination: Consistency properties[END_REF][START_REF] Fix | Discriminatory analysis: Small sample performance[END_REF], Cover and Hart [START_REF] Cover | Nearest neighbor pattern classification[END_REF], Cover [START_REF] Cover | Estimation by the nearest neighbor rule[END_REF][START_REF] Cover | Rates of convergence for nearest neighbor procedures[END_REF]), where a mere rescaling of the coordinate axes has a serious impact on the capabilities of this estimate. This is clearly an undesirable feature, especially in applications where the data measurements represent physically different quantities, such as temperature, blood pressure, cholesterol level, and the age of the patient. In this example, a simple change in, say, the unit measure of the temperature parameter will lead to totally different results, and will thus force the statistician to use a somewhat arbitrary preprocessing step prior to the k-NN estimation process. Furthermore, in several practical implementations, one would like, for physical or economical reasons, to supply the freshly collected data to some machine without preprocessing.

In this paper, we discuss a variation of the k-NN regression estimate whose definition is not affected by affine transformations of the coordinate axes. Such a modification could save the user a subjective preprocessing step and would save the manufacturer the trouble of adding input specifications.

The data set we have collected can be regarded as a collection of independent and identically distributed R d × R-valued random variables D n = {(X 1 , Y 1 ), . . . , (X n , Y n )}, independent of and with the same distribution as a generic pair (X, Y ) satisfying E|Y | < ∞. The space R d is equipped with the standard Euclidean norm . . For fixed x ∈ R d , our goal is to estimate the regression function r(x) = E[Y |X = x] using the data D n . In this context, the usual k-NN regression estimate takes the form

r n (x; D n ) = 1 k n kn i=1 Y (i) (x),
where (X (1) (x), Y (1) (x)), . . . , (X (n) (x), Y (n) (x)) is a reordering of the data according to increasing distances X ix of the X i 's to x. (If distance ties occur, a tie-breaking strategy must be defined. For example, if X ix = X jx , X i may be declared "closer" if i < j, i.e., the tie-breaking is done by indices.) For simplicity, we will suppress D n in the notation and write r n (x) instead of r n (x; D n ). Stone [START_REF] Stone | Consistent nonparametric regression[END_REF] showed that, for all p ≥ 1, E[r n (X)-r(X)] p → 0 for all possible distributions of (X, Y ) with E|Y | p < ∞, whenever k n → ∞ and k n /n → 0 as n → ∞. Thus, the k-NN estimate behaves asymptotically well, without exceptions. This property is called L p universal consistency.

Clearly, any affine transformation of the coordinate axes influences the k-NN estimate through the norm . , thereby illuminating an unpleasant face of the procedure. To illustrate this remark, assume that a nontrivial affine transformation T : z → Az+b (that is, a nonsingular linear transformation A followed by a translation b) is applied to both x and X 1 , . . . , X n . Examples include any number of combinations of rotations, translations, and linear rescalings. Denote by

D ′ n = (T (X 1 ), Y 1 ), . . . , (T (X n ), Y n ) the transformed sample. Then, for such a function T , one has r n (x; D n ) = r n (T (x); D ′ n ) in general, whereas r(X) = E[Y |T (X)]
since T is bijective. Thus, to continue our discussion, we are looking in essence for a regression estimate r n with the following property:

r n (x; D n ) = r n (T (x); D ′ n ). (1.1)
We call r n affine invariant. Affine invariance is indeed a very strong but highly desirable property. In R d , in the context of k-NN estimates, it suffices to be able to define an affine invariant distance measure, which is necessarily data-dependent. With this objective in mind, we develop in the next section an estimation procedure featuring (1.1) which in form coincides with the k-NN estimate, and establish its consistency in Section 3. Proofs of the most technical results are gathered in Section 4.

It should be stressed that what we are after in this article is an estimate of r which is invariant by an affine transformation of both the query point x and the original regressors X 1 , . . . , X n . When the sole regressors are subject to such a transformation, it is then more natural to talk of "affine equivariant" regression estimates rather than of "affine invariant" ones; this is more in line with the terminology used, for example, in Ollila, Hettmansperger, and Oja [START_REF] Ollila | Estimates of regression coefficients based on sign covariance matrix[END_REF] and Ollila, Oja, and Koivunen [START_REF] Ollila | Estimates of regression coefficients based on lift rank covariance matrix[END_REF]. These affine invariance and affine equivariance requirements, however, are strictly equivalent.

There have been many attempts in the nonparametric literature to achieve affine invariance. One of the most natural ones relates to the so-called transformation-retransformation proposed by Chakraborty, Chaudhuri, and Oja [START_REF] Chakraborty | Operating transformation retransformation on spatial median and angle test[END_REF]. That method and many variants have been discussed in texts such as [START_REF] Devroye | A Probabilistic Theory of Pattern Recognition[END_REF] and [START_REF] Györfi | A Distribution-Free Theory of Nonparametric Regression[END_REF] for pattern recognition and regression, respectively, but they have also been used in kernel density estimation (see, e.g., Samanta [START_REF] Samanta | A note on uniform strong convergence of bivariate density estimates[END_REF]). It is worth noting that, computational issues aside, the transformation step (i.e., premultiplication of the regressors by M-1 n , where Mn is an affine equivariant scatter estimate) may be based on a statistic Mn that does not require finiteness of any moment. A typical example is the scatter estimate proposed in Tyler [START_REF] Tyler | A distribution-free M-estimator of multivariate scatter[END_REF] or Hettmansperger and Randles [START_REF] Hettmansperger | A practical affine equivariant multivariate median[END_REF]. Rather, our procedure takes ideas from the classical nonparametric literature using concepts such as multivariate ranks. It is closed in spirit of the approach of Paindaveine and Van Bever [START_REF] Paindaveine | Nonparametric consistent depth-based classifiers[END_REF], who introduce a class of depth-based classification procedures that are of a nearest neighbor nature.

There are also attempts at getting invariance to other transformations. The most important concept here is that of invariance under monotone transformations of the coordinate axes. In particular, any strategy that uses only the coordinatewise ranks of the X i 's achieves this. The onus, then, is to show consistency of the methods under the most general conditions possible. For example, using an L p norm on the d-vectors of differences between ranks, one can show that the classical k-NN regression function estimate is universally consistent in the sense of Stone [START_REF] Stone | Consistent nonparametric regression[END_REF]. This was observed by Olshen [START_REF] Olshen | Comments on a paper by C.J. Stone[END_REF], and shown by Devroye [START_REF] Devroye | A universal k-nearest neighbor procedure in discrimination[END_REF] (see also Gordon and Olshen [15,[START_REF] Gordon | Consistent nonparametric regression from recursive partitioning schemes[END_REF], Devroye and Krzyżak [START_REF] Devroye | New multivariate product density estimators[END_REF], and Biau and Devroye [START_REF] Biau | On the layered nearest neighbour estimate, the bagged nearest neighbour estimate and the random forest method in regression and classification[END_REF] for related works). Rules based upon statistically equivalent blocks (see, e.g., Anderson [START_REF] Anderson | Some nonparametric multivariate procedures based on statistically equivalent blocks[END_REF], Quesenberry and Gessaman [START_REF] Quesenberry | Nonparametric discrimination using tolerance regions[END_REF], Gessaman [START_REF] Gessaman | A consistent nonparametric multivariate density estimator based on statistically equivalent blocks[END_REF], Gessaman and Gessaman [START_REF] Gessaman | A comparison of some multivariate discrimination procedures[END_REF], and Devroye, Györfi, and Lugosi [START_REF] Devroye | A Probabilistic Theory of Pattern Recognition[END_REF]Section 21.4]) are other important examples of regression methods invariant with respect to monotone transformations of the coordinate axes. These methods and their generalizations partition the space with sets that contain a fixed number of data points each.

It would be interesting to consider in a future paper the possibility of morphing the input space in more general ways than those suggested in the previous few paragraphs of the present article. It should be possible, in principle, to define appropriate metrics to obtain invariance for interesting large classes of nonlinear transformations, and show consistent asymptotic behaviors.

An affine invariant k-NN estimate

The k-NN estimate we are discussing is based upon the notion of empirical distance. Throughout, we assume that the distribution of X is absolutely continuous with respect to the Lebesgue measure on R d and that n ≥ d.

Because of this density assumption, any collection

X i 1 , . . . , X i d (1 ≤ i 1 < i 2 < . . . < i d ≤ n) of d points among X 1 , .
. . , X n are in general position with probability 1. Consequently, there exists with probability 1 a unique hyperplane in R d containing these d random points, and we denote it by

H(X i 1 , . . . , X i d ).
With this notation, the empirical distance between d-vectors x and x ′ is defined as

ρ n (x, x ′ ) = 1≤i 1 <...<i d ≤n 1 {segment (x,x ′ ) intersects the hyperplane H(X i 1 ,...,X i d )} .
Put differently, ρ n (x, x ′ ) just counts the number of hyperplanes in R d passing through d out of the points X 1 , . . . , X n , that are separating x and x ′ . Roughly, "near" points have fewer intersections, see Figure 1 that depicts an example in dimension 2. x Figure 1: An example in dimension 2. The empirical distance between x and x ′ is 4. (Note that the hyperplane defined by the pair [START_REF] Chakraborty | Operating transformation retransformation on spatial median and angle test[END_REF][START_REF] Cover | Estimation by the nearest neighbor rule[END_REF] indeed cuts the segment (x, x ′ ), so that the distance is 4, not 3.) This hyperplane-based concept of distance is known in the multivariate rank tests literature as the empirical lift-interdirection function (Oja and Paindaveine [START_REF] Oja | Optimal signed-rank tests based on hyperplanes[END_REF], see also Randles [START_REF] Randles | A distribution-free multivariate sign test based on interdirections[END_REF], Oja [START_REF] Oja | Affine invariant multivariate sign and rank tests and corresponding estimates: A review[END_REF], and Hallin and Paindaveine [START_REF] Hallin | Optimal tests for multivariate location based on interdirections and pseudo-Mahalanobis ranks[END_REF] for companion concepts). It was originally mentioned (but not analyzed) in Hettmansperger, Möttönen, and Oja [START_REF] Hettmansperger | The geometry of the affine invariant multivariate sign and rank methods[END_REF], and independently suggested as an affine invariant alternative to ordinary metrics in the monograph of Devroye, Gÿorfi, and Lugosi [START_REF] Devroye | A Probabilistic Theory of Pattern Recognition[END_REF]Section 11.6]. We speak throughout of distance even though, for a fixed sample of size n, ρ n is only defined with probability 1 and is not a distance measure stricto sensu (in particular, ρ n (x, x ′ ) = 0 does not imply that x = x ′ ). Nevertheless, this empirical distance is invariant under affine transformations x → Ax + b, where A is some arbitrary nonsingular linear map and b any offset vector (see, for instance, Oja and Paindaveine [START_REF] Oja | Optimal signed-rank tests based on hyperplanes[END_REF]Section 2.4]). Now, fix x ∈ R d and let ρ n (x, X i ) be the empirical distance between x and some observation X i in the sample X 1 , . . . , X n . (That is, the number of hyperplanes in R d passing through d out of the observations X 1 , . . . , X n , that are cutting the segment (x, X i )). In this context, the k-NN estimate we are considering still takes the familiar form

r n (x) = 1 k n kn i=1 Y (i) (x),
with the important difference that now the data set (X

p ≥ 1, if k n → ∞ and k n /n → 0, E |r n (x) -r(x)| p → 0 as n → ∞.
The following corollary is a consequence of Theorem 2.1 and the Lebesgue dominated convergence theorem.

Corollary 2.1 (Global L p consistency) Assume that X has a probability density, that Y is bounded, and that the regression function r is µ-almost surely continuous. Then, for all

p ≥ 1, if k n → ∞ and k n /n → 0, E |r n (X) -r(X)| p → 0 as n → ∞.
The conditions of Stone's universal consistency theorem given in [START_REF] Stone | Consistent nonparametric regression[END_REF] are not fulfilled for our estimate. For the standard nearest neighbor estimate, a key result used in the consistency proof by Stone is that a given data point cannot be the nearest neighbor of more than a constant number (say, 3 d ) other points. Such a universal constant does not exist after our transformation is applied. That means that a single data point can have a large influence on the regression function estimate. While this by itself does not imply that the estimate is not universally consistent, it certainly indicates that any such proof will require new insights. The addition of two smoothness constraints, namely that X has a density (without, however, imposing any continuity conditions on the density itself) and that r is µ-almost surely continuous, is sufficient.

The complexity of our procedure in terms of sample size n and dimension d is quite high. There are n d possible choices of hyperplanes through d points. This collection of hyperplanes defines an arrangement, or partition of R d into polytopal regions, also called cells or chambers. Within each region, the distance to each data point is constant, and thus, a preprocessing step might consist of setting up a data structure for determining to which cell a given point x ∈ R d belongs: This is called the point location problem. Meiser [START_REF] Meiser | Point location in arrangements of hyperplanes[END_REF] showed that such a data structure exists with the following properties: (1) it takes space O(n d+ε ) for any fixed ε > 0, and (2) point location can be performed in O(log n) time. Chazelle's cuttings [START_REF] Chazelle | Cutting hyperplanes for divide-and-conquer[END_REF] improve [START_REF] Anderson | Some nonparametric multivariate procedures based on statistically equivalent blocks[END_REF] to O(n d ). Chazelle's processing time for setting up the data structure is O(n d ). Still in the preprocessing step, one can determine for each cell in the arrangement the distances to all n data points: This can be done by walking across the graph of cells or by brute force. When done naively, the overall set-up complexity is O(n 2d+1 ). For each cell, one might keep a pointer to the k nearest neighbors. Therefore, once set up, the computation of the regression function estimate takes merely O(log n) time for point location, and O(k) time for retrieving the k nearest neighbors.

One could envisage a reduction in the complexity by defining the distances not in terms of all hyperplanes that cut a line segment, but in terms of the number of randomly drawn hyperplanes that make such a cut, where the number of random draws is now a carefully selected number. By the concentration of binomial random variables, such random estimates of the distances are expected to work well, while keeping the complexity reasonable. This idea will be explored elsewhere.

Proof of the theorem

Recall, since X has a probability density with respect to the Lebesgue measure on R d , that any collection X i 1 , . . . ,

X i d (1 ≤ i 1 < i 2 < . . . < i d ≤ n) of d points among X 1 , . . . , X n defines with probability 1 a unique hyperplane H(X i 1 , . . . , X i d ) in R d .
Thus, in the sequel, since no confusion is possible, we will freely refer to "the hyperplane H(X i 1 , . . . , X i d ) defined by X i 1 , . . . , X i d " without further explicit mention of the probability 1 event.

Let us first fix some useful notation. The distribution of the random variable X is denoted by µ and its density with respect to the Lebesgue measure is denoted by f . For every ε > 0, we let B x,ε = {y ∈ R d : yx ≤ ε} be the closed Euclidean ball with center at x and radius ε. We write A c for the complement of a subset A of R d . For two random variables Z 1 and Z 2 , the notation

Z 1 ≤ st Z 2 means that Z 1 is stochastically dominated by Z 2 , that is, for all t ∈ R, P{Z 1 > t} ≤ P{Z 2 > t}.
Our first goal is to show that for µ-almost all x, as k n /n → 0, the quantity max i=1,...,kn X (i) (x)x converges to 0 in probability, i.e., for every ε > 0, lim n→∞ P max i=1,...,kn

X (i) (x) -x > ε = 0. (3.1)
So, fix such a positive ε. Let δ be a real number in (0, ε) and γ n be a positive real number (eventually function of x and ε) to be determined later. To prove identity (3.1), we use the following decomposition, which is valid for all x ∈ R d : P max i=1,...,kn

X (i) (x) -x > ε ≤ P    min i=1,...,n X i ∈B c x,ε ρ n (x, X i ) < γ n    + P    max i=1,...,n X i ∈B x,δ ρ n (x, X i ) ≥ γ n    + P {Card {i = 1, . . . , n : X i -x ≤ δ} < k n } := A + B + C. (3.2)
The convergence to 0 of each of the three terms above-from which identity (3.1) immediately follows-are separately analyzed in the next three paragraphs.

Analysis of A.

As for now, taking an affine geometry point of view, we keep x fixed and see it as the origin of the space. Recall that each point in the Euclidean space R d (with the origin at x) may be described by its hyperspherical coordinates (see, e.g., Miller [25, Chapter 1]), which consist of a nonnegative radial coordinate r and d -1 angular coordinates θ 1 , . . . , θ d-1 , where θ d-1 ranges over [0, 2π) and the other angles range over [0, π] (adaptation of this definition to the cases d = 1 and d = 2 is clear). For a (d -1)-dimensional vector Θ = (θ 1 , . . . , θ d-1 ) of hyperspherical angles, we let B x,ε (Θ) be the unique closed ball anchored at x in the direction Θ and with diameter ε (see Figure 2 which depicts an illustration in dimension 2). We also let L x (Θ) be the axe defined by x and the direction Θ, and let as well S x,ε (Θ) be the open segment obtained as the intersection of L x (Θ) and the interior of B x,ε (Θ).

x

Θ ε 2 N 2 x,ε (Θ) N 1 x,ε (Θ) R 2 x,ε (Θ) R 1 x,ε (Θ) S x,ε (Θ) L x (Θ)
Figure 2: The ball B x,ε (Θ) and related notation. Illustration in dimension 2.

Next, for fixed x, ε and Θ, we split the ball

B x,ε (Θ) into 2 d-1 disjoint regions R 1 x,ε (Θ), . . . , R 2 d-1
x,ε (Θ) as follows. First, the Euclidean space R d is sequentially divided into 2 d-1 symmetric quadrants rotating around the axe L x (Θ) (boundary equalities are broken arbitrarily). Next, each region R j

x,ε (Θ) is obtained as the intersection of one of the 2 d-1 quadrants and the ball B x,ε (Θ).

The numbers of sample points falling in each of these regions are denoted hereafter by N 1

x,ε (Θ), . . . , N 2 d-1

x,ε (Θ) (see Figure 2). Letting finally V d be the volume of the unit d-dimensional Euclidean ball, we are now in a position to control the first term of inequality (3.2). Proposition 3.1 For µ-almost all x ∈ R d and all ε > 0 small enough,

P    min i=1,...,n X i ∈B c x,ε ρ n (x, X i ) < γ n    → 0 as n → ∞, provided γ n = n d V d 2 2d+1 ε d f (x) 2 d-1 . Proof of Proposition 3.1 Set p x,ε = min j=1,...,2 d-1 inf Θ µ R j x,ε (Θ) ,
where the infimum is taken over all possible hyperspherical angles Θ. We know, according to technical Lemma 4.1, that for µ-almost all x and all ε > 0 small enough,

p x,ε ≥ V d 2 2d ε d f (x) > 0. (3.3)
Thus, in the rest of the proof, we fix such an x and assume that ε is small enough so that the inequalities above are satisfied.

Let X ⋆ be defined as the intersection of the line (x, X) with B x,ε , and let Θ ⋆ be the (random) hyperspherical angle corresponding to X ⋆ (see Figure 3 for an example in dimension 2).

x X

X ⋆ ε N 2 x,ε (Θ ⋆ ) N 1 x,ε (Θ ⋆ ) Θ ⋆ Figure 3: The ball B x,ε (Θ ⋆ ) in dimension 2.
Denote by N x,ε (Θ ⋆ ) the number of hyperplanes passing through d out of the observations X 1 , . . . , X n and cutting the segment S x,ε (Θ ⋆ ). We have

P    min i=1,...,n X i ∈B c x,ε ρ n (x, X i ) < γ n    ≤ nP {ρ n (x, X ⋆ ) < γ n } = nP {N x,ε (Θ ⋆ ) < γ n } ≤ nP N 1 x,ε (Θ ⋆ ) . . . N 2 d-1 x,ε (Θ ⋆ ) n 2 d-1 -d < γ n ,
where the last inequality follows from technical Lemma 4.2. Thus,

P    min i=1,...,n X i ∈B c x,ε ρ n (x, X i ) < γ n    ≤ n 2 d-1 j=1 P N j x,ε (Θ ⋆ ) < γ n n 2 d-1 -d 1/2 d-1 = n 2 d-1 j=1 P N j x,ε (Θ ⋆ ) < γ 1/2 d-1 n n 1-d/2 d-1 .
Clearly, conditionally on Θ ⋆ , each

N j x,ε (Θ ⋆ ) satisfies Binomial (n, p x,ε ) ≤ st N j x,ε (Θ ⋆ )
and consequently, by inequality

(3.3), Binomial n, V d 2 2d ε d f (x) ≤ st N j x,ε (Θ ⋆ ).
Thus, for each j = 1, . . . , 2 d-1 , by Hoeffding's inequality for binomial random variables (Hoeffding [START_REF] Hoeffding | Probability inequalities for sums of bounded random variables[END_REF]), we are led to

P N j x,ε (Θ ⋆ ) < γ 1/2 d-1 n n 1-d/2 d-1 = E P N j x,ε (Θ ⋆ ) < γ 1/2 d-1 n n 1-d/2 d-1 |Θ ⋆ ≤ exp -2 γ 1/2 d-1 n n 1-d/2 d-1 -n V d 2 2d ε d f (x) 2 /n
as soon as γ

1/2 d-1 n n 1-d/2 d-1 < n V d 2 2d ε d f (x). Therefore, taking γ n = n d V d 2 2d+1 ε d f (x) we obtain P    min i=1,...,n X i ∈B c x,ε ρ n (x, X i ) < γ n    ≤ 2 d-1 n exp -n V d 2 2d ε d f (x) 2 /2 .
The upper bound goes to 0 as n → ∞.

Analysis of B.

Consistency of the second term in inequality (3.2) is established in the following proposition.

Proposition 3.2 For µ-almost all x ∈ R d , all ε > 0 and all δ > 0 small enough,

P    max i=1,...,n X i ∈B x,δ ρ n (x, X i ) ≥ γ n    → 0 as n → ∞, provided γ n = n d V d 2 2d+1 ε d f (x) 2 d-1 . (3.4)
Proof of Proposition 3.2 Fix x in a set of µ-measure 1 such that f (x) > 0 and denote by N x,δ the number of hyperplanes that cut the ball B x,δ . Clearly,

P    max i=1,...,n X i ∈B x,δ ρ n (x, X i ) ≥ γ n    ≤ P {N x,δ ≥ γ n } .
Observe that, with probability 1,

N x,δ = 1≤i 1 <...<i d ≤n 1 {H(X i 1 ,...,X i d )∩B x,δ =∅} ,
whence, since X 1 , . . . , X n are identically distributed,

E[N x,δ ] = n d P {H(X 1 , . . . , X d ) ∩ B x,δ = ∅} ≤ n d d! P {H(X 1 , . . . , X d ) ∩ B x,δ = ∅} .
Consequently, given the choice (3.4) for γ n and the result of technical Lemma 4.3, it follows that

E[N x,δ ] < γ n /2
for all δ small enough, independently of n. Thus, using the bounded difference inequality (McDiarmid [START_REF] Mcdiarmid | On the method of bounded differences[END_REF]), we obtain, still with the choice

γ n = n d V d 2 2d+1 ε d f (x) 2 d-1 , P {N x,δ ≥ γ n } ≤ P {N x,δ -E[N x,δ ] ≥ γ n /2} ≤ exp -2 (γ n /2) 2 n 2d-1 = exp - V d 2 2d+1 ε d f (x) 2 d n/2 .
This upper bound goes to zero as n tends to infinity, and this concludes the proof of the proposition.

Analysis of C.

To achieve the proof of identity (3.1), it remains to show that the third and last term of (3.2) converges to 0. This is done in the following proposition. Proposition 3.3 Assume that k n /n → 0 as n → ∞. Then, for µ-almost all x ∈ R d and all δ > 0,

P {Card {i = 1, . . . , n : X i -x ≤ δ} < k n } → 0 as n → ∞.
Proof of Proposition 3.3 Recall that the collection of all x with µ(B x,τ ) > 0 for all τ > 0 is called the support of µ, and note that it may alternatively be defined as the smallest closed subset of R d of µ-measure 1 (Parthasarathy [START_REF] Parthasarathy | Probability Measures on Metric Spaces[END_REF]Chapter 2]). Thus, fix x in the support of µ and set p x,δ = P{X ∈ B x,δ }, so that p x,δ > 0. Then the following chain of inequalities is valid:

P {Card {i = 1, . . . , n : X i -x ≤ δ} < k n } = P {Binomial (n, p x,δ ) < k n } ≤ P {Binomial (n, p x,δ ) ≤ np x,δ /2} (for all n large enough, since k n /n tends to 0) ≤ exp(-np 2 x,δ /2),
where the last inequality follows from Hoeffding's inequality (Hoeffding [21]). This terminates the proof of Proposition 3.3.

We have proved so far that, for µ-almost all x, as k n /n → 0, the quantity max i=1,...,kn X (i) (x)x converges to 0 in probability. By the elementary inequality

E 1 k n kn i=1
1 { X (i) (x)-x >ε} ≤ P max i=1,...,kn

X (i) (x) -x > ε ,
it immediately follows that, for such an x,

E 1 k n kn i=1 1 { X (i) (x)-x >ε} → 0 (3.5)
provided k n /n → 0. We are now ready to complete the proof of Theorem 2.1.

Fix x in a set of µ-measure 1 such that consistency (3.5) holds and r is continuous at x (this is possible by the assumption on r). Because |a + b| p ≤ 2 p-1 (|a| p + |b| p ) for p ≥ 1, we see that

E |r n (x) -r(x)| p ≤ 2 p-1 E 1 k n kn i=1 Y (i) (x) -r X (i) (x) p + 2 p-1 E 1 k n kn i=1 r X (i) (x) -r(x) p .
Thus, by Jensen's inequality,

E |r n (x) -r(x)| p ≤ 2 p-1 E 1 k n kn i=1 Y (i) (x) -r X (i) (x) p + 2 p-1 E 1 k n kn i=1 r X (i) (x) -r(x) p := 2 p-1 I n + 2 p-1 J n .
Firstly, for arbitrary ε > 0, we have

J n = E 1 k n kn i=1 r X (i) (x) -r(x) p 1 { X (i) (x)-x >ε} + E 1 k n kn i=1 r X (i) (x) -r(x) p 1 { X (i) (x)-x ≤ε} , whence J n ≤ 2 p ζ p E 1 k n kn i=1 1 { X (i) (x)-x >ε} + sup y∈R d : y-x ≤ε |r(y) -r(x)| p (since |Y | ≤ ζ).
The first term on the right-hand side of the latter inequality tends to 0 by (3.5) as k n /n → 0, whereas the rightmost one can be made arbitrarily small as ε → 0 since r is continuous at x. This proves that J n → 0 as n → ∞.

Next, by successive applications of inequalities of Marcinkiewicz and Zygmund [START_REF] Marcinkiewicz | Sur les fonctions indépendantes[END_REF] (see also Petrov [33, pages 59-60]), we have for some positive constant C p depending only on p,

I n ≤ C p E 1 k 2 n kn i=1 Y (i) (x) -r X (i) (x) 2 p/2 ≤ (2ζ) p C p k p/2 n (since |Y | ≤ ζ).
Consequently, I n → 0 as k n → ∞, and this concludes the proof of the theorem.

Some technical lemmas

The notation of this section is identical to that of Section 3. In particular, it is assumed throughout that X has a probability density f with respect to the Lebesgue measure λ on R d . This requirement implies that any collection

X i 1 , . . . , X i d (1 ≤ i 1 < i 2 < . . . < i d ≤ n) of d points among X 1 , . . . , X n define with probability 1 a unique hyperplane H(X i 1 , . . . , X i d ) in R d .
Recall finally that, for x ∈ R d and ε > 0, we set

p x,ε = min j=1,...,2 d-1 inf Θ µ R j x,ε (Θ) ,
where the infimum is taken over all possible hyperspherical angles Θ, and the regions R j x,ε (Θ), j = 1, . . . , 2 d-1 , define a partition of the ball B x,ε (Θ). Recall also that the numbers of sample points falling in each of these regions are denoted by N 1

x,ε (Θ), . . . , N 2 d-1

x,ε (Θ). For a better understanding of the next lemmas, the reader should refer to Figure 2

p x,ε ≥ V d 2 2d ε d f (x) > 0.
Proof of Lemma 4.1 We let x be a Lebesgue point of f , that is, an x such that for any collection A of subsets of B 0,1 with the property that for all A ∈ A, λ(A) ≥ cλ(B 0,1 ) for some fixed c > 0,

lim ε→0 sup A∈A x+εA f (y)dy λ{x + εA} -f (x) = 0, (4.1)
where x + εA = {y ∈ R d : (yx)/ε ∈ A}. As f is a density, we know that µ-almost all x satisfy this property (see, for instance, Wheeden and Zygmund [START_REF] Wheeden | Measure and Integral. An Introduction to Real Analysis[END_REF]). Moreover, since f is µ-almost surely positive, we may also assume that f (x) > 0.

Thus, keep such an x fixed. Fix also j ∈ {1, . . . , 2 d-1 }, and set

p j x,ε = inf Θ µ R j x,ε (Θ) .
Taking for A the collection of regions R j 0,1 (Θ) when the hyperspherical angle Θ varies, that is,

A = R j 0,1 (Θ) : Θ ∈ [0, π] d-2 × [0, 2π) , and observing that λ R j x,ε (Θ) = V d 2 d-1 ε 2 d ,
we may write, for each j = 1, . . . ,

2 d-1 , 2 d-1 p j x,ε V d (ε/2) d -f (x) = inf Θ µ R j x,ε (Θ) λ R j x,ε (Θ) -f (x) = inf A∈A x+εA f (y)dy λ{x + εA} -f (x) ≤ sup A∈A x+εA f (y)dy λ{x + εA} -f (x) .
The conclusion follows from identity (4.1).

Lemma 4.2 Fix x ∈ R d , ε > 0 and Θ ∈ [0, π] d-2 ×[0, 2π). Let N x,ε (
Θ) be the number of hyperplanes passing through d out of the observations X 1 , . . . , X n and cutting the segment S x,ε (Θ). Then, with probability 1,

N x,ε (Θ) ≥ N 1 x,ε (Θ) . . . N 2 d-1 x,ε (Θ) n 2 d-1 -d . Proof of Lemma 4.2 If one of the N j x,ε (Θ) (j = 1, . . . , 2 d-1
) is zero, then the result is trivial. Thus, in the rest of the proof, we suppose that each N j

x,ε (Θ) is positive and note that this implies n

≥ 2 d-1 . Pick sequentially 2 d-1 observations, say X i 1 , . . . , X i 2 d-1 , in the 2 d-1 regions R 1 x,ε (Θ), . . . , R 2 d-1
x,ε (Θ). By construction, the polytope defined by these 2 d-1 points cuts the axe L x,ε (Θ). Consequently, with probability 1, any hyperplane drawn according to d out of these 2 d-1 points cuts the segment S x,ε (Θ). The result follows by observing that there are exactly N 1

x,ε (Θ) . . . N 2 d-1

x,ε (Θ) such polytopes. Lemma 4.3 For 1 ≤ i 1 < . . . < i d ≤ n, let H(X i 1 , . . . , X i d ) be the hyperplane passing through d out of the observations X 1 , . . . , X n . Then, for all x ∈ R d , P {H(X i 1 , . . . , X i d ) ∩ B x,δ = ∅} → 0 as δ ↓ 0.

Proof of Lemma 4.3 Given two hyperplanes H and H ′ in R d , we denote by Φ(H, H ′ ) the (dihedral) angle between H and H ′ . Recall that Φ(H, H ′ ) ∈ [0, π] and that it is defined as the angle between the corresponding normal vectors.

Fix 1 ≤ i 1 < . . . < i d ≤ n. Let E δ be the event

E δ = X i j -x > δ : j = 1, . . . , d -1 ,
and let H(x, X i 1 , . . . , X i d-1 ) be the hyperplane passing through x and the d-1 points X i 1 , . . . , X i d-1 . Clearly, on E δ , the event {H(X i 1 , . . . , X i d )∩B x,δ = ∅} is the same as Φ H(x, X i 1 , . . . , X i d-1 ), H(X i 1 , . . . ,

X i d ) ≤ Φ δ ,
where Φ δ is the angle formed by H(x, X i 1 , . . . , X i d-1 ) and the hyperplane going trough X i 1 , . . . , X i d-1 and tangent to B x,δ (see Figure 4 for an example in dimension 2).

Thus, with this notation, we may write P {H(X i 1 , . . . , X i d ) ∩ B x,δ = ∅} ≤ P{E c n } + P {H(X i 1 , . . . , X i d ) ∩ B x,δ = ∅, E n } ≤ P{E c n } + P Φ H(x, X i 1 , . . . , X i d-1 ), H(X i 1 , . . . , X i d ) ≤ Φ δ , E n .

X i 1 X i 2 δ x Φ δ Figure 4:
The hyperplanes H(x, X i 1 , . . . , X i d-1 ) and H(X i 1 , . . . , X i d ), and the angle Φ δ . Illustration in dimension 2.

Since X has a density, the first of the two terms above tends to zero as δ ↓ 0. To analyze the second term, first note that, conditionally on X i 1 , . . . , X i d-1 , the angle Φ(H(x, X i 1 , . . . , X i d-1 ), H(X i 1 , . . . , X i d )) is absolutely continuous with respect to the Lebesgue measure on R. This follows from the following two observations: i) the random variable X i d has a density with respect to the Lebesgue measure on R d , and ii) conditionally on X i 1 , . . . , X i d-1 , Φ(H(x, X i 1 , . . . , X i d-1 ), H(X i 1 , . . . , X i d )) is obtained from X i d via translations, orthogonal transformations, and the arctan function.

Thus, writing P Φ H(x, X i 1 , . . . , X i d-1 ), H(X i 1 , . . . , X i d ) ≤ Φ δ , E n = E 1 En P Φ H(x, X i 1 , . . . , X i d-1 ), H(X i 1 , . . . ,

X i d ) ≤ Φ δ |X i 1 , . . . , X i d-1
and noting that, on the event E n , for fixed X i 1 , . . . , X i d-1 , Φ δ ↓ 0 as δ ↓ 0, we conclude by the Lebesgue dominated convergence theorem that P Φ H(x, X i 1 , . . . , X i d-1 ), H(X i 1 , . . . , X i d ) ≤ Φ δ → 0 as δ ↓ 0.

  and Figure3.

Lemma 4 . 1

 41 For µ-almost all x ∈ R d and all ε > 0 small enough,

  1 , Y 1 ), . . . , (X n , Y n ) is reordered according to increasing values of the empirical distances ρ n (x, X

i ), not the original Euclidean metric. By construction, the estimate r n has the desired affine invariance property and, moreover, it coincides with the standard (Euclidean) estimate in dimension d = 1. In the next section, we prove the following theorem. The distribution of the random variable X is denoted by µ. Theorem 2.1 (Pointwise L p consistency) Assume that X has a probability density, that Y is bounded, and that the regression function r is µalmost surely continuous. Then, for µ-almost all x ∈ R d and all

d-1 ,

Acknowledgments. We thank two anonymous referees for valuable comments and insightful suggestions.

 [START_REF] Biau | On the layered nearest neighbour estimate, the bagged nearest neighbour estimate and the random forest method in regression and classification[END_REF]Research partially supported by the French National Research Agency under grant ANR-09-BLAN-0051-02 "CLARA". 3 Research carried out within the INRIA project "CLASSIC" hosted by Ecole Normale Supérieure and CNRS. 4 Research sponsored by NSERC Grant A3456 and FQRNT Grant 90-ER-0291. 5 Research sponsored by NSERC Grant RGPIN 402438-2011. 6 Research sponsored by NSERC Grant N00118.

1