Nicolas Dupin
email: nicolas.dupin@univ-angers.fr

Matheuristic variants of DSATUR for the vertex coloring problem

Keywords: Matheuristic, Vertex Coloring, DSATUR, cliques, dual heuristic

If DSATUR heuristic was proven efficient for many instances of the vertex coloring problem, this paper aims to analyze inefficiency causes of DSATUR heuristic, extending DSATUR with matheuristic operators. Using an Integer Linear Programming formulation allows to have larger local greedy optimization in DSATUR construction scheme. Matheuristics allow also to initialize DSATUR heuristics and matheuristics with cliques or a partial optimal coloring. Dual bounds are also obtained, improving lower bounds implied by cliques computed by matheuristics. Computational analyses are provided, highlighting the strengths and weaknesses of DATUR heuristic and matheuristics.

Introduction

The vertex coloring problem (VCP) is one of the most widely studied and popular optimization problems in graph theory, for theoretical insights and practical applications like planning problems or interference avoidance in telecommunications, with many works in exact and heuristic methods to design efficient solvers for VCP [START_REF] Malaguti | A survey on vertex coloring problems[END_REF]. VCP is a NP-hard problem [START_REF] Garey | The complexity of near-optimal graph coloring[END_REF].

Exact methods can solve VCP to optimality with hundreds of vertices for difficult instances with Integer Linear Programming (ILP) techniques [START_REF] Furini | Exact weighted vertex coloring via branch-and-price[END_REF][START_REF] Cornaz | Solving vertex coloring problems as maximum weight stable set problems[END_REF]. First compact ILP formulation, as an assignment problem, is quite inefficient with symmetry issues, and was improved with cutting planes [START_REF] Méndez-Díaz | A cutting plane algorithm for graph coloring[END_REF]. A more recent compact formulation, based on representatives, breaks symmetries and improve the LP relaxation [START_REF] Campêlo | On the asymmetric representatives formulation for the vertex coloring problem[END_REF]. The most efficient ILP resolution techniques use an extended formulation with column generation (CG) algorithms [START_REF] Furini | Exact weighted vertex coloring via branch-and-price[END_REF] or solve VCP after transformation into maximum weight stable set problems [START_REF] Cornaz | Solving vertex coloring problems as maximum weight stable set problems[END_REF] without CG.

Among the best seminal constructive heuristics for VCP, two adaptive greedy algorithms, namely DSATUR [START_REF] Brélaz | New methods to color the vertices of a graph[END_REF] and RLF [START_REF] Leighton | A graph coloring algorithm for large scheduling problems[END_REF], are the most efficient. Some previous works analyzed inefficiency issues of such constructive heuristics [START_REF] Chiarandini | Efficiency issues in the rlf heuristic for graph coloring[END_REF][START_REF] Janczewski | The smallest hard-to-color graph for algorithm dsatur[END_REF]. Note that an exact tree search method was derived from DSATUR, and was proven to be efficient [START_REF] Furini | An improved dsatur-based branch-andbound algorithm for the vertex coloring problem[END_REF][START_REF] Segundo | A new dsatur-based algorithm for exact vertex coloring[END_REF].

Among recent trends to hybridize heuristics, matheuristics rely on exact methods to design heuristics that scales better than exact approaches, may provide results on design of meta-heuristics [START_REF] Boschetti | Matheuristics: using mathematics for heuristic design[END_REF][START_REF] Dupin | Matheuristics to optimize refueling and maintenance planning of nuclear power plants[END_REF], or furnish both lower and upper bounds, with dual bounds provided by varied relaxations [START_REF] Boschetti | Matheuristics: survey and synthesis[END_REF][START_REF] Dupin | Machine learning-guided dual heuristics and new lower bounds for the refueling and maintenance planning problem of nuclear power plants[END_REF]. This paper investigates such methodology to develop matheuristic variants of DSATUR for VCP. The goal is not only to improve DSATUR heuristic, but also to understand better inefficiency issues of DSATUR. Note that matheuristics were used recently for VCP by [START_REF] Chandrasekharan | A constructive matheuristic approach for the vertex colouring problem[END_REF] and for B-coloring [START_REF] Melo | A matheuristic approach for the b-coloring problem using integer programming and a multi-start multi-greedy randomized metaheuristic[END_REF].

2 Problem statements

Definitions and notation

