Angela Felisia

Chiarello

Paola Goatin

Luis Miguel Villada

High-order Finite Volume WENO schemes for non-local multi-class traffic flow models

This paper focuses on the numerical approximation of a class of non-local systems of conservation laws in one space dimension, arising in traffic modeling, proposed by [F. A. Chiarello and P. Goatin. Non-local multi-class traffic flow models. Networks and Heteroge-neous Media, to appear, Aug. 2018]. We present the multi-class version of the Finite Volume WENO (FV-WENO) schemes [C. Chalons, P. Goatin, and L. M. Villada. High-order numerical schemes for one-dimensional non-local conservation laws. SIAM Journal on Scientific Computing, 40(1):A288-A305, 2018.], with quadratic polynomial reconstruction in each cell to evaluate the non-local terms in order to obtain high-order of accuracy. Simulations using FV-WENO schemes for a multi-class model for autonomous and human-driven traffic flow are presented for M = 3 .

Introduction

We consider the following class of non-local systems of M conservation laws in one space dimension, introduced in [START_REF] Chiarello | Non-local multi-class traffic flow models[END_REF] to model multi-class traffic dynamics:

∂ t ρ i (t, x) + ∂ x ρ i (t, x)v i ((r * ω i )(t, x)) = 0, i = 1, ..., M, (1.1)
where

r(t, x) := M i=1 ρ i (t, x), (1.2)
v i (ξ) := v max i ψ(ξ), (1.3)
(r * ω i )(t, x) := x+η i x r(t, y)ω i (y -x) dy , (1.4)
where ρ i is the density of vehicles belonging to the i-th class, v i is the class-specific mean velocity and η i is proportional to the look-ahead distance.

We assume that the following hypotheses hold:

1 Inria Sophia Antipolis -Méditerranée, Université Côte d'Azur, Inria, CNRS, LJAD, 2004 route des Lucioles -BP 93, 06902 Sophia Antipolis Cedex, France. E-mail: {felisia.chiarello, paola.goatin}@inria.fr 2 GIMNAP-Departamento de Matemáticas, Universidad del Bío-Bío, Casilla 5-C, Concepción, Chile and CI 2 MA, Universidad de Concepción, Casilla 160-C, Concepción, Chile. E-mail: lvillada@ubiobio.cl (H1) The convolution kernels ω i ∈ C 1 ([0, η i ]; R + ), η i > 0, are non-increasing functions with interaction strength J i := η i 0 ω i (y) dy. We set W 0 := max i=1,...,M ω i (0). (H2) v max i are the maximal velocities, with 0

< v max 1 ≤ v max 2 ≤ . . . ≤ v max M .
(H3) ψ : R + → R + is a smooth non-increasing function such that ψ(0) = 1 and ψ(r) = 0 for r ≥ 1.

We couple (1.1) with an initial datum

ρ i (0, x) = ρ 0 i (x), i = 1, . . . , M. (1.5)
Model (1.1) is a generalization of the n-population model for traffic flow described in [START_REF] Benzoni-Gavage | An n-populations model for traffic flow[END_REF] and it is a multi-class version of the one dimensional scalar conservation law with non-local flux proposed in [START_REF] Blandin | Well-posedness of a conservation law with non-local flux arising in traffic flow modeling[END_REF]. The term "non-local" refers to the speed functions v i evaluated on a neighborhood of x ∈ R defined by the downstream convolution between the weight functions ω i and the sum of the densities r. This is intended to describe the reaction of drivers that adapt their velocity to the downstream traffic, assigning greater importance to closer vehicles, see also [START_REF] Chiarello | Global entropy weak solutions for general non-local traffic flow models with anisotropic kernel[END_REF][START_REF] Goatin | Well-posedness and finite volume approximations of the LWR traffic flow model with non-local velocity[END_REF]. We consider different anisotropic discontinuous kernels for each equation of the system. The model takes into account the distribution of heterogeneous drivers and vehicles characterized by their maximal speeds and look-ahead visibility in a traffic stream. It is worth to point out that in multi-class dynamic faster vehicles can overtake slower ones and slower vehicles slow down the faster ones, avoiding one of the biggest limitations of the standard LWR traffic flow model [START_REF] Lighthill | On kinematic waves. II. A theory of traffic flow on long crowded roads[END_REF][START_REF] Richards | Shock waves on the highway[END_REF], i.e. the first-in-first-out rule. The computation of numerical solutions for (1.1) is challenging due to the high nonlinearity of the system and the dependence of the flux function on integral terms. First and second order finite volume schemes for (1.1) were proposed and analyzed in [START_REF] Chiarello | Non-local multi-class traffic flow models[END_REF][START_REF] Chiarello | Lagrangian-Antidiffusive Remap schemes for non-local multi-class traffic flow models[END_REF]. In this paper, a high-order finite-volume WENO (FV-WENO) scheme is proposed to solve the nonlocal multi-class system (1.1). The procedure proposed in [START_REF] Chalons | High-order numerical schemes for one-dimensional nonlocal conservation laws[END_REF] is used and extended to the multi-class cases in order to evaluate the non-local term that appears in the flux functions.

The paper is organized as follows. First, in Section 2, we describe the implementation of the high-order FV-WENO scheme for the non-local system (1.1). In Section 3, we provide a couple of numerical test in the case of three populations (M = 3) and convergence studies for third, fifth and seventh accuracy order.

Finite Volume WENO schemes

In this section, we solve the non-local system of conservation laws (1.1) by using a high-order finite volume WENO scheme [START_REF] Shu | Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws[END_REF][START_REF] Shu | Efficient implementation of essentially non-oscillatory shock-capturing schemes[END_REF]. First we consider {I j } M j=1 as a partition of [-L, L] and the points x j are the center of the cells

I j = [x j-1 2 , x j+ 1 2 ], with length |I j | = ∆x = 2
M . We denote the unknowns by ρ i,j (t), the cell average of the exact solution ρ i (t, •) in the cell I j :

ρ i,j (t) := 1 ∆x I j ρ i (t, x)dx.
We extend ω i (x) = 0 for x > η i , and set

ω k i := 1 ∆x k ∆x (k-1) ∆x ω i (x) dx , k ∈ N * , (2.1)
so that ∆x +∞ k=1 ω k i = η i 0 ω i (x) dx = J i (the sum is indeed finite since ω k i = 0 for k ≥ N i sufficiently large). Moreover, we set r j (t) := m i=1 ρ i,j (t) and define the convolution term in the form R i (t, x) := (r * ω i )(t, x). Integrating (1.1) over I j we obtain

d dt ρ i,j (t) = - 1 ∆x f i (t, x j+1/2 ) -f i (t, x j-1/2 ) , i = 1, . . . , M, ∀j ∈ Z,
where

f i (t, x j+1/2 ) := ρ i (t, x j+ 1 2 )v i (R i (t, x j+ 1 2
)). This equation is approximated by the semidiscrete conservative scheme

d dt ρ i,j (t) = - 1 ∆x f i,j+ 1 2 -f i,j-1 2 , i = 1, . . . , M, ∀j ∈ Z, (2.2)
where

f i,j+ 1 2 is a consistent approximation of flux ρ i v i (R i ) at interface x j+1/2 .
Here, we consider the multi-class version of the Godunov scheme [START_REF] Chiarello | Non-local multi-class traffic flow models[END_REF]

f i,j+ 1 2 := f (ρ l i,j+ 1 2 , ρ r i,j+ 1 2 ) = ρ l i,j+ 1 2 v i (R r i,j+1/2 ), (2.3)
where ρ l i,j+ ) obtained from the cell averages ρ i,j (t) j∈Z . In this work, we consider the classical WENO scheme proposed in [START_REF] Shu | Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws[END_REF][START_REF] Shu | Efficient implementation of essentially non-oscillatory shock-capturing schemes[END_REF]. R r i,j+1/2 is the right approximation of R i (t, x) at the interface x j+1/2 . Since R i is defined by a convolution, we naturally set

