Sergei Kucherenko
email: s.kucherenko@imperial.ac.uk

Pietro Belotti
email: belotti@elet.polimi.it

Leo Liberti
email: liberti@lix.polytechnique.fr

Nelson Maculan
email: maculan@cos.ufrj.br

New formulations for the Kissing Number Problem

Keywords: sphere packing, NLP, global optimization, stochastic algorithm, multilevel single linkage, variable neighbourhood search

Determining the maximum number of D-dimensional spheres of radius r that can be adjacent to a central sphere of radius r is known as the Kissing Number Problem (KNP). The problem has been solved for 2, 3 and very recently for 4 dimensions. We present two nonlinear (nonconvex) mathematical programming models for the solution of the KNP. We solve the problem by using two stochastic global optimization methods: a Multi Level Single Linkage algorithm and a Variable Neighbourhood Search. We obtain numerical results for 2, 3 and 4 dimensions.

Introduction

When rigid balls touch each other, in billiard-room terms, they "kiss". In mathematical terms, the kissing number in D dimensions is the number of D- In R 2 the result is trivial: the maximum kissing number is 6 (Fig. 1,a). The situation is far from trivial in R 3 . The problem earned its fame because, according to Newton, the maximum kissing number in 3D was 12, whereas according to his contemporary fellow mathematician David Gregory, the maximum kissing number in 3D was 13 (this fact was stated without proof). This question was settled, at long last, more than 250 years after having been stated, when J. Leech finally proved that the solution in 3D is 12 [START_REF] Leech | The problem of the thirteen spheres[END_REF]. The question for the 4-dimensional case was very recently settled in a yet unpublished paper by O. Musin of Moscow State University [START_REF] Musin | The kissing number in four dimensions[END_REF] which shows that the solution of KNP(4) is 24 spheres. In this paper, we propose a mathematical programming approach based on two nonconvex continuous models, which we solve with two global optimization algorithms for the 2, 3 and 4 dimensional cases. One of these algorithms is a quasi Monte Carlo variant of the Multi Level Single Linkage (MLSL) algorithm called SobolOpt [START_REF] Kucherenko | Application of deterministic low-discrepancy sequences to nonlinear global optimization problems[END_REF], the other is a novel implementation of Variable Neighbourhood Search for constrained NLPs. Previously, the KNP had been used as a validation method for a local optimization code [START_REF] Martínez | Augmented lagrangians and sphere packing problems[END_REF][START_REF] Martínez | Validation of an augmented lagrangian algorithm with a gauss-newton hessian approximation using a set of hard-spheres problems[END_REF].

The models

We propose two models both of which rely on a variable being maximized (or minimized) with an associated threshold, called the feasibility indicator. If the globally optimal value of the feasibility indicator is higher than the threshold, an N -sphere packing in D dimensions is feasible, otherwise it is infeasible. Using a bisection principle, we can quickly pinpoint the optimal solution N * to KNP(D) as the minimum N such that the N + 1 packing is infeasible.

We propose two formulations for the KNP. The first one, which we call the "distance formulation", is a special case of a more general formulation found in the technical report [START_REF] Maculan | Bounds on the kissing numbers in R n : Mathematical programming formulations[END_REF]. Given parameters D (dimension of Euclidean space) and N (number of spheres in the packing), the variables x i = (x i 1 , . . . , x i D ), 1 ≤ i ≤ N determine the position of the center of the i-th sphere around the central one. We maximize the feasibility indicator α ≥ 0 (with threshold 1), which represents the minimum pairwise sphere separation distance in the N -sphere configuration being tested, subject to the necessary geometric constraints. Since the constraints are nonconvex, there may be multiple local minima. If the solution is α ≥ 1, then there is enough space for N spheres, otherwise there are overlapping spheres and

N > N * . max α (1) ∀i ≤ N ||x i || 2 = 4 (2) ∀i < j ≤ N ||x i -x j || 2 ≥ 4α (3) α ≥ 0 (4) ∀i ≤ N x i ∈ R D , (5)
where the norm || • || is taken to be the Euclidean norm. Constraints (2) ensure that the centers of the N spheres all have distance 2 from the center of the central sphere (i.e., the N spheres "kiss" the central sphere). Constraints (3) makes the N spheres non-overlapping. In practice we perform a simplification based on the following:

||x i -x j || 2 = D k=1 (x i k -x j k ) 2 = D k=1 ((x i k ) 2 + (x j k ) 2 -2x i k x j k ) = = ||x i || 2 + ||x j || 2 -2 D k=1 x i k x j k .
By Eq. ( 2) we have

||x i -x j || 2 = 8 -2 D k=1 x i k x j k , thus constraint (3) can be reformulated to ∀i < j ≤ N 2α + D k=1 x i k x j k ≤ 4. ( 6
)
This reduces the number of nonlinear terms in the problem and makes it somewhat faster to solve in a local stage of a global optimization algorithm.

We propose a second formulation, which we call the "angle formulation". As all centres have distance 2 from the centre of the central sphere, their position is uniquely identified by their spherical coordinates ϑ i k , for k ∈ {1, 2, . . . , D -1} and i ∈ {1, 2, . . . , N }. Adopting spherical coordinates ϑ i k , the cartesian coordinates x k i of the sphere centers are defined as

x i k = cos ϑ i k d-k h=1 sin ϑ i h .
We obtain a formulation with polynomial constraints by considering the sines and cosines of each ϑ i k defined as variables σ i k = sin ϑ i k and γ i k = cos ϑ i k , respectively. The angle formulation can be defined as follows (we assume

γ i D = 1 for each i ≤ N ): min α (7) ∀i ≤ N, k < D (σ i k ) 2 + (γ i k ) 2 = 1 (8) ∀i < j ≤ N D k=1 γ i k γ j k D-k h=1 σ i h σ j h ≤ 0.5 + α (9) ∀i ≤ N σ i ∈ [-1, 1] D-1 , γ i ∈ [-1, 1] D . (10)
Constraints ( 8) impose that sin 2 ϑ i k + cos 2 ϑ i k = 1, while constraints [START_REF] Maculan | Bounds on the kissing numbers in R n : Mathematical programming formulations[END_REF] require that the angle between any two vectors must not exceed π 3 (which implies that the spheres will not overlap), where the feasibility indicator α is employed in the minimization direction. If for a given (N, D) the globally optimal α is zero, then we have found a sphere packing.

The methods

A stochastic approach for global optimization, in its simplest form, consists only of random search and it is called Pure Random Search (PRS). In PRS an objective function f (x) is evaluated at P randomly chosen points and the smallest value of f (x) is taken as the global minimum. Advanced stochastic techniques use stochastic methods to search for local minima and then utilize deterministic methods to solve a local minimisation problem. Two phases are considered: global and local. In the global phase, the function is evaluated in a number of randomly sampled points from a uniform distribution over a unit Hypercube H n . In the local phase the sample points are used as starting points for a local minimization search. The efficiency of the multistage methods depends both on the performance of the global stochastic and the local minimization phases.

In the most basic form of the multistage approach a local search is applied to every sample point. Inevitably, some local minima would be found many times. Since the local search is the most CPU-time consuming stage, ideally it should start just once in every region of attraction. This is the idea behind various versions of the so-called clustering methods. Extensive reviews on this subject can be found in [START_REF] Törn | Global Optimization[END_REF][START_REF] Schoen | Two-phase methods for global optimization[END_REF]. One of the most efficient clustering methods is a MLSL algorithm developed by Rinnooy-Kan and Timmer in [START_REF] Rinnooy-Kan | Stochastic global optimization methods; part I: Clustering methods[END_REF][START_REF] Rinnooy-Kan | Stochastic global optimization methods; part II: Multilevel methods[END_REF].

