Odalric-Ambrym Maillard
email: odalricmaillard@ee.technion.ac.il

Robust Risk-averse Stochastic Multi-Armed Bandits

Keywords: Multi-armed bandits, coherent risk measure, cumulant generative function, concentration of measure

We study a variant of the standard stochastic multi-armed bandit problem when one is not interested in the arm with the best mean, but instead in the arm maximising some coherent risk measure criterion. Further, we are studying the deviations of the regret instead of the less informative expected regret. We provide an algorithm, called RA-UCB to solve this problem, together with a high probability bound on its regret.

Introduction

The setting of the stochastic multi-armed bandit problem is an old and wellknown problem (see [START_REF] Thompson | On the likelihood that one unknown probability exceeds another in view of the evidence of two samples[END_REF], [START_REF] Thompson | On the theory of apportionment[END_REF] and [START_REF] Robbins | Some aspects of the sequential design of experiments[END_REF]), with a simple formalization that is nevertheless extremely powerful, leading to a large range of beautiful theoretical developments as well as important practical questions (medical treatment strategies, web advertisement, economy, etc.). The standard stochastic setting considers an agent facing a finite number of distributions (arms) that she can sample one at a time. Each sample is considered as a reward and the goal is to maximize after T trials the cumulative sum of rewards in some sense. Generally one measures performance with the expected regret criterion, that compares the expected cumulative reward of the learner to that of the strategy that constantly pulls the arm with highest mean. However in a number of applications, this criterion is not sufficient. For instance a medical treatment that is very effective on average may still have a high variability and may potentially endanger the patients. Thus, we are interested in this paper in a so-called risk-averse rather than expected measure of performance.

Risk-aversion is an old notion, however with no consensus about its definition (see [START_REF] Neumann | Theory of games and economic behavior[END_REF][START_REF] Markowitz | Portfolio selection[END_REF][START_REF] Ahmadi-Javid | Entropic value-at-risk: A new coherent risk measure[END_REF]). However, any risk measure that is coherent (see [START_REF] Rockafellar | Coherent approaches to risk in optimization under uncertainty[END_REF]) is considered to be a good measure. We here use a standard risk measure that is coherent, defined in equation [START_REF] Auer | Finite-time analysis of the multiarmed bandit problem[END_REF]. In the relevant field of reinforcement learning, risk-aversion is also quite old (see [START_REF] Howard | Risk-sensitive markov decision processes[END_REF]), and a number of works try to solve risk-averse problems [START_REF] Stephen | On terminating markov decision processes with a risk-averse objective function[END_REF][START_REF] Artzner | Coherent multiperiod risk adjusted values and bellman's principle[END_REF][START_REF] Defourny | Risk-aware decision making and dynamic programming[END_REF][START_REF] Liu | An exact algorithm for solving mdps under risksensitive planning objectives with one-switch utility functions[END_REF], although generally on the algorithmic and not on the theoretical side. Risk-aversion has been more closely looked in the online learning setting (where one sees all the rewards after pulling an arm and not only that of the chosen arm), with tight positive and negative results (see e.g. [START_REF] Even-Dar | Risk-sensitive online learning[END_REF][START_REF] Warmuth | Online variance minimization[END_REF]) on what can be done. Now in the bandit literature, the expected regret criterion has been extensively looked over the past, with recent, extremely tight, non-asymptotic results for various algorithms [START_REF] Honda | An asymptotically optimal policy for finite support models in the multiarmed bandit problem[END_REF][START_REF] Honda | Finite-time regret bound of a bandit algorithm for the semi-bounded support model[END_REF][START_REF] Garivier | The KL-UCB algorithm for bounded stochastic bandits and beyond[END_REF][START_REF] Maillard | A finite-time analysis of multi-armed bandits problems with Kullback-Leibler divergences[END_REF][START_REF] Kaufmann | Thompson sampling: An asymptotically optimal finite-time analysis[END_REF]. However, much less attention has been put on the risk-averse problem. One can cite [START_REF] Denardo | Optimal stopping, exponential utility and linear programming[END_REF] about optimality of index policies, for exponential utility and with the Gittins index perspective, which is however quite different than our goal. More recently, the work of [START_REF] Salomon | Robustness of stochastic bandit policies[END_REF] goes in the direction we target as it analyses the deviations of the regret for algorithms that compete with the best expected arm, however no risk-aversion measure is considered. Of special interest is [START_REF] Sani | Risk-aversion in multi-armed bandits[END_REF] that explicitly considers the risk-aversion problem in multi-armed bandits, targeting finite time performance guarantees. They use the very standard risk-measure called the mean-variance (see [START_REF] Markowitz | Portfolio selection[END_REF]), and show that it is possible to get sub-linear regret for such a setting, for a specific definition of regret that they introduce. However, the regret analysis is not completely satisfactory, since their notion of regret takes into account not only the variability of each arms but also the variability of the learning algorithm itself (that is, an algorithm is somehow penalized for switching between arms). Due to that, the considered regret is difficult to interpret and debatable (it is not clear whether penalizing an algorithm for switching is a desired feature).

The present work is inspired by the works of [START_REF] Rockafellar | Coherent approaches to risk in optimization under uncertainty[END_REF][START_REF] Salomon | Robustness of stochastic bandit policies[END_REF][START_REF] Maillard | A finite-time analysis of multi-armed bandits problems with Kullback-Leibler divergences[END_REF] and [START_REF] Sani | Risk-aversion in multi-armed bandits[END_REF]. We consider a notion of regret different from [START_REF] Sani | Risk-aversion in multi-armed bandits[END_REF], that we believe to be more natural and easier to interpret, where only the variability of the distribution of arms defines the regret on the one hand, while a control on the tail of the regret is provided on the other hand, similarly to [START_REF] Salomon | Robustness of stochastic bandit policies[END_REF]. We consider a coherent risk measure that generalizes the mean-variance criterion in the sense that the two measures coincide in the special case of Gaussian distributions, while the former takes into account the entire tail distribution of the random variables, not only the first two moments. We introduce the RA-UCB algorithm, inspired from [START_REF] Maillard | A finite-time analysis of multi-armed bandits problems with Kullback-Leibler divergences[END_REF] and provide a regret analysis in Proposition 1, Theorem 1 that we believe to be tight up to constant factor (note that the focus of this work is not on optimising the leading constants, such as in [START_REF] Honda | Finite-time regret bound of a bandit algorithm for the semi-bounded support model[END_REF][START_REF] Garivier | The KL-UCB algorithm for bounded stochastic bandits and beyond[END_REF][START_REF] Maillard | A finite-time analysis of multi-armed bandits problems with Kullback-Leibler divergences[END_REF], which would require a much more technical and uninformative analysis).

The paper is organized as follows. In Section 2, we introduce the regret and the coherent risk measure that is considered here, together with intuition about its meaning. In Section 3, we provide a generic robust (high-probability) nonasymptotic upper bound on the regret of any algorithm (Proposition 1) that depends only on the risk measure and on a control of the number of pulls of suboptimal arms by the algorithm. Section 4 introduces the RA-UCB algorithm, that is inspired by the K inf strategy of [START_REF] Burnetas | Optimal adaptive policies for sequential allocation problems[END_REF][START_REF] Maillard | A finite-time analysis of multi-armed bandits problems with Kullback-Leibler divergences[END_REF], together with a dual formulation that enables effective implementation. Section 5 concludes the paper with a highprobability bound on the regret of the RA-UCB (Theorem 1, Corollary 1). The analysis makes use of adaptation of concentration tools that are detailed in the Appendix.

Setup and Notations

We consider a standard multi-armed bandit game with A many unknown distributions {ν a } a=1,...,A where for each a, ν a ∈ M + 1 (R). At each time step, a learning algorithm A must choose an arm A A t ∈ {1, . . . , A}, and then receives one new sample (the reward ) from the corresponding distribution ν A A t . We write Y t ∼ ν A A t for the random reward received when the strategy A pulls the arm A A t at time t, X i,a to refer to the i th random variable sampled from arm a from the beginning of the game, and we finally introduce the quantity

N A T,a def = A t=1 I{A A t = a}.
Using these notations, the cumulated reward received by algorithm A up to time T is given by

T t=1 Y t = A a=1 N A T ,a i=1 X i,a .

Measure of risk-aversion