R r i,j+1/2 = R i (t, x j+1/2 ) := R i,j+1/2 (t).
In order to compute the integral R i,j+1/2 , we use the technique proposed in [START_REF] Chalons | High-order numerical schemes for one-dimensional nonlocal conservation laws[END_REF], i.e., we consider a reconstruction of ρ i (x, t) on I j by taking advantage of the high-order WENO reconstructions ρ r i,j-1 2 and ρ l i,j+ 1 2 at the boundaries of I j , as well as the approximation of the cell average ρ n i,j . We consider a quadratic polynomial p i,j (x) defined on I j such that

p i,j (x j-1 2 ) = ρ r i,j-1 2 , p i,j (x j+ 1 2 ) = ρ l i,j+ 1 2 , 1 ∆x I j p i,j (x)dx = ρ n i,j .
In particular, we take

p i,j (x) := a i,j,0 v (0) (ξ j (x)) + a i,j,1 v (1) ξ j (x) + a i,j,2 v (2) ξ j (x) , x ∈ I j , (2.4)
with

v (0) (y) = 1, v (1) (y) = y, v (2) (y) = 1 2 3y 2 -1 , ξ j (x) = x -x j ∆x/2 .
Coefficients in (2.4) can be easily computed as

a i,j,0 = ρ n i,j , a i,j,1 = 1 2 ρ l i,j+ 1 2 -ρ r i,j-1 2 , a i,j,2 = 1 2 ρ l i,j+ 1 2 + ρ r i,j-1 2 -ρ n i,j .
Now, summing for i = 1, . . . , M , we have

P j (x) := M i=1 p i,j (x) = âj,0 v (0) (ξ j (x)) + âj,1 v (1) ξ j (x) + âj,2 v (2) ξ j (x) , x ∈ I j , with âj,0 := M i=1 a i,j,0 = r j , âj,1 := M i=1 a i,j,1 , âj,2 := M i=1 a i,j,2 .
With this polynomial P j (x), we can compute R i,j+

1 2 as R i,j+ 1 2 = M k=1 I j+k P j+k (y)ω i (y -x j+ 1 2 )dy = N k=1 I j+k ω i (y -x j+ 1 2 ) 2 l=0 âj+k,l v (l) (ζ j+k (y))dy = M k=1 2 l=0
âj+k,l

I j+k ω i (y -x j+ 1 2 )v (l) (ζ j+k (y))dy = M k=1 2 l=0 âj+k,l ∆x 2 1 -1 ω i ∆x 2 y + (k - 1 2 )∆x v (l) (y)dy = M k=1 2 l=0 âj+k,l Γ i,k,l ,
where the coefficients Γ i,k,l are computed exactly or using a high-order quadrature approximation.

The utilization of the quadratic polynomial on each cell to evaluate the convolution term suggests the following algorithm to approach the solution of non-local system (1.1):

Algorithm: FV-WENO scheme for non-local multi-class traffic models Given ρ n i,j for j ∈ Z, i = 1, . . . , M , approximation of the cell averages of ρ i (x, t) at t n . 4. Use a high-order accurate Runge-Kutta method to solve the semi-discrete system (2.2), with the CFL condition ∆t ∆x

v max M ψ ∞ ≤ 1 2 . ( 2.5)
In this paper, we use the WENO method of third (WENO3), fifth (WENO5) and seventh (WENO7) accuracy order proposed by [START_REF] Shu | Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws[END_REF][START_REF] Shu | Efficient implementation of essentially non-oscillatory shock-capturing schemes[END_REF]. For the temporal discretization, in order to match the order of spatial accuracy, fifth or seventh explicit Runge-Kutta schemes are used [START_REF] Buchmüller | Improved accuracy of high-order weno finite volume methods on cartesian grids[END_REF].

Numerical tests

In the following numerical tests, we solve (1.1) numerically in the intervals x ∈ [-1, 1] and t ∈ [0, 2]. We propose two tests in order to illustrate the dynamics of the model (1.1) for autonomous and human-driven vehicles, using FV-WENO5 scheme with 1/∆x = 400. For each integration, we set ∆t to satisfy the CFL condition (2.5).

To test the accuracy order of the proposed method, since we cannot compute the exact solution explicitly, we use a reference solution ρ ref obtained using FV-WENO7 on a refined mesh (1/∆x = 6400). The L 1 -error for the cell average is given by

L 1 (∆x) = M i=1   1 N N j=1 |ρ i,j -ρ ref i,j |   ,
where ρ i,j and ρ ref i,j are the cell averages of the numerical approximation and the reference solution respectively. The Experimental Order of Accuracy (E.O.A.) is naturally defined by

γ(∆x) = log 2 L 1 (∆x)/L 1 (∆x/2) .

Test 1, circular road

The aim of this test is to study the possible impact of the presence of Connected Autonomous Vehicles (CAVs) on road traffic performances, as proposed in [START_REF] Chiarello | Global entropy weak solutions for general non-local traffic flow models with anisotropic kernel[END_REF]Section 4.2]. Let us consider a circular road modeled by the space interval [-1, 1] with periodic boundary conditions at x = ±1. The interaction radius of CAVs is much grater than the one of human-driven cars. Moreover, we can assign a constant convolution kernel to CAVs, since we assume that the information they get about surrounding traffic is transmitted through wireless connections and its degree of accuracy does not depend on distance. We consider the following initial data and parameters

ρ 1 (0, x) = αp(x), ω 1 (x) = 1 η 1 , η 1 = 0.3, v max 1 = 0.8, (3.1) ρ 2 (0, x) = β p(x), ω 2 (x) = 1 η 2 , η 2 = 0.3, v max 2 = 1.2, (3.2)
ρ 3 (0, x) = γ p(x), ω 3 (x) = 2 η 3 1 - x η 3 , η 3 = 0.05, v max 3 = 1.2, (3.3)
where p(x) = 0.5 + 0.3 sin(5πx) is the total initial density, α, β, γ ≥ 0 and α + β + γ = 1. Above, ρ 1 represents the density of autonomous trucks, ρ 2 is the density of autonomous cars and ρ 3 is the density of human-driven cars. In Figure 1a we consider the penetration rates

α = 0.5, β = 0.3, γ = 0.2,
and we can compare the total density r = ρ 1 + ρ 2 + ρ 3 with that one in Figure 1b where we have no human-driven cars: α = 0.5, β = 0.5, γ = 0.

We observe that oscillations are reduced if only autonomous vehicles are present. Finally, we compute the E.O.A. for the FV-WENO schemes. We consider parameters α = 0.5, β = 0.3, γ = 0.2 , and compute the L 1 -error at T = 0.2 in Table 1. As expected, we obtain the correct order. Figure 1: (t,x)-plots of the total density r(t, x) = ρ 1 (t, x) + ρ 2 (t, x) + ρ 3 (t, x) computed with the FV-WENO5 scheme, corresponding to different penetration rates of autonomous and nonautonomous vehicles: (a) mixed autonomous / human-driven traffic, (b) fully autonomous traffic.

FV-WENO3

FV-WENO5 FV-WENO7 3), with α = 0.5, β = 0.3, γ = 0.2 and final time T = 0.2. The reference solution is computed with FV-WENO7 scheme for 1/∆x = 6400.

Test 2, stretch of straight road

In this test case, we consider a stretch of road populated by cars and trucks as in the example proposed in [START_REF] Chiarello | Non-local multi-class traffic flow models[END_REF]Section 4.2]. The space domain is given by the interval [-1, 1] and we impose absorbing conditions at the boundaries. The dynamics is described by the equation (1.1) with M = 3, and the following initial conditions and parameter values

ρ 1 (0, x) = 0.5χ [-0.6,-0.1] (x), ω 1 (x) = 2 η 1 1 - x η 1 , η 1 = 0.1, v max 1 = 0.8, (3.4) ρ 2 (0, x) = α 1 χ [-0.9,-0.6] (x), ω 2 (x) = 1 η 2 , η 2 = 0.5, v max 2 = 1.3. (3.5) ρ 3 (0, x) = β 1 χ [-0.9,-0.6] (x), ω 3 (x) = 2 η 3 1 - x η 3 , η 3 = 0.05, v max 2 = 1.3. (3.6)
In this setting, ρ 1 (t, x) describes the density of human-driven trucks, ρ 2 (x, t) the density of autonomous cars and ρ 3 (x, t) is density of human driven cars. We have a red traffic light located at x = -0.1, which turns green at the initial time t = 0.

