D Gauvain

Marquet Bilel

Derbel Arnaud

Liefooghe El

Ghazali Talbi

Shake them all! Rethinking Selection and Replacement in MOEA

In this paper, we build upon the previous efforts to enhance the search ability of Moead (a decomposition-based algorithm), by investigating the idea of evolving the whole population simultaneously at once. We thereby propose new alternative selection and replacement strategies that can be combined in different ways within a generic and problem-independent framework. To assess the performance of our strategies, we conduct a comprehensive experimental study on biobjective combinatorial optimization problems. More precisely, we consider ρMNK-landscapes and knapsack problems as a benchmark, and experiment a wide range of parameter configurations for Moead and its variants. Our analysis reveals the effectiveness of our strategies and their robustness to parameter settings. In particular, substantial improvements are obtained compared to the conventional Moead.

Introduction

Evolutionary multi-objective optimization (EMO) algorithms [START_REF] Deb | Multi-Objective Optimization using Evolutionary Algorithms[END_REF][START_REF] Coello Coello | Evolutionary Algorithms for Solving Multi-Objective Problems[END_REF] have been proved extremely effective in computing a high-quality approximation of the Pareto set, i.e., the set of solutions providing the best compromises between the multiple objectives of an optimization problem. In particular, the so-called decomposition-based (or aggregation-based) EMO is gaining in popularity as an increasing number of studies is being devoted to their development [START_REF] Wagner | Pareto-, Aggregation-, and Indicator-based Methods in Many-objective Optimization[END_REF][START_REF] Zhang | MOEA/D: A Multiobjective Evolutionary Algorithm Based on Decomposition[END_REF][START_REF] Hughes | Multiple Single Objective Pareto Sampling[END_REF]. Recently, Moead [START_REF] Zhang | MOEA/D: A Multiobjective Evolutionary Algorithm Based on Decomposition[END_REF] (Multi-Objective Evolutionary Algorithm based on Decomposition) has attracted a lot of interest; which is due to its simplicity, approximation quality, and computational efficiency. In this paper, we seek for new alternative selection mechanisms for Moead at the aim of enhancing its search quality; and we focus on bi-objective combinatorial problems as a case study.

Generally speaking, Moead builds upon the idea of decomposing the initial multi-objective optimization problem into several single-objective sub-problems by means of scalarizing functions [START_REF] Miettinen | Nonlinear Multiobjective Optimization[END_REF] configured with different weight vectors. The most original part of Moead is to define, for each sub-problem, a neighborhood structure containing the set of the closest sub-problems. Then, Moead iterates over sub-problems and performs the three following basic steps: (i) select parents among the neighbors of the current sub-problem, (ii) generate an offspring by applying problem-specific operators, and (iii) replace neighbors' solutions if the generated offspring is better. We remark that mating selection (Step (i)) is performed exclusively among neighbors. Assuming that nearby sub-problems have similar solutions, the neighborhood size is critical for an accurate exploration/exploitation balance. Moreover, the replacement mechanism (Step (iii)) can lead to a situation where several neighbors are replaced by the same offspring. This can imply a loss of diversity, and likely a loss in performance. These issues have actually been addressed in [START_REF] Li | Multiobjective optimization problems with complicated pareto sets, MOEA/D and NSGA-II[END_REF], where two extra modifications have been introduced when dealing with complicated continuous Pareto sets. The first modification uses an extra probability parameter allowing parents to be selected from the whole population. The second one uses an extra parameter to bound the number of neighbors that can be replaced by a newly generated offspring.

In this paper, we propose new selection strategies to enhance the search ability of Moead for combinatorial problems. The idea behind our strategies stems from the observation that if an offspring is allowed to replace a neighboring solution in Moead, this neighboring solution is then 'lost' and it has no chance to get selected for reproduction in subsequent iterations. To overcome this shortcoming, we investigate an alternative perspective in Moead by optimizing all sub-problems at once. Intuitively, every solution from the population has a more fair chance to participate in the evolution process. This allows us to propose different strategies that can be plugged in the basic version of Moead. Our newly proposed strategies do not distort the basic framework of Moead, neither do they induce a loss in generality nor do they introduce new extra parameters; while being fully compatible with the previous modifications introduced so far. Moreover, they are proven to exhibit substantial improvements in approximation quality when compared with basic Moead and its modifications. Our performance assessment is in fact obtained as the byproduct of a thorough experimental analysis on two bi-objective combinatorial optimization problems, namely, knapsack and ρMNK-landscapes, and by considering a broad range of configurations. In the remainder, we first recall in Sec. 2 some basic definitions as well as a brief description of Moead. In Sec. 3, we describe our algorithmic contribution in designing new selection and replacement strategies for Moead. In Sec. 4, we present the settings of our experimental study. In Sec. 5, we state our main experimental findings. Finally, we conclude the paper in Sec. 6.

Background

Definitions. A multi-objective optimization problem can be defined by an objective function vector f = (f 1 , . . . , f m ) with m 2, and a set X of feasible solutions in the solution space. In the combinatorial case, X is a discrete set. Let Z = f (X) ⊆ IR m be the set of feasible outcome vectors in the objective space. To each solution x ∈ X is assigned an objective vector z ∈ Z, on the basis of the function vector f : X → Z with z = f (x). In a maximization context, a solution x ∈ X is dominated by a solution x ′ ∈ X iff ∀i ∈ {1, . . . , m}, f i (x) f i (x ′ ) and ∃i ∈ {1, . . . , m} such that f i (x) < f i (x ′ ). A solution x ⋆ ∈ X is said to be Pareto optimal (or non-dominated ), if there does not exist any other solution x ∈ X such that x ⋆ is dominated by x. The set of all Pareto optimal solutions is the Pareto set. Its mapping in the objective space is the Pareto front. Decomposition-based EMO. Contrary to existing Pareto-based EMO algorithms, like Nsga-ii or Spea2, which explicitly use the Pareto dominance relation in their selection mechanism, decomposition-based EMO algorithms [START_REF] Giagkiozis | Generalized Decomposition[END_REF] rather seek a good-performing solution in multiple regions of the Pareto front by decomposing the original multi-objective problem into a number scalarized singleobjective sub-problems, which can be solved independently as in Msops [START_REF] Hughes | Multiple Single Objective Pareto Sampling[END_REF], or in a dependent way as in Moead [START_REF] Zhang | MOEA/D: A Multiobjective Evolutionary Algorithm Based on Decomposition[END_REF]. Many different scalarizing functions have been proposed in the literature [START_REF] Miettinen | Nonlinear Multiobjective Optimization[END_REF]. Popular examples are the weighted sum (g ws ) and the weighted Tchebycheff (g te ) functions defined below.

