Xiangjing Lai
email: laixiangjing@gmail.com

Jin-Kao Hao
email: jin-kao.hao@univ-angers.fr

Zhang-Hua Fu
email: fuzhanghua@cuhk.edu.cn

Dong Yue

Neighborhood decomposition based variable neighborhood search and tabu search for maximally diverse grouping

Keywords: Heuristics, graph grouping problem, neighborhood decomposition, local search

The maximally diverse grouping problem (MDGP) is a relevant NP-hard optimization problem with a number of real-world applications. However, solving large instances of the problem is computationally challenging. This work is dedicated to a new heuristic algorithm for the problem, which distinguishes itself by two original features. First, it introduces the rst neighborhood decomposition strategy to accelerate neighborhood examinations. Second, it integrates, in a probabilistic way, two complementary neighborhood decomposition based local search procedures (variable neighborhood descent and tabu search) as well as an adaptive perturbation strategy to ensure a suitable balance between intensication and diversication of the search space. Computational results on 320 benchmark instances commonly used in the literature show that the proposed algorithm competes favorably with the state-ofthe-art MDGP algorithms, by reporting improved best-known results (new lower bounds) of the literature for 220 large instances. Additional experiments are conducted to analyze the main components of the algorithm. The proposed algorithm can help to better solve practical problems that can be formulated by the maximally diverse grouping model.

Introduction

Given a set V of N elements, a distance matrix D = [d ij ] N ×N between the elements, a positive integer m, and m pairs of non-negative integers {L g , U g } (1 ≤ g ≤ m) called the capacity lower and upper limits of groups, the maximally diverse grouping problem (MDGP) is to partition the set V into m disjoint groups such that the size of each group g lies in [L g , U g ] (1 ≤ g ≤ m), while the sum of the distances between the elements in the same groups is maximized. MDGP can also be described as a graph partition problem as follows. We consider an edge-weighted complete graph G = (V, E, D), where V is the set of N vertices, E is the set of N ×(N -1)/2 edges, and D = [d ij ] N ×N denes the set of edge weights. Then MDGP can be considered as a special case of the NP-hard clique partitioning problem (CPP) [START_REF] Grötschel | A cutting plane algorithm for a clustering problem[END_REF]29] with non-negative edge weights and constraints related to the capacity lower and upper limits of groups.

Formally, MDGP can be written as a quadratic binary programming problem [START_REF] Gallego | Tabu search with strategic oscillation for the maximally diverse grouping problem[END_REF][START_REF] Rodriguez | An articial bee colony algorithm for the maximally diverse grouping problem[END_REF]:

Maximize m g=1 N -1 i=1 N j=i+1 d ij X ig X jg (1)
Subject to m g=1

X ig = 1, i = 1, 2, . . . , N (2)
L g ≤ N i=1
X ig ≤ U g , g = 1, 2, . . . , m

X ig ∈ {0, 1}, ∀i ∈ {1, 2, . . . , N }, ∀g ∈ {1, 2, . . . , m},

where X ig is a binary variable that takes 1 if the vertex i locates in the group g and 0 otherwise, the set of constraints [START_REF] Benlic | Breakout local search for the vertex separator problem[END_REF] guarantees that each vertex is located in exactly one group, and the set of constraints [START_REF] Benlic | Breakout local search for the max-cut problem[END_REF] ensures that the size of group g lies in [L g , U g ] (g = 1, 2, . . . , m).

MDGP is a relevant model for formulating many practical problems, such as assignment of students to groups [START_REF] Johnes | Operational Research in education[END_REF][START_REF] Krass | Constrained group balancing: Why does it work[END_REF]28], creation of peer review groups [START_REF] Chen | A hybrid grouping genetic algorithm for reviewier group construction problem[END_REF], and VLSI design [START_REF] Weitz | An empirical comparison of heuristic and graph theoretic methods for creating maximally diverse groups, VLSI design, and exam scheduling[END_REF]. More applications can be found in [START_REF] Gallego | Tabu search with strategic oscillation for the maximally diverse grouping problem[END_REF][START_REF] Lai | Iterated maxima search for the maximally diverse grouping problem[END_REF][START_REF] Palubeckis | Maximally diverse grouping: an iterated tabu search approach[END_REF][START_REF] Rodriguez | An articial bee colony algorithm for the maximally diverse grouping problem[END_REF][START_REF] Uro²evi¢ | Variable neighborhood search for maximum diverse grouping problem[END_REF].

Due to the NP-hardness of MDGP [START_REF] Feo | A class of bounded approximation algorithms for graph partitioning[END_REF], a number of heuristic algorithms have been proposed to nd approximate solutions. Existing heuristic algorithms can be classied into two categories, i.e., trajectory-based local search algorithms and hybrid evolutionary algorithms. As examples of trajectory-based local search algorithms, we mention the multistart algorithm [START_REF] Arani | A three phased approach to nal exam scheduling[END_REF], Lot-Cerveny-Weitz (LCW) algorithm [27], T-LCW method mixing LCW and tabu search [START_REF] Gallego | Tabu search with strategic oscillation for the maximally diverse grouping problem[END_REF], simulated annealing algorithm [START_REF] Palubeckis | Comparative performance of three metaheuristic approaches for the maximally diverse grouping problem[END_REF], variable neighborhood search algorithms [START_REF] Brimberg | Solving the maximally diverse grouping problem by skewed general variable neighborhodd search[END_REF][START_REF] Palubeckis | Comparative performance of three metaheuristic approaches for the maximally diverse grouping problem[END_REF][START_REF] Uro²evi¢ | Variable neighborhood search for maximum diverse grouping problem[END_REF], tabu search algorithms [START_REF] Gallego | Tabu search with strategic oscillation for the maximally diverse grouping problem[END_REF][START_REF] Palubeckis | Maximally diverse grouping: an iterated tabu search approach[END_REF], and iterated maxima search (IMS) algorithm [START_REF] Lai | Iterated maxima search for the maximally diverse grouping problem[END_REF]. Hybrid evolutionary algorithms include hybrid genetic algorithms [START_REF] Fan | A hybrid genetic algorithmic approach to the maximally diverse grouping problem[END_REF][START_REF] Singh | A new hybrid genetic algorithm for the maximally diverse grouping problem[END_REF], hybrid grouping genetic algorithm [START_REF] Chen | A hybrid grouping genetic algorithm for reviewier group construction problem[END_REF], hybrid steady-state genetic algorithm [START_REF] Palubeckis | Comparative performance of three metaheuristic approaches for the maximally diverse grouping problem[END_REF], articial bee colony (ABC) algorithms [START_REF] Rodriguez | An articial bee colony algorithm for the maximally diverse grouping problem[END_REF], and constructive genetic algorithm [START_REF] Antonio | Constructive genetic algorithm for clustering problems[END_REF]. According to the experimental results reported in recent studies such as [START_REF] Brimberg | Solving the maximally diverse grouping problem by skewed general variable neighborhodd search[END_REF][START_REF] Lai | Iterated maxima search for the maximally diverse grouping problem[END_REF][START_REF] Singh | A new hybrid genetic algorithm for the maximally diverse grouping problem[END_REF], the skewed general variable neighborhood search algorithm (SGVNS) [START_REF] Brimberg | Solving the maximally diverse grouping problem by skewed general variable neighborhodd search[END_REF], the iterated tabu search algorithm (ITS) [START_REF] Palubeckis | Maximally diverse grouping: an iterated tabu search approach[END_REF], the iterated maxima search algorithm (IMS) [START_REF] Lai | Iterated maxima search for the maximally diverse grouping problem[END_REF] can be regarded as the current state-of-the-art MDGP algorithms. Finally, the hybrid genetic algorithm NSGGA [START_REF] Singh | A new hybrid genetic algorithm for the maximally diverse grouping problem[END_REF] can also be considered as a state-of-the-art algorithm, but only for the special case where all groups have an equal size.