Let G = (V, E) be an undirected graph on a finite vertex set V and edge set E. Cardinalities of V and E are denoted n = |V | and m = |E| < n(n -1)/2. We write V = {v 1 , . . . , v n }, and denote I = [[1; n]] = [1; n] ∩ Z. An edge e ∈ E linking vertices v i and v j is denoted e = (v i , v j ) with i < j. Note that for VCP there is no sense to consider loop edges (v i , v i ) or multiple edges. For some equations, we consider unordered notation ngb(i, j) = 1 if vertices v i and v j are linked by an edge, i.e. (v min(i,j) , v max(i,j) ) ∈ E, otherwise ngb(i, j) = 0. For i ∈ I, δ i denotes the set of indexes corresponding to neighbors of vertex v i in G:

δ i = {j ∈ I, ngb(i, j) = 1} (1)
The degree of vertex v i is denoted

d i = |δ i |.
A k-coloring of G is an assignment of k colors to vertices V such that adjacent nodes do not share the same color. VCP minimizes the number of color k to have a k-coloring of G, the optimum is the chromatic number of graph G. Formally, a proper k-coloring is denoted (c) where

c i ∈ [[1; k]] denotes the color of vertex v i , fulfilling: ∀i < j, (v i , v j ) ∈ E =⇒ c i ̸ = c j (2)
A clique in G is a subset of vertices C ⊂ V that forms a complete sub-graph of G, i.e. each couple of vertices of C is linked by an edge in G. The cardinality of any clique gives a first lower bound for the chromatic number, an optimal coloring necessarily implies different colors for the clique.

For constructive heuristics, we define a partial k-coloring (c) where

c i ∈ [[1; k]] ∪ {-1}. if c i > 0, c i is the color of vertex v i , otherwise c i = -1 denotes that vertex v i is not colored yet. A partial k-coloring (c) fulfill: ∀i < j, (v i , v j ) ∈ E =⇒ (c i ̸ = c j or c i = c j = -1) (3)
For a given partial k-coloring (c), we consider for each vertex v i the saturation table S i as the set of the assigned colors of the neighbors of v i . The saturation (degree) denotes

s i = |S i |. S i = n j∈δi {c j } \ {-1} (4)
We denote with ≽ the total order among vertices, as the hierarchic order comparing firstly the saturation degrees and then the degree:

∀v i , v j ∈ V, v i ≽ v j ⇐⇒ s i > s j or ( s i = s j and d i ⩾ d j ) (5)

Compact ILP formulations

This section presents two compact ILP formulations for VCP.

Assignment

s.t :

k c=1 z i,c = 1 , ∀i ∈ I (7)
z i,c + z j,c ⩽ y c , ∀(v i , v j ) ∈ E, ∀c ∈ [[1; k]] (8)
Objective function [START_REF] Chiarandini | Efficiency issues in the rlf heuristic for graph coloring[END_REF] counts with variables y the number of colors used. Constraints [START_REF] Cornaz | Solving vertex coloring problems as maximum weight stable set problems[END_REF] ensure that each vertex is colored. Constraints (8) expresses incompatibility of neighbor vertices to have the same color.

Having as initial value k an upper bound of the chromatic number, given by a primal heuristic, or simply k = |V | ensures the feasibility to compute the chromatic number and an optimal assignment of colors solving this last ILP. However, for efficiency issues with ILP resolution, value of k should be as small as possible. Symmetries in this encoding, for instance by permutation of colors, is a bottleneck to solve large VCP with a Branch&Bound (B&B) tree search.

Representatives ILP model A compact formulation, based on representatives, breaks symmetries and improves the LP relaxation [START_REF] Campêlo | On the asymmetric representatives formulation for the vertex coloring problem[END_REF]. Binary variables z i,i ′ ∈ {0, 1}, are defined for all i, i ′ ∈ V with i ⩽ i ′ , z i, ′ = 1 if and only if vertices v i et v i ′ have the same color, and i is the minimum index of its color. It induces following asymmetric ILP formulation:

min z n i=1 x i,i (9)
s.c :

