A Distribution Evolutionary Algorithm for the Graph Coloring Problem

Besides the simulated annealing [9, 10] and the variable neighborhood search [11], the tabu search (TS) is one of the most popular individual-based metaheuristics applied to solve the GCPs [12]. Porumbel et al. [13] improved the performance of TS by evaluation functions that incorporates the structural or dynamic information in addition to the number of conflicting edges. Blöchliger and Zufferey [14] proposed a TS-based constructive strategy, which constructs feasible but partial solutions and gradually increases its size to get the optimal color assignment of a GCP. Hypothesizing that high quality solutions of GCPs could be grouped in clusters within spheres of a specific diameter, Porumbel et al. [15] proposed two improved TS variants using a learning process and a tree-like structure of the connected spheres. Assuming that each vertex only interacts with a limited number of components, Galán [16] developed a decentralized coloring algorithm, where colors of vertexes are modified according to those of the adjacent vertexes to iteratively reduce the number of edge conflicts. Sun et al. [17] established a solution-driven multilevel optimization framework for GCP, where an innovative coarsening strategy that merges vertexes based on the solution provided by the TS, and the uncoarsening phase is performed on obtained coarsened results to get the coloring results of the original graph. To color vertexes with a given color number $k$, Peng et al. [18] partitioned a graph into a set of connected components and a vertex cut component, and combined the separately local colors by an optimized maximum matching based method.

Since a probability model can provide a bird’s-eye view for the landscape of optimization problem, Zhou et al. [19, 20] proposed to enhance the global exploration of individual-based metaheuristics by the introduction of probability models. They deployed a probabilistic model for the colors of vertexes, which is updated with the assistance of a reinforcement learning technology based on discovered local optimal solutions [19]. Moreover, they improved the learning strategy of probability model to develop a three-phase local search, that is, a starting coloring generation phase based on a probability matrix, a heuristic coloring improvement phase and a learning based probability updating phase [20].

Population-based iteration mechanisms are incorporated to the improve exploration abilities of metaheuristics as well. Hsu et al. [21] proposed a modified turbulent particle swarm optimization algorithm for the planar graph coloring problem, where a three-stage turbulent model is employed to strike a balance between exploration and exploitation. Hernández and Blum [22] dealt with the problem of finding valid graphs colorings in a distributed way, and the assignment of different colors to neighboring nodes is asynchronously implemented by simulating the calling behavior of Japanese tree frogs. Rebollo-Ruiz and M. Graña [23] addressed the GCP by a gravitational swarm intelligence algorithm, where nodes of a graph are mapped to agents, and its connectivity is mapped into a repulsive force between the agents corresponding to adjacent nodes. Based on the conflict matrices of candidate solutions, Zhao et al. [24] developed a dimension-by-dimension update method, by which a discrete selfish herd optimizer was proposed to address GCPs efficiently. Aiming to develop an efficient parameter-free algorithm, Chalupa and Nielsen [25] proposed to improve the global exploration by a multiple cooperative searching strategy. For the four-colormap problem, Zhong et al. [26] proposed an enhanced discrete dragonfly algorithm that performs a global search and a local search alternately to color maps efficiently.

The incorporation of probability models are likewise employed to improve the performance of population-based metaheuristics. Bui et al. [27] developed a constructive strategy of coloring scheme based on an ant colony, where an ant colors just a portion of the graph unsing only local information. Djelloul et al. [28] took a collection of quantum matrices as the population of the cuckoo search algorithm, and an adapted hybrid quantum mutation operation was introduced to get enhanced performance of the cuckoo search algorithm.

The population-based metaheuristics can be further improved by hybrid search strategies. Paying particular attention to ensuring the population diversity, Lü and Hao [29] proposed an adaptive multi-parent crossover operator and a diversity-preserving strategy to improve the searching efficiency of evolutionary algorithms, and proposed a memetic algorithm that takes the TS as a local search engine. Porumbel et al. [30] developed a population management strategy that decides whether an offspring should be accepted in the population, which individual needs to be replaced and when mutation is applied. Mahmoudi and Lotfi [31] proposed a discrete cuckoo optimization algorithm for the GCP, where a neighborhood search in radius of the lay egg causes the algorithm hardly trapped in local minimum and producing new eggs. Accordingly, it provides a good balance between diversification and centralizing.