One notices that the best performing MDGP algorithms in the literature rely, with no exception, on a powerful neighborhood search subroutine. Meanwhile, a careful analysis of the underlying neighborhood search procedures indicates that they examine the whole neighborhood used at each iteration, which leads to superuous examinations of many non-promising neighbor solutions and a waste of computation time. To overcome this problem, this work investigates an original neighborhood decomposition method which avoids redundant calculations of the neighborhood examination and thus speeds up the neighborhood search.

The main contributions of this work are summarized as follows.

(1) We propose a novel neighborhood decomposition based heuristic algorithm (NDHA) with two original features. First, NDHA relies a dynamic neighborhood decomposition strategy that allows the algorithm to avoid redundant examination of irrelevant neighbor solutions at each iteration by ignoring uninteresting candidate solutions. As such, the neighborhood decomposition strategy accelerates the neighborhood examination and enables more promising candidate solutions to be examined for a given time budget. Second, the neighborhood decomposition based heuristic algorithm integrates two complementary neighborhood search procedures (i.e., tabu search and variable neighborhood descent) that are applied in a probabilistic way, leading to an enhanced robustness of the algorithm. [START_REF] Benlic | Breakout local search for the vertex separator problem[END_REF] We present computational results of the proposed algorithm on 320 benchmark instances commonly used in the literature and compare our results with those of the state-of-the-art MDGP algorithms. Our comparative studies indicate that the proposed algorithm outperforms signicantly the reference algorithms especially on the large benchmark instances.

The new lower bounds for 220 benchmark instances reported by our algorithm are useful for assessment of other MDGP algorithms. Moreover, the source code of our NDHA algorithm will be made available online, which can be used by researchers and practitioners to better solve various practical problems that can be formulated as MDGP.

(3) The neighborhood decomposition strategy is of general nature and can be advantageously adopted to speed up other neighborhood search algorithms for MDGP and other related clustering problems as well.

In the next section, we describe the proposed algorithm. In Section 3, we evaluate the proposed algorithm by reporting computational results on a large number of benchmark instances and making comparisons with reference algorithms in the literature. In Section 4, we conduct analyses to investigate two essential components of the proposed algorithm as well as one key parameter.

In the last section, we summarize the work and provide research perspectives.

2 Neighborhood decomposition based hybrid heuristic for MDGP

The neighborhood decomposition based heuristic algorithm (NDHA) proposed in this work follows the general iterated local search framework [START_REF] Lourenco | Iterated local search[END_REF] and combines two complementary neighborhood search procedures (i.e., tabu search and variable neighborhood descent) with a perturbation operator to reach a suitable tradeo between intensication and diversication of the search space. Compared to existing algorithms, the proposed algorithm distinguishes itself by two key features, i.e., its neighborhood decomposition strategy aiming to speed up the neighborhood search and a combined use of two local search methods aiming to enhance the robustness of the algorithm. We describe the main framework of the proposed algorithm and its components in this section.

Main Framework

The NDHA algorithm (Algorithm 1) starts with the solution initialization procedure (Section 2.3) to obtain a high-quality initial feasible solution. It then performs a number of iterations to improve the current solution until the given time limit (t max ) is reached (lines 426).

At each iteration, the current solution s is rst perturbed by the perturbation operator (line 5, Section 2.5) and the perturbed solution is then improved by the neighborhood decomposition tabu search (NDTS) procedure (Section 2.4.4) or the neighborhood decomposition variable neighborhood descent (ND-VND) procedure (Section 2.4.3). The decision of applying NDTS or NDVND