The efficiency of MLSL depends on the quality of sampled points. It has been recognized through theory and practice that uniformly distributed deterministic sequences provide more accurate results than purely random sequences. Low-discrepancy sequences (LDS) are designed specifically to place sample points as uniformly as possible. Unlike random numbers, successive low discrepancy points "know" about the position of their predecessors and fill the gaps left previously. Methods based on LDSs are known as quasi Monte Carlo (QMC) methods. In the majority of applications, QMC methods have superior performance compared to that of MC methods. Improvement in time-toaccuracy using QMC methods can be as large as several orders of magnitude. It was shown in [START_REF] Kucherenko | Application of deterministic low-discrepancy sequences to nonlinear global optimization problems[END_REF] that application of LDS can significantly increase the efficiency of MLSL methods. Central to the QMC approach is the choice of LDS. Different principles were used for constructing LDSs by Holton, Faure, Sobol', Niederreiter and others. Many practical studies have proven that Sobol' LDS in many aspects are superior to other LDSs [START_REF] Paskov | Faster evaluation of financial derivatives[END_REF][START_REF] Sobol | On quasi-monte carlo integrations[END_REF].

Variable Neighbourhood Search (VNS) is a relatively recent metaheuristic method which relies on iteratively exploring neighbourhoods of growing size to identify better local optima [START_REF] Hansen | Variable neighbourhood search[END_REF][START_REF] Hansen | Variable neighbourhood search: Principles and applications[END_REF]. More precisely, VNS escapes from the current local minimum x * by initiating other local searches from starting points sampled from a neighbourhood of x * which increases its size iteratively until a local minimum better than the current one is found. These steps are repeated until a given termination condition is met. The search space is defined as the hyper-rectangle given by the set of variable ranges x L ≤ x ≤ x U . For each k ≤ k max we define hyper-rectangular neighbourhoods N k (x * ) with side lengths proportional to those of x L ≤ x ≤ x U , centered at x * , whose sides have been scaled by k kmax . More precisely, we let N k (x * ) be the hyper-rectangle y L ≤ x ≤ y U where, for all i ≤ n:

y L i = x * i - k k max (x * i -x L i ), y U i = x * i + k k max (x U i -x * i ).

Computational results