In Figure 2a and we can compare the space-time evolution of the total density r = ρ 1 + ρ 2 + ρ 3 with the one in Figure 2b, where α 1 = 0, β 1 = 0.5.

In this case, the presence of autonomous cars in a heterogeneous traffic of human-driven vehicles induces higher vehicle densities during the overtaking phase, but for shorter time. In Figure 3 we display the density profiles of ρ 1 , ρ 2 and ρ 3 computed with different FV-WENO schemes at time t = 0.5 in the same setting of Test 2(a). We can appreciate the efficiency of FV-WENO schemes in presence of discontinuities in comparison with the finite volume Godunov type scheme.

Conclusions

In this paper, we applied high-order finite volume WENO schemes to the non-local multi-class traffic flow model proposed in [START_REF] Chiarello | Non-local multi-class traffic flow models[END_REF]. We used quadratic polynomial reconstructions in each cell to evaluate the non-local terms in order to obtain high order of accuracy. The numerical results of the accuracy test show that the proposed schemes maintain the correct order of accuracy. Besides, the considered examples allow to illustrate the interaction dynamics of mixed traffic consisting of both autonomous and human-driven vehicles.

1. Compute ρ l i,j+ 1 2 and ρ r i,j+ 1 2 , 1 2

 121 the left and right high-order WENO approximations for j ∈ Z and i = 1, . . . , M ; 2. Calculate R i,j+ for j ∈ Z and i = 1, . . . , M ; 3. Calculate the Godunov numerical flux (2.3) for j ∈ Z and i = 1, . . . , M ;

  (a) α = 0.5, β = 0.3, γ = 0.2, (b) α = 0.5, β = 0.5, γ = 0.

  we consider the rates α 1 = 0.25, β 1 = 0.25, (a) α 1 = 0.25, β 1 = 0.25, (b) α 1 = 0, β 1 = 0.5.

Figure 2 :Figure 3 :

 23 Figure2: (t, x)-plots of the total density r(t, x) = ρ 1 (t, x) + ρ 2 (t, x) + ρ 3 (t, x) computed with FV-WENO5 scheme, corresponding to different penetration rates of autonomous and non-autonomous vehicles.