g ws (x, λ) = m i=1 λ i • z ⋆ i -f i (x) ; g te (x, λ) = max i∈{1,...,m} λ i • z ⋆ i -f i (x)
where x ∈ X is an element from the solution space, λ = (λ 1 , . . . , λ m ) is a weighting coefficient vector such that λ i 0 for all i, and z ⋆ = (z ⋆ 1 , . . . , z ⋆ m ) is a utopian point, i.e., ∀i, ∀x, z ⋆ i > f i (x). Both functions are to be minimized. MOEA/D in a Nutshell. Let g be a scalarizing function and let (λ 1 , . . . , λ µ ) be a set of µ uniformly distributed weighting coefficient vectors, corresponding to µ sub-problems to be optimized. For each sub-problem i ∈ {1, . . . , µ}, the goal is to approximate the solution x with the best scalarizing function value g(x, λ i ). For that purpose, Moead maintains a population P = (p 1 , . . . , p µ ), each individual corresponding to a good-quality solution for one sub-problem. For each sub-problem i ∈ {1, . . . , µ}, a set of neighbors B(i) is defined with the T closest weighting coefficient vectors. To evolve the population, subproblems are optimized iteratively. At a given iteration corresponding to one sub-problem i ∈ {1, . . . , µ}, two solutions are selected at random from B(i), and an offspring solution x is created by means of variation operators (mutation and crossover). A problem-specific repair or improvement heuristic is potentially applied on solution x to produce x ′ . Then, for every sub-problem j ∈ B(i), if x ′ improves over j's current solution p j then x ′ replaces it; i.e., if g(x ′ , λ j ) < g(p j , λ j ) then set p j = x ′ . The algorithm continues looping over sub-problems, optimizing them one after the other, until a stopping condition is satisfied. We shall also consider the two modifications introduced in [START_REF] Li | Multiobjective optimization problems with complicated pareto sets, MOEA/D and NSGA-II[END_REF] to enhance Moead in the context of continuous complicated Pareto sets. The first one allows to select a parent from the whole population with a small probability parameter (1 -δ). More precisely, when dealing with a sub-problem i, its neighborhood is set to B(i) with probability δ, and to the whole population P with probability (1 -δ). The second one limits by a parameter nr the number of times that an offspring x ′ , created when dealing with a sub-problem i, can replace solutions in the neighborhood of i.

Rethinking Selection and Replacement in MOEA/D

As mentioned in the introduction, Moead could suffer from a lack of diversity due to the locality of its selection and replacement mechanism. We argue that this can also be caused by the fact that in Moead (and its modified variants), subproblems are optimized iteratively. In fact, since parents are selected randomly from the neighborhood of the sub-problem being processed, it might happen that a solution with the potential of producing a good offspring, gets never selected for reproduction. Additionally, because a neighbor's solution might be replaced as soon as a better offspring is found, this solution gets actually no chance to survive in the population. To increase the chance for a solution to survive in the population, we investigate the idea of evolving the whole population simultaneously by optimizing all subproblems in one shot and not iteratively. This idea is depicted in Algorithm 1 and discussed more thoroughly in the following.

o i ← crossover(p k , p ℓ ); o i ← mutation(o i ); 12 else o i ← mutation(p k ); 13 if o i is infeasible then repair(o i ); 14 for i ∈ {1, . . . , µ} do ci ← 0; 15 for i ∈ {1, . . . , µ} do /* Environmental replacement */ 16 if y = s then /* Selfish replacement */ 17 p ′ ← o i ; 18 if g(p ′ , λ i ) < g(p i , λ i ) then p i ← p ′ ; 19 else if y = c then /* Collective replacement */ 20 shuffle(Bi); 21 for j ∈ Bi do 22 p ′ ← o j ; 23 if cj < nr then 24 if g(p ′ , λ i ) < g(p i , λ i ) then p i ← p ′ ; cj ← cj + 1 ;
Algorithm 1 is mainly divided in two stages (lines 2 to 13 and lines 14 to 24). Contrary to Moead where a single offspring is generated at each iteration, our framework is basically a (µ+ µ)-EA where the first stage consists in generating µ offsprings and the second stage consists in updating the whole population for the next round. The first stage corresponds to mating selection where one new offspring is created for every subproblem. Specifically, we consider two alternatives: (i) either the solution of the current subproblem is always selected to be a parent and hence included for variation (x = s), or (ii) parents are picked randomly from neighbors in the usual way Moead proceeds (x = c). Moreover, every offspring is tagged with the identifier of the subproblem where it has been created. Thus, we can identify the subproblem that originated the creation of a given offspring. Only when all subproblems are treated and all µ new offspring solutions are created, the second stage of replacement occurs. In this stage, the subproblems are processed iteratively and we again consider two alternatives: (i) either the solution of a subproblem is compared to the offspring created at this subproblem (y = s), or (ii) the solution of the current subproblem is compared to the offspring created in neighboring subproblems (y = c). In both cases, the solution of a subproblem gets replaced if the considered offspring shows an improvement with respect to the scalarizing function of that particular subproblem.