As mentioned in the introduction, there exists many possible ways to define risk-aversion. From a practical point of view, being risk-averse generally implies avoiding situations when we receive too bad reward (think of a medical treatment strategy, where the actions are the possible treatments, and the reward correspond to the health state of a patient). That is, we want to have a control on the tails, and more specifically on the lower-tail (the mass below the mean). More formally, let us recall that for arbitrary random variable X admitting a finite cumulant generative function around 0, then the two following properties hold (this is by a simple application of Markov's inequality)

P X ≥ inf 1 λ log E exp(λX) + log(1/δ) λ : λ > 0 ≤ δ , (1)
P X ≤ sup - 1 λ log E exp(-λX) - log(1/δ) λ : λ > 0 ≤ δ . (2)
Note that (1) measures the probability that X is big, while (2) measures the probability that X is small, which is what we want to be protected against. Now, for the sake of clarity, it makes sense to introduce the value of the cumulant generative function of the variable X at point λ, rescaled by λ, that we denote

κ λ,ν = 1 λ log E ν exp λX , (3)
and similarly we denote κ -λ,ν the value of κ λ ′ ,ν for λ ′ = -λ. This quantity is at the heart of many key-results and tools of concentration of measure (e.g. the Cramer-Chernoff method, the Chernoff transform, the log-Laplace transform).

More importantly here, κ -λ,ν is a key quantity to control the probability that X is small. We now provide more intuition for people unfamiliar with that quantity.

Example: To understand (1) and ( 2), let us consider t Gaussian random variables {Z k } k=1,...,t i.i.d. from a distribution ν with mean µ and variance σ 2 , then X = t k=1 Z k is Gaussian with mean µt and variance σ 2 t, and simple computations show that κ λ,ν = µt + λσ 2 t 2 , which yields, after optimizing the previous bounds in λ, to the optimal value λ =

2 log(1/δ) σ 2 t
and the familiar concentration bounds for Gaussian random variables

P 1 t t k=1 Z k -µ ≥ σ 2 log(1/δ) t ≤ δ and P µ - 1 t t k=1 Z k ≥ σ 2 log(1/δ) t ≤ δ .
Let us comment on this example. First, the quantity κ -λ,ν = µt -λσ 2 t 2 takes the form of an operator that measures the mean of a random variable, penalized by some higher moment (the variance in that case). This is actually a general property, since by the variational formula for the Kullback-Leibler divergence, we have for a random variable X distributed according to ν ∈ M + 1 (R) that

κ -λ,ν def = inf E ν ′ (X) + 1 λ KL(ν ′ ||ν) : ν ′ ∈ M + 1 (R) ≤ E ν X . (4)
where KL(ν ′ ||ν) denotes the Kullback-Leibler divergence between two distributions ν and ν ′ . Using κ -λ,ν as a measure of risk-aversion is natural for several reasons: Additionally to the formulation (4) and the control (2) that are important for interpretability it is also a standard coherent risk-measure (see [START_REF] Rockafellar | Coherent approaches to risk in optimization under uncertainty[END_REF]). Also, due to its deep link for concentration of measure, it is especially natural for analysis.

Mixability gaps Finally, for completeness, we also introduce the two fundamental quantities m + λ,ν X and m - λ,ν X that we call here the upper (and respectively lower) mixability gap and that are defined by

m + λ,ν = κ λ,ν -E ν X and m - λ,ν = E ν X -κ -λ,ν .
Note that the mixability gaps are always non negative by Jensen's inequality, and that an upper bound on them immediately provides a high probability confidence interval. Indeed, with these notations, the previous equations ( 1) and (2), can thus be rewritten more compactly as

P X -E ν X ≥ inf λ>0 m + λ,ν + log(1/δ) λ ≤ δ , (5)
P E ν X -X ≥ inf λ>0 m - λ,ν + log(1/δ) λ ≤ δ . (6)
2.2 Regrets for risk-averse multi-armed bandits Optimal arm We now naturally define the optimal arm a ⋆ as the one maximizing the risk aversion at some fixed level λ, that is we define

a ⋆ ∈ argmax a=1,...,A κ -λ,ν a ⋆ .
Note again that in the case of Gaussian distributions with mean µ a and variance σ 2 a , we simply have κ -λ,ν a ⋆ = µ a -λσ 2 a 2 , and that in general we always have κ -λ,ν a ⋆ ≤ E νa X .

Regret Now we define the empirical regret R T (λ) of the strategy A with respect to the strategy ⋆ that constantly pulls the same arm a ⋆ ∈ {1, . . . , A} by the difference between the cumulated reward received by algorithm A and the cumulated reward that the strategy ⋆ would have received during the same game, that is, by introducing the fictitious plays

{X i,a ⋆ } N A T ,a ⋆ <i≤T , R T (λ) def = T i=1 X i,a ⋆ - A a=1 N A T ,a i=1 X i,a = T i=N A T ,a ⋆ +1 X i,a ⋆ - a =a ⋆ N A T ,a i=1 X i,a . (7)
Note that we are not interested here in controlling the expected regret R T as it gives no information on the risk of the strategy A and of pulling one arm. Indeed, we have the following standard decomposition

R T = T E ν a ⋆ X -E T s=1 Y s = a∈A E ν a ⋆ [X] -E νa [X] E N T,a , (8)
while one would prefer to have a more informative measure, taking into account for instance the variance of the arms or some control of the tails. For this purpose, an other natural notion of regret is the risk-averse regret R T (λ) defined by

R T (λ) = a∈A κ -λ,ν a ⋆ -κ -λ,νa E N T,a . (9)
In the sequel, we control both ( 7) and ( 9) as they both offer interesting interpretation.

The price for risk-aversion

At a high-level, there is obviously a trade-off between trying to get maximal rewards and being risk-averse. Being too cautious (such as, arguably, Exp3 see [START_REF] Auer | The nonstochastic multiarmed bandit problem[END_REF]) avoids getting linear regret, but prevents from getting high rewards as well. On the other hand, simply targeting the maximal mean (such as UCB see [START_REF] Auer | Finite-time analysis of the multiarmed bandit problem[END_REF]) enables to get close to optimal rewards on average, but possibly very bad rewards in difficult environments (when sub-optimal arms have fat lower tails).

A similar situation appears in the standard expected regret setting for the class of UCB-ρ algorithms as shown by [START_REF] Salomon | Robustness of stochastic bandit policies[END_REF]: for, ρ > ρ ′ , UCB-ρ can compete with a larger class of environments than UCB-ρ ′ . However UCB-ρ ′ will beat UCB-ρ on simpler environments.

Simple and complex environments The risk-averse regret (9) captures the sub-optimality of an algorithm in terms of risk-aversion at some fixed level λ. As such, it is the direct equivalent of the expected regret in multi-armed bandits, and we control this regret for our RA-UCB procedure in Theorem 1. If such control may seem satisfactory for many reasons, it has also some drawback. Namely, the level of risk-aversion is here not related in any way to the actual distribution of rewards, since it is some parameter chosen a priori by the practitioner who wants to be protected against sampling possibly very low rewards. As a result, in easy situations when the rewards distributions have very light tails, a high risk-averse algorithm will be too cautious, and will get lower cumulative rewards than a less risk-averse algorithm, such as UCB. Similarly, if the actual distributions have very fat lower tails, a low risk-averse algorithm may not be cautious enough and thus get bad rewards compared to a more risk-averse algorithm, such as Exp3. See also figure 2.3.

Since such situations, that are of immediate practical interest, are not captured by the risk-averse regret (9) defined for some level λ, this motivates the study of the empirical risk-aversion regret [START_REF] Cappé | Kullback-leibler upper confidence bounds for optimal sequential allocation[END_REF] as this one is able to capture such behaviors (this is because it makes appear the empirical rewards coming from the actual distribution explicitly).

Note that this also raises the question of automatically adapting the level of risk-aversion to some bandit problem, or equivalently getting the best of all RA-UCB-λ algorithms (in terms of cumulated reward), which is very hard, (or even impossible, see [START_REF] Salomon | Robustness of stochastic bandit policies[END_REF] for impossibility results regarding UCB-ρ in the related problem of adaptivity in bandit problems). Since this involves orthogonal ideas that would worsen readability and interpretation, add a difficult layer of complexity, and is little justified in practice (where the level of risk-aversion is often simply fixed), we do not study this question in the present work.

Contribution

The difficult situation for risk-aversion appears when the suboptimal arms produce rewards much lower than their mean (heavy lower tail) while the best arm produces rewards much higher than its mean (heavy upper tail): this creates maximal regret. We introduce in section 4 the RA-UCB algorithm that guarantees a low regret in such difficult environments (contrary to e.g. UCB).

For clarity purpose, we clearly separate in two dedicated sections the analysis that is tighten and intrinsic to the risk-aversion problem (Proposition 1 that holds for any algorithm) from the more standard techniques used in stochastic bandits (Theorem 1). We derive on the way some non-trivial concentration results needed for the proof (Lemma 1, Lemma 2, equation ( 14)). The regret of RA-UCB essentially scales as O(log(T )) with the time horizon T up to a distribution-dependent complexity factor.

A Generic Decomposition of the Empirical Regret

We now introduce a generic decomposition of the regret, valid for any strategy A, that is the direct equivalent of (9) for the empirical regret. Proposition 1. Let us define, for some non negative constants {u a } a=1,...,A the event that sub-optimal arms are pulled too much

Ω def = ∃a = a ⋆ : N A T,a > u a ,
and let us fix some value of λ such that κ -λ,νa exists for all a = 1, . . . , A. Then, for all δ ∈ (0, 1), with probability higher than 1 -δ -P(Ω), the regret of the strategy A is upper bounded by

R T (λ) ≤ a =a ⋆ u a κ -λ,ν a ⋆ -κ -λ,νa + m - -λ,ν a ⋆ a =a ⋆ u a + (A -1) log(2A/δ) λ + inf λ ′ >0 m + λ ′ ,ν a ⋆ a =a ⋆ u a + log(2A/δ) λ ′ . ( 10
)
The first term of (10) makes appear a quantity very similar to that of the optimal regret bounds for the expected regret in the stochastic setting, where the standard optimality gaps E ν a ⋆ X -E νa X are replaced by κ -λ,ν a ⋆ -κ -λ,νa , as expected. Now the second and third terms involve the mixability gaps of the optimal arm. The third term is intuitive: indeed, a regret minimizing algorithm will try to understand κ -λ,νa for each arm, and prevent from large deviations below the mean (bad rewards). However, this does not prevent the optimal arm to have large deviations above the mean (that is, unexpected good rewards), which is precisely captured by the third term. Now the presence of the second term comes from another phenomenon: λ is a parameter of the algorithm that tries to pull the arm with highest risk-aversion at level λ. As such, this goal may be successful or not depending on intrinsic properties of the environment. We say that λ is well-adapted to the environment if it is such that the second term in [START_REF] Dembo | Large Deviations Techniques and Applications[END_REF] is negligible before the first term.

So as to provide some intuition, let us now specialize Proposition 1 to the case of Example 1 for illustration purpose. In this case, the mixability gaps of the optimal arm a ⋆ equal λ 2 σ 2 a ⋆ and λ ′ 2 σ 2 a ⋆ , so that if we introduce for convenience the quantity u def = a =a ⋆ u a , one can rewrite [START_REF] Dembo | Large Deviations Techniques and Applications[END_REF] as

R T (λ) ≤ a =a ⋆ u a κ -λ,ν a ⋆ -κ -λ,νa + uλ 2 σ 2 a ⋆ + (A -1) log(A/δ) λ + 2u log(A/δ)σ a ⋆ . ( 11
)
Thus λ is well-adapted to the environment for instance when λ = Ω(u -1/2 ). Since any reasonable algorithm will pull sub-optimal arms only u a = O(log(T )) times with high probability, this indicates that a well-adapted level of risk aversion for a Gaussian game of length T is of order 1 λ = Ω(log(T ) -1/2 ). A similar reasoning holds for the sub-Gaussian and thus the bounded case as well, since we only need an upper-bound on the mixability gaps rather than an equality here. In the sequel, we consider such a case, disregarding the extremely challenging question of defining and estimating a distribution-dependent optimally-adapted value of λ (it also conveys difficult interpretation since the optimal arm depends on λ). Note finally that contrary to the empirical regret, the risk-averse regret ( 9) is completely blind to such situations, as it basically corresponds to the first term in [START_REF] Dembo | Large Deviations Techniques and Applications[END_REF].

Proof. We begin the proof with the following decomposition. For non negative values λ and λ ′ , using the property that T -N A T,a ⋆ = a =a ⋆ N A T,a , one has by a simple rewriting

R T (λ) = T i=N A T ,a ⋆ +1 X i,a ⋆ -κ λ ′ ,ν a ⋆ + a =a ⋆ N A T,a κ λ ′ ,ν a ⋆ -κ -λ,νa + a =a ⋆ N A T ,a i=1 κ -λ,νa -X i,a .
Now, we have on the one hand that R T (λ) ≤ T under the event Ω, while on the other hand, under its complement Ω c , we have

R T (λ) ≤ a =a ⋆ u a κ λ ′ ,ν a ⋆ -κ -λ,νa + max s≤ a =a ⋆ ua T i=T -s+1 X i,a ⋆ -κ λ ′ ,ν a ⋆ + a =a ⋆ max s≤ua s i=1 κ -λ,νa -X i,a .
In this decomposition, the two last terms are controlled by means of concentration of measure, following suitable adaptations of ( 1) and ( 2), and the first term can be further decomposed due to the inequality

κ λ ′ ,ν a ⋆ -κ -λ,νa ≤ m + λ ′ ,ν a ⋆ + m - -λ,ν a ⋆ + κ -λ,ν a ⋆ -κ -λ,νa
, where we made appear the mixability gaps of the optimal arm a ⋆ . More precisely, the generic bound on the regret of a learning algorithm A now relies on the previous decomposition and on the following two concentration results, whose proof is in the Appendix. Lemma 1. Let τ ⋆ ∈ N be some positive constant, and λ ′ > 0 be such that κ λ ′ ,ν a ⋆ is finite. Then we have the property that for all ǫ > 0, then

P max s≤τ ⋆ T i=T -s+1 X i,a ⋆ -κ λ ′ ,ν a ⋆ ≥ ǫ ≤ exp(-λ ′ ǫ)
Lemma 2. Let τ ∈ N be some positive constant, and λ > 0 be such that κ λ,νa is finite. Then we have the property that for all ǫ > 0, then

P max s≤τ s i=1 κ -λ,νa -X i,a ≥ ǫ ≤ exp(-λǫ)
In order to conclude the proof of Proposition 1, we then apply Lemma 1 to the value τ ⋆ = a =a ⋆ u a with ǫ = log(2A/δ) λ ′ and then Lemma 2 to τ = u a with ǫ = log(2A/δ) λ . We deduce, by a union bound, that with probability higher than 1 -δ -P(Ω), then

R T (λ) ≤ a =a ⋆ u a κ -λ,ν a ⋆ -κ -λ,νa + a =a ⋆ log(2A/δ) λ + inf λ ′ >0 a =a ⋆ u a m + λ ′ ,ν a ⋆ + m - -λ,ν a ⋆ + log(2A/δ) λ ′ .

The Risk-Averse Upper Confidence Bound algorithm

We introduce in this section a strategy A that we call the RA-UCB algorithm.

From now on, we restrict to the case when all distributions belong to M + 1 (R B ), where R B = (-∞, B] for some known value of B. Thus, let us introduce for all a ∈ A, the empirical distribution νt (a) ∈ M + 1 (R B ) associated to ν a , built using the past observations Y 1 , . . . , Y t ; we define νt (a) Further, for clarity purpose, we now use the notation νn,a (with a in subscript) in order to denote the empirical distribution built from the n first samples drawn from ν a , while we reserve the functional notation νt (a) for the empirical distribution built from the samples received from arm a up to time t. Naturally, we have that νt (a) = νN A t (a),a . More generally, for some distribution ν, we also write νn for its empirical distribution built from n samples.

The RA-UCB algorithm is inspired from the strategies introduced by [START_REF] Lai | Asymptotically efficient adaptive allocation rules[END_REF][START_REF] Burnetas | Optimal adaptive policies for sequential allocation problems[END_REF][START_REF] Maillard | A finite-time analysis of multi-armed bandits problems with Kullback-Leibler divergences[END_REF][START_REF] Garivier | The KL-UCB algorithm for bounded stochastic bandits and beyond[END_REF][START_REF] Cappé | Kullback-leibler upper confidence bounds for optimal sequential allocation[END_REF] as it selects at time t + 1 the arm A t+1 = argmax a∈A U t (a), where U t (a) is an upper confidence bound on the risk aversion of arm a at level λ, defined by

U t (a) def = sup κ -λ,ν : K(ν t (a), κ -λ,ν ) ≤ f (t) N t,a , (12)
and where we introduced the following quantity

K(ν t (a), r) def = inf KL(ν t (a)||ν) : ν ∈ M + 1 (R B ), κ -λ,ν ≥ r . (13)
Note that UCB-like algorithms are unnatural in this setting: they are based on empirical means only, while we really need to control the tail distributions here. KL-based algorithm are more suitable, and produce much stronger results. Note also that the parameter λ is here the same that defines the level of risk aversion used in the definition of the regret. The algorithm requires another parameter, that is a non-decreasing function of the time f . A typical choice is such that f (t) = O(log(t)), as mentioned in Theorem 1.

A Useful Formulation with Dual Optimality Conditions

The definition of the bound (12) may seem quite abstract. In order to make it more computable and explicit, we now provide the following result, that is a dual formulation of the optimization problem appearing in the definition of K(ν t (a), r) (see the proof in the appendix). Lemma 3. Let νn denote with a finite number n of atoms. Then the following dual formulation holds

K(ν n , r) = max 1 n n i=1 log 1 - γ ⋆ λ 1 -e -λ(xi-r) : 0 ≤ γ ⋆ ≤ λ 1 -e -λ(B-r) .
This result shows that the optimization problem [START_REF] Even-Dar | Risk-sensitive online learning[END_REF] can actually be solved numerically. and is deeply linked to the numerically efficient dual formulation considered for instance in [START_REF] Borwein | Duality relationships for entropy-like minimization problem[END_REF], [START_REF] Harari-Kermadec | Vraisemblance empirique généralisée et estimation semiparamétrique[END_REF], or re-derived more recently in [START_REF] Honda | An asymptotically optimal bandit algorithm for bounded support models[END_REF] for the related problem of optimal regret bounds in the stochastic multi-armed bandit with expected regret criterion. For completeness, it makes sense to introduce the quantity for distributions

ν ∈ M + 1 (R B ) K(ν, r) = sup E log 1 - γ ⋆ λ 1 -e -λ(X-r) : 0 ≤ γ ⋆ ≤ λ 1 -e -λ(B-r) .

Regret Analysis of the RA-UCB Algorithm

By the generic decomposition result of Proposition 1, we only have to provide a high-probability upper bound on the number of pulls of any sub-optimal arm a, more precisely on the event

Ω def = ∃a = a ⋆ : N A T,a > u a ,
In order to control the probability of such on event, let us introduce, for all a = a ⋆ , the random time t a corresponding to the last round when a is chosen, that is we have N A ta,a = N A T,a -1 and N A ta+1,a = N A T,a . For such a t a , we also have by definition

A A ta+1 = a.
Decomposition of Events (step 1) We start by considering the event A A t+1 = a for a sub-optimal arm. By definition of the algorithm, we have the property that U t (a) ≥ U t (⋆). This event can be decomposed as

A A t+1 = a ⊂ U t (⋆) ≤ κ ⋆ ∪ U t (⋆) > κ ⋆ and A A t+1 = a ⊂ U t (⋆) ≤ κ ⋆ ∪ U t (a) > κ ⋆ and A A t+1 = a ⊂ K(ν t (⋆), κ ⋆ ) ≥ f (t) N A t,⋆ ∪ K(ν t (a), κ ⋆ ) ≤ f (t) N A t,a )
and A A t+1 = a , where we introduced here some quantity κ ⋆ . We now make use of this decomposition in order to show that, for all choice of constant {u a } a=1,...,A with u a > 1, we have

