P Rolet

O Teytaud

Bandit-based Estimation of Distribution Algorithms for Noisy Optimization: Rigorous Runtime Analysis

) that reaches the bound within logarithmic factors. We emphasize the differences with empirical derived published algorithms.

Introduction

Optimization in noisy environments deals with fitness functions for which multiple evaluations do not necessarily give the same result, due to noise in the computation. Noise occurs in many real world applications. In numerous cases, it is due to finite precision of Monte-Carlo sampling, and often the noise is centered to zero, i.e. the goal is to have an optimal expected value.

The simplest solution for reducing noise consists in evaluating several times the same point; by averaging, the variance decreases. Unfortunately, depending on detailed experimental conditions, various points of views can be found in the literature, notably regarding evolutionary strategies. [START_REF] Fitzpatrick | Genetic algorithms in noisy environments[END_REF] concludes in the case of genetic algorithms that averaging does not perform well: when considering the tradeoff between the computational overhead of averaging and the convergence rate, convergences are always slower with averaging. In the case of evolution strategies (see [START_REF] Beyer | The Theory of Evolutions Strategies[END_REF]), [START_REF] Hammel | Evolution strategies on noisy functions: How to improve convergence properties[END_REF] and [START_REF] Arnold | Efficiency and mutation strength adaptation of the (mu/mui,lambda)-es in a noisy environment[END_REF] draw different conclusions: [START_REF] Hammel | Evolution strategies on noisy functions: How to improve convergence properties[END_REF] concludes that strong averaging is required, whereas [START_REF] Arnold | Efficiency and mutation strength adaptation of the (mu/mui,lambda)-es in a noisy environment[END_REF] concludes that increasing the population size is better than averaging. According to [START_REF] Arnold | Efficiency and mutation strength adaptation of the (mu/mui,lambda)-es in a noisy environment[END_REF], results in [START_REF] Hammel | Evolution strategies on noisy functions: How to improve convergence properties[END_REF] are due to a poor mutation strength adaptation schema; however, interestingly, they point out that for various noise models, each usual mutation-strength adaptation schema can lead to poor results: [START_REF] Arnold | Efficiency and mutation strength adaptation of the (mu/mui,lambda)-es in a noisy environment[END_REF] compares mutative self-adaptation [START_REF] Rechenberg | Evolutionstrategie: Optimierung Technischer Systeme nach Prinzipien des Biologischen Evolution[END_REF][START_REF] Schwefel | Numerical Optimization of Computer Models[END_REF] (two variants, one with arithmetic averaging of mutation strength and the other with geometric averaging) and cumulative selfadaptation [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF], and concludes that in both cases, there are simple examples of noise models for which poor behaviors can occur. It might then be that averaging is required to avoid such behaviours: this problem remains open. An improvement of cumulative step-length adaptation is possible with increased population size, but only to a limited extent [START_REF] Arnold | Evolution strategies with cumulative step length adaptation on the noisy parabolic ridge[END_REF]. An analysis using Markov chains has been adapted from the noise-free case [START_REF] Bienvenue | Global convergence for evolution strategies in spherical problems: some simple proofs and difficulties[END_REF][START_REF] Auger | Convergence results for (1,λ)-SA-ES using the theory of ϕ-irreducible markov chains[END_REF] to the noisy case in [START_REF] Jebalia | On multiplicative noise models for stochastic search[END_REF][START_REF] Teytaud | On the adaptation of the noise level for stochastic optimization[END_REF]; reevaluations are suggested, but it is pointed out that in many cases, this would not be sufficient. [START_REF] Heidrich-Meisner | Hoeffding and bernstein races for selecting policies in evolutionary direct policy search[END_REF] proposed the use of "bandits" [START_REF] Lai | Asymptotically efficient adaptive allocation rules[END_REF][START_REF] Auer | Finite time analysis of the multiarmed bandit problem[END_REF] and Bernstein's races [START_REF] Mnih | Empirical Bernstein stopping[END_REF] for choosing the number of function evaluations spent on a given individual (bandits and races will be presented in section 3). This idea, albeit promising, could not be applied as such to the problem; the authors had to add several tricks in the implementation in order to get acceptable results. Furthermore, they could not get convergence rates, convergence proofs, or convincing experimental curves exposing a good behavior. The main goal of this paper is to get convergence proofs and convergence rates using a modified bandit algorithm; the difference between our approach and the approach in [START_REF] Heidrich-Meisner | Hoeffding and bernstein races for selecting policies in evolutionary direct policy search[END_REF] will be cleary emphasized in section 4.

An peculiar case of the noisy optimization setting is when the variance of the noise decreases to zero near the optimum. This has been tackled in [START_REF] Jebalia | On multiplicative noise models for stochastic search[END_REF], when the noise satisfies conditions such as, typically,

f noise (x) = f (x) × (1 + N )
where f (x) is the sphere function f (x) = ||x|| 2 and N is an independent noise variable. Assuming that N has a density in a bounded range [m, M ] with -1 < m < M , [START_REF] Jebalia | On multiplicative noise models for stochastic search[END_REF] shows that the scale-invariant 1 + 1-ES converges. However, there is no information on the convergence rate, and the theoretical algorithm cannot be used in practice (see note 5 in [START_REF] Jebalia | On multiplicative noise models for stochastic search[END_REF]) as it uses some prior information, at each iteration, on the distance to the optimum. In this paper:

-we derive a lower bound on the number of requests necessary for reaching a given precision for a given confidence level; -we consider a different model of noise, including cases in which the noisy measurement can be arbitrarily small with respect to the expected fitness, i.e. f noise (x)/f (x) is not lower bounded; -we propose an explicit algorithm: it does not require information on the position of the optimum, or on the distance towards the optimum; -we prove a tight convergence rate for any fixed dimension.

In all the paper, #E denotes the cardinal of the set E, (x) i denotes the i th coordinate of vector x, and [[a, b]] denotes {a, a + 1, a + 2, . . . , b}.

The paper is organized as follows. Section 2 will show a lower bound on the number of requests necessary for reaching a given precision with a given probability. Section 3 will propose an algorithm for noisy optimization, along with a proof that this algorithm matches the lower bound. This shows the tightness of the lower bound given in section 2, within logarithmic factors. Section 4 will then discuss the difference with existing algorithms and emphasize the importance of subtle elements in the generation of the offspring.

Lower bounds

In this section, we show a lower bound on the number of requests before an optimization algorithm can find the optimum with precision ǫ and confidence 1 -δ. This means that it is not possible to design an algorithm which finds, with a smaller number of requests to the fitness function, the optimum of the fitness function with precision ǫ and with probability 1 -δ.

The optimization framework is described in Algorithm 1. This is a black-box optimization framework: the algorithm can request the fitness value at any point of the domain, and this costs one request. No other information on the fitness function is available.

We consider a fitness function f parameterized by the (unknown) location of its optimum, t. The goal is to find the optimum t of f (., t), by observing noisy measurement of f at x i : 1 with probability 1 -f (x i , t) (i.e. if θ i > f (x i , t)) and 0 otherwise. θ is the random variable accounting for the noise. This fits applications based on highly noisy optimization, such as games: let x be a parameter of a game strategy, that we wish to set at its best value; one noisy observation is a game against a baseline, resulting either in a win or in a loss; the aim is to find the value of x maximizing the probability of winning.

In the following, t is not handled stochastically, i.e. the lower bounds are not computed in expectation w.r.t all the possible fitness functions yielded by different values of t. Rather, we will consider the worst case on t. Therefore the only random variable in this framework is θ, which decides the answers (see Algorithm 1). As a consequence, all probability / expectation operators are w.r.t θ.

Algorithm 1 Noisy optimization framework. Opt is an optimization algorithm; it takes as input a sequence of visited points and their measured fitness values, and outputs a new point to be visited. Fitness values are binary: we are looking for points x of the domain such that f (x, t) is as small as possible. The algorithm Opt is successful on target function f parameterized by t and random noise θ if Reward(t, θ, Opt) is small; details in the text.

Parameters: N , number of fitness evaluations; t, unknown element of X. θ: random state of the nature

∈ [0, 1] N ; each coordinate θi for i ∈ {1, 2, . . . } is uniformly distributed in [0, 1]. for n ∈ [[0, N -1]] do x t,θ n+1 = Opt(x t,θ 1 , . . . , x t,θ n , y t,θ 1 , . . . , y t,θ n ) if n < N -1 then if f (x t,θ n+1 , t) < θn+1 then y t,θ n+1 = 1 else y t,θ n+1 = 0 end if end if end for Loss(t, θ, Opt) = d(t, x t,θ N )
We consider a domain X, a function f : X × X → R, and we define

d(t, t ′ ) = sup x∈X |f (x, t) -f (x, t ′ )|
for t and t ′ in X. In all the paper, B(n, p) is a binomial random variable (sum of n independent Bernoulli variables of parameter p).

The lower-bound theorem can then be stated as follows:

Theorem 1. For any optimization algorithm Opt, let:

-N ∈ N * (a number of points visited); -ǫ 0 > 0 , and 0 < ǫ < ǫ 0 (ǫ is a level of precision of the learned optimum); -D ∈ N * (characterizing a generalized dimension of X, see remark below);