Algorithm 1 is fully compatible with the baseline ideas of Moead; in particular, with the variants in [START_REF] Li | Multiobjective optimization problems with complicated pareto sets, MOEA/D and NSGA-II[END_REF], i.e., parameters δ and nr. Due to lack of space, we omit describing all the standard aspects that are shared with Moead, e.g., weights initialization, neighborhoods, update of the reference point, archiving.

To summarize, Algorithm 1 differs from Moead by essentially the fact that µ offsprings for all subproblems are created at each iteration. Moreover, since two alternatives are designed for mating selection and replacement, four different variants are possible: Moead-xy with x, y ∈ {s, c}s (resp. c) refers to a Selfish (resp. Collective) strategy where a subproblem privileges its own solution (resp. its neighbors' solutions). It is worth to notice that some parameter combinations may not have any impact on some algorithm variants, e.g., nr does not have an impact on Moead-ss and Moead-cs, neither δ on Moead-ss when cr = 0.

Experimental Setup

We analyze our approach on bi-objective ρMNK-landscapes and knapsack problems, with a broad range of instances with different structures and sizes. ρMNK-Landscapes. The family of ρMNK-landscapes constitutes a problemindependent model used for constructing multi-objective multi-modal landscapes with objective correlation [START_REF] Verel | On the structure of multiobjective combinatorial search space: MNK-landscapes with correlated objectives[END_REF][START_REF] Aguirre | Working principles, behavior, and performance of MOEAs on MNK-landscapes[END_REF][START_REF] Kauffman | The Origins of Order[END_REF]. A bi-objective ρMNK-landscape aims at maximizing an objective function vector f : {0, 1}

n → [0, 1] 2 . Solutions are binary strings of size n. The parameter k defines the number of variables that influence a particular position from the bit-string (the epistatic interactions). By increasing the number of variable interactions k from 0 to (n -1), landscapes can be gradually tuned from smooth to rugged. The objective correlation parameter ρ defines the degree of conflict between the objectives. The positive (resp. negative) data correlation allows to decrease (resp. increases) the degree of conflict between the objective function values. This has an impact on the cardinality of the Pareto front [START_REF] Verel | On the structure of multiobjective combinatorial search space: MNK-landscapes with correlated objectives[END_REF]. We investigate six random ρMNK-landscapes for each parameter combination given in Table 1. Knapsack. The knapsack problem is one of the most studied NP-hard problem. Given a collection of n items and a set of 2 knapsacks, the 0 -1 bi-objective bi-dimensional knapsack problem seeks a subset of items subject to capacity constraint based on a weight function vector w : {0, 1} n → N 2 , while maximizing a profit function vector p : {0, 1} n → N 2 . More formally, it can be stated as:

max n j=1 p ij • x j i ∈ {1, 2} s.t. n j=1 w ij • x j c i i ∈ {1, 2} x j ∈ {0, 1}
j ∈ {1, . . . , n} where p ij ∈ N is the profit of item j on knapsack i, w ij ∈ N is the weight of item j on knapsack i, and c i ∈ N is the capacity of knapsack i. We consider the standard instances proposed in [START_REF] Zitzler | Multiobjective evolutionary algorithms: A comparative case study and the strength pareto approach[END_REF], with random uncorrelated profit and weight integer values from [START_REF] Aguirre | Working principles, behavior, and performance of MOEAs on MNK-landscapes[END_REF]100], and where capacity is set to half of the total weight of a knapsack. Thirty different random problem instances are investigated for each parameter combination given in Table 1. Moreover, we use the same advanced weighted repairing procedure to handle constraints as in Moead [START_REF] Zhang | MOEA/D: A Multiobjective Evolutionary Algorithm Based on Decomposition[END_REF]. Parameter Setting. Table 1 shows the parameter settings investigated in our study. We consider the effect of the population size (µ), the neighborhood size (T ), the maximum number of neighboring solutions replaced (nr), the probability to select a parent outside of a neighborhood (1 -δ), the scalarizing function (g), and the crossover probability (cr). The stopping condition is set to 10 6 evaluation (resp. repair) calls for ρMNK-landscapes (resp. knapsack). Standard Moead [START_REF] Zhang | MOEA/D: A Multiobjective Evolutionary Algorithm Based on Decomposition[END_REF][START_REF] Li | Multiobjective optimization problems with complicated pareto sets, MOEA/D and NSGA-II[END_REF] is considered, together with our four variants. We use a bit-flip mutation (where each bit is independently flipped with a rate 1/n) and one-point crossover. The crossover probability parameter (cr) allows us to appreciate the impact of the variation operator, from a pure randomized local search algorithm (cr = 0.0) to a conventional genetic algorithm (cr = 1.0). The initial population is generated randomly. An unbounded archive of all non-dominated solutions is maintained with all the approaches. All algorithms have been executed under comparable conditions and share the same base components. Overall, we tested 15918 different configurations, each one executed 30 times. Due to space limitations, we only highlight a subset of settings allowing us to state our findings.

Experimental Analysis

Algorithm Comparison. We start by studying the approximation quality of competing algorithms. We follow the performance assessment protocol proposed by [START_REF] Knowles | A tutorial on the performance assessment of stochastic multiobjective optimizers[END_REF] using the hypervolume difference and multiplicative epsilon indicators [START_REF] Zitzler | Performance Assessment of Multiobjective Optimizers: An Analysis and Review[END_REF]. The hypervolume difference indicator (I - H ) gives the difference between the portion of the objective space that is dominated by the Pareto set approximation and some reference set. The reference point is set to the worst value obtained over all approximations, and the reference set is the best-found approximation over all tested configurations. The epsilon indicator (I × ε ) gives the minimum multiplicative factor by which the approximation found by an algorithm has to be translated in the objective space to weakly dominate the reference set.