Wu and Hao [32] proposed a preprocessing method that extracts large independent sets by the TS, and the memetic algorithm proposed by Lü and Hao [29] was employed to color the residual graph. For the chromatic problem, Douiri and Elbernoussi [33] initialized the color number by the coloring result of a heuristic algorithm and generated the initial population of genetic algorithm (GA) by finding a maximal independent set approximation of the investigated graph.

Bessedik et al. [34] addressed the GCPs within the framework of the honey bees optimization, where a local search, a tabu search and an ant colony system are implemented as workers and queens are randomly generated. Mirsaleh and Meybodi [35] proposed a Michigan memetic algorithm for GCPs, where each chromosome is associated to a vertex of the input graph. Accordingly, each chromosome is a part of the solution and represents a color for its corresponding vertex, and each chromosome locally evolves by evolutionary operators and improves by a learning automata based local search. Moalic and Gondran [36] integrated a TS procedure with an evolutionary algorithm equipped with the greedy partition crossover, by which the hybrid algorithm can performs well with a population consisting of two individuals. Silva et al. [37] developed a hybrid algorithm iColourAnt, which addresses the GCP using an ant colony optimization procedure with assistance of a local search performed by the reactive TS.

A large number of works have been reported to improve the general performance of EDAs. To improve the general precision of a distribution model, Shim et al. [38] modelled the restricted Boltzmann machine as a novel EDA, where the probabilistic model is constructed using its energy function, and the $k$-means clustering was employed to group the population into small clusters. Approximating the Boltzmann distribution by a Gaussian model, Valdez et al. [39] proposed a Boltzmann univariate marginal distribution algorithm, where the Gaussian distribution obtains a better bias to sample intensively the most promising regions. Considering the multivariate dependencies between continuous random variables, PourMohammadBagher et al. [40] proposed a parallel model of some subgraphs with a smaller number of variables to avoid complex approximations of learning a probabilistic graphical model. Dong et al. [41] proposed a latent space-based EDA, which transforms the multivariate probabilistic model of Gaussian-based EDA into its principal component latent subspace of lower dimensionality to improve its performance on large-scale optimization problems.

To enhance the local exploitation of an EDA, Zhou et al. [42] suggested to combine an estimation of distribution algorithm with cheap and expensive local search methods for making use of both global statistical information and individual location information. Considering that the random sampling of Gaussian EDA usually suffers from the poor diversity and the premature convergence, Dang et al. [43] developed an efficient mixture sampling model to achieve a good tradeoff between the diversity and the convergence, by which it can explore more promising regions and utilize the unsuccessful mutation vectors.

The performance of EDA can also improved by designing tailored update strategies of the probability model. To address the multiple global optima of multimodal problem optimizations, Pẽna et al. [44] introduced the unsupervised learning of Bayesian metwokrs in EDA, which makes it able to model simultaneously the different basins represented by the selected individuals, whereas preventing genetic drift as much as possible. Peng et al. [45] developed an explicit detection mechanism of the promising areas, by which function evaluations for exploration can be significantly reduced. To prevent the Gaussian EDAs from premature convergence, Ren et al. [46] proposed to tune the main search direction by the anisotropic adaptive variance that is scaled along different eigendirections based on the landscape characteristics captured by a simple topology-based detection method. Liang et al. [47] proposed to archive a certain number of high-quality solutions generated in the previous generations, by which fewer individuals are needed in the current population for model estimation. In order to address the mixed-variable newsvendor problem, Wang et al. [48] developed a histogram model-based estimation of distribution algorithm, where an adaptive-width histogram model is used to deal with the continuous variables and a learning-based histogram model is applied to deal with the discrete variables. Liu et al. [49] embedded within the search procedure a learning mechanism based on an incremental Gaussian mixture model, by which all new solutions generated during the evolution are fed incrementally into the learning model to adaptively discover the structure of the Pareto set of an MOP.

