Fabio Daolio
email: fdaolio@shinshu-u.ac.jp

Arnaud Liefooghe
email: arnaud.liefooghe@inria.fr

Sébastien Verel
email: verel@lisic.univ-littoral.fr

Hernán Aguirre

Kiyoshi Tanaka
email: ktanaka@shinshu-u.ac.jp

Global vs Local Search on Multi-objective NK-Landscapes: Contrasting the Impact of Problem Features

Computationally hard multi-objective combinatorial optimization problems are common in practice, and numerous evolutionary multiobjective optimization (EMO) algorithms have been proposed to tackle them. Our aim is to understand which (and how) problem features impact the search performance of such approaches. In this paper, we consider two prototypical dominance-based algorithms: a global EMO strategy using an ergodic variation operator (GSEMO) and a neighborhood-based local search heuristic (PLS). Their respective runtime is estimated on a benchmark of combinatorial problems with tunable ruggedness, objective space dimension, and objective correlation (ρMNK-landscapes). In other words, benchmark parameters define classes of instances with increasing empirical problem hardness; we enumerate and characterize the search space of small instances. Our study departs from simple performance comparison to systematically analyze the correlations between runtime and problem features, contrasting their association with search performance within and across instance classes, for both chosen algorithms. A mixed-model approach then allows us to further generalize from the experimental design, supporting a sound assessment of the joint impact of instance features on EMO search performance.

MOTIVATIONS

Despite the increasing number of available heuristics for evolutionary multi-objective optimization, the understanding of their search dynamics is comparatively less advanced than in the singleobjective case. This knowledge is all the more valuable in a blackbox scenario, where such approaches are mostly adopted. Finding what makes a problem hard, which features have an influence on algorithm performance, and how that changes across different algorithm and problem classes, are relevant and still open challenges in multi-objective optimization.

Following [START_REF] Liefooghe | A feature-based performance analysis in evolutionary multiobjective optimization[END_REF], we precisely investigate the impact of instance features on search performance. We decline this research question in three complementary sub-questions: (i) how strongly are features associated to performance? (ii) what performance change is . to be expected when varying each feature? (iii) what is the relative importance of features in explaining the performance variance? We consider a global and a local EMO heuristics and try to highlight the differences and similarities in the behavior of both algorithm classes, with respect to a set of relevant problem features. As such, the emphasis is more on making inferences, e.g. [START_REF] Chiarandini | Mixed models for the analysis of optimization algorithms[END_REF], rather than making predictions, e.g. [START_REF] Hutter | Algorithm runtime prediction: Methods & evaluation[END_REF]. Hence, this work attempts to provide both insightful understandings and methodological suggestions.

We introduce the experimental setup of our analysis in Section 2, in particular the considered problems, ρMNK-landscapes [START_REF] Verel | On the structure of multiobjective combinatorial search space: MNK-landscapes with correlated objectives[END_REF], and the search algorithms, namely the global simple EMO optimizer (GSEMO) [START_REF] Laumanns | Running time analysis of evolutionary algorithms on a simplified multiobjective knapsack problem[END_REF] and the Pareto local search (PLS) [START_REF] Paquete | Pareto local optimum sets in the biobjective traveling salesman problem: An experimental study[END_REF], of which we measure the estimated runtime to find a (1 + ε)-approximation of the Pareto set. In Section 3, we examine the performance of GSEMO and PLS depending on the ruggedness, the objective space dimension, and the objective correlation of ρMNK-landscapes. In Section 4, we provide an overview of the features that characterize the hardness of problem instances, relating to the difficulties randomized search heuristics might have to face. Via a correlation and regression analysis, we are able to point out how such features influence the performance of global and local search algorithms like GSEMO and PLS. The last section concludes the paper.

INGREDIENTS

Multi-objective Optimization

We are interested in maximizing a black-box objective function vector f : X → Z, that maps any solution from the solution space x ∈ X to a vector in the objective space z ∈ Z, with Z ⊆ IR M , such that z = f (x). We assume that the solution space is a discrete set X = {0, 1} N , where N is the problem size, i.e. the number of binary (zero-one) variables. An objective vector z ∈ Z is dominated by an objective vector z ∈ Z, denoted by z ≺ z , iff for all i ∈ {1, . . . , M } zi z i , and there is a j ∈ {1, . . . , M } such that zj < z j . Similarly, a solution x ∈ X is dominated by a solution x ∈ X iff f (x) ≺ f (x ). An objective vector z ∈ Z is non-dominated if there does not exist any objective vector z ∈ Z such that z ≺ z. A solution x ∈ X is non-dominated, or Paretooptimal, if f (x) is non-dominated. The set of Pareto-optimal solutions is the Pareto set (PS) ; its mapping in the objective space is the Pareto front (PF). One of the main challenges in multi-objective optimization is to identify the Pareto set/front, or a good approximation of it, for large-size and difficult problems.

ρMNK-Landscapes

The family of ρMNK-landscapes constitutes a problem-independent model used for constructing multi-objective multimodal land-scapes with objective correlation [START_REF] Verel | On the structure of multiobjective combinatorial search space: MNK-landscapes with correlated objectives[END_REF]. They extend single-objective NK-landscapes [START_REF] Kauffman | The Origins of Order[END_REF] and multi-objective NK-landscapes with independent objective functions [START_REF] Aguirre | Working principles, behavior, and performance of MOEAs on MNK-landscapes[END_REF]. Candidate solutions are binary strings of size N , i.e. the solution space is X = {0, 1} N . The objective function vector f = (f1, . . . , fi, . . . , fM ) is defined as f : {0, 1} N → [0, 1] M such that each objective function fi is to be maximized. As in the single-objective case, each separate objective function value fi(x) of a solution x = (x1, . . . , xj, . . . , xN ) is an average value of the individual contributions associated with each variable xj. Indeed, for each objective fi, i ∈ {1, . . . , M }, and each variable xj, j ∈ {1, . . . , N }, a component function fij : {0, 1} K+1 → [0, 1] assigns a real-valued contribution for every combination of xj and its K epistatic interactions {xj1, . . . , xjK }. These fij-values are uniformly distributed in the range [0, 1]. Thus, the individual contribution of a variable xj depends on its value and on the values of K < N other variables {xj 1 , . . . , xj K }. The problem can be formalized as follows.

max fi(x) = 1 N N j=1 fij(xj, xj 1 , . . . , xj K ) i ∈ {1, . . . , M } s.t. xj ∈ {0, 1} j ∈ {1, . . . , N }
In this work, the epistatic interactions, i.e. the K variables that influence the contribution of xj, are set uniformly at random among the (N -1) variables other than xj, following the random neighborhood model from [START_REF] Kauffman | The Origins of Order[END_REF]. By increasing the number of epistatic interactions K from 0 to (N -1), problem instances can be gradually tuned from smooth to rugged. In ρMNK-landscapes, fijvalues additionally follow a multivariate uniform distribution of dimension M , defined by an M × M positive-definite symmetric covariance matrix (cpq) such that cpp = 1 and cpq = ρ for all p, q ∈ {1, . . . , M } with p = q, were ρ > -1 M -1 defines the objective correlation degree; see [START_REF] Verel | On the structure of multiobjective combinatorial search space: MNK-landscapes with correlated objectives[END_REF] for details. The positive (respectively negative) data correlation ρ decreases (respectively increases) the degree of conflict between the objective function values. The same correlation coefficient ρ is then defined between all pairs of objectives, and the same epistatic degree K and epistatic interactions are set for all the objectives.

Multi-objective Search Heuristics

In this paper, we consider two randomized search heuristics for ρMNK-landscapes: (i) global SEMO (GSEMO) [START_REF] Laumanns | Running time analysis of evolutionary algorithms on a simplified multiobjective knapsack problem[END_REF], a simple elitist steady-state global EMO algorithm (see Algorithm 1); and (ii) Pareto local search (PLS) [START_REF] Paquete | Pareto local optimum sets in the biobjective traveling salesman problem: An experimental study[END_REF], a population-based multi-objective local search (Algorithm 2). Both algorithms maintain an unbounded archive A of mutually non-dominated solutions. This archive is initialized with one random solution from the solution space. Then, at each iteration, one solution is selected at random from the archive x ∈ A. In GSEMO, each binary variable from x is independently flipped with rate 1 N in order to produce an offspring solution x . The archive is then updated by keeping the non-dominated solutions from A ∪ {x }. In PLS, the solutions located in the neighborhood of x are evaluated. Let N(x) denote the set of solutions located at a Hamming distance 1. The non-dominated solutions from A ∪ N(x) are stored in the archive, and the current solution x is then tagged as visited in order to avoid a useless revaluation of its neighborhood. This process is iterated until a stopping condition is satisfied. While for GSEMO there does not exists any explicit stopping rule [START_REF] Laumanns | Running time analysis of evolutionary algorithms on a simplified multiobjective knapsack problem[END_REF], PLS has a natural stopping condition which is satisfied when all solutions from the archive are tagged as visited.

