Maria Cristina Riff

Elizabeth Montero
email: elizabeth.montero@inf.utfsm.cl

Bertrand Neveu
email: bertrand.neveu@sophia.inria.fr

María Cristina Riff

C-Strategy: A Dynamic Adaptive Strategy for the CLONALG Algorithm

The control of parameters during the execution of bioinspired algorithms is an open research area. In this paper, we propose a new parameter control strategy for the immune algorithm CLONALG. Our approach is based on reinforcement learning ideas. We focus our attention on controlling the number of clones. Our approach provides an efficient and low cost adaptive technique for parameter control. We use instances of the Travelling Salesman Problem. The results obtained are very encouraging.

Introduction

When we design a bio-inspired algorithm to address a specific problem we need to define a representation, and a set of components to solve the problem. We also need to choose parameter values which we expect will give the algorithm the best results. This process of finding adequate parameter values is a timeconsuming task and considerable effort has already gone into automating this process [START_REF] Hinterding | Adaptation in evolutionary computation: A survey[END_REF]. Researchers have used various methods to find good values for the parameters, as these can affect the performance of the algorithm in a significant way. The best known mechanism to do this is tuning parameters on the basis of experimentation, something similar to a generate and test procedure. Taking into account that an immune algorithm run is an intrinsically dynamic adaptive process, we can expect that a dynamic adaptation of the parameters during the search could help to improve the performance of the algorithm, as it has been shown for other metaheuristics [START_REF] Davis | Adapting operator probabilities in genetic algorithms[END_REF], [START_REF] Deb | Understanding interactions among genetic algorithm parameters[END_REF], [START_REF] Smith | Operator and parameter adaptation in genetic algorithms[END_REF], [START_REF] Tuson | Adapting operator settings in genetic algorithms[END_REF], [START_REF] Eiben | Evolutionary algorithms with on-the-fly population size adjustment[END_REF], [START_REF] Montero | Self-calibrating strategies for evolutionary approaches that solve constrained combinatorial problems[END_REF], [START_REF] Montero | Improving MMAS using parameter control[END_REF]. The idea of adapting and self-adapting parameters during the run is not something new, but we need to manage a trade-off between the improvement of the search and the adaptation cost. In this case, the idea is to monitor the search to be able to trigger actions from an adaptive parameter control strategy, in order to improve the performance of the well-known immune algorithm CLONALG [START_REF] De Castro | Artificial Immune Systems: A New Computational Intelligence Approach[END_REF].

We present a brief revision of the related work in the parameter control domain in the following section. The algorithm CLONALG is briefly described in section 3 with remarks on the aspects of parameter control. We propose a method in section 4 which includes a monitoring task. We have tested the algorithm using our strategy with instances of the Travelling Salesman Problem, that is reported in section 5. The conclusions in section 6 resume our experience from this study and the following trends of this research.

Related Work

We can classify the parameter selection process into two different methods:tuning and parameter control, [START_REF] Eiben | Parameter control in evolutionary algorithms[END_REF]. Tuning, as mentioned above, implies in the worst case a generate and test procedure, in order to define which are the "best" parameter values for a bio-inspired algorithm. Usually, these parameter values are fixed for all the runs of the algorithm. As we mentioned before, Revac [START_REF] Nannen | Relevance estimation and value calibration of evolutionary algorithm parameters[END_REF] has been proposed to tune of evolutionary algorithms in an efficient way. The idea is to use an estimation of distribution algorithms in an iterative process to converge upon a set of good values for the parameters. The ParamILS [START_REF] Hutter | Automatic algorithm configuration based on local search[END_REF] technique is also available for tuning, it uses a local search strategy to search the parameter values space in an efficient way. The Racing methods [START_REF] Birattari | A racing algorithm for configuring metaheuristics[END_REF] use statistical information about the performance of the tuned algorithm to carry out a race between different parameter configurations.

We can expect that the "best" parameter values depend on the step of the algorithm we are using. The main disadvantage of the tuning methods is the time they require to perform the calibration process. For that reason, the idea of designing strategies to allow the algorithm to change its parameter values during its execution appears as a promising way to improve the performance of a bio-inspired algorithm. In this context many ideas have been proposed in the evolutionary research community, including self-adaptability methods. In selfadaptive methods, each individual incorporates information which can influence the parameter values or the evolution of these processes can be used to define some criteria to produce changes, [START_REF] Riff | Inheriting parents operators: a new dynamic strategy to improve evolutionary algorithms[END_REF], [START_REF] Gómez | Self adaptation of operador rates in evolutionary algorithms[END_REF], [START_REF] Montero | Self-calibrating strategies for evolutionary approaches that solve constrained combinatorial problems[END_REF].

The parameter control methods can be grouped according to the parameter controlled. For example, we can find strategies which control the population size. The parameter-less algorithm [START_REF] Lobo | The parameter-less genetic algorithm in practice[END_REF] implements a population control performing simultaneous runs of a genetic algorithm with different population sizes, emulating a race among the multiple populations. As the search progresses, parameterless eliminates small populations which present a worse average performance than larger populations. A recent research of Eiben et al. [START_REF] Eiben | Evolutionary algorithms with on-the-fly population size adjustment[END_REF] proposes a strategy to change the population size of the genetic algorithm, taking into account the state of the search. The method is called PRoFIGA. PRoFIGA is able to change the population size according to the convergence and stagnation of the search process.

It is also possible to find strategies which combine methods to control several parameters. In [START_REF] Srinivasa | A self-adaptive migration model genetic algorithm for data mining applications[END_REF], Srinivasa et al. propose the SAMGA method. The SAMGA method incorporates an adaptive parameter control technique to vary the population size, mutation rate and crossover rate of each population of a parallel genetic algorithm. The changes are performed according to the relative performance of one population against the others. Moreover, it dynamically controls the migration rate based on the stagnation of the global search procedure. Recently, Montiel et al. [START_REF] Montiel | Human evolutionary model: A new approach to optimization[END_REF] introduced the HEM algorithm which implements an evolutionary algorithm. The HEM algorithm incorporates an adaptive intelligent system which allows the parameter values of the individuals to be controlled in a self-adaptive way. The intelligent system is able to learn about the evolutionary algorithm performance based on expert knowledge and feedback from the search process.

The parameter setting problem, usually suited for evolutionary algorithms, also has been extended to other research areas. In [START_REF] Mezura-Montes | Parameter control in differential evolution for constrained optimization[END_REF] Mezura and Palomeque proposed a parameter control strategy for the constraint-handling mechanism used by a Differential Evolution Algorithm. In their strategy they self-adapted three parameters and one parameter is deterministically controlled. In [START_REF] Tan | Balancing exploration and exploitation with adaptive variation for evolutionary multi-objective optimization[END_REF], the authors proposed a two-fold strategy to deterministically control the exploration/exploitation balance during the search and to adaptively coordinate the crossover and mutation operation. This strategy is implemented in an evolutionary multi-objective algorithm with binary representation.

One of the first approaches to parameter control in immune systems was proposed by [START_REF] Garrett | Parameter-free, adaptive clonal selection[END_REF]. In this work Garret proposed a self-adaptive control strategy to control the amount of clones generated from each cell. The amount of clones is established according to predictions of the performance of the algorithm based on the performance that the algorithm showed in the previous iterations. In [START_REF] Hu | Adaptive clonal selection with elitism-guided crossover for function optimization[END_REF], Hu et al. propose an adaptive control strategy to control the exploration/exploitation rate in the search. In their approach, they divide the repertoire set in three subsets based on their quality. The cells in the best set are mutated according to a micro-mutation process based on the normalized performance of these cells. This approach implements an elitist crossover operator that operates in the two better set of cells. The strategy proposed is able to control the amount of exploitation/exploration according to the variation of two parameters associated to the micro-mutation and crossover operators here defined. The parameter associated to the micro-mutation operator is controlled in an adaptive way and the parameter associated to the crossover operator is controlled in a deterministic way.

The advantage of changing the parameter values during the execution becomes more relevant when we are tackling NP-complete problems. In this paper we introduce a strategy which is able to implement efficient parameter control by using the idea of taking into account the evolution, while also taking action during the execution in order to accelerate the convergence process.

Description of CLONALG

Artificial Immune Systems (AIS) have been defined as adaptive systems inspired by the immune system and applied to problem solving. In this paper we are interested in CLONALG that is an artificial immune algorithm based on Clonal Selection. CLONALG has successfully been applied to solve various complex problems. It follows the basic theory of an immune system response to pathogens. Roughly speaking, the components of the algorithm are cells named antibodies and an antigen that is an invader attacking the immune system. The immune system reproduces those cells that recognise antigens. Cells that match well are cloned (create copies of themselves). The better the match, the more clones. The clones undergo small mutations which may improve their recognition. A selection process retains the best cells as memory cells. The CLONALG pseudocode is shown in figure 1.

The algorithm starts generating a random population of Cells and computing their fitness related to the problem at hand. The iterative process begins by constructing a sub-population S p of size n, composed by the n best antibodies belonging the population of Cells. A new population of clones P c is constructed by generating C clones of each element on S p . This population of clones follows a mutation process, according to a mutation rate, computed using an exponential distribution with parameter ρ, in order to improve the evaluation value of the clones. The set Cells is updated including the n best antibodies from P c and incorporating new (An) randomly generated antibodies to construct the population of size A.

The standard parameters of CLONALG are: ρ, C, n and A. Our study is focused on dynamically controlling parameter C. The value of C is related to the exploitation of the search space, because at each iteration the algorithm tries to improve these antibodies applying a mutation procedure.

Parameter Control Strategy

The key idea in our approach is to design a low computational cost strategy to control the population size of clones in CLONALG. Our aim is to propose and to evaluate a strategy that allows parameter adaptation in an immune based approach according to the problem at hand. We propose an adaptive reinforcement control evaluated in the experimental comparison section.

Adaptive Control

CLONALG works with a set of antibodies as we mentioned before: A population of C clones. The value of this parameter can be controlled in order to guide the trade-off between the intensification and the diversification process of the algorithm. Our idea for this control comes from reinforcement learning. We clearly distinguish between a positive behaviour and a negative one either to reward (increase) or to penalize (reduce) the number of antibodies that belong to the set of clones.

Before we introduce our strategy we require the following definition:

Definition 1. Given a fitness function F to be maximized, a population P c of C clones and a population S p of n selected antibodies. We define a success measure for the parameter C in the ith iteration, P S i (C) as:

P S i (C) = F (BP c ) -F (BS p ) ( 1 )
where F (BP c ) and F (BS p ) are the respective fitness of the best antibody of the populations P c and S p . The idea is to evaluate if there is an improvement of the best pre-solution found by the algorithm. In this case we reward this situation by increasing the number of C clones, that means more exploitation or intensification of the search f the algorithm. When the algorithm decides to reduce the number of clones it produces that the new population of cells will be completed adding more new cells. From a conceptual point of view it can be interpreted as more exploration or diversification of the search of the algorithm.

Figure 2 shows the structure of the algorithm with the adaptive control. It is important to remark that in lines 3 and 7 the fitness are data available in the original CLONALG code. We explicitly show their evaluations here because we use them in lines 9 and 10. The algorithm computes the P S i (C) value in order to increase (decrease) C (number of clones) value. Obviously, the algorithm manages extreme situations, as C must be greater than 1.