Due to space limitations, we shall not focus on eliciting the best configuration for a specific instance; but give an overview of the differences between algorithms and their robustness to different parameters. To this end, a non-exhaustive set of results is shown in Tables 2 and3. First, notice the strong impact of the scalarizing function (g ws or g te ) on performance and its dependency on the considered problem. Overall, Moead-sc and Moead-cc are highly competitive and exhibit the most appealing behaviors. For knapsack, these two variants perform similarly to Moead. This can be explained by the relative strength of the repair function, and also by the shape of the Pareto front for knapsack problems, which is relatively easy to approximate. However, for ρMNK-landscapes with different structures and correlations, substantial improvements are reported, independently of the parameter setting. Actually, Moead-ss is also found to be competitive, but only when the crossover operator is activated. This is because Moead-ss degenerates to a multiple independent search in this case; and thus it is more likely trapped into independent local optima. At the opposite, Moeadsc and Moead-cc are able to adequately use information from neighbors, even when only a mutation (local search) operator is considered.

The previous discussion is in general valid when conducting an "anytime" analysis. This is illustrated in Fig. 1 rendering the typical convergence of competing algorithms. We see that all algorithms are able to make improvements, with Moead-sc and Moead-cc being consistently better than Moead. These results confirm that shaking many solutions at once can serve the approximation quality till the early stages of the search process. We also remark that Moead-sc and Moead-cc are more systematically improving upon Moead-ss for test instances having conflicting objectives; whereas Moead-ss is able to outperform its competitors as the objective correlation gets higher. Notice in fact that our strategies induce different intensification/diversification trade-offs both at the local level of every single-objective scalarized subproblem; but also at a more global level when considering the whole approximation set. When a selfish (resp. collective) mating selection is considered, the probability that a solution in the population gets selected for reproduction is 1 (resp. 1 -(1 -1/T ) T ). Roughly speaking, this means that all our strategies imply diversified offsprings since no solution in the current population gets replaced before exploring its potential. At the replacement stage, if a collective strategy is adopted, then the single-objective search at every subproblem is intensified since the probability that a locally improving solution can be found is higher. But this might increase the number of copies in the current approximation set. When a selfish replacement is considered, it is more likely that the number of copies is minimized; but at the price of delaying the advance of the population towards the front. For correlated objectives, and since the front is not too large, it is sufficient that only few solutions are able to approach the front in order to get good overall performance. Thus, a selfish replacement can be accurate. This is not the case for anti-correlated objectives where both the local improvements at every subproblem and the global spread of solutions is crucial. This explains the relative performance of our strategies depending on the characteristic of the tackled problem. Impact of Parameters. From Tables 2 and3, we can already extract some interesting observations on the impact of parameters, e.g., notice the differences between g ws and g te . Further observations from our data are sketched in Fig. 2, where only Moead, Moead-sc and Moead-cc are highlighted. First (Fig. 2 left), we confirm the positive impact of small values of parameter nr [START_REF] Li | Multiobjective optimization problems with complicated pareto sets, MOEA/D and NSGA-II[END_REF] on the performance of Moead for combinatorial problems. However, it was not always clear what is the best value to choose independently of the other parameters, e.g., the recommended value of 2 is in fact accurate, but not always optimal. Also, the impact of parameter nr on our strategies is rather mitigated. Although we found that it could often bring improvements, the impact was not pronounced compared to the case of Moead. For neighborhood size T (Fig. 2

middle),

Table 2. Relative performance of a representative subset of configurations w.r.t I - H and I × ε and ρMNK-landscapes (µ = 128, T = 8, δ = 1.0) at termination. For each row, the numbers indicates how many algorithm configurations (over the other 15 configurations given in columns) outperforms the configuration under consideration with a statistical confidence level of 0.05 (the lower, the better). we found that, contrary to knapsack, small values of T are interestingly more accurate for ρMNK-landscapes. We attribute this to the influence of the crossover operator which, combined with the repair mechanism, does enable to find highquality solutions for knapsack. However, for ρMNK-landscapes, it is more likely that the crossover fails in finding good solutions when selecting many often the parents from relatively large neighborhoods. Finally (Fig. 2 right), the impact of parameter δ for neighborhood selection is confirmed to have a positive impact on the performance. However, we find that another parameter has even more effect; namely the population size µ and especially for anti-correlated instances.

nr ∞ 2 - ∞ 2 - ∞ 2 ∞ 2 - ∞ 2 - ∞ 2 cr ρ K I - H 1.0 -0.
We attribute this to the fact that, when the Pareto front gets larger, it is more beneficial to increase the population size in order to avoid loosing promising individuals and to allow the population to be efficiently distributed.

Diversity Issues. We conclude our analysis by illustrating in Fig. 3 the size of the Pareto set approximation extracted at different iterations from the population (without the archive) for Moead and Moead-sc. In fact, we were able to observe that our strategy tends to maintain more spread solutions and to distribute them efficiently over the weight vectors, independently of the population size µ. We argue that this is a key feature of why our variants are able to exhibit better performance over Moead. Of course, this is not the only ingredient for optimal anytime performance; but it contributes much in finding a high-quality approximation, especially for conflicting objectives. This also suggests that other strategies controlling explicitly the number of clones in the whole populationtreated in a global manner-may lead to an even better performance.

Conclusions and Perspectives

We introduced a framework incorporating four strategies to deal with selection and replacement in Moead. Our experimental results show that substantial improvements can be obtained. Moreover, our strategies are found to be rather robust to Moead configuration parameters. We argue that our study opens new possibilities for improving the design of decomposition-based algorithms in several perspectives. Firstly, one can wonder whether more general (µ + λ)-EA can be accurately embedded in our framework. Secondly, we think that our framework opens the road to high-quality local search-based Moead variants for combinatorial optimization problems, e.g., plugging several (1 + λ)-EAs within each subproblem. Finally, our strategies are inherently distributed in the sense that each sub-problem is optimized concurrently in parallel with the others. In this respect, an interesting research issue would be to investigate the effective parallelization of our strategies and to study their running time scalability.