i ′ ⩽i x i ′ ,i ⩾ 1 , ∀i ∈ I (10) x j,i + x j,i ′ ⩽ x j,j ∀(v i , v i ′ ) ∈ E, ∀j ⩽ i, (11)
Objective function [START_REF] Dupin | Machine learning-guided dual heuristics and new lower bounds for the refueling and maintenance planning problem of nuclear power plants[END_REF] counts with variables x i,i the number of colors used. Constraints [START_REF] Dupin | Parallel matheuristics for the discrete unit commitment problem with min-stop ramping constraints[END_REF] ensure that each vertex i is colored: either x i,i = 1 or there exists a representative i ′ < i such that x i ′ ,i = 1. Constraints [START_REF] Dupin | Matheuristics to optimize refueling and maintenance planning of nuclear power plants[END_REF] expresses incompatibility of neighbor vertices to have the same color, and that a variable x i ′ ,i = 1 implies that x i ′ ,i ′ = 1.

Standard DSATUR algorithm

DSATUR colors the vertices one after another, assigning the first color available or adding a new color. The order of traversal follows order ≽. Among the uncolored vertices, the selected vertex to color maximizes firstly the saturation and then the degree in a lexicographic way. Coloring a new vertex updates saturation, the order of traversal is not defined a priori, DSATUR is an adaptive greedy algorithm. Algorithm 1 writes the standard version of DSATUR, iterating with a partial coloring till the graph is fully colored.

DSATUR matheuristic variants

This section provides variants of DSATUR with matheuristic extensions of operators of Algorithm 1.

Initialization

Several initialization strategies can be used before processing the standard DSATUR constructive heuristic. Initialization consists of defining an initial partial coloring and computing the saturation table for the uncolored vertices. Several strategies can be provided:

• maxDeg: only one vertex is colored, one having the highest degree. Only the neighbors of this vertex have a saturation set to 1.

• col-n: one considers n vertices having the highest degrees, and color them solving to optimality ILP ( 9)-( 11) restricted to these n vertices. Note that n is a controlled parameter, that is set so that the ILP is solvable quickly. • clq: one find a large clique in the graph, and color this clique with different colors. Finding a maximum clique being NP-hard, this initialization requires a heuristic. In Appendix, a constructive matheuristic dealing with maximum clique problems of fixed size is presented. • clq-col-n: initialization strategy clq is firstly operated. Then col-n strategy is operated for n vertices having the highest saturation and degrees, considering the clique. Colors of these n vertices are either colors used for the clique, or new colors. If ILP formulations of section 2 can be used, fixing variables of the clique to their assignment or defining them as representatives of old colors, Section 3.2 provides an ILP formulation for VCP with existing colors.

Note that maxDeg is equivalent to standard DSATUR as written in Algorithm 1: no initialization induces that the first node to be selected in the loop is one with the highest degree. Many initialization are possible with maxDeg, there can be many vertices having the maximum degree. With col-n strategies, there is more depth in the initialization, to analyze if first iterations or standard DSATUR algorithm induce bad decisions. Initializing saturation with a clique does not induce a heuristic choice, it is an exact pre-processing, each color of the clique shall be different. The initialization of the saturation table is more advanced with a clique than with a single vertex.

Local optimization with larger neighborhoods

Let (c) be a partial k-coloring. Let k be the current number of colors in c, let C be the set of colored vertices, and U a subset of un-colored vertices:

C = {i ∈ I, c i > 0} (12) U ⊂ {i ∈ I, c i = -1} (13)
In this section, we define an ILP formulation to assign a color for vertices indexed in U while preserving the colors that are assigned in C. Following formulation hybridizes assignment-based ILP formulation for the existing colors, and the representative-based formulation for new colors. Binary variables x u,u ′ are defined here only for u < u ′ ∈ U , considering subset of edges

E U = {(v u , v u ′ )} u<u ′ ∈U .
Binary variables z u,l to assign previous colors are defined for u ∈ U and l ∈ [[1; k]] such that there is no neighbor u that have color l in (c). Variables z u,l are defined for all defined for u ∈ U for l ∈ K u where:

K u = {l ∈ [[1; k]], ∀i ∈ C, c i = l =⇒ ngb(i, j) = 0)} ( 14
)
It induces following ILP formulation to color the vertices indexed in U : Initialization can be any strategy defined in section 3.1. The remaining of the Algorithm simultaneously colors o vertices, solving an ILP [START_REF] Goudet | New bounds and constraint programming models for the weighted vertex coloring problem[END_REF] with o+r vertices and the previously assigned colors. In the standard version of DSATUR heuristic, we have o = 1 and r = 0. Having r > 0 ensures more depth in the local decision making with the possibility to reoptimize these variables after, as in [START_REF] Chandrasekharan | A constructive matheuristic approach for the vertex colouring problem[END_REF][START_REF] Dupin | Matheuristics to optimize refueling and maintenance planning of nuclear power plants[END_REF]. Having r = 0 could lead to threshold effects. To solve efficiently local optimization, there are o + r new vertices to color using ILP formulation [START_REF] Goudet | New bounds and constraint programming models for the weighted vertex coloring problem[END_REF], this parameter should be fixed according to the capability of the ILP solver to solve VCP problems for this size. Note that the r vertices that can be reoptimized are not necessarily chosen for the next iteration as the saturation table is updated with the o fixed colors.

min z u∈U x u,u s.t : z i,l + z i ′ ,l ⩽ 1 ∀(v i , v i ′ ) ∈ E U , ∀l ∈ [[1; k]] x u,i + x u,i ′ ⩽ x u,u ∀(v i , v i ′ ) ∈ E U , ∀u ∈ U, u ⩽ i i ′ ∈U :i ′ ⩽i x i ′ ,i + l∈Ku z i,l ⩾ 1 ∀u ∈ U (15)

General algorithm

Dual bounds

As in [START_REF] Boschetti | Matheuristics: survey and synthesis[END_REF][START_REF] Dupin | Machine learning-guided dual heuristics and new lower bounds for the refueling and maintenance planning problem of nuclear power plants[END_REF], some DSATUR matheuristics allow to have both lower and upper bounds, with dual bounds provided by varied relaxations [START_REF] Boschetti | Matheuristics: survey and synthesis[END_REF][START_REF] Dupin | Machine learning-guided dual heuristics and new lower bounds for the refueling and maintenance planning problem of nuclear power plants[END_REF]. Firstly, we recall that the cardinality of any clique gives a first lower bound, an optimal coloring (as any proper coloring) implies different colors for the clique. Algorithm 3 in Appendix A gives thus first dual bounds for VCP, after the first phase to initialize DSATUR with a clique.

After a clique initialization, any dual bounds of the ILP resolution of (15) assigning n = o + r colors, either in clq-col-n initialization or in the first iteration of in Algorithm 2 after clq initialization, is a dual bound for VCP, relaxing the constraints corresponding to the unoptimized nodes, and without any heuristic reduction of the original problem (which is not true once colors are fixed in a subset that is not a clique). As in [START_REF] Dupin | Machine learning-guided dual heuristics and new lower bounds for the refueling and maintenance planning problem of nuclear power plants[END_REF], dual bounds can be obtained with several relaxations computations of ILP ( 15): an exact ILP resolution of such ILP restricted with small values of n = o + r, larger values of N with computations of LP relaxation, or intermediate dual bounds with truncated ILP resolution with intermediate values of n. Note that such dual heuristics may take advantage of exact reduction techniques as in [START_REF] Goudet | New bounds and constraint programming models for the weighted vertex coloring problem[END_REF] to compute more efficiently dual bounds for smaller and equivalent VCP problems or to have a more relevant selection of the subset of n nodes considered in the relaxation.

Towards randomization and multi-start?

As in [START_REF] Melo | A matheuristic approach for the b-coloring problem using integer programming and a multi-start multi-greedy randomized metaheuristic[END_REF], one may extend Algorithm 2 for a multi-start initialization and use randomization. DSATUR standard algorithm can be randomized with the perturbation of the order of traversal of vertices, with the choice of the color to assign, not necessarily the first one. In Algorithm 2, the color to assign is optimized regarding the depth of traversal, so that such randomization does not make sense anymore. The depth in the order of traversal induces also few impact to have local perturbations of the order of traversal. It makes sense only to break ties for the threshold effects of the algorithm. A more promising randomization is to consider several large cliques to have several initialization.

For the computation of dual bounds, two randomization may be considered: the initial clique and also the selection of the n = o + r nodes for the second phase computation of dual bounds using the ILP solver.

Computational results

