Solving the Latin Square Completion Problem by Memetic Graph Coloring

The Latin square completion (LSC) problem involves completing a partially filled Latin square of order ${n}$ by assigning numbers from 1 to ${n}$ to the empty grids such that each number occurs exactly once in each row and each column. LSC has numerous applications and is, however, NP-complete. In this paper, we investigate an approach for solving LSC by converting an LSC instance to a domain-constrained Latin square graph and then solving the associated list coloring problem. To be effective, we first employ a constraint propagation-based kernelization technique to reduce the graph model and then call for a dedicated memetic algorithm to find a legal list coloring. The population-based memetic algorithm combines a problem-specific crossover operator to generate meaningful offspring solutions, an iterated tabu search procedure to improve the offspring solutions, and a distance-quality-based pool updating strategy to maintain a healthy diversity of the population. Extensive experiments on more than 1800 LSC benchmark instances in the literature show that the proposed approach can successfully solve all the instances, surpassing the state-of-the-art methods. To our knowledge, this is the first approach achieving such a performance for the considered problem. We also report computational results for the related partial Latin square extension problem.


Solving the Latin Square Completion Problem by Memetic Graph Coloring Yan Jin and Jin-Kao Hao
Abstract-The Latin square completion (LSC) problem involves completing a partially filled Latin square of order n by assigning numbers from 1 to n to the empty grids such that each number occurs exactly once in each row and each column.LSC has numerous applications and is, however, NP-complete.In this paper, we investigate an approach for solving LSC by converting an LSC instance to a domain-constrained Latin square graph and then solving the associated list coloring problem.To be effective, we first employ a constraint propagation-based kernelization technique to reduce the graph model and then call for a dedicated memetic algorithm to find a legal list coloring.The population-based memetic algorithm combines a problem-specific crossover operator to generate meaningful offspring solutions, an iterated tabu search procedure to improve the offspring solutions, and a distance-quality-based pool updating strategy to maintain a healthy diversity of the population.Extensive experiments on more than 1800 LSC benchmark instances in the literature show that the proposed approach can successfully solve all the instances, surpassing the state-of-the-art methods.To our knowledge, this is the first approach achieving such a performance for the considered problem.We also report computational results for the related partial Latin square extension problem.Index Terms-Graph coloring, Latin square completion (LSC), list coloring, memetic search, tabu search (TS).

I. INTRODUCTION
A LATIN square L of order n is composed of n × n grids (or cells) such that each grid is filled with a number in {1, . . ., n} (n ∈ N + ) and each number occurs in each row and each column exactly once.If some grids of L remain unfilled (or empty), L is a partial Latin square.The Latin square completion (LSC) problem of order n involves completing the empty grids of a partial Latin square with numbers in {1, . . ., n} to obtain a legal Latin square.
LSC was first studied by Hall [19] and Ryser [37], and was known to be NP-complete in the general case [1], [8], [11].LSC can be considered as a special case of the partial Latin square extension (PLSE) problem, which is to assign numbers in {1, . . ., n} to as many empty grids as possible under the condition that each number has to occur at most once in each row and each column.Both LSC and PLSE arise naturally in a variety of practical applications, such as scheduling, optical routing, error correcting codes as well as combinatorial design [3], [9], [16], [26], [29].
Given their theoretical and practical importance, a number of studies on LSC and PLSE have been reported in the literature.For instance, in 1999, Kumar et al. [26] proposed two approximation algorithms for PLSE with nontrivial worst-case performance guarantees.In 2002, Gomes and Shmoys [17] studied three complete solution methods for solving LSC: 1) a constraint satisfaction-based approach (CSP); 2) a hybrid 0/1 linear programming/CSP-based strategy (LP/CSP); and 3) a Boolean satisfiability-based approach (SAT).In 2004, Ansótegui et al. [2] focused on a systematic comparison of SAT and CSP models for the Latin square (quasigroup) completion problem.In 2004, Gomes et al. [18] presented a natural randomized rounding algorithm based on a packing linear programming relaxation, which yields an e/(e−1)-approximation algorithm.These algorithms are able to solve small LSC instances within a reasonable time, but fail to solve most of the large and hard instances.Recently in 2016, Haraguchi [21] introduced several powerful iterated local search algorithms with multiple neighborhoods to solve PLSE as well as LSC.Assessed on a large set of 1800 instances of various sizes and characteristics, these local search algorithms showed state-ofthe-art performances.In particular, the Trellis-neighborhood search algorithm (Tr-ILS * ) proves to outperform other tested ILS variants and two general optimization solvers (IBM-ILOG CPLEX and LocalSolver).The instances and the associated results presented in [21] will be used as the main references for our computational studies.
Despite much research effort dedicated to LSC and the resulting advances, there are still very few methods that are able to solve the problem effectively.For instance, no existing algorithm can find a solution for some traditional instances tested in [17] and many new instances introduced by Haraguchi [21].It is thus quite useful and challenging to devise a method able to solve large and difficult instances.
In this paper, we investigate for the first time a solution method for solving LSC by converting the problem to a particular graph coloring problem (i.e., precoloring extension [5], then list coloring [12], [28]).With reference to the particular features of the resulting coloring model, we propose a memetic coloring algorithm (MMCOL) to solve it.Note that 1089-778X c 2019 IEEE.Personal use is permitted, but republication/redistribution requires IEEE permission.
as a heuristic, if MMCOL finds a legal coloring for a given LSC instance, then the problem is solved.Otherwise, it says nothing about whether the LSC instance is solvable or not.We summarize the contributions of this paper as follows.
First, from a perspective of solution method, the proposed approach considers the LSC problem as a particular graph coloring problem.In this approach, we start by converting an LSC instance to a domain-constrained Latin square graph (Section II-A).Then, we reduce the graph model by applying a constraint propagation-based kernelization technique (Section II-B), leading to an instance of the list coloring problem.Finally, we seek a legal list coloring of the graph by running a dedicated memetic algorithm (Section III).The kernelization technique recursively uses constraint propagation to remove the vertices with a fixed color (corresponding to filled grids).The memetic algorithm adapts ideas from graph coloring algorithms to effectively solve the underlying list coloring problem.In particular, the algorithm integrates a problem-specific crossover to generate promising offspring solutions, an effective iterated tabu search (ITS) procedure to improve each offspring solution, and a distance-and-qualitybased pool updating strategy to ensure a healthy diversity of the population.
Second, from a perspective of computational performance, we provide experimental results on a large number of LSC benchmark instances available in the literature (over 1800 in total, including 19 traditional benchmark instances from [17] and 1800 new instances from [21]) and show comparisons with various state-of-the-art approaches, including four recent iterated local search algorithms, general ILP and exact CP solvers, and a general heuristic solver.While the reference approaches can only solve a subset of the tested instances, our approach is able to solve all the instances consistently.Such a performance has never been reported in the literature, demonstrating the high effectiveness of considering LSC as a graph coloring problem and using the proposed populationbased memetic algorithm to color Latin square graphs.We also adapt the method to the general PLSE problem and report computational results on additional 1800 PLSE benchmark instances from [21].
Third, and more generally, the proposed method can be used to solve the list coloring and precoloring extension problems, which are relevant graph models both in theory and in practice [28].Indeed, for these two important coloring problems, although the literature offers many theoretical studies on specific graphs, we are not aware of any dedicated and effective algorithm able to handle large graphs.This paper thus fills in this gap.Moreover, since precoloring extension and list coloring are useful models to formulate various applications, our method can be applied in these practical settings as well.
The rest of this paper is organized as follows.Section II describes the converted graph coloring model.Section III presents the proposed MMCOL algorithm.Section IV reports computational results obtained with the proposed method and provides comparisons with state-of-the-art algorithms.Section V shows an analysis of two key components of the method, followed by concluding comments in the last section.
The Appendix reports computational results of the proposed method on the related PLSE problem.