We solved the KNP by using solvers implementing the methods above within the framework of the general-purpose global optimization software framework ooOPS [START_REF] Liberti | ooOPS[END_REF]. Both solvers call a local constrained NLP solver code (SNOPT [START_REF] Gill | User's Guide for SNOPT 5.3[END_REF]) to perform the local descent.

Our computational results were all obtained with a 800MHz Intel Pentium III CPU with 384 MB RAM running Linux. In Table 1 Table 2 Computational results for the angle model ( 7)- [START_REF] Martínez | Augmented lagrangians and sphere packing problems[END_REF]. CPU is in seconds of user time.

Our results show that the distance formulation is more promising, in terms of user CPU time, than the angle formulation, for both SobolOpt and VNS global solvers.

Conclusion

In this paper two nonconvex continuous mathematical programming formulations to solve the decision problem "is N the Kissing Number in Euclidean D-space?" were presented. The Kissing Number was correctly determined in 2, 3 and 4 dimensions, using two different stochastic global optimization algorithms. Research is ongoing for obtaining computational results for the smallest open case, namely 5 dimensions.

Fig. 1 .

 1 Fig. 1. The problem in R 2 (a) and R 3 (b).Figure (b) was taken from [18].

  Fig. 1. The problem in R 2 (a) and R 3 (b). Figure (b) was taken from [18]. spheres of unit radius that can be arranged around a central D-sphere of unit radius so that each of the surrounding spheres touches the central one without overlapping. Determining the maximum kissing number in various dimensions has become a well-known problem in Combinatorial Geometry. Notationally, we indicate the Kissing Number Problem in D dimensions by KNP(D).

  we report results from the first model; n, m are the number of variables and constraints in the problem, rows SobolOpt and VNS report the globally optimal feasibility indicator values and user CPU times (in seconds) taken by the SobolOpt and VNS solvers to globally solve the problem. Each column corresponds to a different (N, D) pair. In 2D, 6 spheres are a feasible packing, and 7 are not; in 3D, 12 spheres are feasible and 13 are not; in 4D, 24 and 25 are known lower and upper bounds to the KNP[START_REF]Sphere Packings, Lattices and Groups[END_REF]. Table2contains the corresponding results for the second model. For each instance and solver, the correct feasibility indicator values were identified.

	Threshold: 1	D = 2	D = 3	D = 4
	Max α	N = 6	N = 7	N = 12	N = 13	N = 24	N = 25
	n	13	15	37	40	97	101
	m	21	28	78	91	300	325
	SobolOpt: α	1.00000	0.75302	1.10558	0.91473	1.00000	0.92537
	(CPU)	5.05	6.24	31.35	40.72	334.36	369.55
	VNS: α	1.00000	0.75302	1.10558	0.9147	1.00000	0.92537
	(CPU)	0.31	0.44	1.88	3.39	30.07	24.35
	Table 1
	Computational results for the distance model (1)-(5). CPU is in seconds of user
	time.
	Threshold: 0	D = 2	D = 3	D = 4
	Min α	N = 6	N = 7	N = 12	N = 13	N = 24	N = 25
	n	13	15	37	40	97	101
	m	21	28	90	104	348	375
	SobolOpt: α	0.00000	0.12349	0.00000	0.04263	0.00000	0.03731
	(CPU)	2.75	5.32	120.58	215.25	4861.53	5565.55
	VNS: α	0.00000	0.12349	0.00000	0.04263	0.00000	0.03731
	(CPU)	0.62	0.69	9.76	20.56	850.22	816.23

Preprint submitted to Elsevier Science


S Cossin

S Malavaud

P Jarno

M Giard

F L'hériteau

L Simon

L Bieler

L Molinier

B Marcheix

A.-G Venier

O Ali-Brandmeyer

C Neels

M Aupée

M Perennec

P Astagneau

F Daniel

C Campion

C Bernet

E Caillat-Vallet

C Bervas

E Reyreaud

P Baillet

J L Yannick

V Jost

L Merle

J Merlo

S Christophe

F Malavaud

J.-M Bruyere

S Thiolet

Barquin- Guichard

Surgical site infection after valvular or coronary artery bypass surgery: 2008 2011 French national SSI-RAISIN surveillance

Keywords:

Background: Multisite information regarding surgical site infection (SSI) rates for cardiac surgery programmes is not widely available. Ward characteristics that may affect outcomes have not been analysed previously. Aim: To determine individual-and ward-level factors associated with SSI occurrence after coronary artery bypass grafting (CABG) and valvular surgery. Methods: A dataset from the French national database SSI-RAISIN 2008 2011 was used.

Only adult patients were included. A standardized questionnaire was completed for each patient who underwent surgery, and patients with and without SSI were characterized.

Patients and ward risk factors for SSI were analysed using a multilevel logistic regression model with SSI as binary outcome (two levels: patient and ward).

Results: Out of 8569 patients from 39 wards, the SSI rate was 2.2%. Micro-organisms were isolated in 144 patients (74%): 35% coagulase-negative staphylococci (N = 51), 23%

Introduction

Surveillance systems (along with infection control programs) can reduce surgical site infections (SSI) rates. 1 3 SSI after cardiac surgery increase the patient morbidity and mortality, the length of hospital stay and hospital costs. [START_REF] Alasmari | Temporal trends in the incidence of surgical site infections in patients undergoing coronary artery bypass graft surgery: a populationbased cohort study, 1993 to 2008[END_REF] Multisite information regarding SSI rates for cardiac surgery programmes is not readily available. Moreover, SSI risk factors after coronary artery bypass graft (CABG) and valvular surgery have rarely been assessed by population-based or multicentre studies. Previously identified risk factors, such as obesity, are generally not amenable to corrective action. 4 7 By identifying the role of the specific risk factors related to the treating facility, such as ward characteristics, it is possible that interventions could be undertaken to reduce SSI rates. [START_REF] Beckmann | Sternal surgical site infection prevention is there any room for improvement?[END_REF] The aim of this study was to describe SSI after CABG and valvular surgery from the French national database SSI-RAISIN 2008 2011, in order to identify significant individual-and ward-level factors associated with SSI occurrence.

Methods

Population

A dataset of adult patients who underwent a valvular or CABG surgical procedure during the four-year period 2008 2011 was analyaed: 8569 patients from 39 wards from 31 healthcare facilities.

Data collection

Yearly surveillance included at least three months during the annual survey period (January to June). Details of data collection and management methods for the SSI-RAISIN The following preoperative characteristics were considered: age, sex, duration of preoperative hospitalization (<48 h vs 48 h), wound class (Altemeier classification), the American Society of Anesthesiologists (ASA) physical status score, duration of surgery ( 75th percentiles versus >75th percentiles; the 75th percentile was 5 h for CABG with the use of internal thoracic arteries, 4 h with the use of saphenous vein and for valvular procedures), emergency versus elective surgery, ambulatory surgery, type of procedure, duration of postoperative follow-up, SSI occurrence during hospitalization or after discharge, depth of infection (superficial, deep, organ/space), average delay between the surgery and SSI diagnosis, revisions surgery, patient outcome (alive/deceased), National Nosocomial Infections Surveillance system (NNIS) risk index of the surgery (based on risk factors for infection: ASA classification >2, duration of surgery >75th percentile, and wound contamination >2). [START_REF] Altemeier | Definitions and classifications of surgical infections. Manual on control of infections in surgical patients[END_REF] The following ward-based characteristics were considered: hospitalization ward, global cardiac surgery SSI rate, valvular and CABG SSI rate, the number of cardiac surgery procedures per month, the average duration of patients' follow-up after surgery, and the proportion of NNIS 0 patients.

Statistical analysis

First, a descriptive analysis of patients with and without SSI was performed with the R ® software version 2.15.3.

Risk factors for SSI were estimated with a multilevel logistic regression model with SSI as a binary outcome. [START_REF] Goldstein | Multilevel statistical models[END_REF] Data structure had two levels: patients (first level) and wards (second level). An 'empty' model (model A) only included a random intercept and allowed detection of a possible contextual dimension. Thereafter, the individual characteristics were

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT
4 included in the model (model B) to investigate the extent to which ward incidence differences were explained by the individual patients of the wards. Ward characteristics were included in a third model (model C). Multivariate analyses were performed by stepwise logistic regression. All variables with P 0.10 in the univariate analysis or which were thought to be associated with SSI from the literature were included in the models. P < 0.05 was considered significant. 12 15 For quantifying heterogeneity between clusters, the median odds ratio (MOR) was calculated and intraclass correlation (ICC) was estimated with the linear threshold model method. [START_REF] Merlo | A brief conceptual tutorial of multilevel analysis in social epidemiology: using measures of clustering in multilevel logistic regression to investigate contextual phenomena[END_REF] MLwiN ® software 2.15 was used for univariate and multivariate analyses.

Cost assessment

The cost assessment adopted the healthcare payer's perspective, the healthcare insurance, and took into account only hospitalization costs. Hospitalization costs (such as nursing, medical care, pharmacy, diagnostic tests, staff, general equipment, administration, security, central supply, dietetics or social services) were assessed using the French diagnosisrelated groups that fund all costs. [START_REF] Ye | Incidence and costs of bleeding-related complications in French hospitals following surgery for various diagnoses[END_REF][START_REF] Yilmaz | The influence of social deprivation on length of hospitalisation[END_REF] The average additional cost of valvular or coronary artery bypass surgery with a deep SSI compared to the same support without SSI was estimated from the database of the University Hospital of Toulouse (€11,000). The total additional cost was calculated by multiplying the total number of deep SSIs by the estimated average incremental cost. Costs were expressed in euros.

Results

The SSI rate was 2.5% in 2008 (52/2071), 1.9% in 2009 (35/1815), 2.2% in 2010 (45/2075), and 2.3% in 2011 (62/2698) with no significant difference between years (P = 0.60).

Ward characteristics

Only five wards participated in the surveillance every year. The number of valvular and CABG surgeries per month ranged from 11 to 111 [median: 44; interquartile range (IQR):

26 77]. The proportion of NNIS-0 patients on the various wards ranged from 0 to 80% (median: 2.5%; IQR: 1 6%).

Six healthcare facilities had two different surgery wards performing cardiac surgery and one had three different wards: as there were substantial differences in SSIs between wards in the same healthcare facility, these wards were considered separately.

Ward SSI rates ranged from 0% to 5% (median: 1.7%; IQR: 0.7 2.9; cf. Figure 1). At ward level, the mean duration of patient's follow-up was 32 days (mean: 30; IQR: 9 41).

Patient characteristics

Out of 8659 CABG or valvular procedures, the overall SSI rate was 2.2%. Patient characteristics are recorded in Table I: 61% of SSIs were diagnosed during hospitalization,

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 5 45%
were superficial SSIs and 48% deep or organ space (sternitis, mediastinitis, and endocarditis) SSI. Depth of infection was not specified for 12 patients (6%). Average (±SD) time from surgery to SSI onset was 14 days (±7.14) without significant difference being noted between CABG and valvular surgery. The mean duration of hospital stay after surgery was 23 days in SSI patients versus 10 days in those without SSI (P < 0.01).

Mortality was 4.6% in SSI patients (9/194) versus 2.4% (205/8465) in patients without infection (P = 0.06); further surgery was required in 45.2% of the SSI patients (only two patients without SSI underwent an additional surgery).

Comparing the patients with deep-organ space sternal SSI versus the patients with superficial SSI, there were no differences with respect to the interval to diagnosis, causative micro-organisms or NNIS score. However, patients with deep SSI tended to have a longer length of hospital stay (32 days versus 19 days, P < 0.001), a longer duration of follow-up (45 days versus 30 days, P < 0.001), a higher mortality (9.6% versus 0%), and a higher risk of revision surgery (67.8% versus 20%, P < 0.001).

Microbiology

Micro-organisms were isolated in 144 (74%) patients: 35% coagulase negativestaphylococci (N = 51), 23% Staphylococcus aureus (N = 33) including five meticillinresistant S. aureus (MRSA), 6% Escherichia coli (N = 8), 5% Pseudomonas aeruginosa (N = 7), 3% Enterococcus faecalis (N = 4), 3% other Gram-positive bacteria (N = 4), 2%

Citrobacter koseri (N = 3), 2% Serratia spp. (N = 2), and 8% other Gram-negative bacteria (N = 11). Five patients had two micro-organisms identified, 15 had three (10%). Cultures were negative in eight patients, not performed in five patients and information was missing in 37.

S. aureus was found in six out of nine patients (67%) who developed SSI and died before discharge, including two patients with MRSA. No trends between years were observed.

Results of multivariate analysis are presented in Table II. Model A was an empty model with no explanatory variable. The part of the residual variance attributed to ward level was 7% and was significantly different from zero (i.e. 7% of total variability of the outcome variable was explained by differences between wards). In a model with patients' characteristics only (model B), a higher probability of SSI was associated with the duration of preoperative hospitalization, the duration of follow-up, the duration of surgery >75th In univariate analysis, only the ward SSI rate was significantly associated with higher probability of SSI. In a multivariate model with both patient and ward characteristics (model C), a higher probability of SSI was associated with the duration of surgery >75th percentile and the ward SSI rate.

Cost

The total additional cost of SSI for facilities was estimated to be €291,000 for 94 patients during the four years studied. If hospital financing did not reimburse a part of the costs of the SSI, the costs would have increased to €1,034,000 for the facilities.

Discussion

This analysis using a national dataset found a 2.2% incidence rate of SSI after valvular and CABG surgery within the first 30 days following surgery. This result is consistent with previous nationwide studies. [START_REF] Berg | National surveillance of surgical site infections after coronary artery bypass grafting in Norway: incidence and risk factors[END_REF][START_REF] Manniën | Surveillance and epidemiology of surgical site infections after cardiothoracic surgery in The Netherlands, 2002 2007[END_REF] A significant heterogeneity of SSI occurrence between wards was found. This is, to our knowledge, the first nationwide study showing the contribution of ward risk factors in cardiac surgery SSI.

Several limitations to the present study should be noted. Whereas the protocol recommended a minimum follow-up of 30 days, methods differed between the wards. [START_REF] Perennec-Olivier | Surveillance des infections du site operatoire, France 2011 [Surveillance of surgical site infections[END_REF] Despite this issue, we believe that underdiagnosis of infection was unlikely since no patient developed an SSI after 22 days. [START_REF] Wilson | Hospitals in Europe Link for Infection Control through Surveillance (HELICS). Inter-country comparison of rates of surgical site infection opportunities and limitations[END_REF][START_REF] Manniën | Impact of postdischarge surveillance on surgical site infection rates for several surgical procedures: results from the nosocomial surveillance network in The Netherlands[END_REF] The study was restricted to the wards that carried out more than 10 cardiac procedures in order to prevent the great variability of small samples which could generate bias in the multilevel analysis. An additional issue which could have affected analysis was the fact that only one type of procedure could be defined on the surveillance form. Surgeons had to choose between CABG or valvular procedure, even when patients underwent a combined procedure. Nevertheless, our findings were in line with literature showing a higher rate of SSI after CABG than after valvular repair. [START_REF] Manniën | Surveillance and epidemiology of surgical site infections after cardiothoracic surgery in The Netherlands, 2002 2007[END_REF][START_REF] Eklund | Mediastinitis after more than 10,000 cardiac surgical procedures[END_REF][START_REF] Mediastinitis | Risk factors for deep sternal wound infection after sternotomy: a prospective, multicenter study[END_REF] This study also confirmed the impact of several previously identified individual risk factors. The impact of the preoperative length of stay was confirmed as noted in previous reports. This phenomenon is probably multifactorial, including increased colonization with hospital pathogens or progressive impairment of the patient's nutritional status and is unlikely to be impacted by interventions. 23 Since this factor disappeared in the model along with ward characteristics, further investigation of the phenomenon should be carried out.

Duration of surgery may reflect either the patient's preoperative morbidity, operating room logistics, or difficulties during the surgery. [START_REF] Mediastinitis | Risk factors for deep sternal wound infection after sternotomy: a prospective, multicenter study[END_REF] Cristofolini et al. pointed out that the ASA risk group may not be applicable in cardiac surgery since almost all patients exhibit an ASA

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

7 score 3, which was the case in our study. 24 The wound contamination class seems irrelevant in the types of cardiac surgery we studied; the NNIS score index should therefore not be used to compare SSI rates after cardiothoracic operations among surgeons or institutions. 25 The logistic EuroSCORE I and II or the Society of Thoracic Surgeons risk score (STS score) should be more appropriate than NNIS for cardiac surgery but was not available at the time of the study. 26 Coagulase-negative staphylococci and S. aureus were implicated in half of the cases of SSIs. These organisms have various virulence factors and are usually found on the skin and the environment which can explain their predominance in SSI after cardiac surgery. 27,28 The importance of these pathogens has led to recommendations for the use of nasal mupirocin for all cardiac surgery cases. 29 We used a multilevel modelling technique to take into account the hierarchical data structure and the correlation of within-group observations. This study showed a significant correlation of infections and the ward location of the patient after adjusting for patient-level variables. In an empty model, we found a between-patient and between-ward heterogeneity.

MOR indicated that the ward of hospitalization would be relevant for understanding variations of the individual probability of acquiring SSI. Healthcare-associated infection prevention practices may explain variation in rates between healthcare facilities. A higher incidence of SSI in the ward was associated with higher probability of SSI for a patient: this relationship could be a reflection of the ward ecology, colonization pressure or suboptimal compliance with infection control measures, none of which was could be evaluated in this survey. Furthermore, organization and patient safety culture differences may also have a role. 30 Although the use of intraoperative antibiotic bone paste and the use of antibioticimpregnated sutures were associated with a lower risk of developing SSIs, the differences were not significant. 22 SSI surveillance is considered to be an effective strategy for reducing SSI incidence, as well as a tool for surgical staff to promote awareness and prevention, but French SSI surveillance was unable to identify specific reasons for ward-related infection risks. Indeed, the findings of our study suggest that more precise and detailed ward characteristics should be included in future cardiac surgery SSI studies in the hope of identifying interventions which could further decrease the risk of these infections.

The present study provided clear evidence of the impact of ward-level variables on the occurrence of SSI in CABG and valvular surgery. Although patient risk factors were more strongly associated with the SSI occurrence, the impact of ward characteristics was not negligible. Ward risk factors should be evaluated by the surgery team and infection control teams by performing root cause analysis to improve patients' quality of care. 31,32

Author queries

1. Please confirm full postal address for corresponding author.

2. Please check the names 'J. Christophe' and 'Yannick' in the study group list.

  previously published.[START_REF] Astagneau | Reducing surgical site infection incidence through a network: results from the French ISO-RAISIN surveillance system[END_REF][START_REF] Rioux | Impact of a six-year control programme on surgical site infections in France: results of the INCISO surveillance[END_REF] SSI diagnostic criteria used by surgical teams accorded with the standardized Centers for Disease Control definitions. This study was given ethical clearance by the French data protection authority.A sternal wound infection included at least one of the following: (i) presence of purulent drainage; (ii) bacterial isolation from the incision; (iii) presence of at least one of the following symptoms: tenderness, swelling, redness, or heat; (iv) diagnosis made by the surgeon or physician.A standardized data collection form was completed for each patient who underwent surgery during the surveillance periods. For cardiac surgery, the classification of the procedure was categorized as follows: closed-heart surgery, open-heart surgery, CABG with saphenous vein graft, CABG with internal thoracic arterial graft, valvular surgery, and intrathoracic arteries surgery as categorized by the French Medical Acts Classification.

  percentile. ASA score, wound contamination class, sex, and age were not significantly associated with SSI. Model B included the separated components of NNIS score and was preferred to a model including NNIS [NNIS score: OR: 1.67 (95% CI 1.23 2.270); MOR: ; ICC: 0.06]. The residual heterogeneity between wards (MOR: 1.53) was equal to the duration of preoperative hospitalization (OR: 1.57).

13 Figure 1 .

 131 Figure 1. Distribution of the 2008 2011 surgical site infection (SSI) rate in valvular and coronary artery bypass graft surgery, French national SSI Surveillance ISO RAISIN (N =39 wards).

3 .

 3 Please note that original ref. 10 was a duplicate of ref. 1 and has been removed please check renumbering.

Table II

 II Multivariate analysis of risk factors for surgical site infection after coronary artery bypass graft or valvular surgery (SSI-RAISIN 2008 2011)

	Variable	Model A	Model B	Model C
	Intercept	3.881 (0.118)	4.278 (0.433)	4.524 (0.463)
	Patient characteristics
	Age 65 years		0.755	0.988
			(0.550 1.035)	(0.974 1.001)
	Male sex		1.023	1.018
			(0.729 1.436)	(0.725 1.429)
	Duration of preoperative		1.573	1.462
	hospitalization 48 h		(1.138 2.173)	(1.062 2.012)
	Duration of postoperative		1.007	0.898
	follow-up		(1.001 1.013)	(0.660 1.219)
	ASA score >2		1.274	1.309
			(0.583 2.778)	(0.593 2.869)
	Wound contamination class >2		0.880	0.738
			(0.209 3.775)	(0.174 3.131)
	Duration of surgery >75th		1.820	1.756
	percentile		(1.299 2.550)	(1.275 2.416)
	Ward characteristics
	Percentage of NNIS-0 patients			1.006
				(0.965 1.048)
	No. of valvular and CABG			1.000
	surgeries per month			(0.992 1.007)
	Ward SSI rates			1.724
				(1.480 2.000)
	Median OR (95% CI)	1.59	1.53	1.05
		(1.431 1.791)	(1.451 1.760)	(1.010 1.08)
	Intraclass correlation coefficient	0.067	0.058	0.001
	NNIS, National Nosocomial Infections Surveillance system; ASA, American Society of

Anesthesiologists; CABG, coronary artery bypass graft; OR, odds ratio; CI, confidence interval; NNIS-0 patient: ASA score 2 + wound contamination class 2 + duration of surgery 75th percentile.

Acknowledgements

We thank the infection control and surgical teams at the healthcare facilities who participated in the RAISIN-ISO surveillance.

Funding sources

The French national surgical site infection surveillance system is supported by the French Ministry of Health.

Conflict of interest statement

None declared.


Jessica Chemali

Alessandro Lazaric

Direct Policy Iteration with Demonstrations

We consider the problem of learning the optimal policy of an unknown Markov decision process (MDP) when expert demonstrations are available along with interaction samples. We build on classification-based policy iteration to perform a seamless integration of interaction and expert data, thus obtaining an algorithm which can benefit from both sources of information at the same time. Furthermore, we provide a full theoretical analysis of the performance across iterations providing insights on how the algorithm works. Finally, we report an empirical evaluation of the algorithm and a comparison with the state-of-the-art algorithms.

Introduction

We study the problem of learning the optimal policy of an unknown Markov decision process (MDP) by reinforcement leaning (RL) when expert samples (limited in number and/or quality) are available. In particular, we identify in approximate policy iteration (API) [START_REF] Bertsekas | [END_REF] the most suitable learning scheme to integrate expert demonstrations. API is an iterative learning scheme where one has access to batch data (i.e., RL samples) obtained by a sampling policy. Starting at an arbitrary policy, API goes through two stages at every iteration: 1) policy evaluation first computes the value function of the current policy, 2) policy improvement computes the greedy policy w.r.t. the value function. At any of the two steps, API may introduce errors that depend on the use of a finite number of samples and the choice of a specific approximation space (e.g., linear regression). Because of such errors, API is not guaranteed to improve policies monotonically at each iteration and, in general, it only converges to a set of policies with bounded performance loss w.r.t. the optimal policy [Munos, 2003]. In this work, we propose to improve such learning scheme by modifying the policy improvement step to bias it towards policies that agree, to some degree, with demonstrated expert actions.