f (s) + δ • d(s , s) > 1) ∧ ( f (s ) f (s * ) + δ • d(s , s * ) > 1 then s ← s end if f (s ) > f (s * ) then s * ← s k ← k min end else k ← k + k step end if k ≥ k max then k ← k min end end return s * depends on a probability: Q × m N for NDTS and 1 -(Q × m N ) for NDVND, where Q (0 ≤ Q ≤ N/m
) is a parameter, N is the number of elements in the problem instance, and m is the number of groups. Subsequently, the improved solution s is accepted, like [START_REF] Brimberg | Solving the maximally diverse grouping problem by skewed general variable neighborhodd search[END_REF], as the current solution based on its objective value f (s ) and its distances to the current solution s and the best solution found so far s * (lines 1315), where the distances between solutions are measured by a partition-based distance function from [START_REF] Brimberg | Solving the maximally diverse grouping problem by skewed general variable neighborhodd search[END_REF][START_REF] Lai | Iterated maxima search for the maximally diverse grouping problem[END_REF]. After that, the best solution found so far (s * ) is accordingly updated if s is better than s * (lines 1619).

The strength k of the perturbation operator is adaptively adjusted during the search process following the strategies of breakout local search [START_REF] Benlic | Breakout local search for the vertex separator problem[END_REF][START_REF] Benlic | Breakout local search for the max-cut problem[END_REF]. k is initially set to the minimum value k min (line 3). Then, k increases by k step if the recorded best solution s * is not updated by s and is reset to k min otherwise (lines 1622). When k reaches the maximum value k max , it is reset to k min as well (lines 2325). In the following subsections, we describe the components of the NDHA algorithm. Given an edge-weighted complete graph G = (V, E, D) and m pairs of capacity limits of groups {L g , U g } (L g ≤ U g , g = 1, 2, . . . , m), the search space Ω explored by the NDHA algorithm contains all m-partitions {G 1 , G 2 , . . . , G m } of the vertex set V satisfying the capacity limits of groups, i.e.,

Search Space and Solution Representation

G 1 G 2 G 3 G 3 4 5 10 1 L 1 =2, U 1 =5 8 11 G 1 6 2 L 2 =3, U 2 =5 G 2
G 1 ∪ G 2 ∪ • • • ∪ G m = V, G i ∩ G j = ∅ (i = j)
, and L g ≤ |G g | ≤ U g for g = 1, 2, . . . , m.

To ensure a high computational eciency, we use a N -dimensional vector x[1 : N ] to represent a candidate solution (i.e., a m-partition of V ), where x[i] (i = 1, 2, . . . , N ) takes its values in {1, 2, . . . , m} and x[i] = g indicates that vertex i is clustered in group g. In addition, to ease the implementation of the neighborhood decomposition strategy, we also maintain a m × U max matrix A, as illustrated in Fig. 1, where U max denotes the largest upper limit of groups, i.e., U max = max 1≤g≤m {U g }.

Initial Solution

Following [START_REF] Lai | Iterated maxima search for the maximally diverse grouping problem[END_REF], the initialization procedure of the NDHA algorithm generates β feasible solutions (typically from ten to a few dozen) and then chooses the best one among them as the initial solution of the algorithm. Specically, the initialization procedure generates a feasible solution as follows. First, it constructs randomly a partial solution satisfying the lower capacity limits of groups, i.e., |G g | ≥ L g (g = 1, 2, . . . , m). Then, the remaining vertices are added into the partial solution one by one to obtain a complete feasible solution such that the upper capacity limits of groups are respected, i.e., |G g | ≤ U g (g = 1, 2, . . . , m). Finally, the quality of the constructed solution is locally improved by the local search procedure described in Algorithm 2.

Neighborhood Decomposition based Local Search Procedures

The NDHA algorithm relies on two key local search procedures which are used to improve solutions generated by the initialization procedure or the perturbation operator. We describe below these two local search procedures.

Neighborhood Structures and their Decomposition

The local search procedures of the NDHA algorithm are based the neighborhood decomposition strategy introduced in this work. Below, we present the two underlying neighborhoods used in this work (the constrained OneM ove neighborhood and swap neighborhood) and their decompositions. Though both neighborhoods have been used in existing studies on MDGP [4,11,18,21 23,25], the neighborhood decomposition strategy is novel, which constitutes one of the key ingredients contributing to the success of our NDHA algorithm.

The constrained OneM ove neighborhood (denoted by N 1 ) can be described as follows. Given a solution s = {G 1 , G 2 , . . . , G m }, the OneM ove move denoted by (v, i, j) transfers a vertex v from its current group G i to another group G j (j = i) such that the resulting neighbor solution (denoted by s ⊕ (v, i, j)) satises the capacity constraints of groups G i and G j . The neighborhood N 1 (s) of solution s is then composed of all feasible solutions which can be obtained by applying OneM ove to s:

N 1 (s) = {s ⊕ (v, i, j) : v ∈ G i , i, j ∈ {1, . . . , m}, i = j, |G i | > L i , |G j | < U j }.
Clearly N 1 (s) has a size bounded by O(N × m), and becomes empty if L g = U g for any group g.

We notice that N 1 (s) can be decomposed into m × (m -1) disjoint subsets that we call neighborhood blocks B 1 [i][j](s) (i, j ∈ {1, . . . , m}, i = j), where each neighborhood block

B 1 [i][j](s) is given by B 1 [i][j](s) = {s ⊕ (v, i, j) : v ∈ G i , |G i | > L i , |G j | < U j }.
As a result, the neighborhood N 1 (s) can be equivalently expressed as

N 1 (s) = ∪ 1≤i =j≤m B 1 [i][j](s).
Now we can use this decomposed neighborhood to accelerate neighborhood examination as follow. At each iteration of the neighborhood search, the current neighborhood N 1 is rst decomposed into m × (m -1) blocks. Then the algorithm skips those non-promising blocks that have been identied in pre-vious iterations and marked in a state matrix (see below), and focuses only on the remaining (promising) blocks. As a result, the neighborhood search is speeded up greatly.

To indicate whether a neighborhood block has been checked or not during the neighborhood search process, we maintain a m × m asymmetric binary state matrix M 1 , where the entry

M 1 [i][j] (1 ≤ i = j ≤ m) corresponds to the neighborhood block B 1 [i][j](s) and the diagonal entries M 1 [i][i] (1 ≤ i ≤ m) are irrelevant since the corresponding B 1 [i][i](s) (1 ≤ i ≤ m) are empty set. M 1 [i][j] takes 0 if the block B 1 [i]
[j](s) has been examined previously without nding any improving solution, and takes 1 otherwise. An illustrative example for the state matrix M 1 is shown in Fig. 2(a).

The swap neighborhood (denoted by N 2 ) is dened by the Swap(v, u) move which generates a neighbor solution by exchanging the group of vertex v and the group of vertex u. Therefore, the neighborhood N 2 (s) of s contains all possible solutions that can be reached by applying the Swap move to s:

N 2 (s) = {s ⊕ Swap(v, u) : v ∈ G i , u ∈ G j , 1 ≤ i < j ≤ m}. This neigh- borhood has a size bounded by O(N 2 ). Similar to N 1 , the neighborhood N 2 (s) can also be decomposed into m × (m - 1)/2 disjoint neighborhood blocks, i.e., N 2 (s) = ∪ 1≤i<j≤m B 2 [i][j](s), where each neighborhood block B 2 [i][j](s) is dened by B 2 [i][j](s) = {s⊕Swap(v, u) : v ∈ G i , u ∈ G j }.
Then, taking advantage of this decomposed neighborhood, we can speed up neighborhood examination during the search process by ignoring the non-promising neighborhood blocks that do not contain any improving solution. For this purpose, we use a m × m symmetric binary state matrix M 2 , where the entry

M 2 [i][j] (= M 2 [j][i]) corresponds to the block B 2 [i][j](s) and takes 0 if B 2 [i]
[j](s) has been examined previously without nding an improving solution, and takes 1 otherwise. An illustrative example of M 2 is given in Fig. 2(b).

The neighborhood decomposition method is based on the fact that the objective function of MDGP is the sum of subunit objectives dened on m groups and thus most neighborhood blocks are mutually independent in terms of move values (i.e., the change of objective values between the current solution and a neighbor solution). With the help of the state matrix of the corresponding neighborhood, the algorithm avoids many redundant neighborhood examinations if the neighborhood is checked in a block-by-block way, as described in Sections 2.4.3 and 2.4.4.

Updating of State Matrices of Neighborhoods

The state matrices M 1 and M 2 are initialized at the beginning of each neighborhood search procedure, where their entries are set to 1 except for the diago- nal elements which are set denitively to 0. Then, M 1 and M 2 are dynamically updated as the search process progresses.

Specically, for the neighborhood

N 1 , M 1 [i][j] is rst set to 0 when the block B 1 [i][j](s) is being examined. For N 2 , both M 2 [i][j] and M 2 [j][i] are set to 0 when the block B 2 [i][j](s) is being examined. After that, M 1 [i][q] (or M 2 [i][q] for N 2 ), M 1 [q][i] (or M 2 [q][i] for N 2 ), M 1 [j][q] (or M 2 [j][q] for N 2 ) and M 1 [q][j] (or M 2 [q][j] for N 2 ) (1 ≤ q ≤ m) all are updated to 1 if a solution in the block B 1 [i][j](s) (or B 2 [i][j](s) for N 2
) is chosen as the current solution, and keeps unchanged otherwise. Two illustrative examples for the update of M 1 and M 2 are given in Fig. 3, where each blue entry corresponds to the neighborhood block being examined, those entries in lilac are those that need to be updated if a neighbor solution is chosen to replace the current solution from the block being examined during the neighborhood search.

As we observe in Fig. 3, if a neighbor solution is chosen to replace the current solution during the neighborhood search, only those blocks corresponding to the lilac entries are impacted in terms of the objective value. Thus, the algorithm only needs to examine those blocks with a state value of 1. In this way, the neighborhood search process can be accelerated signicantly (see the experimental study on this issue presented in Section 4.1).

Neighborhood Decomposition based Variable Neighborhood Descent

The neighborhood decomposition based variable neighborhood descent procedure (NDVND) (Algorithm 2) relies on the general variable neighborhood descent (VND) method [START_REF] Hansen | Variable neighborhood search: Principles and applications[END_REF]. NDVND employs both neighborhoods N 1 and N 2

for i ← 1 to m do 8 for j ← 1 to m do 9 if M 1 [i][j] = 1 then M 1 [i][j] ← 0,f lag ← false for each s ∈ B 1 [i][j](s) do if f (s ) > f (s) then s ← s Improve ← true, f
j ← i + 1 to m do if M 2 [i][j] = 1 then M 2 [i][j] ← 0, M 2 [j][i] ← 0, f lag ← false for each s ∈ B 2 [i][j](s) do if f (s ) > f (s) then s ← s Improve ← true, f
0 1 0 0 0 1 0 0 1 M Updated 1 0 1 0 1 1 1 1 1 0 1 1 0 0 1 1 0 0 1 0 1 0 1 0 1 1 1 1 1 1 0 1 0 1 1 1 1 0 1 1 1 1 1 0 1 1 1 1 1 1 0 1 0 1 0 1 1 0 1 1 0 1 1 0 (a) 0 0 0 1 0 0 0 1 0 0 0 0 1 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 1 1 0 1 0 0 0 0 1 0 0 1 0 0 0 0 1 1 1 0 0 1 0 0 1 0 0 1 2 M 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 0 0 1 1 1 1 0 0 1 1 0 1 0 1 0 1 0 1 0 0 0 1 0 0 1 1 1 1 0 1 0 1 0 1 1 1 1 1 1 1 0 1 1 1 1 1 1 0 1 2 M 1 1 0 1 1 0 1 1 1 0 1 1 0 0 1 Updated (b)
Fig. 3. Illustrative example for the update of the matrices M 1 and M 2 , where an blue entry corresponds to the neighborhood block being examined, and the lilac entries are needed to be updated if a neighbor solution in the neighborhood block being examined is chosen to replace the current solution during the neighborhood search.

to explore candidate solutions. Stating from N 1 , the procedure examines N 1 and N 2 in a token-ring way until no improving solution exists in N 1 (s) and N 2 (s) with respect to the current solution s.

Specically, for a given neighborhood N i (i = 1, 2), it is examined blockby-block, and the current solution s is immediately updated each time an improving neighbor solution is found. Then, the state matrices M 1 and M 2 are updated if at least one improving neighbor solution is found after the associated neighborhood block is completely examined.

Neighborhood Decomposition based Tabu Search

The neighborhood decomposition based tabu search (NDTS) procedure of the NDHA algorithm relies on the tabu search metaheuristic [START_REF] Glover | Tabu search[END_REF] and employs a reduced swap neighborhood N N DT S obtained from the swap neighborhood N 2 . Formally, this reduced swap neighborhood is given by N N DT S (s) = {s ∈

f max ← -∞ for i ← 1 to m do for j ← i + 1 to m do 11 if (M 2 [i][j] = 1) ∨ (rand(0, 1) < µ) then 12 M 2 [i][j] ← 0, M 2 [j][i] ← 0 13 for each s ∈ B 2 [i][j](s) do 14 if ((f (s ) > f max )∧(s is not forbidden))∨(f (s ) > f (s b ) then 15 f max ← f (s ) 16 I ← i, J ← j 17 s nb ← s /
f (s) > f (s b )) then s b ← s N oImprove ← 0 end else N oImprove ← N oImprove + 1 end end return s b B 2 [i][j](s) : M 2 [i][j] = 1 ∨ rand(0, 1) < µ, i < j}
, where rand(0, 1) is a random number in the interval (0, 1) and µ is a parameter which is set to 0.05 in this study (see analysis of µ in Section 4.3). In other words, the reduced neighborhood includes always the promising neighborhood block B 2 [i][j](s) (indicated with M 2 [i][j] = 1). To avoid a too restricted (small) neighborhood at each iteration, the neighborhood block B 2 [i][j](s) of N 2 (s) is always contained in the neighborhood N N DT S (s) with a probability of µ independent of its state (M 2 [i][j] = 1 or 0).

The NDTS procedure described in Algorithm 3 starts from the initialization of the tabu list T (a N × m array) and the state matrix M 2 (line 5), then performs a number of iterations until the best solution s b can not be improved during α consecutive iterations (lines 634), where α is a parameter called the depth of tabu search.

At each iteration, the NDTS procedure examines the neighborhood N N DT S (s) in a block-by-block way and chooses a best neighbor solution (denoted by s nb = s ⊕ Swap(v, u)) that is not forbidden by the tabu list T to replace the current solution s. Then the tabu list T is accordingly updated (line 26), i.e., the corresponding vertices v and u are recorded into T and forbidden to move back to their previous groups for the next tl iterations, where tl = 15 + rand( 5) is the tabu tenure with rand(5) being a random integer between 0 and 4. Moreover, the aspiration criterion is applied, i.e., a neighbor solution s nb always replaces the current solution if the quality of s nb is better than the best solution found so far (s b ) (line 14). After that, the state matrix M 2 is accordingly updated (line 26). To diversify the search process and jump out of local optimum traps, the NDHA algorithm employs a perturbation operator to modify the solutions returned by the local search methods. Specically, starting from the input solution s 0 , the perturbation operator performs a number λ of random swap moves to generate a new solution, where each swap move exchanges the groups of two random vertices located in two distinct groups. As indicated in Section 2.1, the strength of perturbation λ is a changing value which is dynamically adjusted by an adaptive technique. The pseudo-code of the perturbation operator is given in Algorithm 4.

Discussions on the Innovations of the Work

Compared with the existing iterative algorithms in the literature, such as SGVNS [START_REF] Brimberg | Solving the maximally diverse grouping problem by skewed general variable neighborhodd search[END_REF], ITS [START_REF] Palubeckis | Maximally diverse grouping: an iterated tabu search approach[END_REF] and IMS [START_REF] Lai | Iterated maxima search for the maximally diverse grouping problem[END_REF], the proposed NDHA algorithm has two following original features. First, unlike the existing methods where only one local search procedure is used, the proposed NDHA algorithm combines in a probabilistic way two complementary local search procedures (i.e., tabu search and variable neighborhood descent) to reinforce its search robustness on the dierent types of problem instances. Second, both our tabu search and variable neighborhood descent rely on the innovative neighborhood decomposition strategy to speed up the neighborhood evaluation process, which signicantly improves the computational and search eciency of the algorithm (as shown in Section 3).

Actually, the proposed neighborhood decomposition strategy is the most important innovation for this work, and it is able to speed greatly up the neighborhood examination process, without missing improving candidate solutions within the given neighborhood. The neighborhood decomposition strategy follows the general idea of the candidate list strategy to reduce the examined neighborhood size while retaining high-quality solutions. Importantly, it provides a practical and highly eective technique to organize the neighbor solutions into neighborhood blocks and to enable the search procedure to focus on promising neighbor solutions without compromising the quality of the search.

Finally, the idea of neighborhood decomposition is very general and can be advantageously adopted in neighborhood search algorithms for other grouping or clustering problems. Thus, the contribution introduced in this work goes beyond the problem considered here and could potentially benet many heuristic algorithms for dicult combinatorial optimization.

the-art algorithms in the literature.

Benchmark Instances

In our experiments, we test 320 benchmark instances widely used in the literature. These instances belong to ve sets whose main characteristics are summarized as follows 1 .

• RanReal set (20 instances): This set includes 10 instances with dierent group sizes (DGS) and 10 instances with equal group sizes (EGS). For these instances, the number N of elements equals 480 or 960, the number m of groups equals 20 or 24, L g and U g vary between 10 and 50, and the distances d ij (i < j) are a real number generated randomly in the interval (0, 100).

For the EGS instances, both L g and U g are equal to N/m for any group g. These instances were tested in [START_REF] Brimberg | Solving the maximally diverse grouping problem by skewed general variable neighborhodd search[END_REF][START_REF] Gallego | Tabu search with strategic oscillation for the maximally diverse grouping problem[END_REF][START_REF] Lai | Iterated maxima search for the maximally diverse grouping problem[END_REF][START_REF] Palubeckis | Maximally diverse grouping: an iterated tabu search approach[END_REF][START_REF] Rodriguez | An articial bee colony algorithm for the maximally diverse grouping problem[END_REF]. • RanInt set (20 instances): Similar to RanReal set, this set contains 10 DGS instances and 10 EGS instances. The main characteristics of these instances are the same as RanReal instances, while the distances d ij (i < j) between elements are an integer generated randomly between 0 and 100. These instances were tested in [START_REF] Brimberg | Solving the maximally diverse grouping problem by skewed general variable neighborhodd search[END_REF][START_REF] Gallego | Tabu search with strategic oscillation for the maximally diverse grouping problem[END_REF][START_REF] Lai | Iterated maxima search for the maximally diverse grouping problem[END_REF][START_REF] Palubeckis | Maximally diverse grouping: an iterated tabu search approach[END_REF][START_REF] Rodriguez | An articial bee colony algorithm for the maximally diverse grouping problem[END_REF]. • Geo set (20 instances): Similar to RanReal and RanInt sets, this set contains 10 EGS instances and 10 DGS instances whose main characteristics are the same as RanReal and RanInt instances, while the distances d ij (i < j) between elements are Euclidean distances between pairs of points with random coordinates from [0, 10], and the number of coordinates of points varies from 2 to 21. These instances were tested in [START_REF] Brimberg | Solving the maximally diverse grouping problem by skewed general variable neighborhodd search[END_REF][START_REF] Gallego | Tabu search with strategic oscillation for the maximally diverse grouping problem[END_REF][START_REF] Lai | Iterated maxima search for the maximally diverse grouping problem[END_REF][START_REF] Palubeckis | Maximally diverse grouping: an iterated tabu search approach[END_REF][START_REF] Rodriguez | An articial bee colony algorithm for the maximally diverse grouping problem[END_REF]. • MDG-a set (220 instances): This set is composed of 11 subsets, including 6 subsets of DGS instances and 5 subsets of EGS instances, and each subset consists of 20 instances that were generated from 20 edge-weighted complete graphs with N = 2000, where the edge weights d ij (i < j) are an integer generated randomly between 0 and 10. The main characteristics of these 11 subsets are summarized in Table 1. These instances were tested in [START_REF] Brimberg | Solving the maximally diverse grouping problem by skewed general variable neighborhodd search[END_REF][START_REF] Lai | Iterated maxima search for the maximally diverse grouping problem[END_REF][START_REF] Rodriguez | An articial bee colony algorithm for the maximally diverse grouping problem[END_REF]]. • MDG-c set (40 instances): This set is composed of 20 DGS instances and 20 EGS instances with n = 3000 and m = 50, where L g and U g are respectively set to 0.8N/m and 1.2N/m for the DGS instances, and N/m for the EGS instances. The distances d ij (i < j) between elements are an integer generated randomly between 0 and 1000. These instances are the largest instances used in this study and were rst used in [START_REF] Lai | Iterated maxima search for the maximally diverse grouping problem[END_REF]. The NDHA algorithm adopts several parameters whose descriptions and settings are listed in Table 2. The value of each parameter was xed independently according to a preliminary experiment performed on a selection of instances of dierent characteristics. Typically, we tested a number of possible values from a given range to retain the value leading to the best average result. We observe that among the parameters, Q which controls the probability of applying the NDVND or NDTS optimization procedure is one key parameter, for which we provide a detailed analysis in Section 4.2. Notice that the parameter values shown in Table 2 can be considered to be the default parameter setting of the NDHA algorithm, which were used consistently to perform all the experiments reported in this work unless stated otherwise.