A. Partial Latin Square and Latin Square Graph
Let P be a partial Latin square with n × n grids, an associated graph G = (V, E), called Latin square graph, can be conveniently defined with the vertex set V = {{1, . . ., n} × {1, . . ., n}} (|V| = n 2 ) representing the grids and edge set E (|E| = n 2 (n − 1)) where {u, v} ∈ E if and only if u and v represent two grids of the same row or column [6].Then the LSC problem is equivalent to find a legal n-coloring of the associated Latin square graph G by using the colors {1, . . ., n} as follows.Let D(v) denote the color domain of vertex v of the graph.Obviously, if v corresponds to a grid already filled with number k (k ∈ {1, . . ., n}), D(v) is a singleton domain {k}; otherwise, D(v) is initially set to {1, . . ., n}.The above coloring problem is the so-called precoloring extension problem [5], where some vertices have a fixed color and the remaining vertices are to be assigned a color in {1, . . ., n}.
Note that a legal n-coloring of G can also be defined as Basically, in order to legally complete a partial Latin square, each color class must contain exactly n vertices when all the grids are filled.
Fig. 1 shows a partial Latin square of order 3, with two filled grids and seven empty grids [Fig.1(a)] and the corresponding domain-constrained graph G with nine vertices and 18 edges [Fig.1(b)].Let L xy represent the grid with xth row and yth column, then the connection between L xy and its corresponding vertex v i is given by i = (x − 1) × n + y.The objective is to find a legal three-coloring of the associated G by using the colors {1, 2, 3}.The vertices with the blue and red colors (indicated by colors 1 and 2, respectively) represent the filled grids in Fig. 1(a) while the black vertices represent the empty grids.In this example, D(v 3 ) = {1} and D(v 6 ) = {2} while the color domain of each uncolored vertex is D(v i ) = {1, 2, 3} (i = 1, 2, 4, 5, 7, 8, 9).The residual capacities of V 1 -V 3 are 2, 2, and 3, respectively (|R 1 | = 2, |R 2 | = 2, and |R 3 | = 3).Now, completing the partial Latin square is equivalent to finding a legal coloring of the graph by assigning a color in {1, 2, 3} to each uncolored vertex of G.
One notices that unlike the general graph coloring problem, the precoloring extension problem associated to a Latin square graph has a specific feature.That is, if a vertex v of the graph represents a grid already filled with k ∈ {1, . . ., n}, v has a singleton color domain D(v) = {k} and thus receives definitively the unique color k.Moreover, this color is forbidden for any vertex u adjacent to v and should be excluded from the color domain D(u).From a perspective of graph coloring, we can beneficially use this property to perform a preprocessing of the graph to obtain a reduced graph and then color the reduced graph instead of the initial Latin square graph.end for 8: end while 9: return G = (V, E)