Computational experiments were processed using a computer Intel(R) Core(TM) i7-6700, 3.40GHz, running Linux Lubuntu 20.4, using up to 4 threads and 32 Gb of RAM memory. CPLEX version 20.1 was used for ILP resolution. Algorithms were coded in Julia programming language version 1.7.3, using the JuMP library version 1.1.1 to call ILP solvers and LightGraphs version 1.3.5 for graphs. Without specific precision, we use CPLEX 20.1 with its default parameter, except parameters CPX_PARAM_EPAGAP = 0.99999 to stop computation to optimality knowing the objective function is integer, a time limit of 30 seconds maximum for each computation (which was not reached in most of the following results). CPLEX allows to set optimization parameters to sizes n = 123 for clique depth search in Algorithm 3 and size n = o + r = 80 in Algorithm 2, to have partial ILP computations solvable to optimality in at most few seconds. Note that CBC can also be used with JuMP, to have an open source code, for these cases, we set n = 100 and o + r = 60.

Instances

For this study, we consider a subset of 53 DIMACS instances remowing instances that are easy for DSATUR, where DSATUR and the matheuristics gives the BKS that is proven optimal. These instances are highlighted in Appendix B, with their characteristics and their best known lower and upper bounds. For comparing primal heuristics, we used for this paper the instances without the exact pre-processign reduction from [START_REF] Goudet | New bounds and constraint programming models for the weighted vertex coloring problem[END_REF]. For the selected difficult instances, only 13 are reduced by [START_REF] Goudet | New bounds and constraint programming models for the weighted vertex coloring problem[END_REF], which could lead to easy instances, which was used for the computations of dual bounds. Original and reduced instances for VCP are available at https://github.com/Cyril-Grelier/gc wvcp cp.

Standard DSATUR with varied initialization

Table 1 presents for different initialization of DSATUR the total number of colors used for the 53 selected instances, the gaps to the Best Known Solutions (BKS), the number of instances where BKS is equaled, and the other columns are compared with the standard DSATUR Algorithm: #worse and #better counts the number of instances where the considered Algorithms have different values with BKS, respectively worse and better solutions, and quartiles are considered with the absolute gap from DSATUR to the corresponding algorithm, negative values means that the corresponding algorithm has better value than standard DSATUR, quartiles allows to appreciate the dispersion of the results.

Table 1 shows that col-n strategies are disappointing, leading to worse results in average than standard DSATUR. On the contrary, clq initialization improves significantly standard DSATUR. Using col-n strategies after clique initialization improves also significantly standard DSATUR. Note that for instance lr1000.1.col, a BKS is found by DSATUR and clq initialization, not for the other approaches, which explains this instance is considered in the selected pool of difficult instances for DSATUR.

It is interesting that clq-col-80 and clq improve DSATUR solutions on different instances, Considering the best results of both algorithms in the row "Best clq" of Table, as if we consider both approaches in parallel as in [START_REF] Dupin | Parallel matheuristics for the discrete unit commitment problem with min-stop ramping constraints[END_REF][START_REF] Dupin | Matheuristics and Column Generation for a Basic Technician Routing Problem[END_REF], it provides an additional significative improvement. In the row, "Best clq+DSATUR" we consider the best result including also DSATUR, to analyze the complementarity with the original approach. A very slight improvement is observed, as well as considering all the approaches in "Best+DSATUR" row, or removing only the . DSATUR standard approach in row "Best-DSATUR". These last results highlight that for three instances, le450_5b; queen11_11 and queen15_15, none of the other initialization improves or equals standard DSATUR. This section validates to consider both clq-col-80 and clq strategies, in a multi-start of parallel heuristic, and the power of using cliques for DSATUR variants. If the solutions of standard DSATUR have been improved, the gaps from the BKS remain very significant.

DSATUR with larger local optimization

Table 2 has the same shape as Table 1 to compare DSATUR extended matheuristics to the standard Algorithm 1. Parameters o, r are the ones in Algorithm 2, standard version of DSATUR, implemented with Algorithm 1, corresponds to o = 1 and r = 0. This allows to analyze the impact of a larger depth in local optimization and the part of vertices to reoptimize for a better efficiency. (60,[START_REF] Melo | A matheuristic approach for the b-coloring problem using integer programming and a multi-start multi-greedy randomized metaheuristic[END_REF]} improves DSATUR standard algorithm, this is not the case with (o, r) = (80, 0). Coherently with [START_REF] Chandrasekharan | A constructive matheuristic approach for the vertex colouring problem[END_REF][START_REF] Dupin | Matheuristics to optimize refueling and maintenance planning of nuclear power plants[END_REF], it is important to have a significant part of variables that can be reoptimized to avoid bad choices due to threshold effects. A drawback of increasing r value (and decreasing o value to keep value o + r stable), is that computation times are increasing. (o, r) = (40, 40) is a good compromise between solution quality and computation time.

