# How can I speed up the classic GA for graph coloring?

I'm trying to compute the chromatic number of this graph (which is 28):

g = Import@"http://www.info.univ-angers.fr/pub/porumbel/graphs/dsjc250.5.col";


My genetic algorithm is getting stuck at an upper bound of 38 vertex colors:

In[] := Timing @ GAColor[g, 10, 20, 3]
Out= {19.178072, {38, {28, 16, ...., 3, 22}}}


I've written the general GA implementation, but I'm using naive recombination and mutation, and my mathematica code is slow. My question is, how could I improve on this with more clever choices of Combine[] and Mutate[], as well as faster code, in the general? I'm by no means an expert here, so I'm sure there are many possible improvements both theoretically and algorithmically...

GAColor[g_Graph, PopulationSize_Integer:100, NumberOfGenerations_Integer:10, NumberOfMutants_Integer:0, mutationRadius_Integer:Automatic] :=  Module[
{NumberOfVertices = VertexCount @ g, NumberOfBreeders, PermuteColorClasses, MutationRadius, Combine, Mutate, PopulationStep,
InitializePopulation, InitalPopulation, Generations, BestFitness = \[Infinity], BestColoring, GenerationsFitness, Chromatize, Adjacencies, result},

NumberOfBreeders = PopulationSize - NumberOfMutants;

PermuteColorClasses[colors_, n_:1] := Module[{p},
p = {#, Flatten @ Position[colors, #]}& /@ Union[colors];
p[[All, 1]] = RandomSample[p[[All,1]]];
ReplacePart[ConstantArray[0, Length[colors]], Flatten[Rule @@@ Thread[Reverse@#]& /@ p]]
];

Chromatize[colorVector_] := Module[{f, h, min, co = colorVector},
f = Function[{c,v},
ReplacePart[c, v -> With[{ncols = c[[Adjacencies[[v]]]]},
For[min = 1, MemberQ[ncols, min], min++]; min]
]];
h = Function[{c}, Fold[f, c, RandomSample[Range[NumberOfVertices]]]];
FixedPoint[h, co]
];

Combine[colorVector1_, colorVector2_] := MapThread[RandomChoice[{#1, #2}]&, {colorVector1, colorVector2}];
Mutate[colorVector_, mr_] := Permute[colorVector, RandomPermutation[mr]];

PopulationStep[population_, NumberOfBreeders_] := Module[
{fitness = Max /@ population, breeders, children, mutants},
With[{min = Min[fitness]}, If[min < BestFitness, BestFitness = min]];
breeders = RandomChoice[fitness -> population, NumberOfBreeders];
children = Chromatize /@ Table[Combine @@ RandomChoice[breeders, 2], {NumberOfBreeders}
];
mutants = Mutate[#, MutationRadius]& /@ RandomChoice[breeders, NumberOfMutants];
Join[children, mutants]
];

InitalPopulation = With[{color = Chromatize[RandomSample[Range @ NumberOfVertices]]},
Table[PermuteColorClasses[color], {PopulationSize}]
];

Generations = NestList[PopulationStep[#, NumberOfBreeders]&, InitalPopulation, NumberOfGenerations];
GenerationsFitness = Map[Max, Generations, {2}];
BestColoring = Extract[Generations, Position[GenerationsFitness, BestFitness, {2}, 1]][];
If[Or @@ (BestColoring[[First[#]]] == BestColoring[[Last[#]]]& /@ First /@ Most[ArrayRules[AdjacencyMatrix[g]]]),
$Failed, {BestFitness, BestColoring} ] ]  ## For those with Mathematica 7 or Less Here is code that doesn't use the version 8 Graph object, it's pretty much exactly the same: GAColor[adjmatrix_, PopulationSize_Integer:100, NumberOfGenerations_Integer:10, NumberOfMutants_Integer:0, mutationRadius_Integer:Automatic] := Module[ {NumberOfVertices = Length @ adjmatrix, NumberOfBreeders, PermuteColorClasses, MutationRadius, Combine, Mutate, PopulationStep, InitializePopulation, InitalPopulation, Generations, BestFitness = \[Infinity], BestColoring, GenerationsFitness, Chromatize, Adjacencies, result}, MutationRadius = If[mutationRadius === Automatic, NumberOfVertices, mutationRadius]; Adjacencies = Last /@ Transpose /@ GatherBy[First /@ Most[ArrayRules[adjmatrix]], First]; NumberOfBreeders = PopulationSize - NumberOfMutants; PermuteColorClasses[colors_, n_:1] := Module[{p}, p = {#, Flatten @ Position[colors, #]}& /@ Union[colors]; p[[All, 1]] = RandomSample[p[[All,1]]]; ReplacePart[ConstantArray[0, Length[colors]], Flatten[Rule @@@ Thread[Reverse@#]& /@ p]] ]; Chromatize[colorVector_] := Module[{f, h, min, co = colorVector}, f = Function[{c,v}, ReplacePart[c, v -> With[{ncols = c[[Adjacencies[[v]]]]}, For[min = 1, MemberQ[ncols, min], min++]; min] ]]; h = Function[{c}, Fold[f, c, RandomSample[Range[NumberOfVertices]]]]; FixedPoint[h, co] ]; Combine[colorVector1_, colorVector2_] := MapThread[RandomChoice[{#1, #2}]&, {colorVector1, colorVector2}]; Mutate[colorVector_, mr_] := Permute[colorVector, RandomPermutation[mr]]; PopulationStep[population_, NumberOfBreeders_] := Module[ {fitness = Max /@ population, breeders, children, mutants}, With[{min = Min[fitness]}, If[min < BestFitness, BestFitness = min]]; breeders = RandomChoice[fitness -> population, NumberOfBreeders]; children = Chromatize /@ Table[Combine @@ RandomChoice[breeders, 2], {NumberOfBreeders} ]; mutants = Mutate[#, MutationRadius]& /@ RandomChoice[breeders, NumberOfMutants]; Join[children, mutants] ]; InitalPopulation = With[{color = Chromatize[RandomSample[Range @ NumberOfVertices]]}, Table[PermuteColorClasses[color], {PopulationSize}] ]; Generations = NestList[PopulationStep[#, NumberOfBreeders]&, InitalPopulation, NumberOfGenerations]; GenerationsFitness = Map[Max, Generations, {2}]; BestColoring = Extract[Generations, Position[GenerationsFitness, BestFitness, {2}, 1]][]; If[Or @@ (BestColoring[[First[#]]] == BestColoring[[Last[#]]]& /@ First /@ Most[ArrayRules[adjmatrix]]),$Failed, {BestFitness, BestColoring}
]
]


Here is the sample input graph (as a compressed adjacency matrix) to test it on: http://pastebin.com/t7gnTczD

The algorithm should give a chromatic number of 28 in a few seconds. Here are the other benchmarks: http://www.info.univ-angers.fr/pub/porumbel/graphs/

Even in version 8 of Mathematica there still are no tools to compute the chromatic number or index of a graph, let alone a fast upper bound. Here is an illustration of the simulated annealing that's going on inside the algorithm:

g = Uncompress@"1:eJzt...."; (* get this string from pastebin link *)
NumberOfVertices = Length @ g;
color = RandomSample[Range[NumberOfVertices], NumberOfVertices];
n = NumberOfVertices;
A = Last /@ Transpose /@ GatherBy[First /@ Most[ArrayRules[g]], First];

NeighborComplements = Function[c,
Module[{p, n, nc, r},
p = Flatten @ Position[color, c];
n = A[[p]];
nc = Map[color[[#]]&, n, {2}];
Thread @ {p, Complement[Range[c], #, {c}]& /@ nc}
]
];

Chromatic[g_, n_, col_:Range[NumberOfVertices]] := Module[{c=col, f, h, slow, fast},
f = Function[{c, v},
ReplacePart[c,
v -> Module[{i, com = Complement[Range[n], c[[A[[v]]]]]},
RandomChoice[Join[com, {c[[v]]}]]
]
]
];
h = Function[{c}, Fold[f, c, RandomSample[Range[NumberOfVertices], NumberOfVertices]]];
NestWhile[h, c, Max[#]>n&]
];

AbsoluteTiming[Monitor[color = NestWhile[Chromatic[g, n-=1, #]&, color, (color=#;n>1)&],
ListPlot[Sort @ color, PlotRange -> All, PlotLabel -> Max[color]]]]


This is an optimization problem, and I'm sure some of you know this area intimately. When you run this code you will see a plot of the color classes which decrease slowly to around 30 different colors for the 250 vertices, however this is only a local minimum, the global minimum and chromatic number of the graph is actually 28... so my code is inefficient, if you can design a completely new function, and/or use openCL or JavaLink that is ok too...

• Is it necessary to use a GA? Oct 22, 2012 at 23:34
• I don't know that much about GAs but my impression of them is that it's trivial to write something that works after a fashion, but to do the job well requires either a lot of domain-specific knowledge and experience or very many experiments. So, this question might end up being more about GAs than about Mathematica per se. (I don't consider that a bad thing as I like domain-specific questions and find evolutionary computation techniques interesting. Others may or may not agree.) Here's a small hint: Darwinian evolution is very inefficient! Oct 23, 2012 at 2:07
• @DanielLichtblau I'm after a fast upperbound approx algorithm to push stuff, and I'm keen on getting a GA/GP to work. I've tried others, but if a stochastic, deterministic, and annealing method works better who am I to say, just test against the benchmark...
– M.R.
Oct 23, 2012 at 4:11
• Would it be too much trouble for you to give some code and an example of use that does not use Graph? Since I don't have that functionality I'm not sure what your data should look like. Oct 24, 2012 at 17:27
• Thanks for the update. Reading now. Oct 24, 2012 at 18:01

A lot of time is being spent on:

For[min = 1, MemberQ[ncols, min], min++]; min


Replacing it with:

First @ Complement[Range @ NumberOfVertices, ncols]


appears to speed things considerably.

This appears to be about twice as fast as the original. Can you confirm?

GAColor2[adjmatrix_, PopulationSize_Integer : 100,
NumberOfGenerations_Integer : 10, NumberOfMutants_Integer : 0,
PopulationStep, InitializePopulation, InitalPopulation,
Generations, BestFitness = \[Infinity], BestColoring,

Vrange = Range@NumberOfVertices;

Last /@ Transpose /@

NumberOfBreeders = PopulationSize - NumberOfMutants;

PermuteColorClasses[colors_, n_ : 1] :=
Module[{p},
p = ({#1, Flatten[Position[colors, #1]]} &) /@ Union[colors];
p[[All, 1]] = RandomSample[p[[All, 1]]];
ReplacePart[ConstantArray[0, Length[colors]],
Flatten[(Apply[Rule, Thread[Reverse[#1]], {1}] &) /@ p]]];

Chromatize[colorVector_] :=
Module[{f, h},
f = Function[{c, v},
ReplacePart[c,
h = Function[{c}, Fold[f, c, RandomSample[Vrange]]];
FixedPoint[h, colorVector]
];

Combine[colorVector1_, colorVector2_] :=

Mutate[a_, n_Integer] := RandomSample@Take[a, n];

PopulationStep[population_, NumberOfBreeders_] :=
Module[{fitness = Max /@ population, breeders, children, mutants},
With[{min = Min[fitness]},
If[min < BestFitness, BestFitness = min]];
breeders = RandomChoice[fitness -> population, NumberOfBreeders];
children =
Chromatize /@
Table[Combine @@ RandomChoice[breeders, 2], {NumberOfBreeders}];
mutants = (Mutate[#1, MutationRadius] &) /@
RandomChoice[breeders, NumberOfMutants];
Join[children, mutants]];

InitalPopulation =
With[{color = Chromatize[RandomSample[Vrange]]},
Table[PermuteColorClasses[color], {PopulationSize}]];

Generations =
NestList[PopulationStep[#1, NumberOfBreeders] &, InitalPopulation,
NumberOfGenerations];

GenerationsFitness = Map[Max, Generations, {2}];

BestColoring =
Extract[Generations,
Position[GenerationsFitness, BestFitness, {2}, 1]][];

If[Or @@ (BestColoring[[First[#1]]] == BestColoring[[Last[#1]]] &) /@
BestColoring}]
]


You highlighted Combine and Mutate for attention, and though I doubt these are the slow part of your code I shall address them.

These should be a bit faster:

combine[a_, b_] := a # + b (1 - #) & @ RandomInteger[1, Length@a]

mutate[a_, n_Integer] := RandomSample @ Take[a, n]


I am assuming that the second argument of mutate is an integer, not a list.

n = 2000000;
a = RandomInteger[99, n];
b = RandomReal[1, n];

Combine[a, b] // Timing // First

combine[a, b] // Timing // First


1.342

0.593

Mutate[a, 2^#] & ~Array~ 20 // Timing // First

mutate[a, 2^#] & ~Array~ 20 // Timing // First


0.2808

0.04368

• Thanks, see my comments to run the code in version 7.
– M.R.
Oct 24, 2012 at 18:01
• Version 8 has no GraphColor[] functionality... are you able to import the sample inputs from the website? info.univ-angers.fr/pub/porumbel/graphs I believe that you would need v8 Import to get the .col files... I can post them as compressed matrices if you need.
– M.R.
Oct 24, 2012 at 18:02
• @Mr.Wizard Did you profile this? A blog post on profiling would be very welcome... :) Oct 24, 2012 at 20:19
• Yes, this is faster but gets stuck at 40 colors. The annealing gets down to 31...
– M.R.
Oct 24, 2012 at 21:10
• @Mr.Wizard Combining the annealing with some sort of branch and bound partial backtrack could do the trick... The real problem is getting a 28 coloring in a reasonable time, whatever the method, this is the desired result. I'm thinking I might have to write this in c...
– M.R.
Oct 24, 2012 at 21:35

[Not an answer, but too big for a comment]

An alternative approach is to write it as an integer linear program. Below is code for this, with appropriate post-processing omitted to confuse the weak-minded (starting with the author).

n = Length[VertexList[g]];
edges = EdgeList[g];
vars = Array[x, {n, n}];
fvars = Flatten[vars];
tvars = Transpose[vars];
c2 = Map[0 <= # <= 1 &, fvars];
c3 = Table[Map[x[#[], j] + x[#[], j] <= 1 &, edges], {j, n}];
colvars = Array[y, n];
c4 = Map[0 <= # <= 1 &, colvars];
c5 = Table[n*colvars[[j]] >= Total[tvars[[j]]], {j, n}];
obj = Total[colvars];
allvars = Join[fvars, colvars];
constraints =
Flatten[Join[c1, c2, c3, c4, c5, {Element[allvars, Integers]}]];

Timing[min = FindMinimum[{obj, constraints}, allvars]]


This will not run in finite time for the example in question. There may be variations that do better though.

--- edit ---

One possibility for a heuristic approach based on this setup is to change the FindMinimum call to e.g.

NMinimize[{obj, constraints}, allvars,
Method -> {"DifferentialEvolution", "CrossProbability" -> .1,
"SearchPoints" -> 100}, MaxIterations -> 500]


This is not alone sufficient though. I think you will need to play with giving an "InitialPoints" option as well. Have not had time to try that.

--- end edit ---

• Nice idea, but I feel this approach is doomed as this problem is NP-complete. Isn't the point of approximation algorithms to use clever heuristics and shortcuts to get fast OK solutions? The speed is the issue and any sort of backtracking or exactly solving won't ever work on super large graphs.
– M.R.
Oct 24, 2012 at 18:36
• But perhaps you can think of a way to relax these conditions or randomize things to make it faster with some probability of stopping?
– M.R.
Oct 24, 2012 at 21:36
• See edit for an incomplete possible heuristic direction. Oct 25, 2012 at 0:01
• Hmm, it still doesn't terminate...
– M.R.
Oct 25, 2012 at 16:53
• Is FindMinimum guaranteed to find a global minimum in situations like this? Feb 14, 2019 at 15:21