P ∃a = a ⋆ : N A T,a > u a ≤ P ∃a = a ⋆ : N A ta,a > u a -1 and A ta+1 = a ≤ P ∃a = a ⋆ : N A ta,a ≥ u a and K(ν ta (⋆), κ ⋆ ) ≥ f (t a ) N A ta,⋆ +P ∃a = a ⋆ : N A ta,a ≥ u a and K(ν ta (a), κ ⋆ ) ≤ f (t a ) N A ta,a and A A ta+1 = a ≤ P ∃n ≤ T : K(ν n,⋆ , κ ⋆ ) ≥ f (n) n + a =a ⋆ T n=ua P K(ν n,a , κ ⋆ ) ≤ f (T ) n ,
where we used in the last line that f (t a ) ≥ f (N A ta,⋆ ) for the first term, since f is non-decreasing, and similarly that f (t a ) ≤ f (T ) together with a union bound for the second term.

Concentration Inequalities (step 2)

We now make use of concentration inequalities. More precisely, we first use that for the optimal arm it holds that for the value κ ⋆ = κ -λ,ν a ⋆ and for all ǫ > 0, then

P K(ν n,⋆ , κ ⋆ ) ≥ ǫ ≤ e(n + 2) exp(-nǫ) , (14)
(the proof of which is provided in the Appendix, Proposition 2) and that for any suboptimal arm a, for all ǫ > 0 one can resort to an application of nonasymptotic Sanov's lemma. Indeed, under some conditions that we detail below, an easy consequence of [10, Exercise 2.2.38] is that

P K(ν n,a , κ ⋆ ) ≤ ǫ ≤ exp -nχ a (ǫ) , (15)
where the quantity χ a (ǫ) is intrinsic to the complexity of the problem of testing whether a is a good arm and is defined by

χ a ǫ def = inf KL(ν||ν a ) : K(ν, κ ⋆ ) ≤ ǫ . ( 16
)
Note that this term satisfies that χ a K(ν a , κ ⋆ ) = 0 and that for all ǫ < K(ν a , κ ⋆ ) then χ a ǫ > 0. In particular, the inequality ( 15) is non trivial only for such ǫ < K(ν a , κ ⋆ ).

Important remark Note at this point that there may be some topological difficulty in order to meet the conditions needed for [START_REF] Honda | An asymptotically optimal bandit algorithm for bounded support models[END_REF]. As a reminder, these conditions are that χ a (ǫ) < ∞, that the set {ν : K(ν, κ ⋆ ) ≤ ǫ} is convex and moreover that it is close (for the weak topology). The first condition is not restrictive (the case when it is not met is actually even a favorable situation). Now, the convexity of the considered set easily follows from the biconvexity of the Kullback-Leibler divergence (see [START_REF] Cover | Elements of Information Theory[END_REF]) and of the convexity of the set of distributions with high risk-aversion. The latter follows from the concavity of the log function. Now, in order to show that this set is close, it is sufficient to show that K(•, κ ⋆ ) is lower semi-continuous on a set including the limit distribution ν a , which is also not difficult to prove (following [START_REF] Honda | An asymptotically optimal bandit algorithm for bounded support models[END_REF] for instance).

Final Control of the Number of Pulls (step 3) So far, by combining the initial decomposition with ( 14) and ( 15) together with a union bound, we have shown

