Marie-Liesse Cauwet
email: cauwet@lri.fr

Olivier Teytaud

Shih-Yuan Chiu

Kuo-Min Lin

Shi-Jim Yen

David L St-Pierre

Fabien Teytaud
email: fabien.teytaud@lisic.univ-littoral.fr

Parallel Evolutionary Algorithms Performing Pairwise Comparisons

Keywords: G.1.6 [Optimization]: Unconstrained optimization Differential Evolution, Parallelism, Particle Swarm Optimization

We study mathematically and experimentally the convergence rate of differential evolution and particle swarm optimization for simple unimodal functions. Due to parallelization concerns, the focus is on lower bounds on the runtime, i.e upper bounds on the speed-up, as a function of the population size. Two cases are particularly relevant: A population size of the same order of magnitude as the dimension and larger population sizes. We use the branching factor as a tool for proving bounds and get, as upper bounds, a linear speed-up for a population size similar to the dimension, and a logarithmic speed-up for larger population sizes. We then propose parametrizations for differential evolution and particle swarm optimization that reach these bounds.

INTRODUCTION

Evolutionary algorithms are population-based algorithms designed to solve black-box optimization problems. An evolutionary algorithm typically (i) generates a population of search points (ii) selects a candidate solution for the optimization problem. The selection is performed thanks to the ranking/comparisons between the search points. The population makes parallelization highly relevant when the population is large: each point can be evaluated on a single processor [START_REF] Storn | Differential evolution: A simple and efficient heuristic for global optimization over continuous spaces[END_REF][START_REF] Kennedy | Particle swarm optimization[END_REF].

We first introduce parallel black-box optimization in Section 1.1. We then present the convergence rate of an algorithm in Section 1.2 before explaining what is the branching factor in Section 1.3. We then focus on parallel black-box optimization algorithms which use only "paired" comparisons, as explained in Section 1.4. Section 2 applies branching factor analysis to these algorithms, and provides lower bounds on runtimes. Sections 3 and 4 confront bounds with experimental rates for differential evolution and particle swarm optimization respectively.

Throughout the paper, unless stated otherwise, log refers to the natural logarithm.

Parallel black-box optimization

Typical optimization problems consist in searching for the minimum x * of some objective function f : R D → R, i.e. x * is such that ∀x, f (x * ) ≤ f (x). f is also known as the fitness function and D represents the dimension of the problem. A black-box optimization consists in doing so by successive calls to f only. f is seen as an oracle; no internal property of f is used, and the goal of the optimization algorithm consists in finding a good approximation of x * , within a moderate number of calls to the objective function.

Parallel synchronous black-box optimization is the setting in which p simultaneous calls to the objective function are possible. Therefore, p is the number of processors or cores, possibly on GPU architectures, which allow large numbers of cores. When measuring speed-up, i.e. the factor by which speed is improved when using p processors, there are several cases: (i) Cases in which the objective function is fast. Then, the internal cost of the optimization algorithm and possibly communication costs are not negligible. (ii) Cases in which the cost of the objective function is not approximately constant. Then, asynchronous effects matter. (iii) Cases in which the computational cost is mainly in the objective function, and the computational cost of the fitness function is constant. This leads to a synchronous case.

We focus on case (iii). Additionally, we consider that the structure of the optimization algorithm is not modified: we just use an increased population size, so that we can simultaneously evaluate one fitness value on each core or each processor. We will study the speed-up as a function of the population size, assuming that the population size is the number of cores. This is a special case, but a very usual one in practice.

Convergence Rate

We aim to study the effectiveness and limits of some evolutionary algorithms. As such, we define the convergence rate of an algorithm. At each iteration n, the evolutionary algorithm provides an approximation xn of the optimum x * . We consider the following convergence measure:

lim n→∞ 1 n log xn -x * → -C λ,D (1)
where C λ,D is a positive constant, depending on the population size λ and the dimension D only. Both theoretical and empirical results ( [START_REF] Auger | Convergence results for (1,λ)-SA-ES using the theory of ϕ-irreducible Markov chains[END_REF][START_REF] Auger | sigma,eta) : quasi-random mutations for evolution strategies[END_REF][START_REF] Beyer | The Theory of Evolution Strategies[END_REF][START_REF] Rechenberg | Evolutionstrategie: Optimierung Technischer Systeme nach Prinzipien des Biologischen Evolution[END_REF]) show that a wide range of evolutionary algorithms may converge in the sense of Eq. 1 when the objective function is unimodal and reasonably smooth. Our purpose is to exhibit some bounds on this convergence rate C λ,D when it is well defined. We also get lower bounds on the corresponding lim inf.

From the definition of the convergence rate (Eq. 1), it is straightforward that the speed-up of an evolutionary algorithm is equal to the convergence rate up to a multiplicative constant. Hence, in the following, results on the convergence rate also hold for the speed-up.

We have some extra informations available on the convergence rate C λ,D , as follows [START_REF] Fournier | Lower Bounds for Comparison Based Evolution Strategies using VC-dimension and Sign Patterns[END_REF].

(i) Convergence with a population size equal to the dimension, at best:

Linear scalability: lim λ→∞ C λ,λ = C linear . (2)
Eq. 2 is written for C λ,λ . Indeed, Eq. 2 may also hold for C λ,D(λ) where λ → D(λ) is an approximately linear function, i.e. D(λ)/λ = Θ [START_REF] Arabas | Parade: A massively parallel differential evolution template for easea[END_REF]. [START_REF] Fournier | Lower Bounds for Comparison Based Evolution Strategies using VC-dimension and Sign Patterns[END_REF]Section 5] has shown that general classes of comparison-based algorithms can not do better than such a convergence.

This equation means that we "cancel" the curse of dimensionality thanks to

• A parallel implementation with one individual evaluated on each processor.

• A linear number of processors (linear in the dimension).

(ii) Convergence with a large population size:

Logarithmic speed-up:

lim λ→∞ C λ,D = C large,D log(λ) (3) with lim D→∞ D × C large,D = C large (4)
For a wide class of comparison-based algorithms, including evolution strategies and pattern search methods, no algorithm can do better than this bound ([13, Fig. 2 and Prop. 8]).

We here show that these results also hold in the case of differential evolution and particle swarm optimization. Explicit bounds on C linear and C large are available thanks to the branching factor analysis sketched in Section 2.

The convergence rate above (Eq. 1) is a convergence rate in the search space; we can also consider the convergence rate in the fitness space:

lim n→∞ 1 n log f (xn) -inf x f (x) = -C ′ λ,D . (5)
On the sphere function, C ′ λ,D = 2C λ,D . We will use this convergence rate C ′ λ,D (in fitness space) in all our experiments.

Branching factor

Many evolutionary algorithms, after having evaluated the population, keep only a concise representation of the evaluations they have performed. This concise representation is typically:

• The ordered list of the indices of the µ best among the λ generated individuals. There are λ µ × (µ!) possibilities.

• The family (unordered) of the indices of the µ best among the λ generated individuals. There are λ µ possibilities.

A crucial property of such algorithms is that the number of possibilities is finite. The branching factor is this number. For evolution strategies, the branching factor can be strongly reduced, by assuming some regularity of the objective function (VC-dimension assumption on the objective function [START_REF] Vapnik | The Nature of Statistical Learning[END_REF][START_REF] Devroye | A probabilistic Theory of Pattern Recognition[END_REF]).

The branching factor is then used to prove that comparisonbased algorithms have usually convergence rate increasing at most logarithmically with the population size λ as in Eq. 3. This upper bound is proved for evolution strategies in [13, Fig. 2]. Moreover, this bound is tight: ad hoc variants of evolutionary algorithms match it. Thanks to some specific tuning, [START_REF] Teytaud | Log(lambda) Modifications for Optimal Parallelism[END_REF] shows empirically that this logarithmic speedup can be reached by some evolutionary algorithms, namely CMSA, EMNA and CMAES. Similarly, the branching factor leads to finding a lower bound on the convergence rate provided that the population size is limited. Here again, this bound is tight, hence Eq. 2 holds for comparison-based algorithms (see [START_REF] Fournier | Lower Bounds for Comparison Based Evolution Strategies using VC-dimension and Sign Patterns[END_REF]Section 5]).

Differential evolution and particle swarm optimization

Two families of evolutionary algorithms are studied in the present paper, namely Particle Swarm Optimization (PSO) and Differential Evolution (DE). They share several common points. They both are population based optimization algorithms and do not require heavy calculations for gradient or Hessian updates. Furthermore, the run of these algorithms depends on "paired" comparisons only: instead of comparing or sorting most of the population, these two algorithms compare each individual to only one possible mutation. This leads to specific mathematical properties, discussed in Section 2.

Among evolutionary algorithms, Differential Evolution was, from the very beginning [START_REF] Storn | Differential evolution: A simple and efficient heuristic for global optimization over continuous spaces[END_REF], presented as a tool with high relevance for parallel optimization in continuous domains.

Particle Swarm Optimization [START_REF] Kennedy | Particle swarm optimization[END_REF][START_REF] Shi | A Modified Particle Swarm Optimizer[END_REF] is also a general purpose black-box optimization algorithm, aimed at finding quickly approximate solutions to optimization problems.

In the implementation of these algorithms (see Algorithms 1 and 2), we use two loops per generation: one loop for computing fitness values of candidate vectors; one loop for updating "best" points (best search point so far, and best point so far for each given particle for PSO). One can implement variants of PSO and DE in which this is done in one loop only. However, such variants are less parallel, because the mutations, which depend on "best" points, might change during the loop and therefore cannot be computed simultaneously.

Choosing the right parameters for these algorithms is a challenging task. In this paper, we define rules for choosing these parameters in a parallel black-box optimization setting. We restrict this experimental part of the study to the case of unimodal well conditioned fitness functions. We validate our empirically developed formula for the parameters in front of both (i) standard values used for DE and PSO and (ii) theoretically known bounds.

Differential evolution. Differential evolution (DE) [START_REF] Storn | Differential evolution: A simple and efficient heuristic for global optimization over continuous spaces[END_REF] is one of the main evolutionary algorithms for continuous optimization. It performs well on optimization testbeds [START_REF] Storn | Differential evolution: A simple and efficient heuristic for global optimization over continuous spaces[END_REF][START_REF] Ardia | Differential Evolution with DEoptim: An application to non-convex portfolio optimization[END_REF][START_REF] Ardia | Jump-diffusion calibration using Differential Evolution[END_REF][START_REF] Pošík | Jade, an adaptive differential evolution algorithm, benchmarked on the bbob noiseless testbed[END_REF]. At each iteration, each point pi of the population undergo mutations of some of its components, then comparisons are performed between the old point and its corresponding mutated version. The best of the two points is selected for the next iteration. Different parameters are involved in the mutation step. A classical version of DE is presented in Algorithm 1. Many variants exist, including selfadaptive parameters [START_REF] Brest | Self-adapting control parameters in differential evolution: A comparative study on numerical benchmark problems[END_REF][START_REF] Liu | A fuzzy adaptive differential evolution algorithm[END_REF][START_REF] Yang | Self-adapting differential evolution algorithm with chaos random for global numerical optimization[END_REF][START_REF] Price | Differential Evolution -A Practical Approach to Global Optimization[END_REF] and meta-heuristics for choosing parameters [START_REF] Pedersen | Tuning & simplifying heuristical optimization[END_REF]. In some variants [START_REF] Storn | Differential evolution: A simple and efficient heuristic for global optimization over continuous spaces[END_REF] On the technical side, the large scale parallelization of the evaluations is performed with CUDA in [START_REF] Qin | An improved cuda-based implementation of differential evolution on gpu[END_REF]. [START_REF] Olensek | A new asynchronous parallel global optimization method based on simulated annealing and differential evolution[END_REF] focuses on the asynchronous parallelization, which is not considered in