Table 1 :

 1 E.O.A. Test 1, initial condition (3.1)-(3.

	1/∆x	L 1 -err	γ(∆x)	L 1 -err	γ(∆x)	L 1 -err	γ(∆x)
	100	1.51e-03	-	1.09e-04	-	5.64e-05	-
	200	1.38e-04	3.44	9.44e-06	3.53	1.54e-06	5.19
	400	1.20e-05	3.53	4.01e-07	4.56	1.58e-08	6.61
	800	1.27e-06	3.24	1.26e-08	4.99	1.68e-10	6.55
	1600	1.05e-07	3.01	3.60e-10	5.12	4.71e-12	5.15

Acknowledgements

This research was supported by the Inria Associated Team Efficient numerical schemes for non-local transport phenomena (NOLOCO; 2018-2020). LMV is supported by Fondecyt project 1181511 and CONICYT/PIA/Concurso Apoyo a Centros Científicos y Tecnológicos de Excelencia con Financiamiento Basal AFB170001.


Jean-Bastien Grill
email: jbgrill@google.com

Michal Valko
email: michal.valko@inria.frmunos@google.com

Rémi Munos

Optimistic Optimization of a Brownian

 which provided only polynomial rates. Our algorithm is adaptive-each query depends on previous values-and is an instance of the optimism-in-the-face-of-uncertainty principle.

Introduction to sample-efficient Brownian optimization

We are interested in optimizing a sample of a standard Brownian motion on [0, 1], denoted by W. More precisely, we want to sequentially select query points t n ∈ [0, 1], observe W (t n ), and decide when to stop to return a point t and its value M = W t in order to well approximate its maximum M sup t∈[0,1] W (t). The evaluations t n can be chosen adaptively as a function of previously observed values W (t 1 ), ..., W (t n-1 ). Given an ε > 0, our goal is to stop evaluating the function as early as possible while still being able to return t such that with probability at least 1 -ε, M -W t ≤ ε. The number of function evaluations required by the algorithm to achieve this ε-approximation of the maximum defines the sample-complexity.

Motivation There are two types of situations where this problem is useful. The first type is when the random sample function W (drawn from the random process) already exists prior to the optimization.

Either it has been generated before the optimization starts and the queries correspond to reading values of the function already stored somewhere. For example, financial stocks are stored at a high temporal resolution and we want to retrieve the maximum of a stock using a small number of memory queries. Alternatively, the process physically exists and the queries correspond to measuring it.

Another situation is when the function does not exist prior to the optimization but is built simultaneously as it is optimized. In other words, observing the function actually creates it. An application of this is when we want to return a sample of the maximum (and the location of the maximum) of a Brownian motion conditioned on a set of already observed values. For example, in Bayesian optimization for Gaussian processes, a technique called Thomson sampling [START_REF] William | On the likelihood that one unknown probability exceeds another in view of the evidence of two samples[END_REF][START_REF] Chapelle | An empirical evaluation of Thompson sampling[END_REF][START_REF] Russo | A tutorial on Thompson sampling[END_REF][START_REF] Basu | Analysis of Thompson sampling for Gaussian process optimization in the bandit setting[END_REF] requires returning the maximum of a sampled function drawn from the posterior distribution. The problem considered in the present paper can be seen as a way to approximately perform this step in a computationally efficient way when this Gaussian process is a Brownian motion.

Moreover, even though our algorithm comes from the ideas of learning theory, it has applications beyond it. For instance, in order to computationally sample a solution of a stochastic differential equation, [START_REF] Hefter | Optimal strong approximation of the one-dimensional squared Bessel process[END_REF] express its solution as a function of the Brownian motion W and its running minimum. They then need, as a subroutine, an algorithm for the optimization of Brownian motion to compute its running minimum. We are giving them that and it is light-speed fast. Prior work [START_REF] Al-Mharmah | Optimal random non-adaptive algorithm for global optimization of Brownian motion[END_REF] provide a non-adaptive method to optimize a Brownian motion. They prove that their method is optimal among all non-adaptive methods and their sample complexity is polynomial of order 1/ √ ε. More recently, [START_REF] Calvin | Adaptive approximation of the minimum of Brownian motion[END_REF] provided an adaptive algorithm with a sample complexity lower than any polynomial rate showing that adaptability to previous samples yields a significant algorithmic improvement. Yet their result does not guarantee a better rate than a polynomial one.

Our contribution

We introduce the algorithm OOB = optimistic optimization of the Brownian motion. It uses the optimism-in-face-of-uncertainty apparatus: Given n -1 points already evaluated, we define a set of functions U n in which W lies with high probability. We then select the next query point t n where the maximum of the most optimistic function of U n is reached: t n arg max t∈[0,1] max f ∈Un f (t). This begets a simple algorithm that requires an expected number of queries of the order of log 2 (1/ε) to return an ε-approximation of the maximum, with probability at least 1 -ε w.r.t. the random sample of the Brownian motion. Therefore, our sample complexity is better than any polynomial rate.

Solving an open problem [START_REF] Munos | Optimistic optimization of deterministic functions without the knowledge of its smoothness[END_REF] provided sample complexity results for optimizing any function f characterized by two coefficients (d, C) where d is the near-optimality dimension and C the corresponding constant (see his Definition 3.1). It is defined as the smallest d ≥ 0 such that there exists a semi-metric and a constant C > 0, such that, for all ε > 0, the maximal number of disjoint -balls of radius O(ε) with center in {x, f (x) ≥ sup x f (x) -ε} is less than Cε -d . Under the assumption that f is locally (around one global maximum) one-sided Lipschitz with respect to (see his Assumption 2), he proved that for a function f characterized by (d = 0, C), his DOO algorithm has a sample complexity of O(C log(1/ε)), whereas for a function characterized by (d > 0, C), the sample complexity of DOO is O C/ε d . Our result answers a question he raised: What is the near-optimality dimension of a Brownian-motion? The Brownian motion being a stochastic process, this quantity is a random variable so we consider the number of disjoint balls in expectation. We show that for any ε, there exists some particular metric ε such that the Brownian motion W is ε -Lipschitz with probability 1 -ε, and there exists a constant C(ε) = O(log(1/ε)) such that (d = 0, C(ε)) characterizes the Brownian motion. However, there exists no constant C < ∞ independent of ε such that (d = 0, C) characterizes the Brownian motion. Therefore, we solved this open problem. Our answer is compatible with our result that our algorithm has a sample complexity of O(log 2 (1/ε)).

New algorithm for Brownian optimization

Our algorithm OOB is a version of DOO [START_REF] Munos | Optimistic optimization of deterministic functions without the knowledge of its smoothness[END_REF] with a modified upper bound on the function, in order to be able to optimize stochastic processes. Consider the points t 1 < t 1 < ... < t n evaluated so far and t 0 = 0. OOB defines an upper confidence bound B [ti,ti+1] for each interval [t i , t i+1 ] with i ∈ {0, ..., n -1} and samples W in the middle of the interval with the highest upper-confidence bound. Algorithm 1 reveals its pseudo-code. 3 Guarantees: OOB is correct and sample-efficient

Algorithm 1 OOB algorithm 1: Input: ε 2: Init: I ← {[0, 1]}, t 1 = W (1) 3: for i = 2,
Let M sup t∈[0,1] W (t)
be the maximum of the Brownian motion, t ε the output of OOB called with parameter ε > 0, and N ε the number of Brownian evaluations performed until OOB terminates. All are random variables that depend on the Brownian motion W . We now voice our main result. Theorem 1. There exists a constant c > 0 such that for all ε < 1/2,

P M -W tε > ε ≤ ε and E[Nε] ≤ c log 2 (1/ε).
The first inequality quantifies the correctness of our estimator M ε = W t ε . Given a realization of the Brownian motion, our OOB is deterministic. The only source of randomness comes from the realization of the Brownian. Therefore, being correct means that among all possible realizations of the Brownian motion, there is a subset of measure at least 1 -ε on which OOB outputs an estimate M ε which is at most ε-away from the true maximum. Such guarantee is called probably approximately correct (PAC). The second inequality quantifies performance. We claim that the expectation (over W ) of the number of samples that OOB needs to optimize this realization with precision ε is O log 2 (1/ε) .

Corollary 1. We get the classic (δ, ε)-PAC guarantee easily. For any δ > 0 and ε > 0, choose ε = min(δ, ε) and apply Theorem 1 for ε from which we get

P M -W t ε > ε ≤ ε which is stronger than P M -W t ε > ε ≤ δ. Similarly, E[N ε ] ≤ c log 2 (1/ε ) ≤ 4c(log(1/ε) + log(1/δ)) 2 .
Remark 1. Our PAC guarantee is actually stronger than stated in Theorem 1. Indeed, the PAC guarantee analysis can be done conditioned on the collected function evaluations and get

P M -W tε > ε W (t1), ..., W (tN ε ) ≤ ε,
from which taking the expectation on both sides gives the first part of Theorem 1. This means that the unfavorable cases, i.e., the Brownian realizations for which M -M ε > ε, are not concentrated on some subsets of Brownian realizations matching some evaluations in t 1 , ..., t Nε . In other words, the PAC guarantee also holds when restricted to the Brownian realizations matching the evaluations in t 1 , ..., t Nε only. This is possible because N ε is not fixed but depends on the evaluations done by OOB.

One difference from the result of [START_REF] Calvin | Adaptive approximation of the minimum of Brownian motion[END_REF] is that theirs is with respect to the L p norm. For their algorithm, they prove that with n samples it returns

t n ∈ [0, 1] such that ∀r > 1, p > 1, ∃cr,p, E M -W ( tn) p 1/p ≤ cr,p/n r .
To express their result in the same formalism as ours, we first choose to achieve accuracy ε 2 and compute the number of samples n ε 2 needed to achieve it. Then, for p = 1, we apply Markov inequality and get that for all r > 1 there exists c r,1 such that

P M -W ( tn ε 2 ) > ε ≤ ε and Nε ≤ cr,1/ε 1/r .
On the other hand, in our Theorem 1 we give a poly-logarithmic bound for the sample complexity and we are in the business because this is better than any polynomial rate.

Analysis and the proof of the main theorem

We provide a proof of the main result. Let I fin be the set I of intervals tracked by OOB when it finishes.

We define an event C such that for any interval I of the form I = [k/2 h , (k + 1)/2 h ] with k and h being two integers where 0 ≤ k < 2 h , the process W is lower than B I on the interval I.

Definition 1. Event C is defined as

C ∞ h=0 2 h -1 k=0 sup t∈[k/2 h ,(k+1)/2 h ] W (t) ≤ B [k/2 h ,(k+1)/2 h ] •
Event C is a proxy for the Lipschitz condition on W for the pseudo-distance d(x, y) = |y -x| ln(1/|y -x|) because η ε (δ) scales with δ as δ ln (1/δ). We show that it holds with high probability. To show it, we make use of the Brownian bridge which is the process Br(t) (W (t)|W (1) = 0). Lemma 1 follows from the Markov property of the Brownian combined with a bound on the law of the maximum of Br(t) to bound the probability P[sup t∈I W (t) ≥ B I ] for any I of the form [k/2 h , (k + 1)/2 h ] and a union bound over all these intervals. Lemma 1. For any ε, event C from Definition 1 holds with high probability. In particular,

P[C c ] ≤ ε 5 .
Proof. For any interval I,

BI = max(W (a), W (b)) + ηε(b -a) (by definition of BI ) = W (a) + W (b) 2 + |W (a) -W (b)| 2 + ηε(b -a) (max(x, y) = (x+y+|x-y|)/2) ≥ W (a) + W (b) 2 + W (a) -W (b) 2 2 + (ηε(b -a)) 2 ∀x, y > 0, (x+y) 2 ≥ x 2 +y 2 .
We now define

t h,k k 2 h , h,k W (t h,k ) -W (t h,k+1 ) 2 , η h ηε(b -a),
and

A h,k sup t∈[k/2 h ,(k+1)/2 h ] W (t) > B [k/2 h ,(k+1)/2 h ] •
First, for any a < b, the law of the maximum of a Brownian bridge gives us

W (t) > x W (a) = Wa, W (b) = W b = exp - 2(x-Wa)(x-W b ) b -a •
Combining it with the definition of A h,k and the first inequality of the proof we get

P A h,k W (t h , k), W (t h , k + 1) ≤ exp   - 2 W (t h,k+1 )-W (t h,k ) 2 + ∆ 2 h,k + η 2 h W (t h,k )-W (t h,k+1 ) 2 + ∆ 2 h,k + η 2 h t h,k+1 -t h,k   = exp -2 h+1 ∆ 2 h,k + η 2 h -∆ h,k ∆ 2 h,k + η 2 h + ∆ h,k = exp -2 h+1 η 2 h = exp -2 h+1 5 2 • 2 h ln 2 h ε = ε 2 h 5 • By definition, C ∞ h=0 2 h -1 k=0 A c h,k = ∞ h=0 2 h -1 k=0 A h,k
. By union bound on all A h,k we get

P[C c ] ≤ h≥1 2 h -1 k=0 P[A h,k ] ≤ h≥1 2 h -1 k=0 ε 2 h 5 ≤ h≥1 ε 5 2 4h ≤ ε 5 .
Lemma 1 is useful for two reasons. As we bound the sample complexity on event C and the complementary event in two different ways, we can use Lemma 1 to combine the two bounds to prove Proposition 2 in the end. We also use a weak version of it, bounding ε 5 by ε to prove our PAC guarantee. For this purpose, we combine the definition of C with the termination condition of OOB to get that under event C, the best point M ε so far, is close to the maximum M of the Brownian up to ε.

Since C holds with high probability, we have the following PAC guarantee which is the first part of the main theorem.

Proposition 1. The estimator M ε = W t ε is probably approximately correct with

P M -Mε > ε ≤ ε.
Proof. Let I next = [a, b] be the interval that the algorithm would split next, if it was not terminated. Since the algorithm only splits the interval with the highest upper bound then B next = sup I∈Ifin B I . Also let I max ∈ I fin be one of the intervals where a maximum is reached, t max ∈ arg max t∈[0,1] W (t) M and t max ∈ I max . Then, on event C,

M ≤ BI max ≤ BI next = max(W (a), W (b)) + ηε(b -a).
Since the algorithm terminated, we have that η ε (b -a) ≤ ε and therefore,

max(W (a), W (b)) ≥ M -ε,
which combined with Lemma 1 finishes the proof as ε 5 ≤ ε.

In fact, Proposition 1 is the easy-to-obtain part of the main theorem. We are now left to prove that the number of samples needed to achieve this PAC guarantee is low. As the next step, we define the near-optimality property. A point t is said to be η-near-optimal when its value W (t) is close to the maximum M of the Brownian motion up to η. Check out the precise definition below.

Definition 2. When an (h, k, η) verifies W k 2 h ≥ M -η, we say that the point t = (k/2 h ) is η-nearoptimal. We define N h (η) as the number of η-near-optimal points among 0/2 h , 1/2 h , ...,

2 h /2 h , N h (η) k ∈ 0, ..., 2 h , such that W k 2 h ≥ M -η •
Notice that N h (η) is a random variable that depends on a particular realization of the Brownian motion. The reason why we are interested in the number of near-optimal points is that the points the algorithm will sample are η ε 1/2 h -near-optimal. Since we use the principle of optimism in face of uncertainty, we consider the upper bound of the Brownian motion and sample where this upper bound is the largest. If our bounds on W hold, i.e., under event C, then any interval I with optimistic bound B I < M is never split by the algorithm. This is true when C holds because if the maximum of W is reached in I max , then B Imax ≥ M > B I which shows that I max is always chosen over I. Therefore, a necessary condition for an interval [a, b] to be split is that max(W (a), W (b)) ≥ M -η which means that either a or b or both are η-near-optimal which is the key point of Lemma 2. Lemma 2. Under event C, the number of evaluated points N ε by the algorithm verifies

Nε ≤ 2 hmax h=0
N h ηε 1/2 h , with hmax being the smallest h such that ηε 1/2 h ≤ ε.

Lemma 2 explicitly links the near-optimality from Definition 2 with the number of samples N ε performed by OOB before it terminates. Here, we use the optimism-in-face-of-uncertainty principle which can be applied to any function. In particular, we define a high-probability event C under which the number of samples is bounded by the number of near-optimal points N h (η h ) for all h ≤ h max . From the termination condition of the algorithm, we know that any interval that is satisfying I = [k/2 h , (k + 1)/2 h ] with h ≥ h max will not be split during the execution.1 Therefore, another necessary condition for an interval I = [a, b] to be split during the execution is that b -a > 1/2 hmax . Writing η h η ε 1/2 h , we deduce from these two necessary conditions that

N ≤ hmax h=0 2 h -1 k=0 1 k 2 h or k + 1 2 h is η h -near-optimal ≤ hmax h=0 2 h -1 k=0 1 k 2 h is η h -near-optimal + 1 k + 1 2 h is η h -near-optimal ≤ 2 hmax h=0 2 h k=0 1 k 2 h is η h -near-optimal = 2 hmax h=0 N h (η h ).
We now prove a property specific to W by bounding the number of near-optimal points of the Brownian motion in expectation. We do it by rewriting it as two Brownian meanders [START_REF] Durrett | Weak convergence to Brownian meander and Brownian excursion[END_REF], both starting at the maximum of the Brownian, one going backward and the other one forward with the Brownian meander W + defined as

∀t ∈ [0, 1] W + (t) |W (τ + t(1 -τ ))| √ 1 -τ , where τ sup{t ∈ [0, 1] : W (t) = 0}.
We use that the Brownian meander W + can be seen as a Brownian motion conditioned to be positive (Денисов, 1983). This is the main ingredient of Lemma 3. Lemma 3. For any h and η, the expected number of near-optimal points is bounded as

E[N h (η)] ≤ 6η 2 2 h .
This lemma answers a question raised by [START_REF] Munos | Optimistic optimization of deterministic functions without the knowledge of its smoothness[END_REF]: What is the near-optimality dimension of the Brownian motion? We set η h η ε 1/2 h . In dimension one with the pseudo-distance (x, y) = η ε (|y -x|), the near-optimality dimension measures the rate of increase of N h (η h ), the number of η h -near-optimal points in [0, 1] of the form k/2 h . In Lemma 3, we prove that in expectation, this number increases as O η 2 h 2 h = O(log(1/ε)), which is constant with respect to h. This means that for a given ε, there is a metric under which the Brownian is Lipschitz with probability at least 1 -ε and has a near-optimality dimension d = 0 with C = O(log(1/ε)).

The final sample complexity bound is essentially constituted by one O(log(1/ε)) term coming from the standard DOO error for deterministic function optimization and another O(log(1/ε)) term because we need to adapt our pseudo-distance to ε such that the Brownian is -Lipschitz with probability 1 -ε. The product of the two gives the final sample complexity bound of O log 2 (1/ε) .

Proof of Lemma 3. We denote by W, the Brownian motion whose maximum M is first hit at the point defined as t M = inf{t ∈ [0, 1]; W (t) = M } and B + a Brownian meander [START_REF] Durrett | Weak convergence to Brownian meander and Brownian excursion[END_REF]. We also define

B + 0 (t) M -W (tM -t • tM )) √ tM and B + 1 (t) M -W (tM + t(1 -tM )) √ 1 -tM • If L
= denotes the equality in distribution, then Theorem 1 of Денисов (1983) asserts that