In other words, while PLS is based on the exploration of the whole 1-bit-flip neighborhood from x, GSEMO rather uses an ergodic operator, i.e. an independent bit-flip mutation. This means Algorithm 1: GSEMO that there is a non-zero probability of reaching any solution from the solution space at every GSEMO iteration. This makes GSEMO a global optimizer rather than a local optimizer as PLS. In this paper, we are interested in the running time, in terms of a number of function evaluations, until an (1 + ε)-approximation of the Pareto set is identified and is contained in the internal memory A of the algorithm, subject to a maximum budget of function evaluations.

Estimated Expected Running Time (ert)

Let ε be a constant value such that ε 0. The (multiplicative) εdominance relation ( ε) can be defined as follows. For x, x ∈ X, x is ε-dominated by x (x ε x ) iff fi(x)

(1 + ε) • fi(x ), ∀i ∈ {1, . . . , m}. A set X ε ⊆ X is an (1 + ε)-approximation of the Pareto set if for any solution x ∈ X, there is one solution x ∈ X ε such that x ε x . This is equivalent to finding an approximation set whose multiplicative epsilon quality indicator value with respect to the (exact) Pareto set is lower than (1 + ε), see e.g. [START_REF] Zitzler | Performance assessment of multiobjective optimizers: An analysis and review[END_REF]. Interestingly, under some general assumptions, there always exists an (1+ε)-approximation, for any given ε 0, whose cardinality is both polynomial in the problem size and in 1 ε [START_REF] Papadimitriou | On the approximability of trade-offs and optimal access of web sources[END_REF]. Following a conventional methodology from single-objective continuous black-box optimization [START_REF] Auger | Performance evaluation of an advanced local search evolutionary algorithm[END_REF], we measure algorithm performance in the expected number of function evaluations to identify an (1 + ε)-approximation. However, as any heuristic, GSEMO or PLS can either succeed or fail to reach an accuracy of ε in a single run. In case of success, we record the number of function evaluations until an (1 + ε)-approximation was found. In case of failure, we simply restart the algorithm at random. Thus we obtain a "simulated running time" [START_REF] Auger | Performance evaluation of an advanced local search evolutionary algorithm[END_REF] from a set of independent trials on each instance. Such performance measure allows us to take into account both the success rate ps ∈ (0, 1] and the convergence speed of the algorithm with restarts. Precisely, after (t -1) failures, each one requiring T f evaluations, and the final successful run of Ts evaluations, the total runtime is T = t-1 i=1 T f + Ts. By taking the expectation and by considering independent trials as a Bernoulli process stopping at the first success, we have:

E[T ] = 1 -ps ps E[T f ] + E[Ts]
In our case, the success rate ps is estimated with the ratio of suc-

Ti

where ts is the number of successful runs, and Ti is the number of evaluations for successful run i. For more details, we refer to [START_REF] Auger | Performance evaluation of an advanced local search evolutionary algorithm[END_REF].

EXPERIMENTAL ANALYSIS

Experimental Setup

We consider ρMNK-landscapes with an epistatic degree K ∈ {2, 4, 6, 8, 10}, an objective space dimension M ∈ {2, 3, 5}, and ρ ∈ {-0.9, -0.7, -0.4, -0.2, 0.0, 0.2, 0.4, 0.7, 0.9}, such that ρ > -1

M -1 , as objective correlation. The problem size is set to N = 18 in order to enumerate the solution space exhaustively. A set of 30 different landscapes are independently generated at random for each parameter combination ρ, M , and K. They are made available at the following URL: http://mocobench.sf.net.

The target is set at ε = 0.1. The time limit is set to Tmax = 2 N • 10 -1 < 26 215 function evaluations without identifying an (1 + ε)-approximation. Each algorithm is executed 100 times per instance. From these 100 repetitions, the success rate and the expected number of evaluations for successful runs, hence the expected runtime on the given instance, are estimated. For the comparative analysis, we only consider pairwise-complete cases, i.e. instances that have been solved by both algorithms. This brings the total number of available observations to 2 874 per algorithm.

Exploratory Analysis

The expected runtime (ert) distribution across the experimental blocks that are defined by each combination of benchmark parameters is presented in Figure 1. For both algorithms, the ert clearly increases with the ruggedness (K) and the number of objectives (M ), whereas the trend w.r.t. the objective correlation (ρ) is a bit more complex. Indeed, for a small K and a large M , the ert decreases when ρ increases. On the contrary, for large K, problem instances seem to get harder when the objectives are independent (ρ ≈ 0) rather than anti-correlated (ρ < 0), despite the fact that the cardinality of the Pareto set (PS) gets higher in such cases [START_REF] Verel | On the structure of multiobjective combinatorial search space: MNK-landscapes with correlated objectives[END_REF]. This gives rise to the inverted u-shape observed on the right side of the figure, which is particularly pronounced for PLS. Actually, we observed that the ε-value of random approximation sets follows a similar u-shaped trend w.r.t. objective correlation. Also, this ε-value tends to increase with K. This holds for approximation sets containing a constant number of randomly generated solutions. Moreover, we know from [START_REF] Verel | On the structure of multiobjective combinatorial search space: MNK-landscapes with correlated objectives[END_REF] that the number of Pareto local optimal solutions increases with K and decreases with ρ. This could explain the relative advantage of PLS on problem instances with positively correlated objectives. Notice also that the opposite is true for the connectedness of the PS: the smaller K, the larger ρ, the more clustered are Pareto optimal solutions in the solution space.

Figure 2 displays the runtime aggregated over all the instance parameters (ρ, M, K) but one. We clearly see that PLS is significantly outperforming GSEMO overall, and in particular for positively correlated objectives. In fact, the runtime of PLS is shorter than that of GSEMO for 88% of the instances; compared to PLS, on average GSEMO requires more than 17 000 additional function evaluations to identify a 1.1-approximation of the PS. The performance dif- ference between the two algorithms seems to be constant, except for large ρ and w.r.t. K. Notably, the ruggedness of the underlying single-objective objective functions appears to have the highest impact on the search performance. In particular, the ruggedness seems to have more impact on the performance of GSEMO than PLS. In general, finding an (1 + ε)-approximation becomes harder as the number of objectives grows and much harder for highly-rugged instances, whereas the trend w.r.t. objective correlation is less clear, more algorithm-dependent.

In the next section, we list the problem features that are intuitively impacting the performance of dominance-based search heuristics for the class of ρMNK-landscapes, and we explicitly assess their separate and joint effect on the runtime of PLS and GSEMO.

FEATURE-BASED ANALYSIS

Characterizing Problem Difficulty

Besides the benchmark parameters defining ρMNK-landscapes, a total of twelve general-purpose problem features are summarized in Table 1; they can be classified into the following three groups:

(i) Features related to the Pareto set (PS) includes the number of non-dominated solutions, the number of supported solutions [START_REF] Ehrgott | Multicriteria optimization[END_REF], and the hypervolume of the Pareto front, setting the reference point to the origin.

(ii) Features related to the distance and the connectedness of Pareto optimal solutions includes the average and the maximum distance between any pair of solutions from the PS, the number of clusters relative to the PS size, the proportional size of the largest cluster, and the minimum distance to consider such that all solutions belong to the same cluster. The distance measure is here taken as the Hamming distance between binary strings, which is directly related to the bit-flip neighborhood operator.

(iii) Features related to multimodality or ruggedness of the landscape includes the number of Pareto local optimal (PLO) solutions, the length of a single-solution Pareto hill-climbing that iteratively accepts a dominating neighbor until a PLO is found, as well as the autocorrelation of hypervolume of single solutions or of their neighborhood sets (local hypervolume), along a random walk. The random walk length is set to = 10 4 , and the neighborhood is the 1-bit-flip.