According to the computational results reported in [START_REF] Brimberg | Solving the maximally diverse grouping problem by skewed general variable neighborhodd search[END_REF][START_REF] Lai | Iterated maxima search for the maximally diverse grouping problem[END_REF][START_REF] Palubeckis | Maximally diverse grouping: an iterated tabu search approach[END_REF][START_REF] Singh | A new hybrid genetic algorithm for the maximally diverse grouping problem[END_REF], the algorithms ITS [START_REF] Palubeckis | Maximally diverse grouping: an iterated tabu search approach[END_REF], SGVNS [START_REF] Brimberg | Solving the maximally diverse grouping problem by skewed general variable neighborhodd search[END_REF], IMS [START_REF] Lai | Iterated maxima search for the maximally diverse grouping problem[END_REF], and NSGGA [START_REF] Singh | A new hybrid genetic algorithm for the maximally diverse grouping problem[END_REF] (only for instances with equal group sizes) outperform signicantly other algorithms in the literature and can be regarded as the state-of-the-art algorithms for MDGP. Hence, in this work, we use these algorithms as the reference algorithms. Among these reference algorithms, the source code of ITS is available at http://www.proin.ktu.lt/~gintaras/mdgp.html, the source code of IMS is available at http://www.info.univ-angers.fr/pub/hao/mdgp.html, and the executable code of SGVNS was kindly provided by the authors of [START_REF] Brimberg | Solving the maximally diverse grouping problem by skewed general variable neighborhodd search[END_REF]. The proposed NDHA algorithm as well as IMS and ITS were written in C++ and compiled using the g++ compiler with the -O3 option.