Let $n=|V|$ be the vertex number of a graph $G=(V,E)$. An assignment of vertexes with $k$ colors can be represented by an integer vector $\mathbf{x}=(x_{1},\dots,x_{n})$, where $x_{j}$ denotes the assigned color of vertex $v_{j}$. Then, the $k$-coloring problem can be modelled as a minimization problem

While $\delta(j_{1},j_{2})=1$, the adjacent vertexes $v_{j_{1}}$ and $v_{j_{2}}$ are assigned with the same color, and $(v_{j_{1}},v_{j_{2}})$ is called a conflicting edge. Accordingly, the objective values $f_{k}(\mathbf{x})$ represents the total conflicting number of the color assignment $\mathbf{x}$. While $G$ is $k$-colorable, there exists an optimal partition $\mathbf{x}^{*}$ such that $f_{k}(\mathbf{x}^{*})=0$, and $\mathbf{x}^{*}$ is called a legal $k$-color assignment of graph $G$. Thus, the chromatic problem is modelled as

where $\mathbf{x}^{*}$ represents a legal $k$-color assignment that is an optimal color assignment of problem (1).

red

Different from the ACO and the EDA, the QEA employs a quantum matrix for probabilistic modelling of the solution space, modelling the probability distribution of a binary variable $bx$ by a Q-bit $(\alpha,\beta)^{T}$ satisfying $|\alpha|^{2}+|\beta|^{2}=1$ [50]. That is, $|\alpha|^{2}$ and $|\beta|^{2}$ give the probabilities of $bx=1$ and $bx=0$, respectively. Accordingly, the probability distribution of an $n$-dimensional binary vector $\mathbf{bx}=(bx_{1},\dots,bx_{n})$ is represented by

where $|\alpha|_{j}^{2}+|\beta|_{j}^{2}=1$, $j=1,\dots,n$. Then, the value of $bx_{j}$ can be obtained by sampling the probability distribution $(|\alpha|_{j}^{2},|\beta|_{j}^{2})$.

The probability distribution of binary variable is modified in the QEA by the Q-gate, a $2\times 2$ orthogonal matrix

Premultiplying $\vec{q}_{j}$ by $U(\Delta\theta)$, the probability distribution of $bx_{j}$ is modified as

To coloring a graph of $n$ vertexes by $k$ colors, Djelloul et al. [28] modelled the probability distribution of color assignment by a quantum matrix

where

In this way, the Q-gate can be deployed for each unit vector $(q_{2i-1,j},q_{2i,j})^{T}$ to regulate the distribution of color assignment.

Since the unit vector $(q_{2i-1,j},q_{2i,j})^{T}$ models each candidate color $i$ of vertex $j$ independently, the sampling process would lead to multiple assignments for vertex $j$, and an additional repair strategy is needed to get a feasible $k$-coloring assignment [28]. To address this defect, we propose to model the color distribution of vertex $j$ by a $k$-dimensional unit vector $\vec{q}_{j}$, and present the color distribution of $n$ vertexes as

where $\vec{q}_{j}$ satisfies

The distribution model confirmed by (3) and (4) incorporates the advantages of models in EDAs and QEAs.

• 1.

  An feasible $k$-coloring of $n$ vertexes can be achieved by successively sampling $n$ columns of $\mathbf{q}$.

• 2.

  The update of probability distribution can be implemented by both orthogonal transformations performed on column vectors of $\mathbf{q}$ and direct manipulations of components that do not change the norms of columns vectors ¹¹1Details for the update process are presented in Section 4.4..

red As presented in Algorithm 1, DEA-PPM is implemented as two nested loops: the inner loop addressing the $k$-coloring problem and the outer loop decreasing $k$ to get the chromatic number $\chi(G)$. Based on a distribution population $\mathbf{Q}(t)=(\mathbf{q}^{[1]}(t),\dots,\mathbf{q}^{[np]}(t))$ and the corresponding solution population $\mathbf{P}(t)=(\mathbf{x}^{[1]}(t),\dots,\mathbf{x}^{[np]}(t))$, it starts with the initialization of the color number $k$, which is then minimized by the outer loop to get the chromatic number $\chi(G)$.