Related work. While most of past literature focused on learning from either RL or expert samples [START_REF] Ng | [END_REF][START_REF] Abbeel | Apprenticeship learning via inverse reinforcement learning[END_REF]Bagnell and Ross, 2010], there is a recent surge of interest in combining expert demonstrations in RL algorithms [START_REF] Kim | [END_REF][START_REF] Piot | [END_REF] or cost functions in imitation learning [START_REF] Ross | [END_REF]. [START_REF] Kim | [END_REF] proposed API with demonstrations (APID) which modifies regularized-LSPI [START_REF] Farahmand | [END_REF] by adding linear soft max-margin constraints in the policy evaluation step; the constraints require value functions for which expert actions are greedy. [START_REF] Kim | [END_REF] upper bound the Bellman error of APID at any iteration and empirically show that it outperforms both LSPI and imitation learning algorithms. Although the empirical results are a convincing proof of concept, the theoretical analysis is limited to single-iteration errors and does not provide a clear understanding on how effectively RL and expert data are integrated. In fact, while the objective of LSPI is to compute the value function of the current policy, the constraints favor functions with greedy actions equivalent to the expert's. These two requirements may be contrasting (e.g., the greedy policy of the exact value function may not be similar to the expert) and may point towards solutions which are neither good to estimate the target value nor in satisfying the desired constraints. [START_REF] Piot | [END_REF] build on a similar idea but integrate expert constraints directly into the minimization of the optimal Bellman residual. This integration is very appealing since it "regularizes" the global objective of approximating the optimal value function rather than constraining each iteration of API, in which the constraints and the objective over iterations may not always be coherent. Unfortunately, the resulting optimization problem is non-convex, non-differentiable and biased, and so they have to resort to a non-parametric method where the conditional transition matrix is embedded in RKHS and the minimization is carried out by a boosting algorithm. Although the resulting algorithm, called RLED, outperforms APID in practice, the complexity of the optimization prevents from deriving theoretical guarantees on its performance.