In addition, all the computational experiments were carried out on the same computing platform with an Intel E5-2670 processor (2.5 GHz and 2G RAM), running the Linux operating system. Following [START_REF] Lai | Iterated maxima search for the maximally diverse grouping problem[END_REF], the stopping condition of all the algorithms is a cuto time limit t max set to 120, 600, 1200, and 3000 seconds for instances n = 480, n = 960, n = 2000 and n = 3000, respectively. Finally, to assess the average performance of the algorithms, the NDHA algorithm and the reference algorithms (i.e., ITS, SGVNS, IMS) were performed 20 times with dierent random seeds for each run.

Computational Results and Comparison on the Small Instances

The rst experiment is devoted to an assessment of the NDHA algorithm on the 60 small instances with N ≤ 960, from RanInt, RanReal and Geo sets.

The experimental results of the reference algorithms (ITS, SGVNS, IMS) as well as our NDHA algorithm are summarized in In addition, to check the statistical dierence between NDHA and each reference algorithm in terms of f best or f avg , the p-values from the Wilcoxon signed-rank tests are reported in the last row of the tables, and a p-value smaller than 0.05 means that there exists a signicant dierence between the NDHA algorithm and the corresponding reference algorithm.

Table 3 shows that for the instances with N ≤ 960 and dierent group sizes, the NDHA algorithm is very competitive compared to the reference algorithms. In terms of f best , the ITS, SGVNS, IMS and NDHA algorithms obtained respectively the best result for 3, 0, 15 and 12 out of 30 instances. In terms of f avg , these four algorithms obtained the best result for 1, 0, 15 and 14 instances. Moreover, the small p-values (≤ 0.05) mean that the NDHA algorithm outperforms signicantly the ITS and SGVNS algorithms both in terms of f best and f avg . When comparing with IMS and NDHA, we observe that they perform similarly (p-values ≥ 0.05) in terms of f best and f avg .

Table 4 indicates that for the instances with N ≤ 960 and equal group sizes, the NDHA algorithm performs well compared to the reference algorithms. In terms of f best , the ITS, SGVNS, IMS and NDHA algorithms obtained respectively the best result for 11, 0, 8 and 11 out of 30 instances. The large p-values (7.19E-2 and 4.53E-1) imply that there does not exist a signicant dierence between NDHA and ITS (or IMS). However, the small p-value (2.76E-3) means that NDHA performs signicantly better than the SGVNS algorithm. In terms of f avg , the p-values show that the NDHA algorithm outperforms signicantly the ITS and SGVNS algorithms and performs similarly compared to the IMS algorithm.