B. Preprocessing to Simplify the Latin Square Graph
The preprocessing procedure (Algorithm 1) aims to reduce the given Latin square graph by using the colored vertices (i.e., those with a singleton color domain).For this purpose, we apply a kernelization technique based on constraint propagation [36] as follows.We first remove the precolored vertices (corresponding to the filled grids) as well as the edges connected to a colored vertex.Moreover, considering the coloring constraint stating that two adjacent vertices cannot receive the same color, once a vertex v receives color k, k is forbidden for any adjacent vertex u and can be safely removed from its color domain D(u).If the color domain of a vertex u becomes a singleton, vertex u definitively receives the unique color.Since vertex u is now a colored vertex, it can be used to further reduce the graph.This process is repeated until no color domain can be reduced.Notice that if the color domain of a vertex is reduced to the empty set during the preprocessing procedure, then the given LSC instance has no solution, i.e., it cannot be fully completed.
Consider again the example of Fig. 1.After applying the preprocessing to the Latin square graph in Fig. 2(a), we obtain the reduced graph shown in Fig. 2(b).In this particular case, since v 9 is connected to the two colored vertices v 3 and v 6 , the colors 1 and 2 are removed from the color domain of v 9 , causing D(v 9 ) to become a singleton {3}.As a result, v 9 receives the unique color 3.The color domains of other vertices adjacent to v 3 , v 6 , or v 9 are also reduced, leading to the graph of Fig. 2 In terms of graph coloring, a reduced Latin square graph like Fig. 2(b) is a domain-constrained graph because the permissible colors of a vertex are limited to a list of colors in {1, . . ., n} [instead of the whole set {1, . . ., n}].In fact, the underlying coloring problem is the so-called list coloring problem [12], [28], which, like the classic vertex coloring problem, is NP-complete in the general case.Our literature review on list coloring indicates that no practical algorithm is currently able to color large graphs.Meanwhile, it is known that the list coloring problem can be transformed to the vertex coloring problem [28].However, this transformation needs to create an auxiliary graph which is larger than the input graph by adding k ≥ n vertices and k 2 edges.Note that in the case of LSC, the Latin square graphs include already 2500-4900 vertices for n = 50, 60, 70 for the main benchmark instances tested in this paper.To our knowledge, few vertex coloring algorithms are able to effectively color graphs of these sizes given that the benchmark graphs from the well-known DIMACS Challenge (https://mat.gsia.cmu.edu/COLOR) are limited to 1000 vertices.For these reasons, we introduce below a dedicated algorithm specifically designed to solve the list coloring problem of Latin square graphs.

III. MEMETIC ALGORITHM FOR COLORING LATIN SQUARE GRAPHS
We describe in this section the population-based memetic algorithm for coloring domain-constrained Latin square graphs, i.e., solving the associated list coloring problem where each vertex v can only take a color from its given color domain D(v).
The algorithm takes a reduced Latin square graph G as its input and tries to find a legal list coloring in the search space defined in Section III-B.For this purpose, the algorithm starts with an initial population [line 1 (Section III-C) in Algorithm 2].Then, to find a legal n-coloring, MMCOL repeats a number of generations to improve the population until a stopping condition (limited to maxGenerations) is met.At each generation, MMCOL randomly selects two parent colorings from the population [line 5 (Section III-C)

B. Search Space and Evaluation Function
Let G = (V, E) be a Latin square graph with L vertices {v 1 , . . ., v L } and color domains D(v i ) ⊆ {1, . . ., n} (i = 1, . . ., L).Our MMCOL algorithm explores the following space C of candidate list colorings: Given a candidate coloring c in C, if c(u) = c(v) and {u, v} ∈ E (i.e., two adjacent vertices u and v receive the same color), then {u, v} is a conflicting edge in c while u and v are called conflicting vertices.To assess the quality of the coloring c, we use the evaluation or fitness function f given in Equation ( 1), which counts the number of conflicting edges in c Consequently, if f (c) = 0, c is conflict-free and identifies a legal list coloring.Otherwise (f (c) > 0), c is an illegal coloring with conflicting edges.For two candidate solutions c 1 and c 2 , c 1 is considered to be better than c 2 if f (c 1 ) < f (c 2 ) (c 1 contains fewer conflicting edges).
Given the above evaluation function, the objective of our memetic algorithm is to find a legal (conflict-free) list coloring in the search space C by minimizing f .

C. Population Initialization
The MMCOL algorithm applies a randomized coloring strategy to create the initial population P that is composed of p colorings sampled in C (p is the population size and set to 20 in this paper).Let G = (V, E) be the given graph with V = {v 1 , . . ., v L } and D(v i ) ⊆ {1, . . ., n} (i = 1, . . ., L).
To build an initial coloring of G, we iteratively select an uncolored vertex v at random and then assign it a random color k from its color domain D(v).Such an initial solution can be obtained very quickly in O(L), but may involve a high number of conflicting edges.To obtain an initial coloring of reasonable quality, we improve this solution by the local optimization procedure (see Section III-E) and then insert the improved solution into the population if the solution does not exist in P. Otherwise, the solution is discarded and a new solution is generated.This initialization process is repeated until the population is filled up with p different colorings.

D. Crossover Operator
Recombination is an important component of our MMCOL algorithm that aims to transmit meaningful features from parents to offspring solutions [20].For the conventional graph coloring problem, the greedy partition crossover (GPX) [14] is known to be highly effective.However, given that list coloring graphs have restricted color domains (instead of the set {1, . . ., n}), GPX cannot be applied directly in the context of the list coloring problem.On the other hand, the key idea of GPX, i.e., inheriting large color classes, is of interest even in the case of list coloring.As a result, we propose an adaptation of GPX by taking into account the constrained color domains of our graphs.This leads to our maximum approximate group-based crossover (MAGX) for Latin square graph coloring.
The proposed MAGX crossover operator generates one offspring solution from two randomly selected parent solutions (see Algorithm 3).Let P 1 = {V Identify from parents (P 1 and P 2 ) the largest color class Remove the vertices of V b i from P 1 and P 2 6: g ← g + 1 7: end while 8: for each empty color class V o i in o do 9: v is randomly assigned a color from its color domain D(v) 13: end for 14: return o Third, for each vertex v missing in o, v is assigned a random color class in its color domain D(v).
At this stage, a complete offspring solution o is obtained.In case that the offspring is the same as one of the parent solutions (this rarely happens), MAGX applies a slightly different strategy for the first phase such that the largest color class is selected by considering alternatively P 1 and P 2 (instead of considering simultaneously P 1 and P 2 ).Since the three phases have a time complexity of O(n 2 ), O(n), and O(n 2 ), respectively, the time complexity of the MAGX crossover is bounded by O(n 2 ).
Fig. 3 shows an illustration example of the MAGX crossover.This example involves a Latin square graph of order 3 (n = 3) with nine vertices a, b, c, d, e, f , g, h, and i

E. Iterated Tabu Search Procedure
The ITS procedure (Algorithm 4) takes an offspring solution c generated by the MAGX crossover operator as its input and tries to improve its quality in terms of fitness function f (Equation (1), Section III-B).For this purpose, ITS iterates between a tabu search (TS) procedure followed by a relaxation-based perturbation procedure to try to attain a legal coloring by resolving the conflicts (lines 3-13 in Algorithm 4).TS iteratively improves c by recoloring conflicting vertices (see Section III-B).At the end of each TS run, if the conflicts are resolved, then a legal list coloring c * is found, and the whole search terminates immediately.If conflicts remain in the solution, ITS triggers a perturbation procedure to modify the solution and then uses the modified solution as its starting solution for the next TS run (line 9 in Algorithm 4).ITS repeats the above process until a prefixed maximum number of iterations maxLSIters is reached or a legal coloring is obtained.
1) Tabu Search-Based Coloring Procedure: As its key optimization procedure, ITS uses the TS method [15] to improve a given illegal list coloring.Specifically, the TS procedure used in this paper is based on the implementations presented in [10] and [14] of the popular TabuCol algorithm for the conventional graph coloring problem [22].Suppose that the solution c is composed of L vertices {v 1 , v 2 , . . ., v L } and each vertex v i receives a permissible color in its constrained color domain D(v i ) (i ∈ {1, 2, . . ., L}).The TS procedure end if 13: until maxLSIters is reached explores the space C composed of all possible colorings (see Section III-B) to seek a legal list coloring.
To improve the solution, TS iteratively makes transitions from the current solution c to one neighboring solution.To obtain a neighboring solution c from solution c, TS displaces a conflicting vertex v from its current color class V i to another eligible color class V j such that j ∈ D(v) [i.e., the current color i of vertex v is changed to a new permissible color j in v's color domain D(v)].Thus, c and c differ only by the color of a conflicting vertex v. Since the color domains are bounded by n, the size of this neighborhood is bounded by O(n c × n), where n c is the number of conflicting vertices in coloring c.At each iteration, TS selects among the eligible neighboring solutions the best neighbor c b according to the evaluation function f (Equation (1), Section III-B) and uses c b to replace c.A neighboring solution is eligible if it is not forbidden by the tabu list (see the explanation below) or if it is better than the best recorded solution.Suppose that the selected neighboring solution is obtained by changing the color i of conflicting vertex v, (v, i) is recorded in the tabu list, indicating that vertex v is forbidden to receive the color i again for the next β consecutive iterations (β is called the tabu tenure).Following [10] and [14], β is dynamically tuned by β = 0.6 * f (c)+random (10), where random (10) is a random number in {1, . . ., 10}.The TS procedure stops when its iteration counter reaches the given limit α (α is called the TS depth).The best solution c and the number of conflicting edges in c recorded during the search are returned as its output when the procedure terminates.
2) Relaxation-Based Perturbation: It is possible that no legal list coloring is found at the end of a TS run (see line 8 in Algorithm 4).In this case, the search is considered to be trapped in a local optimum and we trigger a relaxation-based perturbation procedure to escape from the trap.The overall procedure of the relaxation-based perturbation is illustrated in Fig. 4. Let X ⊂ V be the set of conflicting vertices (i.e., each vertex of X is involved in at least one conflicting edge in c).The perturbation procedure is performed in three steps.
1) Extract a subgraph G by randomly removing |X|/2 conflicting vertices along with the incident edges.2) Improve the coloring on G using TS.