Remark 1. The previous definition and the procedure can easily be adapted when the function to be optimized must be minimized, as we use in the experimental comparison section with the travelling salesman problem

The main motivation of the work reported in this paper is to address the on-line parameter determination problem. We propose a new method which improves the performance of CLONALG and also allows a self-adaptation of the parameter C without adding a significant overhead and without introducing any major changes to its original algorithmic design. The decisions made during the search are based on information available from a monitoring process. Such information allows the algorithm to trigger changes when deemed necessary, based on the potential advantages that such changes could bring. It is also important to keep a proper balance between the level of improvement achieved and the computational cost required to reach it.

CLONALG-with adaptive control

P Si(C) = F (BPc) -F (BSp) ( 9 ) If (P Si(C) > 0) then C = C + 1 else If C > 2 then C = C -1 ( 10
) i = i + 1 (11

Experimental Comparison

In this section, we experimentally evaluate the CLONALG algorithm using our strategy.

Experimental Set-Up

The hardware platform for the experiments was a PC CPU Intel Core i7-920 4GB RAM, Mandriva 2009.

Travelling Salesman Problem

Travelling Salesman Problem (TSP) is a classical optimization problem, defined as the task of finding the shortest path for visiting N cities and returning to the original point. TSP is a classical combinatorial optimization problem [START_REF] Stützle | A comparison of nature inspired heuristics on the traveling salesman problem[END_REF], [START_REF] Zhang | A novel local search algorithm for the traveling salesman problem that exploits backbones[END_REF], [START_REF] Richter | Improving the efficiency of helsgauns lin-kernighan heuristic for the symmetric tsp[END_REF]. Several approaches have been proposed to face the TSP using artificial immune systems. The more well-known approach was proposed by de Castro and von Zuben [START_REF] De Castro | The clonal selection algorithm with engineering applications[END_REF] that we briefly present in the following section.

TSP-CLONALG:

The original version of CLONALG to solve TSP [START_REF] De Castro | The clonal selection algorithm with engineering applications[END_REF] called TSP-CLONALG, has some particularities that differs of its standard version to be noticed. The number of clones of each antibody is calculated according to equation 2. We can observe that the amount of clones is determined by the parameter β and the ranking of the antibody k. The total amount of clones at each iteration is calculated according to equation 3.

C k = round β • A k (2) C = n k=1 C k (3)
In order to improve the performance of the standard version of the CLONALG to solve TSP the authors have added two new parameters: The number of randomly generated cells to be included (d) after a fixed number of iterations (iter). In the original version of CLONALG the random cells were incorporated at each iteration and the number of random cells was fixed at (An). Thus, this version of the algorithm uses six parameters: ρ, β, A, n, d, iter. Using β and A it computes the total amount of parameters C.

In this paper we compare the performance between TSP-CLONALG using 6 parameters and the standard CLONALG version using our strategy. In our study we have detected that the parameters ρ and β are not relevant for the algorithm when it uses our control strategy. Thus, for our approach we just considered 2 parameters (A, n).

In both versions of CLONALG the representation is a permutation of the cities and the affinity measure is the length of the tour associated to each instantiation or antibody. We have considered as stopping condition when the algorithm is not able to get a better solution during 40 iterations.

Before continuing our description of tests we will describe the tuning method used in our experiments.

Tuning Method: Revac

In our experiments we used the tuning strategy called Relevance Estimation and Value Calibration method, Revac. The Revac method was proposed by Eiben & Nannen in [START_REF] Nannen | Relevance estimation and value calibration of evolutionary algorithm parameters[END_REF]. Revac is defined as an estimation of distribution algorithm [START_REF] Pelikan | A Survey of Optimization by Building and Using Probabilistic Models[END_REF]. It works with a set of parameter configurations, a population. For each parameter, Revac starts the search process with a uniform distribution of values within a given range. As the process advances Revac performs transformation operations with the aim of reducing each parameter distribution to a range of values that performs the best.

Revac method works with a population of 100 parameter configurations and at each iteration generates a new calibration using a crossover and a mutation operator. Revac performs a total amount of 1000 runs of the algorithm to be tuned, that means in average 10 9 number of evaluations.

Evaluation of the Adaptive Strategy

We report the results using the strategy considering Clones Variations (C).

Scenario 1.

In our experiments we have analyzed the results achieved by de Castro and von Zuben in [START_REF] De Castro | The clonal selection algorithm with engineering applications[END_REF]. In their work they study the results obtained by CLONALG method in an instance of 30 cities from [START_REF] Moscato | Stochastic versus deterministic update in simulated annealing[END_REF] displayed in figure 3. For this scenario we have used the information described in [START_REF] De Castro | The clonal selection algorithm with engineering applications[END_REF] and we have interpreted their swap operator as the well-known 2-opt as the move for mutation. We have used two versions of tuned CLONALG. The first configuration uses the parameter values proposed by the authors of the method, we called this configuration hand-tuned configuration. The second one has automatically been obtained by Revac. Table 1 shows the original parameter configuration and this obtained by Revac. In the original version for TSP, the algorithm differs from the classical CLONALG. For TSP two new parameters are included as we mentioned in section 5.2. It considers a set of d = 60 random cells to be included after every 20 iterations. Our adaptive technique works on the standard CLONALG, thus these two parameters are not required for the algorithm using control. Table 2 shows the performance obtained by the tuned parameter configurations and the adaptive strategy proposed for solving the 30 city instance. The three versions of the algorithm found the optimal value of 48873.

Fig. 3. Coordinates of the 30 cities instance

From the results we observe that the C-strategy outperforms all other versions. However it appears using more evaluations than the tuned versions. The algorithm using C-strategy has found the optimal tour 25 times, of the 50 total runs, the original hand-tuned version found 6 times and the version tuned using Revac only found the optimal solution 8 times.

We do not know how many times the algorithm has been run and the corresponding number of evaluations done to obtain the hand-tuned parameters configuration. Revac performed around of 10 8 number of evaluations in order to obtain the parameter values. These number of evaluations are not included in the Table 2. Scenario 2. We have tested our technique using symmetric travelling salesman problem instances from TSPLib 1 . The TSPLib is a set of more than 100 instances of the TSP collected by Professor Gerhard Reinelt in the mid 1990â ȂŹ s and that has been extensively used in the metaheuristic research area in the last years [START_REF] Stützle | Max-min ant system and local search for the traveling salesman problem[END_REF], [START_REF] Sun | An immune optimization algorithm for tsp problem[END_REF], [START_REF] Yang | Solving traveling salesman problems using generalized chromosome genetic algorithm[END_REF].

We have also calibrated the algorithm using Revac to solve other instances of the TSP. For this, we have also estimated the value of parameter d, which corresponds to the number of random cells to be inserted each 20 iterations. The results obtained are shown in Table 3. We can observe that, for each instance, Revac has identified different parameter values. This means that the behaviour of the algorithm strongly depends on both the problem instance and the values of its parameters. The amount of selected antibodies, n, ranges from around 20% to 95% of the population size and the amount of random antibodies, d, incorporated to the population each 20 generations ranges from 3% to 50% of the population size. The size of the population of cells A and ρ are also dependent on the instance at hand. The algorithm using our adaptive strategy does not require to be re-tuned each time it must solve another new instance of the problem. Thus, the time invested on tuning is strongly reduced.

Graphs in figure 4 show the relative distance to optimun for each instance comparing the performance of the C-strategy against to the hand-tuned and Revac-tuned performance. Our strategy helps the algorithm to improve the quality of the solutions found by the tuned versions. The greater is the number of cities the more noteworthy are the results obtained by the C-strategy.

Conclusions

We propose in this paper a strategy for an adaptive parameter control for CLON-ALG.

For intensification, our strategy works with the number of clones to be improved. When the improvement procedure has been applied successfully the algorithm increases the number of clones which will follow the improvement process in the next iteration.

In our design we have taken into account the overhead produced by including our strategy. It is a low cost strategy. The main advantage of using adaptive strategies to improve the behavior of the CLONALG comes from non-existing pre-execution time. The tuning techniques can obtain very good results to find tuned parameters for a given problem, but the parameter configuration step requires a huge amount of time compared to the time needed for the tuned algorithm to find a solution.

The original TSP-CLONALG algorithm uses 6 parameters. Accorditong to the results obtained using Revac, the value of these parameters are different for different instances of the TSP. Thanks to our dynamic control strategy the number of parameters to solve TSP is reduced to only 2. Moreover, these 2 parameters require only one tuning step whatever is the TSP instance to be solved from these used in our tests.

The algorithm using our adaptive strategy does not require to be re-tuned for solving new instances of the problem. Thus, the time invested for tuning is strongly reduced. The strategy that is proposed helps CLONALG to improve its performance when it is solving Travelling Salesman Problem instances. The Cstrategy is able to improve the performance in all instances tested with a relative average improvement of 3.3.

As future work we would like to evaluate our strategy in CLONALG using other hard combinatorial problems. We are also considering to include this strategy in other kinds of artificial immune algorithms based on Immune Network Models. We are working to include this strategy in CD-NAIS, an immune algorithm that solves hard instances of Constraint Satisfaction Problems. A promising research area is the collaboration between various parameter control strategies.

Fig. 1 .

 1 Fig. 1. Standard version of CLONALG

  Begin Cells = Initialise population of A antibodies (1) i = 1; Repeat Sp = Selected n best antibodies from Cells (2) F (BSp)= Fitness of the best antibody in Sp (3) Pc = Generate C clones of each antibody from Sp (4) Pc = Mutation(Pc) ( 5 ) Pc = Selected n best antibodies from Pc (6) F (BPc) = Fitness of the best antibody in Pc (7) New-Cells = Pc + Generate (An) New antibodies (8)

Fig. 2 .

 2 Fig. 2. Clonal Selection Algorithm with Adaptive Control

Fig. 4 .

 4 Fig. 4. Performance: Hand-Tuned, Revac-Tuned and C-strategy

Table 1 .

 1 Parameter configurations for the 30 cities instance

	Instance	ρ β A n d
	Hand-tuned	2.5 2.0 300 150 60
	Revac	1.2 2.8 190 178 97

Table 2 .

 2 Performance evaluation in the 30 cities instance. AS: Average Solution, %ΔAv: Relative distance to optimum of the average of 10 runs, EAv: Average evaluations, E best : Evaluations of the best run, OF: Times the optimum was found. Tuned by Revac for 30 cities instance. Revac required around of 10 8 number of evaluations for tuning this instance. These evaluations are not included in the number of evaluations reported in this table.

	Algorithm	AS %ΔAv EAv	E best OF
	Hand-tuned	50795	3.9 48693	51256 6
	Tuned by Revac *	50976	4.3 45649	47203 8
	C-strategy	50109	2.5 76413	78990 25
	*

Table 3 .

 3 Revac tuning results

	Instance	ρ	β	A	n	d
	burma14	1.8	0.1	87	87	36
	ulysses22	1.4	0.2	84	84	65
	eil51	1.0	0.1	120 120 4
	berlin52	1.8	0.3	112 70	37
	st70	1.1	0.2	109 20	36
	eil76	1.7	0.2	142 92	53
	pr76	0.8	0.1	93	63	28
	gr96	2	0.4	133 94	45
	kroA100	1.7	0.7	104 95	45
	eil101	2.5	0.7	147 103 52
	ch130	1.4	0.8	149 73	52
	ch150	2.4	3.2	219 164 59
	gr202	3.1	0.5	152 150 50
	tsp225	2.4	2.0	215 136 43
	a280	0.9	3.8	222 203 70
	pcb442	0.7	4.1	206 177 59
	gr666	1.1	2.6	194 174 97
	pr1002	0.5	1.4	172 110 66

Table 4 .

 4 Performance of Hand-tuned configuration. AS: Average Solution, %ΔAv: Relative distance to optimum of the average of 10 runs, EAv: Average evaluations, BS: Best solution of 10 runs, %ΔBest: Relative distance to optimum of the best solution found, E best : Evaluations of the best run.

	Instance	AS %ΔAv EAv	BS %ΔBest E best
	burma14	3333	0.3	16747	3323	0.0	9824
	ulysses22	7073	0.8	36631	7013	0.0	30925
	eil51	450	5.5 126175	435	2.1 113980
	berlin52	8093	7.3 135670	7787	3.2 135263
	st70	720	6.6 210339	693	2.7 213953
	eil76	589	9.4 207454	559	3.9 239400
	pr76	115678	7.0 244200	112825	4.3 291478
	gr96	59110	7.1 343631	57276	3.7 346040
	kroA100	22874	7.5 374095	21596	1.5 407959
	eil101	685	8.9 315624	672	6.8 324708
	ch130	6661	9.0 510299	6448	5.5 519051
	ch150	7248 11.0 577297	6906	5.8 603625
	gr202	43353	8.0 884525	42296	5.3 820052
	tsp225	4323 10.4 984360	4234	8.1 1069647
	a280	2965 15.0 1338809	2923	13.3 1356322
	pcb442	58021 14.3 2328181	56343	11.0 2385774
	gr666	338838 15.1 4014526	332335	12.9 4210456
	pr1002	308884 19.2 6568289	299594	15.7 6749630

Table 5 .

 5 Performance of Revac-tuned configurations. AS: Average Solution, %ΔAv: Relative distance to optimum of the average of 10 runs, EAv: Average evaluations, BS: Best solution of 10 runs, %ΔBest: Relative distance to optimum of the best solution found, E best : Evaluations of the best run, Δ + : Improvement of relative distance to optimum compared to the hand-tuned version. The amount of evaluations specified does not include the around of 10 9 evaluations performed by the tuning process.

	Instance	AS %ΔAv EAv †	BS %ΔBest E best †	Δ +
	burma14	3388	2.0	518	3336	0.4	707	-1.7
	ulysses22	7150	2.0	1561	7013	0.0	1823	-1.1
	eil51	473 11.1	4052	451	5.9	4256	-5.5
	berlin52	8305 10.1	9799	7547	0.1	9565	-2.8
	st70	750 11.1	13008	709	5.0	17287	-4.5
	eil76	594 10.4	12049	573	6.5	13426	-1.0
	pr76	119710 10.7	7854	117792	8.9	6036	-3.7
	gr96	61503 11.4	33120	58362	5.7	39705	-4.3
	kroA100	23285	9.4	56371	21924	3.0	50938	-1.9
	eil101	683	8.6 122179	663	5.4 144320	0.2
	ch130	6708	9.8 124983	6455	5.6 131313	-0.8
	ch150	7124	9.1 724714	6964	6.7 749362	1.9
	gr202	43833	9.1 274549	43150	7.4 306989	-1.2
	tsp225	4343 10.9 849222	4264	8.9 819274	-0.5
	a280	2905 12.6 2227123	2819	9.3 2303958	2.3
	pcb442	57221 12.7 3372096	56269	10.8 3344210	1.6
	gr666	333980 13.5 3779875	327663	11.3 3979137	1.7
	pr1002	297017 14.7 4857562	290840	12.3 4965511	4.6
	†

Table 6 .

 6 Performance of C-strategy. AS: Average Solution, %ΔAv: Relative distance to optimum of the average of 10 runs, EAv: Average evaluations, BS: Best solution of 10 runs, %ΔBest: Relative distance to optimum of the best solution found, E best : Evaluations of the best run, Δ + : Improvement of relative distance to optimum compared to the hand-tuned version.

	Instance	AS %ΔAv	EAv	BS %ΔBest	E best	Δ +
	burma14	3326	0.1	13326	3323	0.0	13069	0.2
	ulysses22	7023	0.1	37915	7013	0.0	39049	0.7
	eil51	444	4.2	280361	429	0.7	212817	1.3
	berlin52	7870	4.4	325300	7542	0.0	276394	2.9
	st70	705	4.5	624727	683	1.2	673642	2.1
	eil76	567	5.3	715914	551	2.4	649833	4.1
	pr76	113273	4.7	742192	109504	1.2	838858	2.2
	gr96	57878	4.8	1484345	56145	1.7	1653666	2.2
	kroA100	22545	5.9	1643365	21753	2.2	1629237	1.5
	eil101	667	6.0	1394766	652	3.7	1300641	2.9
	ch130	6444	5.5	2803527	6263	2.5	3069959	3.5
	ch150	6931	6.2	3969191	6722	3.0	4137293	4.9
	gr202	42771	6.5	7422549	41886	4.3	41886	1.4
	tsp225	4186	6.9	9877811	4076	4.1	9599825	3.5
	a280	2809	8.9 17376294	2734	6.0 17240650	6.1
	pcb442	55693	9.7 45624518	54588	7.5 49324914	4.6
	gr666	323210	9.8 117452296	317200	7.8 120874559	5.3
	pr1002	282506	9.1 286994425	276963	6.9 280607399	10.2

http://www.iwr.uni-heidelberg.de/groups/comopt/software/TSPLIB95/

Project 1080110.


Jean-Yves Audibert
email: audibert@imagine.enpc.fr

Sébastien Bubeck
email: sebastien.bubeck@inria.fr

Rémi Munos
email: remi.munos@inria.fr

Best Arm Identification in Multi-Armed Bandits

We consider the problem of finding the best arm in a stochastic multi-armed bandit game. The regret of a forecaster is here defined by the gap between the mean reward of the optimal arm and the mean reward of the ultimately chosen arm. We propose a highly exploring UCB policy and a new algorithm based on successive rejects. We show that these algorithms are essentially optimal since their regret decreases exponentially at a rate which is, up to a logarithmic factor, the best possible. However, while the UCB policy needs the tuning of a parameter depending on the unobservable hardness of the task, the successive rejects policy benefits from being parameter-free, and also independent of the scaling of the rewards. As a by-product of our analysis, we show that identifying the best arm (when it is unique) requires a number of samples of order (up to a log(K) factor) ∑ i 1/∆ 2 i , where the sum is on the suboptimal arms and ∆ i represents the difference between the mean reward of the best arm and the one of arm i. This generalizes the well-known fact that one needs of order of 1/∆ 2 samples to differentiate the means of two distributions with gap ∆.

Introduction

In the multi-armed bandit problem, at each stage, an agent (or forecaster) chooses one action (or arm), and receives a reward from it. In its stochastic version, the reward is drawn from a fixed probability distribution given the arm. The usual goal is to maximize the cumulative sum of rewards, see [START_REF] Robbins | Some aspects of the sequential design of experiments[END_REF]; [START_REF] Auer | Finite-time analysis of the multiarmed bandit problem[END_REF] among many others. Since the forecaster does not know the distributions, he needs to explore (try) the different actions and yet, exploit (concentrate its draws on) the seemingly most rewarding arms. In this paper, we adopt a different viewpoint. We assume that after a given number of pulls, the forecaster is asked to output a recommended arm. He is then only evaluated by the average payoff of his recommended arm. This is the so-called pure exploration problem, [START_REF] Bubeck | Pure exploration in multi-armed bandits problems[END_REF].

The distinguishing feature from the classical multi-armed bandit problem described above is that the exploration phase and the evaluation phase are separated. Thus, there is no explicit trade-off between the exploration and the exploitation while pulling the arms. The target of Hoeffding and Bernstein races, see [START_REF] Maron | Hoeffding races: Accelerating model selection search for classification and function approximation[END_REF]; [START_REF] Mnih | Empirical bernstein stopping[END_REF] among others, is more similar to ours. However, instead of trying to extract from a fixed number of rounds the best action, racing algorithms try to identify the best action at a given confidence level while consuming the minimal number of pulls. They optimize the budget for a given confidence level, instead of optimizing the quality of the recommendation for a given budget size. Another variant of the best arm identification task is the problem of minimal sampling times required to identify an ε-optimal arm, see in particular [START_REF] Domingo | Adaptive sampling methods for scaling up knowledge discovery algorithms[END_REF] and [START_REF] Even-Dar | Action elimination and stopping conditions for the multi-armed bandit and reinforcement learning problems[END_REF].

We now illustrate why this is a natural framework for numerous applications. Historically, the first occurrence of multi-armed bandit problems was given by medical trials, see [START_REF] Robbins | Some aspects of the sequential design of experiments[END_REF]. In the case of a severe disease, ill patients only are included in the trial and the cost of picking the wrong treatment is high. It is important to minimize the cumulative regret, since the test and cure phases coincide. However, for cosmetic products, there exists a test phase separated from the commercialization phase, and one aims at minimizing the regret of the commercialized product rather than the cumulative regret in the test phase, which is irrelevant.

Another motivating example concerns channel allocation for mobile phone communications. During a very short time before the communication starts, a cellphone can explore the set of channels to find the best one to operate. Each evaluation of a channel is noisy and there is a limited number of evaluations Parameters available to the forecaster: the number of rounds n and the number of arms K.

Parameters unknown to the forecaster: the reward distributions ν1, . . . , νK of the arms.

For each round t = 1, 2, . . . , n;

(1) the forecaster chooses It ∈ {1, . . . , K},

(2) the environment draws the reward X It,T I t (t) from νI t and independently of the past given It.

At the end of the n rounds, the forecaster outputs a recommendation Jn ∈ {1, . . . , K}.

Figure 1: The pure exploration problem for multi-armed bandits.

before the communication starts. The connection is then launched on the channel which is believed to be the best. Opportunistic communication systems rely on the same idea. Again the cumulative regret during the exploration phase is irrelevant since the user is only interested in the quality of its communication starting after the exploration phase. More generally, the pure exploration problem addresses the design of strategies making the best possible use of available resources in order to optimize the performance of some decision-making task. That is, it occurs in situations with a preliminary exploration phase in which costs are not measured in terms of rewards but rather in terms of resources that come in limited budget (the number of patients in the test phase in the clinical trial setting and the time to connect in the communication example).

Problem setup

A stochastic multi-armed bandit game is parameterized by the number of arms K, the number of rounds (or budget) n, and K probability distributions ν 1 , . . . , ν K associated respectively with arm 1, . . . , arm K. These distributions are unknown to the forecaster. For t = 1, . . . , n, at round t, the forecaster chooses an arm I t in the set of arms {1, . . . , K}, and observes a reward drawn from ν It independently from the past (actions and observations). At the end of the n rounds, the forecaster selects an arm, denoted J n , and is evaluated in terms of the difference between the mean reward of the optimal arm and the mean reward of J n . Precisely, let µ 1 , . . . , µ K be the respective means of ν 1 , . . . , ν K . Let µ * = max k∈{1,...,K} µ k . The regret of the forecaster is

r n = µ * -µ Jn .
For sake of simplicity, we will assume that the rewards are in [0, 1] and that there is a unique optimal arm. Let i * denote this arm (so, µ i * = µ * ). For i ̸ = i * , we introduce the following suboptimality measure of arm i:

∆ i = µ * -µ i .
For reasons that will be obvious later, we also define ∆ i * as the minimal gap

∆ i * = min i̸ =i * ∆ i .
We introduce the notation (i) ∈ {1, . . . , K} to denote the i-th best arm (with ties break arbitrarily), hence

∆ i * = ∆ (1) = ∆ (2) ≤ ∆ (3) ≤ . . . ≤ ∆ (K) .
Let e n denote the probability of error, that is the probability that the recommendation is not the optimal one:

e n = P(J n ̸ = i * ).
We have Er n = ∑ i̸ =i * P(J n = i)∆ i , and consequently

∆ i * e n ≤ Er n ≤ e n .
As a consequence of this equation, up to a second order term, e n and Er n behave similarly, and it does not harm to focus on the probability e n .

For each arm i and all time rounds t ≥ 1, we denote by T i (t) the number of times arm i was pulled from rounds 1 to t, and by X i,1 , X i,2 , . . . , X i,Ti(t) the sequence of associated rewards. Introduce X i,s = 1 s ∑ s t=1 X i,t the empirical mean of arm i after s pulls. In the following, the symbol c will denote a positive numerical constant which may differ from line to line.

The goal of this work is to propose allocation strategies with small regret, and possibly as small as the best allocation strategy which would know beforehand the distributions ν 1 , . . . , ν K up to a permutation. Before going further, note that the goal is unachievable for all distributions ν 1 , . . . , ν K : a policy cannot perform as well as the "oracle" allocation strategy in every particular cases. For instance, when the supports of ν 1 , . . . , ν K are disjoint, the oracle forecaster almost surely identifies an arm by a single draw of it. As a consequence, it has almost surely zero regret for any n ≥ K. The generic policy which does not have any knowledge on the K distributions cannot reproduce this performance for any K-tuple of disjointly supported distributions. In this work, the above goal of deciding as well as an oracle will be reached for the set of Bernoulli distributions with parameters in (0, 1), but the algorithms are defined for any distributions supported in [0, 1].

We would like to mention that the case K = 2 is unique and simple since, as we will indirectly see, it is optimally solved by the uniform allocation strategy consisting in drawing each arm n/2 times (up to rounding problem), and at the end recommending the arm with the highest empirical mean. Therefore, our main contributions concern more the problem of the budget allocation when K ≥ 3. The hardness of the task will be characterized by the following quantities

H 1 = K ∑ i=1 1 ∆ 2 i and H 2 = max i∈{1,...,K} i∆ -2 (i) .
These quantities are equivalent up to a logarithmic factor since we have (see Section 6.1)

H 2 ≤ H 1 ≤ log(2K)H 2 .
(1) Intuitively, we will show that these quantities are indeed characteristic of the hardness of the problem, in the sense that they give the order of magnitude of the number of samples required to find the best arm with a reasonable probability. This statement will be made precise in the rest of the paper, in particular through Theorem 2 and Theorem 4.

Outline. In Section 3, we propose a highly exploring policy based on upper confidence bounds, called UCB-E (Upper Confidence Bound Exploration), in the spirit of UCB1 Auer et al. (2002). We prove that this algorithm, provided that it is appropriately tuned, has an upper bound on the probability of error e n of order exp ( -c n H1

)

. The core problem of this policy is the tuning of the parameter. The optimal value of the parameter depends on H 1 , which has no reason to be known beforehand by the forecaster, and which, to our knowledge, cannot be estimated from past observations with sufficiently high confidence in order that the resulting algorithm still satisfies a similar bound on e n .

To get round this limitation, in Section 4, we propose a simple new policy called SR (Successive Rejects) that progressively rejects the arms which seem to be suboptimal. This algorithm is parameter-free and its probability of error e n is at most of order exp ( -

n log(2K)H2
)

. Since H 2 ≤ H 1 ≤ log(2K)H 2 , up to at most a logarithmic term in K, the algorithm performs as well as UCB-E while not requiring the knowledge of H 1 .

In Section 5, we prove that H 1 and H 2 truly represent the hardness of the problem (up to a logarithmic factor). Precisely, we consider a forecaster which knows the reward distributions of the arms up to a permutation. When these distributions are of Bernoulli type with parameter in [p, 1 -p] for some p > 0, there exists a permutation of the distributions for which the probability of error of the (oracle) forecaster is lower bounded by exp ( -

cn p(1-p)H2
)

. Section 6 gathers some of the proofs. Section 7 provides some experiments testing the efficiency of the proposed policies and enlightening our theoretical results. We also discuss a modification of UCB-E where we perform a non-trivial online estimation of H 1 . We conclude in Section 8.

Example. To put in perspective the results we just mentioned, let us consider a specific example with Bernoulli distributions. Let

ν 1 = Ber ( 1 2 )
, and

ν i = Ber ( 1 2 -1 K i ) for i ∈ {2, . . . , K}.
Here, one can easily check that H 2 = 2K 2K . Thus, in this case, the probability of missing the best arm of SR is at most of order exp ( -

n 2 log(2K)K 2K )
. Moreover, in Section 5, we prove that there does not exist any forecaster (even with the knowledge of the distributions up to a permutation) with a probability of missing the best arm smaller than exp

( -11n K 2K
) for infinitely many n. Thus, our analysis finds that, for this particular reward distributions, the number of samples required to find the best arm is at least (of order of) K 2K , and SR actually finds it with (of order of) log(K)K 2K samples.

Highly exploring policy based on upper confidence bounds

In this section, we propose and study the algorithm UCB-E described in Figure 2. When a is taken of order log n, the algorithm essentially corresponds to the UCB1 policy introduced in [START_REF] Auer | Finite-time analysis of the multiarmed bandit problem[END_REF] Parameter: exploration parameter a > 0.

For i ∈ {1, . . . , K}, let Bi,s = Xi,s + √ a s for s ≥ 1 and Bi,0 = +∞. For each round t = 1, 2, . . . , n:

Draw It ∈ argmax i∈{1,...,K} B i,T i (t-1) . Let Jn ∈ argmax i∈{1,...,K} X i,T i (n) .
Figure 2: UCB-E (Upper Confidence Bound Exploration) algorithm.

its cumulative regret is of order log n. [START_REF] Bubeck | Pure exploration in multi-armed bandits problems[END_REF] have shown that algorithms having at most logarithmic cumulative regret, have at least a (non-cumulative) regret of order n -γ for some γ > 0. So taking a of order log n is inappropriate to reach exponentially small probability of error. For our regret notion, one has to explore much more and typically use a parameter which is essentially linear in n. Precisely, we have the following result, the proof of which can be found in Section 6.2.

Theorem 1 If UCB-E is run with parameter 0 < a ≤ 25 36 n-K H1 , then it satisfies e n ≤ 2nK exp ( - 2a 25
) .

In particular for a = 25 36 n-K H1 , we have e n ≤ 2nK exp

( -n-K 18H1
) .

The theorem shows that the probability of error of UCB-E is at most of order exp(-ca) for a ≥ log n. In fact, one can easily show a corresponding lower bound. In view of this, as long as a ≤ 25 36 n-K H1 , we can essentially say: the more we explore (i.e., the larger a is), the smaller the regret is. Besides, the smallest upper bound on the probability of error is obtained for a of order n/H 1 , and is therefore exponentially decreasing with n. The constant H 1 depends not only on how close the mean rewards of the two best arms are, but also on the number of arms and how close their mean reward is to the optimal mean reward. This constant should be seen as the order of the minimal number n for which the recommended arm is the optimal one with high probability. In Section 5, we will show that H 1 is indeed a good measure of the hardness of the task by showing that no forecaster satisfies e n ≤ exp ( -cn

H2

) for any distributions ν 1 , . . . , ν K , where we recall that

H 2 satisfies H 2 ≤ H 1 ≤ log(2K)H 2 .
One interesting message to take from the proof of Theorem 1 is that, with probability at least 1 -2nK exp ( -2a

25

)

, the number of draws of any suboptimal arm i is of order a∆ -2 i . This means that the optimal arm will be played at least n -caH 1 , showing that for too small a, UCB-E "exploits" too much in view of our regret target. Theorem 1 does not specify how the algorithm performs when a is larger than 25 36 n-K H1 . Nevertheless, similar arguments than the ones in the proof show that for large a, with high probability, only low rewarding arms are played of order a∆ -2 i times, whereas the best ones are all drawn the same number of times up to a constant factor. The number of these similarly drawn arms grows with a. In the limit, when a goes to infinity, UCB-E is exactly the uniform allocation strategy studied in [START_REF] Bubeck | Pure exploration in multi-armed bandits problems[END_REF]. In general1 , the uniform allocation has a probability of error which can be lower and upper bounded by a quantity of the form exp

( -c n∆ 2 i * K )
. It consequently performs much worse than UCB-E for a = 25 36

n-K H1 , since H 1 ≤ K∆ -2
i * , and potentially H 1 ≪ K∆ -2 i * for very large number of arms with heterogeneous mean rewards.

One straightforward idea to cope with the absence of an oracle telling us the value of H 1 would be to estimate online the parameter H 1 and use this estimation in the algorithm. Unfortunately, we were not able to prove, and do not believe that, this modified algorithm generally attains the expected rate of convergence. Indeed, overestimating H 1 leads to low exploring, and in the event when the optimal arm has given abnormally low rewards, the arm stops being drawn by the policy, its estimated mean reward is thus not corrected, and the arm is finally not recommended by the policy. On the contrary, underestimating H 1 leads to draw too much the suboptimal arms, precluding a sufficiently accurate estimation of the mean rewards of the best arms. For this last case, things are in fact much more subtle than what can be retranscribed in these few lines,

Let A1 = {1, . . . , K}, log(K) = 1 2 + ∑ K i=2 1 i , n0 = 0 and for k ∈ {1, . . . , K -1}, n k = ⌈ 1 log(K) n -K K + 1 -k ⌉ .
For each phase k = 1, 2, . . . , K -1:

(1) For each i ∈ A k , select arm i for n k -n k-1 rounds.

(2) Let A k+1 = A k \ arg mini∈A k Xi,n k (we only remove one element from A k , if there is a tie, select randomly the arm to dismiss among the worst arms).

Let Jn be the unique element of AK . and we notice that keeping track of a lower bound on H 1 would lead to the correct rate only under appropriate assumptions on the decrease of the sequence ∆ (k) , k ∈ {1, . . . , K}. In Section 7 we push this idea and propose a way to estimate online H 1 , however we solely justify the corresponding algorithm by experiments.

In the next section we propose an algorithm which does not suffer from these limitations.

Successive Rejects algorithm

In this section, we describe and analyze a new algorithm, SR (Successive Rejects), see Figure 3 for its precise description. Informally it proceeds as follows. First the algorithm divides the time (i.e., the n rounds) in K -1 phases. At the end of each phase, the algorithm dismisses the arm with the lowest empirical mean.

During the next phase, it pulls equally often each arm which has not been dismissed yet. The recommended arm J n is the last surviving arm. The length of the phases are carefully chosen to obtain an optimal (up to a logarithmic factor) convergence rate. More precisely, one arm is pulled

n 1 = ⌈ 1 log(K) n-K K ⌉ times, one n 2 = ⌈ 1 log(K) n-K K-1
⌉ times, ..., and two arms are pulled

n K-1 = ⌈ 1 log(K) n-K 2 ⌉ times. SR does not exceed the budget of n pulls, since, from the definition log(K) = 1 2 + ∑ K i=2 1 i , we have n 1 + . . . + n K-1 + n K-1 ≤ K + n -K log(K) ( 1 2 + K-1 ∑ k=1 1 K + 1 -k ) = n.
For K = 2, up to rounding effects, SR is just the uniform allocation strategy.

Theorem 2 The probability of error of SR satisfies

e n ≤ K(K -1) 2 exp ( - n -K log(K)H 2
) .

Proof: We can assume that the sequence of rewards for each arm is drawn before the beginning of the game. Thus the empirical reward for arm i after s pulls is well defined even if arm i has not been actually pulled s times. During phase k, at least one of the k worst arms is surviving. So, if the optimal arm i * is dismissed at the end of phase k, it means that X i * ,n k ≤ max i∈{(K),(K-1),...,(K+1-k)} X i,n k . By a union bound and Hoeffding's inequality, the probability of error e n = P(A K ̸ = {i * }) thus satisfies

e n ≤ K-1 ∑ k=1 K ∑ i=K+1-k P( X i * ,n k ≤ X (i),n k ) ≤ K-1 ∑ k=1 K ∑ i=K+1-k P( X (i),n k -µ (i) + µ * -X i * ,n k ≥ ∆ (i) ) ≤ K-1 ∑ k=1 K ∑ i=K+1-k exp ( -n k ∆ 2 (i) ) ≤ K-1 ∑ k=1 k exp ( -n k ∆ 2 (K+1-k)
) .

We conclude the proof by noting that by definition of n k and H 2 , we have

n k ∆ 2 (K+1-k) ≥ n -K log(K) 1 (K + 1 -k)∆ -2 (K+1-k) ≥ n -K log(K)H 2 . ( 2
)
The following theorem provides a deeper understanding of how SR works. It lower bounds the sampling times of the arms and shows that at the end of phase k, we have a high-confidence estimation of ∆ (K+1-k) up to numerical constant factor. This intuition will prove to be useful in Section 7, see in particular Figure 4.

Theorem 3 With probability at least 1 -K 3 2 exp ( -n-K 4log(K)H2
)

, for any arm j, we have

T j (n) ≥ n -K 4log(K)H 2 ∆ 2 j .
(3)

With probability at least 1 -K 3 exp ( - n-K 32log(K)H2
)

, for any k ∈ {1, . . . , K -1}, the dismissed arm

ℓ k = A k+1 \ A k at the end of phase k satisfies 1 4 ∆ (K+1-k) ≤ 1 2 ∆ ℓ k ≤ max m∈A k X m,n k -X ℓ k ,n k ≤ 3 2 ∆ ℓ k ≤ 3∆ (K+1-k) . ( 4
)
Proof: We consider the event E on which for any k ∈ {1, . . . , K -1}, for any arm ℓ in the worst k arms, and any arm j such that 2∆ j ≤ ∆ ℓ , we have

X j,n k -X ℓ,n k > 0.
This event holds with probability at least

1 -K 3 2 exp ( -n-K 4log(K)H2
) , since, from Hoeffding's inequality, a union bound and (2), we have

K-1 ∑ k=1 ∑ ℓ∈{(K),(K-1),...,(K+1-k)} j:2∆j ≤∆ ℓ P ( X j,n k -X ℓ,n k ≤ 0 ) ≤ K-1 ∑ k=1 ∑ ℓ∈{(K),(K-1),...,(K+1-k)} j:2∆j ≤∆ ℓ exp ( -n k (∆ ℓ -∆ j ) 2 ) ≤ K-1 ∑ k=1 kK exp ( -n k ∆ 2 (K+1-k) 4 ) ≤ K 3 2 exp ( - n -K 4log(K)H 2 ) .
During phase k, at least one of the k worst arms is surviving. On the event E, this surviving arm has an empirical mean at the end of the phase which is smaller than the one of any arm j satisfying 2∆ j ≤ ∆ (K+1-k) . So, at the end of phase k, any arm j satisfying 2∆ j ≤ ∆ (K+1-k) cannot be dismissed. Now, for a given arm j, we consider two cases depending whether there exists m ∈ {1, . . . , K} such that ∆ (m-1) ≤ 2∆ j ≤ ∆ (m) . First case. If no such m exists, then we have

∆ 2 j T j (n) ≥ 1 4 ∆ 2 (K) n 1 ≥ n-K 4log(K)H2
, so that (3) holds. Second case. If such m exists, then, from the above argument, the arm j cannot be dismissed before the end of the phase K + 2 -m (since there exists K + 1 -m arms ℓ such that ∆ ℓ ≥ 2∆ j ). From (2), we get

∆ 2 j T j (n) ≥ ∆ 2 j n K+2-m ≥ ∆ 2 j ∆ 2 (m-1) n -K log(K)H 2 ≥ n -K 4log(K)H 2 ,
which ends the proof of (3). We have seen that at the end of phase k, any arm j satisfying 2∆ j ≤ ∆ (K+1-k) cannot be dismissed. Consequently, at the end of phase k, the dismissed arm

ℓ k = A k+1 \ A k satisfies the left inequality of 1 2 ∆ (K+1-k) ≤ ∆ ℓ k ≤ 2∆ (K+1-k) . ( 5
)
Let us now prove the right inequality by contradiction. Consider k such that 2∆ (K+1-k) < ∆ ℓ k . Arm ℓ k thus belongs to the k -1 worst arms. Hence, in the first k -1 rejects, say at the end of phase k ′ , an arm j with ∆ j ≤ ∆ (K+1-k) is dismissed. From the left inequality of (5), we get

∆ (K+1-k ′ ) ≤ 2∆ j < ∆ ℓ k .
On the event E, we thus have X j,n k ′ -X ℓ k ,n k ′ > 0 (since ℓ k belongs to the k ′ worst arms by the previous inequality). This contradicts the fact that j is rejected at phase k ′ . So (5) holds. Now let E ′ be the event on which for any arm j, and any k ∈ {1, . . . , K -1} X j,n k -µ j ≤ ∆ (K+1-k) 8

. Using again Hoeffding's inequality, a union bound and (2), this event holds with probability at least 1 -2K(K -1) exp ( -

n-K 32log(K)H2
)