Computational Results and Comparison on the Large Scale Instances

The second experiment aims to assess the NDHA algorithm on the 260 large scale instances belonging to 11 subsets of the set MDG-a and 2 subsets of the set MDG-c, where each subset has 20 instances with N = 2000 or N = 3000. For the experiment, the NDHA algorithm and the reference algorithms (ITS, SDVNS, and IMS) were respectively run 20 times on each instance, and the detailed results for each subset are reported in Tables A.1A.13 of the appendix, where the same information as in Tables 3 and4 is provided.

A summary of these detailed experimental results is provided in Table 5 shows that the NDHA algorithm performs very well and outperforms signicantly the reference algorithms. The NDHA algorithm obtained the best results for 222 and 220 out of 260 instances in terms of f best and f avg , respectively. Compared to the ITS and SGVNS algorithms, our algorithm obtained a better result for each instance both in terms of f best and f avg . When comparing with the IMS algorithm, the NDHA algorithm yielded a better result both in terms of f best and f avg on most instances except the instances with m = 10, which means that the neighborhood decomposition strategy used in the NDHA algorithm plays an important role in enhancing the performance of the algorithm for the instances with a large number of groups. This experiment indicates that the proposed NDHA algorithm is highly competitive compared to the three reference algorithms especially for the large instances with a large number of elements and groups.

In addition, to make a comparison between the NDHA algorithm and the recent hybrid genetic algorithm NSGGA which is designed for MDGP with equal group sizes [START_REF] Singh | A new hybrid genetic algorithm for the maximally diverse grouping problem[END_REF], we provide the results of the NSGGA and NDHA algorithms in Table 6 for 20 large EGS instances with N = 2000 and m = 200, where the same statistical information is given as in the previous tables. Since the code of the NSGGA algorithm is not available to us, this comparison is based on the results reported in [START_REF] Singh | A new hybrid genetic algorithm for the maximally diverse grouping problem[END_REF], based 20 runs per instance on an Intel Core i5 computer with 4G RAM under the same stopping condition as that used by our NDHA algorithm (specied in Section 3.2). One observes from Table 6 that the NDHA algorithm signicantly dominates the NSGGA algorithm in terms of f best and f avg , which is conrmed by small p-values.

In summary, the experimental results shown in this section indicate that the NDHA algorithm is very competitive compared with the state-of-the-art algorithms in the literature and performs especially well on the large instances.

Analysis and Discussions

In this section, we analyze two essential components of the proposed algorithm, i.e., the neighborhood decomposition strategy employed by the local search methods and the strategy of jointly using two local search methods.

Importance of the Neighborhood Decomposition Strategy

In order to show the eectiveness of the neighborhood decomposition strategy introduced in this study, we carried out an experiment based on 20 large EGS instances with n = 2000, m = 200 and L g = U g = 10. For this study, we created a variant (denoted by NDHA-D) of the NDHA algorithm by disabling the neighborhood decomposition strategy (i.e., setting all entries of state matrices M 1 and M 2 to the value of 1 during the neighborhood search process), while keeping other components unchanged. We ran NDHA-D and NDHA 20 times for each instance according to the experimental protocol of Section 3.2.

The experimental results are summarized in Table 7, where columns 14 give the name and main characteristics of the instances, columns 56 show respectively the best objective value (f best ) over 20 runs for the two compared algorithms, columns 78 indicate the average objective value f avg , and columns 910 present the worst objective value (f worst ). In addition, the row 'Avg.' shows the average result for each column, the row '#better' shows the number of instances for which the associated algorithm obtained a better result than the competing algorithm in terms of f best , f avg , and f worst , respectively, and the last row gives the p-values from the Wilcoxon signed-rank tests for the compared algorithms. The better results between the compared algorithms are indicated in bold.

One observes from Table 7 that the NDHA algorithm dominates the NDHA-D algorithm for each considered performance indicator. Specically, the NDHA algorithm obtained a better result than the NDHA-D algorithm for all tested instances in terms of f best , f avg , and f worst . Moreover, the small p-values imply the dierence between the two compared algorithms is statistically signicant. This experiment shows the neighborhood decomposition strategy used in this study plays an important role for the high performance of the NDHA algorithm.

To further show the eect of the neighborhood decomposition strategy on the tabu search procedure that is one main component of the NDHA algorithm, we carried out another experiment to compare the tabu search procedure with neighborhood decomposition (i.e., NDTS) and a tabu search procedure without neighborhood decomposition (denoted by TS). To ease the presentation, The subgures (a) and (c) of Fig. 4 indicate that NDTS (with neighborhood decomposition) requires much less time than TS (without neighborhood decomposition) to perform the same number of iterations, which means that the neighborhood decomposition technique is able to speed up the neighborhood search process notably. On the other hand, one observes from the subgures (b) and (d) that NDTS yielded much better results in the objective value f (s) than TS within the same computation time, implying the neighborhood decomposition technique is able to enhance signicantly the performance of the tabu search procedure.

Impact of Hybridizing Two Local Search Methods

To enhance its robustness for solving instances with very dierent characteristics, the NDHA algorithm hybridizes two local search procedures (i.e., NDVND and DNTS) in a probabilistic way, where the probability is con- trolled by the parameter Q. To show the rationality of this hybridization and choose an appropriate value for Q, we carried out an additional experiment based on ve representative EGS instances with very dierent numbers (m) of groups. For each considered instance and each Q value in the range {0.0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, N/m}, the NDHA algorithm was performed 20 times using the experimental protocol in Section 3.2. The experimental results are summarized in Fig. 5 using the box and whisker plots, where the X-axis indicates the values of Q and the Y-axis gives the objective values.

Fig. 5 shows that the performance of the algorithm is sensitive to the setting of Q. Specically, when the neighborhood decomposition based tabu search (NDTS) procedure was always applied and the neighborhood decomposition based variable neighborhood descent (NDVND) procedure was disabled (with Q = N/m), the NDHA algorithm yielded the worst results among the considered Q values for the instances with a small number (m ≤ 50) of groups.

On the other hand, when only the NDVND procedure was applied and the NDTS procedure was disabled (with Q = 0), the NDHA algorithm produced the worst results among the considered Q values for the instances with a large number (m ≥ 100) of groups. The experiment shows that the NDTS procedure is more suitable for the instances with a large number of groups, while the NDVND procedure is more suitable for the instances with a small number of groups. This provides the main motivation of combining NDVND and DNTS in a probabilistic way to be able to deal with instances with small and large number of groups. One notices that the setting of Q = 0.1 led globally to a good performance for the considered instances, which was used as the default value of Q in this work. Fig. 6 indicates that the setting of parameter µ signicantly impacts the be-havior of the NDHA algorithm. Specically, For µ < 0.05, the performance of the algorithm gradually increases with the increase of µ, while the performance decreases as the value of µ increases for µ ≥ 0.15. Thus, [0.05, 0.10] is a suitable range for µ and 0.05 was adopted as the default value in this work.

5 Conclusions and Future Work

We presented a new heuristic algorithm (NDHA) for solving the maximally diverse grouping problem (MDGP). The proposed algorithm distinguishes itself from existing algorithms by its speeding-up neighborhood decomposition technique and the joint use of two complementary local search procedures (both are based on neighborhood decomposition). An adaptive perturbation strategy is additionally used to escape local optimum traps. The algorithm was assessed on 320 benchmark instances commonly used in the literature. Our computational results show that the algorithm outperforms signicantly the state-of-the-art algorithms by reporting improved lower bounds for 220 large benchmark instances (i.e., for more than 68% of the tested instances).