3) Construct a new coloring on G by getting it back to G.
The perturbation procedure is based on the consideration that the conflicting vertices of the local optimum are critical vertices for obtaining a legal list coloring.Meanwhile, these are also difficult vertices for conflict resolution.By ignoring some of these difficult vertices, TS has a higher chance to resolve the conflicts of the relaxed subproblem, thus providing new search opportunity when the improved solution of the relaxed subproblem is added back to the ignored partial solution.Notice that, in case that the improved c * p is not a legal coloring after the second step of the perturbation procedure, c * p is still a high quality partial coloring which could be close to a complete solution.Using c * p as its starting point to be extended, TS will explore a new search trajectory and hopefully encounters a legal coloring.

F. Population Updating
In order to avoid premature convergence of our MMCOL algorithm, we apply a population updating strategy similar to those used in [24], [27], [33], and [38].The adopted strategy simultaneously considers solution quality and diversity when using an offspring solution to update the population.
Given two list colorings c i and c j , we use the so-called set-theoretic partition or transfer distance D i,j [34], [35] to measure the dissimilarity of c i and c j , which is defined as the minimum number of vertices that need to be moved between color classes of c i to transform c i to c j .The diversity between one solution and the entire population P is given by D i,P = min j =i {D ij }.Furthermore, we define the goodness score of one n-colorings c i of P in terms of both solution quality and diversity by s(c i ) = f (c i ) + e 0.08n 2 /D i,P , ∀c i ∈ P where f (c i ) is the number of conflicting edges of c i [27].A small (large) s(c i ) value indicates a good (bad) solution with respect to the individuals of P. Given offspring o, the population P is updated with o according to the following procedure.
Step Step 4: If c w is o, remove c w with probability 0.8 and remove c sw with probability 0.2.This updating strategy ensures that the individuals of the population are not only of high quality, but also sufficiently distanced.This property provides a basis for the random strategy used in our algorithm to select the parents for the crossover.