Obviously, since some of those features require the solution space to be enumerated exhaustively, they are not practical for a performance prediction purpose. However, we decided to include them npo Number of solutions from the PS [1, 13] supp Proportion of supported solutions in the PS [START_REF] Knowles | Instance generators and test suites for the multiobjective quadratic assignment problem[END_REF] hv Hypervolume [START_REF] Zitzler | Performance assessment of multiobjective optimizers: An analysis and review[END_REF] of the PS [1] avgd Avg. dist. between solutions from the PS [START_REF] Liefooghe | What makes an instance difficult for black-box 0-1 evolutionary multiobjective optimizers?[END_REF] maxd Max. dist. between solutions from the PS [START_REF] Knowles | Instance generators and test suites for the multiobjective quadratic assignment problem[END_REF] ncomp Relative number of clusters in the PS [START_REF] Paquete | Clusters of non-dominated solutions in multiobjective combinatorial optimization: An experimental analysis[END_REF] lcomp Proportional size of the PS largest cluster [START_REF] Liefooghe | On local search for bi-objective knapsack problems[END_REF] dconn Minimal distance to connect the PS [START_REF] Paquete | Clusters of non-dominated solutions in multiobjective combinatorial optimization: An experimental analysis[END_REF] nplo Number of Pareto local optima [START_REF] Paquete | On local optima in multiobjective combinatorial optimization problems[END_REF] ladapt Length of adaptive walk [START_REF] Verel | On the structure of multiobjective combinatorial search space: MNK-landscapes with correlated objectives[END_REF] corhv Autocorrelation of solution hypervolume [START_REF] Liefooghe | What makes an instance difficult for black-box 0-1 evolutionary multiobjective optimizers?[END_REF] corlhv Autocorrelation of local hypervolume [START_REF] Liefooghe | What makes an instance difficult for black-box 0-1 evolutionary multiobjective optimizers?[END_REF] in order to examine their impact on the algorithm performance. For a more comprehensive explanation of those features and their intercorrelation, we refer to [START_REF] Liefooghe | What makes an instance difficult for black-box 0-1 evolutionary multiobjective optimizers?[END_REF]. In the next section, we relate the value of those features for enumerable ρMNK-landscapes to the performance of both GSEMO and PLS.

Correlation Analysis

A first assessment of the dependency of the search performance on instance features can be done through visual inspection of scatter plots, supported by a correlation analysis. Naturally, correlation does not imply causation, and we do not draw any direct link between each considered feature and the algorithm runtime, even if in our case the eventual link could only go in one direction. Instead, we restrict ourself to measure the association of each feature to the performance metric (ert). We quantify this dependency strength via the Kendall's τ statistic [START_REF] Mcleod | Kendall: Kendall rank correlation and Mann-Kendall trend test[END_REF]. This rank-based nonlinear correlation measure is based on the ordering of all possible pairs, and its value is proportional to the difference between the number of concordant and discordant pairs among all possible pairwise comparisons. As such, when the null hypothesis of mutual independence (H0 : τ = 0) is rejected, τ can be directly interpreted as the prob- The scatter plots and the regressions (with local polynomial fitting) of the average running time of both algorithms as a function of the instance features are provided in Figure 3. In addition, Kendall's τ coefficients are given by the red points in Figure 4. For both algorithms, the average distance between Pareto optimal solution (avgd) is highly positively correlated with ert: the larger this distance, the longer the running time. On the contrary, both ruggedness-related features based on measures of hypervolume autocorrelation (corhv, corlhv), and one feature related to the connectedness, i.e. the size of the largest cluster in the PS (lcomp), are highly negatively correlated. As expected, ruggedness and connectedness play a major role for both algorithms: the runtime decreases with corhv and corlhv, and when a large number of Pareto optimal solutions are clustered in the solution space. Some features have a different impact on the two algorithms, possibly highlighting their respective strengths and weaknesses. In particular, the runtime of PLS increases with the number of Pareto optimal solutions and with the number of PLO solutions. Contrastingly, the scatter plots show that having a high number of PLO has less impact on GSEMO than PLS. Moreover, the runtime of GSEMO is correlated with three others features related to the connectedness of the PS (maxd, ncomp, dconn). Indeed, the relations between Pareto optimal solutions have a large effect on the runtime of GSEMO, especially when the distance between those solutions is large. Surprisingly, the runtime of PLS does not increase when non-dominated solutions are disconnected.

However, we have to be careful when drawing conclusions by aggregating data from different areas of the instance parameters space, since feature values and their range depend, in turn, on the levels of ρ, M , and K. This can be visually appreciated in Figure 3: the autocorrelation measures corhv and corlhv, for example, are clearly clustered around five levels that actually correspond to the different K-values. Similarly, we are able to distinguish three clusters in the hypervolume metric hv, which actually follow the objective space dimension M . Therefore, we deepen the analysis by evaluating the correlation within the instance groups defined by each possible combination of the ρ, M, K -values under study. Black points and lines in Figure 4 display the average τ -value within groups, together with the confidence interval associated to the mean. By comparing them with red points, we can clearly notice how data aggregation slightly enhances the correlation statistic in the corhv and corlhv case, leading to the same inference nonetheless. On the contrary, although hypervolume is very weakly associated with runtime overall, and that its impact is contradictory between GSEMO and PLS, group results are more consistent, showing a strong positive association between ert and hv for both algorithms. Unfortunately, as for features related to the connectedness of the PS, our confidence on the average correlation within groups is too low to make further comments, mainly due to the fact that, in many cases, we could not reject the null hypothesis of mutual independency at the group level. Nevertheless, the previous observations on data blocks and their possible effect motivate the remainder of our analysis.

Regression Analysis

In this section we aim to quantify the impact of instance features, and possibly disentangle their individual contribution to the performance variance, by taking into account the dependency among Points give Kendall's τ statistic for the correlation between runtime and instance feature, evaluated on the whole set of instances (red) or within instance groups (black, see left legend). Group averages, 0.95 confidence intervals, and significance, are estimated with a t-test considering only statistics that were already significant at the group level.

Group and overall significance are based on Kendall's test p-value at the 0.05 level (H0 : τi = 0) [START_REF] Mcleod | Kendall: Kendall rank correlation and Mann-Kendall trend test[END_REF].

measurements that is induced by the experimental plan. Our goal is precisely to generalize from it as much as possible, in order to make inferences about the effect and relative importance of features.

Linear Mixed-Model

Let yi be the log-transformed expected runtime (ert) on the ith problem instance. We treat yi as an observation from a random variable Y with expectation E(Y ) = µ, i.e. yi = µ + εi where εi are taken to be independent, identically distributed, and zero-mean.

In a classical multilinear regression, we would model µ as a linear combination of p predictors, notably (a subset of) the p problem instance features that we can measure. Thus, the performance observation on instance i can be written as:

yi = β0 + p k=1 β k x ki + εi εi ∼ N (0, σ 2 )
where εi is the usual random term, i.e. the regression residual. In this model, performance observations are supposed to be i.i.d. from a normal distribution yi ∼ N (µ, σ 2 ). However, as discussed in the previous sections, our observations are mostly clustered around the different combinations of benchmark parameters; see Figure 1.

In fact, a simple linear regression on a dummy categorical predictor having a different level for each combination of ρ, M , and K, would explain 84.51% and 86.85% of the ert variance of GSEMO and PLS, respectively. Since we rather want to investigate the impact of instance features, we need to decompose that global performance variance into what is due to the grouping of benchmark parameters, from which we would like to generalize, and what is due to the randomness involved in the instance generation process; namely for ρMNKlandscapes, the epistatic interaction links and their contributions to the objective values. That conveys the feature variance within the blocks of our experimental design. To this end, instead of fitting an independent regression model for each instance group, we build a linear mixed-model with random effects for experimental blocks. In such framework, the performance on instance i from group j can then be modelled as:

yij = β0 + p k=1 β k x kij + αj + εij εij ∼ N (0, σ 2 )
where αj are i.i.d random variables, with αj ∼ N (0, σ 2 α ) denoting the group effect. By doing so, we suppose that features have the same impact across groups; in linear terms, constant (fixed) slopes and random intercepts. Notice that estimates β and residuals ε will likely not be the same as in the previous model. Notice also that performance observations are now i.i.d. conditionally on the grouping factor yij|αj ∼ N (µ+αj, σ 2 ), whereas the unconditional model yij ∼ N (µ, σ 2 + σ 2 α ) carries a dependency between measurements on the same group and yields the aforementioned variance decomposition [START_REF] Chiarandini | Mixed models for the analysis of optimization algorithms[END_REF].

The usual approach to estimate the parameters of the unconditional model is the restricted maximum likelihood method; see [START_REF] Bates | Fitting linear mixed-effects models using lme4[END_REF][START_REF] Faraway | Extending the linear model with R: generalized linear, mixed effects and nonparametric regression models[END_REF] for theoretical and implementation details. In the following, we present the results of such estimation on the full model comprising all instance features and for both algorithms. In particular, each estimated βi is tested against the null hypothesis H0 : βi = 0, whereas as for the group effect we only need to check H0 : σ 2 α = 0. In the multilinear framework, the regression coefficients that are statistically significant allow us to assess the runtime effect of each feature conditionally on all the others and, given the random effect formulation, across experimental blocks.

