Cyrille Kone
email: cyrille.kone@inria.fr

Emilie Kaufmann
email: emilie.kaufmann@univ-lille.fr

Laura Richert
email: laura.richert@u-bordeaux.fr

Adaptive Algorithms for Relaxed Pareto Set Identification

In this paper we revisit the fixed-confidence identification of the Pareto optimal set in a multi-objective multi-armed bandit model. As the sample complexity to identify the exact Pareto set can be very large, a relaxation allowing to output some additional near-optimal arms has been studied. In this work we also tackle alternative relaxations that allow instead to identify a relevant subset of the Pareto set. Notably, we propose a single sampling strategy, called Adaptive Pareto Exploration, that can be used in conjunction with different stopping rules to take into account different relaxations of the Pareto Set Identification problem. We analyze the sample complexity of these different combinations, quantifying in particular the reduction in sample complexity that occurs when one seeks to identify at most k Pareto optimal arms. We showcase the good practical performance of Adaptive Pareto Exploration on a real-world scenario, in which we adaptively explore several vaccination strategies against Covid-19 in order to find the optimal ones when multiple immunogenicity criteria are taken into account.

Introduction

In a multi-armed bandit model, an agent sequentially collects samples from several unknown distributions, called arms, in order to learn about these distributions (pure exploration), possibly under the constraint to maximize the samples collected, viewed as rewards (regret minimization). These objectives have been extensively studied for different types of univariate arms distributions [START_REF] Lattimore | Bandit Algorithms[END_REF]. In this paper, we consider the less common setting in which arms are multi-variate distributions. We are interested in the Pareto Set Identification (PSI) problem. In this pure exploration problem, the agent seeks to identify the arms that are (Pareto) optimal, i.e. such that their expected values for all objectives are not uniformly worse than those of another arm.

We formalize this as a fixed-confidence identification problem: in each round t the agent selects an arm A t using an adaptive sampling rule and observes a sample X t ∈ R D from the associated distribution. It further uses an adaptive stopping rule τ to decide when to stop sampling and output a set of arms S τ which is her guess for (an approximation of) the true Pareto set S . Given a risk parameter δ ∈ (0, 1), this guess should be correct with high probability, e.g. satisfy P( S τ = S ) ≥ 1 -δ for exact Pareto set identification, while requiring a small sample complexity τ . This generalizes the well-studied fixed-confidence Best Arm Identification (BAI) problem [START_REF] Even-Dar | Action elimination and stopping conditions for the multi-armed bandit and reinforcement learning problems[END_REF][START_REF] Kalyanakrishnan | PAC subset selection in stochastic multiarmed bandits[END_REF][START_REF] Garivier | Optimal best arm identification with fixed confidence[END_REF] to multiple objectives.

37th Conference on Neural Information Processing Systems (NeurIPS 2023).

Our motivation to study multi-objective adaptive identification stems from the design of adaptive earlyphase clinical trials. In phase I/II trials, the effects of a molecule in humans are explored, and several biological endpoints may be assessed at the same time as indicative markers of efficacy. In particular, in the context of vaccine development, early-phase trials usually assess multiple immunogenicity endpoints (i.e. various markers of the effects of the vaccine on the immune system, such as different facets of antibody responses or other immune parameters). In the absence of a known correlate of protection during early clinical development, these endpoints may not have a clear a priori hierarchy, may not all be correlated, which makes an examination of the Pareto set of different vaccinal strategies particularly relevant. In addition, given the availability of various vaccine platforms (such as mRNA vaccines, viral-vector vaccines, protein vaccines), as exemplified by Covid-19 vaccines, there may be a need to adaptively screen the various resulting vaccine strategies to select the most promising ones. Apart from clinical trials, active Pareto Set Identification can be meaningful in many real-word contexts, and we refer the reader to the various examples given by [START_REF] Zuluaga | Active learning for multi-objective optimization[END_REF], such as in hardware or software design. Other applications include A/B/n testing for marketing or online recommender systems in which it is common to jointly optimize multiple (possibly conflicting) objectives such as user behavioral metrics (e.g. clicks, streams, dwell time, etc), supplier exposure objectives (e.g. diversity) and platform centric objectives (e.g. promotions) [START_REF] Mehrotra | Bandit based optimization of multiple objectives on a music streaming platform[END_REF].

For many applications, the sample complexity of exact PSI can be prohibitive, either when there are many close to optimal arms or when the Pareto set is very large, and different relaxations have been considered in the literature [START_REF] Auer | Pareto front identification from stochastic bandit feedback[END_REF][START_REF] Zuluaga | e-pal: An active learning approach to the multiobjective optimization problem[END_REF]. Going back to our motivation, in an adaptive trial that aims at pre-selecting a certain number of treatments or vaccine strategies for further investigations in clinical trials, practical constraints (the cost and feasibility of the trials) impose a constraint on the maximal number of interesting arms that can be identified. This motivates the introduction of a new setting where the agent is asked to identify at most k Pareto optimal arms. Interestingly the sampling rule that we propose for this setting can be used to solve (some generalizations of) other relaxations considered in the literature.

Related work

The work most closely related to ours is that of Auer et al. [START_REF] Auer | Pareto front identification from stochastic bandit feedback[END_REF], who propose a relaxation, which we refer to as ε 1 -PSI: their algorithm returns a set S that contains w.h.p. all the Pareto optimal arms and possibly some sub-optimal arms, which when increased by ε 1 coordinatewise become Pareto optimal. For arms that have sub-Gaussian marginals, they provide an instancedependent sample complexity bound scaling with some notion of sub-optimality gap for each arm. The work of Zaluaga et al. [START_REF] Zuluaga | Active learning for multi-objective optimization[END_REF][START_REF] Zuluaga | e-pal: An active learning approach to the multiobjective optimization problem[END_REF] studies a structured variant of fixed-confidence PSI in which the means are regular functions of arms' descriptors. They use Gaussian process modeling and obtain worse-case sample complexity bounds. In particular [START_REF] Zuluaga | e-pal: An active learning approach to the multiobjective optimization problem[END_REF] considers the identification of an ε-cover of the Pareto set, which is a representative subset of the (ε)-Pareto set that will be related to our (ε 1 , ε 2 )-PSI criterion. The algorithms of [START_REF] Auer | Pareto front identification from stochastic bandit feedback[END_REF] and those of [START_REF] Zuluaga | Active learning for multi-objective optimization[END_REF][START_REF] Zuluaga | e-pal: An active learning approach to the multiobjective optimization problem[END_REF] in the unstructured setting 1 have the same flavor: they sample uniformly from a set of active arms and remove arms that have been found sub-optimal (or not representative). Auer et al. [START_REF] Auer | Pareto front identification from stochastic bandit feedback[END_REF] further adds an acceptation mechanism to stop sampling some of the arms that have been found (nearly-)optimal and are guaranteed not to dominate an arm of the active set. In this paper, we propose instead a more adaptive exploration strategy, which departs from such accept/reject mechanisms and is suited for different types of relaxation, including our novel k-relaxation.

Adaptive Pareto Exploration (APE) leverages confidence intervals on the differences of arms' coordinates in order to identify a single arm to explore, in the spirit of the LUCB [START_REF] Kalyanakrishnan | PAC subset selection in stochastic multiarmed bandits[END_REF] or UGapEc [START_REF] Gabillon | Best arm identification: A unified approach to fixed budget and fixed confidence[END_REF] algorithms for Top-m identification in (one-dimensional) bandit models. These algorithms have been found out to be preferable in practice to their competitors based on uniform sampling and eliminations [START_REF] Kaufmann | Information complexity in bandit subset selection[END_REF], an observation that will carry over to APE. Besides the multi-dimensional observations, we emphasize that a major challenge of the PSI problem with respect to e.g. Top m identification is that the number of arms to identify is not known in advance. Moreover, when relaxations are considered, there are multiple correct answers. In the one-dimensional settings, finding optimal algorithms in the presence of multiple correct answers is notoriously hard as discussed by the authors of [START_REF] Degenne | Pure exploration with multiple correct answers[END_REF], and their lower-bound based approach becomes impractical in our multi-dimensional setting. Finally, we remark that the k-relaxation can be viewed as an extension of the problem of identifying any k-sized subset out of the best m arms in a standard bandit [START_REF] Chaudhuri | PAC identification of many good arms in stochastic multi-armed bandits[END_REF].

Beyond Pareto set identification, other interesting multi-objective bandit identification problems have been studied in the literature. For example [START_REF] Drugan | Designing multi-objective multi-armed bandits algorithms: A study[END_REF] propose an algorithm to identify some particular arms in the Pareto set through a scalarization technique [START_REF] Miettinen | Nonlinear multiobjective optimization[END_REF]. The idea is to turn the multi-objective pure-exploration problem into a single-objective one (unique optimal arm) by using a real-valued preference function which is only maximized by Pareto optimal arms (see e.g [START_REF] Miettinen | Nonlinear multiobjective optimization[END_REF] for some examples of these functions). In practice, a family of those functions can be used to identify many arms of the Pareto set but it is not always possible to identify the entire Pareto set using this technique (see e.g [START_REF] Ehrgott | Multicriteria Optimization[END_REF] for weighted sum with a family of weights vectors). In a different direction, the authors of [START_REF] Katz-Samuels | Feasible arm identification[END_REF] introduce the feasible arm identification problem, in which the goal is to identify the set of arms whose mean vectors belong to a known polyhedron P ⊂ R D . In a follow up work [START_REF] Katz-Samuels | Top feasible arm identification[END_REF], they propose a fixed-confidence algorithm for finding feasible arms that further maximize a given weighted sum of the objectives. In clinical trials, this could be used to find treatments maximizing efficacy (or a weighted sum of different efficacy indicators), under the constraint that the toxicity remains below a threshold. However, in the presence of multiple indicators of biological efficacy, choosing the weights may be difficult, and an examination of the Pareto set could be more suitable. Finally, some papers consider extensions of the Pareto optimality condition. The authors of [START_REF] Ararat | Vector optimization with stochastic bandit feedback[END_REF] tackle the identification of the set of non-dominated arms of any partial order defined by an R D polyhedral ordering cone (the usual Pareto dominance corresponds to using the cone defined by the positive orthant R D + ), and they provide worst-case sample complexity in the PAC setting. The work of [START_REF] Audiffren | Bandits dueling on partially ordered sets[END_REF] studies the identification of the set of non-dominated elements in a partially ordered set under the dueling bandit setting, in which the observations consists in pairwise comparison between arms.

Outline and contributions First, we formalize in Section 2 different relaxations of the PSI problem: ε 1 -PSI, as introduced by [START_REF] Auer | Pareto front identification from stochastic bandit feedback[END_REF], ε 1 , ε 2 -PSI, of which a particular case was studied by [START_REF] Zuluaga | e-pal: An active learning approach to the multiobjective optimization problem[END_REF] and ε 1 -PSI-k, a novel relaxation that takes as input an upper bound k on the maximum number of ε 1 -optimal arms that can be returned. Then, we introduce in Section 3 Adaptive Pareto Exploration, a simple, adaptive sampling rule which can simultaneously tackle all three relaxations, when coupled with an appropriate stopping rule that we define for each of them. In Section 4, we prove high-probability upper bounds on the sample complexity of APE under different stopping rules. For ε 1 -PSI, our bound slightly improves upon the state-of-the-art. Our strongest result is the bound for ε 1 -PSI-k, which leads to a new notion of sub-optimality gap, quantifying the reduction in sample complexity that is obtained. Then, Section 5 presents the result of a numerical study on synthetic datasets, one of them being inspired by a Covid-19 vaccine clinical trial. It showcases the good empirical performance of APE compared to existing algorithms, and illustrates the impact of the different relaxations.

Problem Setting

In this section, we introduce the Pareto Set Identification (PSI) problem and its relaxations. Fix K, D ∈ N . Let ν 1 , . . . , ν K be distributions over R D with means µ 1 , . . . , µ K ∈ R D . Let A := [K] := {1, . . . , K} denote the set of arms. Let ν := (ν 1 , . . . , ν K ) and X := (µ 1 , . . . , µ K ). We use boldfaced symbols for R D elements. Letting X ∈ R D , u ∈ R, for any d ∈ {1, . . . , D}, X d denotes the d-th coordinate of X and X + u := (X 1 + u, . . . , X D + u). In the sequel, we will assume that ν 1 , . . . , ν K have 1-subgaussian marginals2 . Definition 1. Given two arms i, j ∈ A, i is weakly (Pareto) dominated by j (denoted by

µ i ≤ µ j ) if for any d ∈ {1, . . . , D}, µ d i ≤ µ d j .
The arm i is (Pareto) dominated by j (µ i µ j or i j) if i is weakly dominated by j and there exists d ∈ {1, . . . , D} such that

µ d i < µ d j . The arm i is strictly (Pareto) dominated by j (µ i ≺ µ j or i ≺ j) if for any d ∈ {1, . . . , D}, µ d i < µ d j .
For ε ∈ R D + , the ε-Pareto set S ε (X ) is the set of ε-Pareto optimal arms, that is:

S ε (X ) := {i ∈ A s.t j ∈ A : µ i + ε ≺ µ j }.
In particular, S 0 (X ) is called the Pareto set and we will simply write S (X ) to denote S 0 (X ). When it is clear from the context, we write S (or S ε ) to denote S (X ) (or S ε (X )). By abuse of notation we write S ε when ε ∈ R + to denote S ε , with ε := (ε, . . . , ε).

In each round t = 1, 2, . . . , the agent chooses an arm A t and observes an independent draw X t ∼ ν At with E(X At ) = µ At . We denote by P ν the law of the stochastic process (X t ) t≥1 and by E ν , the expectation under P ν . Let F t := σ(A 1 , X 1 , . . . , A t , X t ) the σ-algebra representing the history of the process. An algorithm for PSI consists in : i) a sampling rule which determines which arm to sample at time t based on history up to time t -1, ii) a stopping rule τ which is a stopping time w.r.t the filtration (F t ) t≥1 and iii) a recommendation rule which is a F τ -measurable random set S τ representing the guess of the learner. The goal of the learner is to make a correct guess with high probability, using as few samples τ as possible. Before formalizing this, we introduce the different notion of correctness considered in this work, depending on parameters

ε 1 ≥, ε 2 ≥ 0 and k ∈ [K].
Our first criterion is the one considered by [START_REF] Auer | Pareto front identification from stochastic bandit feedback[END_REF].

Definition 2. S ⊂ A is correct for ε 1 -PSI if S ⊂ S ⊂ S ε1 .
To introduce our second criterion, we need the following definition.