IV. EXPERIMENTAL RESULTS
In this section, we assess the proposed approach by reporting computational results on the LSC problem and showing comparisons with state-of-the-art methods.We show in the Appendix additional results on the related PLSE problem.

A. Benchmark Instances and Experimental Protocol
To evaluate the performance of the proposed approach for solving LSC, we carry out extensive experiments on the set of 1800 random LSC benchmark instances recently introduced in [21]. 1 These LSC benchmark instances (named as LSCn-r) are evenly divided into 18 types (n ∈ {50, 60, 70} and r ∈ {0.3, 0.4, 0.5, 0.6, 0.7, 0.8}), where n is the order of the partial Latin square and r (r ∈ [0, 1]) denotes the ratio of filled grids over the n×n grids.So each type (n, r) has 100 instances.These instances were generated by randomly removing (1 − r)n2 grids from an arbitrary legal Latin square.As a result, these LSC instances always admit a complete Latin square.
By converting these instances to Latin square graphs (see Section II-A), we obtain graphs with 2500-4900 vertices and 122 500-338 100 edges. 2 To solve each instance, we first apply the preprocessing procedure of Section II-B to obtain a reduced list coloring graph which is then colored with the MMCOL algorithm.The preprocessing step takes typically from several seconds to dozens of seconds.
In addition to these 1800 random instances, we also assess our approach on the set of 19 traditional benchmark instances from the COLOR03 competition3 that were tested, for instance, [17] and [21].
MMCOL was coded in C++ 4 and compiled using g++ with the "−O3" option on a computer running Linux equipped with 2.83 GHz and 8-GB RAM.When running the DIMACS machine benchmark procedure "dfmax.c"5on our machine, we obtain the following results: 0.20, 1.23, and 4.68 s for graphs r300.5, r400.5, and r500.5, respectively.The computational results reported in this section were obtained with the parameter setting shown in Table I.
In the following sections, we first show the results on the 19 traditional instances, and then present a comparative analysis of our computational results on the large set of 1800 benchmark instances with respect to the state-of-the-art results in the literature.Given the stochastic nature of MMCOL, each instance was independently solved 30 times with different random seeds.

B. Results on 19 Traditional Benchmark Instances
The computational results of MMCOL on the 19 traditional Latin square graphs are summarized in Table II.Columns 1-3 of Table II indicate the characteristics of each instance: the name, the Latin square order n, and the ratio r.Columns 4 and 5 present the success rate over 30 trials (SR) and the computation time over the successful runs t(s) in seconds (a successful run means that a legal Latin square is attained for this run).From Table II, one observes that MMCOL can complete the partial Latin square for all the 19 instances.Besides, our MMCOL requires a very short computation time even for the large instances with n ≥ 50.Moreover, the last five instances are critically constrained and fully "balanced," where the number of empty grids is approximately the same over rows and columns.These instances are known to be particularly difficult [25] and only the two smallest ones of these five instances (qwhdec.order33.holes381.bal.1 and qwhdec.order50.holes825.bal.1)can be solved by very few approaches presented in [17] including CSP, hybrid strategy mixing LP/CSP, and SAT-based method.The difficulty of these instances are further confirmed by the most recent study reported in [21], where even the best performing heuristic Tr-ILS* [21] cannot solve any of these balanced instances.We also ran the source code of Tr-ILS* on our computer for a long computation time of 3600 s and still failed to solve any of these five balanced instances.It is thus remarkable that our MMCOL approach solves these instances consistently, even though MMCOL has a low success rate for two instances.We conclude that MMCOL performs very competitively with respect to all of the existing approaches for solving these traditional instances.