. We now work on the event E ∩ E ′ , which holds with probability at least

1 -K 3 exp ( - n-K 32log(K)H2
)

. From (5), the dismissed arm ℓ k at the end of phase k satisfies

X ℓ k ,n k -µ ℓ k ≤ ∆ (K+1-k) 8 ≤ ∆ ℓ k 4 .
Besides, we also have

max m∈A k X m,n k -µ (1) ≤ ∆ (K+1-k) 8 ≤ ∆ ℓ k 4 .
Consequently, at the end of phase k, we have

1 4 ∆ (K+1-k) ≤ 1 2 ∆ ℓ k ≤ max m∈A k X m,n k -X ℓ k ,n k ≤ 3 2 ∆ ℓ k ≤ 3∆ (K+1-k) .

Lower bound

In this section, we provide a general and somewhat surprising lower bound. We prove that, when the reward distributions are Bernoulli distributions with variances bounded away from 0, then for any forecaster, one can permute the distributions on the arms (before the game starts) so that the probability of missing the best arm will be at least of order exp ( -cn

H2

)

. Note that, in this formulation, we allow the forecaster to know the reward distributions up to a permutation of the indices! However, as the lower bound expresses it, whatever Bernoulli distributions with variances bounded away from 0 are considered, the quantity H 2 is a good measure of the hardness of finding the best arm.