Initializing with clq provides again the best results, (o, r) = (40, 40) improves significantly DSATUR with both the standard and the clique initialization. Combining (o, r) = (40, 40) and standard DSATUR construction (o, r) = (1, 0) allows an additional improvement. This highlights that standard DSATUR algorithm has good properties, that can be broken with more depth in local optimization.

Coherently with [START_REF] Dupin | Matheuristics and Column Generation for a Basic Technician Routing Problem[END_REF], using larger neighborhoods in greedy constructive heuristics improves in average the solution quality. However, even with an ensemble of such constructive heuristics that can be computed in parallel, a significant gap remains to the BKS.

Dual bounds

In Appendix B, best known lower bounds reported, mainly after [START_REF] Held | Safe lower bounds for graph coloring[END_REF] and in some cases the chromatic number is known by construction or specific reasoning. To compute the dual bounds, we used the exact reduction from [START_REF] Goudet | New bounds and constraint programming models for the weighted vertex coloring problem[END_REF] for the 13 instances where a reduction is obtained, as shown in Tables 4 and5.

For the DIMACS instances, Algorithm 3 is very efficient to find large cliques, and gives already the maximum clique size which is also the chromatic number for 28 out of the 53 selected instances with parameter n ∈ {100, 125}. This often occurs also for the easy instances that were removed from the dataset for this study as easy instances for DSATUR. We removed thus these instances for results of dual bounds, these case being specific and easy to compute optimal dual bounds. Note that having n = 100 or n = 125 produced the same results with Algorithm 3, computations with n = 125 are slightly longer. For the 25 remaining instances, we report dual bounds and computation times obtained after a clique computation with Algorithm 3 computations with n = 100 with following parameters to analyze the compromise between the number of vertex to consider in the ILP (15):

• n = o + r = 80 and a time limit of 300s, with bounds at the root node of the B&B tree and in truncated resolution time, or the optimal value. • n = o + r = 125 and a time limit of 900s, with bounds at the root node of the B&B tree . • n = o + r = 200 and a time limit of 3600s for B&B tree search. Table 3 shows that the computation of dual bounds with ILP (15) induced improvements of the cliques given in input for most of the instances. For latin_square instance, a clique of size 90 is found easily, it is the actual best Lower Bound known (BLBK) for VCP, the ILP computations of improved dual bounds do not improve this bound. Note that in experiments of [START_REF] Held | Safe lower bounds for graph coloring[END_REF], computations could take around 30 days to have dual bounds for very difficult problems, especially for instances C2000.5 and C4000.5. For such large and difficult instance, where the best cliques knows are far from the BKS and BLBK, our dual bounds are quite limited, which is also the cases for [START_REF] Furini | Exact weighted vertex coloring via branch-and-price[END_REF].

Table 3 shows it is preferable and computable to tackle problems with n = o + r = 200. The higher the value of n is, the higher the optimal solution of ILP [START_REF] Goudet | New bounds and constraint programming models for the weighted vertex coloring problem[END_REF], and truncated ILP resolution remains efficient to computer better dual bounds than the ones computed optimally with smaller subproblems. For three instances, namely dsjr500.5 dsjc125.1 and dsjc125.9, optimal lower bounds are found quickly. Globally, more improvements are observed on dense graphs (suffixed by .5 and .9 indicating the density) than on sparse ones (suffixed by .1 ), the more difficult instances were no improvement is observed are sparse graphs. Note that half of ILP [START_REF] Goudet | New bounds and constraint programming models for the weighted vertex coloring problem[END_REF] with n = 200 are solved to optimality in less than one hour, sometimes very quickly, clique initialization may be helpful to speed up such ILP computations, whereas for other instances some ILP computations of size n = 125 are not solvable in one hour.

Conclusions and perspectives