C. Comparative Results on the Set of 1800 Benchmark Instances
Table III summarizes the computational statistics of our MMCOL algorithm on the set of 1800 benchmark instances, together with the results of seven most recent methods in the literature reported in [21].The reference methods include CPX-IP, CPX-CP, LSSOL, 1-ILS*, 2-ILS, 3-ILS, and Tr-ILS*, where CPX-IP and CPX-CP are integer programming (IP) and constraint programming (CP) solvers from IBM/ILOG CPLEX, LSSOL denotes the LocalSolver, 6 and 1-ILS*, 2-ILS, 3-ILS, and Tr-ILS* are four iterated local search algorithms with (1, ∞)-neighborhood, (2, ∞)-neighborhood, (3, ∞)neighborhood, and Trellis-neighborhood search, respectively.All the reference algorithms are performed on an Intel core i7-4770 processor with 3.90 GHz and 8-GB RAM (which is faster than our computer), with a time limit of 30 s for CPX-IP, CPX-CP, and LSSOL, and 10 s for 1-ILS*, 2-ILS, 3-ILS, and Tr-ILS*.Table IV additionally presents the detailed results of our approach on a subset of 600 difficult benchmark instances.
Columns 1-3 of Table III show the characteristics of the tested instances: the order n of each Latin square, the ratio r, and the number of instances Inst# for each type (n, r).Following [21], columns 4-10 present the results of the seven reference algorithms (CPX-IP, CPX-CP, LSSOL, 1-ILS*, 2-ILS, Tr-ILS*, and 3-ILS), "Suc#" shows for each type of 100 instances the number of instances for which an algorithm can obtain a legal solution.Columns 11 and 12 give the results of our MMCOL algorithm in terms of "Suc#" and the average time t avg (s) in seconds is defined by where t i (s) is the average time over the successful runs for the ith instance.The last row (perfect success times) shows the number of instance types with Suc# = 100, i.e., the number of instance types among the 18 types (n, r) for which all the 100 instances are solved by an algorithm.
From Table III, one observes that both exact methods CPX-IP and CPX-CP solve all 100 instances for only three out of the 18 types.The five ILS heuristics (LSSOL, 1-ILS*, 2-ILS, 3-ILS, and Tr-ILS*) solve all 100 instances for 1, 4, 5, 1, and 6 of 18 types (in bold), respectively.In contrast, our MMCOL algorithm can solve all the instances for all 18 types.The average time to find a solution is less than 11 s except for the 300 instances with r = 0.7 for which MMCOL needs less than 300 s to attain a solution while all reference algorithms fail to solve any of these instances.Besides, we observe that the seven reference algorithms have a worse performance for the types (n ∈ {50, 60, 70}, r ∈ {0.6, 0.7}) which are known to be more difficult [21].On the other hand, MMCOL has no particular difficulty to solve these hard instances.In order to verify if the ILS algorithms can solve more instances by using more computation time, we ran the source code of the best performing algorithm Tr-ILS* under a much relaxed time limit of 3600 s on the instances of types (n ∈ {50, 60, 70}, r ∈ {0.6, 0.7}).One observes that the 100 instances of LSC-50-60 (n = 50, r = 0.6) can be fully completed, 98 and 95 instances of LSC-60-60 (n = 60, r = 0.6) and LSC-70-60 (n = 70, r = 0.6) can be fully completed, respectively.Nevertheless, no instance of the types (n ∈ {50, 60, 70}, r = 0.7) can be fully completed even if these instances have more filled grids for the LSC.In summary, MMCOL competes very favorably with seven most recent methods in the literature and proves to be highly effective in solving the set of 1800 benchmark instances with no exception.

V. ANALYSIS OF TWO KEY COMPONENTS
In this section, we present an analysis of two key components of the proposed method: 1) the constraint propagationbased kernelization and 2) the role of the MAGX crossover operator.

A. Impact of Constraint Propagation-Based Kernelization
The constraint propagation-based kernelization of Section II-B is used to preprocess an initial Latin square graph and thus reduces the search space of the subsequent list coloring task.In order to investigate the impact of this kernelization technique, we evaluate the reduced vertices and the color domains for the 1800 benchmark instances.
Table V summarizes the statistics of the reduced graphs for each type of instances.Columns 1-3 recall the instance characteristics: the order n, the ratio r, and the number of instances Inst#.For each type (n, r) of 100 instances, column 4 "#V avg " indicates the average number of the precolored vertices (i.e., the filled grids) in the initial Latin square graph.Since the color domain of some vertices becomes a singleton during the kernelization process, the graph can be further reduced.Hence, column 5 "#V avg " shows the average number of further reduced vertices.Column 6 "#D avg " presents for each type of 100 instances the average cardinality of the color domains after kernelization, i.e., #D avg = [( j=Inst# j=1 . And column 7 "Solved#" gives the number of instances that are solved during the kernelization process (i.e., all color domains are reduced to a singleton).
From Table V, one observes that the vertices of the converted graphs are dramatically reduced by the kernelization technique.For the types of instances [n ∈ {50, 60, 70}, r = 0.3, 0.4, 0.5 and (n = 60, r = 0.6)], only the precolored vertices (the filled grids) are removed from the graphs, meaning that removing the precolored vertices does not lead to any new singleton domain.However, for the types of instances (n ∈ {50, 60, 70}, r = 0.8), many vertices can be further removed by the kernelization technique.Furthermore, the color domains of the uncolored vertices are obviously reduced from their initial sizes n.In particular, for 32 instances of the type (n = 50, r = 0.8) and nine instances of the type (n = 60, r = 0.8), a solution is found during kernelization (i.e., all color domains are reduced to a singleton), without needing to run the subsequent coloring algorithm.We conclude that the kernelization technique plays an important role in reducing the initial Latin square graphs and helps to ease the subsequent list coloring task.