Theorem 4 (Lower Bound) Let ν 1 , . . . , ν K be Bernoulli distributions with parameters in [p, 1 -p], p ∈ (0, 1/2). For any forecaster, there exists a permutation σ : {1, . . . , K} → {1, . . . , K} such that the probability error of the forecaster on the bandit problem defined by ν1 = ν σ(1) , . . . , νK = ν σ(K) satisfies e n ≥ exp

( - (5 + o(1))n p(1 -p)H 2 ) ,
where the o(1) term depends only on K, p and n and goes to 0 when n goes to infinity (see the end of the proof).

The proof of this result is quite technical. However, it is simple to explain why we can expect such a bound to hold. Assume (without loss of generality) that the arms are ordered, i.e., µ 1 > µ 2 ≥ . . . ≥ µ K , and that all rewards X i,t , t ∈ {1, . . . , n}, i ∈ {1, . . . , K}, are drawn before the game starts. Let i ∈ {2, . . . , K}. If X 1,n/i < X i,n/i ≤ X j,n/i for all j ∈ {2, . . . , i -1}, then it seems reasonable that a good forecaster should not pull arm 1 more than n/i times, and furthermore not select it as its recommendation. One can see that, the probability of the event we just described is of order of exp(-c(n/i)∆ 2 i ). Thus, with probability at least exp(-cn/ max 2≤i≤K i∆ -2 i ), the forecaster makes an error, which is exactly the lower bound we propose. However, note that this argument does not yield a reasonable proof strategy, in particular we assumed a "good" forecaster with a "reasonable" behavior. For instance, it is obvious that the proof has to permute the arms, since a forecaster could, despite all, choose arm 1 as its recommendation, which imply a probability error of 0 as soon as the best arm is in position 1.