Finally, the integration of expert demonstrations has been investigated in Bayesian model-based RL in [Doshi-Velez et al., 2010], where a non-parametric Bayesian approach is employed to combine model knowledge from the expert trajectories and RL samples obtained from independent exploration.

In this paper, we build on classification-based policy iteration [START_REF] Lagoudakis | [END_REF], notably the direct policy iteration (DPI) [START_REF] Lazaric | [END_REF] algorithm, where policy improvement is replaced by a classification problem. Since imitating the expert policy is intrinsically a classifica-Algorithm 1 Direct Policy Iteration with Demonstrations 1: Input: policy space Π, state distribution ρ, expert weight α 2: Initialize: arbitrary policy π0 ∈ Π 3:

for k = 1 to K do 4: Construct the rollout set D k roll = {xi} N RL i=1 , xi iid ∼ ρ 5:
for xi ∈ D k roll and actions a ∈ A do 6:

for j = 1 to M do 7:

Perform a rollout according to policy π k 8: R

π k j (xi, a) = r(xi, a) + H t=1 γ t r(x t , π k (x t )), 9:
end for 10:

Q π k (xi, a) = 1 M M j=1 R π k j (xi, a) 11: Build RL training set D k RL = D k RL ∪ {xi, a * i , wi} 12: end for 13: Build expert training set D k E = {zj, aj, α} N E j=1 14: Build training set D k = D k E ∪ D k RL 15: π k+1 = arg min π∈Π Lπ k (ρ, μ, π) (Classification)
16: end for tion problem, the integration of the expert and RL data sums up to merging the two datasets. Thanks to its simple form, the resulting algorithm (DPID) (Sect. 3) does not increase the complexity of DPI and its performance across iterations can be theoretically analyzed (Sect. 4). Experimental results in two simulated domains show the potential of DPID compared to existing methods (Sect. 5).

Preliminaries

We consider a discounted MDP defined by the tuple (X , A, p, r, γ), where the state space X is bounded in R d , the action space A is finite1 , the kernel p : X × A → X defines the transition probabilities, the reward r : X × A → [-R max ; R max ] is bounded, and γ ∈ [0, 1) is a discount factor. A deterministic policy π is a mapping X → A and Π is a set of policies. We denote by V π : X → R and Q π : X ×A → R the value and action value functions of policy π. V π is the unique fixed-point of the Bellman operator T π defined as