Definition 3. Let ε 1 , ε 2 ≥ 0. A subset S ⊂ A is an (ε 1 , ε 2 )-cover of the Pareto set if : S ⊂ S ε1 and for any i / ∈ S either i / ∈ S or ∃j ∈ S such that µ i ≺ µ j + ε 2 .
The ε-accurate set of [START_REF] Zuluaga | e-pal: An active learning approach to the multiobjective optimization problem[END_REF] is a particular case of (ε 1 , ε 2 )-cover for which ε 1 = ε 2 = ε. Allowing ε 1 = ε 2 generalizes the notion of ε-correct set and can be useful, e.g., in scenarios when we want to identify the exact Pareto set (setting ε 1 = 0) but allow some optimal arms to be discarded if they are too close (parameterized by ε 2 ) to another optimal arm already returned. We note however that the sparse cover of [START_REF] Auer | Pareto front identification from stochastic bandit feedback[END_REF] is an (ε, ε)-cover with and additional condition that the arms in the returned set should not be too close to each over. Identifying a sparse cover from samples requires in particular to identify S ε1 hence it can not be seen as a relaxation of ε 1 -PSI.

Definition 4. S ⊂ A is correct for (ε 1 , ε 2 )-PSI if it is an (ε 1 , ε 2 )-cover of the Pareto set. Definition 5. S ⊂ A is correct for ε 1 -PSI-k if either i) | S|= k and S ⊂ S ε1 or ii) | S|< k and S ⊂ S ⊂ S ε1 holds.
Given a specified objective (ε 1 -PSI, (ε 1 , ε 2 )-PSI or ε 1 -PSI-k), and a target risk parameter δ ∈ (0, 1), the goal of the agent is to build a δ-correct algorithm, that is to guarantee that with probability larger than 1 -δ, her guess S τ is correct for the given objective, while minimizing the number of samples τ needed to make the guess, called the sample complexity.

We now introduce two important quantities to characterize the (Pareto) optimality or sub-optimality of the arms. For any two arms i, j, we let

m(i, j) := min 1≤d≤D (µ d j -µ d i ), and
M(i, j) := max 1≤d≤D (µ d i -µ d j ),
which have the following interpretation. If i j, m(i, j) is the minimal quantity α ≥ 0 that should be added component-wise to µ i so that µ i + α ⊀ µ j , α := (α, . . . , α). Moreover, m(i, j) > 0 if and only if i ≺ j. Then, for any arms i, j, if i ⊀ j, M(i, j) is the minimum quantity α such µ i ≤ µ j + α , α := (α , . . . , α ). We remark that M(i, j) < 0 if and only if i ≺ j. Our algorithms, presented in the next section, rely on confidence intervals on these quantities.

Adaptive Pareto Exploration

We describe in this section our sampling rule, Adaptive Pareto Exploration, and present three stopping and recommendation rules to which it can be combined to solve each of the proposed relaxation. Let T k (t) := t-1 s=1 1(A s = k) be the number of times arm k has been pulled up to round t and

µ k (t) := T k (t) -1 T k (t)
s=1 X k,s the empirical mean of this arm at time t, where X k,s denotes the s-th observation drawn i.i.d from ν k . For any arms i, j ∈ A, we let m(i, j, t) := min

d ( µ d j (t) -µ d i (t)) and M(i, j, t) := max d ( µ d i (t) -µ d j (t)).
The empirical Pareto set is defined as

S(t) := {i ∈ A : j ∈ A : µ i (t) ≺ µ j (t)}, = {i ∈ A : ∀j ∈ A\{i}, M(i, j, t) > 0} .

Generic algorithm(s)

Adaptive Pareto Exploration relies on a lower/upper confidence bound approach, similar to singleobjective BAI algorithms like UGapEc [START_REF] Gabillon | Best arm identification: A unified approach to fixed budget and fixed confidence[END_REF], LUCB [START_REF] Kalyanakrishnan | PAC subset selection in stochastic multiarmed bandits[END_REF] and LUCB++ [START_REF] Simchowitz | The simulator: Understanding adaptive sampling in the moderate-confidence regime[END_REF]. These three algorithms identify in each round two contentious arms: b t : a current guess for the optimal arm (defined as the empirical best arm or the arm with the smallest upper bound on its sub-optimality gap), c t : a contender of this arm; the arm which is the most likely to outperform b t (in all three algorithms, it is the arm with the largest upper confidence bound in [K]\{b t }). Then, either both arms are pulled (LUCB, LUCB++) or the least explored among b t and c t is pulled (UGapEc). The originality of our sampling rule lies in how to appropriately define b t and c t for the multi-objective setting. To define those, we suppose that there exists confidence intervals [L d i,j (t, δ), U d i,j (t, δ)] on the difference of expected values for each pair of arms (i, j) and each objective d ∈ D, such that introducing

E t := K i=1 j =i D d=1 L d i,j (t, δ) ≤ µ d i -µ d j ≤ U d i,j (t, δ) and E = ∞ t=1 E t , (1)
we have P(E) ≥ 1 -δ. Concrete choices of these confidence intervals will be discussed in Section 3.2.

To ease the notation, we drop the dependency in δ in the confidence intervals and further define

M -(i, j, t) := max d L d i,j (t) and M + (i, j, t) := max d U d i,j (t) (2)
m -(i, j, t) := -M + (i, j, t) and m + (i, j, t) := -M -(i, j, t).

(

) Lemma 1. For any round t ≥ 1, if E t holds, then for any i, j ∈ A, M -(i, j, t) ≤ M(i, j) ≤ M + (i, j, t) and m -(i, j, t) ≤ m(i, j) ≤ m + (i, j, t). 3
Noting that S ε1 = {i ∈ A : ∀j = i, M(i, j) + ε 1 > 0}, we define the following set of arms that are likely to be ε 1 -Pareto optimal:

OPT ε1 (t) := {i ∈ A : ∀j ∈ A\{i}, M -(i, j, t) + ε 1 > 0}.

Sampling rule In round t, Adaptive Pareto Exploration samples a t , the least pulled arm among two candidate arms b t and c t given by b t := argmax

i∈A\OPT ε 1 (t) min j =i M + (i, j, t), c t := argmin j =bt M -(b t , j, t)
The intuition for their definition is the following. Letting i be a fixed arm, note that M(i, j) > 0 for some j, if and only if there exists a component d such that µ d i > µ d j i.e i is not dominated by j. Moreover, the larger M(i, j), the more i is non-dominated by j in the sense that there exists d such that µ d i µ d j . Therefore, i is strictly optimal if and only if for all j = i, M(i, j) > 0 i.e α i := min j =i M(i, j) > 0. And the larger α i , the more i looks optimal in the sense that for each arm j = i, there exists a component d j for which i is way better than j. As the α i are unknown, we define b t as the maximizer of an optimistic estimate of the α i 's. We further restrict the maximization to arms for which we are not already convinced that they are optimal (by Lemma 1, the arms in OPT ε1 (t) are (nearly) Pareto optimal on the event E). Then, we note that for a fixed arm i, M(i, j) < 0 if and only if i is strictly dominated by j. And the smaller M(i, j), the more j is close to dominate i (or largely dominates it): for any component d, µ d i -µ d j is small (or negative). Thus, for a fixed arm i, argmin j =i M(i, j) can be seen as the arm which is the closest to dominate i (or which dominates it by the largest margin). By minimizing a lower confidence bound on the unknown quantity M(b t , j), our contender c t can be interpreted as the arm which is the most likely to be (close to) dominating b t . Gathering information on both b t and c t can be useful to check whether b t can indeed be optimal.

Interestingly, we show in Appendix E that for D = 1, our sampling rule is close but not identical to the sampling rules used by existing confidence-based best arm identification algorithms.

Stopping condition

Recommendation

Objective

τ ε1 Z ε1 1 (t) > 0 ∧ Z ε1 2 (t) > 0 O(τ ε1 ) ε 1 -PSI τ ε1,ε2 Z ε1,ε2 1 (t) > 0 ∧ Z ε1,ε2 2 (t) > 0 OPT ε1 (τ ε1,ε2 ) (ε 1 , ε 2 )-PSI τ k |OPT ε1 (t)|≥ k OPT ε1 (τ k ) ε 1 -PSI-k
O(t) := S(t) ∪ {i ∈ S(t) c : j = i : m -(i, j, t) > 0}
consists of the current empirical Pareto set plus some additional arms that have not yet been formally identified as sub-optimal. Those arms should be (ε 1 )-Pareto optimal.

For (ε 1 , ε 2 )-PSI we define a similar stopping rule τ ε1,ε2 where the stopping statistics are respectively replaced with

Z ε1,ε2 1 (t) := min i∈S(t) max(g ε2 i (t), h ε1 i (t)) and Z ε1,ε2 2 (t) := min i∈S(t) c max(g ε2 i (t), h ε1 i (t))
with the convention min ∅ = +∞, and the recommendation is OPT ε1 (τ ε1,ε2 ).

To tackle the ε 1 -PSI-k relaxation, we propose to couple τ ε1 with an additional stopping condition checking whether OPT ε1 (t) already contains k arms. That is, we stop at τ k ε1 := min τ ε1 , τ k where τ k := inf{t ≥ K : | OPT ε1 (t)| ≥ k} with associated recommendation OPT ε1 (τ k ). Depending of the reason for stopping (τ ε1 or τ k ), we follow the corresponding recommendation. Lemma 2. Assume E holds. For ε 1 -PSI (resp. (ε 1 , ε 2 )-PSI , ε 1 -PSI-k), Adaptive Pareto Exploration combined with the stopping rule τ ε1 (resp. τ ε1,ε2 , resp. τ k ε1 ) outputs a correct subset. Remark 1. We decoupled the presentation of the sampling rule to that of the "sequential testing" aspect (stopping and recommendation). We could even go further and observe that multiple tests could actually be run in parallel, for free. If we collect samples with APE (which only depends on ε 1 ), whenever one of the three stopping conditions given in Table 1 triggers, for any values of ε 2 or k, we can decide to stop and make the corresponding recommendation or continue and wait for another "more interesting" stopping condition to be satisfied. If E holds, a recommendation made at any such time will be correct for the objective associated to the stopping criterion (third column in Table 1).

Our instantiation

We propose to instantiate the algorithms with confidence interval on the difference of pair of arms. For any pair i, j ∈ A, we define a function β i,j such that for any d ∈

[D], U d i,j (t) = µ d i (t)-µ d j (t)+β i,j (t) and L d i,j (t) = µ d i (t) -µ d j (t) -β i,j (t)
. We take from [START_REF] Kaufmann | Mixture martingales revisited with applications to sequential tests and confidence intervals[END_REF] the following confidence bonus for timeuniform concentration:

β i,j (t) := 2   C g log K1 δ 2 + a∈{i,j} log(4 + log(T a (t)))     a∈{i,j} 1 T a (t)   , (4)
where K 1 := K(K -1)D/2 and C g ≈ x + log(x) is a calibration function. They result in the simple expressions M ± (i, j, t) = M(i, j, t) ± β i,j (t) and m ± (i, j, t) = m(i, j, t) ± β i,j (t). As an

Algorithm 1: ε 1 -APE-k: Adaptive Pareto Exploration for ε 1 -PSI-k Data: parameter ε 1 ≥ 0, k ∈ [K] initialize : sample each arm once, set t = K, T i (K) = 1 for any i ∈ A for t = K + 1, . . . , do S(t) = {i ∈ A : ∀j ∈ A\{i}, M(i, j, t) > 0}; OPT ε1 (t) = {i ∈ A : ∀j ∈ A\{i}, M(i, j, t) -β i,j (t) + ε 1 > 0}; if |OPT ε1 (t)|≥ k then break and output OPT ε1 (t) if Z ε1 1 (t) > 0 ∧ Z ε1 2 (t) > 0 then break and output O(t) = S(t) {i ∈ S(t) c : j = i : m(i, j, t) -β i,j (t) > 0} b t := argmax i∈A\OPT ε 1 (t) min j =i [M(i, j, t) + β i,j (t)]; c t := argmin i =bt [M(b t , j, t) -β bt,j (t)]; sample a t := argmin i∈{bt,ct} T i (t);
example, we state in Algorithm 1 the pseudo-code of APE combined the stopping rule suited for the k-relaxation of ε 1 -PSI, which we refer to as ε 1 -APE-k.

In Appendix F, we also study a different instantiation based on confidence bounds of the form

U i,j (t) = U i (t) -L j (t) where [L i (t), U i (t)
] is a confidence interval on µ i . This is the approach followed by LUCB for D = 1 and prior work on Pareto identification [START_REF] Auer | Pareto front identification from stochastic bandit feedback[END_REF][START_REF] Zuluaga | e-pal: An active learning approach to the multiobjective optimization problem[END_REF]. In practice we advocate the use of the pairwise confidence intervals defined above, even if our current analysis does not allow to quantify their improvement. For the LUCB-like instantiation, we further derive in Appendix F an upper bound on the expected stopping time of APE for the different stopping rules.

Theoretical analysis

In this section, we state our main theorem on the sample complexity of our algorithms and give a sketch of its proof.

First let us introduce some quantities that are needed to state the theorem. The sample complexity of the algorithm proposed by [START_REF] Auer | Pareto front identification from stochastic bandit feedback[END_REF] for (ε 1 )-Pareto set identification scales as a sum over the arms i of 1/(∆ i ∨ ε 1 ) 2 where ∆ i is called the sub-optimality gap of arm i and is defined as follows. For a sub-optimal arm i / ∈ S (X ),

∆ i := max j∈S m(i, j),
which is the smallest quantity that should be added component-wise to µ i to make i appear Pareto optimal w.r.t {µ i : i ∈ A}. For a Pareto optimal arm i ∈ S (X ), the definition is more involved:

∆ i := min j∈A\{i} ∆ j if S := {i} min(δ + i , δ - i ) else, where δ + i := min j∈S \{i} min(M(i, j), M(j, i)) and δ - i := min j∈A\S {(M(j, i)) + + ∆ j }.
For x ∈ R, (x) + := max(x, 0). We also introduce some additional notion needed to express the contribution of the k-relaxation. Let 1 ≤ k ≤ K. For any arm i, let ω i = min j =i M(i, j) and define

ω k := k max i∈A ω i , S ,k := 1...k argmax i∈A ω i ,
with the k-th max and first to k-th argmax operators. Observe that ω k > 0 if and only if |S (X )| ≥ k. Theorem 1. Fix a risk parameter δ ∈ (0, 1), ε 1 ≥ 0, let k ≤ K and ν a bandit with 1-subgaussian marginals. With probability at least 1 -δ, ε 1 -APE-k recommends a correct set for the ε 1 -PSI-k objective and stops after at most

a∈A 88 ∆ 2 a log 2K(K -1)D δ log 12e ∆ a ,
samples, where for each a ∈ A, ∆ a := max(∆ a , ε 1 , ω k ).