Finally, in order to assess the accuracy of a regression model, the conventional R 2 (ratio of variance explained by the regression) can be extended by taking into account the variance decomposition that is specific to mixed-models [START_REF] Nakagawa | A general and simple method for obtaining r2 from generalized linear mixed-effects models[END_REF]. We obtain a marginal R 2 yielding the proportion of variance explained by instance features, and a conditional R 2 which, despite its name, gives the proportion of variance explained by the entire model, i.e. including the random effect of benchmark parameter combinations. Marginal and conditional R 2 are respectively 0.617 and 0.919 for the regression modelling GSEMO's ert, respectively 0.482 and 0.911 for PLS.

Predictors Effect Size

We report in Figure 5 the values of the regression coefficients estimated from the mixed-model. In a multilinear regression, each coefficient βi predicts the change in the conditional mean of the response variable after a unitary change of the corresponding covariate xi. As such, higher values correspond, visually, to steeper slopes of the partial regression lines where the given predictor is the only covariate and all other predictors remain constant. This simple interpretation motivates the choice of a multilinear model.

For both algorithms, the features having the highest individual impact on the runtime are the hypervolume (hv) and its autocorrelation measures (corhv and corlhv, in order of importance). For GSEMO, the number of Pareto optimal solutions has a significant effect on estimated runtime, the larger npo the shorter ert, whereas we cannot reject the hypothesis that the number of PLO solutions has no effect on GSEMO's ert. Still, the effect size of the length of an adaptive walk (ladapt), which is a good estimator of nplo [START_REF] Verel | On the structure of multiobjective combinatorial search space: MNK-landscapes with correlated objectives[END_REF], is significant for GSEMO. Conversely, the opposite is true for PLS. Once again, the relative size of the largest connected component of the PS (lcomp) impacts the performance of both algorithm. However, when controlling for all other features (i.e. conditionally on all other predictors), we find that an increase in the PS connectedness yields an increase in the expected runtime. Surprisingly also, the number of supported solutions (supp) is a significant predictor for both algorithms, even if none of them explicitly exploits this feature during the search process. Finally, despite being one of the features with the highest correlation to ert, the average distance between non-dominated solutions (avgd) has little impact on GSEMO and no significant impact on PLS.

Relative Importance of Predictors

Variable importance is commonly assessed via feature selection. However, stepwise selection can be misleading: intercorrelated predictors have a confounding effect on each other, but what the multilinear regression tries to measure is precisely the effect of one variable controlling for the others (i.e. fixing all the others), which leads to a poor estimation in the presence of collinearity. In an information theoretical framework, the Akaike Information Crite-rion (AIC) [START_REF] Akaike | Information theory as an extension of the maximum likelihood principle[END_REF] measures the relative quality of a model on a given dataset, not in terms of accuracy, but in terms of likelihood, i.e. the relative distance to the unknown true mechanism generating the observed data. The difference in AIC-values between two models can then be used to estimate the strength of evidence for one model against the other. On a set of alternate models, AIC differences can be transformed into so-called Akaike weights, which can be directly interpreted as the probability for each model to be the "best" one, conditionally on the considered set of models [START_REF] Burnham | Model selection and multimodel inference: a practical information-theoretic approach[END_REF]. In this context, instead of performing feature selection, variable importance can be better assessed by making inference from all candidate models [START_REF] Burnham | Multimodel inference understanding aic and bic in model selection[END_REF]. We perform an exhaustive search in the space of all 2 p models we can build with our p predictors. The sum of Akaike weights for the 2 (p-1) models that contain a given predictor, can give us an estimation of the relative importance of that particular variable.

Results are reported in Figure 6: hv, corhv, and corlhv are strong explanatory features for both models, albeit in different order of relative importance. For GSEMO, npo, ladapt, and lcomp can also be considered as important features, whereas nplo is the second most-important predictor of PLS runtime.

CONCLUSIONS

In this paper we investigated the impact of problem characteristics on the performance of two simple dominance-based multiobjective search heuristics. More particularly, we analyzed the expected runtime of GSEMO and PLS to identify a (1 + ε)-approxi-mation of the Pareto set for enumerable ρMNK-landscapes. Beyond a more conventional performance assessment of competing algorithms, we studied the correlations between algorithm performance and a set of relevant problem features. We contrasted their association with algorithm runtime, within and across instance classes. By building a mixed-model regression taking into account the experimental design, we were able to make more general inferences about the impact of problem features, their relative importance, as well as similarities and differences between both algorithms in this respect. Notably, the global hypervolume indicator value and its local autocorrelation measures, which are related to the ruggedness of the fitness landscape, have the strongest impact for both algorithms. In addition, as expected, our analysis reveals that the runtime of PLS appears to be largely affected by the landscape multimodality, here defined as the number of Pareto local optimal solutions. Moreover, the Pareto set cardinality plays a more important role for GSEMO than for PLS. At last, the features reflecting the connectedness of the Pareto set are important predictors for the runtime of GSEMO, which navigates the search space with an ergodic variation operator. On the contrary, the same features have a nonsignificant effect on the performance of PLS, which is constrained by the exploration of a finite neighborhood.

A similar analysis on other problem classes and large-size instances, with computationally affordable features, would allow us to further understand and gain more knowledge about the main difficulties that EMO algorithms have to face, and eventually to predict their expected performance.

Figure 2 :

 2 Figure2: interaction plots between the average estimated runtime ert (y-axis, log-scale) and the benchmark parameters: objective correlation ρ (left), the number of objectives M (center), and the problem non-linearity K (right, see titles). Results are grouped by algorithm (see legend). The average (lines) and is 0.95 confidence interval (shaded areas) are evaluated through bootstrapping.

Figure 3 :

 3 Figure 3: interaction plots between the average estimated runtime ert (y-axis, fixed logscale) and the instance features (x-axis, see facets titles). Results are grouped by algorithm (see legend). Lines represent a locally-fitted polynomial function (LOESS).

  ability of observing agreement (τ > 0) or disagreement (τ < 0) between the ranks of paired values.

Figure 4 :

 4 Figure4: performance-feature association. Points give Kendall's τ statistic for the correlation between runtime and instance feature, evaluated on the whole set of instances (red) or within instance groups (black, see left legend). Group averages, 0.95 confidence intervals, and significance, are estimated with a t-test considering only statistics that were already significant at the group level. Group and overall significance are based on Kendall's test p-value at the 0.05 level (H0 : τi = 0)[START_REF] Mcleod | Kendall: Kendall rank correlation and Mann-Kendall trend test[END_REF].

Figure 5 :Figure 6 :

 56 Figure 5: conditional impact of instance features on log 10 expected runtime in a linear mixed-model. Points and bars give, respectively, the estimates and 0.95 confidence intervals of model parameters for intercept and fixed effects βi (H0 : βi = 0) [5].

  Create x by flipping each bit of x with a rate 1/N ;

	1 Choose an initial solution x0 uniformly from X;
	2 A ← {x0};
	3 repeat
	4	Select an element x out of A uniformly;
	6	A ← non-dominated solutions from A ∪ {x };
	7 until success ∨ maxeval;
		Algorithm 2: PLS
	1 Choose an initial solution x0 uniformly from X;
	2 A ← {x0};
	3 repeat
	4	Select a non-visited element x out of A uniformly;
	5	Create N(x) by flipping each bit of x in turns;
	6	Flag x as visited;
	7	A ← non-dominated solutions from A ∪ N(x);
	8 until all-visited ∨ success ∨ maxeval;

5

Table 1 :

 1 Summary of ρMNK-landscape benchmark parameters and problem instance features investigated in the paper.

	Benchmark parameters (3)
	K Number of variable (epistatic) interactions
	M Number of objective functions
	ρ Correlation between the objective values
	Problem instance features (12)

Acknowledgments

This work was partially supported by the Japanese-French JSPS project "Global Research on the Framework of Evolutionary Solution Search to Accelerate Innovation" (2013-2016), and by the Japanese-French JSPS-Inria project "Threefold Scalability in Anyobjective Black-Box Optimization" (2015-2017).


Keywords: Joueurs artificiels de Go, Exploration-exploitation, UCT, Monte-Carlo, Motifs