′′ i,j ← p ′ i,j if rand < Cr or if j = R p ′′ i,j ← pi,j otherwise end for Compute f (p ′′ i ) end for for i ∈ {1, 2, . . . , λ} do if f (p ′′ i ) < f (pi) then pi ← p ′′ i if best = 1 & f (pi) < f (g) then g ← pi
U pdate xi ← xi + vi Compute f (xi) end for for each i ∈ P do if f (xi) < f (pi) then pi ← xi if f (pi) < f (g) then g ← pi end if end if end for n ← n + 1 end while Return g
It has several parameters leading to complicated dynamics; for example w (see Algorithm 2) is often thought of as an inertia, but it can be negative or greater than 1. Therefore, many works are devoted to choosing optimally these parameters [START_REF] Parsopoulos | Parameter selection and adaptation in unified particle swarm optimization[END_REF][START_REF] Trelea | The particle swarm optimization algorithm: convergence analysis and parameter selection[END_REF][START_REF] Clerc | The particle swarmexplosion, stability, and convergence in a multidimensional complex space[END_REF].

As explained before, PSO is relevant for parallelization. [START_REF] Gardner | A speculative approach to parallelization in particle swarm optimization[END_REF] uses speculative parallelization. Speculative parallelization is known to ensure an asymptotic logarithmic convergence rate (i.e. the optimal rate as described in Eq. 3), but it does not, in the general case, ensure optimality in large dimension (Eq. 2). They use the topology for defining a parallel variant of PSO, without necessarily increasing the population size. We will here consider the simple increasing of the population size. [START_REF] Mcnabb | Parallel pso using mapreduce[END_REF] focuses on the implementation part of a parallel PSO, in particular taking care of node failure and communication costs, thanks to MapReducewe will here assume that communication costs can be neglected, e.g. due to expensive fitness functions. [START_REF] Mahdad | Fuzzy Controlled Parallel PSO to Solving Large Practical Economic Dispatch[END_REF] uses a fuzzy controller for a parallel PSO, applied to a power dis-patch problem; the fuzzy part generates rules for the simulations. Other experimental works on parallel-PSO include [START_REF] Schutte | Parallel global optimization with the particle swarm algorithm[END_REF] (focusing on parallel global optimization, whereas we focus on unimodal functions) and [START_REF] Chang | A parallel particle swarm optimization algorithm with communication strategies[END_REF] (detailing communication strategies, which are beyond the scope of the present paper, because we focus on the case in which most of the cost is in the fitness evaluations).

Outline of this paper. The branching factor methodology has never been applied to differential evolution and particle swarm optimization. This is the purpose of the theoretical part of this paper. In this paper:

• We show mathematically that, whatever may be the parametrization, differential evolution and particle swarm optimization can not do better than a logarithmic speed-up (Eq. 3) and a linear scalability (Eq.

2).

• We estimate empirically the convergence rate, in the case of the sphere function. These convergence rates, in the easy considered setting, match the lower bounds.

• We propose some variants designed to reach a good convergence rate when the population size is large.

MATHEMATICAL ANALYSIS

In [START_REF] Fournier | Lower Bounds for Comparison Based Evolution Strategies using VC-dimension and Sign Patterns[END_REF], the branching factor was identified as a crucial component of a comparison-based algorithm. We present a general framework of optimization in Alg. 3.

Algorithm 3 General framework of algorithms with bounded branching factor K.

Initialize some state S while not finished do Generate (possibly randomly) a family E of individuals, using S only.

Let I be the extracted information, which depends on (i) the fitness function (ii) S (iii) E; we assume that I has values in a finite set {1, . . . , K}.

Update: S ← update(S, I).

Recommend: let x = recom(S) be the approximation of the optimum proposed by the algorithm. end while Define, for the theorem below, xm the approximation of the optimum that it recommends after m iterations in Algorithm 3. It is known, by the following theorem [START_REF] Teytaud | General lower bounds for evolutionary algorithms[END_REF] that the runtime necessary for hitting the optimum with probability 1 -δ and precision ǫ is lower bounded by a function of the branching factor.

Theorem 1 (Branching factor theorem).

Assume that an optimization algorithm has branching factor K. Let D be the search domain. We define F, n ǫ,δ and N (ǫ) as follow:

• We denote by F the set of objective functions on D, i.e. functions from D to R. Each function in F has a unique minimum x * (f ) in D. We assume that for each x ∈ D, there is at least one objective function f ∈ F with optimum x * (f ) = x.

• Let n ǫ,δ be the minimum number of iterations in the While loop in Algorithm 3 such that, for any f ∈ F , with probability 1 -δ, xn ǫ,δ -x * (f ) ≤ ǫ.

• We define N (ǫ) the maximum integer n such that there exist n points (x1, . . . , xn) ∈ D n satisfying xi -xj ≥ 2ǫ for any i = j.

Then, n ǫ,δ ≥ log(1 -δ) log K + log N (ǫ) log K . (6)
N (ǫ) is the so-called packing number. In [START_REF] Fournier | Lower Bounds for Comparison Based Evolution Strategies using VC-dimension and Sign Patterns[END_REF], bounds on the branching factor are computed for evolution strategies based either on a full ranking of the population or on a selection of the µ best individuals. The bounds are improved when the objective function is simple, e.g. has bounded VCdimension. The purpose of this section is to show a bound on the branching factor in the case of differential evolution (Algorithm 1) and particle swarm optimization (Algorithm 2). Lemma 1 exhibits such a bound.

Lemma 1. If the objective function is f (x) = x -x * 2 over D = (0, 1) D for some x * ∈ D, or any composition x → m( x -x * 2
) with m increasing, then the branching factor of DE/rand (see Algorithm 1, best = 0) is at most:

K ≤ min{(λ + 1) D , 2 λ }. ( 7
)
With the same assumptions, the branching factor of DE/best (see Algorithm 1, best = 1) and PSO (see Algorithm 2) is at most:

K ≤ min{λ(λ + 1) D , λ2 λ }. (8)
Proof. First, let us give the proof when DE/rand is used. The right hand side of the min corresponds to the branching factor of λ comparisons; there are λ comparisons, hence 2 λ possible comparison results, hence the 2 λ bound. We just have to show the left hand side of the min. There are λ comparisons at each iteration, each of them between two points. This does not follow the scheme in [START_REF] Fournier | Lower Bounds for Comparison Based Evolution Strategies using VC-dimension and Sign Patterns[END_REF] because, there, the branching was directly a ranking of the µ best among λ or a selection (without ranking) of the µ best among λ. Here, one of two points is selected, λ times. When p1 and p ′′ 1 are compared, p1 is selected if p1 -x * < p ′′ 1 -x * . Similarly, the result of the comparison of pi and p ′′ i depends on the position of x * w.r.t the median hyperplane of pi and p ′′ i . Therefore, the final comparison result depends on in which cell, among the cells obtained by the arrangement of λ hyperplanes, contains x * . The branching factor is therefore upper-bounded by the number of cells obtained by λ hyperplanes in R D . By Zaslavsky's Theorem (see e. We now extend to DE/best and PSO, as follows:

• When DE/best is used, the best point is one of the λ selected points; hence the additional λ factor: beyond the 2 λ possible results of the pairwise comparisons, we have λ possibilities for the best of the λ selected points. Hence the expected result for DE/best.

• When PSO is used, the "global" best so far g is one of the λ "best so far per particle" pi. Hence, there are λ possibilities; and the expected result for PSO with best particle.

Please note that if φg = 0, then the branching factor of PSO verifies the same bound as DE without best (e.g. DE/rand).

Lemma 2. For an arbitrary family F of objective functions, the branching factor of DE/rand is at most

K ≤ 2 λ , (9)
and the branching factor of PSO and of DE/best is at most

K ≤ λ2 λ .
Proof. The 2 λ and λ2 λ parts of Lemma 1 do not use any property of F .

Property 1 (Runtime of DE/rand). If the objec- tive function is f (x) = x -x * 2 over D = (0, 1) D for some x * ∈ D, or any composition x → m( x -x * 2
) with m increasing, for DE/rand, the runtime n ǫ,δ is lower bounded:

n ǫ,δ ≥ log(1 -δ) D log (λ + 1) + log (⌈1/(2ǫ)⌉) log (λ + 1) . ( 10
)
Again when applying DE/rand, but with an arbitrary family F of objective functions, if λ = D, the runtime n ǫ,δ is lower bounded:

n ǫ,δ ≥ log(1 -δ) λ log(2) + log (⌈1/(2ǫ)⌉) log 2 . ( 11
)
Proof. Just plug the bound (Eq. 7, 8 and 9) on the branching factor above in Theorem 1, Eq. 6, and use log N (ǫ) ≥ D log (⌈1/(2ǫ)⌉) when ǫ → 0 if the domain D is bounded with non-empty interior.

Property 2 (Runtime of PSO and DE/best). If the objective function is

f (x) = x -x * 2 over D = (0, 1) D for some x * ∈ D, or any composition x → m( x -x * 2 )
with m increasing, for PSO and DE /best, the runtime n ǫ,δ is lower bounded:

n ǫ,δ ≥ log(1 -δ) log(λ) + D log (λ + 1) + D log (⌈1/(2ǫ)⌉) log(λ) + D log (λ + 1) . ( 12
)
Again when applying PSO or DE/best, but with an arbitrary family F of objective functions, if λ = D, the runtime n ǫ,δ is lower bounded:

n ǫ,δ ≥ log(1 -δ) log(λ) + λ log(2) + log (⌈1/(2ǫ)⌉) 1 + log 2 . ( 13
)
Proof. See proof of Property 1.

Remarks: Eq. 10, 11, 12, 13 imply that if DE and PSO algorithms converge in the sense of Eq. 1 then