Algorithm 1 : 4 else 6 k ← i; 7 else 8 k 9 ℓ

 146789 Our proposed framework Moead-xy (x, y ∈ {s, c}) Input: λ 1 , . . . , λ µ : weight vectors w.r.t sub-problems; g: a scalarizing function; B(i): the neighbors of sub-problem i ∈ {1, . . . , µ}; P = p 1 , . . . , p µ : the initial population. 1 while Stopping Condition do 2 for i ∈ {1, . . . , µ} do 3 if rand(0, 1) < δ then Bi ← B(i) ; /* Neighborhood Setting */ Bi ← P ; 5 if x = s then /* Selfish mating selection */ if x = c then /* Collective mating selection */ ← rand(Bi); ← rand(Bi); while ℓ = k do ℓ ← rand(Bi); 10 if rand(0, 1) < cr then /* Variation operators */ 11

Fig. 1 .

 1 Fig. 1. Convergence plots w.r.t. hypervolume difference and ρMNK-landscapes. Columns are respectively for ρ ∈ {-0.7, 0.0, 0.7}. Rows are respectively for K ∈ {4, 8}. Results are for µ = 128, T = 8, δ = 1.0, cr = 1.0, nr = 0 and g te .

Fig. 2 .

 2 Fig. 2. Impact of parameters. The two subfigures in the left (resp. middle, right) show the impact of nr (resp. T , µ combined with δ) as shown in the vertical axis for every algorithm. The boxplots are w.r.t. the indicator depicted at the horizontal axis. Results are for a ρMNK-landscape with ρ = -0.7 and K = 4. Whenever not explicit: µ = 128, T = 8, cr = 1.0, δ = 1.0, nr = 2 and g te .

Fig. 3 .

 3 Fig. 3. Solution diversity. The y-axis gives the number of different solutions in the population (in log-scale). Results are for a ρMNK-landscape with ρ = -0.7 and K = 8 (µ = 128, T = 8, cr = 1.0, δ = 1.0, nr = 2 and g te ).

Table 1 .

 1 Parameter setting.

			ρMNK-landscapes		Knapsack
			m = 2, n = 128		m = 2
			ρ ∈ {-0.7, 0.0, 0.7}, K ∈ {4, 8} n = 250 n = 500 n = 750
	pop size	µ	64, 128, 256	150	200	250
	neighborhood size	T	4, 8, 16, 32		10, 20, 30
	max. number of replacements nr	1, 2, 3, 4, ∞		2, 4, 8, 10, ∞
	neighborhood probability	δ	0.9, 1.0
	crossover rate	cr	0.0, 0.9, 1.0
	scalarizing function	g	weighted sum (g ws ), weighted Tchebycheff (g te )
	stopping condition		10 6 evaluation function calls 10 6 repair procedure calls

Table 3 .

 3 Relative performance of a representative subset of configurations for knapsack (T = 20, δ = 0.9). Similar metrics than in Table.2 are reported.

		7 4	2 2	8	0	0	5	0	0			4 5	0	6	2	7	5	6
		-0.7 8	3 1	7	0	0	6	1	0 11 7	6	9	6	9	9	8
		0.0 4	2 1	0	0	0	2	0	0			2 0	0	0	0	0	0	2
		0.0 8	6 0	0	0	1	6	0	0			8 7	3	8	9	12	6	0
		0.7 4 14 10	0	2	2	2	2	1			1 2	0	2	2	2	2	1
		0.7 8	2 2	0	2	1	2	2	2			2 1	0	5	2	1	2	1
		-0.7 4	8 8	1	0	0	1	0	0	I × ε	0 5	0	5	0	8	2	3
		-0.7 8	5 5	0	0	0	5	5	0			7 6	0	9	5	7	7	5
		0.0 4	2 1	0	0	0	1	0	0			0 0	1	1	1	1	0	1
		0.0 8	4 0	0	0	1	3	0	0			9 9	6	9	9	10	7	7
		0.7 4 15 4	0	2	2	2	2	0			2 2	0	2	2	2	3	0
		0.7 8	1 1	0	2	0	1	1	2			2 1	1	6	2	1	4	2
		-0.7 4	0 4	15	0	0	4	0	0	I -H	6 5	13	5	5	12	6	6
		-0.7 8	1 1	15	1	0	6	0	0			9 6	6	6	6	10	6	8
		0.0 4	3 0	15	0	0	0	0	0			3 0	11	0	0	0	0	0
		0.0 8	0 0	15	0	0	4	0	0			6 8	6	8	5	6	6	5
		0.7 4	9 1	9	1	0	1	1	0			3 1	2	0	0	0	0	0
		0.7 8	0 1	0	4	0	1	1	1			2 1	0	1	1	1	4	4
	0.0	-0.7 4	6 8	10	1	1	1	1	1	I × ε	6 2	0	4	5	10	2	4
		-0.7 8	3 5	5	1	1	5	1	1			8 5	0	7	5	10	5	9
		0.0 4	3 1	13	0	0	1	0	0			1 1	1	0	0	1	1	1
		0.0 8	0 0	8	0	0	0	0	0			8 8	6	8	8	8	8	8
		0.7 4	6 0	4	0	0	1	0	0			3 2	3	0	0	0	0	0
		0.7 8	0 1	0	3	0	1	0	0			2 0	0	2	2	0	10 6
						g te										g ws
		M o e a d	M o e a d -s s	M o e a d -s c	M o e a d -c s	M o e a d -c c		M o e a d	M o e a d -s s	M o e a d -s c	M o e a d -c s	M o e a d -c c
		nr ∞ 2	-		∞	2	-	∞	2		∞ 2	-	∞	2	-	∞	2
	cr N 250 8 8	7		8	8	8	8	8	I -H	0 0	0	0	0	0	0	0
		500 8 8	8		8	8	13	8	8			0 0	0	0	0	0	0	0
		750 8 8	8		8	8	14	8	8			0 0	0	0	0	0	0	0
	0.9	250 0 0	0		0	0	4	0	0	I × ε	0 0	0	0	0	0	0	0
		500 8 8	7		7	7	11	7	7			0 5	0	0	4	1	4	5
		750 4 3	6		4	3	7	3	2			5 0	0	5	7	7	0	8