B. Impact of the MAGX Crossover Operator
Within the proposed memetic algorithm, the MAGX crossover described in Section III-D is another key component.To assess its impact, we present an experiment to compare MMCOL (with the crossover MAGX), MMCOL' (with an uniform assignment crossover), and ITS (without MAGX).The uniform assignment crossover builds an offspring solution by inheriting, for each vertex, the color either from parent P 1 or P 2 with an equal probability of 0.5.For a fair comparison, we use the same parameter setting for MMCOL and MMCOL' and set maxLSIters = 10 4 for each ITS run, which corresponds roughly to the same search effort of MMCOL/MMCOL' with 100 generations (see Table I).The initial solutions of all the algorithms are generated by the initialization procedure given in Section III-C.We ran MMCOL, MMCOL', and ITS 30 times to solve each of the 600 difficult benchmark instances (n ∈ {50, 60, 70}, r ∈ {0.6, 0.7}) listed in Table IV.The comparative results are shown in Fig. 5 where we indicate for each type of 100 instances, the number of instances for which an algorithm can find a solution with a perfect success rate 30/30 (i.e., each of its 30 trials finds a legal list coloring).
Additionally, Table VI presents the detailed results on the 22 most difficult instances of LSC-50-70 where we show for each instance and each algorithm (MMCOL, MMCOL', and ITS), the success rate SR over 30 trials and the average computation time t(s) in seconds over the successful runs.This table indicates that MMCOL obtains a higher success rate SR than MMCOL' and ITS for 17 and 22 instances, respectively (in bold), with shorter computation times t(s) except for two cases.This experiment confirms that the population-based evolutionary framework implemented in our memetic algorithm and its crossover operator contribute positively to the performance of the proposed algorithm.

VI. CONCLUSION
We have proposed an approach to solve the LSC problem by converting LSC to a domain-constrained graph coloring problem.By taking advantage of the particular features of Latin square graphs, we have developed a constraint propagation-based kernelization technique to preprocess the given Latin square graph to obtain a reduced graph, for which an associated list coloring problem is defined.To effectively solve the list coloring problem, we have devised a dedicated memetic algorithm MMCOL which integrates a tailored crossover operator to generate new solutions, an ITS procedure to improve each offspring solution and a distance-qualitybased pool updating strategy to ensure a healthy diversity of the population.
Extensive evaluations on a large number of benchmark instances in the literature (19 traditional instances and 1800 random instances) have shown that the proposed approach performs very well with respect to the state-of-the-art methods including those introduced very recently in 2016.In particular, our approach is able to find a solution for all the benchmark instances consistently and effectively, a performance never attained by any existing approach.We have also used a slightly modified version of the method to solve the general PLSE problem and reported computational results on the set of 2018 PLSE instances in the Appendix.
Given that LSC and PLSE have a number of applications, the proposed approach can help to solve these applications.More generally, the method proposed in this paper can be used to approximate the important list coloring and precoloring extension problems, for which few practical algorithms exist in the literature.The method will be particularly useful if large problem instances are considered.
For future work, it would be interesting to identify additional features of Latin square graphs and use them to design effective search strategies and operators.Approaches based on vertex coloring algorithms are also worthy of investigation.

APPENDIX RESULTS ON THE PARTIAL LATIN SQUARE EXTENSION PROBLEM
We now show that the method presented in this paper can be used to solve the related and more general PLSE problem.Recall that PLSE is to assign numbers {1, . . ., n} to as many empty grids as possible under the condition that each number occurs at most once in each row and each column.To apply the proposed method to solve PLSE, we make the following two adjustments.
First, unlike LSC, the color domains of some vertices in the case of PLSE may become empty during the constraint propagation-based kernalization process, implying that the corresponding grids cannot be legally filled by any given number.In terms of graph coloring, if the color domain of a vertex v becomes empty when applying the preprocessing procedure, any color for vertex v is definitively conflicting with at least one of its adjacent vertices.If this happens, we randomly assign a color {1, . . ., n} to vertex v and keep this color unchanged during the search process.
Second, at the end of the MMCOL algorithm, there are two possibilities.If the returned final coloring c * is conflictfree [i.e., f (c * ) = 0], the given partial Latin square is fully completed and an optimal solution is found.Otherwise, some vertices are assigned conflicting colors in c * [i.e., f (c * ) > 0].In this case, we obtain a legal partial solution by dropping from c * some conflicting vertices.The dropped vertices correspond to the grids that cannot be legally filled while the remaining vertices in the legal partial coloring define a solution for the given PLSE instance.To remove conflicting vertices, we first drop any vertex v with empty color domain caused by the preprocessing procedure.Then, if conflicts remain, we repetitively remove the vertex u which is conflicting with the largest number of other vertices in the coloring until we obtain a partial conflict-free coloring.
Table VII summarizes the results of our MMCOL algorithm and seven reference methods on the set of 1800 PLSE benchmark instances introduced in [21].Like the LSC instances of Section IV-A, these 1800 PLSE instances are evenly divided into 18 types (n ∈ {50, 60, 70}, r ∈ {0.3, 0.4, 0.5, 0.6, 0.7, 0.8}) [so 100 instances per type (n, r)].For these instances, n 2 is a trivial upper bound of their optimal solutions.In this experiment, we used the same experimental condition as for solving LSC (Section IV-C).Like [21], we solved each instance once.Columns 1-3 indicate the characteristics of the instances with the same information as in Table III.Columns 4-17 present, for each reference algorithm and for each type (n, r) of 100 instances, the number of fully completed Latin squares "suc#," and the average completed grids over 100 instances "Avg."Columns 18-20 show the results of our MMCOL algorithm in terms of suc#, Avg, and the average computation time "t avg (s)" in seconds.Notice that if a partial Latin square is fully completed, the optimum is attained (so 1186 instances out of 1800 are solved to optimality).Otherwise, the reported result in terms of the filled grids gives a lower bound of the given PLSE instance.
Table VII shows that MMCOL obtains improved or equal average results for 15 out of 18 types (in bold) except the types (n ∈ {50, 60, 70}, r = 0.8).In particular, for the 1000 instances of types (n ∈ {50, 60, 70}, r ∈ {0.3, 0.4, 0.5}), (n = 70, r = 0.6), and 186 instances of types (n ∈ {50, 60}, r = 0.6), (n = 70, r = 0.7), MMCOL attains an optimal solution.Meanwhile, for the instances of types (n ∈ {50, 60, 70}, r = 0.8), MMCOL performs worse than the reference algorithms.We mention that since the instances of types (n ∈ {50, 60, 70}, r = 0.8) are strongly constrained, the color domains of some uncolored vertices are reduced to the empty set during the constraint propagation-based preprocessing of Section II-B, indicating that these instances cannot be fully completed.Let α > 0 be the number of vertices with an empty color domain identified during the preprocessing, n 2 −α defines an upper bound of the given instance, which is strictly tighter than the trivial n 2 bound.

Fig. 2 .Algorithm 2 1 ,
Fig. 2. Latin square graph of Fig. 1(b) and the residual Latin square graph obtained by the preprocessing procedure.
in Algorithm 2] and recombines them to generate an offspring coloring by a dedicated crossover operator [line 6 (Section III-D) in Algorithm 2].This offspring coloring is then improved by an ITS procedure [line 7 (Section III-E) in Algorithm 2].Finally, the improved offspring is used to update the population according to an updating strategy based on a distance-quality criterion [line 11 (Section III-F) in Algorithm 2].During the memetic search process, if a legal coloring is found, MMCOL stops and returns the legal coloring found.
to be assigned to three color classes V 1 -V 3 .Suppose that the color domainsD(a) = D(c) = D(g) = {1, 3}, D(h) = {1, 2}, and D(x) = {1, 2, 3} for x ∈ {b, d, e, f , i}.At the beginning, no color class exists in o, so |R i | = 3 (i = 1, 2, 3).In the first step, V 2 = {d, e, f } of P 1 is identified as the largest color class whose |V 2 | ≤ |R 2 |and V 2 of the offspring o is empty.Thus, this color class {d, e, f } becomes the color class V 2 of the offspring and the vertices d, e, and f are removed from both P 1 and P 2 .Notice that due to the fact that vertices may have different color domains, the vertices of the inherited color class {d, e, f } of o receives the same color as the donor parent (here color 2).Similarly, V 3 = {b, c, i} and V 1 = {a, g} of P 2 are inherited as color classes V 3 and V 1 of o.After these operations, vertex h is still missing in o.Since this vertex belongs to different classes in P 1 and P 2 , h is assigned a random color class from its color domain D(h) = {1, 2}.