This paper studied matheuristic variants of DSATUR, to improve its standard version, and also to help understanding the strengths and weaknesses of this well-known heuristic. Initializing DSATUR with a large clique, using a simple greedy matheuristic, is a very significant improvement of the standard initialization with one vertex of maximal degree. Having larger optimization in the greedy construction is efficient when some vertices can be re-optimized to avoid threshold effects. However, improvement of DSATUR is slight, significant gaps to BKS remain using DSATUR constructive matheuristics. Dual bounds are also provided, highlighting also the interest of cliques for DSATUR. With a newly introduced ILP formulation, dual bounds implied by cliques can be improved in short and long computation times.

These results offer several perspectives. Firstly, exact version of DSATUR [START_REF] Segundo | A new dsatur-based algorithm for exact vertex coloring[END_REF][START_REF] Furini | An improved dsatur-based branch-andbound algorithm for the vertex coloring problem[END_REF] could be improved using cliques for branching in the tree search algorithm, dual and primal heuristics can be used and parametrized to prune some nodes in this tree search. Secondly, dual bound can be improved using other exact techniques for dual bounds, as [START_REF] Furini | Exact weighted vertex coloring via branch-and-price[END_REF][START_REF] Cornaz | Solving vertex coloring problems as maximum weight stable set problems[END_REF][START_REF] Held | Safe lower bounds for graph coloring[END_REF]. Thirdly, perspectives are to extend similarly RLF as matheuristics.

Algorithm 1 :

 1 Standard DSATUR algorithm Input: G = (V, E) a non-empty and non-oriented graph Initialization: define partial coloring c with ci := -1 for all i ∈ I define saturation table S with Si := ∅ for all i ∈ I initialize set U := V , and color k := 0 while U ̸ = ∅ find u ∈ U , a maximum of ≽ in U . if |Su| = k then k:=k + 1 // a new color is added compute ci := min Su // assign color to u remove u from U for all i ∈ δu ∩ U , Si = Si ∪ {ci} // update saturation end while return color k and (c) a k-coloring of G

Algorithm 2 :

 2 Matheuristic DSATUR variants Input: G = (V, E) a non-empty and non-oriented graph Parameters: • an initialization strategy S (from section 3.1) ; • o ∈ N, o > 1 ; • r ∈ N. Initialization: initialize colored set C, and color k with strategy S. initialize W := V \ C. update partial coloring c and saturation table S with strategy S. while W ̸ = ∅ sort W with order ≽. define U1 as the o first elements after sorting. define U2 as the elements of rank o + 1 and min(|W |, o + r) after sorting. solve ILP (15) with C and U = U1 ∪ U2. k := k + OP T where OP T is the optimal value of the last ILP. if o + r ⩽ |W | then U1 = U end if set W := W \ U1 assign colors cu of the ILP for u ∈ U1 end while return color k and (c) a k-coloring of G Algorithm 2 is a general version for an extended DSATUR matheuristic.

  -based ILP model Having k a maximum number of colors, following ILP model minimizes the number of colors to cover G, if the chromatic number is at most k, otherwise the infeasibility of this ILP proves there exists no k-coloring for G. Two types of binary variables are used. On one hand, assignment variables z i,c ∈ {0, 1} are defined for i ∈ I and c ∈ [[1; k]], z i,c = 1 if and only if vertex v i is assigned to color c. On the other hand, availability variables

	y c ∈ {0, 1} are defined c ∈ [[1; k]], y c = 1 if and only if color c is used. It gives
	rise to following ILP:
	k
	min	y c
	c=1

Table 1 .

 1 Comparison of DSATUR matheuristics with different initialization of saturation table. Results parallelizing several strategies are also provided.

		#colors gap #BKS #worse #better Q1 Q2 Q3
	maxDeg	3240 32.03 %	1	0	0	0 0 0
	col-60	3251 32.48 %	1	19	16	-1 0 1
	col-80	3250 32.44 %	2	20	16	-1 0 1
	clq-col-80	3214 30.97 %	2	18	17	-1 0 1
	clq	3209 30.77 %	4	13	19	-1 0 0
	Best clq	3181 29.63 %	6	7	26	-1 0 0
	Best clq+DSATUR 3174 29.34 %	6	0	26	-1 0 0
	Best-DSATUR	3163 28.89 %	6	3	34	-2 -1 0
	Best+DSATUR	3160 28.77 %	6	0	34	-2 -1 0
	BKS	2454 0.00 %	53	0	52	-14 -5 -3