Each iteration of the outer loop begins with the iterative vertex removal (IVR) strategy [51], by which the investigated graph $G$ could be transformed into a reduced graph $G^{\prime}=(V^{\prime},E^{\prime})$, and the complexity of the coloring process could be reduced as well. Then, Lines 5-11 of Algorithm 1 initialize a distribution population $\mathbf{Q}(0)$ and the corresponding solution population $\mathbf{P}(0)$ for $G^{\prime}$. Once $G^{\prime}$ is colored by Lines 12-20 of Algorithm 1, DEA-PPM recovers the obtained color assignment $\mathbf{x}^{*}_{G^{\prime}}$ to get an color assignment $\mathbf{x}^{*}_{G}$ of $G$, and $\mathbf{Q}(t)$ as well as $\mathbf{P}(t)$ is archived for the inherited initialization performed at the next generation. Repeating the aforementioned process until the termination-condition 1 is satisfied, DEA-PPM returns a color number $k$ and the corresponding color assignment $\mathbf{x}^{*}_{G}$.

After the initialization of $\mathbf{x}^{*}_{G^{\prime}}$, $\mathbf{p}_{1}$, $\mathbf{p}_{2}$ and $t$, the inner loop tries to get a legal $k$-coloring assignment for the reduced graph $G^{\prime}$ by evolving both the distribution population $\mathbf{Q}(t)$ and the solution population $\mathbf{P}(t)$. It first performs the orthogonal exploration on $\mathbf{Q}(t)$ to generate $\mathbf{Q}^{\prime}(t)$, and then, generates an intermediate solution population $\mathbf{P}^{\prime}(t)$, which is further refined to get $\mathbf{P}(t+1)$. Meanwhile, $\mathbf{Q}(t+1)$ is generated by refining $\mathbf{Q}^{\prime}(t)$. The inner loop repeats until the termination-condition 2 is satisfied.

The outer loop of DEA-PPM is implemented only once for the $k$-coloring problem. To address the chromatic problem, the termination condition 1 is satisfied if the chromatic number has been identified or the inner loop fails to get a legal $k$-coloring assignment for a given iteration budget. The termination condition 2 is met while a legal $k$-coloring assignment is obtained or the maximum iteration number is reached.

For the $k$-coloring problem, DEA-PPM initializes the color number $k$ by a given positive integer. While it is employed to address the chromatic number problem, we set $k=\Delta G+1$²²2Here $\Delta G$ is the maximum vertex degree of graph $G$. because an undirected graph $G$ is sure to be $(\Delta G+1)$-colorable [52].

To reduce the time complexity of the $k$-coloring algorithm, Yu et al. [51] proposed an iterative vertex removal (IVR) strategy to reduce the size of the investigated graph. By successively removing vertexes with degrees less than $k$, IVR generates a reduced graph $G^{\prime}=(V^{\prime},E^{\prime})$, and put the removed vertexes into a stack $S$. In this way, one could get a graph $G^{\prime}$ where degrees of vertexes are greater than or equal to $k$, and its size could be significantly smaller than that of $G$.

While a $k$-coloring assignment $\mathbf{x}^{*}_{G^{\prime}}$ is obtained for the reduced graph $G^{\prime}$, the inverse recovery (IR) operation is implemented by recovering vertexes in the stack $S$. The IR process is initialized by assigning any legal color to the vertex at the top of $S$. Because the IVR process removes vertexes with degree less than $k$, the IR process can get all recovered vertexes colored without conflicting. An illustration for the implement of the IVR and the IR is presented in Fig. 1.

Depending on the iteration stage of DEA-PPM, the initialization of populations is implemented by the uniform initialization or the inherited initialization.

At the beginning, the uniform initialization generates $np$ individuals of $\mathbf{Q}(0)$ as

$\mathbf{P}(0)$ are generated by sampling model (5) $np$ times.

While $gen>0$, the inherited initialization gets $\mathbf{Q}(0)$ and $\mathbf{P}(0)$ with the assistance of the distribution populations $\mathbf{Q}$ and the solution population $\mathbf{P}$ archived at the last generation. The graph $G$ has been colored with $k$ colors, and it is anticipated to get a legal $k-1$-color assignment. To get an initial color assignment of $k-1$ colors, a color index $l_{m}$ that corresponds to the minimum vertex independent set is identified for $\mathbf{x}^{[i]}=({x}^{[i]}_{1},\dots,{x}^{[i]}_{n})\in\mathbf{P}$. Then, we get the initial color assignment $\mathbf{y}^{[i]}=({y}^{[i]}_{1},\dots,{y}^{[i]}_{n})$ by