Algorithm 4 repeat 4 : 5 :else 11 :
Pseudo-Code of ITS Require: A n-coloring c, depth of tabu search α, color domain D(v) of each vertex v ∈ V Ensure: A legal coloring c * 1: c * ← c; // c * records the best solution found so far 2: f * ← f (c * ); // f * records the smallest number of conflicting edges 3:(c, f ) ← TS(c, α); // Apply the tabu search procedure with search depth α to improve the input coloring c, see Sect.III-E1if f < f * then 6: c * ← c; f * ← f (c); f ) ← Perturbation_Procedure(c); //Apply the perturbation procedure to locate at a promising region, see Sect.III-E2 10: return the legal coloring c * ; 12:

1 :
Insert the offspring solution o into P and compute the score s(c i ) of each individual c i of P. Step 2: Identify the worst individual c w (i.e., with the largest value of the scoring function s) and second worst individual c sw (with the second largest s value).Step 3: If c w is different from o, remove c w from P.
First, MAGX builds a number of color classes of o by inheriting color classes from the parent solutions.To build a new color class, MAGX selects, among the color classes of both P 1 and P 2 , one largest class (call it V b i ) such that its cardinality does not exceed the residual capacity R i of the corresponding color class (line 3 in Algorithm 3).MAGX then uses V b Parent solutionsP 1 = {V 1 1 , . .., V 1 n }, P 2 = {V 2 1 , . .., V 2 n },and color domain D(v) of each vertex v ∈ V Ensure: An offspring solution o = {V o 1 , . . ., V o n } 1: g ← 0 // Count the number of color classes already built in o 2: while g < n do [14] ..., V 1 n } and P 2 = {V21 , ..., V 2 n } be the parent solutions, MAGX generates, in three phases, the offspring solution o = {V o 1 , ..., V o n }, where each V o i (i = 1, ..., n) is initially set to empty. the recombination operation (like GPX of[14]does), for our list coloring problem of Latin square graphs, each color class of the offspring must inherit the color of the selected parent due to the constrained color domains of the vertices.Second, for each color j such thatV o j = ∅ in o, if V 1 j and V 2j share common vertices, these vertices are used to form the color class V o j of the offspring.Algorithm 3 Pseudo-Code of the MAGX Crossover OperatorRequire: For the residual vertices, transmit the vertices that share the same color in both parents 10: end for 11: for each uncolored v ∈ V in o do

TABLE III COMPARATIVE
RESULTS OF MMCOL WITH BEST-PERFORMING ALGORITHMS ON THE SET OF 1800 LSC BENCHMARK INSTANCES

TABLE IV DETAILED
COMPUTATIONAL RESULTS OF MMCOL ON A SUBSET OF 600 INSTANCES Furthermore, we show in TableIVthe detailed results of MMCOL on the 600 difficult instances (n ∈ {50, 60, 70}, r ∈ {0.6, 0.7}) for which most reference algorithms perform badly.For each instance, we present the success rate over 30 trials SR and the average computation time over the successful runs t(s) in seconds.From TableIV, we observe that MMCOL achieves the perfect success rate 30/30 on three types (n ∈ {50, 60, 70}, r = 0.6).MMCOL has a lower success rate only for 28

TABLE V STATISTICS
OF THE REDUCED LATIN SQUARE GRAPHS AFTER KERNELIZATION

TABLE VI COMPARATIVE
RESULTS ON 22 DIFFICULT CASES OF LSC-50-70