Aymeric Stamm

Olivier Commowick

Patrick Pérez

Christian Barillot

FAST IDENTIFICATION OF OPTIMAL FASCICLE CONFIGURATIONS FROM STANDARD CLINICAL DIFFUSION MRI USING AKAIKE INFORMATION CRITERION

Keywords: diffusion MRI, multi-compartment models, model selection, model averaging, semioval center

Analytic multi-compartment models have gained a tremendous popularity in the recent literature for studying the brain white matter microstructure from diffusion MRI. This class of models require the number of compartments to be known in advance. In the white matter however, several non-collinear bundles of axons, termed fascicles, often coexist in a same voxel. Determining the optimal fascicle configuration is a model selection problem. In this paper, we aim at proposing a novel approach to identify such a configuration from clinical diffusion MRI where only few diffusion images can be acquired and time is of the essence. Starting from a set of fitted models with increasing number of fascicles, we use Akaike information criterion to estimate the probability of each candidate model to be the best Kullback-Leibler model. These probabilities are then used to average the different candidate models and output an MCM with optimal fascicle configuration. This strategy is fast and can be adapted to any multicompartment model. We illustrate its implementation with the ball-and-stick model and show that we obtain better results on single-shell low angular resolution diffusion MRI, compared to the state-of-the-art automatic relevance detection method, in a shorter processing time.

INTRODUCTION

Diffusion magnetic resonance imaging (dMRI) is a tool of choice for the analysis of the brain white matter (WM). The diffusion of thermically agitated water molecules in the axons is indeed restricted by the cylindrical geometry of such cells with largest diameters around 20 µm. The typical spatial resolution in clinical dMRI is 2 × 2 × 2 mm 3 cubic voxels, which implies that hundreds of non-collinear axons can presumably populate a given voxel. Axons are usually assumed to be distributed in orientationally homogeneous bundles termed fascicles. Water in a voxel then populates the different fascicle environments and the extra-axonal space in unknown proportions. Multi-compartment models (MCM) rely on finite mix-tures of distributions and are thus well suited to model the voxelwise diffusion-induced molecular displacements. Existing MCMs include the 47 models reviewed in [START_REF] Panagiotaki | Compartment models of the diffusion MR signal in brain white matter: a taxonomy and comparison[END_REF] and Diffusion Directions Imaging [START_REF] Stamm | A new multi-fiber model for low angular resolution diffusion MRI[END_REF]. The number of fascicles however needs to be known in advance, which is not the case in practice. So far, this model selection problem has been solved either by brute-force methods or using Bayesian frameworks.

In brute-force approaches, a set of nested candidate MCMs with increasing number of fascicles is fitted to the data. The best MCM is then identified as the candidate model that best fits the data, where the comparison usually relies on an F-test [START_REF] Alexander | Detection and modeling of non-Gaussian apparent diffusion coefficient profiles in human brain data[END_REF]. Since the more complex the model, the better the fit, the F-test often tends to favor MCMs that overfit the data when the same data is used for estimation and to assess goodness of fit. To limit overfitting, the Bayesian information criterion has been introduced to penalize model complexity that increases with the number of fascicles [START_REF] Schultz | Multidiffusion-tensor fitting via spherical deconvolution: a unifying framework[END_REF]. Recently, the generalization error has been proposed to choose the "optimal" MCM based on its ability to predict new data [START_REF] Scherrer | Reliable selection of the number of fascicles in diffusion images by estimation of the generalization error[END_REF], thus avoiding the overfitting issue. These approaches limit their search of the "optimal" MCM to a predefined candidate set.

Differently, Bayesian frameworks try to estimate the best MCM as the one that maximizes a posterior distribution on the models. They rely on a careful choice of a prior distribution for MCMs. For instance, [START_REF] Demiralp | Generalizing diffusion tensor model using probabilistic inference in Markov random fields[END_REF] uses Markov random fields (MRF) while [START_REF] Behrens | Probabilistic diffusion tractography with multiple fibre orientations: What can we gain?[END_REF] resorts to Automatic Relevance Detection (ARD) in which non-informative priors are assigned to all the MCM parameters except the mixture weights that are assumed to be Beta-distributed. Such priors automatically prune an entire compartment if it is not supported by the data. These methods simultaneously perform model estimation and selection. When translated to clinics, Bayesian methods however have limitations. They are prohibitively computationally expensive: ARD (resp., MRF) requires MCMC simulations (resp., optimization) to come up with a usable posterior distribution. Moreover, dMRI clinical protocols often include a single b-value and a maximum of 30 diffusion-sensitizing gradients (DSG) [START_REF] Jones | The effect of gradient sampling schemes on measures derived from diffusion tensor MRI: a Monte Carlo study[END_REF]. With such small sample sizes, the posterior distribution strongly depends on the prior, making the Bayesian information updating potentially ineffective.

In this paper, we propose to combine the two previous approaches. Starting with a candidate set, we search for the best Kullback-Leibler (KL) MCM (i.e., at minimal KL divergence to the true unknown diffusion model) as a combination of the candidate models. In section 2.1, we show that we can estimate the probability of any MCM to be the best KL MCM using only the Akaike Information Criterion (AIC) and we explain how these probabilities can be used to generate the best KL MCM from the candidate set. In section 2.2, we implement our method for the ball-and-stick model [START_REF] Behrens | Probabilistic diffusion tractography with multiple fibre orientations: What can we gain?[END_REF], which allows us to set up a pilot study (Section 3) to compare our results to the ARD method on standard clinical data. Results are presented in section 4 and discussed in section 5.