P ∃a = a ⋆ : N A T,a > u a ≤ T n=1 e(n + 2)e -f (n) + a =a ⋆ T n=ua exp -nχ a f (T ) n .
The first sum in the right hand side of this inequality is easily controlled. For instance if the parameter function f satisfies that f (n) ≥ log(2e(n + 2)n 2 /δ), then it is less than δ. The last sum can be made more explicit. Let us now define the quantity u a = (1 + ǫ a ) f (T ) K(νa,κ ⋆ ) + 1 for all a = a ⋆ and some ǫ a > 0, and further define ǫ a by

ǫ a = inf ǫ > 0 : 1 + ǫ K(ν a , κ ⋆ ) χ a K(ν a , κ ⋆ ) 1 + ǫ ≥ 1 . (17)
Note that such an ǫ a exists and is finite since χ a is a non increasing function, so that 1+ǫ

K(νa,κ ⋆ ) χ a K(νa,κ ⋆ ) 1+ǫ
is an increasing function of ǫ. With such notation, we deduce

T n=ua exp -nχ a f (T ) n ≤ exp -u a χ a f (T ) u a 1 + ∞ n=1 exp -nχ a f (T ) u a ≤ exp -f (T ) 1 + ǫ a K(ν a , κ ⋆ ) χ a K(ν a , κ ⋆ ) 1 + ǫ a 1 + ∞ n=1 exp -nχ a K(ν a , κ ⋆ ) 1 + ǫ a ≤ δ 2e(T + 2)T 2 1 + 1 1 -exp -χ a K(νa,κ ⋆ ) 1+ǫa
.

where we used in the first line that n > u a and that χ a is non increasing, in the second line that

u a f (T ) = (1 + ǫ a ) 1 K(ν a , κ ⋆ ) + 1 f (T ) > 1 + ǫ a K(ν a , κ ⋆ ) ,
and finally in the third line the bound on f (T ) and the definition of ǫ a .

The previous analysis together with the result of Proposition 1 for the empirical risk-averse regret [START_REF] Cappé | Kullback-leibler upper confidence bounds for optimal sequential allocation[END_REF] on the one hand, and the definition of the risk-averse regret (9) on the other hand enable us to deduce the following bound on the regret of the RA-UCB algorithm, which is the main result of this paper. We further provide Corollary 1 for illustration purpose.

Theorem 1. Assume that for all a = 1, . . . , A, then ν a ∈ M + 1 (R B ), and define ǫ a by equation [START_REF] Honda | Finite-time regret bound of a bandit algorithm for the semi-bounded support model[END_REF]. Let us define f (t) = log(2e(t + e)t 2 /γ) for some γ ∈ (0, 1). Then, the risk-averse regret of the RA-UCB algorithm is upper bounded at time T by

R T (λ) ≤ a =a ⋆ (1 + ǫ a )∆ a K a log 2e(T + e)T 2 /γ + a =a ⋆ (1 + ǫ a )∆ a + a =a ⋆ γ∆ a T + γ∆ a 2e(T + 2)T 1 + 1 1 -exp -χ a Ka 1+ǫa
.

Further, the empirical regret of the RA-UCB algorithm at time T is upper bounded for all δ ∈ [0, 1] by

R T (λ) ≤ a =a ⋆ (1 + ǫ a )∆ a K a log 2e(T + e)T 2 /γ + a =a ⋆ (1 + ǫ a )∆ a + a =a ⋆ log(2A/δ) λ + inf λ ′ >0 a =a ⋆ (1 + ǫ a ) ∆ a K a log 2e(T + e)T 2 /γ + 1 m + λ ′ ,ν a ⋆ + m - -λ,ν a ⋆ + log(2A/δ) λ ′ ,
where we introduced the optimality gaps ∆ a = κ -λ,ν a ⋆ -κ -λ,νa , and K a = K(ν a , κ -λ,ν a ⋆ ), with probability higher than

1 -δ -γ - γ 2e(T + 2)T 2 a =a ⋆ 1 + 1 1 -exp -χ a Ka 1+ǫa
.

Note that the quantity 1 -exp -χ a Ka 1+ǫa

appearing in the high probability bound is problem dependent but is actually a constant, and that similarly ǫ a is also a problem dependent constant. In particular, both quantities are independent on T and on the algorithm. Taking into account these remarks together with the discussion after Proposition 1 regarding well-adapted values of λ leads to the more readable corollary Corollary 1. Under the same assumptions as Theorem 1, for γ = Θ(T -1 ), then

R T (λ) ≤ 5 a =a ⋆ (1 + ǫ a )∆ a K a log T + O 1 . (18)
Further, assuming moreover that the distributions of rewards are all sub-Gaussian, and that the level of risk-aversion is λ = Θ(log(T ) -1/2 ), then for a choice of γ = Θ(T -β ) for some β ≥ 0 the empirical regret of the RA-UCB algorithm at time T is bounded as

R T (λ) ≤ c a =a ⋆ (1 + ǫ a )∆ a K a log T + O log(T ) , (19)
with probability of order 1 -δ -o(1) for some constant c ≤ 4 + β.

Note also that assuming λ = Ω(log(T ) -1/2 ) only makes the second to last terms in Theorem 1 o(log(T )) instead of O(log(T )). Thus even for others value of λ (that is a constant), we still have O(log(T )) empirical regret (with possibly larger constants).

Discussion

The bound [START_REF] Kaufmann | Thompson sampling: An asymptotically optimal finite-time analysis[END_REF] makes appear a first order term scaling with a log(T ), which looks very much like the results for the standard multi-armed bandit with expected regret. Note that that such a dependency is achievable is not obvious since working with risk-aversion is usually considered as much more difficult than working with expectation. This should also be compared to the result of [START_REF] Sani | Risk-aversion in multi-armed bandits[END_REF], although they consider a different, trickier to interpret, setting. We show however that this is possible. The constant before the logarithmic term consists of the ratio ∆a Ka that is also very similar to the known bounds for the expected regret ( [START_REF] Burnetas | Optimal adaptive policies for sequential allocation problems[END_REF][START_REF] Maillard | A finite-time analysis of multi-armed bandits problems with Kullback-Leibler divergences[END_REF]), up to the constant c, that could definitely be reduced by a more careful analysis and parameter tuning (this is not the main focus of this work), and more importantly the constant 1 + ǫ a . Theorem 1 holds for a larger class of distributions than the one considered e.g. in [START_REF] Maillard | A finite-time analysis of multi-armed bandits problems with Kullback-Leibler divergences[END_REF]. The reason for this is precisely because we accept to loose the constant 1 + ǫ a (as opposed to 1 in their work). This term is not entirely intrinsic: one could easily change the threshold 1 in the definition [START_REF] Honda | Finite-time regret bound of a bandit algorithm for the semi-bounded support model[END_REF] to a smaller constant at the price of an increased probability term (in the o(1)), so that one may a priori optimize this term further. It also more complex than the quantities ∆ a and K a whose interpretation is immediate. However, understanding the function x → xχ a (1/x), which is related, by the definition of χ a , to understanding how K(ν, κ -λ,ν a ⋆ ) varies when ν moves from ν a to distributions having higher risk-aversion, is definitely needed, and can not be done simply using ∆ a and K a (one could use for instance the derivative of K(•, κ -λ,ν a ⋆ ) in the direction of high risk-aversion, but this is not a lot more easier to interpret either and still difficult to handle in specific cases). Thus, we here prefer to let the term 1 + ǫ a as it is, since it is anyway fully explicit, and captures such distribution-dependent behavior in a fairly concise way. Note also that a similar difficulty appears in the expected regret setting (see [START_REF] Honda | An asymptotically optimal policy for finite support models in the multiarmed bandit problem[END_REF], [START_REF] Maillard | A finite-time analysis of multi-armed bandits problems with Kullback-Leibler divergences[END_REF], [START_REF] Garivier | The KL-UCB algorithm for bounded stochastic bandits and beyond[END_REF]).

Conclusion

In this work, the variant of the stochastic multi-armed bandit problem when one considers looking for the maximally risk-averse arm for a user-defined level of risk-aversion (instead of the mean) is considered. We first provide a generic decomposition of the regret (Proposition 1) for any algorithm that enables to focus on the number of plays of sub-optimal arms only. We make use of a coherent risk measure based on the cumulant generative function and show that it is possible to achieve optimal performance up to a regret that is logarithmic (Corollary 1) in the time horizon (with distribution dependent constants). This logarithmic regret is achieved by some adaptation of existing algorithms designed for the expected regret, together with new concentration results precisely introduced for the control of risk-aversion rather than of the mean, which are of independent interest.

A Generic Decomposition of the Regret

In this section, we provide a proof for Lemma 1 and for Lemma 2.

Lemma 1 Let τ ⋆ ∈ N be some positive constant, and λ ′ > 0 be such that κ λ ′ ,ν a ⋆ is finite. Then we have the property that for all ǫ > 0, then

P max s≤τ ⋆ T i=T -s+1 X i,a ⋆ -κ λ ′ ,ν a ⋆ ≥ ǫ ≤ exp(-λ ′ ǫ)
Proof. The proof follows by an application of Doob's inequality for submartingales. Let us introduce the quantity M s = T i=T -s+1 X i,a ⋆ -κ λ ′ ,ν a ⋆ , which is a sum of s iid random variables. Note that M s is not centered, and is a supermartingale. Indeed, it satisfies that

E M s+1 M 1 , .., M s = M s + E X T -s,a ⋆ -κ λ ′ ,ν a ⋆ ≤ M s
On the other hand, we have the property that for suitable values of γ > 0, then exp(γM s ) is a submartingale. Indeed, we have the property that

log E exp(γM s+1 ) M 1 , .., M s = γM s + log E exp γX T -s,a ⋆ - γ λ ′ log E exp(λ ′ X) = γM s + γ 1 γ log E exp(γX) - 1 λ ′ log E exp(λ ′ X) = γM s + γ κ γ,ν a ⋆ -κ λ ′ ,ν a ⋆ .
Thus, for all γ > 0 such that κ γ,ν a ⋆ ≥ κ λ ′ ,ν a ⋆ , then exp(γM s ) is a submartingale, and Doob's maximal inequality applies. This shows that for all ǫ > 0, then

P max s≤τ ⋆ M s ≥ ǫ = P max s≤τ ⋆ exp(γM s ) ≥ exp(γǫ) ≤ E exp(γM s ) exp(-γǫ) = E exp(γX) τ ⋆ E exp(λ ′ X) -γ λ ′ τ ⋆ exp(-γǫ) = exp γτ ⋆ κ γ,ν a ⋆ -κ λ ′ ,ν a ⋆ exp(-γǫ) .
We conclude by simply taking γ = λ ′ (though this is a suboptimal choice).

