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

Question

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[293]= {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},

    MutationRadius = If[mutationRadius === Automatic, NumberOfVertices, mutationRadius];
    Adjacencies = Last /@ Transpose /@ GatherBy[First /@ Most[ArrayRules[AdjacencyMatrix[g]]], 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]][[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]][[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...

Answer 1

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, 
  mutationRadius_Integer : Automatic] := 
 Module[{NumberOfVertices = Length[adjmatrix], NumberOfBreeders, 
   PermuteColorClasses, MutationRadius, Combine, Mutate, 
   PopulationStep, InitializePopulation, InitalPopulation, 
   Generations, BestFitness = \[Infinity], BestColoring, 
   GenerationsFitness, Chromatize, Adjacencies, result, Vrange},

  Vrange = Range@NumberOfVertices;

  MutationRadius = 
   If[mutationRadius === Automatic, NumberOfVertices, mutationRadius];

  Adjacencies = 
   Last /@ Transpose /@ 
     GatherBy[ArrayRules[adjmatrix][[;; -2, 1]], First];

  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, 
       v -> First@Complement[Vrange, c[[Adjacencies[[v]]]]]]]; 
    h = Function[{c}, Fold[f, c, RandomSample[Vrange]]];
    FixedPoint[h, colorVector]
    ];

  Combine[colorVector1_, colorVector2_] := 
   MapThread[RandomChoice[{#1, #2}] &, {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]][[1]];

  If[Or @@ (BestColoring[[First[#1]]] == BestColoring[[Last[#1]]] &) /@
      First /@ Most[ArrayRules[adjmatrix]], $Failed, {BestFitness, 
    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

Answer 2

[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];
c1 = Thread[Total[tvars] == 1];
c2 = Map[0 <= # <= 1 &, fvars];
c3 = Table[Map[x[#[[1]], j] + x[#[[2]], 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 ---