(T π V )(x) = r(x, π(x)) + γ X p(dy|x, π(x))V (y), while Q π (x, a) = r(x, a) + γ X p(dy|x, a)V π (y).
The optimal policy π maximizes V π over all states and we denote by V and Q its corresponding functions. V is also the unique fixed-point of the optimal Bellman operator T defined as

(T V )(x) = max a [r(x, a) + γ X p(dy|x, a)V (y)], while Q (x, a) = r(x, a) + γ X p(dy|x, a)V (y).
We say that π is greedy with respect to an action value function Q, if in any state x, π(x) ∈ arg max a∈A Q(x, a). Finally, we introduce the greedy policy operator G mapping policies to greedy policies as (Gπ)(x) = arg max a∈A Q π (x, a).

Direct Policy Iteration with Demonstrations

Building on direct policy iteration (DPI) [START_REF] Lazaric | [END_REF], we introduce DPI with demonstration (DPID), where demonstrations from an expert policy ξ are available. In the original DPI, at each iteration k, a rollout set {x i } N RL i=1 is sampled from a given sampling distribution ρ over X and the value of the current policy π k is estimated on those states using M rollouts truncated to H steps. The resulting estimates Q π k (x i , a) are then used to build a (cost sensitive) training set

D k RL = {x i , a * i , w i } N RL i=1
, where a * i is the best estimated action (i.e., a * i = arg max a Q π k (x i , a)) and

w i = Q π k (x i , a * i ) -Q π k (x i , a +
) is the cost of misclassifying the sample with a + the other action in A. Since A is finite, a * i can be interpreted as the label of x i and thus DPI reduces to solving a (binary) cost-sensitive classification problem using a given policy space Π (e.g., a linear classifier) and it returns a policy π k+1 approximating the greedy policy G( π k ). In DPID we further assume access to a training set

D k E = {z j , ξ(z j )} N E j=1
, where states z j are drawn from a distribution µ over X and ξ is an expert policy (not necessarily optimal). The basic idea of DPID is then to define a suitable integration of D k RL and D k E and solve a classification problem on the joint training set. We introduce the following elements. Definition 1. At each iteration k, we define the state and average empirical RL losses of a policy π as

ˆ RL π k (x; π) = max a∈A Q π k (x, a) -Q π k (x, π(x)), L RL π k ( ρ, π) = 1 N RL N RL i=1 [max a Q π k (x i , a) -Q π k (x i , π(x i ))],
where

Q π k is the rollout estimate of Q π k (see Alg. 1
) and ρ is the empirical counterpart of ρ (i.e., average of Dirac distributions at observed states). We also define the state and average empirical expert losses of a policy π as

E (x; π) = I{ξ(x) = π(x)}, L E ( µ, π) = 1 N E N E j=1 E (z j ; π),
where µ is the empirical counterpart of µ. Finally, the average empirical joint loss is

L π k ( ρ, µ, π) = L RL π k ( ρ, π) + αL E ( µ, π)
, where α ≥ 0 is a parameter of the algorithm.

While the RL loss ˆ RL π k (x; π) depends on the policy at iteration k, the expert loss does not change through iterations since it only measures the "similarity" w.r.t. ξ using a 0/1 loss. In the following we will also use

L E α ( µ, π) = αL E ( µ, π).
As reported in Alg. 1, DPID turns the expert training set to a cost-sensitive set by adding a constant penalty α, which allows to move from DPI (α = 0) to a full imitation learning algorithm (α = ∞). Once the policy π k+1 is computed by minimizing the joint loss L π k , the whole process is repeated over iterations. With DPID we propose a seamless integration of RL and expert data in a common training set, which does not affect the complexity of the original DPI algorithm. This is in contrast with APID and RLED, where the integration of expert data in the form of explicit constraints leads to optimization problems which are more complex and potentially computationally challenging (e.g., RLED needs solving a regularized non-convex optimization problem and APID defines a non-classical constrained convex optimization).

In this section we provide a complete analysis of the performance of DPID. We first define the expected counterparts of the losses introduced in Def. 1. Definition 2. At each iteration k, we define the state and average RL losses of any policy π as

RL π k (x; π) = max a∈A Q π k (x, a) -Q π k (x, π(x)), L RL π k (ρ, π) = X max a∈A Q π k (x, a) -Q π k (x, π(x)) ρ(dx).
Similarly, the expert loss of any policy π is

L E (µ, π) = X E (x; π)µ(dx).
Finally, the joint loss is defined as

L(ρ, µ, π) = L RL π k (ρ, π) + αL E (µ, π) While the difference of RL loss L RL π k (ρ, π) w.r.t. its empiri- cal counterpart L RL π k ( ρ, π
) is in both the sampling distribution and use of estimated Q-values, L E only differs in using samples from µ. We also introduce a series of policies of interest:

π RL k+1 = arg min π∈Π L RL π k (ρ, π), π E = arg min π∈Π L E (µ, π), π k+1 = arg min π∈Π L π k (ρ, µ, π).
These policies are the minimizers within the policy space Π of the losses defined above and thus they are the best policies in approximating the greedy policy G( π k ), the expert policy ξ, and the mixture obtained by averaging the two losses respectively. In the following, we assume that Π has a finite VC-dimension h = V C(Π) and we prove a bound on the loss of the policy returned by DPID at each iteration. Theorem 1. Let π k+1 be the policy returned by DPID and

RL = 20Q max 2 N RL h log eN RL h + log 8/δ , E = 8 N E h log eN E h + log 4 δ .
Then with probability (1 -6δ)

L π k (ρ, µ, πk+1 ) ≤ L π k (ρ, µ, π k+1 ) + 2 all .
with all = RL + α E . Furthermore,

L RL π k (ρ, πk+1 ) ≤ L RL π k (ρ, π RL k+1 ) + 2 all + L E α (µ, π RL k+1 ), L E α (ρ, πk+1 ) ≤ L E α (µ, π E ) + 2 all + L RL π k (ρ, π E ).
Remark (the bound). The first statement refers to the quality of π k+1 in minimizing the joint loss and, as usual, two sources of error appear. The first is the approximation error related to the best policy in Π minimizing the loss L π k , while the second term is the estimation error due to the limited number of samples. To reduce the estimation error, both RL and expert sample sizes should increase. The second and third statements bound the RL and expert losses of π k+1 and compare it to the best possible policy in Π minimizing L RL and L E respectively. In both cases, besides the approximation and estimation errors, π k+1 suffers from an error depending on the errors of π RL k+1 and π E k+1 measured w.r.t. the expert and RL losses respectively. This means that the RL loss of π k+1 will approach the smallest loss possible within Π whenever π RL k+1 is indeed similar to the expert policy ξ, i.e., when

L E α (µ, π RL k+1
) is small (similar for the third statement). As expected, the behavior of DPID strictly depends also on the parameter α. First, the relevance of different estimation errors in all depends on α, so that for, e.g., small α, RL errors are more relevant. Similarly, whenever α is small, π k+1 tends to have a performance very similar to the best approximation of the greedy policy (i.e., π RL k+1 ), while when α is large π k+1 will approach the performance of π E .

Remark (comparison with APID). While the periteration error reported in [START_REF] Kim | [END_REF] is measured w.r.t. the policy Bellman error, here we consider different losses. While the estimation errors are similar in their dependency on N RL and N E and α, the rest of the bounds cannot be easily compared. In fact, in [START_REF] Kim | [END_REF] the impact of terms such as α, the similarity of ξ and π k , and the approximation space Π on the performance of APID is not clearly understandable, while in DPID, thanks to its simple integration of RL and expert data, their impact is relatively straightforward.