B + L = B + 0 L = B + 1 and tM is independent from both B + 0 and B + 1 .
We upper-bound the expected number of η-near-optimal points for any integer h ≥ 0 and any η > 0,

E[N h (η)] = E   2 h k=0 1 W k 2 h > M -η   = 2 h k=0 E 1 W k 2 h > M -η = 2 h k=0 E 1 W k 2 h > M -η ∩ k 2 h ≤ tM ∪ W k 2 h > M -η ∩ k 2 h > tM = 2 h k=0 E 1 B + 0 1- k tM 2 h < η √ tM ∩ k 2 h ≤ tM + 2 h k=0 E 1 B + 1 k/2 h -tM 1 -tM < η √ 1-tM ∩ k 2 h > tM .
If X and Y are independent then for any function f ,

E[f (X, Y )] = E[E[f (X, Y )|Y ]] (law of total expectation) ≤ supy E[f (X, y)|Y = y] (for any Z, E[Z] ≤ supw Z(w)) = supy E[f (X, y)].
(because X and Y are independent)

Since t M is independent from B 0 and B + 1 , then using the above with X = (B + 0 , B + 1 ), Y = t M , and

f : (x0, x1), y → 2 h k=0 1 x0 1- k tM 2 h < η √ y ∩ k 2 h ≤ y +1 x1 k/2 h -tM 1-y < η √ 1-y ∩ k 2 h > y , we can claim that E[N h (η)] = E[f (X, Y )] ≤ sup u∈[0,1] E[f (X, u)] ≤ sup u∈[0,1]    2 h k=0 E 1 B + 0 1 - k u2 h < η √ u ∩ k 2 h ≤ u    + sup u∈[0,1]    2 h k=0 E 1 B + 1 k/2 h -u 1 -u < η √ 1 -u ∩ k 2 h > u    = sup u∈[0,1]    u2 h k=0 P B + 0 1 - k u2 h < η √ u    + sup u∈[0,1]    2 h k= u2 h P B + 1 k/2 h -u 1 -u < η √ 1 -u    = 2 sup u∈[0,1]    u2 h k=0 P B + 0 1 - k u2 h < η √ u    = 2 sup u∈[0,1] {α1 + α2 + α3 + α4}, with α1 = 2 h η 2 k=0 P B + 0 1 - k u2 h < η √ u , α2 = u2 h 2 k= 2 h η 2 P B + 0 1 - k u2 h < η √ u α3 = u2 h -2 h η 2 k= u2 h 2 P B + 0 1 - k u2 h < η √ u , and α4 = u2 h k= u2 h -2 h η 2 P B + 0 1 - k u2 h < η √ u •
Since a probability is always upper-bounded by 1, we bound α 1 and α 4 both by η 2 2 h , to get that α 1 + α 4 ≤ 2η 2 2 h . We now bound the remaining probabilities appearing in the above expression by integrating over the distribution function of Brownian meander [START_REF] Durrett | Weak convergence to Brownian meander and Brownian excursion[END_REF], Equation 1.1),