PROPOSED APPROACH

Model averaging from Akaike information criterion

The AIC is widely used in brute force approaches where the best model is often chosen as the one with minimal AIC. However, by definition, it has a fundamental role in model averaging [START_REF] Posada | Model selection and model averaging in phylogenetics: advantages of Akaike information criterion and Bayesian approaches over likelihood ratio tests[END_REF] that this section aims at summarizing.

The AIC is an asymptotically unbiased estimator the relative expected KL divergence between an estimated model and the true unknown model. It is defined as follows:

AIC = -2 log L + 2K , (1)
where L is the maximized value of the likelihood for the estimated model and K is the number of parameters of the estimated model. When the sample size N is small, the following corrected AIC is recommended to avoid overfitting:

AIC c = -2 log L + 2K + 2K(K + 1) N -K -1 . (2)
Since the true model is unknown, its AIC is approximated by the minimal AIC among those of candidate models. Let {M ℓ } ℓ=0,...,L be the set of candidate models. The AIC difference ∆ ℓ = AIC c (ℓ) -min k∈ 0,L AIC c (k) of model M ℓ can then be used to define its so-called Akaike weight w ℓ :

w ℓ = exp {-∆ ℓ /2} L k=0 exp {-∆ k /2} . (3)
The Akaike weight w ℓ approximates the probability for M ℓ to be the best KL model. Such probabilities happen to be very useful. For instance, if we sort the candidate models from largest to smallest Akaike weight, we can establish a 95% confidence set of models by keeping the first p models such that p ℓ=0 w ℓ = 0.95. An evidence ratio between model M ℓ1 and M ℓ2 can also be estimated as w ℓ1 /w ℓ2 and used to assess how strongly M ℓ1 is more likely to be the best KL model compared to M ℓ2 . Finally, a model-averaged estimate of any parameter θ across all candidate models can be obtained as follows:

θ = L ℓ=0 w ℓ 1 Θ ℓ (θ) θ ℓ L ℓ=0 w ℓ 1 Θ ℓ (θ) , ∀θ ∈ Θ 0 ∪ Θ 1 ∪ • • • ∪ Θ L , (4)
where 1 Θ ℓ is the characteristic function of the parameter set Θ ℓ of model M ℓ and θ ℓ is the estimate of θ under model M ℓ .

S (ℓ) i S 0 =   1 - ℓ j=1 f (ℓ) j   e -bd (ℓ) + ℓ j=1 f (ℓ) j e -bd (ℓ) u T i µ (ℓ) j 2 , (5)
in which S 0 is the expected MR signal in absence of DSG and Θ ℓ = ±µ

(ℓ) 1 , f (ℓ) 1 , . . . , ±µ (ℓ) ℓ , f (ℓ) ℓ , d (ℓ)
is the parameter set of the ℓ-fascicle ball-and-stick MCM where:

• the fascicle orientations ±µ

(ℓ) 1 , . . . , ±µ (ℓ)
ℓ ∈ S 2 ; • the fascicle occupancies f

(ℓ) 1 , . . . , f (ℓ) ℓ ∈ [0, 1]; • the free diffusivity d (ℓ) ≥ 0. When ℓ = 0, it reads S (0) i = S 0 e -bd (0) and Θ 0 = d (0) .
Model averaging. We can observe that, for ℓ ∈ 0, L , only the parameter d appears in more than one MCM with the same interpretation as the free diffusivity (it actually appears in all of them). Indeed, MCMs with an increasing number ℓ of fascicles are nested and there is thus no pairwise matching of the fascicle compartments between two different MCMs, making model-averaged estimates of the fascicle orientations/occupancies hard to define. For example, when L = 3, the unique fascicle of the 1-fascicle MCM M 1 can be averaged either with the first or with the second fascicle of the 2-fascicle MCM M 2 . Each of these two combinations can further be averaged with one of the 3 fascicles of the 3fascicle MCM M 3 , leading to a total of 6 averaged fascicles. In general, a candidate set made of MCMs from 0 to L fascicle compartments can generate up to L! fascicles. To address this issue, we express the candidate MCMs with a same fixed number L! of fascicle compartments, while ensuring a pairwise matching of the compartments between the MCMs. This is achieved by the indexing k = (m-1) L! (ℓ-1)! +(j -1) L! ℓ! +p where p ∈ 1, L! ℓ! and m ∈ 1, (ℓ -1)! , which yields the following reformulation of eq. ( 5):

S (ℓ) i S 0 = 1 - L! k=1 f (ℓ) L,k e -bd (ℓ) + L! k=1 f (ℓ) L,k e -bd (ℓ) u T i µ (ℓ) L,k 2 , (6) where f
(ℓ) L,k = ℓ L! f (ℓ) j and µ (ℓ) L,k = µ (ℓ)
j . Note that this formulation remains valid for ℓ = 0. Finally, we can apply eq. ( 4) to eq. ( 6) to obtain model-averaged estimates of the ball-andstick parameters as follows:

d = L ℓ=0 w ℓ d (ℓ) , f L,k = L ℓ=0 w ℓ f (ℓ) L,k and µ L,k = E 1 [D L,k ] with D L,k = L ℓ=1 w ℓ µ (ℓ) L,k µ (ℓ) T L,k L ℓ=1 w ℓ , (7)
where E 1 [•] is the principal eigenvector operator.

MATERIAL & METHODS

Description of the data

Raw data. We conducted a pilot study in which two healthy volunteers (S1 and S2) underwent a series of 10 diffusion MRI scans on the same MR scanner (Siemens 3T Verio) with the same protocol. This protocol lasted 7 min and comprised a single non-weighted diffusion image B0 and 30 diffusion-weighted (DW) images acquired at b = 1000 s/mm 2 along 30 non-collinear DSG directions u i uniformly spread over the north hemisphere. The following parameters were used: 128 × 128 × 60 image resolution with 2 × 2 × 2 mm 3 voxels, TR= 11 s and TE= 99 ms.