lim λ→+∞ C λ,D = O(log(λ)). ( 14
) lim λ→+∞ C λ,λ = O(1). (15)
We have therefore proved that the differential evolution algorithm can not do better than O(log(λ)) speed-up, when using λ fitness evaluations per iteration. Eq. 11 and 13 do not prevent a linear speed-up until λ = D. In the case of evolution strategies, it is known that such a linear speed-up is possible; for DE and PSO, we are not aware of such a result. We will investigate this empirically in Sections 3 and 4.

EXPERIMENTS: SPEED-UP AND SCALABILITY OF PARALLEL DIF-FERENTIAL EVOLUTION

We have shown mathematical lower bounds on the runtime when using differential evolution; hence, if Eq. 1 holds, this implies lower bounds on CR = -C ′ λ,D (see Eqs. 5, 14 and 15). We now check if these bounds are reached or approximately reached; we will see that this is approximately the case, as far as we can tell on (necessarily finite) experiments.

In all this section, experiments are performed on the sphere function x → x 2 . We use 1000 iterations for evaluating the convergence rate. So now on, the convergence rate CR is estimated as (see Eq. 5):

CR = 1 1000
(log(best f itness at 1000 iterations)

-log(initial best f itness)) .

Differential evolution depends on two parameters that take place in the mutation step of the algorithm: the differential weight F and the crossover probability Cr. We are interested in finding the parameters F and Cr which enable to reach these lower bounds. We first do experiments with the baseline differential evolution DE/rand (Algorithm 1 with best = 0) in Section 3.1. We then perform additional experiments with large population sizes and extreme values of the parameters (Section 3.2). Finally, we test adaptive variants of differential evolution (Section 3.3). The purpose of this section is to experimentally check Eq. 3 by considering a fixed dimension D = 5. We increase the population size λ, and check if the convergence rate reaches the optimal C ′ λ,D = Θ(log(λ)) rate. Section 2 has shown that it can not be better (see Eq. 14) than this log(λ) rate. Results presented in Figure 1 are averaged over 15 runs and the number of generations N is fixed to 1000. The x-axis is the base-10 logarithm of the population size, and the yaxis shows the convergence rate CR, multiplied by D, as suggested in Eq. 4. Hence, if Eq. 3 holds, we should observe curves decreasing linearly with log(λ).

Experiments with standard differential evolution

We first experiment differential evolution with various constant parametrizations (i.e. parametrizations which are independent from the population size λ). Different parameters F and Cr are tested in Figure 1a. It appears that (i) the parameter Cr seems to have no impact on the convergence rate CR; (ii) when the population size increases, a smaller parameter F is better. For example, when log 10(λ) ∈ [1, 1.3], F = 0.41 is a good parameter choice. Then when log 10(λ) ∈ [1.3, 1.5], F = 0.24 becomes a better choice. Finally, when log 10(λ) ∈ [2.7, 4.4], F = 0.07 appears to be a better parametrization. That is why from now on we experiment with some differential evolution parametrizations with F decreasing when the population size λ increases. In Figure 1b, we investigate the case of parameter F equal to .5 log(λ) -A for various parameters A. Parameter A fixed to 0.76 is seemingly a good choice, as the corresponding curve is linear for log 10(λ) large enough; whereas other choices of parametrizations lead to plateauing curves. Figure 1c displays results for various values of F when F = .5λ -A . For several values of A (A = 0.4, A = 0.11), the curves have the expected behavior: quasi-linear, when log 10(λ) is large enough. The slopes seem also better than in the case of F = .5 log(λ) -0.76 . Specifically, F = .5λ -0.4 reaches the best slope of these experiments. Hence, a detailed look at results shows that, for population sizes as in these experiments,

Cr = .44, F = .5λ -0.4 (16)
is a reasonable tuning -though only mathematics could validate the (expected) logarithmic asymptotic convergence rate (≃ log(λ)) for λ large. The purpose of this section is to experimentally get Eq. 2.

We set D = λ/2, and we observe what happens when λ goes to infinity. Results are presented in Figure 2. Each experiment is launched 15 times and the number of generations N is set to 1000. The x-axis is the base-10 logarithm of the population size, and the y-axis shows the convergence rate CR. We expect a negative constant convergence rate if Eq. 2 holds in the case of differential evolution (see Eq. 15).

First, various constant parameters F and Cr are tested (see Figure 2a); in this case, none of the results are satisfactory as the curves converge toward 0: the algorithm does not converge. Figure 2b exhibit experiments when F is inversely proportional to an exponent of log(λ). All but one converge toward 0 and the last one decreases linearly as a function of log(λ), which is also not the expected behavior. Last, parameter F decreasing as an exponent of λ is investigated (see Figure 2c). A seemingly good behavior happened when F = .5λ -0.4 ; with this value, the curve approximately converges toward a constant value -0.005. With other tunings, either the curves converge toward 0 or decrease linearly. Crossover probability Cr has seemingly no impact on the convergence rate CR. Experiments suggest that for λ large and 1000 iterations the best performing tested variant is Eq. 16 (Cr = .44 and F = .5λ -0.4 ), independently of the dimension. 3.2 Extension: the case F = 0 and larger population sizes Section 3.1 shows some surprisingly good results with F very small. We investigate then the specific case F = 0. When we set F = 0, DE can still perform optimization, through the use of coordinate-wise mutations. This optimization does not converge asymptotically, because it is restricted to a finite set; yet, unless the number of iterations is huge, this is quite effective on the empirically measured convergence rate. This is shown in Figure 3, with a number of iterations N set to 1000. Each result is the average of 15 runs. We present here experiments with larger population size, including F = 0.

Figure 3a displays the Convergence Rate in y-axis and the population size λ in x-axis. Dimension D is fixed to 3. Figure 3b presents experiments of Figure 3a with the Convergence Rate divided by log(λ) in y-axis and the population size λ in x-axis. Hence, if Eq. 3 holds, we should observe a negative plateau if λ is large enough. We notice that F = 0 is a bad parameter choice -even with a moderate number of iterations (1000)-as the convergence rate remains in 0 (the algorithm does not converge). Two parametrizations appear to perform well when the population size is large: (F = .1, Cr = .5) and (F = .5λ -0.4 , Cr = 0.44).

Figure 3c exhibits results when population size and dimension are linked. Dimension D is here λ/2. We see that F = 0 is very satisfactory in such a case, and F = .5λ -0.4 (which is also satisfactory in Figures 3a and3b with fixed dimension) has the same performance in this D = λ/2 setting.

As a result, these experiments show that (i) F = 0 is suitable in the important case λ proportional to the dimension and a fixed number of iterations. (ii) However, F = 0 is strongly suboptimal for a fixed dimension and λ goes to infinity. (iii) F decreasing such as F = .5λ -0.4 independently of the dimension is seemingly a reasonably good idea in all cases. It is not clear whether it succeeds asymptotically for our two parallelization criteria (Eqs. 2 and 3 for λ → ∞), but the behavior was better than other methods as far as we have seen on a simple sphere function problem.

We see on these experiments that for a fixed total budget, when the dimension and/or the population size increase, the effect of the F parameter is much more important than the Cr value.

Adaptive DE

We here test adaptive variants of DE, aimed at choosing adaptively the Cr and F parameters. Experiments are performed with F = 0.5, Cr = 0.9, for DE/rand/1, DE/rand/2, DE/best/1, DE/best/2. The variants JADE and MDE pBX are adaptive and choose their parameters themselves. We refer to [START_REF] Zhang | Jade: adaptive differential evolution with optional external archive[END_REF][START_REF] Islam | An adaptive differential evolution algorithm with novel mutation and crossover strategies for global numerical optimization[END_REF][START_REF] Mallipeddi | Differential evolution algorithm with ensemble of parameters and mutation strategies[END_REF] for more on the parametrization of differential evolution and the adaptive variants. Experimental results are presented in Figure 4 with 50 iterations. The adaptive variants of DE fail in preserving the log(λ) speedup in case of large population size for a fixed D: the convergence rate decreases and is lower than for our non-adaptive variant (see Figure 4a). For D = λ/2, JADE succeeds reasonably well (see Figure 4b). large population sizes, including a parametrization with F = 0, some classical parametrizations, and our proposed formula F = 0.5λ -0.4 . Each point is averaged over 15 runs. These experiments performed in different settings are aimed at validating our proposed formula for large population sizes and checking that it is better than (i) some usual parametrizations (ii) the limit value F = 0. We see that F = 0 is not efficient for a fixed dimension: the algorithm does not converge. When the dimension is huge, then the case F = 0 performs well, whereas one can show mathematically that it does not converge. F prescribed by .5λ -0.4 performs well. Figure 4a shows that adaptive variants of DE can not compete with the standard DE (e.g. DE/best/1), which could not compete with our variants with F = 1 2 λ -0.4 in previous experiments. Figure 4b shows that JADE performs reasonably well when the population size is proportional to the dimension.

EXPERIMENTS: SPEED-UP AND SCALABILITY OF PARALLEL PAR-TICLE SWARM OPTIMIZATION

We now try to reach the previous theoretical bounds in the case of PSO. The problem investigated in this section is the possibility to find parameters of PSO such that we get approximately these two good properties, namely linear scalability (Eq. 2) and logarithmic speed-up (Eq. 3). Since asymptotically everyone dies, we focus on the following nonasymptotic versions:

P RD(λ) = D 100 (log(best f itness at 100 iterations) -log(initial best f itness)) (17) P R lin (λ) = 1 100
(log(best f itness at 100 iterations)

-log(initial best f itness)) with D = λ/2 [START_REF] Mahdad | Fuzzy Controlled Parallel PSO to Solving Large Practical Economic Dispatch[END_REF] where P R stands for Progress Rate. Implicitly, here, fitness values depend on λ, which is the population size.

We do experiments on simple unimodal functions, the sphere function f (x) = x 2 , plus a validation on the cigar function.

The parameters of PSO are w, φg and φp (see Algorithm 2). We evaluate the impact of each parameter. For the sake of comparison, we need a standard PSO [START_REF] Zambrano-Bigiarini | Standard particle swarm optimisation 2011 at cec-2013: A baseline for future pso improvements[END_REF][START_REF] Clerc | Beyond standard particle swarm optimisation[END_REF][START_REF] Bratton | Defining a standard for particle swarm optimization[END_REF] as a baseline. We choose ω = 1/(2 log( 2)) ≃ 0.72, φp = φg = 0.5 + log(2) ≃ 1.2.

In Section 4.1, we investigate the case of large population size and logarithmic speed-up. Population size of the same order as the dimension and linear scalability are studied in Section 4.2.

Asymptotic scalability: large population size, fixed dimension

In this section we check if Eq. 3 holds. The fitness function used in these experiments shown in Figure 5 is

