Cristina Silvano

Giovanni Agosta

Andrea Bartolini

Andrea R Beccari
email: andrea.beccari@dompe.it

Luca Benini
email: lbenini@iis.ee.ethz.ch

Loïc Besnard

João Bispo

Radim Cmar
email: rcmar@sygic.com

Joao M R Cardoso

Carlo Cavazzoni

João M P Cardoso

Stefano Cherubin

Davide Gadioli

Martin Golasowski

Imane Lasri

Jan Martinovič

Gianluca Palermo

Pedro Pinto

Erven Rohou

Nico Sanna

Kateřina Slaninová

Emanuele Vitali

ANTAREX: A DSL-based Approach to Adaptively Optimizing and Enforcing Extra-Functional Properties in High Performance Computing

Keywords: High Performance Computing, Autotuning, Adaptivity, DSL, Compilers, Energy Efficiency

published or not. The documents may come

L'archive ouverte

I. INTRODUCTION

Designing and optimising applications for energy-efficient High Performance Computing systems up to the Exascale era is an extremely challenging problem. Exascale supercomputers (reaching billions of billions of floating point operations per second) cannot be built by simply expanding the number of processing nodes and leveraging technology scaling, as power demand would increase unsustainably (up to hundreds of MW). To reach the DARPA target of 20MW of Exascale supercomputers for the year 2020, current supercomputers (reaching up to 93 PetaFlop/s) must achieve an energy efficiency "quantum leap". The Green500 list 1 looks at the GigaFlops per Watt as energy efficiency metric to rank supercomputers by their energy efficiency. According to the latest Green500 list published in November 2017, the "most green" supercomputer SHOUBU SystemB installed in Japan achieves 17 GigaFlops/W during its 842-TeraFlop/s Linpack performance run. The top positions in Green500 are all occupied by heterogeneous systems based on high-performance processors and co-processors such as the latest NVIDIA Volta GV100 GPU and PEZY SC2 accelerator to further accelerate the computation. The dominance of heterogeneous systems in the Green500 list is expected to continue for the next coming years to reach the target of 20MW Exascale supercomputers. To this end, European efforts have recently been focused on building supercomputers out of the less power-hungry ARM cores and GPGPUs [START_REF] Rajovic | Supercomputing with Commodity CPUs: Are Mobile SoCs Ready for HPC?[END_REF].

Designing and implementing HPC applications are difficult and complex tasks, which require mastering several specialized languages and tools for performance tuning. This is incompatible with the current trend to open HPC infrastructures to a wider range of users. The current model where the HPC center staff directly supports the development of applications will become unsustainable in the 1 www.green500.org, November 2017 long term. Thus, the availability of effective standard programming languages and APIs is crucial to provide migration paths towards novel heterogeneous HPC platforms as well as to guarantee the ability of developers to work effectively on these platforms.

To fulfil the 20MW target, energy-efficient heterogeneous supercomputers need to be coupled with radically new software stacks to exploit the benefits offered by heterogeneity at all levels (supercomputer, job, node).

The ANTAREX [START_REF] Silvano | AutoTuning and Adaptivity appRoach for Energy efficient eXascale HPC systems: the ANTAREX Approach[END_REF][START_REF] Silvano | The ANTAREX approach to autotuning and adaptivity for energy efficient HPC systems[END_REF][START_REF] Silvano | The ANTAREX Tool Flow for Monitoring and Autotuning Energy Efficient HPC Systems[END_REF] project aims at providing a holistic approach spanning all the decision layers composing the supercomputer software stack and exploiting effectively the full system capabilities, including heterogeneity and energy management. The main goal of ANTAREX is to express by means of a DSL the application selfadaptivity and to runtime manage and autotune applications for green heterogeneous HPC systems up to the Exascale level. The use of a DSL allows the introduction of a separation of concerns, where selfadaptivity and energy efficient strategies are specified separately from the application functionalities. The DSL is based on previous efforts regarding the LARA language [START_REF] Cardoso | LARA: An Aspect-oriented Programming Language for Embedded Systems[END_REF][START_REF] Cardoso | Performance-driven instrumentation and mapping strategies using the LARA aspect-oriented programming approach[END_REF] and makes possible to express at compile time the adaptivity/energy/performance strategies and to enforce at runtime application autotuning and resource and power management. The goal is to support the parallelism, scalability and adaptivity in a dynamic workload by exploiting the full system capabilities (including energy management) for emerging large-scale and extreme-scale systems, while reducing the Total Cost of Ownership (TCO) for companies and public organizations.

The project is driven by two use cases taken from highly relevant HPC application scenarios: (1) a biopharmaceutical application for drug discovery deployed on the 1.21 PetaFlops heterogeneous NeXtScale Intel-based IBM system at CINECA and (2) a selfadaptive navigation system for smart cities deployed on the serverside on the 1.46 PetaFlops heterogeneous Intel® Xeon Phi™ based system provided by IT4Innovations National Supercomputing Center.

The ANTAREX Consortium comprises a wealth of expertise in all pertinent domains. Four top-ranked academic and research partners (Politecnico di Milano, ETHZ Zurich, University of Porto and INRIA) are complemented by the Italian Tier-0 Supercomputing Center (CINECA), the Tier-1 Czech National Supercomputing Center (IT4Innovations) and two industrial application providers, one of the leading biopharmaceutical companies in Europe (Dompé) and the top European navigation software company (Sygic). The complementarity and deep expertise of the Consortium partners has the potential to generate a breakthrough innovation from the ANTAREX project. Moreover, the presence of leading edge industrial partners will ensure The ANTAREX approach and related tool flow, as shown in Figure 1, operate both at design-time and runtime. The application functionality is expressed through C/C++ code (possibly including legacy code), whereas the extra-functional aspects of the application, including parallelisation, mapping, and adaptivity strategies, are expressed through DSL code (based on LARA) developed in the project. As a result, the expression of such aspects is fully decoupled from the functional code. The Clava tool is the centerpoint of the compiletime phase, performing a refactoring of the application code based on the LARA aspects, and instrumenting it with the necessary calls to other components of the tool flow.

The ANTAREX compilation flow leverages a runtime phase with compilation steps, through the use of partial dynamic compilation techniques enabled by libVC. The application autotuning, performed via the mARGOt tool, is delayed to the runtime phase, where the software knobs (application parameters, code transformations and code variants) are configured according to the runtime information coming from application self-monitoring as well as from system monitoring performed by the ExaMon tool. Finally, the runtime power manager, PowerCapper, is used to control the resource usage for the underlying computing infrastructure given the changing conditions. At runtime, the application control code, thanks to the design-time phase, now contains also runtime monitoring and adaptivity strategy code derived from the DSL extra-functional specification. Thus, the application is continuously monitored to guarantee the required Service Level Agreement (SLA), while communication with the runtime resourcemanager takes place to control the amount of processing resources needed by the application. The application monitoring and autotuning is supported by a runtime layer implementing an application level collect-analyse-decide-act loop.

The rest of this paper is organized as follows. In Section II we review the technology portfolio provided by the ANTAREX tool flow. In Section III we provide an assessment of the impact of the proposed DSL on application specifications, while Section IV gives an overview of how the Tool Flow has been used in one of the use cases. Finally, in Section V we draw some conclusions.

II. ANTAREX TECHNOLOGY PORTFOLIO

A. The ANTAREX DSL HPC applications might profit from adapting to operational and situational conditions, such as changes in contextual information (e.g., workloads), in requirements (e.g., deadlines, energy), and in availability of resources (e.g., connectivity, number of processor nodes available). A simplistic approach to both adaptation specification and implementation (see, e.g., [START_REF] Floch | Using architecture models for runtime adaptability[END_REF]) employs hard coding of, e.g., conditional expressions and parameterizations. In our approach, the specification of runtime adaptability strategies relies on a DSL implementing key concepts from Aspect-Oriented Programming (AOP) [START_REF] Irwin | Aspect-oriented Programming[END_REF], mainly specifying adaptation concerns, targeting specific execution points, separately from the primary functionality of the application, with minimum or no changes to the application source code.

Our approach is based on the idea that certain application/system requirements (e.g., target-dependent optimizations, adaptivity behavior and concerns) should be specified separately from the source code that defines the main functionality. Those requirements are expressed as DSL aspects that embody strategies. An extra compilation step, performed by a weaver, merges the original source code and aspects into the intended program [START_REF] Elrad | Aspect-oriented Programming: Introduction[END_REF]. Using aspects to separate concerns from the core objective of the program can result in cleaner programs and increased productivity (e.g., reusability of strategies). As the development process of HPC applications typically involves two types of experts (application-domain experts and HPC system architects) that split their responsibilities along the boundary of functional description and extra-functional aspects, our DSL-aided toolflow provides a suitable approach for helping to express their concerns.