Given that MDGP is a general model able to formulate a number of realworld applications, the proposed algorithm provides a valuable tool to better solve the related practical problems. The availability of the source code of our algorithm further facilitates such applications.

Additional analyses indicate that the eectiveness of the algorithm is mainly attributed to two essential ideas, i.e., the neighborhood decomposition strategy which speeds up the neighborhood search process and the stochastic combination strategy of two local search procedures which enhances the robustness of the algorithm for dierent types of problem instances.

The ideas of neighborhood decomposition and combination of dierent local search procedures are rather general. It is interesting to verify their usefulness for solving other related grouping or clustering problems, such as the capacitated p-median problem [START_REF] Díaz | Hybrid scatter search and path relinking for the capacitated p-median problem[END_REF], the normalized cut clustering problem [START_REF] Dhillon | Weighted graph cuts without eigenvectors: a multilevel approach[END_REF], and the balanced k-means clustering problem [START_REF] Costa | Less is more: basic variable neighborhood search heuristic for balanced minimum sum-of-squares clustering[END_REF]. Moreover, the proposed algorithm can be further reinforced by following two directions. First, the current algorithm does not explicitly deal with the issue of solution symmetry for the case where some or all groups have an equal group size. It is meaningful to investigate dedicated search strategies and techniques exploiting the symmetry property. Second, population-based methods such as memetic computing are known to be general frameworks with the potential of surpassing local search algorithms [START_REF] Hao | Memetic algorithms in discrete optimization[END_REF]. Then it is worth designing such hybrid algorithms where the proposed algorithm plays the key role of local optimization for intensication.

[27] Weitz R., Lakshminarayan S., 1998, An empirical comparison of heuristic methods for creating maximally diverse groups. Journal of the Operational Research Society, 49, 635646.

[28] Yeoh H.K., Nor M.I.M, 2011, An algorithm to form balanced and diverse groups of students. Computer Applications in Engineering Education, 19(3), 582590.

[29] Zhou Y., Hao J.K., Goëon A., 2016, A three-phased local search approach for the Clique Partitioning problem. Journal of Combinatorial Optimization, 32(2), 469491.

A Appendix

Detailed results of the proposed NDHA algorithm and three main reference algorithms (ITS [START_REF] Palubeckis | Maximally diverse grouping: an iterated tabu search approach[END_REF], SGVNS [START_REF] Brimberg | Solving the maximally diverse grouping problem by skewed general variable neighborhodd search[END_REF], and IMS [START_REF] Lai | Iterated maxima search for the maximally diverse grouping problem[END_REF]) on 260 large scale instances with N = 2000 or 3000 are summarized in Tables A.1A.13, where each table corresponds to a subset of benchmarks and the statistical information is the same as in the tables of Section 3. The dominating values are indicated in bold.

One observes from the tables that the NDHA algorithm outperforms significantly the reference algorithms for most instances. Specically, for the 40 instances with m = 10, NDHA outperforms signicantly ITS and SGVNS, but performs worse than IMS. For the remaining 220 instances with m ≥ 25, NDHA outperforms consistently the reference algorithms on all the instances both in terms of f best and f avg .

These outcomes show that the proposed algorithm is particularly suitable to solve large scale instances with a high number of groups, which can be attributed to the neighborhood decomposition strategy.

Fig. 1 .

 1 Fig. 1. An illustrative example for solution representation with N = 11, m = 3, and U max = 6. We use a 3 × 6 matrix A (right gure) to indicate a solution (left gure).

N 2 Fig. 2 .

 22 Fig. 2. State matrices of the constrained OneM ove neighborhood N 1 and the swap neighborhood N 2 , where M 1 and M 2 both are a m × m (m = 9) 01 matrix. For the current solution s, the entries M 1 [i][j] and M 2 [i][j] correspond respectively to the neighborhood blocks B 1 [i][j](s) and B 2 [i][j](s), which take 0 if the corresponding block has been checked without nding an improving solution.

  lag ← true end end if f lag = true then Update M 1 and M 2 /* Section 2.4.

2. 5

 5 Perturbation Operator Algorithm 4: Perturbation Operator Function Perturbation(s 0 , η) Input: Input solution s 0 , strength of perturbation λ Output: The perturbed solution s s ← s 0 for l ← 1 to λ do Randomly select two vertices v and u locating at dierent groups in s s ← s ⊕ Swap(v, u) /* Swap the groups of vertices v and u */ end return s

  experiment with the results obtained by NDTS and TS on two large instances (i.e., MDG-a_21 with N=2000, m = 200 and L g = U g = 10 for any g, and MDG-a_24 with N = 2000, m = 50, L g = 26, and U g = 54 for any g). For the experiment, NDTS and TS were performed one time per instance, starting from the same initial solution. The experimental results are shown in Fig.4, where the subgures (a) and (c) show the running times of NDTS and TS as a function of the number of iterations, and the subgures (b) and (d) show their best objective values f (s) found as a function of running time.

Fig. 4 .

 4 Fig. 4. Comparison between NDTS and TS in terms of the objective function value f (s) and the running time.

-a_21 with m = 200 Fig. 5 .

 2005 Fig.5. Sensitivity analysis for the parameter Q used to control the probability that the NDTS procedure is applied. The X-axis indicates the settings of Q and the Y-axis indicates the objective value.