Before deriving a full propagation analysis for DPID, we report an intermediate lemma showing how the performance loss of π k+1 (i.e., the difference w.r.t. the optimal policy) can be upper-bounded simultaneously by the quality of the expert policy and the standard propagation of error. Lemma 1. Let V E be the expert value function and P π the transition kernel induced by any policy π, then at the end of iteration k + 1, we have the state-wise bounds

V -V π k+1 ≤ γP π (V -V π k )+(I -γP π k+1 ) -1 RL π k ( π k+1 ), V -V π k+1 ≤ V -V E + 2R max 1 -γ (I -γP ξ ) -1 E ( π k+1 ).
This lemma shows that the performance loss of DPID is related to two terms. The first (standard) term links the performance of π k+1 with the previous iteration and accounts for the RL loss, measuring how well π k+1 approximates the greedy policy G( π k ). The second term reveals that the performance is also related to the performance of the expert and the actual difference between π k+1 and ξ. Thus, the actual performance loss of π k+1 will be the minimum between these bounds, suggesting that at any iteration k DPID will benefit from either the quality of the expert or the policy improvement depending on which of the two performs the best. This means that whenever the expert is nearly optimal, π k+1 is likely to perform as well as ξ plus the expert loss. On the other hand, if the expert is suboptimal, the standard approximate policy iteration guarantee holds.

To bound the performance loss bound over iterations, we introduce standard assumptions on the concentrability coeffi-cients of the MDP and distributions ρ and µ, w.r.t. to the target distribution η (see [START_REF] Lazaric | [END_REF] for a discussion).

Assumption 1. For any policy π and non-negative integers s and t, there exists a constant C η,ρ (s, t) such that η(P ) s (P π ) t ≤ C η,ρ (s, t)ρ and C

η,ρ = (1 - γ) 2 ∞ s=0 ∞ t=0 γ s+t C η,ρ (s, t) < ∞.
For the expert policy ξ, there exists C E η,µ (s, t) such that η(P )

s (P ξ ) t ≤ C E η,µ (s, t)µ and C E η,µ = ∞ t=0 γ s+t C E η,µ (s, t) < ∞. Theorem 2. Let ∆ k = V -V π k and ∆ E = V -V E
be the performance loss of the algorithm and the expert respectively. Under Assumption 1, if DPID is run over K iterations then

∆ K 1,η ≤ min A K , B K , C K , with A K = ∆ E 1,η + 2R max (1 -γ) 2 C η,µ L E (µ, π K ), B K = min j γ K-j ∆ E 1,η + 2γ K-j (1 -γ) 2 C E η,µ L E (µ, π j ) + C η,ρ (1 -γ) 2 max l=j,...,K-1 L RL π l (ρ, π l+1 ) , C K = γ K 1 -γ + C η,ρ (1 -γ) 2 max k=1,...,K-1 L RL π l (ρ, π k+1 ).
Remark (propagation bound). The performance loss is the minimum over three terms. The term C K is the same as in the propagation error of DPI and it depends on the largest policy improvement error over iterations. The term A K depends on the quality of the expert and the expert loss of the last policy π K . The most interesting term is B K , which gives an important insight on how DPID can actually benefit from the repeated balancing of expert and RL data. In fact, whenever there is an iteration j in which the combined RL and expert losses are small, this will lead to an overall small performance loss. Finally, combining Thm. 1 and 2, we can obtain the global finite-sample bound. The resulting bound (omitted for the sake of space) shows the critical impact of α and the quality of ξ in the final performance. In fact, whenever the ξ is optimal (i.e. ∆ E = 0) and α is large, the second statement in Thm. 1 guarantees that DPID performs almost as well as the policy π E which approximates ξ the best within the policy space Π (term A K is the bound). On the other hand, if ξ is highly suboptimal and α is small, then DPID would approach the performance of DPI (term C K in the bound). In all other intermediate cases, DPID could still effectively take advantage of both the DPI process and the imitation of ξ to achieve a nearly optimal performance (term B K in the bound).

Experiments

We compare the performance of DPID 2 to DPI, which uses only RL data, and to DAGGER [START_REF] Ross | [END_REF] when the expert can be queried at arbitrary states or pure imitation learning when batch expert data are available. We also compare to 2 The implementation of DPID is available at https://www.dropbox.com/s/jj4g9ndonol4aoy/dpid.zip?dl=0 APID (implemented using CVX [START_REF] Grant | [END_REF]) and RLED (implemented using decision trees as weak learners).

Garnet domain. The Garnet framework [START_REF]On the generation of markov decision processes[END_REF]Bhatnagar et al., 2009;[START_REF] Piot | [END_REF] allows to generate random finite MDPs characterized by 3 parameters: number of states N s , number of actions N a , and maximum branching factor N b . We construct the transition function p(x, a, x ) to mimic dynamical systems where an action moves the state to states that are close by. Given a state x, next states are generated by sampling N b times from a Gaussian distribution centered at x. The probabilities of going to each of the next states are chosen uniformly at random from [0, 1] and then normalized to 1. For every action a, we construct the reward r(x, a) by sampling at random 20% of the states and randomly assigning them rewards of 1 or -1. Everywhere else the reward is 0. For such MDPs, we compute π and V using exact policy iteration and we measure the performance of the policy π alg returned by any of the algorithms as the normalized average error

Err alg = E[V * -V π alg ] E[V * ]
. In the following, we use N s = 15, N a = 3, N b (s) ∈ [3, 6] and we estimate Err alg over 100 independent runs, while error bars are computed as 95% Gaussian confidence intervals.