First, when k = K, observing that ε 1 -APE-K provides a δ-correct algorithm for ε 1 -PSI, our bound improves the result of [START_REF] Auer | Pareto front identification from stochastic bandit feedback[END_REF] for the ε 1 -PSI problem in terms of constant multiplicative factors and log log ∆ -1 terms instead of log ∆ -2 . It nearly matches the lower bound of Auer et al. [START_REF] Auer | Pareto front identification from stochastic bandit feedback[END_REF] for the ε 1 -PSI problem (Theorem 17 therein). It also shows the impact of the k-relaxation on the sample complexity. In particular, we can remark that for any arm i ∈ S \S ,k , max(∆ i , ω k ) = ω k . Intuitively, it says that we shouldn't pay more than the cost of identifying the k-th optimal arm, ordered by the ω i s. A similar result has been obtained for the any k-sized subset of the best m problem [START_REF] Chaudhuri | PAC identification of many good arms in stochastic multi-armed bandits[END_REF]. But the authors have shown the relaxation only for the best m arms while our result shows that even the sub-optimal arms should be sampled less.

In Appendix D, we prove the following lower bound showing that in some scenarios, ε 1 -APE-k is optimal for ε 1 -PSI-k , up to D log(K) and constant multiplicative terms. We note that for ε 1 -PSI a lower bound featuring the gaps ∆ a and ε 1 was already derived by Auer et al. [START_REF] Auer | Pareto front identification from stochastic bandit feedback[END_REF]. Theorem 2. There exists a bandit instance ν with

|S |= p ≥ 3 such that for k ∈ {p, p -1, p -2} any δ-correct algorithm for 0-PSI-k verifies E ν (τ δ ) ≥ 1 D log 1 δ K a=1 1 (∆ k a ) 2 ,
where ∆ k a := ∆ a + ω k and τ δ is the stopping time of the algorithm.

In Appendix C, we prove that Theorem 1 without the ω k terms also holds for (ε 1 , ε 2 )-APE. This does not justifies the reduction in sample complexity when setting ε 2 > 0 in (ε 1 , ε 2 )-PSI observed in our experiments but it at least guarantees that the ε 2 -relaxation doesn't make things worse.

Furthermore, since our algorithm allows ε 1 = 0, it is also an algorithm for BAI when D = 1, ε 1 = 0. We prove in Appendix E that in this case, the gaps ∆ i s matches the classical gaps in BAI [START_REF] Audibert | Best arm identification in multi-armed bandits[END_REF][START_REF] Kaufmann | On the complexity of best-arm identification in multi-armed bandit models[END_REF] and we derive its sample complexity from Theorem 1 showing that it is similar in theory to UGap [START_REF] Gabillon | Best arm identification: A unified approach to fixed budget and fixed confidence[END_REF], LUCB [START_REF] Kalyanakrishnan | PAC subset selection in stochastic multiarmed bandits[END_REF] and LUCB++ [START_REF] Simchowitz | The simulator: Understanding adaptive sampling in the moderate-confidence regime[END_REF] but have better empirical performance.

Sketch of proof of Theorem 1

Using Proposition 24 of [START_REF] Kaufmann | Mixture martingales revisited with applications to sequential tests and confidence intervals[END_REF] we first prove that the choice of β i,j in (4) yields P(E) ≥ 1 -δ for the good event E defined in [START_REF] Ararat | Vector optimization with stochastic bandit feedback[END_REF]. Combining this result with Lemma 2, yields that ε 1 -APE-k is correct with probability at least 1 -δ.

The idea of the remaining proof is to show that under the event E, if APE has not stopped at the end of round t, then the selected arm a t has not been explored enough. The first lemma showing this is specific to the stopping rule τ k ε1 used for

ε 1 -PSI-k. Lemma 3. Let ε 1 ≥ and k ∈ [K]. If E t holds and t < τ k ε1 then ω k ≤ 2β at,at (t).
The next two lemmas are more general as they apply to different stopping rules.

Lemma 4. Let ε 1 ≥ 0. Let τ = τ k ε1 for some k ∈ [K] or τ = τ ε1,ε2
for some ε 2 ≥ 0. If E t holds and t < τ then ∆ at ≤ 2β at,at (t).

Lemma 5. Let ε 1 ≥ 0 and τ be as in Lemma 4. If E t holds and t < τ then ε 1 ≤ 2β at,at (t).

As can be seen in Appendix B, the proofs of these three lemmas heavily rely on the specific definition of b t and c t . In particular, to prove Lemma 4 and 5, we first establish that when t < τ , any arm j ∈ A satisfies m(b t , j, t) ≤ β bt,j (t). The most sophisticated proof is then that of Lemma 4, which relies on a case distinction based on whether b t or c t belongs to the set of optimal arms. These lemmas permit to show that, on

E t if t < τ k ε1 then ∆ at < 2β at,at (t) ≤ 2β Ta t (t)
, where we define β n to be the expression of β i,j (t) when T i (t) = T j (t) = n. Then we have

τ k ε1 1{E} ≤ ∞ t=1 a∈A 1 {a t = a} ∧ ∆ a ≤ 2β Ta(t) a ≤ a∈A inf n ≥ 2 : ∆ a > 2β n . (5)
A careful upper-bounding of the RHS of ( 5) completes the proof of Theorem 1.

Experiments

We evaluate the performance of Adaptive Pareto Exploration on a real-world scenario and on synthetic random Bernoulli instances. For a fair comparison, Algorithm 1 of [START_REF] Auer | Pareto front identification from stochastic bandit feedback[END_REF], referred to as PSI-Unif-Elim and APE are both run with our confidence bonuses β i,j (t) on pairs of arms, which considerably improve single-arm confidence bonuses 3 . As anytime confidence bounds are known to be conservative, we use K 1 = 1 in (4) instead of its theoretical value coming from a union bound. Still, in all our experiments, the empirical error probability was (significantly) smaller than the target δ = 0.1.

Real-world dataset COV-BOOST [START_REF] Munro | Safety and immunogenicity of seven COVID-19 vaccines as a third dose (booster) following two doses of ChAdOx1 nCov-19 or BNT162b2 in the UK (COV-BOOST): a blinded, multicentre, randomised, controlled, phase 2 trial[END_REF] is phase 2 trial which was conducted on 2883 participants to measure the immunogenicity of different Covid-19 vaccines as third dose (booster) in various combinations of initially received vaccines (first two doses). This resulted in a total of 20 vaccination strategies being assessed, each of them defined by the vaccines received as first, second and third dose. The authors have reported the average responses induced by each candidate strategy on cohorts of participants, measuring several immunogenicity markers. From this study, we extract and process the average response of each strategy to 3 specific immunogenicity indicators: two markers of antibody response and one of the cellular response. The outcomes are assumed to have a log-normal distribution [START_REF] Munro | Safety and immunogenicity of seven COVID-19 vaccines as a third dose (booster) following two doses of ChAdOx1 nCov-19 or BNT162b2 in the UK (COV-BOOST): a blinded, multicentre, randomised, controlled, phase 2 trial[END_REF]. We use the average (log) outcomes and their variances to simulate a multivariate Gaussian bandit with K = 20, D = 3. We give in Appendix H.2 some additional details about the processing of the data, and report the means and variances of each arm. In Appendix H.1 we further explain how APE can be simply adapted when the marginals distributions of the arms have different variances. In this experiment, we set ε 1 = 0, δ = 0.1 and compare PSI-Unif-Elim to 0-APE-k (called As |S | = 2, we first observe that, without the relaxation (i.e. for k > 3), APE outperforms its state-of-the-art competitor PSI-Unif-Elim . Moreover for k = 1 or k = 2, the sample complexity of APE-k is significantly reduced. For k = 2 when the stopping time τ k is reached some sub-optimal arms have possibly not yet been identified as such, while for k = 3, even if the optimal arms have been identified, the remaining arms have to be sampled enough to ensure that they are sub-optimal before stopping. This explains the gap in sample complexity between k = 2 and k = 3. In Appendix H.3, we compare APE to an adaptation of PSI-Unif-Elim for the k-relaxation, showing that APE is always preferable.

APE[k=1] APE[k=2*] APE[k=3] APE[
Experiments on random instances To further illustrate the impact of the k-relaxation and to strengthen the comparison with PSI-Unif-Elim , we ran the previous algorithms on 2000 randomly generated multi-variate Bernoulli instances, with K = 5 arms and different values of the dimension D. We set δ = 0.1 and ε 1 = 0.005 (to have reasonable running time). The averaged sample complexities are reported in Table 2. We observe that APE (with k = K) uses 20 to 25% less samples than PSI-Unif-Elim and tends to be more efficient as the dimension increases (and likely the size of the Pareto set, since the instances are randomly generated). We also note that identifying a k-sized subset of the Pareto set requires considerably less samples than exact PSI. In Appendix H.3 we also provide examples of instances for which APE takes up to 3 times less samples than PSI-Unif-Elim . To illustrate the impact of the ε 2 relaxation, setting ε 1 = 0 we report the sample complexity of APE associated with the stopping time τ 0,ε2 for 20 equally spaced values of ε 2 ∈ [0.01, 0.05], averaged over 2000 random Bernoulli instances. Figure 2 shows the decrease of the average sample complexity when ε 2 increases (left) and the average ratio of the size of the returned set to the size of the Pareto set (right). Note that for ε 1 = 0, we have O(τ 0,ε2 ) ⊂ S . The average sample complexity reported decreases up to 86% for the instance with K = 5, D = 2 while the returned set contains more than 90% of the Pareto optimal arms. In Appendix H.3, we further illustrate the behavior of APE with the ε 2 relaxation on a fixed instance in dimension 2.

ε 1 -APE-1 ε 1 -APE-2 ε 1 -APE-3 ε 1 -APE-4 ε 1 -APE-5 ε 1 -PSI-Unif-Elim D =

Conclusion and perspective

We proposed and analyzed APE, an adaptive sampling rule in multi-variate bandits models that when coupled with different stopping rules can tackle different relaxations of the fixed-confidence Pareto Set Identification problem. Our experiments revealed the good performance of the resulting algorithms compared to the state-of-the-art PSI algorithm as well as the great reductions in sample complexity brought by the relaxations. In future work, we intend to make our algorithms more practical for possible applications to clinical trials. For this purpose, as measuring efficacy takes time, we will investigate its adaptation to a batch setting, following, e.g. the work of [START_REF] Jun | Top arm identification in multi-armed bandits with batch arm pulls[END_REF] for BAI. We will also investigate the use of APE beyond the fixed-confidence setting, to the possibly more realistic fixed-budget [START_REF] Audibert | Best arm identification in multi-armed bandits[END_REF] or anytime exploration [START_REF] Jun | Anytime exploration for multi-armed bandits using confidence information[END_REF] settings. To the best of our knowledge, no algorithm exists in the literature for PSI in such settings. Finally, following the works of [START_REF] Zuluaga | Active learning for multi-objective optimization[END_REF][START_REF] Auer | Pareto front identification from stochastic bandit feedback[END_REF], we defined the ε 1 , ε 2 relaxations with scalar values, so that the same slack applies to all components. Although we could easily modify our algorithms to tackle vectorial values ε 1 , ε 2 , so far we could only prove a dependence on min d ε d 1 in the sample complexity. We intend to study the right quantification in the sample complexity when ε 1 and ε 2 are vectorial.

Outline and Notation

In this section, we provide an outline of the supplemental material and define some additional notation. Appendix A proves the correctness of ours algorithms and some concentration lemmas. In Appendix B, we prove Theorem 1 and the lemmas used in its proof. In Appendix C we analyze the correctness and sample complexity of APE associated to the stopping time τ ε1,ε2 . Appendix D describes our worst-case lower bound and in Appendix E we relate our algorithm to other algorithms for BAI. In Appendix F we derive an upper-bound on the expectation of the sample complexity of ε 1 -APE-k with a LUCB1-like instantiation and in Appendix G we recall or prove some technical lemmas that are used in the main proofs. Finally, in Appendix H we give further details about the experiments together with additional experimental results.

Notation

Type

Description A Set of arms K N Number of arms D N Dimension or number of attributes [n] {1, . . . , n} m A 2 → R m(i, j) := min{µ d j -µ d i : d ∈ [D]} M A 2 → R M(i, j) := max{µ d i -µ d j : d ∈ [D]} (x) + R → R + max(0, x) µ a R D Mean of arm a ∈ A S ε {i ∈ A : j = i ∈ A : µ i + ε ≺ µ j } Table 3: Table of notation.

A Correctness for different stopping rules

In this section, we gather and prove results that are related to the correctness of our algorithms, either for their generic form (Lemma 1 and Lemma 2) or some specific calibration. We recall the definition of the events

E t = K i=1 j =i D d=1 L d i,j (t) ≤ µ d i -µ d j ≤ U d i,j (t) and E = ∞ t=1 E t .
A.1 Proof of Lemma 1 Lemma 1. For any round t ≥ 1, if E t holds, then for any i, j ∈ A, M -(i, j, t) ≤ M(i, j) ≤ M + (i, j, t) and m -(i, j, t) ≤ m(i, j) ≤ m + (i, j, t).

Proof of Lemma 1. This result simply follows from the definition of E t . Since

E t := K i=1 j =i D d=1 L d i,j (t) ≤ µ d i -µ d j ≤ U d i,j (t) , if E t holds, then for any i, j M -(i, j, t) := max d L d i,j (t) ≤ M(i, j) := max d (µ d i -µ d j ) ≤ max d U d i,j (t)
:= M + (i, j, t), and the second point follows by noting that m(i, j) = -M(i, j) and m + (i, j, t) := -M -(i, j, t); m -(i, j, t) := -M + (i, j, t) for any pair of arms.

We remark that when the algorithm uses confidence bonus of form (

µ d i (t) -µ d j (t)) ± β i,j (t), M + (i, j, t) := max d U d i,j (t) = max d ( µ d i (t) -µ d j (t)) + β i,j (t) = M(i, j, t) + β i,j (t), M -(i, j, t) := max d L d i,j (t) = max d ( µ d i (t) -µ d j (t)) -β i,j (t) = M(i, j, t) -β i,j (t)
, and the previous lemma implies that on E t , |M(i, j) -M(i, j, t)|≤ β i,j (t) and |m(i, j) -m(i, j, t)|≤ β i,j (t), which is extensively used in our sample complexity analyses.

A.2 Proof of Lemma 2 Lemma 2. Assume E holds. For ε 1 -PSI (resp. (ε 1 , ε 2 )-PSI , ε 1 -PSI-k), Adaptive Pareto Exploration combined with the stopping rule τ ε1 (resp. τ ε1,ε2 , resp. τ k ε1 ) outputs a correct subset.

We show the correctness of ε 1 -PSI-k (for any k) and we derive the correctness for ε 1 -PSI which is equivalent to ε 1 -PSI-K. The correctness of (ε 1 , ε 2 )-PSI is shown separately in Lemma 11 (see Appendix C).

Proof of Lemma 2. Assume E holds. Let t = τ k ε1 and i ∈ OPT ε1 (t). Since i ∈ OPT ε1 (t), for any j = i,

M(i, j) + ε 1 E ≥ M -(i, j, t) + ε 1 > 0, that is i ∈ S ε1
. Therefore, on the event E, OPT ε1 (t) ⊂ S ε1 . Thus, if the stopping has occurred because |OPT ε1 (t)|≥ k, since in this case O(t) ⊂ OPT ε1 (t) ⊂ S ε1 , all the recommended arms will be (ε 1 )-Pareto optimal. On the contrary, if |OPT ε1 (t)|< k, then from the definition of τ k ε1 it holds that Z ε1 1 (t) > 0 and Z ε1 2 (t) > 0, and the recommended set is then

O(t) := S(t) ∪ {i ∈ S(t) c : j = i : m -(i, j, t) > 0}.
For any i ∈ O(t) c , by the definition of the recommended set and since Z 2 (t) > 0,

∃j ∈ A such that m(i, j) E ≥ m -(i, j, t) > 0, so i is a sub-optimal arm. Therefore, S ⊂ O(t). Moreover, for any i ∈ O(t) ∩ S(t), since Z ε1 1 (t) > 0 we have h ε1 i (t) > 0, that is min j∈A\{i} M(i, j) + ε 1 E ≥ min j∈A\{i} M -(i, j, t) + ε 1 > 0. (6)
If i ∈ O(t) ∩ S(t) c , by definition of O(t), we have g i (t) < 0. However, since Z ε1 2 (t) > 0, max(g i (t), h ε1 i (t)) > 0 so we also have h ε1 i (t) > 0 and (6) applies. Thus, for any i ∈ O(t), min j∈A\{i} M(i, j) + ε 1 > 0, that is i ∈ S ε1 , so S ⊂ O(t) ⊂ S ε1 . Finally we can conclude that ε 1 -APE-k and ε 1 -APE output a correct subset on E.

A.3 Calibration of the confidence intervals

In Section 3.2 we proposed to instantiate our algorithms with the confidence intervals

U d i,j (t) = µ d i (t) -µ d j (t) + β i,j (t) and L d i,j (t) = µ d i (t) -µ d j (t) -β i,j (t) . (7)
We prove below that E is indeed a high-probability event for a suitable choice of β i,j (t). Lemma 6. Let ν be a bandit with 1-subgaussian marginals. For the confidence intervals defined in [START_REF] Ehrgott | Multicriteria Optimization[END_REF], with

β i,j (t) = 2   C g log K1 δ 2 + a∈{i,j} log(4 + log(T a (t)))     a∈{i,j} 1 T a (t)   .
the event

E = ∞ t=1 E t with E t = K i=1 j =i D d=1 L d i,j (t) ≤ µ d i -µ d j ≤ U d i,j (t)
is such that P(E) ≥ 1 -δ.

Proof. By observing that for any pair of arm β i,j = β j,i , E t can be rewritten as

E t = {i,j}∈Γ D d=1 L d i,j (t, δ) ≤ µ d i -µ d j ≤ U d i,j (t, δ) , (8)
= {i,j}∈Γ D d=1 |( µ d i (t) -µ d j (t)) -(µ d i -µ d j )|≤ β i,j (t) , (9)
where Γ := 2 [K] is the set of pair of 2 elements of [K], which satisfies |Γ|= K(K -1)/2. Therefore, using a union bound,

P(E c ) = P(∃t ≥ 1 : E c t holds ), = P ∃t ≥ 1, d ∈ [D], {i, j} ∈ Γ : |( µ d i (t) -µ d j (t)) -(µ d i -µ d j )|> β i,j (t) , ≤ {i,j}∈Γ D d=1 P ∃t ≥ 1 : |( µ d i (t) -µ d j (t)) -(µ d i -µ d j )|> β i,j (t) , ≤ {i,j}∈Γ D d=1 δ K 1
(by Proposition 24 of [START_REF] Kaufmann | Mixture martingales revisited with applications to sequential tests and confidence intervals[END_REF] which we recall below in Lemma 7),

= δ, since K 1 := K(K-1)D 2
and |Γ|= K(K -1)/2. Lemma 7 (Proposition 24 of [START_REF] Kaufmann | Mixture martingales revisited with applications to sequential tests and confidence intervals[END_REF]). Let X, Y be centered 1-subgaussian random variables and δ ∈ (0, 1). Let X, X 1 , X 2 , . . . be i.i.d random variables and Y, Y 1 , Y 2 , . . . be i.i.d random variables. With probability at least 1 -δ, for all p, q ≥ 1,

1 p p s=1 X s - 1 q q s=1 Y s ≤ 2 C g log(1/δ) 2
+ log log(e 4 p) + log log(e 4 q)

1 p + 1 q
where C g (x) ≈ x + log(x).

B Sample complexity analysis

In this section we prove Theorem 1 which is restated below. Theorem 1. Fix a risk parameter δ ∈ (0, 1), ε 1 ≥ 0, let k ≤ K and ν a bandit with 1-subgaussian marginals. With probability at least 1 -δ, ε 1 -APE-k recommends a correct set for the ε 1 -PSI-k objective and stops after at most

a∈A 88 ∆ 2 a log 2K(K -1)D δ log 12e ∆ a ,
samples, where for each a ∈ A, ∆ a := max(∆ a , ε 1 , ω k ).

Proof of Theorem 1. The correctness follows from Lemma 2 and the fact that P(E) ≥ 1 -δ (Lemma 6). The upper-bound on the sample complexity is a direct consequence of Lemma 3, Lemma 4, Lemma 5 which are proved later in this section. Indeed, using these lemmas we have that, if ε 1 -APE -k has not stopped during round t i.e t < τ k ε1 and the event E t holds, then a)

ω k ≤ 2β at,at (t) (Lemma 3), b) ∆ at ≤ 2β at,at (t)(Lemma 4), c) ε 1 ≤ 2β at,at (t) (Lemma 5)
hold simultaneously. Then, if we do not count the first K rounds due to initialization, and letting