The ANTAREX DSL relies on the already existing DSL technology LARA [START_REF] Cardoso | LARA: An Aspect-oriented Programming Language for Embedded Systems[END_REF][START_REF] Cardoso | Performance-driven instrumentation and mapping strategies using the LARA aspect-oriented programming approach[END_REF]. In particular, the LARA technology provides a framework that we adopted to implement the ANTAREX aspects and APIs. Moreover, we developed other LARA-related tools such as the Clava2 weaver to leverage the rest of the ANTAREX tool flow.

LARA is a programming language that allows developers to capture non-functional requirements and concerns in the form of strategies, which are decoupled from the functional description of the application. Compared to other approaches that usually focus on code injection (e.g., [START_REF] Spinczyk | AspectC++: An Aspect-oriented Extension to the C++ Programming Language[END_REF]), LARA provides access to other types of actions, e.g., code refactoring, compiler optimizations, and inclusion of additional information, all of which can guide compilers to generate more efficient implementations. Additional types of actions may be defined in the language specification and associated weaver, such as software/hardware partitioning [START_REF] Cardoso | Controlling a complete hardware synthesis toolchain with LARA aspects[END_REF] or compiler optimization sequences [START_REF] Nobre | Use of Previously Acquired Positioning of Optimizations for Phase Ordering Exploration[END_REF]. One important feature of the LARA-aided sourceto-source compiler developed in ANTAREX is the capability to refactor the code of the application in order to expose adaptivity behavior and/or adaptivity design points that can be explored by the ANTAREX autotuning component. In the following sections we show illustrative examples3 of some of the strategies that can be specified using LARA in the context of a source-to-source compiler and currently used for one of the use cases.

B. Precision Tuning

Error-tolerating applications are increasingly common in the emerging field of real-time HPC, allowing to trade-off precision for performance and/or energy. Thus, recent works investigated the use of customized precision in HPC as a way to provide a breakthrough in power and performance. We developed a set of LARA aspects enabling mixed precision tuning on C/C++ and OpenCL kernels. In our precision tuning we combine an adaptive selection of floating and fixed point arithmetic, targeting HPC applications.

Figure 2 presents part of a LARA strategy that changes all declarations of a certain type to a target type (e.g., from double to float) for a given function. We note, however, that a practical and reusable aspect needs to deal with further issues, such as the cloning of functions whose types we want to change but are also called by other unrelated functions in the code, assignments of constants, casts, recursion, changing functions definitions and library functions to the ones related to the type used (e.g., sqrtf vs sqrt in Math.h), etc. In this example, changeType is a function that analyzes and changes compound types, such as double * and double[]. If the type described in $old is found inside the type of the declaration, it is replaced with the type described in $new. To be more specific, if $old is double, $new is float and $decl.type is double * , the type of the declaration will be changed to float * . If the original declaration type does not contain the $old type, it is not changed. A LARA aspect consists of three main steps. Firstly, one captures the points of interest in the code using a select statement, which in this example selects variable declarations. Then, using the apply statement, one acts over the selected program points. In this case, it will define the types of the captured declared variables, using the type attribute. Finally, we can then specify conditions to constrain the execution of the apply (i.e., only if the declared variable has a specific type). This can be done via conditional statements (ifs) as well as via special condition blocks that constrain the entire apply. LARA promotes modularity and aspect reuse, and supports embedding JavaScript code, to specify more sophisticated strategies. As shown in [START_REF] Nobre | Aspect-driven mixed-precision tuning targeting gpus[END_REF], we support exploration of mixed precision OpenCL kernels by using half, single, and double precision floating point data types. We additionally support fixed point representations through a custom C++ template-based implementation for HPC systems, which has already been used in [START_REF] Cherubin | Implications of Reduced-Precision Computations in HPC: Performance, Energy and Error[END_REF]. In both cases the LARA aspects automatically insert code for proper type conversion before and after the critical section that has been converted to exploit a reduced precision data type.

The LARA aspect in Figure 3 shows the generation of different mixed-precision versions to be dynamically evaluated. It is possible to specify -as input of the aspect -the number of mix combination to generate, and a rule set to filter out precision mix combinations which are very likely to lead to useless and/or not efficient results. We exploit programmer's application domain knowledge by relying on them to define test cases to evaluate the different code versions at runtime. LARA automatically inserts code to dynamically perform the exploration over the space of the generated versions with different precision mix.

C. Code Versioning

One of the strategies supported in the ANTAREX toolflow is the capability to generate versions of a function and to select the one that satisfies certain requirements at runtime. Figure 4 shows an aspect that clones a set of functions and changes the types of the newly generated clones. Each clone has the same name as the original with the addition of a provided suffix. We start with a single user-defined function which is cloned by the aspect CloneFunction (called in line 13). Then, it recursively traverses calls to other functions inside the clone and generates a clone for each of them. Inside the clones, calls to the original functions are changed to calls to the clones instead, building a new call tree with the generated clones. At the end of the aspect CreateFloatVersion (lines 16-17,) we use the previously defined ChangePrecision aspect to change the types of all generate clones.

The aspect Multiversion -in Figure 5 -adapts the source code of the application in order to call the original version of a function or a generated cloned version with a different type, according to the value of a parameter given by the autotuner at runtime. The main aspect calls the previously shown aspect, CreateFloatVersion, which clones the target function and every other function it uses, while also changing their variable types from double to float (using the aspects presented in Figure 4 and Figure 2). This is performed in lines 8-9 of the example. From lines 13 to 34, the Multiversion aspect generates and inserts code in the application that is used as switching mechanism between the two versions. It starts by declaring a variable to be used as a knob by the autotuner, then it generates the code for a switch statement and replaces the statement containing 1 import ChangePrecision; 2 import clava.ClavaJoinPoints; the original call with the generated switch code. Finally, in lines 36-38, the aspect surrounds both calls (original and float version) with timing code. An excerpt of the resulting C code can bee seen in Figure 6.

1 import CreateFloatVersion; 2 import lara.code.Timer; 3 import clava.ClavaJoinPoints; In the ANTAREX toolflow, the capability of providing several versions of the same function is not limited to static features. LIBVERSIONINGCOMPILER [START_REF] Cherubin | libVersioningCompiler: An easy-to-use library for dynamic generation and invocation of multiple code versions[END_REF] versioning of multiple versions of the same compute kernel in a HPC scenario. It can be used to support continuous optimization, code specialization based on the input data or on workload changes, or to dynamically adjust the application, without the burden of a full just-in-time compiler. LIBVC allows a C/C++ compute kernel to be dynamically compiled multiple times while the program is running, so that different specialized versions of the code can be generated and invoked. Each specialized version can be versioned for later reuse. When the optimal parametrization of the compiler depends on the program workload, the ability to switch at runtime between different versions of the same code can provide significant benefits [START_REF] Chen | Evaluating iterative optimization across 1000 datasets[END_REF][START_REF] Tartara | Continuous learning of compiler heuristics[END_REF]. While such versions can be generated statically in the general case, in HPC execution times can be so long that exhaustive profiling may not be feasible. LIBVC instead enables the exploration and tuning of the parameter space of the compiler at runtime. Figure 7 shows an example of usage of LIBVC through LARA, which demonstrates how to specialize a function. The user provides this aspect with a target function call and a set of compilation options. These include compiler flags and possible compiler definitions, e.g., data discovered at runtime, which is used as a compile-time constant in the new version. Based on the target function call, the aspect finds the function definition which is passed to the library. After the options are set, the original function call is replaced with a call of the newly compiled and loaded specialized version of the kernel.

1 import antarex.libvc.LibVC; It is worth noting that the combination of LARA and LIBVC can also be used to support compiler flag selection and phase-ordering both statically and dynamically [START_REF] Ashouri | Cobayn: Compiler autotuning framework using bayesian networks[END_REF][START_REF] Ashouri | Micomp: Mitigating the compiler phase-ordering problem using optimization sub-sequences and machine learning[END_REF].

D. Memoization