Lemma 2 Let τ ∈ N be some positive constant, and λ > 0 be such that κ λ,νa is finite. Then we have the property that for all ǫ > 0, then

P max s≤τ s i=1 κ -λ,νa -X i,a ≥ ǫ ≤ exp(-λǫ)
Proof. The proof is very similar do that of Lemma 1. We have the property that for s = 1, . . . , τ , M s = s i=1 κ -λ,νa -X i,a is a supermartingale since κ -λ,νa ≤ E νa X , and then that exp(γM s ) is a submartingale for values of γ such that κ -λ,νa ≥ κ -γ,νa .

B Dual Formulation of the Minimization Algorithm

Lemma 3 Let νn denote the empirical distribution build from ν ∈ M + 1 (R B ) with n i.i.d samples. We have the following rewriting

K(ν n , r) = max 1 n n i=1 log 1- γ ⋆ λ 1-exp -λ(x i -r) : 0 ≤ γ ⋆ ≤ λ 1 -exp(-λ(B -r))
.

Proof. Let x 1 < • • • < x n denote the support of the empirical distribution νn and note that the optimization ( 13) can be reduced to distributions with support in {x i } i≤n up maybe to one extra point x n+1 that will receive extra weight if needed. For ν ≪ νn , let p i be the probability weight associated to x i , and let q be the one associated to x n+1 . The optimization problem can be rewritten in the following form minimize over {p i } i≤n , q :

n i=1 1 n log 1/n p i subject to - 1 λ log n i=1 exp(-λx i )p i + exp(-λx n+1 )q ≥ r ,
and it is immediate to see that x n+1 < x n does not help. Now by introducing appropriate Lagrange-multipliers, this corresponds to minimizing the following quantity over {p i } i≤n , q ∈ R + , γ, α, {ξ i } i≤n , ξ ∈ R:

V def = n i=1 1 n log 1/n p i + γ r + 1 λ log n i=1 exp(-λx i )p i + exp(-λx n+1 )q +α 1 - n i=1 p i -q - n i=1 ξ i p i -ξq ,
with Karush-Kuhn Tucker optimality conditions corresponding to a distribution ν ⋆ being:

-γ ⋆ r + 1 λ log n i=1 exp(-λx i )p ⋆ i + exp(-λx n+1 )q ⋆ = 0, with the condi- tions γ ⋆ ≥ 0 and -1 λ log n i=1 exp(-λx i )p ⋆ i + exp(-λx n+1 )q ⋆ ≥ r , - n i=1 p ⋆ i + q ⋆ = 1 , -ξ ⋆ i p ⋆ i = 0, ξ ⋆ i ≥ 0, p ⋆ i ≥ 0 and ξ ⋆ q ⋆ = 0, ξ ⋆ ≥ 0, q ⋆ ≥ 0 , --1 np ⋆ i + γ ⋆ λ n j =i=1 exp(λ(x i -x j ))p ⋆ j +p ⋆ i +exp(λ(x i -x n+1 ))q ⋆ -1 -ξ ⋆ i -α ⋆ = 0 , -γ ⋆ λ n j=1 exp(λ(x n+1 -x j ))p ⋆ j + q ⋆ -1 -ξ ⋆ -α ⋆ = 0 .
Now, in order for K(ν n , ν ⋆ ) to be finite, we need p ⋆ i > 0 i.e. ξ ⋆ i = 0 for all i ≤ n. If γ ⋆ = 0, then we deduce that α ⋆ = -1 np ⋆ i = -ξ ⋆ for all i ≤ n. In particular, this means that ξ ⋆ > 0 and thus that q ⋆ = 0. Thus, since p ⋆ i = 1 nξ ⋆ and that n i=1 p ⋆ i = 1, then ⋆ = 1 and p ⋆ i = 1 n for all i ≤ n. This correspond to the cases when V = 0. Now let us consider that γ ⋆ > 0. In that case, we deduce that we must have for all i ≤ n exp(λ(x i -r)) = n j =i=1 exp(λ(x i -x j ))p ⋆ j + p ⋆ i + exp(λ(x i -x n+1 ))q ⋆ and also exp(λ(x n+1 -r)) = n j=1 exp(λ(x n+1 -x j ))p ⋆ j + q ⋆ , from which we deduce that we must have

--1 np ⋆ i + γ ⋆ λ exp(-λ(x i -r)) -α ⋆ = 0 , -γ ⋆
λ exp(-λ(x n+1 -r)) -ξ ⋆ -α ⋆ = 0 . Thus, on the one hand we have p If q ⋆ > 0, then ξ ⋆ = 0 and thus α ⋆ = γ ⋆ λ exp(-λ(x n+1 -r)). Thus, this entails that exp(λ(x n+1 -x i )) exp(-λ(x i -r)) -exp(-λ(x n+1 -r)) ,

i.e. that the optimal value of γ ⋆ is given by

γ ⋆ = λ 1 -exp(λ(x n+1 -r)) 1 n n i=1 1 -exp(-λ(x i -x n+1 )) exp(-λ(x i -r)) -exp(-λ(x n+1 -r)) = λ 1 -exp(-λ(x n+1 -r))
.

We thus have shown that the optimal weights are given by

p ⋆ i = 1 n
1 -exp(-λ(x n+1 -r)) exp(-λ(x i -r)) -exp(-λ(x n+1 -r)) ,

Proof. Using the change of variable γ = This follows from the fact that for all u ∈ [0, 1] there exists u ′ ∈ U c such that for all x ≤ x n+1 , then

log 1 -u 1 -exp -λ(x -r) 1 -exp -λ(x n+1 -r) ≤ c + log 1 -u ′ 1 -exp -λ(x -r) 1 -exp -λ(x n+1 -r) . ( 20
)
More precisely we use the fact that for all u, u ′ ∈ [0, 1] such that u ≤ u ′ ≤ 1/2 or u ≥ u ′ ≥ 1/2, then for all y ≤ 1 log 1 -uy ≤ log 1 -u ′ y + 2|u ′ -u| .

The proof of (20) then follows from simple algebra. We then apply this result with the following grid that has at most 2 + 1/c elements U c def = {1/2, 1} ∪ {1/2 + c, . . . , 1/2 + ⌊1/(2c)⌋c} ∪ {1/2 -c, . . . , 1/2 -⌊1/(2c)⌋c} .

We then prove the following proposition, that is the main result of this section.

Proposition 2. Let νn be the empirical distribution built from ν with n i.i.d samples. Then, for all ǫ > 0, it holds that P K(ν n , κ -λ,ν ) ≥ ǫ ≤ e(n + 2) exp(-nǫ) .

Proof. Let r = κ ν . By Lemma C, we know that for all c > 0, there exist a finite set Γ c ⊂ [0, γ ⋆ ≤

Fig. 1 .

 1 Fig.1. Plot of arms' densities and their mean: left) an environment where no arm has fat lower tail. right) an environment where for some λ, the best arm (green) does not have best mean, and sub-optimal arms (red, blue) have fat lower tails.

δ

  Ys I{A s = a} where N A

  exp(-λ(xi-r))-α ⋆ and q ⋆ = 1exp(-λ(xi-r))-α ⋆ , and on the other hand, we have q ⋆ = exp(λ(x n+1 -r))-1 n n i=1 exp(λ(xn+1-xi)) γ ⋆ λ exp(-λ(xi-r))-α ⋆ ,and it remains to determine α ⋆ and γ ⋆ .

  λ(x i -r)) -exp(-λ(x n+1 -r)) = exp(λ(x n+1 -r)) -

λ 1 - 1 -u 1 -

 111 exp(-λ(xn+1-r)) u, with x n+1 = B this corresponds exp -λ(x i -r) 1 -exp -λ(x n+1 -r) ≤ c + max exp -λ(x i -r) 1 -exp -λ(x n+1 -r) .

λ 1 - 1 - 1 - 1 -≤ log exp βc γ∈Γc n i=1 E ν 1 - γ λ 1 -E ν 1 - γ λ 1 -

 1111i=111 exp(-λ(B-r)) ] with at most 2 + 1/c points such thatK(ν n , r) ≤ c + max γ∈Γc exp -λ(x i -r).Then, we thus have for anyβ > 0 log E exp βK(ν n , r) ≤ log E exp βc + β max γ∈Γc exp -λ(x i -r) exp -λ(x i -r) exp -λ(x i -r) β n.Thus, for β = n, we get, since r = -1 λ log E ν exp -λX , exp -λ(X -r) exp -λX exp(λr) = 1 , from which we deduce thatlog E exp βK(ν n , r) ≤ inf c>0 log exp nc |Γ c | .In particular, for c = 1/n, we get the boundlog E exp βK(ν n , r) ≤ log e|Γ 1/n | ≤ log e(n + 2) ,and the result then follows via a simple Markov's inequality.

Such (weak) dependency with T is intuitive: if we only have 10 trials do to something, we would be much more risk-averse (big λ) than with 1000 trials.

and that the corresponding value is

. On the other hand, we must have q ⋆ = 0 and thus

Thus combining these two equations we deduce that γ ⋆ λ -α ⋆ = 1 and thus that

The value is thus given in that case by

.

C Concentration Inequalities for K(ν n , r)

Let us first show the following key intermediate lemma

Lemma

We have the property that there exists some discrete set Γ c with at most 2 + 1/c elements such that


Emmanuel Caruyer
email: emmanuel.caruyer@inria.fr.

Christophe Lenglet

Guillermo Sapiro

Rachid Deriche

Design of multi-shell sampling schemes with uniform coverage in diffusion MRI

Purpose In diffusion MRI, a technique known as Diffusion Spectrum Imaging reconstructs the propagator with a discrete Fourier transform, from a Cartesian sampling of the diffusion signal. Alternatively, it is possible to directly reconstruct the orientation distribution function in Q-ball imaging, providing so-called high angular resolution diffusion imaging. In between these two techniques, acquisitions on several spheres in Q-space offer an interesting trade-off between the angular resolution and the radial information gathered in diffusion MRI. A careful design is central in the success of multi-shell acquisition and reconstruction techniques.

Methods The design of acquisition in multi-shell is still an open and active field of research however. In this work, we provide a general method to design multi-shell acquisition with uniform angular coverage. This method is based on a generalization of electrostatic repulsion to multi-shell.

Results

We evaluate the impact of our method using simulations, on the angular resolution in one and two bundles of fibers configurations. Compared to more commonly used radial sampling, we show that our method improves the angular resolution, as well as fiber crossing discrimination.

Discussion

We propose a novel method to design sampling schemes with optimal angular coverage, and show the positive impact on angular resolution in diffusion MRI.

Introduction