Preprocessing. We preprocessed the data with FSL [10]: (i) For each scan, we performed (a) a rigid registration of the DW images on the B0 to correct for subject motion and rotated the gradient tables accordingly and (b) an affine registration of the DW images on the B0, guided by the previously estimated rigid transformation, to correct for distortions. (ii) We performed a rigid registration of the B0 images of the different scans on the B0 image of the first scan and applied the corresponding transformation to the subsequent DW images. We rotated the gradient tables accordingly. (iii) For each scan, we reduced the noise in the images using the Rician-adapted non-local means filter [START_REF] Wiest-Daesslé | Rician noise removal by nonlocal means filtering for low signal-to-noise ratio MRI: applications to DT-MRI[END_REF]. (iv) We extracted the brain using the BET algorithm. (v) We computed a WM mask using the FAST algorithm.

Experimental setup

Maximal number of candidate MCMs. According to known anatomy [START_REF] Kretschmann | Cranial neuroimaging and clinical neuroanatomy: magnetic resonance imaging and computed tomography[END_REF], the semioval center is an area where association, commissural and projection fascicles cross. The number L of fascicle compartments in the candidate models thus needs to be at least 3 (L ≥ 3). We aim at showing that it is reasonable to set L = 3. To this end, we observed that the estimate of free water occupancy is highly biased for low L and gets better and stabilizes as L increases. Hence, we performed the model averaging procedure for L = 3, 4 and compared the resulting free water occupancy maps by means of Dice coefficient for different occupancy thresholds ranging in ]0, 1[: if these two maps are highly similar, the contribution of the 4-fascicle MCM to the best MCM is negligible.

Comparison to ARD. Next, we propose to compare our approach to the ARD method as implemented in FSL [START_REF] Jenkinson | FSL[END_REF] using L = 3. ARD alters the estimation of the fascicle occupancies in such a fashion that f ARD L estimates the probability that the L-fascicle MCM contributes to the best unknown MCM. With our approach, this probability is given by the Akaike weight w L of the L-fascicle MCM. In order to compare f ARD 3 and w 3 ,we computed these indices on 100 bootstrap DW images generated out of the original 10. For both indices, we then performed an FDR-corrected z-test of nullity in all WM voxels. The aim is to compare the areas where the methods detect a statistically significant contribution of the 3-fascicle MCM to the best MCM, with respect to known anatomy. To this end, we focused on the semioval center and computed the mean of both indices wherever significantly non null. We also provide a qualitative visualization of the estimated fascicles.

RESULTS

Maximal number of candidate MCMs. We computed the Dice coefficient between binarized free water occupancy maps obtained with L = 3, 4. Results are summarized in table 1. For all thresholds employed, the Dice coefficient remains greater than 0.95, which means that including a 4-fascicle MCM in the averaging does not bring much additional information. We can thus reasonably conclude that the contribution of the 4-fascicle MCM to the best MCM is negligible in the WM. 1. Dice similarity index between free water occupancy (f 0 ) maps obtained assuming L = 3, 4 fascicle compartments. Different thresholds were used to binarize the maps.

f
Comparison to ARD. Figure 1 shows the results of the model selection obtained with our approach (left) and with ARD (right) on S1. A global coronal view of the statistically significant mean probabilities w 3 (Figure 1(a)) and f ARD 3 (Figure 1(b)) with same window level reveals that our approach offers a more robust detection of statistically significant contributions of the 3-fascicle MCM to the best MCM. We also zoomed on the semioval center (Figures 1(c) and 1(d)) where one expects to see predominantly 3 distinct fascicles. Results show that ARD hardly identifies voxels with 3 fascicles. It also requires to set a threshold on the fascicle occupancies (here, 0.05) to discard fascicles with too low occupancy, on which the resulting fascicle configurations highly depend. In contrast, our approach does not need to set any parameters and, most of the time, correctly identifies 1 fascicle in the corpus callosum and 3 fascicles in the semioval center. Comparative performances were identical for S2.

DISCUSSION

We presented a new approach to the model selection problem that enables the use of MCMs in a clinical setting. We combined ideas from brute-force approaches starting with a set of candidate models and ideas from Bayesian estimation quantifying the probability of each candidate to be the best one. We focused on the best KL model, which only requires AIC computation which is instantaneous. The use of a bias-corrected AIC makes the procedure robust to small sample sizes thus suited to clinical data. On such data, our method proved to be more efficient than ARD in correctly identifying 3-fascicle crossing areas. It is worth noting that no user-defined parameters are required contrary to ARD. From a computational point of view, the averaging procedure is very fast (a few seconds) but requires the estimation of MCMs up to 3 fascicles, which is a longer step. Nonetheless, estimation and averaging took about 5 hours for ARD and only 1 hour for our approach on an 8-core computer, which makes it closer to a clinically acceptable running time. Future works will focus on the accuracy of the estimated fascicle orientations. This is a tricky problem since (i) no ground truth is available and (ii) the op-timal fascicle configurations might be model-dependent.

2. 2 .

 2 Ball-and-stick model averaging The ℓ-fascicle ball-and-stick model. It assumes that the expected MR signal S (ℓ) i induced by a DSG with b-value b and direction u i has the following parametric form for ℓ ≥ 1 [7]:

Fig. 1 . 3 (

 13 Fig. 1. Model selection results in the semioval center. Coronal view of the mean probabilities w 3 (left) and f ARD 3 (right) masked with the result of FDR-corrected z-test of nullity. Second row focuses on the semioval center and estimated fascicles are overlaid as equally shaped cones. Fascicle orientations are color-coded and match the cone orientations.