Optimising applications for energy-efficiency is a challenge of the ANTAREX project. We introduce in this section a memoization technique integrated in the ANTAREX toolflow. Performance can be improved by caching results of pure functions (i.e. deterministic functions without side effects), and retrieving them instead of recomputing a result. We have implemented the work of [START_REF] Suresh | Compile-Time Function Memoization[END_REF] generalized for C++ and aided with extensions regarding user/developer flexibility. We describe here only the principles of this technique and more details can be found in [START_REF] Suresh | Intercepting Functions for Memoization: A Case Study Using Transcendental Functions[END_REF] [START_REF] Suresh | Compile-Time Function Memoization[END_REF]. Consider a memoizable C function foo as shown in Figure 8. The memoization consists in:

1) the insertion of a wrapper function foo_wrapper and an associated table.

2) The substitution of the references to foo by foo_wrapper in the application. This technique has been extended for C++ memoizable methods and takes into account the mangling, the overloading, and the references to the objects. Memoization is proposed in the ANTAREX project by relying on aspects programmed using the DSL. The advantage of these aspects is that the memoization is integrated into the application without requiring user modifications of the source code. The code generated by Clava is then compiled and linked with the associated generated memoization library.

An example of a LARA aspect for memoization is shown in Figure 9. It defines the memoization (lines 1-13) of a method (aMethod) of a class (aClass) with nbArg parameters of same type as the returned type (Type). Note that the inputs nbArg and Type are required to manage the overloading of the object-oriented languages such that C++. Other parameters (from line 4) are provided to improve several memoization approaches. For examples, the user can specify (1) the policy in case of conflicts regarding the same table entry (line 11): replacement or not in case of conflict to the same entry of the table for different parameters of the memoized function, and (2) the size of the table (line 15). After some verifications, not detailed here, on the parameters (lines 14-15), the method is searched (lines [START_REF] Tartara | Continuous learning of compiler heuristics[END_REF][START_REF] Ashouri | Cobayn: Compiler autotuning framework using bayesian networks[END_REF][START_REF] Ashouri | Micomp: Mitigating the compiler phase-ordering problem using optimization sub-sequences and machine learning[END_REF][START_REF] Suresh | Compile-Time Function Memoization[END_REF][START_REF] Suresh | Intercepting Functions for Memoization: A Case Study Using Transcendental Functions[END_REF][START_REF] Gadioli | Application autotuning to support runtime adaptivity in multicore architectures[END_REF][START_REF] Brun | Engineering Self-Adaptive Systems through Feedback Loops[END_REF][START_REF] Gadioli | Socrates -a seamless online compiler and system runtime autotuning framework for energy-aware applications[END_REF]. Then, in case of success, the code of the wrapper is added (line 28) to produce the memoization library, and (line 30) the code of the application is modified for calling the created "wrapper", this wrapper is also declared as a new method of the class.

Moreover, some variables are exposed for autotuning in the memoization library. For each function or method, a variable that manages the dynamical "stop/run" of the memoization is exposed, as well as Fig. 9. An example of LARA aspect defined for the memoization.

the variable that manages the policy to use in case of conflict to the table. To be complete about the memoization, a LARA aspect is proposed to automatically detect the memoizable functions or methods. Then the user may decide to apply or not the memoization on these selected elements.

E. Self-Adaptivity & Autotuning

In ANTAREX, we consider each application's function as a parametric function that elaborates input data to produce an output (i.e., o = f (i, k1, . . . , kn) ), with associated extra-functional requirements. In this context, the parameters of the function (k1, . . . , kn) are software-knobs that modify the behavior of the application (e.g., parallelism level or the number of trials in a MonteCarlo simulation). The main goal of mARGOt4 [START_REF] Gadioli | Application autotuning to support runtime adaptivity in multicore architectures[END_REF] is to enhance an application with an adaptive layer, aiming at tuning the software knobs to satisfy the application requirements at runtime. To achieve this goal, the mARGOt dynamic autotuning framework developed in ANTAREX is based on the MAPE-K feedback loop [START_REF] Brun | Engineering Self-Adaptive Systems through Feedback Loops[END_REF]. In particular, it relies on an application knowledge, derived either at deploy time or at runtime, that states the expected behavior of the extra-functional properties of interest. To adapt, on one hand mARGOt uses runtime observations as feedback information for reacting to the evolution of the execution context. On the other hand, it considers features of the actual input to adapt in a more proactive fashion. Moreover, the framework is designed to be flexible, defining the application requirements as a multi-objective constrained optimization problem that might change at runtime.

To hide the complexity of the application enhancement, we use LARA aspects for configuring mARGOt and for instrumenting the code with related API. Figure 10 provides a simple example of a LARA aspect where mARGOt has been configured (lines 5-20) to actuate on a software knob Knob1 and target error and throughput metrics [START_REF] Gadioli | Socrates -a seamless online compiler and system runtime autotuning framework for energy-aware applications[END_REF]. In particular, the optimization problem has been defined as the maximization of the throughput while keeping the error under a certain threshold. The last part of the aspect (lines 23-27) is devoted to the actual code enhancement including the needed mARGOt call for initializing the framework and for updating the application configuration. The declarative nature of the LARA library developed for integrating mARGOt simplifies its usage hiding all the details of the framework.

F. Monitoring

Today processing elements embed the capability of monitoring their current performance efficiency by inspecting the utilization of the micro-architectural components as well as a set of physical parameters (i.e., power consumption, temperature, etc). These metrics are accessible through hardware performance counters which in x86 systems can be read by privilege users, thus creating practical problems for user-space libraries to access them. Moreover, in addition to sensors which can be read directly from the software running on the core itself, supercomputing machines embed sensors external to the computing elements but relevant to the overall energy-efficiency. These elements include the node and rack cooling components as well as environmental parameters such as the room and ambient temperature. In ANTAREX, we developed ExaMon [START_REF] Beneventi | Continuous learning of hpc infrastructure models using big data analytics and in-memory processing tools[END_REF] (Exascale Monitoring) to virtualise the performance and power monitoring access in a distributed environment. ExaMon decouples the sensor readings from the sensor value usage. Indeed, ExaMon uses a scalable approach were each sensor is associated to a sensing agent which periodically collects the metrics and sends the measured value with a synchronized time-stamp to an external data broker. The data broker organises the incoming data in communication channels with an associated topic. Every new message on a specific topic is then broadcast to the related subscribers, according to a list kept by the broker. The subscriber registers a callback function to the given topic which is called every time a new message is received. To let LARA take advantage of this monitoring mechanism we have designed the Collector API, which allow the initialization of the Collector component associated with a specific topic that keeps an internal state of the remote sensor updated. This internal state can then be queried asynchronously by the Collector API to gather its value. LARA aspects have been designed to embed the Collector API and to make the application code self-aware.

Figure 11 shows a usage example of ExaMon through LARA, which subscribes to a topic on a given broker and inserts a logging message in the application. To define the connection information, the user needs to provide the address to connect to, as well as the name of the topic to subscribe. As for the integration in the original application code the user needs to provide a target function, where the collector will be managed, and a target statement, where the query of the data and logging will be performed.

1 import antarex.examon.Examon; 2 import lara.code.Logger;

G. Power Capping

Today's computing elements and nodes are power limited. For this reason, state-of-the-art processing elements embed the capability of fine-tuning their performance to control dynamically their power consumption. This includes dynamic scaling of voltage and frequency, and power gating for the main architectural blocks of the processing elements, but also some feedback control logic to keep the total power consumption of the processing element within a safe power budget. This logic in x86 systems is named RAPL [START_REF] David | Rapl: Memory power estimation and capping[END_REF]. Demanding the power control of the processing element entirely to RAPL may not be the best choice. Indeed, it has been recently discovered that RAPL is application agnostic and thus tends to waist power in application phases which exhibit IOs or memory slacks. Under these circumstances there are operating points that proved to be more energy efficient than the ones selected by RAPL while still respecting the same power budget [START_REF] Cesarini | Benefits in relaxing the power capping constraint[END_REF]. However, these are only viable if the power capping logic is aware of the application requirements. To do so, we have developed a new power capping run-time based on a set of user space APIs which can be used to define a relative priority for the given task currently in execution on a given core. Thanks to this priority, the run-time is capable of allocating more power to the higher priority process [START_REF] Bartolini | Self-aware thermal management for high performance computing processors[END_REF][START_REF] Cesarini | Prediction horizon vs. efficiency of optimal dynamic thermal control policies in hpc nodes[END_REF]. In ANTAREX, these APIs can be inserted by LARA aspects in the application code.

III. EVALUATION

Tables I and II show static and dynamic metrics collected for the weaving process of the presented examples. In Table I, we can see the number of logical lines of source code for the LARA strategies, as well as for the input code and generated output code (the SLoC-L columns). In the last two columns we report the difference in SLoC and functions between the input and output code (the delta columns). Note the woven and delta results for the HalfPrecisionOpenCL strategy are the sum of all generated code, totaling 31 versions.

An inspection of columns LARA SLoC-L and Delta SLoC-L reveals that, in most examples, there is a large overhead in terms of LARA SLoC-L over application SLoC-L. While this may seem a problem, we need to consider that a large part of the work being performed by these strategies is code analysis, which does not translate directly to SLoC-L in the final application. Furthermore, the Delta SLoC-L metric does not account for removed application code and for these cases a metric based on the similarity degree among code versions could be of more interest. Also, in real-world applications, the ratio of LARA SLoC-L to application SLoC-L would be definitely more favorable, thanks to aspect reuse.

To better understand the impact of analysis, we report in the first two columns of Table II the number of code points and of their attributes analysed, which can be compared with the last three column of the same table, which instead report the corresponding effects, in terms of the number of modified points and lines of code inserted. To understand the impact of removed lines of code, we look at the Inserts and Actions columns, which show that circa one half of the actions do not insert code. The end line is that the analysis work exceeds the transformation work by an order of magnitude, and the insertions only underestimate significantly the work performed.

Another benefit for user productivity when using LARA is how the techniques presented in the examples can scale into large-scale applications and scenarios. Most of the presented strategies are parameterized by function, i.e., they receive a function join point or name and act on the corresponding function. This could be performed manually, albeit crudely, using a search function of an IDE. Consider the case where we instead want to target a set of functions, whose names we may not know, based on their function signatures, or based on the characteristics of the variables declared inside their scope. This kind of search and filtering based on syntactic and semantic information available in the program is one of the key features of LARA and it cannot be easily attained with other tools. As the aspects presented here illustrate, LARA strategies can be made reusable and applied over large applications, greatly out scaling the effort needed to develop them.

IV. CASE STUDY: SELF-ADAPTIVE NAVIGATION SYSTEM

In this section, we provide an overview of the application of the ANTAREX tool flow to the Self-Adaptive Navigation System developed in Use Case 2. The system is designed to process large volumes of data for the global view computation and to handle dynamic loads represented by incoming routing requests from users of the system. Both disciplines require HPC infrastructure in order to operate efficiently while maintaining contracted SLA. Integration of the ANTAREX self-adaptive holistic approach can help the system to meet the mentioned requirements and pave the way to scaling its operation to future Exascale systems. Core of the system is a routing pipeline with several stages which uses our custom algorithm library written in C++. The library provides an API for the individual routing algorithms and for data access layer, which provides abstraction of a graph representation of the road network. The graph is stored in a HDF5 file, which is a well known and convenient storage format for structured data on HPC clusters.

As an example, we are using LARA aspects to generate C++ code for mapping native data types to types defined by the HDF5 API. The aspects are applied using the Clava tool which is a C++ frontend for the LARA toolchain. The Clava tool is integrated in our CMakebased build process as a custom build step, which parses the C++ structures representing the routing graph in memory and produces part of the HDF5 data access API. Details of the implementation can be found in [START_REF] Golasowski | Expressing and applying c++ code transformations for the hdf5 api through a dsl[END_REF]. Using the same process, other LARA aspects can be easily applied on the source code of the library, which greatly simplifies integration of other tools of the ANTAREX toolchain, such as mARGOt [START_REF] Gadioli | Application autotuning to support runtime adaptivity in multicore architectures[END_REF].

Furthermore, the mARGOt [START_REF] Gadioli | Application autotuning to support runtime adaptivity in multicore architectures[END_REF] autotuner is used in the Probabilistic Time-Dependent routing (PTDR) algorithm [START_REF] Golasowski | Performance evaluation of probabilistic time-dependent travel time computation[END_REF] to dynamically adjust the number of Monte Carlo samples used for the particular routing request. This parameter directly affects load generated by the PTDR stage and precision of its output. The autotuner uses operation point lists generated by a Design Space Exploration phase. The operation points in the context of PTDR are represented by a number of MC samples as an adjustable algorithm parameter and expected values of various metrics. The autotuner then dynamically selects the operation point according to the current request input. This approach can significantly reduce computational load generated by the PTDR phase, contributing to the overall efficient operation of the system.

Currently, our codebase is ready to use the DSL to integrate other tools from the ANTAREX tool flow. The autotuner is manually integrated in the routing pipeline, while verification of its correct function is ongoing. The next step is to use LARA to integrate the autotuner to the target application and evaluate its impact.

We have developed a server-side routing dashboard web application which is used to monitor the current status of the routing service. The application also provides a consistent environment for testing the service performance. It provides a way to execute a benchmark of the service by adjusting its parameters and sending a pre-defined set of routing requests. The service performance is then measured and results of the testing are stored for further analysis. The application also provides a consistent visualisation of the results which can be used for further analysis. This infrastructure will be used for validation of the ANTAREX tools integrated in the routing service [START_REF] Zadnik Stirn | Server-side navigation service benchmarking tool[END_REF].

V. CONCLUSIONS

To fully exploit the heterogeneous resources of future Exascale HPC systems, new software stacks are needed to provide power management, optimization, and autotuning to the parallel applications deployed on such systems. The ANTAREX project provides a holistic system-wide adaptive approach for next generation HPC systems, centered around a domain specific language that allows a full decoupling of functional and extra-functional specifications for each application, providing integration with a wide range of support tools. We have shown how the ANTAREX tool flow allows developers to control the precision of a computation, to manage dynamic code specialization, monitoring, power capping, and dynamic autotuning. The impact and benefits of such technology are far reaching, beyond traditional HPC domains.

Fig. 1 .

 1 Fig. 1. The ANTAREX Tool Flow

1 aspectdef ChangePrecision 2 3 5 / 9 varFig. 2 .

 12592 Fig. 2. Example of LARA aspect to change the types of variables declared inside a given function.

17 /Fig. 3 .

 173 Fig. 3. Example of LARA aspect that generates different precision mix versions of the same OpenCL kernel.

3 4 aspectdef CreateFloatVersion 5 input $func, suffix end 6 output $clonedFunc end 7 8 $ 11 / 15 / 19 $Fig. 4 .

 5681115194 Fig. 4. Example of LARA aspect to clone an existing function and change the type of the clone.

4 5 aspectdef Multiversion 6 input $func, knobName end 7 8 12 / 15 / 23 / 30 /Fig. 5 .

 67121523305 Fig. 5. Example of LARA aspect that generates an alternative version of a function and inserts a mechanism in the code to switch between versions.

2 caseFig. 6 .

 26 Fig. 6. Excerpt of the C code resulting from the generation of alternative code versions.

var 12 var

 12 $function = $target.definition; 10 var lvc = new LibVC($function, {logFile:"log.txt"}, name); 11 lvcOptions = new LibVCOptions(); 13 for (var o of options) { 14 lvcOptions.addOptionLiteral(o.name, o.value, o.value);

Fig. 7 .

 7 Fig. 7. Example of LARA aspect to replace a function call to a kernel with a call to a dynamically generated version of that kernel.

Fig. 8 .

 8 Fig. 8. A memoizable C function and its wrapper.

1 aspectdef mARGOt_Aspect 2 /

 12 * Input: TargetFunctionCall * / 3 input targetCallName end 4 / * mARGOt configuration * / 5 var config = new MargotConfig(); 6 var targetBlock = config.newBlock($targetCallName);

7 13 / 18 /Fig. 10 .

 131810 Fig. 10. Example of a LARA aspect for autotuner configuration and code enhancement.

var 14 /

 14 broker = new ExamonBroker(ip); 11 12 var exa = new ExamonCollector(name, topic); 13 / manage the collector on the target function 15 select $manageFunction.body end 16 apply 17 exa.init(broker, $body); 18 exa.start($body);

27 /Fig. 11 .

 2711 Fig. 11. Example of a LARA aspect integrate an ExaMon collector into an application.

  FileToSave, // filemane to save the table, or 'none'

	1 aspectdef Memoize_Method_overloading_ARGS
	2 input
	3	aClass,	// Name of a class
	4	aMethod,	// Name of a method of the class aClass
	5	pType,	// Name of the selected type
	6	nbArgs,	// Number of parameters of the method
	10	Replace,	// Always replace in case of collisions
	11	approx,	// Number of bits to delete for approximation.
	12	tsize	// Size of the internal table.
	13 end
	14 // Control on the parameters of the aspect: nbArgs in [1,3]
	15	...
	16 // Searching the method.
	17 var MethodToMemoize, found=false;
	18 select class{aClass}.method{aMethod} end
	19 apply
	20 if (! found) {
	21	found = isTheSelectedMethod($method, nbArgs, pType);
	22	if (found) MethodToMemoize=$method;
	23 }
	24 end
	25 if (!found)
	26 { / * message to the user * /}
	27 else {
	28 GenCode_CPP_Memoization(aClass, aMethod, pType, nbArgs,
	29 fileToLoad, FullOffLine, FileToSave, Replace, approx, tsize);
	30 call CPP_UpdateCallMemoization(aClass, aMethod, pType, nbArgs);
	31 }
	32 end

7 fileToLoad, // filename for init of the table, or 'none' 8 FullOffLine,// yes for a fully offline strategy 9

https://github.com/specs-feup/clava

Complete working versions for all examples can be found in https://github.com/specs-feup/specs-lara/tree/master/2018%20DSD

https://gitlab.com/margot_project/core

ACKNOWLEDGEMENTS

The ANTAREX project is supported by the EU H2020 FET-HPC program under grant 671623.


Cédric Colas

Olivier Sigaud

Pierre-Yves Oudeyer

How Many Random Seeds? Statistical Power Analysis in Deep Reinforcement Learning Experiments

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Perhaps the most surprising thing is this: running the same algorithm 10 times with the same hyper-parameters using 10 different random seeds and averaging performance over two splits of 5 seeds can lead to learning curves seemingly coming from different statistical distributions. Notably, all the deep RL papers reviewed by Henderson et al.. (theirs included) used 5 seeds or less. Even worse, some papers actually report the average of the best performing runs. As demonstrated in [START_REF] Henderson | Deep reinforcement learning that matters[END_REF], these methodologies can lead to claim that two algorithms performances are different when they are not. A solution to this problem is to use more random seeds, to average more different trials in order to obtain a more robust measure of the algorithm performance. But how can one determine how many random seeds should be used? Shall we use 5, 10 or 100, as in [START_REF] Mania | Simple random search provides a competitive approach to reinforcement learning[END_REF]?

This work assumes one wants to test a difference in performance between two algorithms. Section 2 gives definitions and describes the statistical problem of difference testing while Section 3 proposes two statistical tests to answer this problem. In Section 4, we present standard guidelines to choose the sample size so as to meet requirements in the two types of statistical errors. Finally, we challenge the assumptions made in the previous section and propose guidelines to estimate error rates empirically in Section 5. The code is available on Github at https://github.com/flowersteam/rl-difference-testing. We might consider that Algo1 outperforms Algo2 because there is not much overlap between the 95% confidence intervals. But is it sufficient evidence that Algo1 really performs better? Below, we show that the performances of these algorithms is actually the same, and explain which methods should be used to have more reliable evidence of the (non-)difference among two algorithms.

2 Definition of the statistical problem

First definitions

Two runs of the same algorithm often yield different measures of performance. This might be due to various factors such as the seed of the random generators (called random seed or seed thereafter), the initial state of the agent, the stochasticity of the environment, etc. Formally, the performance of an algorithm can be modeled as a random variable X and running this algorithm in an environment results in a realization x i . Repeating the procedure N times, one obtains a statistical sample x = (x 1 , .., x N ). A random variable is usually characterized by its expected value or mean µ and its standard deviation, noted σ. While the mean characterizes the expected value of a realization, the standard deviation evaluates the square root of the squared deviations to this mean, or in simpler words, how far from the mean the realization are expected to fall. Of course, the values of µ and σ are unknown. The only thing one can do is to compute their unbiased estimations x and s:

x = n i=1 x i , s = N i+1 (x i -x) 2 N -1 , (1)
where x is called the empirical mean, and s is called the empirical standard deviation. The larger the sample size N , the more confidence one can be in the estimations.

Here, two algorithms with respective performances X 1 and X 2 are compared. If X 1 and X 2 follow normal distributions, the random variable describing their difference (X diff = X 1 -X 2 ) also follows a normal distribution with parameters σ diff = (σ 2 1 + σ 2 2 ) 1/2 and µ diff = µ 1 -µ 2 . In this case, the estimator of the mean of X diff is x diff = x 1 -x 2 and the estimator of σ diff is s diff = s 2 1 + s 2 2 . The effect size can be defined as the difference between the mean performances of both algorithms: = µ 1 -µ 2 .

Testing for a difference between the performances of two algorithms (µ 1 and µ 2 ) is mathematically equivalent to testing a difference between their difference µ diff and 0. The second point of view is considered from now on. We draw a sample x diff from X diff by subtracting two samples x 1 and x 2 obtained from X 1 and X 2 .

Example 1. To illustrate difference testing, we use two algorithms (Algo1 and Algo2) and compare them on the Half-Cheetah environment from the OpenAI Gym framework [START_REF] Brockman | Openai gym[END_REF]. The algorithms implemented are not so important here and will be revealed later. First, we run a study with N = 5 random seeds for each. Figure 1 shows the average learning curves with 95% confidence intervals. Each point of a learning curve is the average cumulated reward over 10 evaluation episodes. The measure of performance X i of Algo i is the average performance over the last 10 points (i.e. last 100 evaluation episodes). From Figure 1, it seems that Algo1 performs better than Algo2. Moreover, the confidence intervals do not seem to overlap much. However, we need to run statistical tests before drawing any conclusion.

Comparing performances with a difference test

In a difference test, statisticians define the null hypothesis H 0 and the alternate hypothesis H a . H 0 assumes no difference whereas H a assumes one:

• H 0 : µ diff = 0 • H a : µ diff = 0
These hypotheses refer to the two-tail case. When an a priori on which algorithm performs best is available, (say Algo1), one can use the one-tail version:

• H 0 : µ diff ≤ 0 • H a : µ diff > 0
At first, a statistical test always assumes the null hypothesis. Once a sample x diff is collected from X diff , one can estimate the probability p (called p-value) of observing data as extreme, under the null hypothesis assumption. By extreme, one means far from the null hypothesis (x diff far from 0). The p-value answers the following question: how probable is it to observe this sample or a more extreme one, given that there is no true difference in the performances of both algorithms? Mathematically, we can write it this way for the one-tail case:

p-value = P (X diff ≥ x diff | H 0 ), (2)
and this way for the two-tail case:

p-value = P (X diff ≥ x diff | H 0 ) if x diff > 0 P (X diff ≤ x diff | H 0 ) if x diff ≤ 0. ( 3
)
When this probability becomes really low, it means that it is highly improbable that two algorithms with no performance difference produced the collected sample x diff . A difference is called significant at significance level α when the p-value is lower than α in the one-tail case, and lower than α/2 in the two tail case (to account for the two-sided test2 ). Usually α is set to 0.05 or lower. In this case, the low probability to observe the collected sample under hypothesis H 0 results in its rejection. Note that a significance level α = 0.05 still results in 1 chance out of 20 to claim a false positive, to claim that there is a true difference when there is not. It is important to note that, when one is conducting N E experiments, the false positive rate grows linearly with the number of experiments. In this case, one should use correction for multiple comparisons such as the Bonferroni correction α Bon = α/N E [START_REF] Rice | Analyzing tables of statistical tests[END_REF]. This controls the familywise error rate (F W ER), the probability of rejecting at least one true null hypothesis (F W ER < α). Its use is discussed in [START_REF] Cabin | To bonferroni or not to bonferroni: when and how are the questions[END_REF].

Another way to see this, is to consider confidence intervals. Two kinds of confidence intervals can be computed:

• CI 1 : The 100 • (1 -α) % confidence interval for the mean of the difference µ diff given a sample x diff characterized by x diff and s diff .

• CI 2 : The 100 • (1 -α) % confidence interval for any realization of X diff under H 0 (assuming µ diff = 0).

Having CI 2 that does not include x diff is mathematically equivalent to a p-value below α. In both cases, it means there is less than 100 • α% chance that µ diff = 0 under H 0 . When CI 1 does not include 0, we are also 100 • (1 -α) % confident that µ = 0, without assuming H 0 . Proving one of these things leads to conclude that the difference is significant at level α.

Statistical errors

In hypothesis testing, the statistical test can conclude H 0 or H a while each of them can be either true or false. There are four cases:

False positive α

True positive 1 -β

This leads to two types of errors:

• The type-I error rejects H 0 when it is true, also called false positive. This corresponds to claiming the superiority of an algorithm over another when there is no true difference. Note that we call both the significance level and the probability of type-I error α because they both refer to the same concept. Choosing a significance level of α enforces a probability of type-I error α, under the assumptions of the statistical test.

• The type-II error fails to reject H 0 when it is false, also called false negative. This corresponds to missing the opportunity to publish an article when there was actually something to be found.

Message 1:

• In the two-tail case, the null hypothesis

H 0 is µ diff = 0. The alternative hypothesis H a is µ diff = 0. • p-value = P (X diff ≥ x diff | H 0 ).
• A difference is said statistically significant when a statistical test passed. One can reject the null hypothesis when 1) p-value< α; 2) CI 1 does not contain 0; 3) CI 2 does not contain x diff .