f (x) = x 2 .
We limited the search space to an upper bound bu = 1 and a lower bound b l = -1. The number of generations N is fixed to 100 and the dimension D to 30. In Figure 5, the x-axis shows the natural log of the population size. Figure 5a displays P R30(λ) in order to show progress rate whereas Figure 5b displays P R30(λ)/ log(λ) to check if PSO reaches the optimal asymptotic behavior P R30(λ) = Θ(log(λ)). If it is the case, the curves should show a plateau. Each experiment is repeated 10 times. Figure 5 seems to indicate that there exist different parameter values for which PSO reaches the optimal convergence rate (within a constant factor), at least as far as we can see in dimension 30. This comes a bit as a surprise as other algorithms struggle with such convergence rates (see [START_REF] Teytaud | Log(lambda) Modifications for Optimal Parallelism[END_REF] for efforts aimed at making some classical algorithms verify this optimal behavior) -standard PSO appears to be more naturally parallel than most existing evolutionary algorithms, from the point of view of the convergence rate as the population size grows to infinity, for simple unimodal functions. However, we will see below that for lower dimensions, things are not so nice for standard PSO.

Choosing φg and φp in PSO: Figures 6a, 6b and 6c look at the impact of φg and φp for different ω around these standard values. The x-axis still shows the population size expressed as log(λ) and the y-axis displays the progress rate.

From Figure 6, it appears that φp and φg can only slightly change the convergence rate in our simple unimodal setting, which is in line with the current literature. From several empirical tests, in the following experiments we choose 0.9 instead of the recommended 0.5 + log(2) as it seems to give better results regardless of the ω for the unimodal function under study in this paper.

Choosing w in PSO: Figure 7 evaluates the impact of different ω for a given dimension D = {2, 3, 5, 10}, respectively Figures 7a, 7b, 7c ,7d, to see whether a dynamic ω could yield better results. As such, we are looking for curves that crosses one another. The x-axis still shows the population size expressed as log(λ) and the y-axis displays the progress rate. The first conclusion we can infer from Figure 7 is that the parameter ω clearly varies in relation to λ. Especially ω:0.6 φ p :0.9 φ g :0.9 ω:0.8 φ p :0.9 φ g :0.9 (a) Progress rate P R 30 (λ) ω:0.6 φ p :0.9 φ g :0.9 ω:0.8 φ p :0.9 φ g :0.9 (b) Progress rate divided by log(λ), i.e. P R 30 (λ)/ log(λ)

Figure 5: The progress rate, expressed as variants of Eq. ( 17), for different set of parameters ω, φp, φg, including the values of standard PSO (black dash line with values ω = 0.72 and φp = φg = 1.2). We observe that standard PSO performs well and even seems to reach the optimal convergence rate as λ grows. For small λ a higher value of ω is slightly more appropriate and eventually (around a population λ = 80k) smaller values of ω outperform standard PSO. Standard deviations are of order 0.01. on lower dimensions (Figures 7a,7b, 7c and7d), we observe patterns that emerge. For instance, in dimension D = 5, when the population size log(λ) is between [START_REF] Ardia | Jump-diffusion calibration using Differential Evolution[END_REF][START_REF] Ardia | Differential Evolution with DEoptim: An application to non-convex portfolio optimization[END_REF], the best ω is 0.7. As the population size increases to log(λ) ∈ [START_REF] Ardia | Differential Evolution with DEoptim: An application to non-convex portfolio optimization[END_REF][START_REF] Auger | Convergence results for (1,λ)-SA-ES using the theory of ϕ-irreducible Markov chains[END_REF], the best ω is equal to 0.6. At log(λ) ∈ [START_REF] Auger | Convergence results for (1,λ)-SA-ES using the theory of ϕ-irreducible Markov chains[END_REF][START_REF] Auger | sigma,eta) : quasi-random mutations for evolution strategies[END_REF], the best ω becomes 0.5, etc. The rate at which the different curves crosses one another gives an important indication that there is a dependency of ω over both λ and the dimension D. 8a,8b,8c,8d. The first set of parameters are those provided by standard PSO. As we observe in Figure 7, as the population size λ grows, smaller ω yields better performance. Thus, we can suppose that when λ → ∞, ω = 0 could be the best parameter. This is the second set of parameters with φp = φg = 0.9. The third one ω λ is a set of parameters that depends on both We see that ω smaller (than standard PSO) is better for λ large and ω larger is better for λ small. The impact is, however, moderate in this case (dimension 30). Standard deviations are of order 0.01. λ and the dimension D; this formula is extrapolated from the empirical results of Figure 7. It is given by:

ω λ = max 0.025 + 4 * exp - λ 70 , 0.9 -log λ 2.5 * D . (19)
From Figure 8, it appears that using ω λ instead of standard ω:0 φ p ,φ g :0.9 ω:0.025 φ p ,φ g :0.9 ω:0.05 φ p ,φ g :0.9 ω:0.075 φ p ,φ g :0.9 ω:0.1 φ p ,φ g :0.9 ω:0.2 φ p ,φ g :0.9 ω:0.3 φ p ,φ g :0.9 ω:0.4 φ p ,φ g :0.9 ω:0.5 φ p ,φ g :0.9 (a) Dimension 2, φg = φp = 0.9 ω:0 φ p ,φ g :0.9 ω:0.1 φ p ,φ g :0.9 ω:0.2 φ p ,φ g :0.9 ω:0.3 φ p ,φ g :0.9 ω:0.4 φ p ,φ g :0.9 ω:0.5 φ p ,φ g :0.9 ω:0.6 φ p ,φ g :0.9 ω:0.7 φ p ,φ g :0.9 ω:0.8 φ p ,φ g :0.9 ω:0.9 φ p ,φ g :0.9 ω:1 φ p ,φ g :0.9 (c) Dimension 5, φg = φp = 0.9 ω:0 φ p ,φ g :0.9 ω:0.1 φ p ,φ g :0.9 ω:0.2 φ p ,φ g :0.9 ω:0.3 φ p ,φ g :0.9 ω:0.4 φ p ,φ g :0.9 ω:0.5 φ p ,φ g :0.9 ω:0.6 φ p ,φ g :0.9 ω:0.7 φ p ,φ g :0.9 ω:0.8 φ p ,φ g :0.9 ω:0.9 φ p ,φ g :0.9 ω:1 φ p ,φ g :0.9 (d) Dimension 10, φg = φp = 0.9 In all cases, as λ becomes large, smaller ω yield better results. The effect is more exacerbate on lower dimensions {2, 3, 5, 10}. The different values of ω produce curves that cross one another which indicate that there is a dependency over the population size λ for the parameter ω. We use these curves for designing rules for choosing w as a function of the population size and the dimension (Eq. [START_REF] Mallipeddi | Differential evolution algorithm with ensemble of parameters and mutation strategies[END_REF]. Standard deviations are of order 0.01.

PSO yields results more robust in relation to different dimensions. In lower dimensions D = {2, 3, 5, 10}, the factor of improvement oscillates between 7 and 8. This means that using ω λ instead of standard PSO in a parallel environment is up to 8 times more efficient. The parameters of standard PSO are essentially optimal for a dimension D = 30. It is important to note that ω λ performs at least as good as standard PSO. As the dimension increases D = 50, ω λ becomes again more efficient than standard PSO by a small margin. The special case ω = 0 gives very good results in small dimensions {2, 3, 5}, improving over standard PSO by an average factor of 5, but never equals to ω λ and rapidly becomes irrelevant in higher dimension D > 10 for the tested population sizes λ.

Figure 9 presents the same study as Figure 8 but using another unimodal fitness function, the Cigar. This fitness is defined as f (xi) = x 2 i,0 + 10 6 D d=1 x 2 i,d . The use of ω λ also appears to yield better results in relation to different dimensions. We observe the same behavior for ω = 0 as in Figure 8.

Linear scalability: population size linear in the dimension

In this section we check if, for different variants of PSO, Eq. 2 holds. As a reminder, the fitness function used in these experiments is still f (x) = x 2 . We limit the search space to an upper bound bu = 1 and a lower bound b l = -1. The number of generations N is fixed to 100. The dimension D is given by λ/2. Each experiment is repeated 10 times. Figure 10 presents the results for different parametrizations of ω, φp and φg.

From Figure 10 it appears that again ω λ yields the best progress rate. Standard PSO is following closely (in this setting; not in other settings as mentioned previously). The case where ω = 0 is by far the worst setting of the 3 tested. Thus, we conclude that ω λ seems to be a good parametrization for a parallel PSO environment, in terms of convergence rates.

CONCLUSION

We have studied the speed-up of PSO and DE, as a function of the number of simultaneous function evaluations -a model of synchronous parallelism in which the main computational cost is in the objective function.

We have investigated (i) theoretically the lower bound on runtime for a large class of objective functions and (ii) empirically on the sphere function ( i.e. a model of unimodal well conditioned function) the convergence rate of differential evolution and particle swarm optimization when the population size λ is large. Roughly speaking, the present paper extends [START_REF] Fournier | Lower Bounds for Comparison Based Evolution Strategies using VC-dimension and Sign Patterns[END_REF] in the sense that it gives more precise bounds for DE and PSO, and it extends [START_REF] Teytaud | Log(lambda) Modifications for Optimal Parallelism[END_REF] in the sense that it provides (unproved) formulas for approximately optimal speed-ups for DE and PSO. The improvement compared to standard parametrizations grows seemingly indefinitely for DE, whereas it is rather limited for PSO where the standard parametrizations perform well even with huge values of λ.

We have identified a similarity between PSO and DE: both The performance (progress rate) for different population size λ where the dimension is given by λ/2. Three sets of parameters are tested. The first set is the values given by standard PSO. The second test is given by ω = 0, φg,φp = 0.9. The third and last set is the formula given by ω λ , φg,φp = 0.9. Again, we observe that ω λ gives the best progress rate -though standard PSO is almost equivalent. Standard deviations are of order 0.01. are based on pairwise comparisons. We have shown that DE and PSO:

• can not be faster than -log( xn-x * ) n = Θ(log(λ)), with xn the best individual at iteration n, a fixed dimension, and x * the optimum;

• can not do better than -log( xn-x * ) n = Θ(1) when the dimension is proportional to λ, e.g. D ≃ λ/2. This mathematical lower bound on the runtime is for arbitrary families of fitness functions, provided that the optimum can be anywhere in the domain (i.e. for all x in the domain, the family of functions contains at least one function with optimum x). The detailed bounds are improved when we consider the sphere function.

Interestingly, these rates are exactly the known limits for evolution strategies [START_REF] Fournier | Lower Bounds for Comparison Based Evolution Strategies using VC-dimension and Sign Patterns[END_REF]. In the case of DE or PSO, we have no mathematical proof that these limits are reached. We checked experimentally whether some parametrization can be as fast as allowed by the mathematical bounds. Runtimes are estimated experimentally.

Parametrization of differential evolution for large population sizes. Experimentally, for differential evolution, a good parametrization might be F ≃ .5λ -0.4 and some fixed Cr (Eq. 16). Maybe such a parametrization reaches the bounds above. The improvement compared to standard parametrizations seemingly (on experiments) grows indefinitely as the population size goes to infinity.

Parametrization of particle swarm optimization for large population sizes. The ω parameter should decrease, whereas other parameters can be chosen close to those of standard PSO; for large population sizes, we get a factor 8 on the progress rate in dimension 5 thanks to the parametrization proposed in Eq. 19. Furthermore, standard PSO performs well on the important case λ proportional to the dimension D, as well as when the dimension is 30; we had only minor improvements in these cases. Overall, PSO, under its standard forms, except for moderate dimensions (roughly from 2 to 10), performs well in a parallel setting with a large population size.

We certainly do not claim that Eq. 19 (for PSO) and Eq. 16 (for DE) are universal solutions for choosing parameters in PSO and DE respectively. They are relevant formulas for a clearly defined setting, namely a unimodal well conditioned setting with λ large. Let us discuss limitations. First, this work is limited to unimodal well conditioned functions. Another limitation of this work is the choice of the number of iterations. We checked the convergence rates on limited numbers of iterations, which might be misleading. However, we believe that in a practical world, our population sizes (corresponding to parallelization on hundreds of thousands of cores -quite a big number even for clusters of machines equipped with GPU) and iteration numbers make sense.

Further work. We have shown runtime lower bounds for a class of algorithms including differential evolution and particle swarm optimization, for large population sizes. We have shown empirically that the runtime is close to that bound, for simple problems. This suggests that the runtime lower bounds are tight. Proving a corresponding upper bound on the runtime, matching the lower bound, is the main further work. Designing adaptive rules which provide good rates for a wide class of fitness functions, beyond the simple unimodal scenario in the present paper, is also part of the agenda. Parametrization of ω (resp. F ) proposed for PSO (resp. DE) might be improved, using some information on the budget, or maybe parameters depending on the iteration index; mathematical analysis of these formulas and/or improved variants is an interesting challenge. Another natural extension is to adapt the bound for differential evolution variants with topology (e.g. [START_REF] Weber | Parallel random injection differential evolution[END_REF]).

  g. [28, Th 4, p. 27]), this number of cells is upper bounded by D i=0 λ i ; this is upper bounded by (λ + 1) D , hence the expected result for DE/rand.