Diffusion MRI is an imaging technique discovered in the mid 1980's [START_REF] Bihan | Imagerie de diffusion it in vivo par rsonnance magntique nuclaire[END_REF][START_REF] Taylor | The spatial mapping of translational diffusion coefficients by the nmr imaging technique[END_REF][START_REF] Merboldt | Self-diffusion nmr imaging using stimulated echoes[END_REF] that maps the local displacement probability of water molecules within biological tissues. This has provided a great tool over the past 20 years for the early diagnosis of brain disorders and the study of brain structural connectivity, and a precious help for pre-operative planning in neurosurgery [START_REF] Jones | Diffusion MRI: Theory, Methods, and Applications[END_REF][START_REF] Johansen-Berg | Diffusion MRI: From quantitative measurement to in-vivo neuroanatomy[END_REF]. Using a pulsedgradient spin-echo sequence [START_REF] Stejskal | Spin diffusion measurements: spin echoes in the presence of a timedependent field gradient[END_REF], the MR echo magnitude is attenuated due to proton diffusion, and by repeating acquisitions with various diffusion encoding gradients, one can characterize and recover information on the water molecules displacement. This in turns reflect the organization of the underlying tissue microstructure, which restricts the displacement of water molecules.

Formally, the diffusion of water molecules is well described by the ensemble average propagator (EAP), P (r; τ ), which is the full 3D probability density function of water molecules displacement during a given observation time τ . Under the narrow-pulse assumption, which is verified whenever the diffusion time τ is much larger than the pulse duration δ, the signal attenuation, E(q), is related to the EAP through a Fourier transform [START_REF] Callaghan | Principles of Nuclear Magnetic Resonance Microscopy[END_REF]:

E(q) = R 3 P (r)e 2ıπq • r dr, (1)
where the wavevector is q = γδG/2π, γ is the gyromagnetic ratio, and G is the magnetic field gradient.

A technique, known as q-space imaging (QSI) or diffusion spectrum imaging (DSI) [START_REF] Callaghan | Principles of Nuclear Magnetic Resonance Microscopy[END_REF][START_REF] Wedeen | Mapping complex tissue architecture with diffusion spectrum magnetic resonance imaging[END_REF], consists in sampling the diffusion weighted MR attenuation E(q) on a Cartesian grid in the reciprocal space. The EAP is then reconstructed through a discrete Fourier transform, and other quantities of interest can be computed such as the orientation distribution function (ODF), ψ(u), which is the angular marginal probability density of water molecules motion. This technique has been applied successfully for the in vivo imaging of human brain white matter structures. However, the long acquisition time required for such technique limits its application in a clinical context. Several approaches [START_REF] Wu | Hybrid diffusion imaging[END_REF][START_REF] Khachaturian | Boosting the sampling efficiency of q-ball imaging using multiple wavevector fusion[END_REF][START_REF] Aganj | Multiple q-shell ODF reconstruction in q-ball imaging[END_REF][START_REF] Descôteaux | Diffusion propagator imaging: Using laplace's equation and multiple shell acquisitions to reconstruct the diffusion propagator[END_REF][START_REF] Assemlal | Efficient and robust computation of pdf features from diffusion MR signal[END_REF] use a hybrid sampling method, where the diffusion attenuation E(q)

is sampled on several spheres in the q-space, rather than on a Cartesian lattice. Multiple q-shell sampling offers an interesting trade-off between classical high angular resolution diffusion imaging (HARDI) techniques in q-ball imaging [START_REF] Tuch | Q-ball imaging[END_REF][START_REF] Descoteaux | Regularized, fast, and robust analytical q-ball imaging[END_REF], and DSI.

The design of sampling protocols in single shell has received much attention from the community [START_REF] Jones | Optimal strategies for measuring diffusion in anisotropic systems by magnetic resonance imaging[END_REF][START_REF] Papadakis | Minimal gradient encoding for robust estimation of diffusion anisotropy[END_REF][START_REF] Jansons | Persistent angular structure: new insights fom diffusion magnetic resonance imaging data[END_REF]. However, the design of multi-shell sampling is still an active subject of research. There are many degrees of freedom in the construction of a multi-shell design. If we try to determine all the parameters that fully describe a sampling protocol on several shells, we have: the total number of samples, K; the number of spheres in the q-space, S; the shell radii (related to the b-values), q s ; the number of points on each shell, K s , such that s≤S K s = K, and the directions on each sampling shell, u s,k ∈ S 2 . With the large and increasing number of reconstruction techniques in the literature [START_REF] Khachaturian | Boosting the sampling efficiency of q-ball imaging using multiple wavevector fusion[END_REF][START_REF] Wu | Hybrid diffusion imaging[END_REF][START_REF] Assemlal | Efficient and robust computation of pdf features from diffusion MR signal[END_REF][START_REF] Aganj | Reconstruction of the ODF in single and multiple shell q-ball imaging within constant solid angle[END_REF][START_REF] Descoteaux | Multiple q-shell diffusion propagator imaging[END_REF][START_REF] Ye | Efficient interlaced multi-shell sampling scheme for reconstruction of diffusion propagators[END_REF][START_REF] Caruyer | A computational framework for experimental design in diffusion MRI[END_REF], it is not clear whether a preferred, unified strategy would fit the needs of every method. In particular we believe that the selection of parameters K, S, q s and K s should be driven by the needs of the reconstruction method, the time limit for the acquisition and the physical limitations of the imaging system. Still, the choice of sampling directions should follow the same unified strategy, commonly adopted in q-ball imaging [START_REF] Jones | Optimal strategies for measuring diffusion in anisotropic systems by magnetic resonance imaging[END_REF][START_REF] Papadakis | Minimal gradient encoding for robust estimation of diffusion anisotropy[END_REF][START_REF] Jansons | Persistent angular structure: new insights fom diffusion magnetic resonance imaging data[END_REF], in order to achieve an angular coverage as uniform as possible. As for the choice of shell radii, number of shells and number of points per shell, we will follow the strategy recently proposed in [START_REF] Caruyer | A computational framework for experimental design in diffusion MRI[END_REF], which is optimized for reconstruction in the continuous and orthonormal basis called modified spherical polar Fourier (mSPF) [START_REF] Caruyer | Diffusion MRI signal reconstruction with continuity constraint and optimal regularization[END_REF].

In particular, recent studies have shown that for a better angular resolution, the sampling directions should be as different as possible from one shell to another [START_REF] Zhan | Differential information content in staggered multiple shell HARDI measured by the tensor distribution function[END_REF][START_REF] Ye | Efficient interlaced multi-shell sampling scheme for reconstruction of diffusion propagators[END_REF][START_REF] Koay | Sparse and optimal acquisition design for diffusion mri and beyond[END_REF]. The constructions in [START_REF] Zhan | Differential information content in staggered multiple shell HARDI measured by the tensor distribution function[END_REF] and [START_REF] Koay | Sparse and optimal acquisition design for diffusion mri and beyond[END_REF] start from a uniform set of directions on the unit sphere, which is then separated into S subsets of quasi-uniform directions on the sphere, that correspond to the S spheres under consideration for sampling. This method leads to a discrete search problem, that is solved at the cost of a timely optimization procedure in [START_REF] Koay | Sparse and optimal acquisition design for diffusion mri and beyond[END_REF], or at the cost of a sub-optimal greedy approach in [START_REF] Zhan | Differential information content in staggered multiple shell HARDI measured by the tensor distribution function[END_REF]. Another related method, presented in [START_REF] Santis | Improved precision in the charmed model of white matter through sampling scheme optimization and model parsimony testing[END_REF] and, in deeper details, in a private communication with the authors, starts from independent sets of directions uniformly spread on the sphere, and then finds the optimal composition of rigid rotations of each set with respect to the others, so that the projection of the rotated sets of points onto the same unit spheres gives the most uniform coverage possible. Finally, the method in [START_REF] Ye | Efficient interlaced multi-shell sampling scheme for reconstruction of diffusion propagators[END_REF] is based on a geometrical construction, using the concept of uniform and dual polyhedra. The vertices of two polyhedra dual to each other give two complementary sets of directions, that are used to design multi-shell sampling strategy. The drawback with this technique is that it only provides two separate sets of directions, which is not enough for more than 2 shells sampling. Moreover, the choice of the number of points per shell is limited by the finite number of uniform polyhedra (see [START_REF] Harel | Uniform solution for uniform polyhedra[END_REF] and the related software Kaleido for a mathematical construction).

In this work, we present a general method to design multi-shell sampling protocols, with a uniform coverage on each shell, and a global uniform angular coverage. We propose a cost function that generalizes electrostatic repulsion to multi-shell, which is presented in the first section, together with an efficient method to minimize this cost function. We also show that the same objective function can be used to design incremental acquisition schemes, that are compatible with prematurely aborted scans. In the results section, we compute some geometrical properties of the optimal point sets. We demonstrate the properties of our new sampling scheme through Monte-Carlo simulations.

A preliminary version of this work was presented in [START_REF] Caruyer | Optimal design of multiple q-shells experiments for diffusion mri[END_REF][START_REF] Caruyer | Incremental gradient table for multiple q-shells diffusion mri[END_REF].

Methods

The cost function we present in this section is a generalization of the electrostatic repulsion in high angular resolution diffusion imaging [START_REF] Jones | Optimal strategies for measuring diffusion in anisotropic systems by magnetic resonance imaging[END_REF][START_REF] Papadakis | Minimal gradient encoding for robust estimation of diffusion anisotropy[END_REF][START_REF] Jansons | Persistent angular structure: new insights fom diffusion magnetic resonance imaging data[END_REF], adapted to the problem of multi-shell sampling.

The cost function is made of two terms: the first term controls the repartition on each shell, while the second controls the angular uniformity of the set of directions as a whole. A point set minimizing this cost function can be constructed by direct minimization, and we report how to minimize this cost function. Alternatively, to construct incremental acquisition schemes, that are compatible with prematurely stopped scans (see [START_REF] Dubois | Optimized diffusion gradient orientation schemes for corrupted clinical dti data sets[END_REF] for instance in q-ball imaging), the same cost function can be minimized incrementally, following the method in [START_REF] Deriche | Optimal real-time q-ball imaging using regularized Kalman filtering with incremental orientation sets[END_REF].

Multi-shell q-space sampling with uniform coverage