The main idea of our proposed proof goes as follows. A bandit problem is defined by a product distribution ν = ν 1 ⊗ • • • ⊗ ν K . One can consider that at the beginning of the game, n K-tuples of rewards are sampled from this product distribution. This defines a table of nK rewards. A forecaster will explore a subpart of this table. We want to find a permutation σ of {1, . . . , K} so that the indices of the best arm for ν and ν = ν σ(1) ⊗ • • • ⊗ ν σ(K) are different and such that the likelihood ratio of the explored part of the table of nK rewards under ν and ν is at least of order exp(-cn/H 2 ) with probability with respect to ν ⊗n lower bounded by a positive numerical constant. This would imply the claimed bound. Remark that, the "likelihood cost" of moving distribution ν i to arm j depends on both the (Kullback-Leibler) distance between the distributions ν i and ν j , and the number of times arm j is pulled. Thus, we have to find the right trade-off between moving a distribution to a "close" distribution, and the fact that the target arm should not be pulled too much. To do this, we "slice" the set of indices in a non-trivial (and non-intuitive) way. This "slicing" depends only on the reward distributions, and not on the considered forecaster. Then, to put it simply, we move the less drawn arm from one slice to the less drawn arm in the next slice. Note that the preceding sentence is not well defined, since by doing this we would get a random permutation (which of course does not make sense to derive a lower bound). However, at the cost of some technical difficulties, it is possible to circumvent this issue.

To achieve the program outlined above, as already hinted, we use the Kullback-Leibler divergence, which is defined for two probability distributions ρ, ρ ′ on [0, 1] with ρ absolutely continuous with respect to ρ ′ as:

KL(ρ, ρ ′ ) = ∫ 1 0 log ( dρ dρ ′ (x) ) dρ(x) = E X∼ρ log ( dρ dρ ′ (X)
) .

Another quantity of particular interest for our analysis is

KL i,t (ρ, ρ ′ ) = ∑ t s=1 log ( dρ dρ ′ (X i,s )
)