3. 1 . 1

 11 Testing D = 5 and large population size λ.

3. 1 . 2

 12 Testing D = λ/2, large population size λ.

FF=0. 5 /F=0. 5 /

 55 independent of λ. Seemingly, the convergence rate reaches a plateau. pow(log(lambda),0.76) Cr=0.86 F=0.5/pow(log(lambda),3.6) Cr=0.68 F=0.5/pow(log(lambda),0.4) Cr=0.44 F=0.5/pow(log(lambda),1.7) Cr=0.96 F=0.5/pow(log(lambda),0.11) Cr=0.81 (b) F inversely proportional to an exponent of log(λ). Results are better than above. pow(lambda,0.76) Cr=0.86 F=0.5/pow(lambda,3.62) Cr=0.68 F=0.5/pow(lambda,0.4) Cr=0.44 F=0.5/pow(lambda,1.67) Cr=0.96 F=0.5/pow(lambda,0.11) Cr=0.81 (c) F inversely proportional to an exponent of λ. Results are better, in particular with Cr = .44 and F = 1 2 λ -0.4 .

Figure 1 :F=0. 5 /

 15 Figure 1: Experiments in dimension 5.Here we plot the convergence rates CR (Eq. 5) of various parametrizations of differential evolution (estimated on 1000 iterations). All experiments are averaged over 15 runs. Standard deviations are very small and not presented. We essentially see that F should decrease when λ increases. Figure3will experiment with a few variants with larger population sizes.

Figure 2 :

 2 Figure 2: Here we plot the convergence rates (Eq. 5) of various parametrizations of differential evolution in the case D = λ/2 with D the dimension. All experiments are averaged over 15 runs. Standard deviations are very small and not presented.Figure 3 presents experiments with a few selected variants with larger population sizes.

Figure 3 presents

 3 experiments with a few selected variants with larger population sizes.

  (a) Dimension D = 3, D × CR as a function of λ. (b) Dimension D = 3, D × CR/ log(λ) as a function of λ.(c) Dimension D = λ/2, D × CR as a function of λ.

Figure 3 :

 3 Figure 3: Convergence rate in the fitness space (Eq. 5) for

  (a) Variants of DE with D = 10 (b) Variants of DE with D = λ/2

Figure 4 :

 4 Figure 4: Experiments aimed at comparing our proposed variants with small values for F as λ → ∞ to adaptive variants of DE. In these experiments, variants of DE are tested with D = 10 and with D = λ/2 respectively, as the population size λ goes to infinity. These results are convergence rates in the fitness space.Figure4ashows that adaptive variants of DE can not compete with the standard DE (e.g. DE/best/1), which could not compete with our variants with F = 1 2 λ -0.4 in previous experiments. Figure4bshows that JADE performs reasonably well when the population size is proportional to the dimension.

  042 φ p :0.32 φ g :0.99 ω:0.72 φ p :1.2 φ g :1.2

  rate:D = 30) / log(λ) ω:0.7 φ p :0.16 φ g :072 φ p :1.2 φ g :1.2

Figure 8

 8 Figure 8 compares 3 different sets of parameters across dimensions D = {2, 3, 5, 50}, respectively Figures8a, 8b, 8c, 8d. The first set of parameters are those provided by standard PSO. As we observe in Figure7, as the population size λ grows, smaller ω yields better performance. Thus, we can suppose that when λ → ∞, ω = 0 could be the best parameter. This is the second set of parameters with φp = φg = 0.9. The third one ω λ is a set of parameters that depends on both