∆ a := max(ω k , ε 1 , ∆ a ), τ k ε1 1{E} -1 ≤ ∞ t=1 1{E}1{τ k ε1 > t}, ≤ ∞ t=1 1{max(ω k , ε 1 , ∆ at ) ≤ 2β at,at (t)} = ∞ t=1 1{ ∆ at ≤ 2β at,at (t)} = ∞ t=1 K a=1 1{{a t = a} ∧ { ∆ a ≤ 2β a,a (t)}} = K a=1 ∞ t=1 1{{a t = a} ∧ { ∆ a ≤ 2β a,a (t)}} ≤ K a=1 inf{n ≥ 2 : ∆ a > 2β n },
where β n is the expression of β i,j (t) when T i (t) = T j (t) = n, that is

β n = 2 C g log K1 δ 2 + 2 log(4 + log(n)) 2 n .
Then, an inversion result given in Lemma 19 yields

inf{s ≥ 2 : 2β s < ∆ a } ≤ 88 ∆ 2 a log 2K(K -1)D δ log 12e ∆ a .
Therefore,

τ k ε1 1{E} ≤ a∈A 88 ∆ 2 a log 2K(K -1)D δ log 12e
∆ a .

We will now prove the lemmas involved in the proof of the main theorem. Two of them (Lemma 4, Lemma 5) rely on the following result, which is an important consequence of the definition of the APE sampling rule.

Lemma 8. Let ε 1 ≥ 0, ε 2 ≥ 0 and k ∈ [K]. If τ = τ k ε1 or τ = τ ε1,ε2
the following holds. If t < τ then for any j ∈ A, m(b t , j, t) ≤ β bt,j (t).

Proof. The proof is split in two steps.

Step 1 If t < τ k ε1 then for any j ∈ A, m(b t , j, t) ≤ β bt,j (t). First, note that t < max(τ k ε1 , τ ε1 ) implies that Z ε1 1 (t) ≤ 0 or Z ε1 2 (t) ≤ 0. By definition of b t and noting that M(i, j, t) = -m(i, j, t), we have

b t ∈ argmin i∈OPT ε 1 (t) c max j =i m(i, j, t) -β i,j (t). ( 10
)
so that if there exists j such that m(b t , j, t) > β bt,j (t), then

max j =bt m(b t , j, t) -β bt,j (t) > 0, therefore, ∀i ∈ OPT ε1 (t) c , max j =i m(i, j, t) -β i,j (t) > 0 i.e g i (t) > 0. (11)
Furthermore, for any i ∈ OPT ε1 (t), h ε1 i (t) > 0. Putting things together, if there exists j such that m(b t , j, t) > β bt,j (t) then, Z ε1 1 (t) > 0 and Z ε1 2 (t) > 0.

Step 2 If t < τ ε1,ε2 then for any j ∈ A, m(b t , j, t) ≤ β bt,j (t).

Recall that by definition t < τ ε1,ε2 implies that Z ε1,ε2 1 (t) ≤ 0 or Z ε1,ε2 2 (t) ≤ 0. Using [START_REF] Garivier | Optimal best arm identification with fixed confidence[END_REF], if there exists j such that m(b t , j, t) > β bt,j (t), then max j =bt m(b t , j, t) -β bt,j (t) > 0.

Combining this with

g ε2 i (t) := max j∈A\{i} m(i, j, t) -β i,j (t) + ε 2 1{j ∈ OPT ε1 (t)}, yields ∀i ∈ OPT ε1 (t) c , 0 < max j =i m(i, j, t) -β i,j (t) ≤ g ε2 i (t). (12)
Furthermore, since we have ∀i ∈ OPT ε1 (t) , h ε1 i (t) > 0,

the initial assumption would yield that for any arm i, max(h ε1 i (t), g ε2 i (t)) > 0, so Z ε1,ε2

1 (t) > 0 and Z ε1,ε2 2 (t) > 0.
We conclude that if τ = τ k ε1 or τ ε1,ε2 , t < τ implies that for any j ∈ A, m(b t , j, t) ≤ β bt,j (t).

B.1 Proof of Lemma 3 Lemma 3. Let ε 1 ≥ and k ∈ [K]. If E t holds and t < τ k ε1 then ω k ≤ 2β at,at (t).
Proof of Lemma 3. First, note that if k > |S |, then the lemma holds trivially since ω k < 0. In the sequel, we assume E t holds and k

≤ |S |. If t < τ k ε1 then it holds that |OPT ε1 (t)|< k. So S ,k ∩ OPT ε1 (t) c = ∅. Let i ∈ S ,k ∩ OPT ε1 (t) c , we have ω k ≤ ω i = min j∈A\{i} M(i, j), (a)
≤ min j∈A\{i} M(i, j, t) + β i,j (t),

(b) ≤ min j∈A\{bt} M(b t , j, t) + β bt,j (t), ≤ M(b t , c t , t) + β bt,ct (t), (c)
≤ 2β bt,ct (t), ≤ 2β at,at (t),
where (a) uses that E t holds and Lemma 1, (b) uses the definition of b t and (c) follows from the definition of c t and the fact that b t / ∈ OPT ε1 (t),which yields M(b t , c t , t) ≤ β bt,ct (t). The last inequality follows since a t is the least sampled among b t , c t and β is decreasing.

B.2 Proof of Lemma 4 Lemma 4. Let ε 1 ≥ 0. Let τ = τ k ε1 for some k ∈ [K] or τ = τ ε1,ε2
for some ε 2 ≥ 0. If E t holds and t < τ then ∆ at ≤ 2β at,at (t).

Before proving the Lemma 4, we state the following lemma which is used to derive an upper bound on the gap of an optimal arm. Its proof is postponed to the end of the section. Lemma 9. For any Pareto optimal arm i, ∆ i ≤ min j =i M(i, j).

Proof of Lemma 4. Assume that E t holds. We consider four different cases depending on whether b t and c t are optimal or sub-optimal.

Case 1.1 b t is a Pareto optimal arm. From the definition of the gap of an optimal arm and using Lemma 9 it follows ∆ bt ≤ M(b t , c t ) which on E t and using Lemma 1 yields

∆ bt + ε 1 ≤ M(b t , c t , t) + β bt,ct (t) + ε 1 (14)
then, noting that there exists j ∈ A \ {b t } such that M(b t , j, t) + ε 1 ≤ β bt,j (t), by definition of c t , we have

M(b t , c t , t) + ε 1 ≤ β bt,ct (t), (15) therefore,
∆ bt + ε 1 ≤ 2β bt,ct (t).
Case 1.2 b t is a sub-optimal arm. By definition of c t and using M = -m, we have

c t ∈ argmax j∈A\{bt} m(b t , j, t) + β bt,j (t), (16)
then, from the definition of the gap of a sub-optimal arm and since E t holds, we know that there exists an arm b t such that

∆ bt = m(b t , b t ) ≤ m(b t , b t , t) + β bt,b t (t), (a)
≤ m(b t , c t , t) + β bt,ct (t),

(b) ≤ 2β bt,ct (t).
where (a) uses the definition of c t and (b) uses Lemma 8.

Case 2.1 c t is a Pareto optimal arm. If b t is also an optimal arm, it follows that

∆ ct ≤ M(b t , c t ) which on E t yields ∆ ct ≤ M(b t , c t , t) + β bt,ct (t)
≤ max(2β bt,ct (t), m(b t , c t , t) + β bt,ct (t))

(b) ≤ 2β bt,ct (t).
where (a) follows from (x -y) + + (x + y) ≤ max(x + y, 2x) and (b) follows from m(b t , c t , t) ≤ β bt,ct (t) (Lemma 8).

Case 2.2 c t is a sub-optimal arm. We know that there exists an arm c t such that ∆ ct = m(c t , c t ). If c t = b t then, since m(j, i) ≤ M(i, j) (follows from the definition), we have

∆ ct = m(c t , c t ) = m(c t , b t ), ≤ M(b t , c t ), (a)
≤ then recalling that β i,j = β j,i ,

µ d ct -µ d bt + 2β bt,ct (t) Et ≥ ( µ d ct (t) -µ d bt (t) -β bt,ct (t)) + 2β bt,ct (t) ≥ µ d c t -µ d bt .
Put together, there exists d ∈ [D] such that

µ d c t -µ d ct ≤ 2β bt,ct (t), so ∆ ct = min d (µ d c t -µ d ct ) ≤ 2β bt,ct (t),
Putting the four cases together, we have proved that if t < max(τ k ε1 , τ ε1,ε2 ) then both

∆ bt ≤ 2β bt,ct (t) and ∆ ct ≤ 2β bt,ct (t) (17)
holds. Further noting that a t is the least sampled among among b t , c t and β is non-increasing, β bt,ct (t) ≤ β at,at (t), (17) yields ∆ at ≤ 2β at,at (t), which achieves the proof.