• statistically significant does not refer to the absolute truth. Two types of error can occur. Type-I error rejects H 0 when it is true. Type-II error fails to reject H 0 when it is false.

• The rate of false positive is 1 out of 20 for α = 0.05. It grows linearly with the number of experiment N E . Correction procedures can be applied to correct for multiple comparisons.

Choice of the appropriate statistical test

In statistics, a difference cannot be proven with 100% confidence. To show evidence for a difference, we use statistical tests. All statistical tests make assumptions that allow them to evaluate either the p-value or one of the confidence intervals described in the Section 2. The probability of the two error types must be constrained, so that the statistical test produces reliable conclusions. In this section we present two statistical tests for difference testing. As recommended in [START_REF] Henderson | Deep reinforcement learning that matters[END_REF], the two-sample t-test and the bootstrap confidence interval test can be used for this purpose3 .

T-test and Welch's t-test

We want to test the hypothesis that two populations have equal means (null hypothesis H 0 ). A 2-sample t-test can be used when the variances of both populations (both algorithms) are assumed equal. However, this assumption rarely holds when comparing two different algorithms (e.g. DDPG vs TRPO).

In this case, an adaptation of the 2-sample t-test for unequal variances called Welch's t-test should be used [START_REF] Welch | The generalization ofstudent's' problem when several different population variances are involved[END_REF]. Both tests are strictly equivalent when the standard deviations are equal. T -tests make a few assumptions:

• The scale of data measurements must be continuous and ordinal (can be ranked). This is the case in RL.

• Data is obtained by collecting a representative sample from the population. This seem reasonable in RL.

• Measurements are independent from one another. This seems reasonable in RL.

• Data is normally-distributed, or at least bell-shaped. The normal law being a mathematical concept involving infinity, nothing is ever perfectly normally distributed. Moreover, measurements of algorithm performances might follow multi-modal distributions. In Section, 5, we investigate the effects of deviations from normality.

Under these assumptions, one can compute the t-statistic t and the degree of freedom ν for the Welch's t-test as estimated by the Welch-Satterthwaite equation, such as:

t = x diff s 2 1 N1 + s 2 2 N2 , ν ≈ s 2 1 N1 + s 2 2 N2 2 s 4 1 N 2 1 (N1-1) + s 4 2 N 2 2 (N2-1) , (4)
with x diff = x 1 -x 2 ; s 1 , s 2 the empirical standard deviations of the two samples, and N 1 , N 2 their sizes. Sample sizes are assumed equal (N 1 = N 2 = N ) thereafter. The t-statistics are assumed to follow a t-distribution, which is bell-shaped and whose width depends on the degree of freedom. The higher this degree, the thinner the distribution.

Figure 2 helps making sense of these concepts. It represents the distribution of the t-statistics corresponding to X diff , under H 0 (left distribution) and under H a (right distribution). H 0 assumes µ diff = 0, the distribution is therefore centered on 0. H a assumes a (positive) difference µ diff = , the distribution is therefore shifted by the t-value corresponding to , t . Note that we consider the one-tail case here, and test for a positive difference.