P B + 0 (t) < x = 2 √ 2π x 0 y exp -y 2 /(2t) t √ 2πt y 0 exp -w 2 /(2(1 -t)) 2π(1 -t) dw dy ≤ 2 t 2πt(1 -t) x 0 y 2 exp - y 2 2t dy ≤ 2x 3 3t 2πt(1 -t) ≤ 2x 3 3t(1 -t) 2πt(1 -t) = 2 3 √ 2π x t(1 -t) 3 ,
where the first two inequalities are obtained by upperbounding the terms exp(•) by one. Now, we use the above bound to bound α 2 + α 3 ,

α2 + α3 = u2 h 2 k= 2 h η 2 P B + 0 1 - k u2 h < η √ u + u2 h -2 h η 2 k= u2 h 2 P B + 0 1 - k u2 h < η √ u ≤ u2 h 2 k= 2 h η 2 2 3 √ 2π   η √ u 1 -k u2 h k u2 h   3 + u2 h -2 h η 2 k= u2 h 2 2 3 √ 2π   η √ u 1 -k u2 h k u2 h   3 ≤ u2 h 2 k= 2 h η 2 1 6 √ π   η √ u k u2 h   3 + u2 h -2 h η 2 k= u2 h 2 1 6 √ π   η √ u 1 -k u2 h   3 ≤ u2 h 2 k= 2 h η 2 1 6 √ π   η √ u k u2 h   3 + u2 h -u2 h 2 k= 2 h η 2 1 6 √ π   η √ u u2 h u2 h + k u2 h  
Changing the indexation as k = -k + u2 h , we get

α2 + α3 ≤ 2 h η 2 3/2 6 √ π    u2 h 2 k= 2 h η 2 1 k + u2 h -u2 h 2 k= 2 h η 2 1 k 3/2    ≤ 2 h η 2 3/2 3 √ π ∞ k= 2 h η 2 1 k 3/2 ≤ 2 h η 2 3/2 3 √ π 3 2 h η 2 ≤ 1 √ π η 2 2 h ≤ η 2 2 h ,
where in the last line we used that for any k

0 ≥ 1, ∞ k=k 0 1 k 3/2 ≤ 1 k 3/2 0 + ∞ k=k 0 +1 k k-1 1 u 3/2 du = 1 k 3/2 0 + ∞ k 0 1 u 3/2 du = 1 k 3/2 0 + 2 √ k0 ≤ 3 √ k0 •
We finally have

∀u, α1 + α2 + α3 + α4 ≤ 3η 2 2 h
and therefore

E[N h (η)] ≤ 6η 2 2 h .
To conclude the analysis, we put Lemma 2 and Lemma 3 together in order to bound the sample complexity conditioned on event C.

Lemma 4. There exists c > 0 such that for all ε

≤ 1/2, E[N ε | C ] ≤ c log 2 (1/ε).
Proof. By definition of h max ,

ε 2 ≤ 5 2 × 2 hmax ln 2 hmax ε ≤ 5 2 2 hmax /ε 2 hmax
, from which we deduce that Furthermore, using Lemma 1, we get that for any ε ≤ 1/2,

E[N h (η h )] = E[N h (η h )| C ]P[C] + E[N h (η h )| C c ]P[C c ] ≥ E[N h (η h )| C ](1 -ε) ≥ E[N h (η h )| C ] 2 •
We now use Lemma 2 and Lemma 3 to get

E[Nε| C ] ≤ 2 hmax h=0 E N h 5 2 ln(2 h+1 /ε) 2 h C ≤ 60 hmax h=0 ln 2 h+1 /ε = 60 hmax+1 h=1 (ln(1/ε) + h ln 2) = O h 2 max + hmax log(1/ε) = O log 2 (1/ε) .
We also bound the sample complexity on C c , the event complementary to C, by the total number of possible intervals [k/2 h , (k + 1)/2 h ] with η ε 1/2 h ≤ ε. Then, we combine it with Lemma 1, Lemma 2, and Lemma 3 to get Proposition 2 which is the second part of the main theorem. Proposition 2. There exists c > 0 such that for all ε < 1/2,

E[Nε] ≤ c log 2 (1/ε).
Proof. From the law of total expectation,

E[Nε] = E[Nε| C ]P[C] + E[Nε|C c ]P[C c ] ≤ E[Nε| C ] + ε 5 E[Nε|C c ].
By Lemma 4, we have that E[N ε | C ] = O log 2 (1/ε) . Now let h max be the maximum h at which points are evaluated. Then,

Nε ≤ hmax h=0 2 h = 2 hmax+1 -1 ≤ 2 hmax+1 .
As in the proof of Lemma 4, we get 2 = O ε -5 . We finally obtain that

E[Nε] ≤ O log 2 (1/ε) + O(1) = O log 2 (1/ε) .
Proposition 2 together with Proposition 1 establish the proof of Theorem 1, the holy grail.

Numerical evaluation of OOB

For an illustration, we ran a simple experiment and for different values of ε, we computed the average empirical sample complexity N ε on 250 independent runs that you can see on the left plot. We also plot one point for each run of OOB instead of averaging the sample complexity, to be seen on the right. The experiment indicates a linear dependence between the sample complexity and log 2 (1/ε).

Conclusion and ideas for extensions and generalizations

We presented OOB, an algorithm inspired by DOO [START_REF] Munos | Optimistic optimization of deterministic functions without the knowledge of its smoothness[END_REF] that efficiently optimizes a Brownian motion. We proved that the sample complexity of OOB is of order O log 2 (1/ε) which improves over the previously best known bound [START_REF] Calvin | Adaptive approximation of the minimum of Brownian motion[END_REF]. As we are not aware of a lower bound for the Brownian motion optimization, we do not know whether O log 2 (1/ε) is the minimax rate of the sample complexity or if there exists an algorithm that can do even better.

What would be needed to do if we wanted to use our approach for Gaussian processes with a different kernel? The optimistic approach we took is quite general and only Lemma 3 would need additional work as this is the ingredient most specific to Brownian motion. Notice, that Lemma 3 bounds the number of near-optimal nodes of a Brownian motion in expectation. To bound the expected number of near-optimal nodes, we use the result of Денисов (1983) which is based on 2 components:

1 A Brownian motion can be rewritten as an affine transformation of a Brownian motion conditioned to be positive, translated by an (independent) time at which the Brownian motion attains its maximum. 2 A Brownian motion conditioned to be positive is a Brownian meander. It requires some additional work to prove that a Brownian motion conditioned to be positive is actually properly defined.

A similar result for another Gaussian process or its generalization of our result to a larger class of Gaussian processes would need to adapt or generalize these two items in Lemma 3. On the other hand, the adaptation or generalization of the other parts of the proof would be straightforward. Moreover, for the second item, the full law of the process conditioned to be positive is actually not needed, only the local time of the Gaussian process conditioned to be positive at points near zero.

  More formally, let ε be the required precision, the only given argument of the algorithm. For any 0 ≤ a < b ≤ 1, the interval [a, b] is associated with an upper bound B [a,b] defined by B [a,b] max(W (a), W (b)) + ηε(b -a), where ∀δ > 0 s.t. εδ ≤ of a set I of intervals [a, b] with W already being sampled at a and b. The algorithm first samples W (1), W (1) ∼ N (0, 1), in order to initialize the set I to the singleton {[0, 1]}. Then, OOB keeps splitting the interval I ∈ I associated with the highest upper bound B I quam necessarium.

  ∀x ≥ max(W (a), W (b)) : P sup t∈[a,b]

Proof.

  Let I = [a, b] be an interval of I t such that max(W (a), W (b)) + η ε (b -a) < M . Let I next ∈ I t be the interval that the algorithm would split after t function evaluations.. Since the algorithm only splits the interval with the highest upper bound, then B Inext = sup I∈It B I . Moreover, if we let I max ∈ I t be one of the intervals where a maximum is reached, t max ∈ arg max t∈[0,1] W (t) M and t max ∈ I max , then on event C, max(W (a), W (b)) + ηε(b -a) BI < M ≤ BI max ≤ BI next . Therefore, under C, a necessary condition for an interval I = [a, b] to be split during the execution of OOB is that max(W (a), W (b)) ≥ M -η ε (b -a), which means that either a or b or both are η ε (b -a)-near-optimal.

  Output: location t ε ← arg max ti W (t i ) and its value W t ε

		3, 4, . . . do
	4:	[a, b] ∈ arg max I∈I B I {break ties arbitrarily}
	5:	if η ε (b -a) ≤ ε then
	6:	break
	7:	end if
	8: 9: 10: end for t i ← W a+b 2 I ← {I ∪ [a, a+b 2 ] ∪ [ a+b 2 , b]}\{[a, b]}
	11:

This holds despite ηε(•) is not always decreasing.

•

Acknowledgements This research was supported by European CHIST-ERA project DELTA, French Ministry of Higher Education and Research, Nord-Pas-de-Calais Regional Council, Inria and Otto-von-Guericke-Universität Magdeburg associated-team north-European project Allocate, French National Research Agency projects ExTra-Learn (n.ANR-14-CE24-0010-01) and BoB (n.ANR-16-CE23-0003), FMJH Program PGMO with the support to this program from Criteo, a doctoral grant of École Normale Supérieure, and Maryse & Michel Grill.


C.-H Pham

C Tor-Díez

H Meunier

N Bednarek

R Fablet

N Passat

F Rousseau

SIMULTANEOUS SUPER-RESOLUTION AND SEGMENTATION USING A GENERATIVE ADVERSARIAL NETWORK: APPLICATION TO NEONATAL BRAIN MRI

Keywords: Super-resolution, segmentation, 3D generative adversarial networks, neonatal brain MRI

The analysis of clinical neonatal brain MRI remains challenging due to low anisotropic resolution of the data. In most pipelines, images are first re-sampled using interpolation or single image super-resolution techniques and then segmented using (semi-)automated approaches. Image reconstruction and segmentation are then performed separately. In this paper, we propose an end-to-end generative adversarial network for simultaneous high-resolution reconstruction and segmentation of brain MRI data. This joint approach is first assessed on the simulated low-resolution images of the high-resolution neonatal dHCP dataset. Then, the learned model is used to enhance and segment real clinical low-resolution images. Results demonstrate the potential of our proposed method with respect to practical medical applications.

INTRODUCTION

Long-term studies of the outcome of prematurely born infants have clearly documented that the majority of such infants may have significant motor, cognitive, and behavioral deficits. However, there is a limited understanding of the nature of the cerebral abnormality underlying these adverse neurologic outcomes. Magnetic Resonance Imaging (MRI) provides unique opportunities for in vivo investigation of the early developing human brain. However, the analysis of clinical neonatal brain MRI data remains challenging mainly due to low anisotropic image resolution. Improving morphological data processing such as image resolution enhancement