. In particular, note that, if arm i has distribution ρ, then this quantity represents the (non re-normalized) empirical estimation of KL(ρ, ρ ′ ) after t pulls of arm i. Let P ν and E ν the probability and expectation signs when we integrate with respect to the distribution ν ⊗n . Another important property is that for any two product distributions ν, ν ′ , which differ only on index i, and for any event A, one has:

P ν (A) = E ν ′ 1 A exp ( -KL i,Ti(n) (ν ′ i , ν i ) ) , (6)
since we have

∏ Ti,n s=1 dνi dν ′ i (X i,s ) = exp ( -KL i,Ti(n) (ν ′ i , ν i )
) . Proof: First step: Notations. Without loss of generality we can assume that ν is ordered in the sense that

µ 1 > µ 2 ≥ . . . ≥ µ K . Moreover let L ∈ {2, . . . , K} such that H 2 = L/∆ 2 L , that is for all i ∈ {1, . . . , K}, i/∆ 2 i ≤ L/∆ 2 L . ( 7
)
We define now recursively the following sets. Let

k 1 = 1, Σ 1 = { i : µ L ≤ µ i ≤ µ L + ∆ L L 1/2 k 1 } ,
and for j > 1,

Σ j = { i : µ L + ∆ L L 1/2 k j-1 < µ i ≤ µ L + ∆ L L 1/2 k j } ,
where k j is the smallest integer (if it exists, otherwise set

k j = +∞) such that |Σ j | > 2|Σ j-1 |. Let ℓ = max{j : k j < +∞}.
We define now the random variables Z 1 , . . . , Z ℓ corresponding to the indices of the less sampled arms of the respective slices Σ 1 , . . . , Σ ℓ : for j ∈ {1, . . . , ℓ},

Z j ∈ argmin i∈Σj T i (n).
Finally let Z ℓ+1 ∈ argmin i∈{1,...,L}\{Jn} T i (n).

Second step: Controlling T Zj (n), j ∈ {1, . . . , ℓ + 1}. We first prove that for any j ∈ {1, . . . , ℓ},

3|Σ j | ≥ L 1- 1 2 k j+1 -1 . ( 8
)
To do so let us note that, by definition of k j+1 , we have

2|Σ j | ≥ { i : µ L + ∆ L /L 1/2 k j < µ i ≤ µ L + ∆ L /L 1/2 k j+1 -1 } ≥ { i : µ i ≤ µ L + ∆ L /L 1/2 k j+1 -1 } -(|Σ 1 | + . . . + |Σ j-1 |).
Now remark that, by definition again, we have

|Σ 1 | + . . . + |Σ j-1 | ≤ (2 -(j-1) + . . . + 2 -1 )|Σ j | ≤ |Σ j |.
Thus we obtain

3|Σ j | ≥ { i : µ i ≤ µ L + ∆ L /L 1/2 k j+1 -1 } .
We finish the proof of ( 8) with the following calculation, which makes use of (7). For any v ≥ 1,

|{i : µ i ≤ µ L + ∆ L /v}| = |{i : ∆ i ≥ ∆ L (1 -1/v)}| ≥ { i : √ i L ∆ L ≥ ∆ L (1 -1/v) } = { i : i ≥ L(1 -1/v) 2 } ≥ L ( 1 -(1 -1/v) 2 ) ≥ L/v.
Now (8) directly entails (since a minimum is smaller than an average), for j ∈ {1, . . . , ℓ},

T Zj (n) ≤ 3L 1 2 k j+1 -1 -1 ∑ i∈Σj T i (n). (9)
Besides, since Z ℓ+1 is the less drawn arm among L -1 arms, we trivially have

T Z ℓ+1 (n) ≤ n L -1 . (10) Third step: A change of measure. Let ν ′ = ν L ⊗ ν 2 ⊗ • • • ⊗ ν K be
a modified product distribution where we replaced the best distribution by ν L . Now let us consider the event

C n = { ∀t ∈{1, . . . , n}, i ∈ {2, . . . , L}, j ∈ {1, . . . , L}, KL i,t (ν i , ν j ) ≤ t KL(ν i , ν j ) + o n and KL 1,t (ν L , ν j ) ≤ t KL(ν L , ν j ) + o n } ,
where o n = 2 log(p -1 ) √ n log(2L). From Hoeffding's maximal inequality, one can prove that we have

P ν ′ (C n ) ≥ 1/2. We thus have ∑ 1≤z1,...,z ℓ+1 ≤L P ν ′ ( C n ∩ {Z 1 = z 1 , . . . , Z ℓ+1 = z ℓ+1 } ) ≥ 1/2. Moreover note that Z 1 , . . . , Z ℓ are all distinct. Thus there exist ℓ + 1 constants z 1 , . . . , z ℓ+1 such that, for A n = C n ∩ {Z 1 = z 1 , . . . , Z ℓ+1 = z ℓ+1 }, we have P ν ′ (A n ) ≥ 1 2L × L! . ( 11
)
Since, by definition Z ℓ+1 ̸ = J n , we have

A n ⊂ {J n ̸ = z ℓ+1 }. ( 12
)
In the following we treat differently the cases z ℓ+1 = 1 and z ℓ+1 ̸ = 1. First, let us assume that z ℓ+1 = 1.

Then, an application of ( 6) and ( 12) directly gives, by definition of A n ,

e n (ν) = P ν (J n ̸ = 1) = E ν ′ 1 Jn̸ =1 exp ( -KL 1,T1(n) (ν L , ν 1 )
)

≥ E ν ′ 1 An exp ( -KL 1,T1(n) (ν L , ν 1 )
)

≥ E ν ′ 1 An exp ( -o n -T Z ℓ+1 (n)KL(ν L , ν 1 ) ) ≥ 1 2L × L! exp ( -o n - n L -1 KL(ν L , ν 1 )
) ,