A t-distribution is defined by its probability density function T ν distrib (τ ) (left curve in Figure 2), which is parameterized by ν. The cumulative distribution function CDF H0 (t) is the function evaluating the area under T ν distrib (t) from τ = -∞ to τ = t. This allows to write: In Figure 2, t α represents the critical t-value to satisfy the significance level α in the one-tail case. When t = t α , p-value= α. When t > t α , the p-value is lower than α and the test rejects H 0 . On the other hand, when t is lower than t α , the p-value is superior to α and the test fails to reject H 0 . As can be seen in the figure, setting the threshold at t α might also cause an error of type-II. The rate of this error (β) is represented by the dark blue area: under the hypothesis of a true difference (under H a , right distribution), we fail to reject H 0 when t is inferior to t α . β can therefore be computed mathematically using the CDF :

p-value = 1 -CDF H0 (t) = 1 - t -∞ T ν distrib (τ ) • dτ. (5)
β = CDF Ha (t α ) = tα -∞ T ν distrib (τ -t ) • dτ. (6)
Using the translation properties of integrals, we can rewrite β as:

β = CDF H0 (t α -t ) = tα-t -∞-t =-∞ T ν distrib (τ ) • dτ. (7)
The procedure to run a Welch's t-test given two samples (x 1 , x 2 ) is:

• Computing the degree of freedom ν and the t-statistic t based on s 1 , s 2 , N and x diff .

• Looking up the t α value for the degree of freedom ν in a t-table 4 or by evaluating the inverse of the CDF function in α.

• Compare the t-statistic to t α . The difference is said statistically significant (H 0 rejected) at level α when t ≥ t α .

Note that t < t α does not mean there is no difference between the performances of both algorithms. It only means there is not enough evidence to prove its existence with 100 • (1 -α)% confidence (it might be a type-II error). Noise might hinder the ability of the test to detect the difference. In this case, increasing the sample size N could help uncover the difference.

Selecting the significance level α of the t-test enforces the probability of type-I error to α. However, Figure 2 shows that decreasing this probability boils down to increasing t α , which in turn increases the probability of type-II error β. One can decrease β while keeping α constant by increasing the sample size N . This way, the estimation x diff of µ diff gets more accurate, which translates in thinner distributions in the figure, resulting in a smaller β. The next section gives standard guidelines to select N so as to meet requirements for both α and β.

Bootstrapped confidence intervals

Bootstrapped confidence interval is a method that does not make any assumption on the distribution of X diff . It estimates the confidence interval CI 1 for µ diff , given a sample x diff characterized by its empirical mean x diff . It is done by re-sampling inside x diff and by computing the mean of each newly generated sample. The test makes its decision based on whether the confidence interval of x diff contains 0 or not. It does not compute a p-value as such.

Without any assumption on the data distribution, an analytical confidence interval cannot be computed. Here, X diff follows an unknown distribution F . An estimation of the confidence interval CI 1 can be computed using the bootstrap principle.

Let us say we have a sample x diff made of N measures of performance difference. The empirical bootstrap sample x * diff of size N is obtained by sampling with replacement inside x diff . The bootstrap principle then says that, for any statistic u computed on the original sample and u * computed on the bootstrap sample, variations in u are well approximated by variations in u * 5 . Therefore, variations of the empirical mean such as its range can be approximated by variations of the bootstrapped samples. The bootstrap confidence interval test assumes the sample size is large enough to represent the underlying distribution correctly, although this might be difficult to achieve in practice. Deviations from this assumption are discussed in Section 5. Under this assumption, the bootstrap test procedure looks like this:

• Generate B bootstrap samples of size N from the original sample x 1 of Algo1 and B samples from from the original sample x 2 of Algo2.

• Compute the empirical mean for each sample:

µ 1 1 , µ 2 1 , ..., µ B 1 and µ 1 2 , µ 2 2 , ..., µ B 2 • Compute the differences µ 1:B diff = µ 1:B 1 -µ 1:B 2
• Compute the bootstrapped confidence interval at 100 • (1 -α)%. This is basically the range between the 100 • α/2 and 100 • (1 -α)/2 percentiles of the vector µ 1:B diff (e.g. for α = 0.05, the range between the 2.5 th and the 97.5 th percentiles).

The number of bootstrap B should be chosen large (e.g. > 1000). If the confidence interval does not contain 0, it means that one can be confident at 100 • (1 -α)% that the difference is either positive (both bounds positive) or negative (both bounds negative), thus, that there is a statistically significant difference between the performances of both algorithms6 .