Table 2 .

 2 Comparison of DSATUR matheuristics with different initialization of saturation, and values of optimization parameters o and r in Algorithm 2

	Init satur	o r #colors gap #BKS #worse #better Q1 Q2 Q3
	maxDeg	1 0 3240 32,03 %	1	0	0	0 0 0
	col-80	1 0 3250 32,44 %	2	20	16	-1 0 1
	col-80	20 60 3181 29,63 %	6	12	30	-3 -1 0
	col-80	40 40 3218 31,13 %	5	20	26	-2 0 1
	col-80	80 0 3322 35,37%	2	35	13	0 1 2
	clq	1 0 3209 30,77 %	4	13	19	-1 0 0
	clq	40 40 3155 28,57% 10	9	32	-3 -1 0
	Best Clq	3134 27,71 % 10	4	37	-3 -1 0
	Best-DSATUR	3125 27,34 % 10	3	40	-3 -2 -1
	Best+DSATUR	3122 27,22 % 10	0	40	-3 -2 -1
	BKS	2454 0,00 %	53	0	52	-14 -5 -3

Table 2

 2 does not provide results with maxDeg initialization, first iteration of Algorithm 2 induce a similar saturation than using col-n initialization.

	Using clq initialization, results are very stable considering parameters value
	(o, r) ∈ {(20, 60); (40, 40); (60, 20); (80, 0)}. Using col-80 initialization, setting
	parameter values (o, r) ∈ {(20, 60); (40, 40);

Table 3 .

 3 Comparaison of the lower bounds obtained by initial clique computation and ILP refinements to the BKS and BKLB, LB are reported as well as computation times. For ILP refinements, the time to compute the clique is not counted, these times are additional time to improve the lower bounds given by clique

		UB LB		LB			t
		BKS BKLB clq n = 125 n = 200 clq n = 125 n = 200
	C2000.5	145	99	15	20	21	185 163	3600
	C4000.5	259 107 17	21	22	252	99	3600
	dsjc125.1	5	5	4	5	5	0,2 81,8	118
	dsjc125.5	17	17	10	14	14	14 131,4	3600
	dsjc125.9	44	44	34	43	44	33	1	1
	dsjc250.1	8	7	4	6	5	0,7 186,4	3600
	dsjc250.5	28	26	12	16	17	160 206	3600
	dsjc250.9	72	71	41	56	70	53	1,5	105
	dsjc500.1	12	9	5	5	5	5	4	3600
	dsjc500.5	48	43	13	17	19	167	61	3450
	dsjc500.9	126 123 51	65	79	50	0,4	274
	dsjc1000.1	20	10	6	6	6	38	8,6	3600
	dsjc1000.5	83	73	14	19	20	175 172,6	3600
	dsjc1000.9	222 215 59	73	86	80	2,3	15
	dsjr500.1c	85	85	76	77	79	47	3	11
	dsjr500.5	122 122 114 122	122	5	0,6	10
	flat300 26 0 26	26	11	15	16	167 217	3600
	flat300 28 0 28	28	12	15	16	160 259	3600
	flat1000 50 0 50	50	13	17	19	175 186	3600
	flat1000 60 0 60	60	13	17	19	178 135	3600
	flat1000 76 0 76	76	14	18	19	166 179	3600
	latin square 97	90	90	90	90	32	0,2	12
	r1000.1c	98	96	87	88	88	143	73	9
	r1000.5	234 234 213 214	220	81	8,5	19
	TOTAL	1965 1716 928 1039	1101
	Average						95	87	2033

Appendix A: matheuristic to find large cliques and stables

This appendix present the matheuristic that computes a large clique as initialization of Algorithm 2. To ease presentation, we present the matheuristic in Algorithm 3 for the Maximum Independent Set (MIS) problem applied to the complementary graph of G. Indeed, it is equivalent for a subset V to be a clique in the graph G and an independent (or stable) set in the complementary graph. Algorithm 3 computes iteratively an independent set based on MIS of fixed size n. Defining with U 1 a subset of V , a maximum independent set in U 1 can be computed using the following ILP formulation, where binary variables z v ∈ {0, 1}, are defined with

Algorithm 3 is an adaptive greedy algorithm: once vertices are added in the current independent set, the next candidate vertices are chosen with the minimum degrees in the updated graph, removing neighbors of selected points that cannot be added in the current stable.

Appendix B: Selected instances and their characteristics