B.3 Proof of Lemma 5

The following lemma holds for each of the stopping times τ ε1 , τ ε1,ε2 and τ k ε1 . Lemma 5. Let ε 1 ≥ 0 and τ be as in Lemma 4. If E t holds and t < τ then ε 1 ≤ 2β at,at (t).

Proof of Lemma 5. By Lemma 8, we have m(b t , c t , t) ≤ β bt,ct (t) or equivalently

M(b t , c t , t) ≥ -β bt,ct (t) . (18)
Then, knowing that b t / ∈ OPT ε1 (t), there exists an arm j such that ε 1 + M(b t , j, t) ≤ β bt,j (t). Using further the definition of c t , it follows that ε 1 + M (b t , c t , t) ≤ β bt,ct (t). Combining this with inequality [START_REF] Kaufmann | On the complexity of best-arm identification in multi-armed bandit models[END_REF] and noting that a t is the least sampled among b t , c t yields

β at,at (t) ≥ β bt,ct (t) ≥ ε 1 /2.

B.4 Auxiliary results

We state the following lemma which is used to prove Lemma 9. Lemma 10. For any sub-optimal arm a, there exists a Pareto optimal arm a such that µ a ≺ µ a and ∆ a = m(a, a ) > 0. Moreover, For any i

∈ A \ S , j ∈ S , i) max k∈S m(i, k) = max k∈A m(i, k),
ii) If i ∈ argmin k∈A\{j} M(j, k) then j is the unique arm such that µ i ≺ µ j Proof. Assume there are p < n dominated arms. Without loss of generality, we may assume they are µ 1 , . . . , µ p . Let i 1 ≤ p. Suppose that no Pareto-optimal arm dominates µ i1 . Since µ i1 is not optimal, by the latter assumption, there exists i 2 ≤ p such that µ i1 ≺ µ i2 . If µ i2 is dominated by a Pareto optimal arm, this arm also dominates µ i1 (strict dominance is transitive) which contradicts the initial assumption. If not, there exits i 3 ≤ p such that µ i1 ≺ µ i2 ≺ µ i3 . Again we can use the same reasoning as before for i 3 . In any case we should stop in at most p steps, otherwise we would have µ i1 ≺ µ i2 ≺ • • • ≺ µ ip and µ ip should be dominated by a Pareto-optimal arm, otherwise it would be itself Pareto-optimal, which is not the case. Therefore, for any a ∈ A \ S , there exists a ∈ S such that a ≺ a and ∆ a = m(a, a ) > 0.

Letting i be a sub-optimal arm, since for any a ∈ A \ S , there exists a ∈ S such that a ≺ a , it follows that ∀d ∈

[D], µ d a -µ d i < µ d a -µ d i , which leads to m(i, a) ≤ m(i, a ), so max j∈A m(i, j) = max j∈S m(i, j) > 0,
which achieves the proof of the first point i). For the second point, let q ∈ A \ S and q such that q ≺ q and q ∈ argmin a∈A\{j} M(j, a).

By direct algebra, since q ≺ q , we have

M(j, q ) < M(j, q),
which is impossible if q = j (because q belongs to the argmin). Therefore, if

q ∈ argmin a∈A\{j} M(j, a)
is a sub-optimal arm, then j is the only arm such that q ≺ j (i.e µ q ≺ µ j ).

We now prove Lemma 9 which follows from the previous lemma.

Proof of Lemma 9. If argmin j =i M(i, j) ⊂ S , then the lemma follows from the definition of the gap of an optimal arm recalled in Section 4. If min j =i M(i, j) = M(i, a), a / ∈ S , then, from Lemma 10, i is the unique arm which dominates a so ∆ a = m(a, i) and using the definition of the gap of an optimal arm,

∆ i ≤ M(a, i) + + ∆ a , = 0 + m(a, i) ≤ M(i, a),
where we have used the the fact that m(p, q) ≤ M(q, p) for any pair of arms p, q (which follows from the definition). Therefore, for an optimal arm i, we always have

∆ i ≤ min j =i M(i, j).
C Algorithm for finding an (ε 1 , ε 2 )-cover

In this section, we analyse the sample complexity of APE when it is associated to the stopping time τ ε1,ε2 for identifying an (ε 1 , ε 2 )-cover of the Pareto set. The sampling rule remains unchanged an we prove that the algorithm does not require more samples to find an (ε 1 , ε 2 )-cover than to solve the ε 1 -PSI problem.

Algorithm 2: (ε 1 , ε 2 )-APE Data: parameter ε 1 ≥ 0 initialize :sample each arm once, set t = K, T i (t) = 1 for any i ∈ A for t = K + 1, . . . , do b t := argmax i∈A\OPT ε 1 (t) min j =i M + (i, j, t); c t := argmin j =bt M -(b t , j, t); if Z ε1,ε2 1 (t) > 0 ∧ Z ε1,ε2 2
(t) > 0 then break and output OPT ε1 (t); sample a t := argmin i∈{bt,ct} T i (t);

We recall the stopping time τ ε1,ε2 .

Stopping rule

Let ε 1 , ε 2 ≥ 0 and 0 < δ < 1. Then τ ε1,ε2 := inf {t ≥ K : Z ε1,ε2 1 (t) > 0 ∧ Z ε1,ε2 2 (t) > 0} , (19)
where,

Z ε1,ε2 1 (t) := min i∈S(t) max(g ε2 i (t), h ε1 i (t)) Z ε1,ε2 2 (t) := min i∈S(t) c max(g ε2 i (t), h ε1 i (t)), and
g ε2 i (t) := max j∈A\{i} m -(i, j, t) + ε 2 1{j ∈ OPT ε1 (t)} h ε1 i (t) := min j∈A\{i} M -(i, j, t) + ε 1
Recommendation rule When it is associated to the stopping time τ ε1,ε2 , APE recommends

O(τ ε1,ε2 ) := OPT ε1 (τ ε1,ε2 ),
which can be understood as follows. When τ ε1,ε2 is reached, the arms that are not yet identified as (nearly) optimal are either ε 2 -dominated by an arm in OPT ε1 (τ ε1,ε2 ) or sub-optimal, which is proven formally in Lemma 11. Lemma 11. Fix δ ∈ (0, 1), ε 1 , ε 2 ≥ 0 then (ε 1 , ε 2 )-APE recommends an (ε 1 , ε 2 )-cover of the Pareto set on the event E.

Proof of Lemma 11. Assume E holds. Let t = τ ε1,ε2 and i ∈ OPT ε1 (t). Since i ∈ OPT ε1 (t),

for any j = i, M(i, j) + ε 1 E ≥ M -(i, j, t) + ε 1 > 0 that is i ∈ S ε1 . Therefore, on the event E, OPT ε1 (t) ⊂ S ε1 . When the stopping time τ ε1,ε2 δ is reached, Z ε1,ε2 1 (t) > 0 and Z ε1,ε2 2 (t) > 0. Under this condition, OPT ε1 (t) = ∅.
Indeed, since Z ε1,ε2 1 (t) > 0 and Z ε1,ε2 2 (t) > 0, if OPT ε1 (t) = ∅ then, by the stopping rule and since OPT ε1 (t) = ∅, for any arm i, we would have h ε1 i (t) < 0 and g ε2 i (t) > 0. That is, for any arm i ∈ A,

∃j = i such that m(i, j) E > m -(i, j, t) > 0,
so every arm would be strictly dominated, which is impossible since the Pareto set cannot be empty. Then, OPT ε1 (t) = ∅ and for any i ∈ O(t) c = OPT ε1 (t) c , by the stopping rule it holds that max(g ε2 i (t), h ε1 i (t)) > 0. Further noting that for such arm i ∈ OPT ε1 (t) c , h ε1 i (t) < 0 , we thus have

g ε2 i (t) > 0, that is m -(i, j, t) + ε 2 1{j ∈ O(t)} > 0,
which on the event E yields

m(i, j) + ε 2 1{j ∈ O(t)} > 0.
Therefore, for such arm i, either

i) ∃j ∈ A such that m(i, j) > 0 that is µ i ≺ µ j or ii) ∃j ∈ O(t) such that m(i, j) + ε 2 > 0 that is µ i ≺ µ j + ε 2 with ε 2 := (ε 2 , . . . , ε 2 ).
Put together, O(t) ⊂ S ε1 and for any i / ∈ O(t), either i / ∈ S (i is a sub-optimal arm) or there exists

j ∈ O(t) such that µ i ≺ µ j + ε 2 . Thus O(t) is an (ε 1 , ε 2 )-cover of the Pareto set and (ε 1 , ε 2 )-APE is correct for (ε 1 , ε 2 )-cover identification.
The two lemmas restated below are used to prove identically to Theorem 1, the main theorem of this section. Lemma 4. Let ε 1 ≥ 0. Let τ = τ k ε1 for some k ∈ [K] or τ = τ ε1,ε2 for some ε 2 ≥ 0. If E t holds and t < τ then ∆ at ≤ 2β at,at (t).

The following lemma holds for each of the stopping times τ ε1 , τ ε1,ε2 and τ k ε1 . Lemma 5. Let ε 1 ≥ 0 and τ be as in Lemma 4. If E t holds and t < τ then ε 1 ≤ 2β at,at (t). Theorem 3. Fix δ ∈ (0, 1), ε 1 , ε 2 ≥ 0. Then (ε 1 , ε 2 )-APE outputs an (ε 1 , ε 2 )-cover of the Pareto set with probability at least 1 -δ using at most

a∈A 88 (∆ ε a ) 2 log 2K(K -1)D δ log 12e ∆ ε a ( 20
)
samples, where for all a ∈ A, ∆ ε a := max(∆ a , ε 1 ). This is the first problem-dependent sample complexity upper-bound for the (ε 1 , ε 2 )-cover of the Pareto set. In particular, this result holds for the ε-accurate Pareto set identification [START_REF] Zuluaga | e-pal: An active learning approach to the multiobjective optimization problem[END_REF] which corresponds to the particular case ε 1 = ε 2 = ε of the Pareto set cover. Therefore, (ε, ε)-APE could be compared to ε-PAL for ε-accurate Pareto set, which however relies on a Gaussian process modeling assumption.

While this sample complexity result upper-bound does not clearly show the dependence in ε 2 , we note that for some problems, we have a nearly matching lower bound that does not depend on ε 2 . In particular, consider the case D = 1, ε 1 = 0, ε 2 > 0 and assume there is a unique best arm (classical assumption in BAI) a . For this setting, an algorithm for (ε 1 , ε 2 )-cover identification is required to output a set S such that S ⊂ S = {a } and for any i = a either µ i < µ a or µ i ≤ µ a + ε 2 which trivially holds as long as S ⊂ S . Therefore, this problem is equivalent to (exact) Best Arm Identification. Almost matched lower bounds for BAI are known and does not depend on ε 2 ([18, 26, 10]). This observation can be generalized to any configuration where there is a unique (Pareto) optimal arm. Letting D ≥ 1, ε 1 = 0, ε 2 > 0 and ν a bandit with one Pareto optimal arm a , any algorithm for (ε 1 , ε 2 )-covering is required to output a set S ⊂ S = {a }. And for any i = a either µ i ≺ µ a or µ i ≺ µ a + ε 2 which trivially holds as long as S ⊂ S = {a }. So, on theses instances, (0, ε 2 )-covering is equivalent 0-PSI and the nearly matched lower of [START_REF] Auer | Pareto front identification from stochastic bandit feedback[END_REF] for 0-PSI does not depend on ε 2 (Theorem 17 therein).

In our experiments (see subsection H.3), we will see that in configurations with multiple Pareto optimal arms, the parameter ε 2 can still help to empirically reduce the sample complexity. Quantifying its precise impact on the sample complexity is left as future work.

D Lower Bound

In this section, we give a gap-dependent lower-bound for the k-relaxation in some configurations. We use the change of distribution lemma of [START_REF] Kaufmann | On the complexity of best-arm identification in multi-armed bandit models[END_REF] (Lemma 1 therein). Theorem 2. There exists a bandit instance ν with |S |= p ≥ 3 such that for k ∈ {p, p -1, p -2} any δ-correct algorithm for 0-PSI-k verifies

E ν (τ δ ) ≥ 1 D log 1 δ K a=1 1 (∆ k a ) 2 ,
where ∆ k a := ∆ a + ω k and τ δ is the stopping time of the algorithm.

Proof. Let p = K = |S |. w.l.o.g assume S = {1, . . . , p} and S ,k = {1, . . . , k}. Let µ 0 ∈ R D and for (p -2) ≤ i ≤ p, define

µ d i :=    2 p-i ω if d = 1 -2 p-i ω else if d = 2 µ d 0 else.
,

for 1 ≤ i ≤ p -3,
µ d i :=    (4 + 2i)ω if d = 1 -(4 + 2i)ω else if d = 2 µ d 0 else. for which P ν (E i ) ≥ 1 2 and P ν (i) (E i ) ≤ δ, it comes that a∈A E ν (T a (τ δ )) KL(ν a , ν (i) a ) ≥ d(P ν (E i ), P ν (i) (E i )), hence E ν (T i (τ δ )) KL(ν i , ν (i) i ) ≥ d(P ν (E i ), P ν (i) (E i )), (22)
where d(x, y) = x log(x/y) + (1 -x) log((1 -x)/(1 -y)) is the binary relative entropy. Since

P ν (E i ) ≥ 1 2 and P ν (i) (E i ) ≤ δ, (22)
yields (see [START_REF] Kaufmann | On the complexity of best-arm identification in multi-armed bandit models[END_REF]),

E ν (T i (τ δ )) ≥ 1 KL(ν i , ν (i) i ) 1 2 log 1 2δ + log 1 2(1 -δ) = 1 2 KL(ν i , ν (i) i ) log 1 δ(1 -δ) ≥ 1 2 KL(ν i , ν (i) i )
log(1/δ).

By direct algebra, we compute (independent marginals since the covariance is diagonal I D ),

KL(ν i , ν (i) i ) = 1 2 µ j -δ + i -ω (k) -γ -µ j 2 2 = D 2 (-δ + i -ω (k) -γ) 2 .
Noting that on this instance all the arms are optimal, we have for any arm i, ∆ i = δ + i . Finally, letting γ -→ 0 proves that for any arm i,

E ν (T i (τ δ )) ≥ 1 D(∆ k i ) 2 log(1/δ), further noting that E(τ δ ) = K i=1 E(T i (τ δ )
) achieves the proof. We have chosen a diagonal matrix matrix simplicity, we believe that choosing carefully correlated objectives like in [START_REF] Auer | Pareto front identification from stochastic bandit feedback[END_REF] could give a tighter lower bound especially regarding the dependence in the dimension D.

E Best Arm Identification