( 2 (Figure 6 :

 26 Figure 6: P R 30 (λ) for different values of ω. In this figure, dif-ferent values of φg and φp for values around those proposed in standard PSO are applied to the study of the (unimodal) sphere function x → x 2 . We see that ω smaller (than standard PSO) is better for λ large and ω larger is better for λ small. The impact is, however, moderate in this case (dimension 30). Standard deviations are of order 0.01.

  .4 φ p ,φ g :0.9 ω:0.5 φ p ,φ g :0.9 (b) Dimension 3, φg = φp = 0.9

Figure 7 :

 7 Figure 7: In dimension D = {2, 3, 5, 10} we check the progress rate P R D (λ) for different values of w.In all cases, as λ becomes large, smaller ω yield better results. The effect is more exacerbate on lower dimensions {2, 3, 5, 10}. The different values of ω produce curves that cross one another which indicate that there is a dependency over the population size λ for the parameter ω. We use these curves for designing rules for choosing w as a function of the population size and the dimension (Eq.[START_REF] Mallipeddi | Differential evolution algorithm with ensemble of parameters and mutation strategies[END_REF]. Standard deviations are of order 0.01.

Figure 8 :

 8 Figure 8: In dimension D = {2, 3, 5, 50} respectively, we check the progress rate P R D (λ) for different values of w. At dimension D = 2, the use of ω λ over standard PSO for parallelization yield an improvement by a factor of 7. At dimension D = 3, ω λ improves over standard PSO by a factor of 8. At dimension D = 5, the factor is also 8. At dimension D = 10 (unpresented), the factor is 7.5. The case of dimension D = 30 (unpresented) provides no improvement, standard PSO is equal to ω λ . At dimension D = 50, there is a small advantage to use ω λ over standard PSO. Standard deviations are of order 0.01.

Figure 9 :

 9 Figure 9: In dimension D = {2, 3, 5, 50} respectively, we check the progress rate P R D (λ) for different values of w. At dimension D = 2, the use of ω λ over standard PSO for parallelization yield an improvement by a factor of more than 5. At dimension D = 3, ω λ improves over standard PSO by a similar factor. At dimension D = 5, the factor is around 9. At dimension D = 10 (unpresented), the factor is 2.5. The case of dimension D = 30 (unpresented) provides no improvement, standard PSO is essentially equal to ω λ . At dimension D = 50, there is a small advantage to use ω λ over standard PSO. Standard deviations are of order 0.01.

Figure 10 :

 10 Figure 10: The performance (progress rate) for different popu-

  Algorithm 1 Pseudo-code of basic differential evolution. The first loop on the population contains the fitness evaluations, supposed to be the key part of the computational cost; it is fully parallel. For j ∈ {1, . . . , D}, xi,j denotes the j th coordinate of a vector xi ∈ R D . U(a, b) D stands for an uniform random variable in the domain (a, b) D .

	Require: D: Dimension, N : Number of generations, F ∈
	[0, 2]: Differential weight, Cr ∈ [0, 1]: Crossover probabil-
	ity, λ: Population size, best ∈ {0, 1}
	for each i ∈ {1, . . . , λ} do
	Initialize pi ∼ U (-500, 500) D
	Compute f (pi)
	end for
	if best = 1 then
	g ← arg min	f (pi)
	p 1 ,...,p λ
	end if
	n = 1
	while n<N do
	for i ∈ {1, . . . , λ} parallel loop do
	Randomly draw a, b and c distinct
	in {1, . . . , i -1, i + 1, . . . λ} Define p ′ i ← pa + F (p b -pc)
	Randomly draw R ∈ {1, . . . , D}
	for each j ∈ {1, . . . , D} do
	p

  ∈ P, xi ∈ R D is in the search space and the parti-cle has a velocity vi ∈ R D . Let pi be the best known position of the particle i and g be the best found solution so far. Algorithm 2 presents a basic PSO following this notation.Algorithm 2 Basic PSO. The first loop on the population contains the fitness evaluations, supposed to be the key part of the computational cost; it is fully parallel. Notations are the same as in Algorithm 1, ie. xj (resp. xi,j) is the j th component of vector x (resp. xi).

	end if end if end for n ← n + 1 end while if best = 0 then g ← arg min p 1 ,...,p λ f (pi) end if Return g the present paper; they combine DE and simulated anneal-ing. [36] considers a ring topology and focuses on avoiding premature convergence. [38] considers the problem of multi-ple populations from the point of view of global convergence and maintaining diversity, while we, in the present paper, consider mainly (at least for experimental rates -our lower bounds on runtime are for arbitrary functions with unique global optimum) local convergence. Consistently with [1], we consider the performance as a function of the population size; we also prove lower bounds and discuss the optimal parametrization of DE with large population size. Particle Swarm Optimization. Particle Swarm Opti-mization (PSO) is another classical evolutionary algorithm for continuous domain. Let P be a population of particles; for each i Require: D: Dimension, N : Number of generations, P: Population, b l : Lower bound of the search space, u l : Up-per bound of the search space, w,φp,φg: Parameters for each particle i ∈ P do Initialize its position xi from a uniform distribution: xi ∼ U (b l , bu) D Initialize its best known position: pi ← xi Initialize its velocity: vi ∼ U (-|bu -b l |, |bu -b l |) D end for Initialize best solution: g ← arg min x i , i∈P f (xi) n = 1 while n < N do for each i ∈ P parallel loop do Pick rp, rg ∼ U(0, 1) for each j ∈ {1, . . . , D} do vi,j ← wvi,

j + φprp(pi,j -xi,j) + φgrg(gj -xi,j) end for


David Auger
email: david.auger@prism.uvsq.fr

Alcaap

Jialin Liu
email: jialin.liu@inria.fr

Sylvie Ruette
email: sylvie.ruette@math.u-psud.fr

David L Saint-Pierre
email: davidlstp@gmail.com

Olivier Teytaud
email: olivier.teytaud@inria.fr

Soon Cheng

Tu Bao Ong

Ho

Sparse Binary Zero-Sum Games

Keywords: Sparsity, bandit algorithms, zero-sum matrix games

Solving zero-sum matrix games is polynomial, because it boils down to linear programming. The approximate solving is sublinear by randomized algorithms on machines with random access memory. Algorithms working separately and independently on columns and rows have been proposed, with the same performance; these versions are compliant with matrix games with stochastic reward. (Flory and Teytaud, 2011) has proposed a new version, empirically performing better on sparse problems, i.e. cases in which the Nash equilibrium has small support. In this paper, we propose a variant, similar to their work, also dedicated to sparse problems, with provably better bounds than existing methods. We then experiment the method on a card game.

Introduction

Solving Nash equilibria (NE) of matrix games is important by itself, e.g. for some economical models; as a building block for Monte-Carlo Tree Search algorithms [START_REF] Flory | Upper confidence trees with short term partial information[END_REF] when simultaneous actions are involved; and for robust stochastic optimization [START_REF] Lambert | A fictitious play approach to large-scale optimization[END_REF]. Bandits algorithms [START_REF] Lai | Asymptotically efficient adaptive allocation rules[END_REF][START_REF] Michael | A sublinear-time randomized approximation algorithm for matrix games[END_REF][START_REF] Auer | Gambling in a rigged casino: the adversarial multi-armed bandit problem[END_REF][START_REF] Audibert | Minimax policies for adversarial and stochastic bandits[END_REF] are tools for handling exploration vs exploitation dilemma that are in particular useful for finding NE [START_REF] Michael | A sublinear-time randomized approximation algorithm for matrix games[END_REF][START_REF] Audibert | Minimax policies for adversarial and stochastic bandits[END_REF], with applications to two-player games [START_REF] Kocsis | Bandit based Monte-Carlo planning[END_REF][START_REF] Flory | Upper confidence trees with short term partial information[END_REF]. To the best of our knowledge, the only paper in which sparsity in NE is used for accelerating bandits is [START_REF] Flory | Upper confidence trees with short term partial information[END_REF].

Section 2 presents the framework of sparse bandits and related algorithms. Section 3 mathematically analyzes sparse bandits. Section 4 presents experiments on a card game and Section 5 concludes.

Algorithms for sparse bandits

In this section we define the problem and our proposed approach. Section 2.1 presents below the notion of NE in matrix games. Section 2.2 defines a classical bandit algorithm framework that aims at solving such a problem and Section 2.3 introduces two proposed algorithms adapted for sparse problems.

Throughout the paper, [[a, b]] denotes {a, a + 1, a + 2, . . . , b -1, b} and lcm(x 1 , . . . , x n ) denotes the least common multiple of x 1 , . . . , x n , where a, b, x 1 , . . . , x n are integers.

Matrix games and Nash equilibria

Consider a matrix M of size K × K with values in {0, 1} (we choose a square matrix for short notations, but the extension is straightforward). Player 1, the row player, chooses an action i ∈ [[1, K]] and player 2, the column player, chooses an action j ∈ [[1, K]]; both actions are chosen simultaneously. Then, player 1 gets reward M i,j and player 2 gets reward 1 -M i,j . The game therefore sums to 1 (we consider games summing to 1 for convenience, but 0-sum games are equivalent).

A NE is a pair (x * , y * ) (both in [0, 1] K and summing to 1) such that if i is distributed according to the distribution x * (i.e. i = k with probability x * k ) and if j is distributed according to the distribution y * (i.e. j = k with probability y * k ) then none of the players can expect a better average reward by changing unilaterally its strategy, i.e. ∀x, y, (x * ) T M y ≥ (x * ) T M y * ≥ x T M y * .

There might be several NE, but the value of the game M , given by v = x * M y * , remains the same. Moreover, we define an -NE as a pair (x, y) such that inf y x T M y > (x * ) T M y * -and sup

x (x ) T M y < (x * ) T M y * + .

Bandit algorithms

A state-of-the-art bandit algorithm to approximate a NE (which includes matrix games) is EXP3 [START_REF] Auer | Gambling in a rigged casino: the adversarial multi-armed bandit problem[END_REF]. Here we present the version used in [START_REF] Audibert | Minimax policies for adversarial and stochastic bandits[END_REF].

At iteration t ∈ [[1, T ]],
our version of EXP3 proceeds as follows (this is for one of the players; the same is done, independently, for the other player):

• At iteration 1, S is initialized as a null vector.

• Action i is chosen with probability p(i) = α t /K + (1 -α t ) × exp(η t S i )/ j exp(η t S j )
for some sequences (α) t≥1 and (η t ) t≥1 , e.g. in [START_REF] Bubeck | Regret analysis of stochastic and nonstochastic multi-armed bandit problems[END_REF]) α t = 0 and η t = log(K)/(tK) or in [START_REF] Audibert | Minimax policies for adversarial and stochastic bandits[END_REF]

η t = min( 4 5 log(K) T K , 1 K ) and α t = Kη t .
• Let r be the received reward.

• Update S i : S i ← S i + r/p(i) (and S j for j = i is not modified). This algorithm, as well as its variants, converges to the NE as explained in [START_REF] Audibert | Minimax policies for adversarial and stochastic bandits[END_REF][START_REF] Bubeck | Regret analysis of stochastic and nonstochastic multi-armed bandit problems[END_REF]) (see also [START_REF] Michael | A sublinear-time randomized approximation algorithm for matrix games[END_REF][START_REF] Auer | Gambling in a rigged casino: the adversarial multi-armed bandit problem[END_REF]).

2.3. Rounded bandits. [START_REF] Flory | Upper confidence trees with short term partial information[END_REF] proposes a generic way of adapting bandits to the sparse case. Basically, the bandit runs as usual, and then the solution is pruned. We here propose a variant of their algorithm, as explained in Alg. 1. Our definition of sparsity does not assume that the matrix M is sparse; rather, we assume that the two vectors of the NE x * and y * have support of moderate size k, i.e. k << K. This assumption certainly holds for many games.

Algorithm 1 RBANDIT, the rounded bandit algorithm. δ ∈]0, 1[ is a parameter.

Input: A K × K matrix M , defined by a mapping (i, j) → M i,j .

Sparsity assumption: We assume that the NE (x * , y * ) is unique and sparse, i.e. x * and y * have support of size ≤ k. Define p = 1/(4ck k ), where c = lcm(1, 2, . . . k k /2 ).

Let q be such that on K × K matrices the bandit algorithm with q iterations provides a p-NE with probability 1 -δ. Run the bandit algorithm for q iterations. Let (x, y) be the resulting approximate NE.

Approximate by multiples of 1/c as follows:

x i = cx i + 1 2 /c y i = cy i + 1 2 /c. Renormalize: x = x / i x i ; y = y / i y i .
Output (x , y ) as an approximate NE.

The version in [START_REF] Flory | Upper confidence trees with short term partial information[END_REF] is based on truncations (components less than a given are removed), instead of roundings; such a version is presented in Alg. 2 (the exact solving after truncation was not in [START_REF] Flory | Upper confidence trees with short term partial information[END_REF]). The drawback is that for an exact solving we need a second step of exact solving, but we will see that the resulting complexity is better than the complexity of the rounded bandit (Alg. 1) because in case of rounding we need a higher precision in the bandit part. [START_REF] Flory | Upper confidence trees with short term partial information[END_REF] gives no proof and does not provide any tool for choosing c. The aim of this paper is to show that our versions (Alg. 1 and 2, i.e. rounded/truncated bandits) find exact NE, faster than existing methods when the solution is sparse.

Algorithm 2 TBANDIT, the truncated bandit algorithm. δ ∈]0, 1[ is a parameter.

Input: A K × K matrix M , defined by a mapping (i, j) → M i,j .

Sparsity assumption: We assume that the Nash equilibrium (x * , y * ) is unique and sparse, i.e. x * and y * have support of size ≤ k. Define p = 1/(4ck k ), where c = k k/2 . Let q be such that on K × K matrices the bandit algorithm with q iterations provides a p-NE with probability 1 -δ.

Run the bandit algorithm for q iterations. Consider the game restricted to rows i such that x i > 1/(2c) and to columns j such that y j > 1/(2c). Let (x , y ) be the exact Nash equilibrium of this restricted game (computed in polynomial time by linear programming). Output (x , y ) as approximate Nash equilibrium of the complete game (complete with 0's for missing coordinates).

Analysis

This section is devoted to the mathematical analysis of sparse bandit algorithms. Section 3.1 introduces the useful notations. Section 3.2 shows some properties of supports of Nash equilibria. Section 3.3 gives some useful results on denominators of rational probabilities involved in Nash equilibria. Section 3.4 presents stability results (showing that in sparse bandits, good strategies are also close to the Nash equilibrium). Section 3.5 concludes by properties on sparse bandits algorithms.

Terminology

We use the classical terminology of game theory. We consider a Matrix Game M of size

K × K as above. A pure strategy is an index i ∈ [[1, K]].
A mixed strategy is a vector of size K with non-negative coefficients summing to 1. Let e i be the vector (0, 0, . . . , 0, 1, 0 . . . , 0) T with a one only at the i th position; by a slight abuse of notation we will use this notation independently of the dimension of the vector (i.e. e i can be used both in R 10 and R 50 ).

Let ∆ denote the set of probability vectors, that is, ∆ = {y : j y j = 1 and ∀j, y j ≥ 0}; this implicitly depends on the dimension of the vectors, we do not precise it since there will be no ambiguity. The support of a vector y ∈ ∆ is the set of indices j such that y j > 0. For short, sup x (or sup y , inf x , inf y equivalently) means supremum on

x ∈ ∆. The value of y ∈ ∆ for M is V (y) = sup x x T M y.
Recall that v denotes the value of the game M , that is, it satisfies

∀x, y ∈ ∆, x T M y * ≤ v ≤ (x * ) T M y, (1)
and v = (x * ) T M y * = V (y * ) if (x * , y * ) is a Nash equilibrium.

Supports of Nash equilibria

Here we consider general matrices A, M with real coefficients. The following lemma is well known (see [START_REF] Gale | Linear programming and the theory of games[END_REF], Lemma 1 page 318)).

Lemma 1 (Farkas Lemma) There exists y ≥ 0 satisfying Ay = b if and only if there is no x such that A T x ≥ 0 and x T b < 0.

The following lemma is adapted from [START_REF] Dantzig | Linear Programming 2: Theory and Extensions[END_REF]. We give the proof for the sake of completeness.

Lemma 2 Let I be the set of column indices i such that for all optimal solutions x * of the row player we have (x * ) T M e i = v. Then there exists an optimal solution y * for the column player whose support is exactly I.

Proof First, we show that it is sufficient to prove that for any i ∈ I there is an optimal solution y i ∈ ∆ such that y i i > 0. Then one considers any strictly convex combination of the y i , for instance y * = 1 |I| i∈I y i . The vector y * ∈ ∆ has support including I (because y i i > 0). On the other hand, y * has support included in I (because by construction any optimal solution has support included in I).

So, it is sufficient to show that for any i ∈ I there is an optimal y i such that y i i > 0. We now prove this. Without loss of generality fix i = 1 ∈ I. Let us suppose that no optimal solution y of the column player has a positive coordinate y 1 > 0. In other words, the system

       i y i = 1 M y ≤ v1 y 1 > 0 y ≥ 0
has no solution, where 1 is the vector with all coefficients equal to 1. Equivalently, this means that the following system has no solution

       i y i = 1 (M -v1 K×K )y ≤ 0 y 1 > 0 y ≥ 0
where 1 K×K is the K × K matrix with all coefficients equal to 1. Introducing the variable z = y y 1 and a slack variable w of size K × 1, this is equivalent to saying that the system

   (M -v1 K×K )z + w = 0 z 1 = 1 z, w ≥ 0
has no solution. By Lemma 1, applied with the concatenation (z, w) as y, b = (0, 0, . . . , 0, 1) and

A = M -v1 K×K Id K 1 0 • • • 0 0 ,
we deduce the existence of a vector x and a real number such that

       x T M •1 -v i x i + ≥ 0 x T (M -v1 K×K ) ≥ 0 x ≥ 0 < 0
where M •1 denotes the first column of M . By the first equation above, x is not zero. Thus we can normalize x to get a vector in ∆, and we infer the existence of an optimal strategy x * ∈ ∆ for the row player such that

(x * ) T M •1 > v;
this implies that we cannot have 1 ∈ I, a contradiction.

Corollary 3 If M admits a unique Nash Equilibrium (x * , y * ) with support J × I then:

• for all i ∈ I and j ∈ J we have

(x * ) T M e i > v and e T j M y * < v; (2)
• the submatrix M of M with rows and columns respectively in J and I has a unique Nash Equilibrium which is the projection of x * , y * on J × I.

Proof The first part is a consequence of Lemma 2. Indeed, if (x * ) T M e i = v with i / ∈ I, then there exists an optimal solution y' whose support contains i, which contradicts the uniqueness of y * . Thus (x * ) T M e i > v by Eq. ( 1). The statement for j / ∈ J is symmetric. For the second part, the projection is clearly a Nash equilibrium. Then, suppose that there is another Nash equilibrium for M and let (x , y ) be the only 2K vector whose projection on J × I is equal to this other equilibrium (in other words add zero coordinates for i ∈ I and j ∈ J).

Consider now (1 -t)y * + ty with t > 0 and a row index j; then if j ∈ J e T j M ((1 -t)y * + ty ) = (1 -t)v + tv = v and if j ∈ J, then by the first part of this corollary the left part is at most v for t small enough. Since we have a finite number of rows this implies that by choosing t small enough we obtain a vector (1 -t)y * + ty which is another optimal solution for the column player in M , which contradicts the uniqueness.

Denominators of Nash equilibria

Consider a matrix M , with coefficients in {0, 1}, of size

k 1 × k 2 with k 1 ≥ 2 and k 2 ≥ 2.
Lemma 4 Assume that the Nash equilibrium (x * , y * ) is unique and that ∀i, j, x * i > 0 and y * j > 0. a) Then k 1 = k 2 and x * and y * are rational vectors which can be written with common denominator at most