where we used (10) and ( 11) for the last equation. Now, for any p, q ∈ [0, 1], the KL divergence between Bernoulli distributions of parameters p and q satisfies KL(Ber(p), Ber(q)) ≤ (p -q) 2 q(1 -q) . (

This can be seen by using log u ≤ u -1 on the two logarithmic terms in KL(Ber(p), Ber(q)). In particular, it implies KL(ν L , ν 1 ) ≤ ∆ 2 L p(1-p) , which concludes the proof in the case z ℓ+1 = 1. Assume now that z ℓ+1 ̸ = 1. In this case we prove that the lower bound holds for a well defined permuted product distribution ν of ν. We define it as follows. Let m be the smallest j ∈ {1, . . . , ℓ + 1} such that z m = z ℓ+1 . Now we set ν as follows: νzm = ν 1 , νzm-1 = ν zm , . . ., νz1 = ν z2 , ν1 = ν z1 , and νj = ν j for other values of j in {1, . . . , K}. Remark that ν is indeed the result of a permutation of the distributions of ν. Again, an application of ( 6) and ( 12) gives, by definition of A n ,

e n (ν) = P ν (J n ̸ = z m ) = E ν ′ 1 Jn̸ =zm exp ( -KL 1,T1(n) (ν L , ν z1 ) - m-1 ∑ j=1 KL zj ,Tz j (n) (ν zj , ν zj+1 ) -KL zm,Tz m (n) (ν zm , ν z1 )
)

≥ E ν ′ 1 An exp ( -(m + 1)o n -T 1 (n)KL(ν L , ν Z1 ) - m-1 ∑ j=1 T Zj (n)KL(ν Zj , ν Zj+1 ) -T Zm (n)KL(ν Zm , ν Z1 ) ) . ( 14
)
From ( 13), the definition of Σ j , and since the parameters of the Bernoulli distributions are in [p, 1 -p],

we have

KL(ν L , ν Z1 ) ≤ 1 p(1-p) ∆ 2 L L , KL(ν Zm , ν Z1 ) ≤ ∆ 2 L p(1-p)
, and for any j ∈ {1, . . . , m -1},

KL(ν Zj , ν Zj+1 ) ≤ 1 p(1 -p) ( ∆ L L 1/2 k j+1
) 2 .

Reporting these inequalities, as well as ( 9), ( 10) and ( 11) in ( 14), we obtain:

e n (ν) ≥ E ν ′ 1 An exp ( -(m + 1)o n -3 ∆ 2 L p(1 -p)L ( T 1 (n) + m-1 ∑ j=1 ∑ i∈Σj T i (n) + nL 3(L -1) )) ≥ 1 2L × L! exp ( -L o n -3n ∆ 2 L p(1 -p)L ( 1 + L 3(L -1) )) Since L ≤ K and 2K × K! ≤ exp ( 2K log(K) )
and from the definitions of o n and L, we obtain

e n (ν) ≥ exp ( -2K log(K) -2K log(p -1 ) √ n log(2K) -5 n p(1 -p)H 2 ) ,
which concludes the proof.

6 Proofs 6.1 Proof of Inequalities (1)

Let log(K) = 1 2 + ∑ K i=2 1 i . Remark that log(K +1)-1/2 ≤ log(K) ≤ log(K)+1/2 ≤ log(2K).
Precisely, we will prove

H 2 ≤ H 1 ≤ log(K) H 2 ,
which is tight to the extent that the right inequality is an equality when for some 0 < c ≤ 1/ √ K, we have

∆ (i) =
√ ic for any i ̸ = i * , and the left inequality is an equality if all ∆ i 's are equal. Proof: The left inequality follows from: for any i ∈ {1, . . . , K}, H

1 = ∑ K k=1 ∆ -2 (k) ≥ ∑ i k=1 ∆ -2 (i) ≥ i∆ -2 (i) . The right inequality directly comes from ∑ K i=1 ∆ -2 (i) = ∆ -2 (2) + ∑ K i=2 1 i i∆ -2 (i) ≤ log(K) max i∈{1,...,K} i∆ -2 (i) .
6.2 Proof of Theorem 1

First step. Let us consider the event

ξ = { ∀i ∈ {1, . . . , K}, s ∈ {1, . . . , n}, | X i,s -µ i | < 1 5 √ a s } .
From Hoeffding's inequality and a union bound, we have

P(ξ) ≥ 1 -2nK exp ( -2a 25
) . In the following, we prove that on the event ξ we have J n = i * , which concludes the proof. Since J n is the empirical best arm, and given that we are on ξ, it is enough to prove that

1 5 √ a T i (n) ≤ ∆ i 2 , ∀i ∈ {1, . . . , K},
or equivalently:

T i (n) ≥ 4 25 a ∆ 2 i , ∀i ∈ {1, . . . , K}. ( 15
)
Second step. Firstly we prove by induction that

T i (t) ≤ 36 25 a ∆ 2 i + 1, ∀i ̸ = i * . ( 16
)
It is obviously true at time t = 1. Now assume that the formula is true at time t -1. Parameter: exploration rate c > 0.

If I t ̸ = i then T i (t) = T i (t -
Definitions:

For k ∈ {1, . . . , K -1}, let n k = ⌈ 1 log(K) n-K K+1-k ⌉ , t0 = 0, t1 = Kn1, and for k > 1, t k = n1 + . . . n k-1 + (K -k + 1)n k .
For i ∈ {1, . . . , K} and a > 0, let Bi,s(a) = Xi,s + √ a s for s ≥ 1 and Bi,0 = +∞. Algorithm: For each phase k = 0, 1, . . . , K -1: Let H 1,k = K if k = 0, and otherwise

H 1,k = max K-k+1≤i≤K i ∆ -2 <i> , where ∆i = ( max 1≤j≤K X j,T j (t k ) ) -X i,T i (t k ) and < i > is an ordering such that ∆<1> ≤ . . . ≤ ∆<K>. For t = t k + 1, . . . , t k+1 : Draw It ∈ argmax i∈{1,...,K} B i,T i (t-1) (cn/ H 1,k ). Recommendation: Let Jn ∈ argmax i∈{1,...,K} X i,T i (n) .
Figure 4: Adaptive UCB-E algorithm. Its intuitive justification goes as follows: The time points t k correspond to the moments where the Successive Rejects algorithm would dismiss an arm. Intuitively, in light of Theorem 3, one can say that at time t k a good algorithm should have reasonable approximation of the gaps between the best arm and the k worst arms, that is the quantities ∆ (K-k+1) , . . . , ∆ (K) . Now with these quantities, one can build a lower estimate of H 2 and thus also of H 1 . We use this estimate between the time points t k and t k+1 to tune the parameter a of UCB-E. Now we prove an other useful formula:

T i (t) ≥ 4 25 min ( a ∆ 2 i , 25 36 (T i * (t) -1) ) , ∀i ̸ = i * . ( 17
)
With the same inductive argument as the one to get equation ( 16), we only need to prove that this formula holds when I t = i * . By definition of the algorithm, and since we are on ξ, when I t = i * we have for all i:

µ * + 6 5 √ a T i * (t -1) ≥ µ i + 4 5 √ a T i (t -1) , which implies T i (t -1) ≥ 16 25 a ( ∆ i + 6 5 √ a T i * (t-1) ) 2 .
We then obtain (17) by using u + v ≤ 2 max(u, v), T i (t) = T i (t -1) and T i * (t -1) = T i * (t) -1. Third step. Recall that we want to prove equation (15). From (17), we only have to show that

25 36 (T i * (n) -1) ≥ a ∆ 2 i *
, where we recall that ∆ i * is the minimal gap ∆ i * = min i̸ =i * ∆ i . Using equation ( 16) we obtain:

T i * (n) -1 = n -1 - ∑ i̸ =i * T i (n) ≥ n -K - 36 25 a ∑ i̸ =i * ∆ -2 i ≥ 36 25 a∆ -2 i * ,
where the last inequality uses 36 25 H 1 a ≤ n -K. This concludes the proof.

Experiments

We propose a few simple experiments to illustrate our theoretical analysis. As a baseline comparison we use the Hoeffding Race algorithm, see [START_REF] Maron | Hoeffding races: Accelerating model selection search for classification and function approximation[END_REF], and the uniform strategy, which pulls equally often each arm and recommend the arm with the highest empirical mean, see [START_REF] Bubeck | Pure exploration in multi-armed bandits problems[END_REF] for its theoretical analysis. We consider only Bernoulli distributions, and the optimal arm always has parameter 1/2. Each experiment corresponds to a different situation for the gaps, they are either clustered in few groups, or distributed according to an arithmetic or geometric progression. In each experiment we choose the number of samples (almost) equal to H 1 (except for the last experiment where we run it twice, the second time with 2H 1 samples). If our understanding of the meaning of H 1 is sound, in each experiment the strategies SR and UCB-E should be able to find the best arm with a reasonable probability (which should be roughly of the same order in each experiment). We report our results in Figure 5. The parameters for the experiments are as follows:

-Experiment 1: One group of bad arms, K = 20, µ 2:20 = 0.4 (meaning for any j ∈ {2, . . . , 20}, µ j = 0.4) -Experiment 2: Two groups of bad arms, K = 20, µ 2:6 = 0.42, µ 7:20 = 0.38.

-Experiment 3: Geometric progression, K = 4, µ i = 0.5 -(0.37) i , i ∈ {2, 3, 4}.

-Experiment 4: 6 arms divided in three groups, K = 6, µ 2 = 0.42, µ 3:4 = 0.4, µ 5:6 = 0.35.

-Experiment 5: Arithmetic progression, K = 15, µ i = 0.5 -0.025i, i ∈ {2, . . . , 15}.

-Experiment 6: Two good arms and a large group of bad arms, K = 20, µ 2 = 0.48, µ 3:20 = 0.37.

-Experiment 7: Three groups of bad arms, K = 30, µ 2:6 = 0.45, µ 7:20 = 0.43, µ 21:30 = 0.38.

The different graphics should be read as follows: Each bar represents a different algorithm and the bar's height represents the probability of error of this algorithm. The correspondence between algorithms and bars is the following:

-Bar 1: Uniform sampling strategy.

-Bar 2-4: Hoeffding Race algorithm with parameters δ = 0.01, 0.1, 0.3.

-Bar 5: Successive Rejects strategy.

-Bar 6-9: UCB-E with parameter a = cn/H 1 where respectively c = 1, 2, 4, 8.

-Bar 10-14: Adaptive UCB-E (see Figure 4) with parameters c = 1/4, 1/2, 1, 2, 4.

Conclusion

This work has investigated strategies for finding the best arm in a multi-armed bandit problem. It has proposed a simple parameter-free algorithm, SR, that attains optimal guarantees up to a logarithmic term (Theorem 2 and Theorem 4). A precise understanding of both SR (Theorem 3) and a UCB policy (Theorem 1) lead us to define a new algorithm, Adaptive UCB-E. It comes without guarantee of optimal rates (see end of Section 3), but performs better than SR in practice (for c = 1, Adaptive UCB-E outperformed SR on all the experiments we did, even those done to make it fail). One possible explanation is that SR is too static: it does not implement more data driven arguments such as: in a phase, a surviving arm performing much worse than the other ones is still drawn until the end of the phase even if it is clear that it is the next dismissed arm. Extensions of this work may concentrate on the following problems. (i) What is a good measure of hardness when one takes into account the (empirical) variances? Do we have a good scaling with respect to the variance with the current algorithms or do we need to modify them ? (ii) Is it possible to derive a natural anytime version of Successive Rejects (without using a doubling trick)? (iii) Is it possible to close the logarithmic gap between the lower and upper bounds? (iv) How should we modify the algorithm and the analysis if one is interested in recommending the top m actions instead of a single one?

Figure 3 :

 3 Figure 3: SR (Successive Rejects) algorithm.

  1) and the formula still holds. On the other hand, if I t = i, then in particular it means that B i,Ti(t-1) ≥ B i * ,T i * (t-1) . Moreover, since we are on ξ, we have B i * ,T i * (t-1) ≥ µ * and B i,Ti(t-1) ≤

	µ i + 6 5	√	a Ti(t-1) . Thus, we have 6 5	√

a Ti(t-1) ≥ ∆ i . By using T i (t) = T i (t -1) + 1, we obtain (16).

We say "in general" to rule out some trivial cases (like when the reward distributions are all Dirac distributions) in which the probability of error en would be much smaller.

: These results support our theoretical findings in the following sense: Despite the fact that the experiments are very different, one can see that since we use a number of samples (almost) equal to the hardness H 1 , in all of them we get a probability of error of the same order, and moreover this probability is small enough to say that we identified the best arm. Note that the Successive Rejects algorithm represents in all cases a substantial improvement over both the naive uniform strategy and Hoeffding Race. These results also justify experimentally the algorithm Adaptive UCB-E. Indeed one can see that with the constant c = 1, we obtain better results than SR in all experiments, even in experiment 6 which was designed to be a difficult instance of Adaptive UCB-E.


R Grompone Von Gioi

P Monasse

J.-M Morel

Z Tang

TOWARDS HIGH-PRECISION LENS DISTORTION CORRECTION

Keywords: Planar textured pattern, virtual pinhole camera, nonparametric, high-precision, lens distortion

This paper points out and attempts to remedy a serious discrepancy in results obtained by global calibration methods:

The re-projection error can be rendered very small by these methods, but we show that the optical distortion correction is far less accurate. This discrepancy can only be explained by internal error compensations in the global methods that leave undetected the inadequacy of the distortion model. This fact led us to design a model-free distortion correction method where the distortion can be any image domain diffeomorphism. The obtained precision compares favorably to the distortion given by state of the art global calibration and reaches a RMSE of 0.08 pixels. Nonetheless, we also show that this accuracy can still be improved.

INTRODUCTION

This paper presents a small step forward in a research programme whose aim is to devise a highly accurate camera calibration method. By highly accurate, we mean that the residual error between the camera and its numerical model obtained by calibration should be far smaller than the pixel size. At first sight, this problem seemed to have been solved adequately by recent global calibration methods. The celebrated Lavest et al. method [START_REF] Lavest | Do We Really Need an Accurate Calibration Pattern to Achieve a Reliable Camera Calibration?[END_REF] measures the non-flatness of a pattern and yields a remarkably small re-projection error of about 0.02 pixels, which outperforms the precision of other methods. The experiments described below will actually confirm this figure. For the goals of computer vision, this precision would be more than sufficient. Yet, this paper describes a seriously discrepant accuracy measurement contradicting this hasty conclusion. As in [START_REF] Alvarez | An algebraic approach to lens distortion by line rectification[END_REF] we shall use the following fact: Theorem 1 [START_REF] Faugeras | Straight lines have to be straight[END_REF] a camera follows the pinhole model if and only if the projection of every straight line in space onto the camera is a line. Following this theorem, the accuracy criterion used herewith directly measures the straightness of straight lines. We shall see that this straightness criterion gives a RMSE in the order of 0.2 pixel, which contradicts the 0.02 pixel re-projection accuracy. This significant discrepancy means that, in the global optimization process, errors in the external and internal camera parameter are being compensated by opposite errors in the distortion model. Thus, an inaccurate distortion model can pass undetected. Such facts raise a solid objection to global calibration methods, which estimate simultaneously the lens distortion and the camera parameters. This paper reconsiders the whole calibration chain and examines an alternative way to guarantee a high accuracy. A useful tool toward this goal will be proposed and carefully tested. It is a direct nonparametric distortion correction method. By nonparametric, we mean that the distortion model allows for any diffeomorphism.

We shall follow the usual assumption that a real camera deviates from the ideal pinhole model [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF] by a lens distortion [START_REF] Brown | Close-range camera calibration[END_REF]. Thus, distortion correction is a fundamental step in multiview geometry applications such as 3D reconstruction. The above mentioned error measurement discrepancy may explain why three categories of distortion correction methods still coexist: classic pattern-based methods [START_REF] Tsai | A versatile camera calibration technique for highaccuracy 3d machine vision metrology using off-the-shelf tv cameras and lenses[END_REF][START_REF] Zhang | Flexible camera calibration by viewing a plane from unknown orientations[END_REF][START_REF] Lavest | Do We Really Need an Accurate Calibration Pattern to Achieve a Reliable Camera Calibration?[END_REF], plumb lines methods [START_REF] Faugeras | Straight lines have to be straight[END_REF][START_REF] Hlaváč | Correcting radial lens distortion without knowledge of 3-d structure[END_REF][START_REF] Prescott | Line-based correction of radial lens distortion[END_REF][START_REF] Alvarez | An algebraic approach to lens distortion by line rectification[END_REF] whose goal is to correct the distortion only, and methods based on enlarged epipolar geometry [START_REF] Li | A non-iterative method for correcting lens distortion from nine point correspondences[END_REF][START_REF] Fitzgibbon | Simultaneous linear estimation of multiple view geometry and lens distortion[END_REF][START_REF] Kukelova | Two minimal problems for cameras with radial distortion[END_REF][START_REF] Stein | Lens distortion calibration using point correspondences[END_REF]]. Yet, these methods are all parametric and depend on the a priori choice of a distortion model with a fixed number of parameters. This per se is a drawback: such calibration methods require several trials and a manual model selection. Most methods assume a radial distortion modeled as a low-order even polynomial [START_REF] Tsai | A versatile camera calibration technique for highaccuracy 3d machine vision metrology using off-the-shelf tv cameras and lenses[END_REF][START_REF] Ma | Implicit and explicit camera calibration: theory and experiments[END_REF] (models with more parameters are needed for wide-angle or fish-eye lens [START_REF] Brandt | A generic camera model and calibration method for conventional, wide-angle, and fish-eye lenses[END_REF][START_REF] Faugeras | Straight lines have to be straight[END_REF]).