In this section, we discuss the sample complexity and the performance of APE associated to the stopping rule τ 1 0 for BAI. Noting that when D = 1, the Pareto set is just the argmax over the means, BAI and PSI are the same for uni-dimensional bandits. For this setting we should expect algorithms for PSI to be competitive with existing algorithms for BAI. We will show that it is actually the case for APE . Let D = 1 and ν be a one-dimensional K-armed bandit. Letting a denote the unique optimal arm of the bandit ν, i.e S = {a }, one can easily check that the gaps defined for PSI matches the common notion of gaps for BAI. Indeed, for any a = a ,

∆ a := max j∈S m(a, j), = m(a, a ) = µ a -µ a ,
and

∆ a = min j =a {M(j, a ) + + ∆ j } = min j =a ∆ j ,
which matches the definition of the gap in the one-dimensional bandit setting ( [START_REF] Audibert | Best arm identification in multi-armed bandits[END_REF][START_REF] Kaufmann | On the complexity of best-arm identification in multi-armed bandit models[END_REF]). Therefore, the sample complexity of APE for BAI can be deduced from Theorem 1. Theorem 4. Let δ ∈ (0, 1), K ≥ 2 and ν a K-armed bandit with a unique best arm a and 1-subgaussian distributions. APE associated to be stopping time τ 1 0 identifies the best arm a with probability at least 1 -δ using at most the following number of samples

K a=1 88 ∆ 2 a log 2K(K -1) δ log 12e ∆ a .
In particular, the k relaxation is not meaningful in this setting. Under the unique optimal arm assumption, the algorithm will always stop when the best arm has been identified. And we remark that from the definition of ω i 's

ω 1 = min j =a M(a , j) = min j =a ∆ j = ∆ a and ∀i = a , ω i < 0, so for any k ≤ K, max(ω k , ∆ a ) = ∆ a .
Remark 2. Theorem 4 could be slightly improved. On the event E we consider that for any pair of arms the difference of their empirical mean does not deviate too much from its actual value. For BAI, since we know that there is a unique optimal arm (enforced by assumption), it is sufficient to control the difference between the best arm and any other arm, therefore we could replace the K(K -1)/2 term due to union bound in the confidence bonus by K -1 and we could show that this will reflect in the sample complexity by replacing K(K -1) by 2(K -1). However, this cannot be done in general for PSI since we do not know in advance the number of optimal arms.

When D = 1, APE reduces to sample at each round t, the least sampled among

b t := argmax i min j =i U i,j (t) , (23)
c t := argmin j =bt L bt,j (t), (24)
where U i,j (t) := µ i (t) -µ j (t) + β i,j (t) and L i,j (t) := µ i (t) -µ j (t) -β i,j (t) are upper and lower bounds on the difference µ i -µ j . To be in the same setting as LUCB and UGapEc which uses confidence interval on single arms, we would have β i,j (t) := β i (t) + β j (t), where β i s are confidence bonuses on single arms such that L i (t) := µ i (t) -β i (t) and U i (t) := µ i (t) + β i (t) are lower and upper confidence bounds on µ i . Then ( 23) and ( 24) rewrite as

b t := argmax i U i (t) -max j =i L j (t) , c t := argmax j =bt U j (t),
This resembles the sampling rule of UGap, which defines

b UGap t := argmax i L i (t) -max j =i U j (t) , c UGap t := argmax j =bt U i (t),
and also pulls the least sampled so far. We note that a variant of our algorithm in which both b t , c t would be sampled (in the spirit of LUCB [START_REF] Kalyanakrishnan | PAC subset selection in stochastic multiarmed bandits[END_REF]) could also be analyzed using the same arguments employed in the proof of Theorem 1.

Note that when ε 1 = 0, for any i ∈ S(t) c , g i (t) > h 0 i (t). Indeed, by definition,

h 0 i (t) = min j =i (M(i, j, t) -β i,j (t)) = min j =i (-m(i, j, t) -β i,j (t))
and since i ∈ S(t) c , there exists i such that m(i, i , t) > 0 (i.e µ i (t) ≺ µ i (t)) and so

-m(i, i , t) -β i,i (t) < m(i, i , t) -β i,i (t).
Therefore,

min j =i (-m(i, j, t) -β i,j (t)) := h 0 i (t) < max j =i (m(i, j, t) -β i,j (t)) := g i (t).
Thus for

ε 1 = 0, Z 0 2 (t) = min i∈S(t) c g i (t). (25)
In the sequel, for this section, we remove the dependence on ε 1 to write Z i (t) instead of Z 0 i for i = 0 and i = 1. In particular, when D = 1, ε 1 = 0, the stopping time τ 0 can be simplified to

τ 0 = inf{t ≥ K : Z 1 (t) > 0},
which is a consequence of the following lemma.

Lemma 12. For D = 1, ε 1 = 0, inf{t ∈ N : Z 1 (t) > 0} = inf{t ∈ N : Z 1 (t) > 0 ∧ Z 2 (t) > 0}.
Proof of Lemma 12. Let S(t) = { a t }. Using the definition of h 0 i , g i and ( 25), Z 1 (t) and Z 2 (t) simplifies to

Z 1 (t) = min i = at { µ at (t) -µ i (t) -β at,i (t)} , (26)
Z 2 (t) = min i = at max j =i [ µ j (t) -µ i (t) -β i,j (t)] . (27)
We have :

Z 1 (t) > 0 =⇒ ∀i = a t , µ at (t) -µ i (t) -β at,i (t) > 0, =⇒ ∀i = a t , max j =i [ µ j (t) -µ i (t) -β j,i (t)] > 0, =⇒ Z 2 (t) > 0.
Thus, Z 1 (t) > 0 =⇒ (Z 1 (t) > 0 ∧ Z 2 (t) > 0) and the reverse holds trivially. So

Z 1 (t) > 0 ⇐⇒ (Z 1 (t) > 0 ∧ Z 2 (t) > 0).
Letting a t denote the empirical best arm after t rounds, the stopping rule of APE (with the instantiation proposed in Section 3.2 based on confidence intervals on pairs of arms) reduces to

τ 0 = inf    t ≥ K : ∀i = a t , ( µ at (t) -µ i (t)) 2 2 1 T a t (t) + 1 Ti(t) ≥ 2C g log(K 1 /δ) 2 + 2 a∈{ at,i} log(4 + log(T a (t)))
   which is very close to a Generalized Likelihood Ratio (GLR) stopping rule assuming Gaussian distributions with variance 1 for the rewards (which is known to be also correct for sub-Gaussian rewards) [START_REF] Garivier | Optimal best arm identification with fixed confidence[END_REF][START_REF] Kaufmann | Mixture martingales revisited with applications to sequential tests and confidence intervals[END_REF]. This modified stopping rule compared to LUCB1 and UGapEc can partially explains the empirical improvement observed in Section H.3.

F LUCB1-like instantiation of APE

In this section we derive an upper bound on the expectation of the sample complexity τ k ε1 when APE is run with confidence bonuses similar to LUCB1 [START_REF] Kalyanakrishnan | PAC subset selection in stochastic multiarmed bandits[END_REF]. This is different from Theorem 1 for which the sample complexity is bounded only on the high-probability event E but as, for many algorithms in pure-exploration [START_REF] Simchowitz | The simulator: Understanding adaptive sampling in the moderate-confidence regime[END_REF][START_REF] Gabillon | Best arm identification: A unified approach to fixed budget and fixed confidence[END_REF][START_REF] Auer | Pareto front identification from stochastic bandit feedback[END_REF] we do not control what happens on E c . Therefore, our goal here is to upper-bound E(τ ) instead of E(1{E}τ ) which we did in Theorem 1. To adapt the strategy employed in [START_REF] Kalyanakrishnan | PAC subset selection in stochastic multiarmed bandits[END_REF], we use similar confidence bonuses, thus we define for any arm i,

β i (t) = 2 T i (t) log 5KDt 4 2δ , (28)
and for any pair i, j ∈ A, β i,j (t) = β i (t) + β j (t). Recalling the definition of E and E t introduced in Section 3.1,

E t := K i=1 j =i D d=1 L d i,j (t, δ) ≤ µ d i -µ d j ≤ U d i,j (t, δ) , and E = ∞ t=1 E t ,
the lemma hereafter shows that with the choice of β i 's in [START_REF] Zuluaga | Active learning for multi-objective optimization[END_REF] and for 1-subgaussian marginals, P(E) ≥ 1 -δ.

Lemma 13. It holds that P(E) ≥ 1 -δ.

Proof. Letting

E t := K i=1 D d=1 | µ d i (t) -µ d i |≤ β i (t) ,
we have E t ⊂ E. Indeed, on E t , for any i, j ∈ A and d ≤ D,

µ d i (t) -µ d j (t) -β i (t) -β j (t) ≤ µ d i -µ d j ≤ µ d i (t) -µ d j (t) + β i (t) + β j (t),
which combined with β i,j (t) = β i (t) + β j (t) yields E t ⊂ E t so

P(E c ) ≤ ∞ t=1 P(E c t ) ≤ ∞ t=1 P( E c t ). (29)
Applying Hoeffding's inequality to the 1-subgaussian marginals yields

P( E c t ) ≤ K i=1 D d=1 P | µ d i (t) -µ d i |> β i (t) , ≤ K i=1 D d=1 t s=1 P(| µ d i,s -µ d i |> β t,s ) where β t,s = 2 s log 5KDt 4 2δ , ≤ K i=1 D d=1 t s=1 4δ 5KDt 4 , = 4δ 5t 3 .
Finally,

P(E c ) ≤ ∞ t=1 P( E c t ), ≤ 4δ 5 ∞ t=1 1 t 3 , ≤ δ.
We can now state the main theorem of this section. Theorem 5. Let ε 1 ≥ 0, k ≤ K and ν a bandit with 1-subgaussian marginals. APE run with the β i s of (28) and associated to the stopping time τ k ε1 outputs a valid set and its expected sample complexity is upper-bounded as follows :

E ν (τ k ε1 ) ≤ 64 √ eH log 5KD 2δ + 256 √ eH log (256H) + 8π 2 15 + 1, with H := K a=1 max(∆ at , ε 1 , ω k ) -2 .
Proof of Theorem 1. The correctness follows from Lemma 2 combined with Lemma 13. It remains to upper-bound E(τ k ε1 ). Note that this proof technique has been already used in [START_REF] Kalyanakrishnan | PAC subset selection in stochastic multiarmed bandits[END_REF][START_REF] Kaufmann | Information complexity in bandit subset selection[END_REF][START_REF] Jun | Anytime exploration for multi-armed bandits using confidence information[END_REF] for LUCB-like algorithms. Let n ≥ 1 to be specified later and

E(n) = t∈[ 1 2 n,n] E t . (30)
Remark that

E(n) ∩ {τ k ε1 > n} holds =⇒ n t=1 1{{τ k ε1 > t} ∩ E(n)} = n. (31)
We will show that for some choice of n, the RHS of (31) will be strictly less than n so the LHS does not hold. We proceed by upper-bounding the RHS

n t=1 1{{τ k ε1 > t} ∩ E(n)} ≤ n 2 + n t= n 2 1{{τ k ε1 > t} ∩ E(n)}, (32)
≤ n 2 + n t= n 2 1{{τ k ε1 > t} ∧ E t }. (33)
From Lemma 4, Lemma 5 and Lemma 3, we have that for any t ∈ [ n 2 , n], {τ k ε1 > t} ∩ E t =⇒ max(∆ at , ε 1 , ω k ) ≤ 2β at,at (t), with β at,at (t) = 2β at (t). Therefore, using this result back in (32) and letting c δ := (5KD/(2δ))

1/4 , ∆ a := max(∆ a , ε 1 , ω k ) yields n t=1 1{{τ k ε1 > t} ∩ E(n)} ≤ n 2 + n t= n 2 1{ ∆ at ≤ 4β at (t)} ≤ n 2 + n t= n 2 K a=1 1{(a t = a) ∧ ∆ a ≤ 4β a (t)}, ≤ n 2 + K a=1 n t= 1 2 n 1{{a t = a} ∧ {T a (t) ≤ 128 ∆ 2 a log(c δ t)}} ≤ n 2 + K a=1 n t= n 2 1{{a t = a} ∧ {T a (t) ≤ 128 ∆ 2 a log(c δ n)}} ≤ n 2 + K a=1 128 ∆ 2 a log(c δ n) ≤ n 2 + 128H log(c δ n),
where H := a ∆ -2 a . Then, choosing n such that