We start with a set of experiments where the exact state representation is used (i.e., no approximation error). This means that all algorithms have the potential to reach a zero error as the number of samples increases thus avoiding confounding effects due to the specific choice of feature design. In Fig. 1A we compare DPID with DPI and pure imitation learning (implemented with SVM [Chang and [START_REF] Chang | LIBSVM: A library for support vector machines[END_REF]) when the expert is optimal and number of demonstrations is small (N E = 15) and fixed, while the RL samples are increased at every iteration by 50. While DPI and DPID both perform better than imitation learning, DPID achieves the best performance thanks to the positive bias provided by the expert demonstrations. As expected, the greater the weight of ex- pert demonstrations the stronger the bias in initial iterations, although they later converge to the same performance. In Fig. 1B and C, we replace the optimal expert with a suboptimal one by sampling random non-optimal actions 25% and 50% of the time respectively. Since α rules the trade-off between RL and expert losses, we can see that depending on α, DPID performs either better (B, α small), similarly (C, α small), or worse (B and C, α large) than DPI. In Fig. 2A we fix the total number of samples N = N RL + N E = 1500, but change the proportion of expert to RL samples, starting from N E = 50 up to 1000. As the number of expert samples increases, DPID and DAGGER converge to the optimal policy, however DPI reaches a plateau due to the variance caused by the limited number of samples. In Fig. 2B and C we replace the optimal expert with a suboptimal one. All algorithms have errors that decrease with increasing expert samples. However, DAGGER performs the worst and DPID with a small α converges the best. In B the expert is good enough to guide towards the optimal policy, while in C it induces a bias and slows down the convergence. We now introduce an approximate representation such that for every state, we construct a binary feature vector of length d = 6 < N s . The number of ones in the representation is set to l = 3 and their locations are chosen randomly as in [Bhatnagar et al., 2009]. Fig. 3A compares the performance of DPI and DPID with or without features in the same setting as Fig. 1A. We notice that DPI is significantly affected by the bias introduced by the features, but DPID is barely so and converges to the optimal policy. In B we compare LSPI (which corresponds to APID with α = 0) and APID. Here both algorithms are affected by the bias, but APID is still better than LSPI in both cases. In C we show the performance of RLED, which does not use our representation since it constructs features automatically even in the exact representation case. We note the different effect that feature design had on the different algorithms, with DPID being the most robust. Vehicle Brake Control domain. In this domain, the state space is continuous and consists of a 4-dimensional vector representing the current velocity, the acceleration pedal value, the brake pedal value, and the target velocity. The goal of the agent is to reach the target velocity and maintain it. There are 5 actions available: acceleration up, acceleration down, brake up, brake down, do nothing. The reward equals -10 times the absolute difference between the current and the target velocity. As in [START_REF] Kim | [END_REF] the expert is implemented using the dynamics of the pedal pressure and the velocity output with the addition of random noise. We compare the performances of DPID, APID, and RLED. As in [START_REF] Kim | [END_REF] we use Radial Basis Functions to represent the state in APID. Fig. 4A and D are the DPID experiments and show the results of using an optimal and suboptimal expert respectively. This domain turns out to be easy for our classificationbased algorithm which converges quickly to a very good policy. In D, a larger α hurts the performance possibly because the expert biases the sampling process and renders RL useless. B and E show the results of APID, which seems to be less suited to this problem and converge very slowly. Expert samples help in both cases, with larger α performing better. Finally, C and F are the results of RLED. The performance is similar to that of APID which is unsurprising because they are both estimating the same quantity. Surprisingly, a larger α leads to worse performance in both cases. Although the comparison among these methods is partially biased by the difference in the underlying learning schemes (classificationbased vs projected Bellman residual minimization vs optimal Bellman residual minimization), these results show that DPID is very effective in integrating RL and expert data and it is relatively robust to expert sub-optimality and poor choice of α.

Conclusion

In this paper, we introduced DPID that naturally integrates RL and expert data. We provided a full analysis of the propa- gation error which resulted in insightful bounds. In particular, the impact of α suggests that a natural improvement on DPID would be to adapt α over iterations depending on the quality of the current policy compared to the expert policy. A principled solution may take inspiration from conservative policy iteration [START_REF] Kakade | Sham Kakade and John Langford. Approximately optimal approximate reinforcement learning[END_REF], where monotonic policy improvement is guaranteed at each iteration. Finally, we showed empirically that even a few or suboptimal expert demonstrations can lead to significant improvements and that DPID compares favorably to APID and RLED.

A Proofs

We first introduce two technical lemma whose proofs follow exactly the same steps as in Lemma 1 of [START_REF] Lazaric | [END_REF].

Lemma 2. Let {xi} N RL i=1 be a rollout set sampled from ρ. If in each state we simulate M truncated rollouts, then (2) where all = 1 RL + 1 E + 2 RL . The last inequality holds for any policy π, in particular π RL k+1 or π E k+1 which completes the proof. Proof of Lemma 1. The second statement is proved in Section 4.2 of [START_REF] Lazaric | [END_REF], while the first statement is obtained from

P sup π∈Π 1 NRL N RL i=1 1 M M j=1 R π k j (xi, π(xi)) - 1 NRL N RL i=1 Q π k (xi, π(xi)) > 2 RL ≤
V ξ -V π k+1 = T ξ V ξ -T ξ V π k+1 + T ξ V π k+1 -V π k+1 = γP ξ (V ξ -V π k+1 ) + T ξ V π k+1 -V π k+1 = (I -γP ξ ) -1 (T ξ V π k+1 -V π k+1 ) = (I -γP ξ ) -1 (Q π k+1 (ξ) -Q π k+1 (π k+1 )) ≤ 2Rmax 1 -γ (I -γP ξ ) -1 I(ξ = π k+1 ) = 2Rmax 1 -γ (I -γP ξ ) -1 E (π k+1 ).
Proof of Theorem 2. We rewrite the right-hand side of the first statement of Lem. 1 as

A k = (V -V E ) + 2Rmax 1 -γ (I -γP ξ ) -1 E ( π k ),
and we introduce E k = (I -γP π k+1 ) -1 and ∆ k+1 = V -V π k+1 . Then Lem. 1 can be written in a more compact way as

∆ k+1 ≤ min A k+1 ; γP ∆ k + E k RL π k ( π k+1
) . Now we can apply Lem. 1 recursively and obtain

∆ k+1 ≤ min A k+1 ; γP * min A k ; γP ∆ k-1 + E k-1 RL π k-1 ( π k ) + E k RL π k ( π k+1 ) = min A k+1 ; γP * A k + E k RL π k ( π k+1 ); (γP * ) 2 ∆ k-1 + γP * E k-1 RL π k-1 ( π k ) + E k RL π k ( π k+1 ) .
At each application of Lem. 1, we generate additional elements similar to the second term in the previous equation, and thus As a direct application of the previous inequality, we bound the absolute value of the performance loss of DPID after K iterations and then take the η-norm, thus concluding the proof.

Figure 1 :

 1 Figure 1: Panel A: Fixed NE 15 optimal expert demonstrations and increasing NRL by 50 at each iteration, starting with NRL = 50-NE (NE = 0 for DPI). Panel B and C: Same setting but expert with non-optimal actions with probability 0.25 (B) or 0.5 (C).

Figure 2 :

 2 Figure 2: N = 1500. Panel A: NRL = 1500 -NE for DPID with NE optimal demonstrations increases by 50 at each iteration. NRL = 1500 for DPI. Panel B and C: Same setting but expert with non-optimal action with probability 0.25 (B) or 0.5 (C).

Figure 3 :

 3 Figure 3: Fixed NE = 50 optimal expert samples and NRL increasing over iterations by 500, starting with NRL = 500-NE (NE = 0 for DPI and API). Panel A: DPI and DPID with exact and approximate representation. Panel B-C: Same for API, APID, and RLED.

Figure 4 :

 4 Figure 4: Panel A: DPI versus DPID when fixing NE = 15 optimal expert demonstrations and increasing NRL by 500 at every iteration. Panel D: DPI versus DPID when fixing NE = 100 suboptimal expert demonstrations and increasing the RL data by 500 at iteration (0.5 prob. of wrong action). Panels B/E and C/F: Same as A and D for API/APID and RLED respectively.

∆

  k+1 ≤ min A k+1 ; min j=1,...,k (γP * ) k+1-j Aj + k l=j (γP * ) k+1-l E l RL π l (π l+1 ) ; (γP * ) k+1 ∆0 + k l=0 (γP * ) k-l E l RL π l (π l+1 ) .

  δ,with 2 RL = 8(1-γ H ) √ RL (h log ( eN RL h ) + log ( 8 δ )). Lemma 3. Let {zj} N Ej=1 be a set sampled from µ, thenP sup π∈Π |L E (μ, π) -L E (µ, π)| > E ≤ δ,(1)with E = 8 N E (h log( eN E h ) + log( 4 δ )). LRL π k (ρ, πk+1 ) + αL E π k (μ, πk+1 ) + 1 RL + 1 E + 2 RL ≤ LRL π k (ρ, π) + αL E π k (μ, π) + 1 RL + 1 E + 2 RL ≤ L RL π k (ρ, π) + αL E π k (µ, π) + 2 all

	Proof of Theorem 1.
	L RL π k (ρ, πk+1 ) + αL E π k (µ, πk+1 )
	≤

M N RL

Qmax 2(h log( eM N RL h )+log(

8

δ )) and

P sup π∈Π |L RL π k (ρ, π) -L RL π k (ρ, π)| > 1 RL ≤ δ, with 1 RL = 16Qmax 2 N

In order to simplify the theoretical analysis we will consider the case |A| =

throughout the rest of the paper.

Acknowledgments

This work was supported by the French Ministry of Higher Education and Research and the French National Research Agency (ANR) under project ExTra-Learn n.ANR-14-CE24-0010-01.