In the design of q-ball imaging acquisition, the most popular approach is to spread out the sampling directions as uniformly as possible on the sphere [START_REF] Jones | Optimal strategies for measuring diffusion in anisotropic systems by magnetic resonance imaging[END_REF][START_REF] Papadakis | Minimal gradient encoding for robust estimation of diffusion anisotropy[END_REF][START_REF] Jansons | Persistent angular structure: new insights fom diffusion magnetic resonance imaging data[END_REF]. Similarly, we propose to construct the sampling directions on each shell, with the following principles: First, for the sake of rotational invariance, it is desirable that the coverage of each sphere be as uniform as possible. Second, we require the whole set of directions (i.e. the set of all sampling points in q-space, reprojected onto the unit sphere) to be as uniform as possible. We therefore propose a split cost function

V = αV 1 + (1 -α)V 2
, where V 1 and V 2 address, respectively, the two principles stated above.

The cost function, is a natural generalization of the method in [START_REF] Jones | Optimal strategies for measuring diffusion in anisotropic systems by magnetic resonance imaging[END_REF][START_REF] Jansons | Persistent angular structure: new insights fom diffusion magnetic resonance imaging data[END_REF] to multi-shell,

V 1 = 1 S S s=1 1 K 2 s i<j v(u s,i , u s,j ), ( 2
)
where S is the number of shells, K s is the number of points on shell s, and {u s,k , k = 1 . . . K s } are the unit directions on shell s. The magnitude of the elementary electrostatic repulsion force, v, between two points u and v on the unit sphere is "according to Coulomb's law, inversely proportional to the square of the distance between the charges" (see [START_REF] Jones | Optimal strategies for measuring diffusion in anisotropic systems by magnetic resonance imaging[END_REF]),

v(u, v) = 1 ||u -v|| 2 + 1 ||u + v|| 2 .
(

) 3
This elementary repulsion force ensures that a sampling direction, u, and its opposite -u play the same role [START_REF] Jones | Optimal strategies for measuring diffusion in anisotropic systems by magnetic resonance imaging[END_REF][START_REF] Cook | Optimal acquisition orders of diffusion-weighted mri measurements[END_REF]. Minimizing the cost function V 1 , corresponds to minimizing the electrostatic cost function on each shell, separately. Note that the cost of shell s is weighted by 1/K 2 s ; indeed, the cost function increases as K 2 [START_REF] Cook | Optimal acquisition orders of diffusion-weighted mri measurements[END_REF], and this weighting ensures equal importance of the uniformity for each shell. So that the whole set of sampling directions be uniform, we also consider this second cost function

V 2 = 1 K 2 s =t Ks i=1 Kt j=1 v(u s,i , u t,j ), ( 4
)
where K is the total number of sampling points in Q-space. This cost function, V 2 , is the sum of magnitude of all electrostatic forces between any two directions coming from two different shells.

This will ensure that the optimal directions are different from one shell to another.

Finally, we construct our set of points by minimizing the cost function

V = αV 1 + (1 -α)V 2 , (5)
under the constraint that each direction lies on the unit sphere, arg min This minimization is a non-convex optimization problem under quadratic equality constraint.

u s,k ∈R 3 V (u s,k ), subject to ∀s ∈ {1 . . . S}, ∀k ∈ {1 . . . K s }, ||u k,s || 2 = 1. (6)
We implemented the minimization through sequential least squares quadratic programming, with analytically pre-calculated gradient to speedup the optimization. The gradient of the functional V is derived from the partial derivatives of v(u, v),

∂v ∂x u = -2 x u -x v ||u -v|| 4 + x u + x v ||u + v|| 4 , (7)
where the directions are parameterized by their Cartesian coordinates, u = [x u y u z u ] T . Besides, the equality constraint is also easily implemented, as the normal to the surface defined by the equality constraint is known: it is simply u s,k . Our implementation, in Python TM , makes use of the fmin_slsqp routine from the SciPy optimization package [START_REF] Jones | SciPy: Open source scientific tools for Python[END_REF]. An example of optimal configuration on 3 shells, with 30 points per shell, and a weighting parameter α = 0.5 is plotted on As a result, the global coverage is non-uniform, and so we get almost "aligned" sampling on both shells. (right) α = 1: the global angular coverage is uniform, but if we consider the repartition of points on shell 1 (green dots) or shell 2 (red dots) separately, the angular coverage is not uniform. each single shell. These two extremal cases are illustrated on Fig. 2.

We plot on Fig. 3 the value of both energies V 1 and V 2 for the configuration minimizing V , while the α parameter varies from 0 to 1. Except near these extremal values of α, the solution is not much sensitive to the choice of α. We chose α = 0.5 throughout our experiments.

Incremental construction of an optimal arrangement

Some recent work on single shell acquisition design [START_REF] Cook | Optimal acquisition orders of diffusion-weighted mri measurements[END_REF][START_REF] Dubois | Optimized diffusion gradient orientation schemes for corrupted clinical dti data sets[END_REF][START_REF] Deriche | Optimal real-time q-ball imaging using regularized Kalman filtering with incremental orientation sets[END_REF] have shown the practical interest of incrementally optimal sampling sequences. In short, these acquisition sequences are carefully designed and ordered, so that an interruption at any point of the acquisition would result in an approximately optimal design. Such designs are compatible with acquisitions that have been stopped early because of subject motion or request. In q-ball imaging, the greedy optimization method proposed in [START_REF] Deriche | Optimal real-time q-ball imaging using regularized Kalman filtering with incremental orientation sets[END_REF] was shown to be very computationally attractive. Besides, in Q-ball imaging, point sets generated with this method well compare to results obtained by a modification of the cost function, introducing weighting in the sum of magnitude of electrostatic elementary repulsion forces, as in [START_REF] Dubois | Optimized diffusion gradient orientation schemes for corrupted clinical dti data sets[END_REF][START_REF] Cook | Optimal acquisition orders of diffusion-weighted mri measurements[END_REF]. For these reasons, we adapt the approach in [START_REF] Deriche | Optimal real-time q-ball imaging using regularized Kalman filtering with incremental orientation sets[END_REF] to the problem of multi-shell design.

To find a nearly optimal design for K points, we construct the point set incrementally. At each iteration k, we first select the sphere s k on which to add the k-th point, so that the ratio of samples on each shell is approximately respected. Then we find the direction u minimizing V , while the k -1 previously defined directions remain fixed. This minimization algorithm is very efficient, as at each step, we need to optimize a single direction (two parameters). Furthermore, as we show in the results section, the geometrical properties of this incrementally constructed point set are close to the optimal one, provided that enough sample points are considered.

Results

Geometrical properties of multi-shell point sets

We present in this section some geometrical measures on the generated point sets, and compare to optimal configurations on the sphere. We also compare the incremental point sets configurations to the optimal configurations. As a measure of the uniformity of a set of directions (right) globally. We compare the sets generated using generalized electrostatic analogy, for S = 2, 3 and 4 shells, to the optimum electrostatic point sets in one shell.

X = {u k , k ∈ K}, 0 10
we compute the minimum angular distance between any two points, defined as

d(X ) = min i =j arccos(|u i • u j |).
The larger d(X ), the more uniform is the set of directions X . For reference, we compare to the minimum angular distance of the optimal electrostatic distribution with the same number of points.

We generated sampling protocols on S = 2, 3 and 4 shells, with the same number of points per shell. The total number of points ranges from K = 12 to K = 240. We report on Fig. 4 (left)

the minimum angular distance between any two points on the same shell. Similarly, we report the minimum angular distance between any two points of the whole set of points, on Fig. 4 (right).

The minimum angular distance remains close to that of the optimal point sets with same number of points, which suggests that the method offers a very good compromise between uniformity per shell, and global uniformity. This means there is very little price to pay on the uniformity per shell, to get a globally angular uniform configuration. This validate our approach, which consists in optimizing both the angular distribution of the sampling protocol, and the distribution on each single shell.

We also evaluate the geometrical properties of incrementally constructed multi-shell point sets.

We report on Fig. 5 the minimum angular distance between any two points, per shell and as a whole. We observe that, provided that enough measurements are considered, the values of the minimum angular distance are reasonably close to that of the optimal configuration on a single shell, for the same number of points.

Finally, we compare our method to three other state-of-the-art approaches that were briefly presented in the introduction, and aim at constructing acquisition designs on multiple shells with uniform angular coverage. The first method, called sparse and optimal acquisition (SOA) designs [START_REF] Koay | Sparse and optimal acquisition design for diffusion mri and beyond[END_REF], starts from a uniform arrangement of K points on the unit sphere. Then it searches for the optimal partition of this set of points into S subsets (corresponding to S spheres), with K s points on shell s, with respect to an objective function. We use the Java implementation provided in the HI-SPEED software, publicly available on the first author's website. The second method [START_REF] Zhan | Differential information content in staggered multiple shell HARDI measured by the tensor distribution function[END_REF] also starts from a set of K points uniformly spread on the unit sphere, though obtained with a different method [START_REF] Wong | A strategy for sampling on a sphere applied to 3d selective rf pulse design[END_REF]. Then it constructs S subsets which are nearly-optimal, minimizing electrostatic energy in a greedy fashion [START_REF] Zhan | How does angular resolution affect diffusion imaging measures?[END_REF]. We will refer to this method as discrete greedy. The last method was presented in [START_REF] Santis | Improved precision in the charmed model of white matter through sampling scheme optimization and model parsimony testing[END_REF] and, in deeper details, in a private communication with the authors. This method starts from S independent sets of directions uniformly spread on the sphere, and then finds the optimal composition of rigid rotations of each set with respect to the others, so that the projection of the rotated sets of points onto the same unit spheres gives the most uniform coverage possible. We will Shell 1

Shell 2 Shell 3 global

K 1 = 28 K 2 = 28 K 3 = 28 K = 84
SOA design [START_REF] Koay | Sparse and optimal acquisition design for diffusion mri and beyond[END_REF] 15.9

K s = 28.
On the last row, we also report the minimum angular distance of the optimal single shell point set with the same number of points (28 points for the first three rows, and 84 for the last row).

refer to this last method as solid rotation optimization. The optimization is done using simulated annealing, and implemented in Python TM , with the anneal routine from the SciPy optimization package [START_REF] Jones | SciPy: Open source scientific tools for Python[END_REF]. The minimum angular distance per shell and globally are compared and reported in Table 1.

Minimization of the cost function

As already stated above, the minimization of cost function V is a non-convex optimization problem, and we have no mathematical clue on the convergence of the optimization. We decided to test this point by repeating the optimization process, with random initial guess. The difference observed in the cost function at the end of the optimization is of the order of 10 -2 , which is negligible.