The research leading to these results has been supported by the ANR MAIA project, grant ANR-15-CE23-0009 of the French National Research Agency (http://recherche.imt-atlantique.fr/maia); INSERM and Institut Mines Télécom Atlantique (Chaire "Imagerie médicale en thérapie interventionnelle"); the Fondation pour la Recherche Médicale (grant DIC20161236453); and the American Memorial Hospital Foundation. We also gratefully acknowledge the support of NVIDIA Corporation with the donation of the Titan Xp GPU used for this research. and brain segmentation, is a key-point to provide robust morphometry analysis tools to the community.

One of the first key components of the processing pipeline of clinical MRI data is the upsampling image estimation. Super-resolution (SR) is a post-processing technique that aims at enhancing the resolution of an imaging system [START_REF] Greenspan | Super-resolution in medical imaging[END_REF]. SR is a challenging inverse problem; in particular the estimation of texture and details remains difficult. Recently, supervised deep learning-based techniques have shown great improvement over model-based approaches. In this context, applying 3D convolutional neural networks (CNNs) yields promising results for MRI data [START_REF] Pham | Brain MRI super-resolution using deep 3D convolutional networks[END_REF][START_REF] Chen | Brain MRI super resolution using 3D deep densely connected neural networks[END_REF]. However, the use of 2 -norm loss leads to smooth, unrealistic high resolution images [START_REF] Johnson | Perceptual losses for real-time style transfer and super-resolution[END_REF][START_REF] Ledig | Photo-realistic single image superresolution using a generative adversarial network[END_REF]. Generative adversarial networks (GANs) have thus been proposed to estimate textured and sharper images [START_REF] Ledig | Photo-realistic single image superresolution using a generative adversarial network[END_REF][START_REF] Chen | Efficient and accurate MRI super-resolution using a generative adversarial network and 3D multilevel densely connected network[END_REF].

Once the high resolution image reconstruction is performed, the implementation of an automatic segmentation robust approach is crucial for fine brain structure analysis [START_REF] Makropoulos | A review on automatic fetal and neonatal brain MRI segmentation[END_REF]. Segmenting thin structures such as the neonatal cortical gray matter remains difficult and the segmentation step is always considered separately from image reconstruction.

In this article, we propose an end-to-end GAN-based approach which can generate both the perceptually superresolved image and a cortical segmentation map from a single low-resolution (LR) image. The proposed approach called SegSRGAN is both assessed on simulated data and real clinical data.

METHOD

In this work, we focus on joint reconstruction and cortical segmentation of clinical LR neonatal brain MRI data. Section 2.1 describes the mathematical formulation of SR and Section 2.2 focuses on image segmentation. The proposed joint method is presented in Section 2.3, while technical details on the induced network architecture are provided in Section 2.4. Experiments on simulated and clinical MRI data are described and discussed in Section 3.

Formulation of single image super-resolution

The objective of single image SR is to estimate a highresolution (HR) image X ∈ R m from one observed LR image Y ∈ R n . SR problem can be formulated using the following linear observation model:

Y = H ↓ BX + N = ΘX + N ( 1
)
where N is the additive noise, B ∈ R m×m is a blur matrix (depending on the point spread function),

H ↓ : R m → R n is a downsampling decimation and Θ = H ↓ B ∈ R n×m (m > n).
A popular approach that solves SR problem defines the matrix Θ -1 as the combination of a restoration operator F ∈ R m×m and an upscaling interpolation operator

S ↑ : R n → R m computing the interpolated LR image Z ∈ R m (Z = S ↑ Y).
In the context of supervised learning, given a set of HR images X i and their corresponding LR images Y i , the restoration operator F can be estimated by minimizing the following loss function:

F = arg min F i X i -F (Z i ) 2
2 . However, it is known that the use of 2 -norm may lead to oversmoothing high resolution images. In order to provide realistic HR images, perceptual loss function [START_REF] Johnson | Perceptual losses for real-time style transfer and super-resolution[END_REF] have been used in a GAN [START_REF] Ledig | Photo-realistic single image superresolution using a generative adversarial network[END_REF]. This is a paradigm shift since it is no longer a question of minimizing only the reconstruction error but of estimating a realistic image, i.e. a high resolution image that corresponds to the observation model with realistic texture details.

A perceptual loss can be formulated as the weighted sum of the content loss (based, e.g., on pixel-wise mean squared error loss) and an adversarial loss component. In GAN-based approaches, the purpose is to train a generating network G that estimates for a given LR input image a corresponding HR image. The goal of the discriminator network D is to classify real and simulated HR images.

Formulation of image segmentation

In this work, image segmentation is viewed as a supervised regression problem: S X = R (X) where R denotes a non-linear mapping from the upscaled image X to the segmentation map S X . Similarly to the SR problem, assuming that we have a set of images X i and corresponding segmentation maps S Xi , a general approach for solving this segmentation problem is to find the mapping R by minimizing the following loss function: R = arg min R i S Xi -R(X i ) 2 2 . Unlike the SR problem, the use of 2 -norm is less critical as it is expected to estimate smooth segmentation maps.

Joint mapping by generative adversarial networks

We propose the use of a GAN-based approach to estimate jointly a HR image and its corresponding segmentation map from one LR image. To this end, a convolution-based generator network G takes as input an interpolated LR image Z and computes a HR image X and a segmentation map S X by minimizing the following reconstruction loss:

L rec = min G i ρ ((X, S X ) i -G(Z i )) (2)
where (X, S X ) i are concatenated along the feature channel.

In this work, we use a robust loss as Charbonnier loss [START_REF] Charbonnier | Deterministic edge-preserving regularization in computed imaging[END_REF]:

ρ(x) = √ x 2 + ν 2
where ν is set to 10 -3 . The discriminator network D attempts to distinguish the real data (X, S X ) and the generated ones G (Z). The game between the generator G and the discriminator D is usually modeled with a minimax objective. However, using such loss function, GAN may be unstable or can suffer from mode collapse during training. Thus, in this work, we propose to use Wasserstein GAN loss described in [START_REF] Gulrajani | Improved training of Wasserstein GANs[END_REF]:

L adv = min G max D E X∼PX,SX∼P S X [D((X, S X ))]- E Z∼PZ [D(G (Z))] + λ gp E XS [( (∇ XS D( XS) 2 -1) 2 ]
(3) where XS is the interpolation of the true data and the generated one as (1-)(X, S X )+ G (Z), ∼ U[0, 1] (uniform distribution). λ gp and ∇ denote the gradient penalty coefficient and gradient operator, respectively. The images X, S X and Z are extracted randomly from the data distributions of HR images P X , HR segmentation maps P SX and LR images P Z . The terms D((X, S X )), D(G (Z)) and D( XS) are the responses of the discriminator with respect to the real data, the generated data and the interpolated data, respectively. The full objective function is expressed as: L total = L rec + λ adv L adv where λ adv is a trade-off parameter between reconstruction loss and adversarial loss.

Architecture of generator and discriminator networks

The generator network is a convolution-based network with residual blocks. Let C i j -S k be a block consisting of the following layers: a convolution layer of j filters of size i 3 with stride of k, an instance normalization layer (InsNorm) [START_REF] Ulyanov | Improved texture networks: Maximizing quality and diversity in feed-forward stylization and texture synthesis[END_REF] and a rectified linear unit (ReLU). R k denotes a residual block as Conv-InsNorm-ReLU-Conv-InsNorm that contains 3 3 convolution layers with k filters. U k denotes layers as Upsampling-Conv-InsNorm-ReLU layer with k filters of 3 3 and stride of 1. After the last layer, we apply a sigmoid activation for the channel of segmentation map and an element-wise sum of the intensity channel and the interpolated LR image. The generator architecture is:

C 7 16 -S 1 , C 3 32 -S 2 , C 3 64 -S 2 , R 64 , R 64 , R 64 , R 64 , R 64 , R 64 , U 32 , U 16 , C 7
2 -S1 . The discriminator network contains five convolutional layers with an increasing number of filter kernels, increasing by a factor of 2 from 32 to 512 kernels. Let C k be a block consisting of the following layers: a convolution layer of k filters of size 4 3 with stride of 2 and a Leaky ReLU with a negative slope of 0.01. The last layer C 2 1 is a 2 3 convolution filter with stride of 1. No activation layer is used after the last layer. The discriminator consists of C 32 , C 64 , C 128 , C 256 , C 512 , C 2 1 .

EXPERIMENTS AND RESULTS

Datasets and network training

To assess the ability to reconstruct HR volume and segment the cerebral cortex , we applied the proposed method on T2weighted (T2w) MR images of the developing Human Connectome Project 1 (dHCP). 40 T2w images were acquired using a 3T Achieva scanner with a 0.5 × 0.5 × 0.5 mm 3 resolution with TR = 12 000 ms TE = 156 ms, respectively. 30 images were used for training networks, whereas the other 10 were used as testing images. As in [START_REF] Greenspan | Super-resolution in medical imaging[END_REF], LR images were generated by using a Gaussian blur with the full-width-at-halfmaximum (FWHM) set to slice thickness before a downsampling step to obtain a 0.5 × 0.5 × 1.5 mm 3 resolution. We have also applied the proposed method onto clinical neonatal MRI data acquired in the neonatology service of Reims Hospital. These LR images have a resolution of 0.446 × 0.446 × 3 mm 3 . 40 HR images of the dataset dHCP were filtered and downsampled as in [START_REF] Greenspan | Super-resolution in medical imaging[END_REF] in order to generate LR images with a same resolution as clinical data. The network was trained using 40 pairs of simulated data and then applied to real LR images for visual evaluations. All data had bias correction and for network training, they were normalized between 0 and 1. No subjects nor image patches appear twice in the different subsets.

The 3D network was trained over 200 epochs with batch size of 16, using Adam method with learning rate of 0.0001 and updates the discriminator 5 times before training the generator as in [START_REF] Gulrajani | Improved training of Wasserstein GANs[END_REF]. The parameters λ adv and λ gp were set to 0.001 and 10 respectively. The training patch size is 64 3 . At test time, the whole HR image and segmentation volume were reconstructed by the weighted predictions of patches. A thresholding at 0.5 has been performed to obtain binary segmentation maps.

Results

Peak signal-to-noise ratio (PSNR) and structural similarity index (SSIM) have been used to evaluate the performance of SR reconstructions. Table 1 provides a summary of quantitative evaluations for the following methods: cubic spline interpolation, a 20-layers CNN-based SR approach (20L-SRReCNN) [START_REF] Pham | Multi-scale brain MRI super-resolution using deep 3D convolutional networks[END_REF] and our proposed SegSRGAN. It can be seen that 20L-SRReCNN provides highest PSNRs as in [START_REF] Johnson | Perceptual losses for real-time style transfer and super-resolution[END_REF][START_REF] Ledig | Photo-realistic single image superresolution using a generative adversarial network[END_REF] since this approach minimizes a 2 -norm-based loss. However, while the two CNN-based approaches (20L-SRReCNN and SegSRGAN) lead qualitatively to similar realistic results on dHCP dataset (see Fig. 1), the proposed approach provides the best reconstructed HR image on clinical data with better contrast on cortical gray matter (see Fig. 2). The Dice index is used to evaluate the cortical segmentation maps obtained by the following state-of-the-art methods: iterative multi-atlas patch-based approach (IMAPA) [START_REF] Tor-Díez | An iterative multi-atlas patch-based approach for cortex segmentation from neonatal MRI[END_REF], DrawEM [START_REF] Makropoulos | Automatic whole brain MRI segmentation of the developing neonatal brain[END_REF] and the proposed SegSRGAN. As in a typical clinical setting, the three methods have been applied on interpolated images. Table 2 shows that quantitatively the proposed approach lead to the best cortical segmentation results with significant improvement with respect to the two other methods. Moreover, as mentioned in [START_REF] Tor-Díez | An iterative multi-atlas patch-based approach for cortex segmentation from neonatal MRI[END_REF], the use of IMAPA applied on original HR dHCP images leads to a mean DICE of 0.887 (std. dev. of 0.011) that is very similar to the results obtained with SegSRGAN (applied on interpolated images). Results on real LR data (see Fig. 2) confirm the potential of the proposed approach for fine analysis of clinical neonatal brain MRI.

CONCLUSION

In this paper, we have presented a simultaneous superresolution and segmentation method for 3D brain MR images using a generative adversarial network. Our experiments on both simulated and clinical data have shown that better performance can be achieved by this joint approach compared to state-of-the-art techniques, opening up new perspectives in the processing of clinical LR neonatal brain MRI data.

Fig. 1 .

 1 Fig. 1. SR and segmentation results for one dHCP subject: (a) original HR image; (b-d) SR reconstruction of the LR image generated from (a); (e) segmentation ground-truth of (a); (f,g) segmentation of (b); (h) HR segmentation from the LR image using the joint SegSR-GAN method.

Fig. 2 .

 2 Fig. 2. Reconstruction (b-d) and segmentation results (e) on a real LR neonatal brain image (a) with voxel size of 0.446 × 0.446 × 3 mm 3 , re-sampled to 0.5 × 0.5 × 0.5 mm 3 .

Table 1 .

 1 Quantitative evaluation of SR methods on dHCP dataset.

		Interpolation 20L-SRReCNN [11] SegSRGAN
	PSNR	30.70	35.84	31.75
	SSIM	0.9492	0.9739	0.9624

Table 2 .

 2 Quantitative evaluation of segmentation methods on dHCP dataset.

		IMAPA [12] DrawEM [13] SegSRGAN
	Dice (mean)	0.788	0.818	0.886
	Dice (std. dev.)	0.061	0.014	0.011

http://www.developingconnectome.org