[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 ---