n 2 + 128H log(c δ n) < n, that is n > T := inf s ∈ N : 128H log(c δ s) s < 1 2 , ( 34
) would yield n t=1 1{{τ k ε1 > t} ∩ E(n)} < n, (35) so
E(n) ∩ {τ k ε1 > n} = ∅, which means {τ k ε1 > n} ⊂ E(n) c . Therefore, for any n > T , {τ k ε1 > n} ⊂ E(n) c . (36) Thus, E ν (τ k ε1 ) = E ν τ k ε1 1{τ k ε1 ≤ T } + τ k ε1 1{τ k ε1 > T } ≤ T + E ν (τ k ε1 1{τ k ε1 > T } ≤ T + ∞ n=T +1 P ν (τ k ε1 > n) ≤ T + ∞ n=T +1 P ν (E(n) c ),
using (30) and union bound yields,

P(E(n) c ) ≤ n t= n 2 4δ 5t 3 , ≤ 4δ 5
( 1 /2)n ( 1 /2) 3 n 3 , = 16δ 5
1 n 2 . Then, E ν (τ k ε1 ) ≤ T + 16δ 5 π 2 6 ≤ T + 8π 2 15 .
Upper-bounding T will conclude the proof.

Lemma 14. It holds that

T -1 ≤ 1 c δ exp -W -1 - 1 256c δ H ≤ 256 √ eH log (256c δ H) .
Finally,

E ν (τ δ ) ≤ 256 √ eH log 256(5KD/(2δ)) 1/4 H + 8π 2 15 + 1, (37)
≤ 64 √ eH log 5KD 2δ + 256 √ eH log (256H) + 8π 2 15 + 1, (38)
which achieves the proof.

Remark 3. The same technique could be applied to upper-bound E ν (τ ε1,ε2 ).

Now we prove Lemma 14.

Proof of Lemma 14. We have

128H log(c δ s) s < 1 2 ⇐⇒ log(c δ s) s < 1 256H (39)
then, using Lemma 15 yields

(39) =⇒ s > 0 if 1 256H < c δ /e 0 < s ≤ 1 c δ or s ≥ N else, with N = 1 c δ exp -W -1 - 1 256c δ H .
Therefore,

T = inf s ∈ N : 128H log(c δ s) s < 1 2 ≤ 1 if 1 256H < c δ /e N else.
Using Corollary 1, to upper bound N yields

T -1 ≤ 256 √ eH log + (256c δ H) , (40)
where log + (x) = max(0, log(x)).

G Technical Lemmas

Lemma 15. Let a, b > 0. If b < a/e then log(ax)

x < b =⇒ 0 < x ≤ 1 a or x ≥ 1 a exp -W -1 - b a .
Moreover, if b ≥ a/e, then for any x > 0, log(ax)/x ≤ b.

Proof. We have

log(ax) x < b =⇒ - 1 ax log 1 ax < b a =⇒ y log(y) > - b a , y = 1 ax =⇒ 1 ax ≥ 1 or - b a < y log(y) < 0
since -b/a > -1/e and the negative branch W -1 of the Lambert function is decreasing on [-1/e, 0],

log(ax) x < b =⇒ 0 < x ≤ 1 a or W -1 (y log(y)) ≤ W -1 (-b/a) =⇒ 0 < x ≤ 1 a or log(y) ≤ W -1 (-b/a) =⇒ 0 < x ≤ 1 a or ax ≥ exp(-W -1 (-b/a)) =⇒ 0 < x ≤ 1 a or x ≥ 1 a exp(-W -1 (-b/a)).
Proving the second part of the lemma just follows from log(x) ≤ x/e.

The following lemma is taken from [START_REF] Jourdan | Top two algorithms revisited[END_REF] Lemma 16 ([12]). For any x ∈ [0, -e -1 ],

-log(-x)+log(-log(-x)) ≤ -W -1 (x) ≤ -log(-x)+log(-log(-x))+min 1 2 , 1

-x log(-x) Corollary 1. Let 0 < a < 1/e. It holds that exp(-W -1 (-a)) ≤ e 1/2 a log 1 a .
We recall the following lemma which is taken from [START_REF] Kaufmann | Mixture martingales revisited with applications to sequential tests and confidence intervals[END_REF].

Proof. Using Lemma 16 yields,

-W -1 (-a) ≤ -log(a) + log(-log(a)) + 1
and by definition of W 0 , we have

W 0 (x) exp(W 0 (x)) = x, so exp W 0 e 4 ∆ 2 = e 4 ∆ 2 1 W 0 (e 4 ∆ -2 ) ,
and therefore,

t 0 = 1 ∆ 2 1 W 0 (e 4 ∆ -2 )
.

We will show that t ∆ > t 0 . Indeed, since W 0 is increasing,

1 ∆ 2 > 3/e =⇒ W 0 e 4 ∆ 2 ≥ W 0 (3e 3 ) = 3 =⇒ 1 ∆ 2 1 W 0 (e 4 ∆ -2 ) ≤ 1 3 1 ∆ 2 , that is t 0 ≤ 1 3 1 ∆ 2 . ( 44
)
On the other side,

1 ∆ 2 > 3/e =⇒ log 2 log 3e 2 2∆ 2 > log(2 log(9e/2)) =⇒ t ∆ ≥ 1 ∆ 2 log(2 log(9e/2)) > 1 3 1 ∆ 2
Therefore,

t 0 ≤ 1 3 1 ∆ 2 < log(2 log(9e/2)) ∆ 2 ≤ t ∆ .
Thus, we have shown that t 0 ≤ t ∆ and for any t ≥ t ∆ , g (t) ≥ 0 so

∀t ≥ t ∆ , g(t) ≥ g(t ∆ ). (45)
Showing that g(t ∆ ) > 0, will conclude the proof. Letting a = 3e 2 /2, we have where g is a non-negative function. Then, for any α ∈ (0, 1) and t ≥ 2,

g(t ∆ ) > 0 ⇐⇒ 1 ∆ 2 log(2 log(a/∆ 2 )) - 1
t ≥ 1 ∆ 2 1 α g(δ) + c 1 -α log + 2 log c (1 -α)∆ 2 =⇒ f (t) < ∆. (54)
Proof. Letting t ≥ 2, we have

t ≥ t 1 := 1 α 1 ∆ 2 g(δ) =⇒ g(δ) t ≤ α∆ 2 . (55)
Furthermore, using Lemma 17 yields

t ≥ t 2 := c (1 -α)∆ 2 log + 2 log 3e 2 2(1 -α)∆ 2 =⇒ log log(e 4 t) t ≤ (1 -α)∆ 2 /c. (56)
Combining ( 55) and (56) yields for t ≥ 2,

t ≥ max(t 1 , t 2 ) =⇒ f (t) 2 < ∆ 2 , so t ≥ t 1 + t 2 ≥ max(t 1 , t 2 ) =⇒ f (t) < ∆. (57

Time and memory complexity

The time complexity of ε 1 -APE-k is O(K 2 D) and its memory complexity is O(K 2 ). The main computational bottleneck is the computation of the M (i, j, t) for each (i, j) ∈ K × K, which requires a triple-nested for-loop over

[K] × [K] × [D].
To give an idea of the runtime, a single run on a random Bernoulli instance with K = 1000, D = 10 takes around 4 minutes for 0.1-APE-1000 on a personal computer (a single 3GHz ARM core used, 8 GB RAM, 256 GB disk storage) with δ = 0.01 and no particular optimization. Due its fully sequential nature, our algorithm may have a higher computational cost compared to uniform sampling strategies. However, in our implementation the most time-consuming operation was actually collecting a sample from the selected arm(s), especially for multivariate Gaussians. So that finally, in the experiments, our algorithm which ultimately require less samples had in practice a similar computational cost compared to PSI-Unif-Elim which uses uniform sampling.

Adaptation to bandits with marginals of different scaling We have presented the algorithm and the results specialized to the case where all the marginals are 1-subgaussian. Indeed our results can be simply extended where the marginals are instead all σ-subgaussian. Furthermore, there is a simple way to adapt the algorithm to the case where the marginals have different known subgaussianity parameter (i.e different scaling) but they are the same for every arm. The idea is to rescale each observation with the subgaussianity parameters. Let σ := σ 1 , . . . σ D , σ i > 0. Assuming that the marginal distributions of each arm are respectively σ 1 , . . . σ D -subgaussian, each observation X At , s from arm A t will be rescaled component-wise to X d At,s /σ d before being given to the algorithm. It is easy to see that this rescaling does not change the Pareto set since all the means are divided by the same values coordinate-wise. Furthermore, by defining

M σ (i, j) := max d µ d i -µ d j σ d ,
and m σ (i, j) =: -M σ (i, j), all the results proved for 1-subgaussian distributions still holds using m σ and M σ in the definition of the gaps (Section 4).

H.2 Data processing

Dataset The dataset is extracted from [START_REF] Munro | Safety and immunogenicity of seven COVID-19 vaccines as a third dose (booster) following two doses of ChAdOx1 nCov-19 or BNT162b2 in the UK (COV-BOOST): a blinded, multicentre, randomised, controlled, phase 2 trial[END_REF] and some processing steps are applied to compute the covariance matrix of the distribution. First, as observed in [START_REF] Munro | Safety and immunogenicity of seven COVID-19 vaccines as a third dose (booster) following two doses of ChAdOx1 nCov-19 or BNT162b2 in the UK (COV-BOOST): a blinded, multicentre, randomised, controlled, phase 2 trial[END_REF], the 3 immunogenicity indicators extracted are weakly correlated, therefore, we assume the covariance matrix to be diagonal. To compute the variance of the marginals, we use the log-normal assumption as assumed for the data reported in [START_REF] Munro | Safety and immunogenicity of seven COVID-19 vaccines as a third dose (booster) following two doses of ChAdOx1 nCov-19 or BNT162b2 in the UK (COV-BOOST): a blinded, multicentre, randomised, controlled, phase 2 trial[END_REF]. Using this log-normal assumption the authors have provided for each arm and each indicator: the geometrical mean, the sample size and a 95% confidence interval on the geometrical mean based on the central limit theorem.

For each of the K = 20 arms (combination of three doses), we use these information to compute the sample variance of each immunogenicity indicator. Moreover, we compute the arithmetic average of the log outcomes which is obtained by taking the log (base e) of the geometrical empirical mean:

x = log(x geometrical ), = log   n i=1 x i 1/n   , = n -1 n i=1 log(x i ),
where x 1 , . . . , x n are the observations which are assumed to be log-normal. x represents by assumption the empirical mean of a Gaussian distribution, which we use as a proxy for its true, unknown mean. From this, we built a bandit model where each arm is a 3-dimensional Gaussian distribution with independent coordinates, whose means are given by the corresponding mean estimates (reported in Table 4) and in which the variance of each indicator is the pooled variance over the different arms (given in In this section we show that for some instances our algorithm can require up to 3 times less samples compared to PSI-Unif-Elim . This is due to the strategy of PSI-Unif-Elim which continue sampling arms identified as optimal until there are shown not to dominate any arm in the active set. For example on Figure 3a, the optimal arm 2 is "easy" to identify as such. However, since it slightly dominates the sub-optimal arm 1, PSI-Unif-Elim should continue sampling arm 2 until arm 1 is removed from the active set ( likely this will happen when the algorithm "sees" that arm 1 is dominated by arm 3). We would expect our adaptive sampling rule to avoid this behaviour.

Figure 3b shows that APE takes nearly half the average sample complexity of PSI-Unif-Elim on this instance. In particular, Table 6 shows the average number of pulls taken by PSI-Unif-Elim divided by the average number of pulls taken by 0-APE-K for each arm. We can observe that the major difference in sample complexity is due to arm 2 being pulled nearly 6 times more by PSI-Unif-Elim than APE . 6: Average number of pulls taken by PSI-Unif-Elim divided by the average number of pulls taken by 0-APE-K for each arm.

By increasing the number of arms and the dimension we can generate instances similar to Figure 3a where the gap between our algorithm and PSI-Unif-Elim is even larger. We chose a specific instance where K = 12, D = 10 and there are 11 optimal arms. On this instance (Figure 4), we can see that our algorithm uses 3 times less samples than PSI-Unif-Elim . Finally, combining these additional experiments with the results of Table 2 we observe that on average ε 1 -APE-k performs nearly 20% better than APE but there are some instances where the gap can be even larger. Of course, this also means that there should exist instances in which the improvement is smaller than 20% to compensate for instances like Figure 3a. But we note that instances like Figure 3a are very unlikely to be generated randomly so they should only be a few in the 2000 instances used for Figure 3a. So that 20% is fairly representative of the average improvement on "normal instances" (excluding instances like Figure 3a where the improvement can be way larger). We investigate the empirical behavior of (ε 1 , ε 2 )-APE for identifying an (ε 1 , ε 2 )-cover. We set ε 1 = 0, δ = 0.01 and we test different values of ε 2 ∈ {0, 0.05, 0.1, 0.2, 0.25}. We average the results over 2000 independent trials with different seeds on the same instance (Figure 5a). We use multi-variate Bernoulli with independent marginals. The instance of Figure 5a is a toy example where (ε 1 , ε 2 )-covering can be meaningful and reduce the sample complexity. The 3 Pareto optimal vectors are chosen by hand and the last 2 vectors are randomly uniformly generated.

We can observe on Figure 5b that the sample complexity decreases as ε 2 increases. This is further confirmed in Figure 6 which shows the empirical sample complexity and the average size of the recommended cover versus ε 2 for 50 equally-spaced values of ε 2 between 0 and 1/2. The drops observed in Figure 6b correspond to the values of ε 2 for which APE removes an optimal arm from the cover to save some samples (Figure 6a). A major decrease in the sample complexity corresponds more or less to an arm being removed from the recommended set. We observe in Figure 7 the histogram of occurrence of each arm in the recommended set for 3 values of ε 2 corresponding more or less to the middle of each plateau. We can see that for ε 2 = 0.15, arm 0 is always recommended, but the others are recommended on half of the runs. For ε 2 = 0.4, the algorithm nearly always recommend arm 0, which as the largest ω i term (i.e the easiest to identify as optimal).

H.3.4 Comparison to some BAI algorithms

We evaluate the performance of APE for Best Arm Identification (D = 1) on two randomly generated instances: one with K = 5 and means (rounded) X 1 := (0.25, 0.16, 0.87, 0.22, 0.98), and the second one with K = 10 and means (rounded) X 2 := (0.43, 0.33, 0.56, 0.85, 0.20, 0.93, 0.70, 0.82, 0.56, 0.78) We use the instantiation of APE with confidence bonuses on pair of arms but without the possible improvement invoked in Remark 2. UGap and LUCB are implemented with the tightest known (to our knowledge) confidence bonus taken from [START_REF] Katz-Samuels | Top feasible arm identification[END_REF] (in spirit of the finite-time law of the iterated logarithm [START_REF] Jamieson | lil' ucb : An optimal exploration algorithm for multi-armed bandits[END_REF]). LUCB++ is used with the improved scheme given in [START_REF] Simchowitz | The simulator: Understanding adaptive sampling in the moderate-confidence regime[END_REF]. We set δ = 0.01 but the empirical error was way smaller. The results are averaged over 1000 independent trials. On the y-axis is the sample complexity in units of (an approximation of) the lower bound of BAI for Gaussian rewards with σ = 1 /24 ( [START_REF] Kaufmann | On the complexity of best-arm identification in multi-armed bandit models[END_REF][START_REF] Simchowitz | The simulator: Understanding adaptive sampling in the moderate-confidence regime[END_REF]) :

Figure 1 :

 1 Figure 1: On the left is the log of the empirical sample complexity of PSI-Unif-Elim and APE on the real-world scenario plot (right) for 2 out of the 3 immunogenicity indicators.

  (a) Average sample complexity vs ε2 (b) Average ratio |O(τ 0,ε 2 )| /|S | vs ε2

Figure 2 :

 2 Figure 2: APE with τ 0,ε2 averaged over 2000 random Bernoulli instance with K = 5 arms.

1

 1 Generic algorithm(s) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Our instantiation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . perspective A Correctness for different stopping rules A.1 Proof of Lemma 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.2 Proof of Lemma 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.3 Calibration of the confidence intervals . . . . . . . . . . . . . . . . . . . . . . . . B Sample complexity analysis B.1 Proof of Lemma 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.2 Proof of Lemma 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.3 Proof of Lemma 5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.4 Auxiliary results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C Algorithm for finding an (ε 1 , ε 2 )-cover D Lower Bound E Best Arm Identification F LUCB1-like instantiation of APE G Technical Lemmas H Implementation and Additional Experiments H.1 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . H.2 Data processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . H.3 Additional experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . H.3.1 Additional experiments for ε 1 -APE-k . . . . . . . . . . . . . . . . . . . . H.3.2 (ε 1 , ε 2 )-APE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . H.3.3 Comparing ε 1 -APE-k to an adaptation of PSI-Unif-Elim . . . . . . . . . . H.3.4 Comparison to some BAI algorithms . . . . . . . . . . . . . . . . . . . . .

Figure 3 :

 3 Figure 3: For this instance (left) S = {0, 2, 3} and on the right is the average sample complexity over the trials on the same instance.

Figure 4 :

 4 Figure 4: Average sample complexity on a specific instance with K = 12, D = 10 and |S | = 11

Figure 5 :

 5 Figure 5: For this instance (left) S = {0, 1, 2}. The difference in on x and y axis is 0.1 between arm 0 and 1 and 0.05 between arm 1 and 2. The rightmost figure is the sample complexity averaged over 2000 trials.

2 (

 2 Average |O(τε 1 ,ε 2 )| vs ε2 b) Sample complexity vs ε2

Figure 6 :

 6 Figure 6: On the left is the average length of O(τ ε1,ε2 ) (over the different runs) versus ε 2 and on the right is the empirical sample complexity (averaged over the runs) versus ε 2 . The empirical probability of error was equal to zero.

  K = 50, D = 2, ε1 = 0.1

Figure 8 :

 8 Figure 8: Average sample complexity on 2000 random instances (left) and 500 random instances (right).

Figure 9 :

 9 Figure 9: Empirical sample complexity expressed in units of the lower bound H log(1/2.4δ) (blue line) on a random instance with K = 5 (left) and K = 10 (right), δ = 0.01. The blue line has a coordinate of 1 on the y-axis (lower-bound).

Table 1 :

 1 Stopping conditions and associated recommendationStopping and recommendation rule(s) Depending on the objective, Adaptive Pareto Exploration can be plugged in with different stopping rules, that are summarized in Table1with their associated recommendations. To define those, we define for all i ∈ A, ε 1 , ε 2 ≥ 0,

	g ε2 i (t) := max j =i	m -(i, j, t) + ε 2 1{j ∈ OPT ε1 (t)} and h ε1 i (t) := min j =i	M -(i, j, t) + ε 1 .
	and let g i (t) := g 0 i (t). Introducing
	Z ε1 1 (t) := min i∈S(t)	h ε1 i (t), and Z ε1 2 (t) := min i∈S(t) c max(g i (t), h ε1 i (t)),
	for ε 1 -PSI, our stopping rule is τ ε1 := inf{t ≥ K : Z ε1 1 (t) > 0 ∧ Z ε1 2 (t) > 0} and the associated
	recommendation is O(τ ε1 ) where

Table 2 :

 2 Average sample complexity over 2000 random Bernoulli instances with K = 5 arms. On average the size of the Pareto set was (2.295, 4.0625, 4.931) respectively for the dimensions 2, 4, 8.

	2	811	39530	109020	145777	150844	190625
	D = 4	214	6410	19908	68061	124001	157584
	D = 8	119	204	405	1448	20336	27270

  M(b t , c t , t) + β bt,ct (t),

	µ d ct (t) -µ d bt (t) + β bt,ct (t) ≥ µ d c t (t) -µ d bt (t) + β bt,c t (t)	Et ≥ µ d c t -µ d bt ,
	(b)

≤ 2β bt,ct (t), where (a) follows from E t and (b) has been already justified in the case 1.1. If b t = c t , then by definition of c t , we have m(b t , c t , t) + β bt,ct (t) ≥ m(b t , c t , t) + β bt,c t (t), which implies that there exists d ∈ [D] such that

  ≥ 0 for t such that ∆ 2 t log(e 4 t) ≥ 1. Using the Lambert function W 0 , which is increasing on [0, ∞), ∆ 2 t log(e 4 t) ≥ 1 ⇐⇒ e 4 t log(e 4 t) ≥

	Therefore, in the sequel, we assume ∆ 2 < e/3. Let
		t ∆ :=	1 ∆ 2 log 2 log	3e 2 2∆ 2	,
	and
		g(t) = t -	1 ∆ 2 log(log(e 4 t)).
	Then,
		g (t) = 1 -	1 ∆ 2 t log(e 4 t)	,
	and g (t) e 4 ∆ 2	(41)
			⇐⇒ log(e 4 t) ≥ W 0	e 4 ∆ 2	(42)
			⇐⇒ t ≥ t 0 :=	1 e 4 exp W 0	e 4 ∆ 2
							2	,
	and taking exp on both sides gives the result.
	Lemma 17. Let ∆ 2 > 0. Then, for t ≥ 2,
	t ≥	1 ∆ 2 log 2 log	3e 2 2∆ 2	=⇒	log log e 4 t t	< ∆ 2 .
	Proof. We note that if ∆ 2 ≥ e 3 , then the result follows trivially since it can be easily checked that for t ≥ 2, log log(e 4 t) ≤ e 3 t.

  ∆ 2 log(log(e 4 t ∆ )) > 0 (46) ⇐⇒ log(2 log(a/∆ 2 )) -log(log(e 4 t ∆ )) > 0 (47) ⇐⇒ 2 log(a/∆ 2 ) -log(e 4 t ∆ ) > 0 (48) ⇐⇒ log(a/∆ 2 ) -log(e 4 t ∆ ∆ 2 /a) > 0

						(49)
	⇐⇒	a ∆ 2 -	e 4 a	∆ 2 t ∆ > 0	(50)
	⇐⇒	a ∆ 2 -	e 4 a	log(2 log(a/∆ 2 )) > 0,	(51)
	then, observing that for x ≥ 12,	log(2 log(x)) ≤	x e 2 ,
	and since	a/∆ 2 > (3e 2 /2) × (3/e) > 12,
	we have	log(2 log(a/∆ 2 )) ≤	1 e 2	a ∆ 2	(52)
	so, using (52) yields that the LHS of (51) is larger than
		3 2	e 2 1 ∆ 2 -e 2 1 ∆ 2 ,
	which is always positive. Therefore,
		∀t ≥ t ∆ , g(t ∆ ) > 0,
	that is	∀t ≥ t ∆ ,	log log(e 4 t) t	< ∆ 2 .	(53)

Lemma 18

. Let δ ∈ (0, 1), ∆ > 0 and c > 0. Let f : t → g(δ)+c log log(e 4 t) t

  )

	Lemma 19. Let ∆ > 0 and δ ∈ (0, 1). Let
	f (t) := 4	2C g (log(1/δ)/2) + 4 log log(e 4 t) t	.
	Then
	inf{n ≥ 2 : f (n) < ∆} ≤	88 ∆ 2 log	4 δ	log	12e ∆	.

Table 5 )

 5 . Sampling an arm in this bandit simulates the measurement of the (log of the) 3 immunogenicity criteria in consideration on a new patient.The 20 arms are classified into two groups. Each three/four-letters acronym denotes a vaccine candidate. Prime BNT/BNT corresponds to giving BNT as first and second dose and similarly for Prime ChAd/ChAd. For example ChAd in the group Prime BNT/BNT means to give BNT as first and second dose and ChAd as third dose (booster).

Table 4 :

 4 Table of the empirical arithmetic mean of the log-transformed immune response for three immunogenicity indicators. Each acronym corresponds to a vaccine. There are two groups of arms corresponding to the first 2 doses: one with prime BNT/BNT (BNT as first and second dose) and the second with prime ChAd/ChAd (ChAd as first and second dose). Each row in the table gives the values of the 3 immune responses for an arm (i.e. a combination of three doses).

	Dose 1/Dose 2	Dose 3 (booster)	Immune response Anti-spike IgG NT 50 cellular response
		ChAd	9.50	6.86	4.56
		NVX	9.29	6.64	4.04
		NVX Half	9.05	6.41	3.56
		BNT	10.21	7.49	4.43
	Prime BNT/BNT	BNT Half VLA	10.05 8.34	7.20 5.67	4.36 3.51
		VLA Half	8.22	5.46	3.64
		Ad26	9.75	7.27	4.71
		m1273	10.43	7.61	4.72
		CVn	8.94	6.19	3.84
		ChAd	7.81	5.26	3.97
		NVX	8.85	6.59	4.73
		NVX Half	8.44	6.15	4.59
		BNT	9.93	7.39	4.75
	Prime ChAd/ChAd	BNT Half VLA	8.71 7.51	7.20 5.31	4.91 3.96
		VLA Half	7.27	4.99	4.02
		Ad26	8.62	6.33	4.66
		m1273	10.35	7.77	5.00
		CVn	8.29	5.92	3.87

Table 5 :

 5 Pooled variance of each group.

			Immune response
		Anti-spike IgG NT 50 cellular response
	Pooled sample variance	0.70	0.83	1.54
	H.3 Additional experiments
	H.3.1 Additional experiments for ε 1 -APE-k

PAL relies on confidence intervals that follow from Gaussian process regression, but can also be instantiated with simpler un-structured confidence intervals as those used in our work and in Auer's

A random variable X is σ-subgaussian if for any λ ∈ R, E(exp(λ(X -E(X))) ≤ exp( λ 2 σ 2 2 ).

In their experiments,[START_REF] Auer | Pareto front identification from stochastic bandit feedback[END_REF] already proposed the heuristic use of confidence bonuses of this form

as Bernoulli distributions are 1 /2-subgaussian

Acknowledgments and Disclosure of Funding

Cyrille Kone is funded by an Inria/Inserm PhD grant. Emilie Kaufmann acknoweldges the support of the French National Research Agency under the projects BOLD (ANR-19-CE23-0026-04) and FATE (ANR22-CE23-0016-01).

Let ν be a bandit where each arm i is a multivariate Gaussian with mean µ i and covariance matrix I D i.e ν i ∼ N (µ i , I D ) (with I D the identity matrix in dimension D). By direct calculation, for 1 ≤ i, j ≤ p -3, M(i, j) = M(j, i) = 2ω|i -j|, and for p -2 ≤ i, j ≤ p, M(i, j) = M(j, i) = 2 p ω|2 -i -2 -j |, for i ≤ p -3 and (p -2) ≤ j ≤ p, M(i, j) = M(j, i) = (4 + 2i -2 p-j )ω ≥ 2ω.

Therefore, computing ω i and δ + i for any i ∈ [p] yields δ + i := min j∈[p]\{i} min(M(i, j), M(j, i)),

additionally, for any i ≤ p,

Thus,

Let γ > 0. For any optimal arm i, since

Moreover, it is easy to observe that for k ∈ {p -2, p -1, p} and any i ∈ [p],

we define the alternative bandit ν (i) which is also Gaussian with the same covariance matrix I D and means given by

Therefore, since A is δ-correct, and by what precedes,

The first point follows simply from the definition of δ + i and the fact that by design M(i, j) = M(i, j) for i, j ∈ [p], For the second point, if k ∈ {p, p -1}, in the bandit ν (i) , at least one arm of S (ν) is no longer optimal, then |S (ν

Letting KL denote the Kulback-Leiber divergence and using Lemma 1 of [START_REF] Kaufmann | On the complexity of best-arm identification in multi-armed bandit models[END_REF], on F τ -measurable event

Proof. We have

Therefore, letting g(δ) := 32C g (log(1/δ)/2) and c = 64 and further using Lemma 18 yields for any α ∈ (0, 1) and t ≥ 2,

where

Since C g (x) ≈ x + log(x) [START_REF] Kaufmann | Mixture martingales revisited with applications to sequential tests and confidence intervals[END_REF], and log(x) ≤ x/e we have

Taking α = α such that 16 + 16 e

that is setting

By numerical evaluation,

Therefore, putting these results together, for t ≥ 2

H Implementation and Additional Experiments

In this section, we give additional details about the experiments and additional experimental results.

H.1 Implementation

Setup We have implemented the algorithms mainly in C++17 compiled with GCC12 and interfaced with python through the cython package. The experiments are run on an ARM64 8GB RAM/8 core/256GB disk storage computer. For the function C g we have used the approximation C g (x) ≈ x + log(x) which is usually admitted [START_REF] Kaufmann | Mixture martingales revisited with applications to sequential tests and confidence intervals[END_REF]. For the experiments on real-world scenario we generate a certain number of seeds (usually 2000) and we use a different seed for each run on the same bandit. This procedure is identical for every experiment where we report the average sample complexity on the same bandit. To assess the robustness of our algorithm, the experiments on the synthetic dataset consisted in randomly uniformly sampling some bandit means for each configuration. For each sampled bandit, the algorithms compared are run once on the same instance and we note their empirical sample complexity. Finally, we report the average sample complexity across all the bandits of the same configuration.

The plateau in the sample complexity for large values of ε 2 (> 0.3) is explained by the fact the algorithm needs to identify at least one optimal arm (which is reflected in the size of the returned set Figure 6a). Indeed, for ε 1 = 0 fixed, an algorithm for (ε 1 , ε 2 )-covering still need to assert that the arms in the recommended set are truly optimal which will require some samples even when ε 2 is very large. Thus, for ε 2 > 0.3 the algorithm need to identify at least one optimal arm and we can see on Figure 6a that for these values of ε 2 , the recommended set contains only one optimal arm. Actually we can observe empirically that the "limit" sample complexity observed in Figure 6b is close to the average sample complexity of 0-APE-1 on the same instance (4073 samples).

H.3.3 Comparing ε 1 -APE-k to an adaptation of PSI-Unif-Elim

In this section, we compare ε 1 -APE-k to an adaptation of PSI-Unif-Elim which stops earlier if at least k optimal arms have been identified. The pseudo-code of the algorithm is given in algorithm 3. As shown in [START_REF] Auer | Pareto front identification from stochastic bandit feedback[END_REF], arms in P 1 (t) are already identified as optimal but when the goal is to identify the Pareto set, some of them (namely P 1 (t)\P 2 (t)) need to be sampled again until all the arms they potentially dominate are removed from A(t) and only then those arms will belong to P 2 (t) and will be removed from the active set. However, for the k-relaxation, identifying k optimal arms is enough so the algorithm can stop as soon as |P (t -1) ∪ P 1 (t)|≥ k. It this never occurs, then the algorithm will follow the initial stopping condition, that is to stop when A(t) = ∅.

Algorithm 3: Adaptation of PSI-Unif-Elim for the ε 1 -APE-k objective Data: We set δ = 0.1 and we compare both algorithms on 2 type of randomly uniformly generated Bernoulli instances. For the first type we set K = 10, D = 2 and for the second one, we set K = 50, D = 2.

For the instances with K = 10 we set ε 1 = 0.05 and run both algorithms on 2000 random Bernoulli instances with. For the second type of instances K = 50, we set ε 1 = 0.1 and we benchmark the algorithms on 500 random instances. The average size of the Pareto set was 2.90 (for K = 10, D = 2) and 4.51 (K = 50, D = 2). Figure 8 shows the average sample complexity of the algorithms for different values of k ∈ {1, . . . , 5}. We can observe that the difference between ε 1 -APE-k and APE for the 5 values reported is more important for K = 50 than for K = 10. Put together, these experiments show that our algorithm is still preferable for the k-relaxation.