Another concern is about the computational cost of the minimization. While we cannot precisely measure the implementation of the SOA method [START_REF] Koay | Sparse and optimal acquisition design for diffusion mri and beyond[END_REF], as the HI-SPEED software packages come with a GUI, we observed a time of several hours for the computation of a point set with 100 points on a modern desktop computer. Our method, the solid rotation method [START_REF] Santis | Improved precision in the charmed model of white matter through sampling scheme optimization and model parsimony testing[END_REF] and the discrete greedy method [START_REF] Zhan | Differential information content in staggered multiple shell HARDI measured by the tensor distribution function[END_REF] give a result within minutes on the same hardware, and therefore are more attractive from a computational perspective.

Monte-Carlo simulations

We provide in this section an experimental validation of our sampling approach. We test the sampling schemes, and evaluate the angular resolution of generalized electrostatic, compared to the other, above mentioned approaches.

We simulated a single fiber by a cylinder with impermeable wall. The diffusion signal was simulated following [START_REF] Soderman | Restricted Diffusion in Cylindrical Geometry[END_REF]. The cylinder dimensions were set to L = 5 mm, and ρ = 5 µm; and the diffusion time is set to τ = 20.8 ms. The signal is simulated for different acquisitions schemes, with a maximum wavevector norm q max = 60 mm -1 , which corresponds to a b value b max = 3000 s • mm -2 .

The diffusion weighted signal was corrupted by Rician noise, with corresponding SNR = 25 with respect to the baseline, non diffusion-weighted signal. This corresponds to a SNR range from 0.1 to 7.0 on the diffusion-weighted signal.

From these noisy measurements, we fit a diffusion tensor model, using the weighted linear least squares method [START_REF] Basser | Estimation of the effective self-diffusion tensor from the nmr spin echo[END_REF]. This process is repeated 100 times for a given cylinder model, to get an estimate of the average angular error. To test the rotational invariance of acquisition schemes, this in turn is repeated 2000 times, for cylinders with different axis orientations. Acquisition protocols are defined on S = 3 shells in q-space, where the shell radii q 1 , q 2 and q 3 are evenly distributed between q min and q max . The number of points per shell is the same on each shell K s = 8, which corresponds to a total of K = 24 measurements, typical in diffusion tensor imaging. For reference, we also compare to a single shell experiment with K = 24 measurements, q = q max . Results are reported on Fig. 6. and in Table 2.

We also compare the accuracy of multiple fiber detection and discrimination. The signal is now simulated using a mixture of 2 Gaussian distributions, with equal compartment size, E(q) = 0.5(exp(-q T D 1 q) + exp(-q T D 2 q))).

The corresponding diffusion tensors D Table 3: Crossing angle for a reconstructed fiber crossing (ground truth is 60 • ). See legend in Fig. 7 for more details. SNR = 25 with respect to the baseline, non diffusion-weighted signal. This corresponds to a SNR range from 1.0 to 20.0 on the diffusion-weighted signal.

The signal is estimated in the continuous, modified spherical polar Fourier (mSPF) basis with

Laplace regularization [START_REF] Caruyer | Diffusion MRI signal reconstruction with continuity constraint and optimal regularization[END_REF]. This is an orthonormal basis of continuous functions to represent the diffusion signal. The coefficients of the signal are estimated analytically, with a Laplace regularization penalty to promote smoothness. From the coefficients of the signal in this basis, the orientation distribution function (ODF) in constant solid angle is estimated analytically using the SPFI method in [START_REF] Cheng | Model-free, regularized, fast, and robust analytical orientation distribution function estimation[END_REF]. The angular and radial truncation orders were set to L = 6 and N = 2, respectively. The maxima of the ODF are computed with a discrete search on a set of 10000 points on the sphere, leading to an accuracy of 0.1 • .

We compare several acquisition protocols with K = 100 points, and a maximum b-value of 2500 s • mm -2 . For radial sampling, discrete greedy sampling [START_REF] Zhan | Differential information content in staggered multiple shell HARDI measured by the tensor distribution function[END_REF] and solid rotation sampling [START_REF] Santis | Improved precision in the charmed model of white matter through sampling scheme optimization and model parsimony testing[END_REF],

we use 4 shells with radii linearly spaced between q min and q max , and K s = 25 points per shell. For the SOA design, we use 10 shells, each with K s = 10 points per shell, as advocated in [START_REF] Koay | Sparse and optimal acquisition design for diffusion mri and beyond[END_REF]. Finally, for the generalized electrostatic, we use the sampling strategy proposed in [START_REF] Caruyer | A computational framework for experimental design in diffusion MRI[END_REF] and optimized for reconstruction in mSPF with radial truncation order N = 2. This corresponds to a design on S = 2

Crossing angle

Figure 7: Crossing angle for a reconstructed fiber crossing. The q-space signal was reconstruction in the mSPF basis, with Laplace regularization (see [START_REF] Caruyer | Diffusion MRI signal reconstruction with continuity constraint and optimal regularization[END_REF]), fixed regularization weight λ = 0.5. The ODF was further estimated using the SPFI method in [START_REF] Cheng | Model-free, regularized, fast, and robust analytical orientation distribution function estimation[END_REF]. The boxplot represents min and max values, as well as the first and third quartiles. We also compare to the ODF reconstructed in constant solid angle [START_REF] Aganj | Reconstruction of the ODF in single and multiple shell q-ball imaging within constant solid angle[END_REF], for a single shell experiment with the same number of points K = 100.

See also Table 3.

struction with our method gives the lowest dispersion in angle estimation, and a very small bias, when it comes to multiple shell detection.

Discussion

We have proposed a novel, efficient method to generate sampling schemes for multi-shell acquisitions. Our method is dedicated to the placement of acquisition directions on each sphere, with a good angular distribution per shell and as a whole. The cost function we propose is an extension of the electrostatic repulsion to multi-shell. With an adequate minimization algorithm, we show that this cost function can be used to create acquisition with incremental angular distribution, that is compatible with prematurely stopped scans.

The focus of the study is more on the angular distribution of points. As for the choice of number of shells, shell radii, number of points per shell, we suggest to follow the strategy optimized for reconstruction in mSPF, recently presented in [START_REF] Caruyer | A computational framework for experimental design in diffusion MRI[END_REF].

The geometrical properties of the generated point sets are satisfactory: the angular distribution per shell is indeed very close to the optimal arrangements in q-ball imaging. This is also the case for the angular distribution as a whole. The compromise between the global angular distribution and the angular distribution on each shell is better than for the other methods to which we compared:

we are able to increase the global uniformity, at a very small cost on the uniformity on each shell.

In simulations, for the reconstruction of single and multiple fiber compartments, this sampling strategies gives a better angular precision than prior arts methods. Single fiber orientation is determined with a lower angular error than simply radial sampling, and the crossing angle in two fibers case is also more accurate.

Figure 1 :

 1 Figure 1: Optimal point set on 3 shells, 30 points per shell. The weighting parameter is set to α = 0.5. (left) sample points in the q-space, (right) sampling directions, i.e. the sampling points reprojected onto the unit sphere.

Fig. 1 .

 1 Fig. 1.The weighting parameter α ∈ [0, 1] balances the trade-off between uniformity on each shell (measured by cost function V 1 ), and angular uniformity of the sampling scheme as a whole (measured by function V 2 ). The case α = 0 is somehow underdetermined, as any rotation of one shell with respect to the others does not change the cost function. At the opposite, the case α = 1 would put all the weight on the global angular uniformity, but nothing would prevent poor coverage of

Figure 2 :

 2 Figure 2: Extremal values of the weighting factor α: minimum cost configurations for S = 2 shells, K 1 = K 2 = 50 points per shell, points reprojected on the unit sphere. (left) α = 0: all the importance is given to the uniformity of the configuration on each shell, independently of the other.

2 Figure 3 :

 23 Figure 3: (left) Cost function V 1 , measuring uniformity on each shell, and (right) V 2 , measuring uniformity of the set of directions as a whole. Cost function evaluated for the optimum configurations, for various α. Except near 0 and 1, the solution does not depend strongly on the choice of α.

Figure 4 :

 4 Figure 4: Minimum angular distance ( • ) between any two points (left) within the same shell and

Figure 5 :

 5 Figure 5: Geometrical properties of the incrementally constructed point sets on S shells. Minimum angular distance ( • ) between any two points (left) within the same shell and (right) globally. We compare the incremental sets generated using generalized electrostatic analogy, for S = 2, 3 and 4 shells, to the optimum electrostatic point sets in one shell.

1 and D 2 Figure 6 :

 26 Figure 6: Average angular error between the principal direction of the estimated tensor and the axis of the cylinder. Signal simulation and reconstruction were repeated with 2000 different cylin-der orientations following a random uniform distribution on the sphere, and 100 repetitions with random, Rician noise for each orientation. We plot here the median, first and third quartiles, and extremal values of the angular error (in degrees). The height of each boxplot represents how rotationally invariant each method performs.. See Table2for quantitative results.

Table 1 :

 1 Minimum angular distance for a 3-shell design, with equal number of points per shell

	•	15.8 •	16.7 •	15.1 •

Table 2

 2 ±0.22 2.71 • ±0.09 2.44 • ±0.10 2.44 • ±0.10 2.47 • ±0.11 2.88 • ±0.34

	for quantitative results.

Table 2 :

 2 Angular error between the principal direction of the estimated tensor and the axis of the cylinder. We report the median and angular deviation. See Fig.6for more details.

	radial	CSA	generalized	soa	solid	discrete
		ODF	electrostatic		rotation	greedy
	estimated crossing 58.5 • ±3.8 66.1 • ±2.6	61.5 • ±2.7	59.4 • ±4.0 59.2 • ±3.5 59.2 • ±3.7

  shells, with b-values 900 and 2500 s • mm -2 , and K 1 = 24, k 2 = 76 points per shell, respectively. Results of the reconstructed crossing angle are reported on Fig 7 and in Table3. As for the previous experiments, the simulation was repeated 1000 times, with Rician noise and random rotation of the fibers to test the rotational invariance of the reconstruction.

	40 45 50 55 60 65 70 75 80	radial	CSA ODF	generalized electrostatic	soa	solid rotation	greedy discrete

While for the reconstruction of a single fiber, either method give comparable results, the recon-

Acknowledgements

Sampling schemes can be downloaded from the author's webpage. This work was partly supported by Inria CD-MRI Associate Team program, ANR NucleiPark, France-Parkinson Foundation, and NIH Grants P41 RR008079 (NCRR), P41 EB015894 (NIBIB), P30 NS057091, and P30 NS076408. C. Lenglet and G. Sapiro were partly supported by the NIH Human Connectome Project U54 MH091657 and Grant R01 EB008432. G. Sapiro received additional support from ONR, NGA, NSF, DARPA, NSSEFF, and ARO. The authors thank the authors of SOA design method [25] for their publicly available implementation, in the HI-SPEED software (https://sites.google.com/site/hispeedpackets).