Meanwhile, delete the $l_{m}$-th row of $\mathbf{q}^{[i]}$, and normalize its columns to get an initial distribution $\mathbf{r}^{[i]}$ corresponding to $k-1$ colors. Details of the inherited initialization are presented in Algorithm 2.

Based on the distribution model defined by (3) and (4), DEA-PPM performs the orthogonal exploration on individuals of $\mathbf{Q}(t)$ to explore the probability space. Moreover, distribution individuals of $\mathbf{Q}(t)$ are refined by an exploitation strategy or a disturbance strategy.

To search the solution space efficiently, DEA-PPM generates a solution population by sampling with inheritance, and then, refines it using a multi-parent crossover operation followed by the TS search proposed in Ref. [55].

By setting $k$ as the best known color numbers of the benchmark problems, preliminary experiments for the $k$-coloring problem show that the performance of DEA-PPM is significantly influenced by the population size $np$, the regulation parameter $\alpha$ and the maximum iteration budget $Iter_{max}$ of the TS. Then, we first demonstrate the univariate influence of parameters by the one-way analysis of variance (ANOVA), and then, perform a descriptive comparison to get a set of parameter for further numerical investigations. \colorblue The benchmark instances selected for the parameter study are the $k$-coloring problems of DSJC500.5, flat300_28_0, flat1000_50_0, flat1000_76_0, le450_15c, le450_15d.

red

In DEA-PPM, the evolution of distribution population is implemented by the orthogonal exploration strategy and the exploitation strategy. We try to validate the positive effects of these strategies in this section.

To validate the efficiency of the orthogonal exploration strategy, we compare two variants, the DEA-PPM with orthogonal exploration (DEA-PPM-O) and the DEA-PPM without orthogonal exploration (DEA-PPM-N), and show in Fig. 4 the statistical results of running time for $k$-coloring of the easy benchmark problems (DSJC500.1 ($k=12$), le450_15c ($k=15$), led450_15d ($k=15$)) and the hard benchmark problems(DSJC500.5 ($k=48$), DSJC1000.1 ($k=20$), DSJC1000.9 ($k=226$)). The box plots imply that with the employment of the orthogonal exploration strategy, DEA-PPM-O performs generally better than DEA-PPM-N, resulting in smaller values of the median value, the quantiles and the standard deviations of running time.

red The positive impact of exploitation strategy is verified by comparing the DEA-PPM with exploitation (DEA-PPM-E) with the variant without exploitation (DEA-PPM-W), and the box plots of running time are included in Fig. 5. It is demonstrated that DEA-PPM-E generally outperforms DEA-PPM-W in terms of the median value, the quantiles and the standard deviation of running time.

blue Besides the statistical comparison regarding the exact values of running time, we perform a further comparison by the Wilcoxon rank sum test with a significance level of 0.05, where the statistical test is based on the sorted rank of running time instead of its exact values. The results are included in Tab. 4, where “P” is the p-value of hypothesis test. For the test conclusion “R”, “+”, “-” and “$\sim$” indicate that the performance of DEA-PPM is better than, worse than and incomparable to that of the compared variant, respectively. The results demonstrate that DEA-PPM outperforms DEA-PPM-N and DEA-PPM-W on two instances, and is not inferior to them for all benchmark problems. It further validates the conclusion that both the exploration strategy and the exploitation strategy significantly improve the performance of DEA-PPM.

To demonstrate the competitiveness of DEA-PPM, we perform numerical comparison for the chromatic problem and the $k$-coloring problem with SDGC [16], MACOL [29], SDMA [17], PLSCOL [20], and HEAD [36], the parameter settings of which are presented in Tab. 5. If an algorithm cannot address the chromatic problem or the $k$-coloring problem in 3600 seconds, a failed run is recorded by the running time of 3600 seconds.