Example 1 (continued). Here, the type-I error requirement is set to α = 0.05. Running the Welch's t-test and the bootstrap confidence interval test with two samples (x 1 , x 2 ) of 5 seeds each leads to a pvalue of 0.031 and a bootstrap confidence interval such that P µ diff ∈ [259, 1564] = 0.05. Since the p-value is below the significance level α and the CI 1 confidence interval does not include 0, both test passed. This means both tests found a significant difference between the performances of Algo1 and Algo2 with a 95% confidence. There should have been only 5% chance to conclude a significant difference if it did not exist.

In fact, we did encounter a type-I error. We know this for sure because Algo 1 and Algo 2 were the exact same algorithm. They are both the canonical implementation of DDPG [START_REF] Lillicrap | Continuous control with deep reinforcement learning[END_REF] from the OpenAI baselines [START_REF] Dhariwal | Openai baselines[END_REF]. The first conclusion was wrong, we committed a type-I error, rejecting H 0 when it was true. We knew this could happen with probability α = 0.05. Section 5 shows that this probability might have been under-evaluated because of the assumptions made by the statistical tests.

Message 2:

• T -tests assume t-distributions of the t-values. Under some assumptions, they can compute analytically the p-value and the confidence interval CI 2 at level α.

• The Welch's t-test does not assume both algorithms have equal variances but the t-test does.

• The bootstrapped confidence interval test does not make assumptions on the performance distribution and estimates empirically the confidence interval CI 1 at level

• Selecting a test with a significance level α enforces a type-I error α when the assumptions of the test are verified.

4 In theory: power analysis for the choice of the sample size

In the Section 3, we saw that α was enforced by the choice of the significance level in the test implementation. The second type of error β must now be estimated. β is the probability to fail to reject H 0 when H a is true. When the effect size and the probability of type-I error α are kept constant, β is a function of the sample size N . Choosing N so as to meet requirements on β is called statistical power analysis. It answers the question: what sample size do I need to have 1 -β chance to detect an effect size , using a test with significance level α? The next paragraphs present guidelines to choose N in the context of a Welch's t-test.

As we saw in Section 3.1, β can be analytically computed as:

β = CDF H0 (t α -t ) = tα-t -∞-t =-∞ T ν distrib (τ ) • dτ, (8)
where CDF H0 is the cumulative distribution function of a t-distribution centered on 0, t α is the critical value for significance level α and t is the t-value corresponding to an effect size . In the end, β depends on α, , (s 1 , s 2 ) the empirical standard deviations computed on two samples (x 1 , x 2 ) and the sample size N .

Example 2. To illustrate, we compare two DDPG variants: one with action perturbations (Algo1) [START_REF] Lillicrap | Continuous control with deep reinforcement learning[END_REF], the other with parameter perturbations (Algo2) [START_REF] Plappert | Parameter space noise for exploration[END_REF]. Both algorithms are evaluated in the Half-Cheetah environment from the OpenAI Gym framework [START_REF] Brockman | Openai gym[END_REF].

Step 1 -Running a pilot study

To compute β, we need estimates of the standard deviations of the two algorithms (s 1 , s 2 ). In this step, the algorithms are run in the environment to gather two samples x 1 and x 2 of size n. From there, we can compute the empirical means (x 1 , x 2 ) and standard deviations (s 1 , s 2 ).

Example 2 (continued). Here we run both algorithms with n = 5. We find empirical means (x 1 , x 2 ) = (3523, 4905) and empirical standard deviations (s 1 , s 2 ) = (1341, 990) for Algo1 (blue) and Algo2 (red) respectively. From Figure 3, it seems there is a slight difference in the mean performances

x diff = x 2 -x 1 > 0.
Running preliminary statistical tests at level α = 0.05 lead to a p-value of 0.1 for the Welch's t-test, and a bootstrapped confidence interval of CI 1 = [795, 2692] for the value of x diff = 1382. The Welch's t-test does not reject H 0 (p-value> α) but the bootstrap test does (0 ∈ CI 1 ). One should compute β to estimate the chance that the Welch's t-test missed an underlying performance difference (type-II error).

Step 2 -Choosing the sample size

Given a statistical test (Welch's t-test), a significance level α (e.g. α = 0.05) and empirical estimations of the standard deviations of Algo1 and Algo2 (s 1 , s 2 ), one can compute β as a function of the sample size N and the effect size one wants to be able to detect.

Example 2 (continued). For N in [2,50] and in [0.1, .., 1] × x 1 , we compute t α and ν using the formulas given in Section 3.1, as well as t for each . Finally, we compute the corresponding probability of type-II error β using Equation 8. Figure 4 shows the evolution of β as a function of N for the different . Considering the semi-dashed black line for = x diff = 1382, we find β = 0.51 for N = 5: there is 51% chance of making a type-II error when trying to detect an effect = 1382. To meet the requirement β = 0.2, N should be increased to N = 10 (β = 0.19). In our example, we find that N = 10 was enough to be able to detect an effect size = 1382 with a Welch's t-test, using significance level α and using empirical estimations (s 1 , s 2 ) = (1341, 990). However, let us keep in mind that these computations use various approximations (ν, s 1 , s 2 ) and make assumptions about the shape of the t-values distribution. Section 5 investigates the influence of these assumptions.

Step 3 -Running the statistical tests

Both algorithms should be run so as to obtain a sample x diff of size N . The statistical tests can be applied.

Example 2 (continued). Here, we take N = 10 and run both the Welch's t-test and the bootstrap test. We now find empirical means (x 1 , x 2 ) = (3690, 5323) and empirical standard deviations (s 1 , s 2 ) = (1086, 1454) for Algo1 and Algo2 respectively. Both tests rejected H 0 , with a p-value of 0.0037 for the Welch's t-test and a confidence interval for the difference µ diff ∈ [732, 2612] for the bootstrap test. Both tests passed. In Figure 5, plots for N = 5 and N = 10 can be compared. With a larger number of seeds, the difference that was not found significant with N = 5 is now more clearly visible. With a larger number of seeds, the estimate x diff is more robust, more evidence is available to support the claim that Algo2 outperforms Algo1, which translates to tighter confidence intervals represented in the figures. Message 3: Given a sample size N , a minimum effect size to detect and a requirement on type-I error α the probability of type-II error β can be computed. This computation relies on the assumptions of the t-test. The sample size N should be chosen so as to meet the requirements on β.

In practice: influence of deviations from assumptions

Under their respective assumptions, the t-test and bootstrap test enforce the probability of type-I error to the selected significance level α. These assumptions should be carefully checked, if one wants to report the probability of errors accurately. First, we propose to compute an empirical evaluation of the type-I error based on experimental data, and show that: 1) the bootstrap test is sensitive to small sample sizes; 2) the t-test might slightly under-evaluate the type-I error for non-normal data. Second, we show that inaccuracies in the estimation of the empirical standard deviations s 1 and s 2 due to low sample size might lead to large errors in the computation of β, which in turn leads to under-estimate the sample size required for the experiment.

Empirical estimation of the type-I error

Remember, type-I errors occur when the null hypothesis (H 0 ) is rejected in favor of the alternative hypothesis (H a ), H 0 being correct. Given the sample size N , the probability of type-I error can be estimated as follows:

• Run twice this number of trials (2 × N ) for a given algorithm. This ensures that H 0 is true because all measurements come from the same distribution.

• Get average performance over two randomly drawn splits of size N . Consider both splits as samples coming from two different algorithms.

• Test for the difference of both fictive algorithms and record the outcome.

• Repeat this procedure T times (e.g. T = 1000)

• Compute the proportion of time H 0 was rejected. This is the empirical evaluation of α.

Example 3 We use Algo1 from Example 2. From 42 available measures of performance, the above procedure is run for N in [2,21]. Figure 6 presents the results. For small values of N , empirical estimations of the false positive rate are much larger than the supposedly enforced value α = 0.05. In our experiment, the bootstrap confidence interval test should not be used with small sample sizes (< 10). Even in this case, the probability of type-I error (≈ 10%) is under-evaluated by the test (5%). The Welch's t-test controls for this effect, because the test is much harder to pass when N is small (due to the increase of t α ). However, the true (empirical) false positive rate might still be slightly under-evaluated. In this case, we might want to set the significance level to α < 0.05 to make sure the true positive rate stays below 0.05. In the bootstrap test, the error is due to the inability of small samples to correctly represent the underlying distribution, which impairs the enforcement of the false positive rate to the significance level α. Concerning the Welch's t-test, this might be due to the non-normality of our data (whose histogram seems to reveal a bimodal distribution). In Example 1, we used N = 5 and encountered a type-I error. We can see on the Figure 6 that the probability of this to happen was around 10% for the bootstrap test and above 5% for the Welch's t-test.