Modification d'UCT avec motifs dans le Monte-Carlo Go Résumé : L'algorithme UCB1 pour le problème du bandit-manchot a récemment été étendu en l'algorithme UCT (Upper bound Confidence for Tree) pour la recherche arborescente min-max. Nous avons développé un joueur artificiel de Go, MoGo, basé sur des simulations Monte-Carlo, qui est le premier programme de Go utilisant UCT. Nous exposons notre modification de l'algorithme UCT pour l'application au jeu de Go, ainsi que l'utilisation de motifs dans les simulations aléatoires qui ont permis d'augmenter significativement le niveau de MoGo. Nous introduisons d'autre part des techniques d'élagage dans l'algorithme UCT pour les grands Goban, ainsi que de la parallélisation d'UCT. MoGo est maintenant un joueur artificiel de premier plan sur les Gobans de taille 9 × 9 et 13 × 13.

Introduction

The history of Go stretches back some 4000 years and the game still enjoys a great popularity all over the world. Although its rules are simple (see http://www.gobase.org for a comprehensive introduction), its complexity has defeated the many attempts done to build a good Computer-Go player since the late 70's [START_REF] Bouzy | Computer go: An AI oriented survey[END_REF]. Presently, the best Computer-Go players are at the level of weak amateurs; Go is now considered one of the most difficult challenges for AI, replacing Chess in this role.

Go differs from Chess in many respects. First of all, the size and branching factor of the tree are significantly larger. Typically the Go board ranges from 9 × 9 to 19 × 19 (against 8 × 8 for the Chess board); the number of potential moves is a few hundred against a few dozen for Chess. Secondly, no efficient evaluation function approximating the minimax value of a position is available. For these reasons, the powerful alpha-beta search used by Computer-Chess players (see [START_REF] Newborn | Computer Chess Comes of Age[END_REF]) failed to provide good enough Go strategies.

Recent progress has been done regarding the evaluation of Go positions, based on Monte-Carlo approaches [START_REF] Bruegmann | Monte carlo go[END_REF] (more on this in section 2). However, this evaluation procedure has a limited precision; playing the move with highest score in each position does not end up in winning the game. Rather, it allows one to restrict the number of relevant candidate moves in each step. Still, the size of the (discrete) search space makes it hardly tractable to use some standard Reinforcement Learning approach [START_REF] Sutton | Reinforcement learning: An introduction[END_REF], to enforce the exploration versus exploitation (EvE) search strategy required for a good Go player.

Another EvE setting originated from Game Theory, the multi-armed bandit problem, is thus considered in this paper. The multi-armed bandid problem models the gambler, choosing the next machine to play based on her past selections and rewards, in order to maximize the total reward [START_REF] Auer | Gambling in a rigged casino: the adversarial multi-armed bandit problem[END_REF]. The UCB1 algorithm proposed by Auer et al. in the multiarmed bandit framework [START_REF] Auer | Finite-time analysis of the multiarmed bandit problem[END_REF] was recently extended to tree-structured search space by Kocsis et al. (algorithm UCT) [START_REF] Kocsis | Bandit based monte-carlo planning[END_REF].

The main contributions of the player we present (named MoGo) are: (i) modification of UCT algorithm for Go, (ii) original use of sequence-like simulations in Monte-Carlo evaluation function. Several algorithmic (dynamic tree structure [START_REF] Coulom | Efficient selectivity and backup operators in monte-carlo tree search[END_REF], parallelized implementation) or heuristic (simple pruning heuristics) issues were also tackled. MoGo has reached a comparatively good Go level: MoGo has been ranked as the first Go program out of 142 on 9 × 9 Computer Go Server (CGOS1 ) since August 2006; and it won all the tournaments (9x9 and 13x13) on the international Kiseido Go Server2 on October and November 2006.

This paper is organized as follows. Section 2 briefly introduces related work, assuming the reader's familiarity with the basics of Go. Section 3 describes MoGo, focussing on our contributions: the implementation of UCT in large sized search spaces, and the use of prior, pattern-based, knowledge to bias the Monte-Carlo evaluation. Experiment results are reported and discussed in Section 4. The paper concludes with some knowledge-and computer-intensive perspectives for improving MoGo.

Previous Related Work

Our approach is based on the Monte-Carlo Go and multi-armed bandit problems, which we present respectively in Section 2.1 and 2.2. UCT, which applies multi-armed bandit techniques to minimax tree search, is presented in Section 2.3. We suppose minimax tree and alpha-beta search is well known for the reader.

Monte-Carlo Go

Monte-Carlo Go, first appeared in 1993 [START_REF] Bruegmann | Monte carlo go[END_REF], has attracted more and more attention in the last years. Monte-Carlo Go has been surprisingly efficient, especially on 9 × 9 game; CrazyStone, developed by Rémi Coulom [START_REF] Coulom | Efficient selectivity and backup operators in monte-carlo tree search[END_REF], a program using stochastic simulations with very little knowledge of Go, is the best known 3 .

Two principle methods in Monte-Carlo Go are also used in our program. First we evaluate Go board situations by simulating random games until the end of game, where the score could be calculated easily and precisely. Second we combine the Monte-Carlo evaluation with minimax tree search. We use the tree structure of CrazyStone [START_REF] Coulom | Efficient selectivity and backup operators in monte-carlo tree search[END_REF] in our program.

Remark 1

We speak of a tree, in fact what we have is often an oriented graph. However, the terminology "tree" is widely used. As to the Graph History Interaction Problem (GHI) explained in [START_REF] Kishimoto | A general solution to the graph history interaction problem[END_REF], we ignore this problem considering it not very serious, especially compared to other difficulties in Computer-Go.

Bandit Problem

A K-armed bandit, is a simple machine learning problem based on an analogy with a traditional slot machine (one-armed bandit) but with more than one arm. When played, each arm provides a reward drawn from a distribution associated to that specific arm. The objective of the gambler is to maximize the collected reward sum through iterative plays 4 . It is classically assumed that the gambler has no initial knowledge about the arms, but through repeated trials, he can focus on the most rewarding arms.

The questions that arise in bandit problems are related to the problem of balancing reward maximization based on the knowledge already acquired and attempting new actions to further increase knowledge, which is known as the exploitation-exploration dilemma in reinforcement learning. Precisely, exploitation in bandit problems refers to select the current best arm according to the collected knowledge, while exploration refers to select the suboptimal arms in order to gain more knowledge about them.

A K-armed bandit problem is defined by random variables X i,n for 1 ≤ i ≤ K and n ≥ 1, where each i is the index of a gambling machine (i.e., the "arm" of a bandit). Successive plays of machine i yield rewards X i,1 ,X i,2 ,... which are independent and identically distributed according to a certain but unknown law with unknown expectation µ i . Here independence holds also for rewards across machines; i.e., X i,s and X j,t are independent (probably not identically distributed) for each 1 ≤ i < j ≤ K and each s, t ≥ 1. Algorithms choose the next machine to play depending on the obtained results of the previous plays. Let T i (n) be the number of times machine i has been played after the first n plays. Since the algorithm does not always make the best choice, its expected loss is studied. Then the regret after n plays is defined by

µ * n - K j=1 µ j E[T j (n)] where µ * = max 1≤i≤K µ i
E[ ] denotes expectation. In the work of Auer and Al. [START_REF] Auer | Finite-time analysis of the multiarmed bandit problem[END_REF], a simple algorithm UCB1 is given, which ensures the optimal machine is played exponentially more often than any other machine uniformly when the rewards are in

[0, 1]. Note Xi,s = 1 s s j=1 X i,j , Xi = Xi,Ti(n) ,
then we have:

Algorithm 1 Deterministic policy: UCB1

• Initialization: Play each machine once.

• Loop: Play machine j that maximizes Xj + 2 log n Tj (n) , where n is the overall number of plays done so far.

One formula with better experimental results is suggested in [START_REF] Auer | Finite-time analysis of the multiarmed bandit problem[END_REF]. Let

V j (s) = 1 s s γ=1 X 2 j,γ -X2 j,s + 2 log n s
be an estimated upper bound on the variance of machine j, then we have a new value to maximize:

Xj + log n T j (n) min{1/4, V j (T j (n))} . (1)
According to Auer and Al., the policy maximizing (1) named UCB1-TUNED, considering also the variance of the empirical value of each arms, performs substantially better than UCB1 in all his experiments. This corresponds to our early results and then we use always the policy UCB1-TUNED in our program5 .

UCT: UCB1 for Tree Search

UCT [START_REF] Kocsis | Bandit based monte-carlo planning[END_REF] is the extension of UCB1 [START_REF] Auer | Finite-time analysis of the multiarmed bandit problem[END_REF] to minimax tree search. The idea is to consider each node as an independent bandit, with its child-nodes as independent arms. Instead of dealing with each node once iteratively, it plays sequences of bandits within limited time, each beginning from the root and ending at one leaf. The algorithm UCT is defined in Table 1 6 . The program continues playing one sequence each time, which is defined from line 1 to line 8. Line 9 to line 21 are the function using UCB1 for choosing one arm (one child-node in the UCT case). Line 15 ensures each arm be selected once before further exploration. Line 16 applies the formula of UCB1. After each sequence, the value of played arm of each bandit is updated7 iteratively from the father-node of the leaf to the root by formula UCB1, described in function updateV alue from line 22 to line 29. Here the code deals with the minimax case. In general, the value of each node converges to the real max (min) value as the number of simulations increases. In the problems of minimax tree search, what we are looking for is often the optimal branch at the root node. It is sometimes acceptable if one branch with a score near to the optimal one is found, especially when the depth of the tree is very large and the branching factor is big, like in Go, as it is often too difficult to find the optimal branch within short time.

In this sense, UCT outperforms alpha-beta search. Indeed we can outlight three major advantages. First, it works in an anytime manner. We can stop at any moment the algorithm, and its performance can be somehow good. This is not the case of alpha-beta search. Figure 2 shows if we stop alpha-beta algorithm prematurely, some moves at first level has even not been explored. So the chosen move can be far from optimal. Of course iterative deepening can be used, and solve partially this problem. Still, the anytime property is stronger for UCT and it is easier to finely control time in UCT algorithm. Second, UCT is robust as it automatically handles uncertainty in a smooth way. At each node, the computed value is the mean of the value for each child weighted by the frequency of visits. Then the value is a smoothed estimation of max, as the frequency of visits depends on the difference between the estimated values and the confidence of this estimates. Then, if one child-node has a much higher value than the others, and the estimate is good, this child-node will be explored much more often than the others, and then UCT selects most of the time the 'max' child node. However, if two child-nodes have a similar value, or a low confidence, then the value will be closer to an average. Third, the tree grows in an asymmetric manner. It explores more deeply the good moves. What is more, this is achieved in an automatic manner. Figure 1 gives an example.

Figure 1 and Figure 2 compares clearly the explored tree of two algorithms within limited time. However, the theoretical analysis of UCT is in progress [START_REF] Kocsis | Improved monte-carlo search[END_REF]. We just give some remarks on this aspect at the end of this section. It is obvious that the random variables involved in UCT are not identically distributed nor independent. This complicates the analysis of convergence. In fact we can define the bias for the arm i by:

δ i,t = µ * i - 1 t t s=1 X i,s ,
where µ * i is the minimax value of this arm. It is clear that at leaf level δ i,t = 0. We can also prove that

δ i,t ≤ K D log t t ,
with K constant and D the depth of the arm (counted from the root down). This corresponds to the fact that the bias is amplified when passing from deep level to the root, which prevents the algorithm from finding quickly the optimal arm at the root node. An advantage of UCT is that it adapts automatically to the 'real' depth. For each branch of the root, its 'real' depth is the depth from where δ i,t = 0 holds true. For these branches, the bias at the root is bounded by K d log t t with the real depth d < D. The values of these branches converging faster than the other, UCT spends more time on other interesting branches.

Main Work

In this section we present our program MoGo using UCT algorithm. Section Carlo program: the quality of simulations (then the estimation of score) and the depth of the tree, we show in the two following sections our corresponding improvements. Section 3.2 presents the sequence-like random simulation with patterns. Section 3.3 presents ideas for tree search pruning on large Go board. Section 3.4 presents the modification on the exploring order of non-visited nodes. At last, Section 3.5 presents parallelization.

Application of UCT for Computer-Go

MoGo contains mainly two parts, namely the tree search part and the random simulation part, as shown in Figure 3. Each node of the tree represents a Go board situation, with child-nodes representing next situations after corresponding move.

The application of UCT for Computer-Go is based on the hypothesis that each Go board situation is a bandit problem, where each legal move is an arm with unknown reward but of a certain distribution. We suppose that there are only two kinds of arms, the winning ones and the losing ones. We set respectively reward 1 and 0. We ignore the case of draw, which is too rare in Go.

In the tree search part, we use a parsimonious version of UCT by introducing the same dynamic tree structure as in CrazyStone [START_REF] Coulom | Efficient selectivity and backup operators in monte-carlo tree search[END_REF] in order to economize memory. The tree is then created incrementally by adding one node after each simulation as explained in the following. This is different from the one presented in [START_REF] Kocsis | Bandit based monte-carlo planning[END_REF], and is more efficient because less nodes are created during simulations. In other words, only nodes visited more than twice are saved, which economizes largely the memory and accelerates the simulations. The pseudocode is given in Table 2. Again we do not talk about optimization.

During each simulation, MoGo starts from the root of the tree that it saves in the memory. At each node, MoGo selects one move according to the UCB1 formula 1. MoGo then descends to the selected child node and selects a new move (still according to UCB1) until such a node has not yet been created in the tree. This part corresponds to the code from line 1 to line 5. The tree search part ends by creating this new node (in fact one leaf) in the tree. This is finished by createN ode. Then MoGo calls the random simulation part, the corresponding function getV alueByM C at line 7, to give a score of the Go board at this leaf.

In the random simulation part, one random game is played from the corresponding Go board till the end, where score is calculated quickly and precisely according to the rules. The nodes visited during this random simulation are not saved. The random simulation done, the score received, MoGo updates the value at each node passed by the sequence of moves of this simulation8 .

Improving simulation with domain knowledge

In this section we introduce our sequence-like random simulation with patterns. Its advantage is obvious compared with the classical random simulation, which we call pure random mode. In the following we talk about our improved random mode as well as our implementation of patterns.

In the random simulation part, it is very important to have clever simulations giving credible scores. Using some simple 3 × 3 patterns inspired by Indigo [START_REF] Bouzy | Associating domain-dependent knowledge and monte carlo approaches within a go program[END_REF] (similar patterns can also be found in [START_REF] Ralaivola | Svm and pattern-enriched common fate graphs for the game of go[END_REF]), our random simulation is likely to have more meaningful sequences in random simulations than before, which has improved significantly the level of MoGo.

The essential difference on the usage of patterns between the ours and other patternbased Monte-Carlo programs ( [START_REF] Bouzy | Associating domain-dependent knowledge and monte carlo approaches within a go program[END_REF]) is that we use patterns to create meaningful sequences in simulations by finding local answers. The moves played are not necessarily globally better moves. It is not obvious that is more important to get better sequences rather that better moves to make the Monte-Carlo evaluation more accurate. However our experiments showed that the main improvement came from the use of local answers. If the same patterns are used to find interesting moves everywhere in the board instead of near the previous moves, the accuracy decreases. We believe that this claim is not obvious, and one of the main contribution of MoGo. We also don't use patterns for pruning in the tree. We have not investigated the more sophisticated patterns equipped by other programs like GnuGo.

In our pure random mode, legal moves are played on the Go board uniformly randomly, with few rules preventing the program from filling its own eyes. We also privilege the moves capturing some stones. On CGOS the bot named MoGo using exactly this mode has achieved rank score 1647 ELO 9 . Currently the rank of MoGo is close to 2200 ELO.

Then, since we were not satisfied by the pure random simulations which gave meaningless games most of the time, local patterns are introduced in order to have some more reasonable moves during random simulations. Our patterns are defined as 3 × 3 intersections, centered on one free intersection, where one move is supposed to be played. Our patterns consist of several functions, testing if one move in such a local situation (3 × 3) is interesting. More precisely, we test if one move satisfies some classical forms in Go games, for example cut move, Hane move, etc.

Moreover, we look for interesting moves only around the last played move on the Go board. This is because that local interesting moves look more likely to be the answer moves of the last moves, and thus local sequence appears when several local interesting moves are tested and then played continuously in random simulations.

We describe briefly how the improved random mode generates moves. It first verifies whether the last played move is an Atari; if yes, and if the stones under Atari can be saved (in the sense that it can be saved by capturing stones or increasing liberties), it chooses one saving move randomly; otherwise it looks for interesting moves in the 8 positions around the last played move and plays one randomly if there is any; otherwise it looks for the moves capturing stones on the Go board, plays one if there is any. At last, if still no move is found, it plays one move randomly on the Go board. Surely, the code of MoGo is actually complicated in details, with many small functions equipped with hand-coded Go knowledges. However, we believe the main frame given here is the most essential to have sequence-like simulations.

Figure 4 shows the first 30 moves of two random games using different modes. Moves generated by the improved random mode are obviously much more meaningful. We now give the detailed information on our patterns. The patterns are 3×3 intersections centered on an empty position, say p, where is supposed to play the next move. Each pattern is a boolean function, answering the question whether the next move playing on p is an interesting move. True is returned (when pattern is matched), if and only if the state of each position on the Go board is the same as the one on the corresponding position of the pattern, or there is a cross on the corresponding position (which means the situation of this position is ignored). Normally there is no constraints on the color of the next move (one move good for black is supposed to be also good for white). Some special cases are explained when mentioned. The symmetry, rotations and exchange of stone colors of patterns are considered. Moves are tested by patterns only if they are neither illegal moves nor self-Atari moves.

We have tried several patterns during the development of MoGo and implemented finally the ones shown in Figure 5, 6, 7 and8, where the position with a square is where the next move is supposed to be played. We used hand-coded patterns in our implementation. However, it will be more interesting if this can be achieved by a learning system. Another approach using Bayesian generation can be found in Bouzy's work [START_REF] Bouzy | Bayesian generation and integration of k-nearest-neighbor patterns for 19x19[END_REF].

Remark 3

We believe that it is not always better to have more 'good' patterns in the random modes, meanwhile what is more important is whether the random simulation can have some meaningful often. This claim needs more experiments.

The Table 3 shows clearly how patterns improve the overall performance.

INRIA

Random mode

Win

UCT with pruning ideas

In this section we show our ideas (heuristics) to reduce the huge tree size, which makes MoGo relatively strong on large Go board. Thus we gain a larger local depth in the search tree by losing the global view of UCT. Obviously pruning heuristics may lead to a sub-optimal solution. First we define group by Go knowledge to reduce the branching factor in tree search. Then zone division is derived from group, which helps to have a more precise score. We use group and zone mode for 13 × 13 and 19 × 19 Go board. Figure 9 and 10 will give two convincing examples.

Remark 4 As we are not very experienced for Go-knowledge-based programming and we had little time working on it, we believe other programs like GnuGo and AyaGo, or Monte-Carlo programs as [START_REF] Cazenave | Abstract proof search[END_REF][8] have more clever pruning techniques. Then, due to the space limitation, we do not give the detailed pseudo code of this part. However, our experimental results of combining UCT with pruning techniques are already encouraging.

First we define one group as a set of strings and free intersections on a Go board according to certain Go knowledge, which gathers for example one big living group and its close enemies. We have implemented Common Fate Graph (CFG) [START_REF] Graepel | Learning on graphs in the game of go[END_REF] in our program to help the calculation of groups. The method starts from one string and recursively adds close empty intersections and strings close to these empty intersections until no more close strings are found within a distance controled by a parameter.

INRIA

In group mode, in the tree search part we search only the moves in the group instead of all over the Go board. In random simulation part there is no more such restriction. Using groups, we reduce the branching factor to less than 50 at the opening period. Then, depth of MoGo's tree could be around 7-8 on large Go board. Table 4 shows MoGo becomes competitive on large Go board by using group pruning technique. However, sophisticated pruning techniques are undoubtedly necessary to improve the level of Computer-Go programs.

Opponents

Win Then we developed the zone mode. As explained above, group mode limits the selection of moves in the tree search part. It has however no restriction on the random simulation. In zone mode, the random moves is generated only in a certain zone instead of on the whole Go board. This is because on the large Go board, random simulations are too long (of an average length more than 400 for 19 × 19 Go board) to be able to return credible scores. The score obtained in this zone is then tuned for the global Go board.

The tuning idea is very crucial in our zone mode, otherwise the global evaluation of situation by random simulations, which is very imprecise most of the time for the large Go board, biases too much the local decision. More precisely, before UCT we estimate, by a certain number (5000 in our experiments) of Monte-Carlo simulations, the score s0 for the zone decided by the current situation (including last move of opponent), and the score s1 for the area outside the zone. Then we tune the s 1 by

ŝ1 = α(s 1 + s0 -komi) -(s 0 -komi) (2)
and then after each simulation in the zone, with local score s 0 we return s 0 + ŝ1 as the final score of the simulation. In this way, the score is always closer to a close game score than it should have been by using s 0 + s1 . This helps MoGo play more normally locally.

Modification of exploring order for non-visited nodes

UCT works very well when the node is frequently visited as the trade-off between exploration and exploitation is well handled by UCB1 formula. However, for the nodes far from the root, whose number of simulations is very small, UCT tends to be too much exploratory. This is due to the fact that all the possible moves in one position are supposed to be explored before using the UCB1 formula. Thus, the values associated to moves in deep nodes are not meaningful, since the child-nodes of these nodes are not all explored yet and, sometimes even worse, the visited ones are selected in fixed order. This can leads to bad predicted sequences. Two modifications are made to have a better order.

First-play urgency

UCB1 algorithm begins by exploring each arm once, before using the formula (1). This can sometimes be unefficient especially if the number of trials is not large comparing to the number of arms. This is the case for numerous nodes in the tree (number of visits is small comparing to the number of moves). For example if an arm keeps returning 1 (win), there is no good reason to explore other arms. We have set a fixed constant named firstplay urgency (FPU) in the algorithm. For each move, we name its urgency by the value of formula [START_REF] Auer | Finite-time analysis of the multiarmed bandit problem[END_REF]. The FPU is by default set to ∞ for each legal move before first visit (see line 15 in Table 1). Any node, after being visited at least once, has its urgency updated according to UCB1 formula. We play the move with the highest urgency. Thus, the FPU ∞ ensures the exploration of each move once before further exploitation of any previously visited move.

On the other way, smaller FPU ensures earlier exploitations if the first simulations lead to an urgency smaller than FPU (in this case the other unvisited nodes are not selected). This improved the level of MoGo according to our experiment as shown in Table 7.

INRIA

Use information of the parents

One assumption that can be made in go game is that given a situation, good moves may sometimes still be good ones on the following move. When we encounter a new situation, instead of exploring each move m in any order, we can use the value estimation of m in an earlier position to choose a better order. We typically use the estimated values of the grandfather of the node. We believe this helps MoGo on the large Go board, however we do not have enough experiments to claim significant results.

Parallelization

As UCT scales well with time, we made MoGo run on a multi-processors machine with shared memory. The modifications to the algorithm are quite straightforward. All the processors share the same tree, and the access to the tree is locked by mutexes. As UCT is deterministic, all the threads could take exactly the same path in the tree, except for the leaf. The behavior of the multithreaded UCT as presented here is then different from the monothreaded UCT. Two experiments has then to be done. First, we can compare the level of MoGo using the monothreaded or the multithreaded algorithms while allowing the same number of simulations per move. All such experiments showed non significant differents in the play level 10 . Second, we can compare the level using the same time per move (the multithreaded version will then make more simulations per move). As UCT benefits from the computational power increase, the multithreaded UCT is efficient (+100 ELO on CGOS with 4 processors).

Results

We list in this section several experiment results who reflect characteristics of the algorithm.

All the tests are made by letting MoGo play against GnuGo 3.6 with default mode. Komi are set to 7.5 points. In the tables, are given, the winning rate when MoGo plays black and white are given along the number of games played in each color (between parenthesis). The confidence interval is 95%.

Dependence of Time

The performance of our program depends on the given time (equally the number of simulations) for each move.

Parametrized UCT

We parametrize the UCT implemented in our program by two new parameters, namely p and F P U . First we add one coefficient p to formula UCB1-TUNED (1), which by default is 1. This leads to the following formula: choose j that maximizes:

Xj + p log n T j (n) min{1/4, V j (n j )}
p decides the balance between exploration and exploitation. To be precise, the smaller the p is, the deeper the tree is explored. According to our experiment shown in 33% ± 9.4% 0.80 33% ± 9% (100) 39% ± 10% (50) 36% ± 6.7% 1.0 40% ± 8% (150) 38% ± 8% (150) 39% ± 5.6% 1.1 39% ± 8% (150) 41% ± 8% (150) 40% ± 5.6% 1.2 40% ± 8% (150) 44% ± 8% (150) 42% ± 5.7% 1.5 30% ± 13% (50) 26% ± 12% (50) 28% ± 9% 3.0 36% ± 13% (50) 24% ± 12% (50) 30% ± 9% 6.0 22% ± 11% (50) 18% ± 10% (50) 20% ± 8% The second is the first-play urgency (FPU) as explained in 3.4. Some results are shown in Table 7. We believe that changing exploring order of non-visited nodes can bring further improvement. We finally use 1.1 as the initiate FPU for MoGo on CGOS.

Results On CGOS

MoGo is ranked as the first program on 9 × 9 Computer Go Server since August.

Conclusion

The success of MoGo shows the efficiency of UCT compared to alpha-beta search in the sense that nodes are automatically studied with better order, especially in the case when search time is too limited. We have discussed the advantages of UCT relevant to Computer-Go. It is worthy to mention that a growing number of top level Go programs now use UCT.

FPU

Winning Rate Winning rate Total for Black Games for White Games Winning Rate 1.4 37% ± 9% (100) 38% ± 10% (100) 37.5% ± 7% 1.2 46% ± 10% (100) 36% ± 10% (100) 41% ± 7% 1.1 45% ± 6% (250) 41% ± 6% (250) 43.4% ± 4.4% 1.0 49% ± 6% (300) 42% ± 6% (300) 45% ± 4% 0.9 47% ± 8% (150) 32% ± 8% (150) 40% ± 5.6% 0.8 40% ± 14% (50) 32% ± 13% (50) 36% ± 9.6% Table 7: Influence of FPU (70000 simulations/move).

We have discussed improvements that could be made to UCT algorithm. In particular, UCT does not help to choose a good ordering for non-visited moves, nor is it so effective for little explored moves. We proposed some methods adjusting the first-play urgency to solve this problem, and futher improvements are expected in this direction.

We have proposed the pattern-based (sequence-like) simulation which has improved significantly the level of MoGo. We implemented 3 × 3 patterns in random simulations in order to have more meaningful sequences. We believe that significant improvements can still be made in this direction, for example by using larger patterns, perhaps automatically generated ones. It is also possible to implement some sequence-forced patterns to improve the quality of simulations.

We have also shown the possibilities of combining UCT with pruning techniques in order to have a strong program on large Go board. Having had some encouraging results, we believe firmly further improvements in this direction.

A straightforward parallelization of UCT on shared-memory computer is made and has given some positive results. Parallelization on a cluster of computers can be interesting but the way to achieve that is yet to be found.

Figure 1 :

 1 Figure 1: UCT search. The shape of the tree enlarges asymmetrically. Only updated values (node[i].value) are shown for each visited nodes.

Figure 2 :

 2 Figure 2: Alpha-beta search with limited time. The node with ' ?' is the node not explored yet. This happens often during the large-sized tree search where entire search is impossible. Iterative deepening solves partially this problem.

Figure 3 :

 3 Figure 3: MoGo contains the tree search part using UCT and the random simulation part giving scores. The numbers on the bottom correspond to the final score of the game (win/loss). The numbers in the nodes are the updated values of the nodes (node[i].value)

Figure 4 :

 4 Figure 4: Left: beginning of one random game simulated by pure random mode. Moves are sporadically played with little sense. Right: beginning of one random game simulated by the pattern-based random mode. From move 5 to move 29 one complicated sequence is generated.

Figure 5 :

 5 Figure 5: Patterns for Hane. True is returned if any pattern is matched. In the right one, a square on a black stone means true is returned if and only if the eight positions around are matched and it is black to play.

Figure 6 :

 6 Figure 6: Patterns for Cut1. The Cut1 Move Pattern consists of three patterns. True is returned when the first pattern is matched and the next two are not matched.

Figure 7 :

 7 Figure 7: Pattern for Cut2. True is returned when the 6 upper positions are matched and the 3 bottom positions are not white.

Figure 8 :

 8 Figure 8: Patterns for moves on the Go board side. True is returned if any pattern is matched. In the three right ones, a square on a black (resp. white) stone means true is returned if and only if the positions around are matched and it is black (resp. white) to play.

Figure 9 :

 9 Figure 9: From one game played by MoGo (White) on KGS. In the left figure, it found almost the exact joseki response in the left-up corner using group and zone mode. In the right figure, the union of the two groups related to the last two moves (move 4 and move 5) are marked by squares (only on the free positions). The other area of zone (besides the stones and the group area) at this moment is marked by shadow.

Figure 10 :

 10 Figure 10: The opening of one game between MoGo and Indigo in the 18th KGS Computer Go Tournament. MoGo (Black) was in advantage at the beginning of the game, however it lost the game at the end.

Table 1 :

 1 Pseudocode of UCT for minimax tree.

	1: function playOneSequence(rootNode);
	2:	node[0] := rootNode; i = 0;
	3:	while(node[i] is not leaf) do
	4:	node[i+1] := descendByUCB1(node[i]);
	5:	i := i + 1;
	6:	end while ;
	7:	updateValue(node, -node[i].value);
	8: end function;
	9: function descendByUCB1(node)
	10:	nb := 0;
	11:
		1;
	26:	value := -value;
	27:	end for;
	28: end function;

for i := 0 to node.childNode.size() -1 do 12: nb := nb + node.childNode[i].nb; 13: end for; 14: for i := 0 to node.childNode.size() -1 do 15: if node.childNode[i].nb = 0 do v[i] := ∞; 16: else v[i] := -node.childNode[i].value /node.childNode[i].nb +sqrt(2*log(nb)/(node.childNode[i].nb) 17: end if ; 18: end for; 19: index := argmax(v[j]); 20: return node.childNode[index]; 21: end function; 22: function updateValue(node,value) 23: for i := node.size()-2 to 0 do 24: node[i].value := node[i].value + value; 25: node[i].nb := node[i].nb +

Table 2 :

 2 3.1 presents our application of UCT. Then, considering two important aspects for having a strong Monte-Pseudocode of UCT for MoGo

	INRIA

Table 3 :

 3 Different modes with 70000 random simulations/move in 9x9.

		. Rate	Win. rate	Total
		for B. Games for W. Games	Win. Rate
	Pure	46% (250)	36% (250) 41.2% ± 4.4%
	Sequence-like	77% (400)	82% (400)	80% ± 2.8%

Table 4 :

 4 MoGo with 70000 simulations per move, on 13 × 13 Go board, using or not the group mode heuristic.

		. Rate	Win. rate	Total
		for B. Games for W. Games	Win. Rate
	No group vs GnuGo level 0	53.2%(216)	51.8% (216)	52% ± 4.8%
	No group vs GnuGo level 8	24.2%(300)	30% (300)	27% ± 3.6%
	group vs GnuGo level 0	67.5% (80)	61.2% (80) 64.3% ± 7.5%
	group vs GnuGo level 8	51.9% (160)	60% (160)	56% ± 5.5%

Table 5

 5

		shows its level improves as this number increases. The
	outstanding performance of MoGo on double-processors and quadri-processors also supports
	this claim.
	Seconds	Winning Rate	Winning rate	Total
	per move for Black Games for White Games Winning Rate
	5 26% ± 12% (50) 26% ± 12% (50)	26% ± 8.7%
	20 41% ± 6% (250) 42% ± 6% (250) 41.8% ± 4.4%
	60 53% ± 7% (200) 50% ± 7% (200)	51.5% ± 5%

Table 5 :

 5 Pure random mode with different times.

Table 6 ,

 6 UCB1-TUNED is almost optimal in this sense.

	p	Winning Rate	Winning rate	Total
		for Black Games for White Games Winning Rate
	0.05	2% ± 4% (50)	4% ± 5% (50)	3% ± 3.4%
	0.55 30% ± 13% (50) 36% ± 13% (50)

Table 6 :

 6 Coefficient p decides the balance between exploration and exploitation. (Pure random Mode)

http://cgos.boardspace.net/

http://www.weddslist.com/kgs/past/index.html RR n • 6062

CrazyStone won the gold medal for the 9 × 9 Go game during the 11th Computer Olympiad at Turin 2006, beating several strong programs including GnuGo, Aya and GoIntellect.

We will use "play an arm" when refering to general multi-armed problems, and "play a move" when refering to Go. In Go application, the "play" will not refer to a complete game but only one move.INRIA

We will however say UCB1 for short.RR n •

Giving the pseudocode in order to have a clear explanation, we do not talk about the optimization here.

Here we use the original formula in Algorithm 1.INRIA

RR n • 6062

It is possible to arrive at one end game situation during the tree search part. In this case, one score could be calculated immediately and there is no need to create the node nor to call the random simulation part.INRIA

The ELO (http://en.wikipedia.org/wiki/Elo rating system) is a rating system where the probability of winning is related to the difference between the ranks.RR n • 6062

we had only access to a 4 processors computer, the behavior can be very different with many more processors.INRIA

Acknowledgments

We would like to thank Pierre-Arnaud Coquelin for the help during the development of MoGo. We specially thank Rémi Coulom for sharing his experiences of programming Crazy-Stone. We also appreciate the discussions on the Computer-Go mailing list.