-δ ∈]0, 1[ (desired confidence level).
Assuming that:

-H(ǫ 0 , D):

∀ǫ 1 < ǫ 0 ∃(t 1 , . . . , t D ) ∈ X D , ∀(i, j) ∈ [[1, D]] 2 , i = j ⇒ d(t i , t j ) = ǫ 1 (generalized dimension) -H P AC (ǫ, N, δ): ∀t, P (d(x t,θ N , t) < ǫ/2) ≥ 1 -δ. Then, if δ < 1/2D, P (B(N, ǫ) ≥ ⌈log 2 (D)⌉) ≥ 1 -Dδ. (1)
This theorem shows that if, with probability 1 -δ, in N iterations, we can find the solution with precision ǫ, then Eq. 1 holds; as discussed below, this is a lower bound on the number of iterations N such that precision ǫ is reached.

Remarks:

-Eq. 1 implies a clearer expression of the lower bound (using Chebychev inequality):

N = Ω(log 2 (D)/ǫ)
where N is the number of iterations required for reaching precision ǫ with confidence 1 -δ for δ < 1/2D. -As stated above, H(ǫ 0 , D) is analogous to a dimensionality assumption, and more general. In dimension D and for simple examples of function f (see examples below), there always exists Θ(D) points with the same distance to each other. -H P AC (ǫ, N, δ) states that the algorithm reaches precision ǫ with probability 1 -δ after N visited points.

Examples.

-A first natural example is X = [0, 1] D0 , ǫ 0 = 1/10, and f (x, t) = ||x -t||.
In this simple cases, the theorem applies with D = D 0 + 1, and the distance d boils down to the considered norm (i.e. d(t, t ′ ) = ||t -t ′ ||): reaching a precision ǫ = ||x θ N -t|| with confidence 1 -δ for any δ < 1/(2D) will require at least Ω(log 2 (D)/ǫ) points.

-The fitness at the optimum needs not have a variance of 0 (consider for instance f (x, t) = min(1,1 3 + ||x -t||)). 1 -Interestingly, the theorem has the exact same interpretation in the less trivial case of f (x, t) = min(1, ||x -t|| p ), p > 0. Indeed, for ǫ 0 sufficiently small, H(ǫ 0 , D) holds for D = D 0 + 1 and applying the theorem yields N = Ω(log 2 (D)/ǫ). It means that, for any fixed p, the distance between the N th iterate (distance for d(., .)) and the optimum that can be guaranteed with confidence 1 -δ for δ < 1/2D, is θ(1/N ). Since in this case d(x, y) = Θ(||x -y||), a somewhat surprising conclusion can be drawn: the convergence is lower bounded by Θ(1/N ) for any value of p used in the family of fitness functions (the constant in Θ might change, but the rate, N -1 , does not improve).

Besides, for any p > 0; there is no Lipschitz condition on the considered family of functions, and in all those cases a lower bound on N can be derived such that H P AC (ǫ, N, δ) holds.

Proof of Theorem 1:

Consider a fixed ǫ < ǫ 0 , D > 0 and δ < 1/2D. We assume H(ǫ 0 , D) and H P AC (ǫ, N, δ). We consider t 1 , . . . , t D as proposed by hypothesis H(ǫ 0 , D):

∀i = j, d(t i , t j ) = ǫ. (2)
We now show that Eq. 1 holds.

To simplify the proof, let us first introduce some relevant definition. First, let us define

C n i (θ) = {j ∈ [[1, D]]; (y tj ,θ 1 , . . . , y tj ,θ n ) = (y ti,θ 1 , . . . , y ti,θ n )}.
C n i (θ) is the set of t in t 1 , . . . , t D such that all answers given by the oracle are the same as on t i . It will be useful to consider the range of values of f at the t i 's such that C n i is large. For that, we define k n (θ) as

k n (θ) = min{i ∈ [[1, D]]; #C n i (θ) is maximal}.
The important point here is that k n (θ) is an index i such that C n i (θ) is of maximal cardinal; the min is only here in case of equality. So, k n (θ) is the index of some t i such that "many" (as many as possible) t j 's raise the same labels y tj ,θ 1 , . . . , y tj ,θ n . Then we define:

s min,n (θ) = inf i∈C n k (θ) f (t i , x ti,θ n ). s max,n (θ) = max i∈C n k (θ) f (t i , x ti,θ n ).
[s min,n (θ), s max,n (θ)] has the important property that it is "small"; precisely,

∀(n, θ), s min,n (θ) ≥ s max,n (θ) -ǫ (3)
thanks to assumption H(ǫ 0 , D) and thanks to ǫ < ǫ 0 , D > 0 and δ < 1/2D. The proof of Eq. 1 is now as follows:

1. Define the property Error(t, θ) : d(x t,θ N , t) ≥ ǫ/2. Define also G, the "good" set defined by

θ ∈ G iff ∀i, #C N i (θ) = 1. (4)
Then,

θ ∈ G ⇒ ∃i = j, x ti,θ N +1 = x tj ,θ N +1
⇒ Error(t i , θ) or Error(t j , θ) thanks to Eq. 2.

(5)

2. Also,

H P AC (ǫ, N, δ) ⇒ P (G) ≥ 1 -Dδ. (6)
Proof of Eq. 6:

We assume H P AC (ǫ, N, δ), and show that P (G) ≥ 1 -Dδ.

If θ ∈ G, then ∃i, Error(t i , θ) (by Eq. 5 above). Therefore,

P (¬G) ≤ i∈[[1,D]] P (Error(t i , θ)) ≤ Dδ by H P AC (ǫ, N, δ).
This concludes the proof of Eq. 6. 3. Equation 6and the definition of G (Eq. 4) immediately lead to (under assumption H P AC (ǫ, N, δ)):

P (∀i, #C N i (θ) = 1) ≥ 1 -Dδ. (7)
4. We have seen with Eq. 7 that if H P AC is verified, then C N i (θ) must be small for all i and with high probability on θ. We will now, in order to complete the proof, show that C N i (θ) can't be small if N is not large. More precisely, we now show:

P (∀i, #C N i (θ) ≤ 1) ≤ P (B(N, ǫ) ≥ ⌈log 2 (D)⌉). (8)
Proof of Eq. 8 by induction:

-Initialization of the induction:

#C 0 k0(θ) (θ) = D; -With k = k n (θ) for short: • P (#C n k (θ) = #C n-1 k (θ)) ≥ 1 -ǫ (by Eq.
3), and

• if #C n k (θ) = #C n-1 k (θ)) then C n k (θ) ≥ ⌈C n-1 k (θ)/2⌉ (the set C n k θ is at most divided by 2).
This implies by immediate induction that P (#C N i (θ) ≤ 1) ≤ P (B(N, ǫ) ≥ ⌈log 2 (D)⌉). This concludes the proof of Eq. 8. 5. Equations 7 and 8 together conclude the proof of Eq. 1. This concludes the proof of Theorem 1.

3 Tightness: an algorithm matching the lower bound

We have shown in the previous section a lower bound on the number of fitness evaluations for the noisy case. We now show that a bandit-based algorithm matches the bound. We consider optimization in the domain

X = [x - 0 , x + 0 ] = [(x - 0 ) 1 , (x + 0 ) 1 ] × [(x - 0 ) 2 , (x + 0 ) 2 ] × • • • × [(x - 0 ) D , (x + 0 ) D ]; x - 0 and x + 0 are elements of R X and therefore X ⊂ R D . We assume that ||x - 0 - x + 0 || ≤ 1.
Sketch of Algorithm 2. We will use Algorithm 2 for showing the tightness of Theorem 1 within logarithmic factors. It proceeds by iteratively splitting the domain in two (not necessarily equal) halves, and retaining the one that most probably contains the optimum. At iteration n, from the nth domain

[x - n , x + n ], the (n + 1)th domain [x - n+1 , x + n+1 ] is obtained by: -Finding the coordinate c such that δ max n = (x + n ) c -(x - n ) c
is maximal; -Selecting 3 regularly spaced points along this coordinate; -Repeatingly assessing those 3 points until we have confidence that the optimum is closer to one point x ). This idea of using bad points as well as good points has already been used in [START_REF] Arnold | Cumulative step length adaptation for evolution strategies using negative recombination weights[END_REF], and in the optimisation heuristic Breda [START_REF] Gelly | Comparison-based algorithms are robust and randomized algorithms are anytime[END_REF] (similarly to Algorithm 2, but in the deterministic case).

The remaining part of this section will be devoted to the analysis of the runtime of our algorithm. There are two points to prove regarding Algorithm 2:

-Algorithm 2 finds the optimum with precision ǫ and confidence 1 -δ; -This is done with computational cost O( 1ǫ log( 1 ǫ )). The following upper-bounding theorem summarizes the results of this section. 2 Then, Algorithm 2 ensures the following properties:

Theorem 2 (Bandit-based noisy optimization). Consider a domain

[x - 0 - x + 0 ] ⊂ R D with ||x - 0 -x + 0 || ≤ 1. Consider a noisy fitness function f noise with values in [0, 1] such that Ef noise (x) = ||x -x * || and x * ∈ [x - 0 , x + 0 ].
-Soundness: The n th domain has diameter ≤ √ D( 3 4 ) ⌊n/D⌋ ; -Consistency: With probability at least 1 -δ, for all n > 0 the optimum is in the n th domain [x n -, x n + ]; -Rate: For a fixed D, the number of function evaluations before the domain

has diameter ≤ ǫ is O( 1 ǫ log(1/ǫ)).
Remark: the multi-armed bandit framework. This theorem is termed "banditbased" because it is inspired from the multi-armed bandit framework (see e.g. [START_REF] Auer | Finite time analysis of the multiarmed bandit problem[END_REF]). The multi-armed bandit framework is named after the famous "one-armed bandit" machines. There are multiple arms that can be tried by a player, each giving a stochastic reward with unknown mean and variance. The goal is to play the machine that yields the best mean reward most of the time.

In this paper, points of the domain X are similar to arms: they yield a random "reward" (between 0 and 1) and the goal is to find the one that has the best mean (the optimum). This is the spirit of the Berstein Race algorithm, Algorithm 3.

Remark: Bernstein races. Hoeffding's bound [START_REF] Hoeffding | Probability inequalities for sums of bounded random variables[END_REF] is a classical concentration inequality: it shows that the difference between the averages and the expectations in "small", if the sample size is enough. It can be applied to so-called Hoeffding's races [START_REF] Maron | Hoeffding races: Accelerating model selection search for classification and function approximation[END_REF][START_REF] Maron | The racing algorithm: Model selection for lazy learners[END_REF]: given a sample of candidates, one might evaluate all of them, get the stochastic rewards, and keep only for further evaluations the candidates for which the Hoeffding's bound is not too bad. After several iterations, only good candidates are kept.

Bernstein's bound [START_REF] Chernoff | A measure of asymptotic efficiency for tests of a hypothesis based on the sum of observations[END_REF][START_REF] Bernstein | On a modification of chebyshev's inequality and of the error formula of laplace[END_REF][START_REF] Bernstein | The Theory of Probabilities[END_REF] is an improvement of Hoeffding's bound which takes into account the variance of candidates, and not only their range. Their are variants of Bernstein's bounds which can even take into account the empirical variance, and not only the variance itself; these bounds can therefore directly be used instead of Hoeffding's bound in races-this is Bernstein races. The results about Bernstein races used in this paper can be found in [START_REF] Mnih | Empirical Bernstein stopping[END_REF] and [11, p124].

We first give an intuitive idea of our proof, and then the complete proof of Theorem 2.Proof sketch. The algorithm geometrically decreases the domain, hence the soundness property. Thanks to the specific family of functions (the sphere), the optimum is closer to the "good" arm than to the "bad arm" with high probability: this provides the consistency. The number of requests to the fitness functions will then be upper bounded by classical Bernstein bounds; this is only possible thanks to the specific mutation operator in Algorithm 2: it ensures that arms are all "sufficiently different".

Proof of Theorem 2:

In the following,Ef noise (x) will be referred to as f (x), and f (x) will stand for an empirical estimate of f (x) by averaging the rewards given by f noise (x).

First of all, notice that

f noise (x) ∈ [0, 1] implies V ar f noise (x) ≤ E[f noise (x)] = ||x -x * ||.
A (quick and easy) proof of this is provided in appendix.

A step towards proving theorem 2 is stated by the following lemma: Algorithm 2 Algorithm for optimizing noisy fitness functions. Bernstein denotes a Bernstein race, as defined in Algorithm 3. The algorithm is illustrated in Fig. 1.

n ← 0 K0 = √ 1 + D - q 1 4 + D while True do c = arg maxi(x + n )i -(x - n )i //
= (x + n )c -(x - n )c for i ∈ [[1, 3]] do x ′i n = 1 2 (x - n + x + n ). // we consider the middle point (x ′i n )c = (x - n )c + i-1 2 (x + n -x - n )c. //
= Bernstein(x ′1 n , x ′2 n , x ′3 n , 6δ π 2 (n+1) 2 , ǫ1, ǫ2). Let H be the halfspace {x ∈ R D ; ||x -good|| ≤ ||x -bad||. Split the domain: [x - n+1 , x + n+1 ] = H ∩ [x - n , x + n ]. n ← n + 1 end while
Lemma 1 (Convergence of Alg. 2). If Algorithm 2 reaches the n th 0 iteration and x * is in the n th 0 domain, then for all n ≤ n 0 , δ n = x + n -x - n verifies: .

δ c n ≤ δ c n-1 × 3 4 where c = arg max i (x + n ) i -(x - n ) i ; δ max n+D ≤ δ max n × 3
Remark: Since it hasn't yet been proved that Algorithm 3 terminates with x * in the remaining domain, the result is stated under those conditions. It will be shown in propositions 1 and 2 that those conditions are satisfied with high probability.

Proof of Lemma 1:

This is an immediate property of Algorithm 2: at each iteration, the part of the domain which is on the wrong side of H (meaning that x * is in the other part,see Fig. 2) is removed, and its size is at least δ max n /4.

We will now prove that provided Algorithm 2 reaches iteration n, the assumptions made in Algorithm 3 guaranteeing its proper behavior are met. At each iteration, a main axis is selected (the one on which the domain has maximum range). Three equally spaced points are generated in the domain on this axis (this is the offspring). Then, a Bernstein race is applied for choosing a "good" and a "bad" arm among these points. The domain is reduced thanks to this knowledge, removing one fourth or one half of the domain (depending on the position of the good arm and of the bad arm -the best case is when the good and the bad arm are diametrically opposed: see Fig. 2).

Lemma 2 (The conditions of the Bernstein race are met). Assume that

x * ∈ [x - n , x + n ]. Then max (i,j)∈[[1,3]] 2 f (x ′ j n ) -f (x ′ i n ) ≥ 2K 0 δ n max where K 0 = 1 4 √ 1 + D - 1 4 + D ; ( 10
) max i∈[[1,3]] f (x ′ i n ) ≤ √ Dδ n max . ( 11
)
Proof of Lemma 2: Eq. 11 is immediate by definition of f . Let us show why eq. 10 holds. Let x * n be the projection of x * on the line on which lie x

′ 1 n , x ′ 2 n , x ′ 3
n . The result will now be proved for 2). First of all, we have

( x * n ) c ∈ [(x ′ 1 n ) c , (x ′ 2 n ) c ]. The proof for the case ( x * n ) c ∈ [(x ′ 2 n ) c , (x ′ 3 n ) c ] is symmetric (see Figure
∆ n . = max i,j∈[[1,3]] 2 f (x ′ i n ) -f (x ′ j n ) ≥ f (x ′ 3 n ) -f (x ′ 2 n )
Algorithm 3 Bernstein comparison between 3 arms. Eq. 9 is Bernstein's inequality for estimating the risk, for arms with variance ≤ ǫ 1 and precision ǫ (see e.g. [11, p124]).

Bernstein(a1, a2, a3, δ ′ , ǫ1, ǫ2) //We assume that arm a1 (resp. a2, a3) has (expected) fitness p1 (resp. p2, p3). //We assume that p1, p2, p3 ≥ 0. //We assume that ǫ1 ≥ max(p1, p2, p3). //We assume that ǫ2 ≤ maxi pi -mini pii. t = 0, ǫ = ǫ2/2 repeat t ← t + 1 Evaluate the fitness of arms a1, a2, a3 once. Evaluate the precision:

ǫ (t) = sup{ǫ > 0; 3 exp " - tǫ 2 2ǫ1 + 2ǫ/3 « ≤ 6δ ′ π 2 t 2 }. ( 9
)
until There are two arms good and bad with empirical difference ≥ 2ǫ: f (good) -

f (bad) ≥ 2ǫ By Pythagora's theorem, ∀i ∈ [[1, 3]], ||x ′ i -x * || 2 = ||x ′ i n -x * n || 2 -|| x * n -x * || 2 . Thus, ∆ n ≥ ||x ′ 3 n -x * n || 2 -|| x * n -x * || 2 -||x ′ 2 n -x * n || 2 -|| x * n -x * || 2 Note that ||x ′ 3 -x * n || = ||x ′ 2 n -x * n || + δ max n /2 . Define d = || x * n -x * || 2 and a = ||x ′ 3 -x * n ||. Then, observing that δ max n ≥ a ≥ δ max n /2, we have ∆ n ≥ a 2 + d -(a -δ max n /2) 2 + d ≥ a   1 + d/a 2 - 1 - δ max n 2a 2 + d a 2   ≥ δ max n 2 1 + d/a 2 - 1 4 + d a 2 (12)
By setting u = d/a 2 , it is clear that ∆ n is greater than the minimum of u →

√ 1 + u -1/4 + u on the interval [0, D] (since √ d = || x * n -x * || ≤ √ Dδ max n /2
). It remains to find the minimum of that function. Its derivative is

1 2 √ 1 + u - 1 2 1/4 + u
which is always negative, since the left-hand side denominator is always bigger than the right-hand side denominator. Therefore, the function is decreasing, and is minimized for u = √ D. Finally, setting K 0 = 1/4( √ 1 + D -1/4 + D) and injecting in eq. 12 yields ∆ n ≥ 2K 0 δ max n , as stated by Eq. 10. n , x ′2 n , x ′3 n ; the left arm is the "bad" arm, whereas the arm in the center is the "good" arm, i.e. the one which proved to be closer to the optimum than the left arm, with confidence 1 -6δ/(π 2 n 2 ).

"good" arm "bad" arm

Discarded part of the domain

Lemma 1 shows that the domain gets smaller and smaller at a given fixed rate. Lemma 2 shows that, if the optimum is in [x - n , x + n ], then the conditions of Algorithm 3, as called in Algorithm 2, are met, i.e. all arms have a probability ≤ ǫ 1 and there are two arms with parameters p i and p j such that p i > p j + ǫ 2 .

It implies the following proposition:

Proposition 1 (Consistency of Algorithm 3). In the previously defined setting in which algorithms 2 and 3 operate, with probability at least 1 -δ, for all calls to Algorithm 3:

∀n, in the call to Alg. 3 in iteration n of Alg. 2, ∀t, ∀i ∈ {1, 2, 3},

|p i -p i | ≤ ǫ (t) (13)
and therefore

∀n, x * ∈ [x - n , x + n ]. (14)
and

δ max n ≤ 4 5 ⌊n/D⌋ (15)
Proof of Proposition 1:

The fact that δ max n ≤ 4 5 ⌊n/D⌋ is already proved in Lemma 1. We have to show Eq. 13 and Eq. 14.

This can be shown by induction on n.

(i) The case n = 0 holds. This is because in our setting, it is stated that x * ∈ [x - 0 , x + 0 ] and δ 0 max < 1; this shows Eq. 14 for n = 0. By Lemma 2, we can apply Bernstein's race (see e.g. [11, p124] or [START_REF] Mnih | Empirical Bernstein stopping[END_REF]), which state Eq. 13. Equations 13 and 14 are therefore proved for n = 0.

(ii) Let's assume Equations 13 and 14 for a given n, and let's prove them for n + 1. Eq. 14 and the halting condition of Alg. 3 ensure that Algorithm 3 terminates and is right (and, by definition of the Bernstein race, it will terminate and be right with probability at least 1 -δ ′ ), its output is 2 arms (good, bad) such that f (good) < f (bad).

Thus, the split that occurs in Algorithm 2 after the call to the Berntein Race ensures with probability 1-δ ′ that x * remains in the (n+1)th domain; this is Eq. 14. Besides, Lemma 1 guarantees that the domain's edge of length δ max n decreases as wished (Eq. 15). We have now to show Eq. 13. From Eq. 14, Lemma 2 states that the conditions of the Bernstein race are met; This implies Eq. 13 by usual properties of Bernstein races.

(iii) Each induction step n holds with probability 1 -δ ′ n if the preceeding one holds, so the proposition is true for all n with probability at least

∞ n=0 (1 -δ ′ n ) = ∞ n=0 (1 - 6δ π 2 (n + 1) 2 ) = δ
which concludes the proof.

It remains to prove that, at each iteration, Algorithm 3 requires a finite number of requests to the fitness function small enough so that in the end the lower bound is reached: Proposition 2 (Number of iterations of the Bernstein race). Consider Algorithm 3 called by Algorithm 2, for a given iteration n with

x * ∈ [x - n , x + n ] (16)
Assume that ∀t, ∀i ∈ {1, 2, 3},

|p i -p i | ≤ ǫ (t) (17)
Then, with probability 1 -δ ′ , the number of iterations t before the halting condition in Algorithm 3 is O(-log(ǫ 1 )/ǫ 1 ).

Remarks. We have assumed Eq. 17 and Eq. 16. These two assumptions are verified, for all n, with probability at least 1 -δ, as shown in Prop. 1; we can thereforeapply Proposition 1 and Proposition 2 simultaneously, and for all n, with probability 1 -δ.

Proof of Proposition 2: Let's bound the number of fitness evaluations of a call to the Bernstein race. We consider a fixed D, and we will use the following facts:

-Algorithm 3 is called with ǫ 2 = Θ(ǫ 1 ); -The iteration n of Algorithm 2 verifies n = Θ(log(ǫ 1 )); -The Bernstein race 3 stops as soon as ∃i, j; ǫ (t) ≤ 1 4 (p i -p j ), and therefore necessarily the last iteration t n of Algorithm 3 when called at iteration n of algorithm 2 is such that

ǫ (tn-1) > ǫ 2 /4 (18)
Plugging Eq. 18 in Eq. 9 and using ǫ 2 = Θ(ǫ 1 ), n = Θ(log(ǫ 1 )),

δ ′ = θ(1/n 2 ), yields exp (-t n Θ(ǫ 1 )) ≥ Θ(1/n 2 t 2 ). (19)
Eq. 19 yields

t n = O(-log(ǫ 1 )/ǫ 1 ). ( 20
)
This concludes the proof of Proposition 2.

We have, with Eq. 20, a bound on the number of fitness evaluations in terms of ǫ 1 . Consider a precision ǫ.

L(ǫ) = O(-log(ǫ)/ǫ).
We now have to bound the complexity of the whole run of Algorithm 2, and not only the last call to the Bernstein race; we'll see that the cost is indeed essentially in the last iteration. The complexity of the last D iterations (iterations N, N -1, . . . , N +1-D) are also bounded by L(ǫ), because δ n max is non-increasing and our bound (Eq. 20) on the cost of one iteration is decreasing as a function of In the same way, the cost of iteration N -kd-i, for some k ∈ N and i

∈ [[0, D-1]], is bounded from above by L(ǫ)( 3 4 ) k if N -kd -i ≥ K 2 (D)
. The overall cost is therefore bounded from above by

DL(ǫ) × (1 + 3 4 + ( 3 4 ) 2 + ( 3 4
) 3 + . . . ).

Therefore the computational cost for ensuring

√ Dδ max N ≤ ǫ is O(DL(ǫ))
, and is therefore O(-log(ǫ)/ǫ). The tightness within logarithmic factors is proved, and Theorem 2 is proved.

Further remarks:

-We have not optimized the dependency in D. Our algorithm does not need that the fitness values at a given point is Bernoulli; it works for any distribution with values between 0 and 1 and variance upper bounded by the expectation. -Our algorithm uses a specific form of Bernstein race, in which ǫ 1 and ǫ 2 are known in advance (i.e. from δ max n we can guess an upper bound on the fitness values in [x - n , x + n ] and a lower bound on the difference between the best and the worst arm). It means that the algorithm works e.g. for f noise (x) a Bernoulli random variable B(||x -x * ||) with parameter ||x -x * ||, but not for f noise (x) = B(2||x -x * ||) or f noise (x) = B(k||x -x * ||) for arbitrary k > 0. Fortunately, [START_REF] Mnih | Empirical Bernstein stopping[END_REF] provides other Bernstein-based races, bearing no significant increase of the computational cost without this prior knowledge; therefore, replacing our simple Bernstein race by one of these results in the same soundness, and the same consistency; the convergence rate is only modified in the sense that it holds for a fixed D and a fixed k. This is a simple modification; a less trivial further work consists in analyzing this algorithm with a wider family of fitness functions; we guess that the soundness and consistency are preserved for any fitness function with values in [0, 1] with Ef noise = g(||x -x * ||) with g some increasing function from [0, ∞[ to [0, 1], and optimal complexity (within logarithmic factors) O(1/ǫ 2 ) for e.g. g(t) = c + t for some c > 0.

Conclusion

The key result of this paper is that including a bandit as in [START_REF] Heidrich-Meisner | Hoeffding and bernstein races for selecting policies in evolutionary direct policy search[END_REF] in the evaluation step of an evolution strategy is a good idea. The underlying subtlety is that we must not find all the µ best points as in [START_REF] Heidrich-Meisner | Hoeffding and bernstein races for selecting policies in evolutionary direct policy search[END_REF]; when trying to find all the µ best points (even if µ = 1!) one might spend of huge computational effort for separing two points with very close fitness. This happened with several tested variants. It is much more stable to generate an offspring in a controlled manner-so that there's no risk of having all points with almost equal fitness-and then to try and ensure that one of the points is better than another one, nothing more. Further work. A first extension to the results presented in this paper consists in generalizing the result for including "natural" algorithms; our algorithm is somehow reasonnable, but implies specifically derandomized offspring; practitionners might prefer more classical mutation algorithms, with an ad hoc bandit algorithm for avoiding troubles.

Another further work has been outlined after Theorem 2: considering a Bernstein race without prior on the shape of bandits, with empirical variance estimates instead. This is easy to do in our setting, but a non trivial extension would be the use for cases in which the variance of the noise around the optimum does not converge to zero.

We have shown a lower bound on the complexity of optimization of noisy fitness functions, with assumptions on the variance of noise around the optimum. There is no log-linear convergence (for convergence rate normalized by the number of function evaluations), but instead a linear decrease of the precision as a function of the number of iterations. We guess that the lower bound becomes the square root of the inverse of the number of evaluations if the variance does not decrease to zero around the optimum-proving it is another necessary further work.

Parallelization should provide a much better speed-up in the case of noisy optimization than in deterministic frameworks. It is therefore quite appealing to extend the above results for more than 1 computation unit.

A

We will show that given a random variable X whose probability distribution's support is in [0,1], its variance is always bounded by its mean. Let us assume that the distribution admits a density p. Then, define α = E[X] = Since ∀x ∈]0, 1[, x > x 2 and p(x) ≥ 0, E[X] -E[(X -E[X]) 2 ] is positive, which proves the original claim for variables admitting a density. Furthermore, it is clear that the result is also true for random variables who do not admit a density, by seeing their probability distribution as a limit of a sequence of probability distributions admitting densities, and applying the above reasoning.

n

  than to another x ′ j n (by Bernstein race); -Splitting the domain by the hyperplane in the middle of these points and normal to the line they define, and keeping only the side of the domain containing x ′ i n . Related works. Our "domain reduction" step relies on a good point (x

Fig. 1 .

 1 Fig.1. Noisy optimization algorithm (cf Algorithm 2). At each iteration, a main axis is selected (the one on which the domain has maximum range). Three equally spaced points are generated in the domain on this axis (this is the offspring). Then, a Bernstein race is applied for choosing a "good" and a "bad" arm among these points. The domain is reduced thanks to this knowledge, removing one fourth or one half of the domain (depending on the position of the good arm and of the bad arm -the best case is when the good and the bad arm are diametrically opposed: see Fig.2).

Fig. 2 .

 2 Fig. 2. The large rectangle is the domain [x - n , x + n ]. The three circles are arms x ′1n , x ′2 n , x ′3 n ; the left arm is the "bad" arm, whereas the arm in the center is the "good" arm, i.e. the one which proved to be closer to the optimum than the left arm, with confidence 1 -6δ/(π 2 n 2 ).

.

  Therefore, the cost for these D iterations is at most DL(ǫ). Each of the D previous iterations (iterations N -D, N -1 -D, . . . , N + 1 -2D) has computational cost bounded by ( 3 4 )L(ǫ), because -for n ′ ≤ n -D, δ n ′ max ≥ 4 3 ǫ; -the cost (Eq. 20) is superlinear as a function of ǫ 1 = √ Dδ n max .

1 0 1 0

 11 xp(x)dx. E[X] -E[(X -E[X] 2 ] = (x -(x -α) 2 )p(x)dx = (1 + 2α)α -α 2x 2 )p(x)dx

  we pick up the coordinate with highest uncertainty δ max

	n

  except that we modify the c th coordinate

	end for ǫ1 = √ Dδ max n	.	// this is a bound on the expected fitness of the x ′i n
	ǫ2 = 2K0δ max n	.	// this is a lower bound on the difference between
				// the best and the worst arm
	(good, bad)

  If x * ∈ [x - Consider the case in which consistency holds (we have proved, in Proposition 1, that this holds with probability 1 -δ.Then, if n is the last iteration in which the Bernstein race is applied, δ max This equation, combined with Eq. 20, gives an upper bound L(ǫ) as follows on the complexity (in terms of the number of fitness evaluations) of the last call to Bernstein(.):

	N , x + N ] and N x * || ≤ ǫ; therefore, we have reached precision ǫ if δ max √ Dδ max N ≤ ǫ, then necessarily ||x -N √ D and the ≤ ǫ/ consistency holds. n ≥ ǫ/ √ D.

Note that in our problem setting, the function g = 1 3 +||x-t|| should in all generality range in [0, 1], which is not the case here. We can either set the measurement to 0 when g(x, t) > 1, or consider f (x, t) = min(1, g(x, t)) instead, which is what will be done from now on.

Note that this formulation is equivalent to the problem setting described at the beginning of section 2, with x * in the role of t; it is more suited to describe our algorithm and proof.

Acknowledgements This work was supported by the French National Research Agency (ANR) grant No. ANR-08-COSI-007-12, and through COSINUS program (project EXPLO-RA nANR-08-COSI-004). We are very grateful to Hervé Fournier for his valuable comments on a preliminary version of this paper, and to Anne Auger for fruitful discussions around noisy optimization. We are also grateful to the people of the "MoGo" project for joint work around noisy optimization.


Ege Gursoy

Akansel Cosgun

Dana Kulic

Andrea Cherubini

PICKBOT: Dual-arm Harvesting Robot

Objective

Develop an efficient robotic system that uses multiple tools to explore the potential of dual-arm robots in agricultural applications.

Challenges

Different types of fruits need different approaches Fruits are unique in color, shape, weight etc.

Highly unstructured and dynamic environment

Important metrics to evaluate the system are speed, accuracy and damage to the fruit/plant

Future work

Explore more fruits and tools Avoid self-occlusions Validate in outdoor setups Improve planning: approaching angles, branch avoidance

Fruit

Localization

Use case: Orange harvesting [START_REF] Wang | Improved YOLOv5 network for real-time multi-scale traffic sign detection[END_REF] LIRMM UM/CNRS, Montpellier, France 2 Monash University, Melbourne, Australia YOLO algorithm [START_REF] Wang | Improved YOLOv5 network for real-time multi-scale traffic sign detection[END_REF] trained on 300 images/fruit [START_REF] Pawara | Comparing local descriptors and bags of visual words to deep convolutional neural networks for plant recognition[END_REF] performs 2D detection Projection using an RGB-D camera gives us fruit 3D locations Real-time tracking algorithm [START_REF] Wojke | Simple online and realtime tracking with a deep association metric[END_REF] improves robustness of the position

We send the position of the fruit which is closest to the robot

The online trajectory generator [START_REF] Kröger | Opening the door to new sensor-based robot applications-The Reflexxes Motion Libraries[END_REF] guides the robot to the fruit position in real-time while avoiding self collisions [START_REF] Benallegue | Fast C 1 proximity queries using support mapping of sphere-torus-patches bounding volumes[END_REF] The absolute -Relative task description ensures safe dual-arm manipulation Hierarchical Quadratic Programming based controller [START_REF] Tarbouriech | An admittance based hierarchical control framework for dual-arm cobots[END_REF] let us define priorities between tasks

Control

The BAZAR robot [START_REF] Cherubini | A collaborative robot for the factory of the future: Bazar[END_REF] has two 7 DOF arms, 1 RGB-D camera mounted on the head and a mobile base Two primary tools for harvesting are electric scissors and a fruit basket, which can be placed on any arm according to the needs.

After receiving the position of the fruit, the cutting arm aims at the peduncle while the collecting arm positions itself relatively to the cutting one

We repeat the harvesting loop each time a new fruit position is received

Absolute task

Relative task


Vu-Linh Nguyen
email: vu-linh.nguyen@hds.utc.fr

Haifei Zhang
email: haifei.zhang@hds.utc.fr

Sébastien Destercke
email: sebastien.destercke@hds.utc.fr

Learning Sets of Probabilities Through Ensemble Methods

Keywords: Ensemble learning, Credal sets approximation, Set-valued prediction, Quantile-based approach

A possible approach to obtain set-valued predictions is to learn for each query instance a probability set (a.k.a. credal set) representing its associated uncertainty. Theoretically founded decision rules extending classical expectation and inducing a partial order between predictions can the be used to derive set-valued predictions. However, obtaining such a credal set by imprecisiating a given learning algorithm is usually computationally challenging, except for simple models such as decision trees or naive Bayes classifiers. In this paper, we propose a simple, easy to use quantile-based framework for estimating credal sets using output of ensemble methods, that can also cope with complex types of data, such as images and mixed/multimodal data, etc. Experiments are conducted to highlight the usefulness of the proposed framework.

Introduction

Classification algorithms are usually designed to produce, for each instance, a prediction in the form of a unique element of the set of possible outputs. Under the presence of uncertainty, which is often a consequence of model inadequacy and/or data imperfections (in terms of quality and/or quantity), the model can however be uncertain about its predictions and make unreliable precise predictions. In such a case, it might be more desirable to provide imprecise (or indeterminate) set-valued predictions which aims to balance correctness (the true output is an element of the set-valued prediction) and precision (the cardinality of the set-valued prediction) in some appropriate manner [START_REF] Del Coz | Learning nondeterministic classifiers[END_REF][START_REF] Mortier | Efficient set-valued prediction in multi-class classification[END_REF][START_REF] Troffaes | Decision making under uncertainty using imprecise probabilities[END_REF][START_REF] Zaffalon | Evaluating credal classifiers by utilitydiscounted predictive accuracy[END_REF].

Learning with a reject option is the simplest case of learning set-valued predictions, in which the classifier is allowed to either produce a singleton prediction or refuse to make a prediction for a given query instance. Threshold-based classifiers have been proposed for that purpose, in which a (global/local) threshold will be employed to decide whether a query instance should be rejected or predicted and then a conventional classifier is called only if the instance should be classified [START_REF] Bartlett | Classification with a reject option using a hinge loss[END_REF][START_REF] Chow | On optimum recognition error and reject tradeoff[END_REF][START_REF] Cortes | Learning with rejection[END_REF][START_REF] Franc | On discriminative learning of prediction uncertainty[END_REF][START_REF] Grandvalet | Support vector machines with a reject option[END_REF][START_REF] Hellman | The nearest neighbor classification rule with a reject option[END_REF]. Threshold-based classifiers have been developed for multi-class classification (MCC) [START_REF] Del Coz | Learning nondeterministic classifiers[END_REF][START_REF] Mortier | Efficient set-valued prediction in multi-class classification[END_REF], when the classifier is allowed to return top (locally/globally) ranked classes. While such classifiers are intuitive and easy to implement, they often require reliable estimates of the class probabilities to be performant, which is hard to ensure when information is lacking.

By considering more expressive uncertainty representations, imprecise probabilistic classifiers [START_REF] Corani | Learning reliable classifiers from small or incomplete data sets: The naive credal classifier 2[END_REF][START_REF] Cozman | Credal networks[END_REF][START_REF] Mantas | Credal-c4. 5: Decision tree based on imprecise probabilities to classify noisy data[END_REF][START_REF] Zaffalon | The naive credal classifier[END_REF] can provide, at least in theory, more reliable outputs. They are developed based on the assumption that uncertainty is described by a (not necessarily convex) set of probabilities, i.e., a credal set [START_REF] Levi | The enterprise of knowledge: An essay on knowledge, credal probability, and chance[END_REF], a description to which can then be applied theoretically justified decision rules [START_REF] Jansen | Quantifying degrees of e-admissibility in decision making with imprecise probabilities[END_REF][START_REF] Troffaes | Decision making under uncertainty using imprecise probabilities[END_REF] to produce set-valued predictions. Moving from a single distribution to a credal set is a natural way to model the lack of information, an aspect that unique probabilities can hardly capture. Unfortunately, imprecise probabilistic classifiers often suffer from the limited use to certain types of (tabular) data, as well as from the high computational cost that represent a credal extension of a given learning method. A solution might be to consider the credal set as a neighbourhood of the initial estimated distribution [START_REF] Montes | Unifying neighbourhood and distortion models: part i-new results on old models[END_REF][START_REF] Rahimian | Distributionally robust optimization: A review[END_REF], yet ensuring the quality of the initial estimated distribution is a challenge itself.

In this paper, we propose a quantile-based framework for estimating credal sets from the output of ensembles [START_REF] Dietterich | Ensemble methods in machine learning[END_REF]. We specifically seek a correctness-precision trade-off when constructing estimates of credal sets, i.e., the estimates are expected to be informative and at the same time not very large. This shall be done by defining "median" of set of distributions and use the "median" to filter out a proportion of "extreme" distributions before forming credal sets. Moreover, we only require the availability of an ensemble of probabilistic classifiers. Thus, the base learner (ensemble) can be freely chosen according to our needs. This flexibility of the proposed approach is remarkably different from existing imprecise probabilistic classifiers. Therefore, we hope to broaden the use of generalized decision rules [START_REF] Jansen | Quantifying degrees of e-admissibility in decision making with imprecise probabilities[END_REF][START_REF] Troffaes | Decision making under uncertainty using imprecise probabilities[END_REF] to applications with complex types of data, such as mixed data [START_REF] De Leon | Classification with discrete and continuous variables via general mixed-data models[END_REF], image/video [START_REF] Yang | Tensor-train recurrent neural networks for video classification[END_REF][START_REF] Yin | Towards efficient tensor decompositionbased dnn model compression with optimization framework[END_REF] and multimodal data [START_REF] Ngiam | Multimodal deep learning[END_REF][START_REF] Xu | Mufasa: Multimodal fusion architecture search for electronic health records[END_REF].

We provide in Section 2 a minimal description of MCC with sets of probabilities. Our main contribution which is a quantile-based approach for estimating credal sets is presented in Section 3. The inference problem with sets of probabilities is summarized in Section 4. Section 5 presents some preliminary experiments on tabular data sets to motivate the use of the proposed framework. Section 6 concludes this work and sketches out future work.

Preliminary

We shall recall basics of classification with sets of probabilities and notations.

Probabilistic Classification

Let X denote an instance space, and let Y = {y 1 , . . . , y K } be a finite set of classes. We assume that an instance x ∈ X is (probabilistically) associated with members of Y. We denote by p(Y | x) the conditional distribution of Y given X = x. Given training data D = {(x n , y n )|n = 1, . . . , N } drawn independently from p(X, Y ), the goal in MCC is to learn a classifier h, which is a mapping X -→ Y that assigns to each instance x ∈ X a class ŷ . .= h(x) ∈ Y.

To evaluate the performance of a classifier h, a loss function : Y ×Y -→ R + is needed, which compares a prediction ŷ with a ground-truth y. Each classifier h is evaluated using its expected loss

R(h) . .= E (Y, h(X)) = (y, h(x)) d P(x, y) ,
where P is the joint probability measure on X × Y characterizing the underlying data-generating process. Therefore, the Bayes-optimal classifier is given by

h * ∈ argmin h∈H R(h) , (1)
where H ⊆ Y X is the hypothesis space. When H is probabilistic, we can follow maximum likelihood estimation and define the Bayes-optimal classifier as the classifier which optimizes the conditional log likelihood (CLL) function:

ĥ . .= p ∈ argmax p∈H CLL(p | D) . .= argmax p∈H 1 N N n=1 log p(y n | x n ) . (2)
To avoid overfitting, the CLL is often augmented by a regularization term [START_REF] Murphy | Machine learning: a probabilistic perspective[END_REF][START_REF] Nguyen | Probabilistic multi-dimensional classification[END_REF].

Once the classifier (2) is learned from D, we can in principle find an optimal prediction of any loss function at the prediction time [START_REF] Elkan | The foundations of cost-sensitive learning[END_REF][START_REF] Mortier | Efficient set-valued prediction in multi-class classification[END_REF]. More precisely, assume the classifier (2) is made available, and predicts for each query instance x a probability distribution p(• | x) on the set of labelings Y. The Bayes-optimal prediction (BOP) of any is then given by the expected loss minimizer

ŷ = ŷ(x) ∈ argmin ȳ∈Y E (y, ȳ) = argmin ȳ∈Y y∈Y (y, ȳ) p(y | x) .
(3)

Classification with Set of Probabilities

Under this setting, we assume that our uncertainty is described by a (not necessarily convex) set of probabilities P(Y | x), i.e., a credal set [START_REF] Levi | The enterprise of knowledge: An essay on knowledge, credal probability, and chance[END_REF]. Clearly, the decision rule ( 3) is no longer directly applicable. Therefore, it is necessary to use some generalized decision rule such as the ones benefiting from strong theoretical justifications [START_REF] Jansen | Quantifying degrees of e-admissibility in decision making with imprecise probabilities[END_REF][START_REF] Troffaes | Decision making under uncertainty using imprecise probabilities[END_REF].

Credal sets can arise in different ways, either as a native result of the learning method [START_REF] Augustin | Introduction to imprecise probabilities[END_REF], as the result of an agnostic (with respect to the missingness process) estimation in presence of imprecise data, or as a neighbourhood taken over an initial estimated distribution p(Y | x) [START_REF] Montes | Unifying neighbourhood and distortion models: part i-new results on old models[END_REF][START_REF] Rahimian | Distributionally robust optimization: A review[END_REF]. These approaches seems to introduce some inconvenience. Native credal classifiers can be hard to learn, and are unavailable for complex inputs such as such as mixed data and images. Approximating P(Y | x) as a neighbourhood taken over an initial estimated distribution p(Y | x) does not face this inconvenience, but requires that the initial estimated distribution is well-estimated, a hard to ensure quality.

In the next section, we propose a simple, flexible and easy to use quantilebased framework for estimating credal sets using output of ensemble methods [START_REF] Dietterich | Ensemble methods in machine learning[END_REF]. This is especially designed to make use of the current and future development of both probabilistic classification and generalized decision rules in a unified framework to broaden the application of imprecise probability (IP) to real-world applications with complex data types.

Credal Sets Approximation

We assume an ensemble

H . .= {h m | m ∈ [M ] . .= {1, . . . , M }} of M probabilistic classifiers h m , m ∈ [M ]
is made available and provides, for each instance x, a set of M probabilistic predictions

H(x) . .= {h m (x) | m ∈ [M ]} = {p m . .= (p m 1 , p m 2 , . . . , p m K ) | m ∈ [M ]} . ( 4
)
Our goal is to aggregate this set of probabilistic predictions into a credal set P(Y | x) in some meaningful way.

A Quantile-Based Approach

The intention of this approach is to seek a correctness-precision trade-off, i.e., the estimations of P(Y | x) are expected to be informative and at the same time not very large. We define the reference point of H(x) as follows:

p * = argmin p: K k=1 p k =1 M m=1 d(p, p m ) . ( 5
)
where d is some distance defined for pairs of probability distributions.

Once the reference point p * is made available, it allows us to define a preference order, reflecting how common/weird each distribution in H(x) is:

p p if d(p * , p) < d(p * , p ) . ( 6
)
Such a preference order in turn allows us to "discard" a given percentage of outliers among elements of H(x).

Let α ∈ [0, 1] be some threshold. We define H α (x) as the set of (1 -α) * 100 % of closest distributions in H(x) with respect to the preference order [START_REF] Corani | Learning reliable classifiers from small or incomplete data sets: The naive credal classifier 2[END_REF]. We approximate the credal set p(Y | x) of x by the convex hull of

H α (x). Let H α (x) . .= {p m | m ∈ [M α ]}. The convex hull is defined as CH α (x) . .= p . .= Mα m=1 γ m p m | γ m ≥ 0, m ∈ [M α ], Mα m=1 γ m = 1 . (7)
The computational complexity of the problem of determining the reference point (5) can greatly depend on the nature of the distance d. In the next section, we recall commonly used distances. Due to page length limit, we only mention few convex distances and refer to [START_REF] Cha | Comprehensive survey on distance/similarity measures between probability density functions[END_REF][START_REF] Gibbs | On choosing and bounding probability metrics[END_REF][START_REF] Lee | Measures of distributional similarity[END_REF][START_REF] Sriperumbudur | Hilbert space embeddings and metrics on probability measures[END_REF] for more distances.

The Cases of Convex Distances

For completeness, we shall start with few definitions and remarks, which are quite basic and would have appeared in textbooks and papers (see, e.g., [START_REF] Boyd | Convex optimization[END_REF][START_REF] Datta | Numerical linear algebra and applications[END_REF][START_REF] Pugh | Real mathematical analysis[END_REF]).

Definition 1. A function f : R K -→ R is convex if for every p, p ∈ R K and every λ 1 , λ 2 ∈ [0, 1] such that λ 1 + λ 2 = 1, we have the inequality f (λ 1 p + λ 2 p ) ≤ λ 1 f (p) + λ 2 f (p ) . ( 8
) Remark 1. Let z ∈ R K . Let • be a norm on R K . f (p) . .= p -z is convex.
Proof. The convexity of f (p) follows consequently from the triangle inequality of norms:

f (λ 1 p + λ 2 p ) = λ 1 p + λ 2 p -z = λ 1 (p -z) + λ 2 (p -z) ≤ λ 1 (p -z) + λ 2 (p -z) = λ 1 p -z + λ 2 p -z = λ 1 f (p) + λ 2 f (p ) .
Remark 2. Conical combinations of convex functions are also convex.

Proof. The proof is trivial. It is enough to multiply the inequalities, one per convex function, by non-negative scalars and sum them up.

In the following, we show that if

f m (p) . .= d(p, p m ) is convex, m ∈ [M ]
, then the problem of finding a reference point (5) of H(x) can be straightforwardly formulated as a convex optimization problem. This is indeed computationally advantageous because with recent advances, convex programming is nearly as straightforward as linear programming [START_REF] Boyd | Convex optimization[END_REF][START_REF] Rockafellar | Lagrange multipliers and optimality[END_REF].

Definition 2. A standard convex optimization problem is of the form

minimize p f (p) subject to g i (p) ≤ 0 , i ∈ [I] , h j (p) = 0 , j ∈ [J] (9)
where: p ∈ R K is the optimization variable; The objective function f : R K -→ R is convex; The inequality constraint functions

g i : R K -→ R, i ∈ [I] are convex;
The equality constraint functions h j : R K -→ R, j ∈ [J], are of the form:

h i (p) = a j p -b j
, where a j is a vector and b j is a scalar.

We can encode the condition that the reference point must be a valid probability distribution by using K inequality constraint functions g i and 1 equality constraint function h 1 :

g k (p) . .= -p k ≤ 0 , k ∈ [K] , h 1 (p) . .= 1 K p -1 = 0 , (10)
where

1 K = (1, . . . , 1). The constraints p k ≤ 1, k ∈ [K],
are implicitly enforced by the K constraints g k (i.e., p k ≥ 0, k ∈ [K]) and h 1 (i.e.,

K k=1 p k = 1, k ∈ [K]
). Therefore, we can use any existing package to find p * (5).

Using Remark 1-2, we can verify that different distances (See [START_REF] Cha | Comprehensive survey on distance/similarity measures between probability density functions[END_REF][START_REF] Gibbs | On choosing and bounding probability metrics[END_REF][START_REF] Lee | Measures of distributional similarity[END_REF][START_REF] Sriperumbudur | Hilbert space embeddings and metrics on probability measures[END_REF] and elsewhere) are convex. Examples are members of the L p Minkowski family

f p (p) . .= L p (p, z) . .= p K k=1 |p k -z k | p , p ≥ 1 , (11)
and Chebyshev distance

f cheb (p) . .= L ∞ (p, z) . .= max k∈[K] |p k -z k | . (12)
Moreover, a closer look at Definition 1 is enough to verify the convexity of some other distances (discussed in [START_REF] Cha | Comprehensive survey on distance/similarity measures between probability density functions[END_REF][START_REF] Gibbs | On choosing and bounding probability metrics[END_REF][START_REF] Lee | Measures of distributional similarity[END_REF][START_REF] Sriperumbudur | Hilbert space embeddings and metrics on probability measures[END_REF] and elsewhere). Examples are the Squared Euclidean distance (whose square function allows triangle inequality)

f sqe (p) . .= d sqe (p, z) . .= K k=1 (p k -z k ) 2 , (13)
and KL divergence (inequality (8) can be verified using the log sum inequality):

f KL (p) . .= d KL (p, z) . .= K k=1 p k log (p k /z k ) , . (14)
To solve the problem (9) efficiently, one should carefully look at the nature of the given convex distance. For example, for any given K, closed-form solution for the f Sqe (13) can be derived (See Proposition 1). This is also a special case where the additional constraints (i.e., K k=1 p k = 1 and p k ≥ 0, k ∈ [K]) do not change the minimizer. However, it is not always the case. For example, these additional constraints can change the minimizer of f 1 (11) (See Proposition 2). Also, different distances may seek for the same minimizer. Examples of such distances are Topsør and Jensen-Shannon [START_REF] Cha | Comprehensive survey on distance/similarity measures between probability density functions[END_REF]. Moreover, for some distance, such as Inner Product [START_REF] Cha | Comprehensive survey on distance/similarity measures between probability density functions[END_REF], the problem ( 9) is reduced to a linear program. Proposition 1. The reference point p * (5) under Squared Euclidean distance f sqe [START_REF] Elkan | The foundations of cost-sensitive learning[END_REF] is uniquely defined as

p * k = 1 M M m=1 p m k , k ∈ [K] . ( 15
)
Proof. The proof is trivial and is given for completeness. For any k ∈ [K], the partial derivative of

f (p) = M m=1 f m Sqe (p) = K k=1 M m=1 (p k -p m k ) 2 (16)
with respect to the variable p k is

∂f ∂p k (p) = 2 M m=1 (p k -p m k ) = 2 M p k - M m=1 p m k . (17)
Since f sqe (p) (13) is strictly convex, its unique minimizer is attained when the partial derivatives are zeros, i.e., p * is defined in [START_REF] Gibbs | On choosing and bounding probability metrics[END_REF]. p * is a valid distribution because the set of possible distributions is a convex set.

Proposition 2. Except for K = 2, the reference point p * (5) under f 1 [START_REF] Del Coz | Learning nondeterministic classifiers[END_REF] may not be the minimizer of the relaxed optimization problem

p ∈ argmin p M m=1 L 1 (p, p m ) = argmin p K k=1 M m=1 |p k -p m k | . ( 18
)
Proof. Without enforcing the probability axioms (i.e.,

K k=1 p k = 1 and p k ≥ 0, k ∈ [K]
), a minimizer p of the relaxed optimization problem ( 18) is defined as

pk . .= median(p 1 k , . . . , p M k ) , k ∈ [K] . (19)
This can be verified by showing that, for any p = p, we have

f 1 (p k ) . .= M m=1 |p k -p m k | ≥ M m=1 |p k -p m k | := f 1 (p k ) , k ∈ [K] , (20)
which implies the relation f 1 (p) ≥ f 1 (p). Let L k be the number of p m k which is larger than pk . Let S k be the number of p m k which is smaller than pk . By definition of "median", we have L k = S k . p k > pk : We have the following relations

|p k -p m k | = |p k -p m k | + |p k -pk | if p m k ≤ pk , |p k -p m k | ≥ |p k -p m k | -|p k -p k | if p m k ≥ pk .
Therefore, we have

f 1 (p k ) = M m=1 |p k -p m k | ≥ M m=1 |p k -p m k | + |p k -pk |S k -|p k -pk |L k = M m=1 |p k -p m k | + |p k -pk |(S k -L k ) = M m=1 |p k -p m k | = f 1 (p k ) .
p k < pk : We have the following relations

|p k -p m k | = |p k -p m k | + |p k -pk | if p m k ≥ pk , |p k -p m k | ≥ |p k -p m k | -|p k -p k | if p m k ≤ pk .
Therefore, we have

f 1 (p k ) = M m=1 |p k -p m k | ≥ M m=1 |p k -p m k | + |p k -pk |L K -|p k -pk |S k = M m=1 |p k -p m k | + |p k -pk |(L k -S k ) = M m=1 |p k -p m k | = f 1 (p k ) .
For K > 2, p may not satisfy the probability axioms (see next Table ). When K = 2, the probability axioms of p are ensured by the fact that the total rank of each distribution p m , m ∈ [M ], on the first and the second classes is always M +1 (as the masses should sum up to 1). Thus, p is either one element of H(x) or the average of two elements of H(x). Let us illustrate this property using an example where M = 9: p 1 p 2 p 3 p 4 p 5 p 6 p 7 p 8 p 9 p 1 Value 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 Rank 1 2 3 4 5 6 7 8 9

K = 3 K > 3 p 1 0.8 0.1 0.1 p 1 0.4 0.2 0.4/(K-3) . . . 0.4/(K-3) p 2 0.
p 2 Value 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Rank 9 8 7 6 5 4 3 2 1

In this example, the total rank is 10 and p is p 5 .

In the next section, we recall the inference problem with credal sets [START_REF] Jansen | Quantifying degrees of e-admissibility in decision making with imprecise probabilities[END_REF][START_REF] Troffaes | Decision making under uncertainty using imprecise probabilities[END_REF].

Inference Problem

As said, when our uncertainty is described by a credal set P(Y | x), instead of a single probability p(Y | x), it is necessary to make predictions using some theoretically founded decision rule extending classical expectation [START_REF] Jansen | Quantifying degrees of e-admissibility in decision making with imprecise probabilities[END_REF][START_REF] Troffaes | Decision making under uncertainty using imprecise probabilities[END_REF]. For any p ∈ P(Y | x) and any loss function , we shall denote by

ŷp ∈ argmin ȳ∈Y y∈Y (y, ȳ) p(y | x) . (21)
Definition 3. An optimal set-valued prediction under the E-admissibility rule is

ŶE ,P = {y ∈ Y | ∃ p ∈ P s.t. y = ŷp } . (22)
Definition 4. An optimal set-valued prediction ŶM ,P under the Maximality rule is the set of the maximal, non-dominated elements of the partial order π P such that ȳ ,P ȳ if

inf p∈P Ep ( (y, ȳ ) -(y, ȳ)) > 0 . ( 23
)
In other words, we have

ŶM ,P = {ȳ ∈ Y | ∃ ȳ s.t. ȳ ,P ȳ} . (24)
It is known that the set-valued prediction given by the E-admissibility rule is a subset of the one given by the Maximality rule [START_REF] Troffaes | Decision making under uncertainty using imprecise probabilities[END_REF].

In the following, we discuss the computational complexity of the inference problem when is the 0/1 loss, i.e., (y, ȳ) = y = ȳ , where A = 1 if the predicate A is true and equals 0 otherwise Let us start with the case of Maximality rule. For any p ∈ CH α (x), we have

Ep ( (y, ȳ ) -(y, ȳ)) = p(ȳ | x) -p(ȳ | x) . (25)
Thus, the relation ȳ ,P ȳ holds if the maximum of the linear program

maximize p f (p) . .= p(ȳ | x) -p(ȳ | x) (26)
subject to p - Mα m=1 γ m p m = 0, γ m ≥ 0, Mα m=1 γ m = 1 , (27)
is negative. Note that if f (p) has a maximum value on the feasible region, then it has this value on (at least) one of the extreme points, i.e., elements of H α (x) [START_REF] Murty | Linear programming[END_REF][Theorem 3.3]. Thus, a naive algorithmic solution is to compute f (p) for the extreme p and compare it with 0. This requires time O(K 2 M α ) because in the worst case, one needs to check all the K(K -1) relation ȳ ,P ȳ , ȳ = ȳ ∈ Y.

We now tackle the case of the E-admissibility rule. Reminding that, ∀ y ∈ ŶE ,P , there must exist at least one p ∈ CH α (x) such that y = ŷp . This is equivalent to having at least one p ∈ CH α (x) such that p(y | x) ≥ p(y | x), y = y. Thus, given any outer approximation Y O ,P of ŶE ,P we can follow the suggestion of [START_REF] Jansen | Quantifying degrees of e-admissibility in decision making with imprecise probabilities[END_REF]

γ m = 1 , (29)
p(y | x) -p(y | x) ≥ 0 , y ∈ Y \ {y, y } , (30)
where y = y, and comparing it with 0. Hence, finding ŶE ,P requires solving |Y O ,P | linear programs, one per y ∈ |Y O ,P |. The naive algorithmic solution, i.e., iterating over all the extreme points, can not be applied here because a class y may be optimal only for probabilities in the interior of CH α (x).

Experiment

To motivate the potential use of the proposed framework, we perform some experiments on 9 tabular datasets from the UCI repository (cf. the left part of Table 1), following a 10-fold cross-validation procedure.

We employ random forests (RFs) [START_REF] Ho | Random decision forests[END_REF] (with default setting of scikit-learn) as the base learner. RFs are compared to an instantiation of our framework, where H α (x) is constructed under the f sqe [START_REF] Elkan | The foundations of cost-sensitive learning[END_REF] and used to produce ŶM ,P . For each train test split, we follow a 10-fold nested cross-validation procedure to choose α optimizing u 65 . The RF is then retrained using the entire training dataset and the chosen α is used to construct H α (x) during the inference phase. The source code has been made public at https://github.com/Haifei-ZHANG/Probability-Sets-Model. Overall results (accuracy, u 65 and u 80 scores [START_REF] Zaffalon | Evaluating credal classifiers by utilitydiscounted predictive accuracy[END_REF] and cardinality | ŶM ,P |) show that our proposal may provide a promising correctness-precision trade-off, compared to RFs itself. Ideally, a reliable classifier should be more cautious on difficult cases, on which the conventional classifier is likely to fail [START_REF] Nguyen | Reliable multi-class classification based on pairwise epistemic and aleatoric uncertainty[END_REF][START_REF] Yang | Nested dichotomies with probability sets for multi-class classification[END_REF]. To verify this ability of our proposal, for each dataset, we report the correctness (i.e., the percentage of times the true class is in ŶM ,P ), given the prediction was imprecise, versus the accuracy of RF on those instances. The results (in the right part of Table 1) strongly support our proposal.

This also suggests the use of the E-admissibility rule (listed as future work) may improve the overall results because, under the f sqe , predictions of RFs should belong to ŶE ,P ⊂ ŶM ,P [START_REF] Troffaes | Decision making under uncertainty using imprecise probabilities[END_REF]. More precisely, under the f sqe , our proposal should always gain in the term of correctness1 and the use of the E-admissibility rule may help to produce smaller (reliable) imprecise predictions.

To gain more insights about the influence of α, we consider u 65 and u 80 scores on the test set as functions of the value of α. The results in Figure 1 are indeed in agreement with our expectations. The P(Y | x) induced by H α (x) with small α may contain extreme/noisy distributions and produce large ŶM ,P (resulting in low u 65 and u 80 scores). Moderate α may provide a nice correctness-precision trade-off (reflected via promising u 65 and u 80 scores). For large α, P(Y | x) is shrunken as (small) neighborhood of the p * (5) and our proposal (under f sqe ) becomes similar to RFs, which use the p * to make predictions. The results also suggest that, in practice, nested cross-validation procedure can help us to find some good value of α (even if the ideal gain provided by the optimal α is small).

Conclusion

We propose a simple, easy to use quantile-based framework for estimating credal sets using output of ensemble methods, that can also cope with complex types of data. Preliminary experiments suggest that our proposal may provide a promising correctness-precision trade-off, compared to ensemble methods. To seek for a complete picture on the usefulness of our proposal, we envision the following works: (1) implement our proposal with other distances and the E-admissibility rule and analyze (dis)advantages provided by different combinations of distance and decision rule, (2) include threshold-based classifiers as competitors, and (3) include complex types of data (such as images) into our empirical studies.

Our theoretical results also inform that voting ensembles, such as RFs, use the p * under f sqe to make predictions. It would be interesting to investigate whether using the p * under other distances to make predictions may bring significant difference, thought our primary goal is not to study the problem of how to aggregate the probabilistic predictions provided by ensemble members into the final singleton predictions.

  2 0.5 0.3 p 2 0.2 0.7 0.1/(K-3) . . . 0.1/(K-3) p 3 0.1 0.4 0.5 p 3 0.1 0.6 0.3/(K-3) . . . 0.3/(K-3) p 0.2 0.4 0.3 p 0.2 0.6 0.3/(K-3) . . . 0.3/(K-3)

  and formulate the problem of checking whether a given y ∈ Y O ,P satisfies the relation y ∈ ŶE ,P as finding the maximum value of a linear program maximize p f (p) . .= p(y | x) -p(y | x) (28) subject to p -Mα m=1 γ m p m = 0, γ m ≥ 0, Mα m=1

Fig. 1 .

 1 Fig. 1. u65 and u80 scores on the test set as functions of the value of α

Table 1 .

 1 Statistics of data sets (P is the number of features) and experimental results.

		Statistics of data sets	RF	Overall results Ours		Cases of abstention RF Ours
	#	Name	N P K Acc. (%) u65 (%) u80 (%) | ŶM ,P | Acc. (%) Corr. (%)
	1	ecoli	336 7 8	78.35	77.77 79.38 2.05	69.84	93.59
	2 balance scale 625 4 3	80.50	82.17 83.15 2.02	26.75	67.75
	3	vehicle	846 18 4	74.46	78.16 82.63 2.04	47.31	90.24
	4	vowel	990 10 11 65.35	65.89 68.71 2.05	41.05	71.80
	5 wine quality 1599 11 6	57.91	61.67 68.54 2.02	49.69	86.73
	6	optdigits	1797 64 10 96.95	97.03 97.19 2.03	50.74	80.19
	7	segment	2300 19 7	98.05	98.02 98.22 2.09	50.12	78.93
	8 waveform 5000 21 3	85.52	85.81 88.33	2	62.06	99.91
	9	letter	20000 16 26 96.57	96.58 96.64 2.03	34.33	81.71

Its predictions are identical to the ones provided by RFs in the cases of singleton/precise predictions. In the cases of imprecise predictions, its predictions cover the predictions provided by RFs.

Acknowledgement

This work was funded/supported by the Junior Professor Chair in Trustworthy AI (Ref. ANR-R311CHD).