Influence of the empirical standard deviations

The Welch's t-test computes t-statistics and the degree of freedom ν based on the sample size N and the empirical estimations of standard deviations s 1 and s 2 . When N is low, estimations s 1 and s 2 under-estimate the true standard deviation in average. Under-estimating (s 1 , s 2 ) leads to smaller ν and lower t α , which in turn leads to lower estimations of β. Finally, finding lower β leads to the selection of smaller sample size N to meet β requirements. Let us investigate how big this effect can be. In Figure 7, one estimates the standard deviation of a normally distributed variable N (0, 1). The empirical estimation s is quite variable and underestimates σ = 1 in average. We consider estimations of the false negative rate as a function of N when comparing two normal distributions (σ = 1), one centered on 3, the other on 3 + . When we select n = 5 for a preliminary study and compute estimations (s 1 , s 2 ) from this sample, our average error is mean(s n=5 )= -0.059 (see above Figure 7). One could also make larger errors: mean(s n=5 )-std(s n=5 )= -0.40 from the same figure.

Figure 8 shows the effect of an error of 0.40 on the evaluation of β. One can see that, if we want to detect an effect size = 0.9 (green curve) and meet a requirement β = 0.2, one would choose N = 17 when standard deviations are correctly estimated (left) and N = 7 when they are under-evaluated. When the number of samples n available in the preliminary study to compute (s 1 , s 2 ) grows, the under-estimation reduces in average and in the worst case. This, in turn, reduces the inaccuracy in the estimation of β and therefore in the required N . Another solution is to systematically choose the sample size larger than what is prescribed by the computation of β. Message 4:

• One should not blindly believe in statistical tests results. These tests are based on assumptions that are not always reasonable.

• α must be empirically estimated, as the statistical tests might underestimate it, because of wrong assumptions about the underlying distributions or because of the small sample size.

• The bootstrap test evaluation of type-I error is strongly dependent on the sample size. A bootstrap test should not be used with less than 20 samples.

• The inaccuracies in the estimation of the standard deviations of the algorithms (s 1 , s 2 ), due to small sample sizes n in the preliminary study, lead to under-estimate the sample size N required to meet requirements in type-II errors.

Figure 1 :

 1 Figure 1: Algo1 versus Algo2 are two famous Deep RL algorithms, here tested on the Half-Cheetah benchmark. The mean and confidence interval for 5 seeds are reported.We might consider that Algo1 outperforms Algo2 because there is not much overlap between the 95% confidence intervals. But is it sufficient evidence that Algo1 really performs better? Below, we show that the performances of these algorithms is actually the same, and explain which methods should be used to have more reliable evidence of the (non-)difference among two algorithms.

Figure 2 :

 2 Figure 2: Representation of H0 and Ha under the t-test assumptions. Areas under the distributions represented in red, dark blue and light blue correspond to the probability of type-I error α, type-II error β and the statistical power 1 -β respectively.

Figure 3 :

 3 Figure 3: DDPG with action perturbation versus DDPG with parameter perturbation tested in Half-Cheetah. Mean and 95% confidence interval computed over 5 seeds are reported. The figure shows a small difference in the empirical mean performances.

Figure 4 :

 4 Figure 4: Evolution of the probability of type-II error as a function of the sample size N for various effect sizes , when (s1, s2) = (1341, 990) and α = 0.05. The requirement 0.2 is represented by the horizontal dashed black line. The curve for = x diff is represented by the semi-dashed black line.

Figure 5 :

 5 Figure 5: Performance of DDPG with action perturbation (Algo1) and parameter perturbation (Algo2) with N = 5 seeds (left) and N = 10 seeds (right). The 95% confidence intervals on the right are smaller, because more evidence is available (N larger). The underlying difference appears when N grows.

Figure 6 :

 6 Figure 6: Empirical estimations of the false positive rate on experimental data (Example 3) when N varies, using the Welch's t-test (blue) and the bootstrap confidence interval test (orange).

Figure 7 :

 7 Figure 7: Empirical standard deviation of X ∼ N (0, 1). The true standard deviation σ = 1 is represented in red. Mean +/-std are shown.

Figure 8 :

 8 Figure 8: Evolution of the probability of type-II error as a function of the sample size N and the effect size , when (s1, s2) = (1 -error, 1 -error) and α = 0.05. Left: error = 0, this is the ideal case. Right: error = 0.40, a large error that can be made when evaluating s over n = 5 samples. The compared distributions are normal, one is centered on 3, the other on 3 + .

Table 1 :

 1 Hypothesis testing

	predicted/true	H 0	H a
	H 0	True negative 1 -α	False negative β
	H a

Reproducibility in MachineLearning and Deep Reinforcement Learning (RL) in particular has become a serious issue in the recent years. As pointed out in[START_REF] Islam | Reproducibility of benchmarked deep reinforcement learning tasks for continuous control[END_REF] and[START_REF] Henderson | Deep reinforcement learning that matters[END_REF], reproducing the results of an RL paper can turn out to be much more complicated than expected. Indeed, codebases are not always released and scientific papers often omit parts of the implementation tricks. Recently, Henderson et al. conducted a thorough investigation of various parameters causing this reproducibility crisis. They used trendy deep RL algorithms such as DDPG[START_REF] Lillicrap | Continuous control with deep reinforcement learning[END_REF], ACKTR[START_REF] Wu | Scalable trust-region method for deep reinforcement learning using kronecker-factored approximation[END_REF], TRPO[START_REF] Schulman | Trust region policy optimization[END_REF] and PPO[START_REF] Schulman | Proximal policy optimization algorithms[END_REF] with OpenAI Gym[START_REF] Brockman | Openai gym[END_REF] popular benchmarks such as Half-Cheetah, Hopper and Swimmer, to study the effects of the codebase, the size of the networks, the activation function, the reward scaling or the random seeds. Among other results, they showed that different implementations of the same algorithm with the same set of hyper-parameters led to drastically different results.1 Available on github at https://github.com/flowersteam/rl-difference-testing

See Wikipedia's article for more details on one-tail versus two-tail tests: https://en. wikipedia.org/wiki/One-_and_two-tailed_tests

Henderson et al. also advised for the Kolmogorov-Smirnov test which tests whether two samples comes from the same distribution. This test should not be used to compare RL algorithms because it is unable to prove any order relation.

Available at http://www.sjsu.edu/faculty/gerstman/StatPrimer/t-table.pdf.

More explanations and justifications can be found in https://ocw.mit.edu/courses/ mathematics/18-05-introduction-to-probability-and-statistics-spring-2014/ readings/MIT18_05S14_Reading24.pdf.

An implementation of the bootstrap confidence interval test can be found at https: //github.com/facebookincubator/bootstrapped.

Acknowledgement

This research is financially supported by the French Ministère des Armées -Direction Générale de l'Armement.

Conclusion

In this paper, we outlined the statistical problems that arise when comparing the performance of two RL algorithms. We defined type-I, type-II errors and proposed appropriate statistical tests to test for performance difference. Finally and most importantly, we detailed how to pick the right number of random seeds (the sample size) so as to reach the requirements in both error types.

The most important part is what came after. We challenged the hypotheses made by the Welch's t-test and the bootstrap test and found several problems. First, we showed significant difference between empirical estimations of the false positive rate in our experiment and the theoretical values supposedly enforced by both tests. As a result, the bootstrap test should not be used with less than N = 20 samples and tighter significance level should be used to enforce a reasonable false positive rate (< 0.05). Second, we show that the estimation of the sample size N required to meet requirements in type-II error were strongly dependent on the accuracy of (s 1 , s 2 ). To compensate the under-estimation of N , N should be chosen systematically larger than what the power analysis prescribes.

Final recommendations

• Use the Welch's t-test over the bootstrap confidence interval test.

• Set the significance level of a test to lower values (α < 0.05) so as to make sure the probability of type-I error (empirical α) keeps below 0.05.

• Correct for multiple comparisons in order to avoid the linear growth of false positive with the number of experiments.

• Use at least n = 20 samples in the pilot study to compute robust estimates of the standard deviations of both algorithms.

• Use larger sample size N than the one prescribed by the power analysis. This helps compensating for potential inaccuracies in the estimations of the standard deviations of the algorithms and reduces the probability of type-II errors.