The proposed distortion correction method does not belong to any of the three categories. Indeed, it is nonparametric, non-iterative and model-free. Like most methods in the second and third category, the method decouples the distortion estimation from the calibration of camera internal parameters, thus avoiding any error compensation between them [START_REF] Herniou | Camera calibration with distortion models and accuracy evaluation[END_REF][START_REF] Li | A non-iterative method for correcting lens distortion from nine point correspondences[END_REF].

Our plan is as follows. Section 2 gives the necessary definitions of the real camera and the pinhole model. It explains why a distortion correction up to a homography is sufficient for 3D applications, and defines the concept of virtual pinhole camera. The proposed nonparametric distortion correction is detailed in section 3, and is followed by experimental results in section 4. The last section 5 discusses how the high accuracy quest could be pursued.

THE VIRTUAL PINHOLE CAMERA

The classic camera model is

Ĉ := DKR[I | -T ] ( 1
)
where T is a 3D point representing the optic center of camera, R is a 3D rotation representing camera orientation, K is a 3 × 3 upper triangular matrix containing the camera internal parameters, and D is a diffeomorphism of the image domain representing the non-linear lens distortion. By pinhole camera model we mean a distortion-free camera model

C := KR[I | -T ]. (2)
The nonparametric method will estimate the distortion up to an arbitrary invertible homography H: D = DH. The correction precision evaluation will be based on the straightness of corrected lines, which is preserved by any homography. Applying D-1 on Ĉ, yields

C = D-1 DKR[I | -T ] = H -1 D -1 DKR[I |-T ] = H -1 KR[I | -T ].( 3
)
Thus, inverting the distortion on all images produced by the camera yields a new camera model which becomes pinhole. H, K being invertible, the decomposition H -1 K = KR is unique by QR decomposition. So after distortion correction (up to a homography) we have

C = KR R[I | -T ] = K R[I | -T ],
which we call the virtual pinhole camera obtained after distortion correction. The orientation and internal parameters of this virtual model do not match the physics of the actual camera, but yield a virtual pinhole camera that can be used to the very same purposes. Indeed, consider several positions of the physical camera inducing as many camera models

C i = DKR i [I | -T i ].
Applying the correction D-1 to all images obtained from these camera positions yields virtual pinhole cameras Ci = K Ri [I | -T i ], which maintains the same relative orientations. From these cameras the whole 3D scene can be reconstructed by standard methods, up to a 3D similarity.

NONPARAMETRIC DISTORTION CORRECTION

The experimental set up

A nonparametric method requires the use of a highly textured planar pattern obtained by printing a textured image and pasting it on a very flat object (a mirror was used in the experiments). Two photos of the pattern are taken by the camera in a fixed lens configuration (e.g., with fixed zoom and focus). Ideally, the whole captor must be covered by the whole pattern.

Feature points

The distortion is estimated (up to a homography) as the diffeomorphism mapping the original digital pattern to a photograph of it. This requires a dense registration, which is obtained by an improved version of the SIFT [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF] method providing high-precision localization about 0.03 pixels [START_REF] Buades | A note on multiimage denoising[END_REF]. Denote by I the original digital pattern, by P the printed pattern, and by v the photo of P . The set of SIFT points matching from I to v is denoted by I v and the corresponding points in v by v I . Since the pattern is planar, there is a planar homography H such that DHI v = v I . Knowing DH permits to synthesize a virtual pinhole camera by applying (DH) -1 on Ĉ, as shown in Eq. (3).

Triangulation and affine interpolation

The correspondences (I v , v I ) actually only define the distortion field DH on the SIFT points I v . The distortion being very smooth and the SIFT points dense enough, an affine interpolation is sufficient. This interpolation is performed after the image domain has been partitioned by a Delaunay triangulation of the SIFT points in v and I respectively (Fig. 2). Fig. 1 shows the texture pattern, selected to yield a maximal density of reproducible SIFT points at fine scales.

Outliers elimination: a loop validation

The few wrong SIFT matches (outliers) are nonetheless a serious problem. In our case, precisely because of the lens distortion, matching points are not related by a homography, and directly applying RANSAC [START_REF] Fischler | Random sample consensus: A paradigm for model fitting with applications to image analysis and automated cartography[END_REF] would not work. The problem is solved by a loop validation (Fig. 3). Consider two similar photos of the pattern, u and v (obtained by moving slightly the camera between two successive snapshots). With straightforward notation we have u I = DH u I u and v I = DH v I v (since the same camera and configuration are used, D does not change). The points v I can be projected back on I by the distortion field from u to I, obtaining I uv = (DH u ) -1 v I . It follows that I v and I uv are related by a homography (without distortion) because

I uv = (DH u ) -1 DH v I v = H -1 u H v I v . (4)
This homography can be estimated by RANSAC algorithm and all the outliers not compatible with the homography are eliminated.

v I = H H I v v u I DH v DH u DH v (DH ) u -1 I v I uv v I v I uv u -1
Fig. 3: The loop validation used to remove outliers.

3.5. Algorithm summary 1. Take two slightly different photos of a textured planar pattern with constant camera settings; 2. apply high precision SIFT between the original digital pattern and both photographs; 3. eliminate outliers by the loop validation step; 4. interpolate the remaining matches to get a dense reverse distortion field; 5. by applying the reverse distortion field to all images produced by the real camera, the camera is converted into a virtual pinhole camera.

EXPERIMENTS

The experiments were made with a Canon EOS 30D reflex camera and an EFS 18-55mm lens. The minimal focal length (18mm) was chosen to produce a fairly large distortion. The RAW images were demosaicked by summing up the four pixels of each 2×2 Bayer cell, obtaining a half-size image. Fig. 1 shows the digital pattern I that was used. Fig. 4a shows a subsampling of the resulting distortion field, after the validation loop, and Fig. 4b shows the modulus of the interpolated distortion field on the discrete image domain. The distortion field is not circular symmetric, which is natural, the distortion being estimated up to an unknown homography. To check the quality of the correction, we built a pattern with tightly stretched strings, that guarantees straightness, see 5b shows the image corrected by our nonparametric method. The lines numbered in Fig. 5a were used to evaluate the distortion error and compare it with the error left by the Lavest et al. algorithm [START_REF] Lavest | Do We Really Need an Accurate Calibration Pattern to Achieve a Reliable Camera Calibration?[END_REF]. On each corrected line, subpixel precision edge points were obtained by Devernay's algorithm [START_REF] Devernay | A Non-Maxima Suppression Method for Edge Detection with Sub-Pixel Accuracy[END_REF]. Then, their regression line was computed and the RMS (root-mean-square) distance from each edge point to the line was used as the error measure. Table 1 shows the results. The proposed nonparametric method shows an improvement over the Lavest et al. method. But this difference is not just quantitative. The Lavest et al. result is somewhat final. Indeed, it already includes a correction of the non-flatness of the pattern while the nonparametric method does not. Figure 6 plots the straightness error of the proposed method along line number one (that is, the distance between edge points to the regression line). One can see small oscillations easily explained by noise and aliasing, but also a global tendency which can only be due to the non-flatness of the pattern. A parallel deterministic tendency observed on the other lines confirms this explanation. Fig. 7 shows that a flatness error of 100 µm (the thickness of a normal paper sheet) can produce the observed tendency. This non-flatness effect is stronger near the border of the image because the angle-of-view is larger. Simple physical measurements by the classic ruler method confirmed that the pattern showed a non-flatness of this amount. resulting "distortion" error flatness error pattern camera α Fig. 7: A flatness error in the pattern can be mistaken for a "distortion". In the experiments, with α ≈ 65 • and a flatness error about 100 µm, the produced "distortion" error is about 64 µm (flatness error × tan α 2 ). Our pattern has the size 406 × 271 mm and produces a 1761 × 1174 image, then one pixel corresponds to 230 µm. Thus the observed error would be approximately 0.3 pixel.

DISCUSSION

The above experimental setting suggests two ways to eventually reach a still higher precision. The first way would simply be to use very flat patterns. But this raises the problem of making smart patterns. An alternative strategy was suggested in [START_REF] Alvarez | An algebraic approach to lens distortion by line rectification[END_REF], using accurate straight objects, like the tightly stretched strings we already used here. Yet, it seems advisable to try to keep the strength and beauty of the Lavest et al. method, which is to estimate and correct the pattern's shape by the calibration process itself. An iterative method could be envisaged where, first, the distortion is corrected by the proposed nonparametric method and, second, the physical shape of the pattern is computed by the Lavest et al. method with no distortion model. Using this correction the distortion would be recomputed, and so on. This is, however, a complex process, which will require a heavier procedure and a mathematical analysis.

Fig. 1 :Fig. 2 :

 12 Fig. 1: Digital pattern: 1761 × 1174 pixels.

Fig. 4 :

 4 Fig. 4: (a) The (subsampled) distortion field directly defined on "inlier" correspondences after the loop validation. (b) Modulus of the dense distortion field obtained by affine interpolation.

Fig. 5 : 1 :

 51 Fig. 5: Distorted lines marked by numbers in (a) will be used to evaluate the precision after distortion correction.

Fig. 5a .

 5a Fig.5a. The distortion is visible near the border of the image. Fig.5bshows the image corrected by our nonparametric method. The lines numbered in Fig.5awere used to evaluate the distortion error and compare it with the error left by the Lavest et al. algorithm[START_REF] Lavest | Do We Really Need an Accurate Calibration Pattern to Achieve a Reliable Camera Calibration?[END_REF]. On each corrected line, subpixel precision edge points were obtained by Devernay's algorithm[START_REF] Devernay | A Non-Maxima Suppression Method for Edge Detection with Sub-Pixel Accuracy[END_REF]. Then, their regression line was computed and the RMS (root-mean-square) distance from each edge point to the line was used as the error measure. Table1shows the results. The proposed nonparametric method shows an improvement over the Lavest et al. method. But this difference is not just quantitative. The Lavest et al. result is somewhat final. Indeed, it already includes a correction of the non-flatness of the pattern while the nonparametric method does not.

Fig. 6 :

 6 Fig.6: The distance in pixels from the edge points to their regression line on line 1 in Fig.5a, after correction by the proposed method.

ACKNOWLEDGMENTS

The authors wish to thank Jean-Marc Lavest and Eric Royer (Université de Clermont-Ferrand) for lending us their calibration software and their gracious help. Part of this work was funded by the Agence Nationale de la Recherche, Callisto project (ANR-09-CORD-003).