k k/2 with k = k 1 = k 2 . b) Moreover, for all x, y ∈ ∆, x T M y * = (x * ) T M y = (x * ) T M y * = v.
Proof The Nash equilibrium y * verifies the following properties:

• the sum of the probabilities of all strategies for the "column" player is 1, i.e. i y * i = 1.

(3)

• the expected reward for the "row" player playing strategy i against y * is independent of i, i.e. ∀i,

j M i,j y * j = j M 1,j y * j . ( 4
)
If there is another solution y = y * to Eqs. ( 3), ( 4), then y = y * -α(y -y * ) is another strategy for the column player. If α is small enough, then y ≥ 0 and y i = 1; it is a correct mixed strategy, and its value is x * M y * -αx * M (y -y * ) which is less than or equal to x * M y * if α has the same sign as x * M (y -y * ). This contradicts the uniqueness of y * as a Nash strategy for the column player.

As a consequence, Eqs. ( 3) and ( 4) are a characterization of the unique Nash equilibrium. Thus y * can be computed by solving Eqs. ( 3) and (4); this is a linear system Zy * = (1, 0, 0, . . . , 0) with one single solution and Z a matrix with k 1 rows and k 2 columns, and where Z has values in {-1, 0, 1}. The solution is unique, therefore k 1 = k 2 and Z is invertible.

Z -1 can be computed as

Z -1 = 1 det(Z) (cofactor(Z)) T where (cofactor(Z)) ij = (-1) i+j det (Z i j ) i =i, j =j .
The matrix Z has coefficients in {-1, 0, 1}; therefore, by Hadamard's maximum determinant problem: [START_REF] Hadamard | Résolution d'une question relative aux déterminants[END_REF], see e.g. (Brenner and Cummings, 1972, p 626)). Moreover, the matrix cofactor(Z) has integer coefficients. This concludes the proof of the fact that y * is rational with denominator D = | det(Z)| at most k k/2 , where k = k 1 = k 2 . The same arguments using x * instead of y * show that x * can also be written as a rational with the same denominator D.