4. 3 Fig. 6 .

 36 Fig. 6. Sensitivity analysis of parameter µ based on the 20 EGS instances with N = 2000 and m = 200, where the X-axis indicates the parameter values and the Y-axis shows the average objective values (f avg ) obtained over 20 independent runs.

  number of groups m, time limit t max , and parameters δ, Q, k min , k step , k max Output: The best feasible solution found (s * )

	Algorithm 1: Neighborhood decomposition based heuristic algorithm
	(NDHA) for MDGP
	s ← InitialSolution(G, m)	/* Section 2.3 */
	s * ← s	/* s * records the best solution found */
	k ← k min	/* k denotes the current perturbation strength */
	while Time() ≤ t max do
	s ← P erturbation(s, k)	/* Algorithm 4 */
	r ← rand(0, 1) /* rand(0, 1) denotes a random number in (0, 1)
	*/
	if r < Q × m N then s ← N DT S(s )	/* Tabu search, Algorithm 3 */
	end
	else
	s ← N DV N D(s )	/* Variable neighborhood descent,
	Algorithm 2 */
	end
	if (

Input: An edge-weighted complete graph G = (V, E, D), f (s )

Table 1

 1 Main characteristics of the instances in the set MDG-a, where '#' indicates the number of instances in the corresponding subset.

			DGS			EGS
	N	m	Lg	Ug	#	Lg = Ug	#
	2000	50	32	48	20
	2000	10	173	227	20	200	20
	2000	25	51	109	20	80	20
	2000	50	26	54	20	40	20
	2000	100	13	27	20	20	20
	2000	200	6	14	20	10	20
	3.2 Parameter Setting and Experimental Protocol

Table 2

 2

	Setting of important parameters
	Parameters Section	Description	Values
	δ	2.1	parameter used in acceptance criterion	0.01
	Q	2.1	probability of applying NDTS and NDVND	0.1
	k min	2.1	minimum perturbation strength	0.2 × N/m
	kmax	2.1	maximum perturbation strength	2 × N/m
	kstep	2.1	incremental value of perturbation strength	0.2 × N/m
	α	2.4.4	depth of tabu search	500
	µ	2.4.4	parameter used in the neighborhood N N DT S	0.05

Table 3 (

 3 with dierent group sizes) andTable 4 (with equal group sizes). Columns 25 of each table give the best objective value (f best ) over 20 runs respectively for the compared algorithms, and columns 69 show the average objective value (f avg ). The row 'Avg.' indicates the average result for each column, and the row '#best' indicates the number of instances for which the corresponding algorithm obtained the best result among the compared algorithms in terms of f best and f avg . The best f best and f avg values among the compared algorithms are indicated in bold for each instance.

Table 3

 3 Comparison of the proposed NDHA algorithm with three state-of-the-art algorithms in the literature on the small instances with dierent group sizes. The best results between the compared algorithms in terms of f best and f avg are indicated in bold.

			f best				favg
	Instance	ITS	SGVNS	IMS	NDHA	ITS	SGVNS	IMS	NDHA

Table 5 ,

 5 where each row represents one subset. Columns 15 of the table give the main characteristics of the instances in the corresponding subset, columns 69 indicate the average results (i.e., the Avg. value in Tables A.1A.13) in terms of f best respectively for each subset and each algorithm, and columns 1013

Table 4

 4 Comparison of the proposed NDHA algorithm with three state-of-the-art algorithms in the literature on the small instances with equal group sizes. The best results between the compared algorithms in terms of f best and f avg are indicated in bold.

	f best	favg
	Instance Geo_n480_ss_01 Geo_n480_ss_02 Geo_n480_ss_03 Geo_n480_ss_04 Geo_n480_ss_05 RanInt_n480_ss_01 379408.00 379532.00 380127.00 377481.62 ITS SGVNS IMS NDHA 552206.89 552040.71 552073.45 552194.91 1047462.35 1047245.92 1047228.47 1047433.23 633855.88 633574.99 633626.24 633740.92 789891.16 789621.27 789657.63 789767.10 945974.03 945667.90 945782.66 945865.98	ITS 552165.47 552020.33 552045.83 552167.92 SGVNS IMS NDHA 1047405.08 1047155.51 1047182.54 1047322.44 633746.62 633544.22 633590.50 633713.2514 789791.95 789544.51 789610.72 789731.7703 945898.63 945643.47 945697.53 945819.7054

Table 5

 5 Summary comparison of the NDHA algorithm with three state-of-the-art algorithms on 260 large instances with n = 2000 (220 instances) and n = 3000 (40 instances).The best results between the compared algorithms in terms of f best and f avg are indicated in bold. results in terms of f avg for each subset and each algorithm. The last row of the table indicates the number of instances for which the associated algorithm obtained the best result in terms of f best and f avg among the compared algorithms.

	Set	Instance sets m Lg Ug T ype	ITS	f best SGVNS	IMS	NDHA	ITS	SGVNS	favg	IMS	NDHA
						0		0	40		220
	shows the average

Table 6

 6 Comparison of the NDHA algorithm with the NSGGA algorithm on 20 EGS instances with N = 2000 and m = 200. The best results between the compared algorithms in terms of f best and f avg are indicated in bold.

	Instances m Lg Ug MDG-a_21 200 10 10 Graph MDG-a_22 200 10 10 MDG-a_23 200 10 10 MDG-a_24 200 10 10 MDG-a_25 200 10 10 MDG-a_26 200 10 10 MDG-a_27 200 10 10 MDG-a_28 200 10 10 MDG-a_29 200 10 10 MDG-a_30 200 10 10 MDG-a_31 200 10 10 MDG-a_32 200 10 10 MDG-a_33 200 10 10 MDG-a_34 200 10 10 MDG-a_35 200 10 10 MDG-a_36 200 10 10 MDG-a_37 200 10 10 MDG-a_38 200 10 10 MDG-a_39 200 10 10 MDG-a_40 200 10 10 Avg. #Better p-value	NSGGA 77610 77671 77567 77401 77536 77442 77510 77670 77442 77575 77557 77470 77480 77538 77684 77519 77630 77493 77461 77544 77540.00 0 8.86E-5	f best	NDHA 78193 78423 78253 78300 78266 78324 78220 78208 78271 78187 78380 78252 78234 78193 78332 78348 78335 78189 78290 78340 78276.90 20	NSGGA 77299.80 77290.50 77271.30 77213.35 77317.30 77263.25 77241.30 77290.85 77242.60 77272.60 77323.55 77245.05 77271.85 77327.65 77334.60 77273.10 77373.85 77269.00 77228.35 77321.50 77283.57 0 8.86E-5	favg NDHA 78101.00 78098.35 78111.00 78075.35 78143.55 78107.90 78085.00 78107.75 78104.90 78092.05 78255.45 78117.95 78085.05 78082.05 78094.25 78158.55 78126.90 78100.00 78122.50 78155.55 78116.26 20

Table 7

 7 Comparison between the NDHA algorithm and its variant NDHA-D on 20 EGS instances with N = 2000 and m = 200.

Table A .

 A 1 Comparison of the proposed NDHA algorithm with three best performing algorithms in the literature on the 20 DGS instances with N = 2000 and m = 10. Table A.2 Comparison of the proposed NDHA algorithm with three best performing algorithms in the literature on the 20 DGS instances with N = 2000 and m = 25.

	Graph	Instance m Lg Ug	ITS	f best SGVNS	IMS	NDHA	ITS	SGVNS	favg	IMS	NDHA
	p-value		8.86E-5 8.86E-5 1.20E-4	1	0 8.86E-5 8.86E-5 8.86E-5 0 20	0
	Graph	Instance m Lg Ug	ITS	f best SGVNS	IMS	NDHA	ITS	SGVNS	favg	IMS	NDHA
	p-value		8.86E-5 8.86E-5 8.86E-5		0 8.86E-5 8.86E-5 8.86E-5 0 0	20

Table A .

 A [START_REF] Brimberg | Solving the maximally diverse grouping problem by skewed general variable neighborhodd search[END_REF] Comparison of the proposed NDHA algorithm with three best performing algorithms in the literature on the 20 DGS instances with N = 2000, m = 50, L g = 32, and U g = 48.

	Graph	Instance m Lg Ug	ITS	f best SGVNS	IMS	NDHA	ITS	SGVNS	favg	IMS	NDHA

Table A .

 A 7 Comparison of the proposed NDHA algorithm with three best performing algorithms in the literature on the 20 EGS instances with N = 2000 and m = 10.

	Graph	Instance m Lg Ug	ITS	f best SGVNS	IMS	NDHA	ITS	SGVNS	favg	IMS	NDHA
	p-value		8.86E-5 8.86E-5 8.86E-5		0 8.86E-5 8.86E-5 8.86E-5 0 0	20

Computational Experiments and AssessmentsIn this section, we assess the NDHA algorithm by performing large experiments on 320 benchmark instances and making comparisons with state-of-

The benchmark instances as well as the source code of our algorithm will be available at http://www.info.univ-angers.fr/pub/hao/NDHA.html

Acknowledgments

We are grateful to the reviewers for their valuable comments and suggestions which helped us to improve the paper. This work was partially supported by the National Natural Science Foundation of China (Grant No. 61703213), the Natural Science Foundation of Jiangsu Province of China (Grant No. BK20170904), six talent peaks project in Jiangsu Province (Grant No. RJFW-011), and Shenzhen Science and Technology Innovation Commission under Grant Nos. JCYJ20180508162601910 and 2019-INT003.