| det(Z)| ≤ k k/2 ((
To show b), notice that Eq. ( 4) can be rewritten as follows: ∀i, (M y

* ) i = (M y * ) 1 . If x ∈ ∆, then x T M y * = i x i (M y * ) i = i x i (M y * ) 1 = (M y * ) 1 because i x i = 1. Thus x T M y * is independent of x ∈ ∆, which implies that x T M y * = (x * ) T M y * = v. By symmetry, one also has: ∀y ∈ ∆, (x * ) T M y = (x * ) T M y * = v.
Please note that k k/2 is known nearly optimal for matrices with coefficients in {0, 1} (by Hadamard's work) for any k of the form 2 m . Also there are examples of matrices for which V (y) = V (y * ) + 1 | det(Z)| y -y * for y arbitrarily close to y * .

Stability of Nash equilibria

In the general case of a zero-sum game, two mixed strategies can be far from each other, whenever both of them are very close, in terms of performance, to the performance of the (assumed unique) Nash equilibrium. However, with a matrix M with values in {0, 1}, this is not true anymore, as explained by the two lemmas below.

Lemma 5 Let k ≥ 2. Consider a k × k matrix M with elements in {0, 1} such that the Nash equilibrium (x * , y * ) is unique and no pure strategy has a null weight. Then for all y ∈ ∆ we have

V (y) ≥ V (y * ) + 1 k k/2 y -y * ∞ .
(5)

Proof By convexity of V , it is sufficient to prove that min

u; i u i =0, u ∞=1 lim t→0, t>0 V (y * + tu) -V (y * ) t ≥ 1/k k/2 .
This is equivalent to

min u; i u i =0, u ∞=1 max i M i .u ≥ 1/k k/2 , ( 6
)
with M i the i th row of M as previously.

Let u be a vector in which the minimum is reached (it exists by compactness). Since u ∞ = 1, there exists i 0 such that | u i 0 | = 1. Let us assume without loss of generality, that u i 0 = 1. The proof is the same if u i 0 = -1. Thus, in Eq. ( 6), we can restrict to the vectors u such that

u i 0 = 1, i u i = 0 and ∀i ∈ {1, . . . , k}, -1 ≤ u i ≤ 1.
This is indeed a linear programming problem, as follows:

min u∈R k , w∈R w under constraints ∀i ∈ {1, . . . , k}, -1 ≤ u i ≤ 1, ∀i ∈ {1, . . . , k}, M i .u ≤ w, u i 0 = 1, 1≤i≤k u i = 0.
It is known that when a linear problem in dimension k + 1 has a non infinite optimum, then there is a solution (u, w) with k + 1 linearly independent active constraints.

Let us pick such a solution u. It is solution of k + 1 linearly independent equations of the form either

u i = 1, or u i = -1, or M i .u = w, or i u i = 0. Let us note the system ∀i ∈ P, u i = 1, ∀i ∈ N, u i = -1, ∀i ∈ H, M i .u = w, 1≤i≤k u i = 0,
where P, N, H are the subsets of {1, . . . , k} where the corresponding constraints are active.

• u < 1/(4k k/2 k k ) for the TBANDIT algorithm.

• u < 1/(4Ck k ) for the RBANDIT algorithm. By Theorem 7, this implies that

• x -x * ∞ ≤ 2uk k < 1/(2k k/2 ) (idem for y -y * ∞ ) for TBANDIT; • x -x * ∞ ≤ 1/(2C) (idem for y -y * ∞ )
for RBANDIT. Then:

• Truncated algorithm: all non-zero coordinates of x * are at least 1/k k/2 and |x * i -x i | < 1/(2k k/2 ) with probability ≥ 1 -δ (and the same for y * , y); so with probability 1 -δ the Nash equilibrium (x , y ) of the reduced game is the solution (x * , y * ) (after filling missing coordinates with 0).

• Rounded algorithm: the denominator of the coordinates of x * , y * is a divisor of C, so with probability 1 -δ,

Cx -Cx * ∞ < 1/2, Cy -Cy * ∞ < 1/2,
and Cx * and Cy * are integers. So x * = Cx + 1 2 /C; RBANDIT finds the exact solution with probability ≥ 1 -δ.

For example, if using the Grigoriadis & Khachiyan algorithm (Grigoriadis and Khachiyan, 1995), or variants of EXP3 [START_REF] Audibert | Minimax policies for adversarial and stochastic bandits[END_REF], one can ensure precision u with fixed probability in time

• O K log K 1 u = O K(log K)(k 3k
) for the truncated version;

• O K log K 1 u = O(K(log K)k 2k exp(2k k/2 (1 + o(1))))
for the rounded version. Then, we get, after rounding (rounded version) or after truncating and polynomial-time solving (truncated version), the exact solution y * with fixed probability and time

• O K log K • k 3k + poly(k) for the truncated algorithm (Alg. 2); • O K log K • k 2k exp(2k k/2 (1 + o(1)
)) for the rounded algorithm (Alg. 1).

The truncated algorithm (Alg. 2) is therefore better.

Experiments

We work on the Pokemon card game. More precisely, we work on the metagaming part, i.e. the choice of the deck; the ingaming is then handled by a simulator with exact solving. The source code is freely available at http://www.lri.fr/ ~teytaud/games.html.

At first a normal EXP3 is executed using our empirically tuned formula to compute the probability of arm i and we normalize the probabilities if needs be. After T iterations, the TEXP3 takes the decision whether an arm is part of the NE is based upon a threshold ζ as explained in Alg. 3. Alg. 3 is based on Alg. 2, with constants adapted empirically and without the exact solving at the end.

p(i) = 1 + c -1 √ t -1 ×   1/(c × √ t) + (1 -1/ √ t) × exp(S i / √ t)/ j exp(S j / √ t)   (11)
Algorithm 3 TEXP3, an algorithm proposed in [START_REF] Flory | Upper confidence trees with short term partial information[END_REF] used in these experiments.

Input: A K × K matrix M , defined by a mapping (i, j) → M i,j . A number T of rounds.

Run EXP3 with Eq. ( 11), which provides an approximation (x, y) of the Nash equilibrium.

Define:

ζ = max a∈K (T x a ) 0.7 T x i = x i if x i ≥ ζ and x i = 0 otherwise. x i = x i / j∈{1,...,K} x j . Define: ζ = max a∈K (T y a ) 0.7 T y i = y i if y i ≥ ζ and y i = 0 otherwise. y i = y i / j∈{1,...,K} y j .
Output x and y as approximate Nash equilibrium of the complete game.

Figures 1(a), 1(c) and 1(e) show the performance of TEXP3 playing against EXP3 for different values of c and Figures 1(b), 1(d ) and 1(f ) present the performance of TEXP3 and EXP3 when playing against the random uniform baseline; the probability distribution obtained by EXP3 and TEXP3 after T iterations of EXP3 and (for TEXP3) after truncation and renormalization are used against random. To ensure that no player gains from being the first player, we make them play as both the row player and the column player and we display the result. Each point in the Figure consists in the means from 100 independent runs.

In all Figures, TEXP3 provides a consistent improvement over EXP3. Even in Figure 1(b), where EXP3 seems relatively weak against the random baseline, TEXP3 manages to maintain a performance similar to the ones in Figure 1(d ) and 1(f ).

Conclusion

(Grigoriadis and [START_REF] Michael | A sublinear-time randomized approximation algorithm for matrix games[END_REF][START_REF] Auer | Gambling in a rigged casino: the adversarial multi-armed bandit problem[END_REF][START_REF] Audibert | Minimax policies for adversarial and stochastic bandits[END_REF] are great steps forward in zero-sum matrix games, and beyond. They provide algorithms solving K × K matrix games, with precision and for a fixed confidence, in time O(K log K/ 2 ).

• Another extension is the case of stochastic versions of matrix games: i.e. if player 1 plays i and player 2 plays j, then the reward is a random variable M i,j , with expectation M i,j , independently sampled at each trial. This does not change the result provided that M i,j verifies a condition as above, i.e. a common denominator bounded by some known integer.

Second, we point out that in the real world, using sparsity provides substantial benefit. For example, [START_REF] Flory | Upper confidence trees with short term partial information[END_REF][START_REF] Chou | Strategic choices: Small budgets and simple regret[END_REF] get benefits on a real-world internet card game using sparsity; their algorithms are similar to ours (for finding approximate results), but without the final step for finding an exact solution. We reproduce here their experiments in yet another game (Section 4) and get improvements far better than the theoretical bound. This suggests that, beyond the guaranteed improvement, there is a large improvement in many practical cases.

Further work

K is usually huge; an algorithm linear in K might be not practical. Algorithms with complexity O( √ K), O(log(K)) might be useful; we know that such algorithms can not do any approximation of a Nash equilibrium without additional assumptions, and maybe some regularity assumptions on the strategies of the row player and on the strategies of the column player are required for this. Such a framework (K huge but regularity assumptions on row strategies and regularity assumptions on column strategies) looks like a relevant framework for applications.

A distinct further work is the case in which we have no upper bound on the sparsity parameter k.

Finally, two assumptions might be partially removed: uniqueness of the Nashequilibrium, and 0-sum nature of the game. In particular, extending to the case of a unique subgame perfect equilibrium looks like a promising direction.

Figure 1 :

 1 Figure 1: Performance (%) in terms of budget T for the game of Pokemon using 2 cards. The left column shows TEXP3 playing against EXP3 for different values of c.The right column shows EXP3 and TEXP3 playing against the random uniform baseline. We tested a wide range of values for c and TEXP3 performs better than EXP3 regardless of c.

c 2014 D. Auger, J. Liu, S. Ruette, D.L. Saint-Pierre & O. Teytaud.

Acknowledgment

We are grateful to the BIRS seminar on Combinatorial Game Theory, to the Dagstuhl seminar on the Theory of Evolutionary Algorithms, and to the Bielefeld seminar on Search Methods, to National Science Council (Taiwan) for NSC grants 99-2923-E-024-003-MY3 and NSC 100-2811-E-024-001. We are grateful to C. Teytaud for implementing the Pokemon deck.

We can remove w by setting M j .u = M i .u for some fixed i ∈ H and all j ∈ H \ {i}. Then, u is solution of a system of k equations in dimension k, with coefficients in {-1, 0, 1}.

We use the same trick as in the proof of Lemma 4; u is solution of a system of k linear equations with all coefficients in {-1, 0, 1}. Therefore all coordinates of u are rational numbers with a common denominator D ≤ k k/2 . Then M i . u = M i .u has a denominator D ≤ k k/2 and is positive; therefore M i . u ≥ 1/k k/2 . This proves the expected result.

Lemma 6 Consider M a K × K matrix with coefficients in {0, 1} and assume that the Nash equilibrium (x * , y * ) is unique. Let J be the support of y * and k = #J. Then

Proof According to Lemma 4, v and the coefficients of x * , y * are multiple of some constant c ≥ 1

. Thus, for every j, (x * ) T M e j is also a multiple of c. Fix j ∈ J. By Corollary 3,

Combining Lemmas 5 and 6 yields the following Theorem 7 Consider a matrix M of size K × K and with coefficients in {0, 1}. Assume that there is a unique Nash equilibrium (x * , y * ). Let k be the size of the supports of x * , y * . Then

. Let J be the support of y * . For every y ∈ ∆, one can write y = ay + by , with a = j∈J y i ∈ [0, 1], a + b = 1 and y , y ∈ ∆ satisfying: ∀j ∈ J, y j = 0 and ∀j ∈ J, y j = 0. For every index i, one has

By Lemma 4(b), (x * ) T M y = v; and by Lemma 6, (x * ) T M e j ≥ v + c for all j / ∈ J. thus

By Eq. ( 8), either δ = bδ 2 ≤ b, or δ = a(y -y * )-by * ∞ . If δ ≤ b, the result is given by Eq. ( 9) because c < 1 and hence c 2 ≤ c. From now on, assume that δ = a(y

We split the end of the proof into two cases.

Case 1: b ≥ cδ/2. Then V (y) -v ≥ cb by Eq. ( 9)

≥ c 2 δ/2 by assumption on b which gives the expected result in case 1.

Case 2: b < cδ/2. Since y has the same support as y * , there exists x ∈ ∆ with the same support as x * such that x M y -v ≥ cδ 1 by Lemma 5. Moreover,

Hence,

-vb (using Eq. ( 10))

Since 1 -c ≥ c, we get the expected result in case 2.

Application to sparse bandit algorithms

Consider a matrix M as in Theorem 7. By Lemma 4, (x * , y * ) can be written with a common denominator at most k k/2 . Define C = lcm(1, 2, 3, . . . , k k/2 ). By the prime number theorem, it is known

We discuss in parallel the truncated bandit algorithm (Alg. 2) and the rounded bandit algorithm (Alg. 1), as follows.

By construction of the algorithms, with probability 1 -δ, the bandit algorithm finds a u-Nash equilibrium (x, y), for 1. See details in http://mathworld.wolfram.com/LeastCommonMultiple.html.

As noticed in [START_REF] Michael | A sublinear-time randomized approximation algorithm for matrix games[END_REF], this has the surprising property that the complexity is sublinear in the size of the matrix, with a fixed risk.

We here show that, with coefficients in {0, 1}, if there is a unique sparse Nash equilibrium with support of size k for each player, then this bound can be reduced to K log K • k 3k , with no precision parameter (we provide an exact solution), with a fixed confidence 1 -δ:

• the dependency in K is the same as in (Grigoriadis and Khachiyan, 1995);

• there is no dependency in .

Practical relevance of this work

We want here to discuss the practical relevance of our results; two aspects are (i) the existence of very sparse problems, and (ii) the possible implementation of real world algorithms inspired by our results.

The first point is the existence of very sparse problems. We have seen that the sparsity level that we need, for our algorithm to outperform the state of the art for exact solutions, is k << log(K)/ log(log(K)). Obviously, one can design arbitrarily sparse zero-sum matrix games. Also, in real games, the quantity of moves worth being analyzed is usually much smaller than the number of legal moves; it is difficult to quantify this numerically for existing games as finding the meaningful strategies requires a lot of expertise and is difficult to quantify. For theoretical games, the classical centipede [START_REF]Games of perfect information, predatory pricing, and the chain store[END_REF] game does not have a unique Nash equilibrium; but if we remove strategies which make no sense (i.e. we do not consider the variants of a strategy after a move which finishes the game) the centipede game is highly sparse -only one strategy is in the Nash equilibrium, the one which immediately defects. The centipede game is not a zero-sum game, but zero-sum variants exist with the same sparsity property [START_REF] Fey | An experimental study of constant-sum centipede games[END_REF].

We also provide experimental results which show the relevance of this work in a realworld case. In particular, the TEXP3 modification, with its default parametrization from [START_REF] Flory | Upper confidence trees with short term partial information[END_REF], performs better than the EXP3 counterpart for all tested values of the parameters. TEXP3 is not new, but the formal proof given here is new.

Extensions

The algorithm that we propose is based on the rounding of the solution given by EXP3 [START_REF] Auer | Gambling in a rigged casino: the adversarial multi-armed bandit problem[END_REF], or Grigoriadis' algorithm [START_REF] Michael | A sublinear-time randomized approximation algorithm for matrix games[END_REF], or INF [START_REF] Audibert | Minimax policies for adversarial and stochastic bandits[END_REF]. Our algorithm works thanks to Theorem 7; any matrix game such that Theorem 7 is valid and such that Grigoriadis' algorithm, or EXP3 or INF works properly can be tackled similarly. In particular:

• Our bound in Theorem 7 (based on the constant in Lemma 4) can be adapted to the case of rational coefficients in M (instead of just 0 and 1);

• EXP3 and INF have no problem for stochastic cases as discussed below.

This leads to the following extensions:

• A first natural extension is the case of rational coefficients with a given denominator. The bound can be adapted to this case.