Graph Coloring

Graph Coloring

You are an expert in vertex coloring: exact MIP and CP models, construction heuristics (DSATUR, RLF), tabu-search-based k-coloring (tabucol), Kempe-chain moves, and clique-based lower bounds. This skill covers the minimum-coloring problem, the k-coloring decision problem, and the application patterns that reduce to them (register allocation, frequency assignment, exam timetabling). Use the framework below to pick the right model and method for the instance size at hand, and always pair an upper bound (a coloring) with a lower bound (a clique or LP bound) so you can state the optimality gap.

Initial Assessment

Establish these facts before proposing a model or algorithm:

• Objective type. Minimum number of colors (chromatic number), or a fixed color budget k where you only need feasibility (k-coloring decision)? Frequency-style problems often fix k and minimize interference instead.
• Graph size and density. Vertices n, edges m, density 2m / (n(n-1)). Exact methods are realistic up to roughly n = 80-100 on dense random graphs; sparse structured graphs can be far larger. Heuristics handle millions of vertices.
• Graph structure. Random, geometric, interval, planar, or derived from an application (interference graph, conflict graph)? Interval graphs and chordal graphs are colorable optimally in polynomial time — check before deploying heavy machinery.
• Hard vs soft constraints. Pure coloring has only hard "endpoints differ" constraints. If the user mentions preferences, spread requirements, or penalties (exams close together, adjacent channels), the problem is a coloring-flavored timetabling/assignment problem; the coloring core still applies but the objective changes.
• Precoloring or list constraints. Are some vertices fixed to specific colors (precoloring extension)? Does each vertex have its own allowed color list (list coloring)? Both are easy to add to MIP/CP and to construction heuristics, but they invalidate symmetry-breaking tricks based on color interchangeability.
• Solver availability. Gurobi license available? If not, OR-Tools CP-SAT is free and is usually the stronger exact tool for coloring anyway; HiGHS/CBC handle the MIP variant.
• Quality requirement. Proof of optimality required (publication, exact benchmark), or is a good coloring with a reported gap acceptable (engineering use)?
• Time budget. Seconds (greedy/DSATUR only), minutes (tabucol descent + clique bound), hours (exact attempt with CP-SAT or branch-and-price)?
• Data format. Adjacency matrix, edge list, DIMACS .col file? Standard benchmarks (DIMACS challenge graphs: DSJC, flat, le450 families) come as DIMACS edge lists.
• Reproducibility. Fix seeds for instance generation and for every randomized heuristic; report them.

Problem Definition and Model Landscape

Formal definition

Given an undirected graph $G = (V, E)$ with $|V| = n$ and $|E| = m$, a k-coloring is a map $c : V \to \{1, \dots, k\}$. It is proper if $c(u) \neq c(v)$ for every edge $\{u, v\} \in E$. The chromatic number $\chi(G)$ is the smallest $k$ admitting a proper k-coloring. Each color class $c^{-1}(i)$ is an independent set, so coloring is equivalent to partitioning $V$ into the fewest independent sets.

Key bounds, with $\omega(G)$ the clique number, $\Delta(G)$ the maximum degree, and $d(G)$ the degeneracy (largest minimum degree over all subgraphs):

$$\omega(G) \;\le\; \chi_f(G) \;\le\; \chi(G) \;\le\; d(G) + 1 \;\le\; \Delta(G) + 1$$

where $\chi_f$ is the fractional chromatic number (LP bound of the set-covering formulation). Brooks (1941) sharpens the upper bound: $\chi(G) \le \Delta(G)$ unless $G$ is complete or an odd cycle. Greedy coloring in a smallest-last (degeneracy) order achieves $d(G)+1$ colors (Matula & Beck 1983).

Complexity: deciding k-colorability is NP-complete for every fixed $k \ge 3$; $\chi(G)$ is NP-hard to approximate within $n^{1-\varepsilon}$ (Zuckerman 2007). Exceptions worth checking: bipartite graphs ($\chi = 2$, test via BFS), interval and chordal graphs (perfect elimination ordering gives $\chi = \omega$ in linear time), planar graphs ($\chi \le 4$).

Assignment MIP model

With a color budget $H$ (any valid upper bound, e.g. the DSATUR value), binary $x{vc} = 1$ iff vertex $v$ takes color $c$, and $yc = 1$ iff color $c$ is used:

$$ \begin{aligned} \min \; & \sum{c=1}^{H} yc \\ \text{s.t.} \; & \sum{c=1}^{H} x{vc} = 1 \qquad && \forall v \in V \\ & x{uc} + x{vc} \le yc && \forall \{u,v\} \in E,\; c = 1,\dots,H \\ & yc \le y{c-1} && c = 2,\dots,H \\ & x{vc},\, y_c \in \{0,1\} \end{aligned} $$

Two structural weaknesses to know in advance. First, the LP relaxation is nearly useless: $x{vc} = 1/H$, $yc = 2/H$ is feasible whenever $E \neq \emptyset$, so the LP bound is at most 2 regardless of $\chi(G)$. Second, colors are interchangeable, so the model has $H!$ symmetric copies of every solution; the ordering constraints $yc \le y{c-1}$ remove only part of this. Always pre-fix a clique (vertex $i$ of a clique $Q$ gets color $i$) — it breaks symmetry much harder and injects the bound $\chi \ge |Q|$. See integer-programming-techniques for the general theory of symmetry and relaxation strength.

Stronger exact alternatives:

• Representatives formulation (Campêlo, Campos & Corrêa 2008): binary $x_{uv}$ = "u represents the color class of v", defined on non-adjacent pairs. No color indices, hence no color symmetry; tighter LP bound.
• Set covering over independent sets (Mehrotra & Trick 1996): $\min \sumS \lambdaS$ over independent sets $S$ covering every vertex. Its LP value is $\chi_f(G)$ — the strongest practical lower bound — solved by column generation with a maximum-weight independent set pricing problem; branch-and-price closes many DIMACS instances.
• CP model: one integer variable per vertex with domain $\{0,\dots,H-1\}$, a binary disequality per edge, minimize the maximum color. Posting AllDifferent on extracted cliques strengthens propagation substantially. CP-SAT with clique fixing is usually the best off-the-shelf exact method for coloring.

Method selection

| Situation | Method | |---|---| | Need proof of optimality, n up to ~80-100 (dense) | CP-SAT with clique fixing; MIP as fallback | | Need proof, larger sparse or structured graph | Branch-and-price over independent sets (Mehrotra & Trick 1996) | | Fixed k, want a proper coloring fast | tabucol; PARTIALCOL for hard instances | | Quick upper bound, any size | DSATUR; RLF when density > ~0.3 | | Huge sparse graph (millions of vertices) | Greedy in degeneracy order, no matrix storage | | Lower bound | Multi-start greedy clique; exact max clique on small graphs; $\chi_f$ via column generation | | Bipartite / interval / chordal suspected | Test the structure first; polynomial algorithms apply |

Applications map

| Application | Graph construction | Notes | |---|---|---| | Register allocation | Interference graph: variables live at the same time are adjacent | k = number of registers; spill code when k-coloring fails (Chaitin 1982) | | Exam timetabling | Conflict graph: exams sharing a student are adjacent | Colors = time slots; soft constraints push beyond pure coloring — see timetabling-and-rostering (de Werra 1985) | | Frequency assignment | Interference graph between transmitters | Distance/bandwidth variants: $|c(u) - c(v)| \ge d{uv}$ (Aardal et al. 2007, frequency assignment survey) | | Sports league scheduling | Edge coloring of the match graph | Edge coloring of $Kn$ = round-robin schedule | | Sudoku, Latin squares | Precoloring extension on a structured graph | CP handles these naturally |

Instance Generation and Validation

Group C contract: every coloring study needs a seeded generator and an independent validator. The validator recomputes feasibility and the objective from the adjacency matrix alone — never trust the solver's own bookkeeping.

import numpy as np

def random_graph(n: int, p: float, seed: int) -> np.ndarray:
    """Erdos-Renyi G(n, p) as a symmetric boolean adjacency matrix (no self-loops)."""
    rng = np.random.default_rng(seed)
    upper = np.triu(rng.random((n, n))  tuple[np.ndarray, np.ndarray]:
    """Random graph that is k-colorable by construction.

    Vertices are split into k classes uniformly at random; edges appear only
    between different classes, each with probability p. Returns (adjacency,
    planted coloring). Useful for testing: chi(G)  list[tuple[int, int]]:
    """Edge list with u  dict[str, int | bool]:
    """Independent feasibility + objective check for a vertex coloring.

    Returns conflict count (violated edges), the number of distinct colors used,
    and whether the color indices are compact (0..K-1 with no gaps). Recomputes
    everything from the adjacency matrix; shares no code with any solver.
    """
    n = adj.shape[0]
    if colors.shape != (n,):
        raise ValueError(f"colors has shape {colors.shape}, expected ({n},)")
    iu, jv = np.nonzero(np.triu(adj, k=1))
    conflicts = int(np.sum(colors[iu] == colors[jv]))
    used = np.unique(colors)
    return {
        "feasible": conflicts == 0,
        "conflicts": conflicts,
        "num_colors": int(used.size),
        "compact_indices": bool(np.array_equal(used, np.arange(used.size))),
    }

if __name__ == "__main__":
    c5 = np.zeros((5, 5), dtype=bool)
    for u, v in [(0, 1), (1, 2), (2, 3), (3, 4), (4, 0)]:
        c5[u, v] = c5[v, u] = True
    print(validate_coloring(c5, np.array([0, 1, 0, 1, 2])))
    # Expected: feasible=True, conflicts=0, num_colors=3, compact_indices=True
    print(validate_coloring(c5, np.array([0, 1, 0, 1, 0])))
    # Expected: feasible=False, conflicts=1 (edge (4,0) has both endpoints colored 0)

Exact Models: MIP and CP-SAT

The assignment MIP below follows the Group C constraint-builder pattern: each constraint family lives in its own named function, so families can be unit-tested and reused independently. Pass a clique from the lower-bound routine in Advanced Techniques to activate the clique-fixing family — on symmetric coloring models this is the single highest-impact modeling decision.

import gurobipy as gp
from gurobipy import GRB
import numpy as np

def add_assignment_constraints(model: gp.Model, x: gp.tupledict, data: dict) -> None:
    """Each vertex receives exactly one color."""
    for v in range(data["n"]):
        model.addConstr(
            gp.quicksum(x[v, c] for c in range(data["H"])) == 1, name=f"assign[{v}]"
        )

def add_conflict_constraints(
    model: gp.Model, x: gp.tupledict, y: gp.tupledict, data: dict
) -> None:
    """Edge endpoints cannot share color c; also links x to the color-use variable y."""
    for u, v in data["edges"]:
        for c in range(data["H"]):
            model.addConstr(x[u, c] + x[v, c]  None:
    """Color c may be used only if color c-1 is used (partial symmetry breaking)."""
    for c in range(1, data["H"]):
        model.addConstr(y[c]  None:
    """Pre-color a known clique: the i-th clique vertex takes color i.

    Valid because clique vertices need pairwise distinct colors and colors are
    interchangeable. Kills most color symmetry and implies the bound chi >= |Q|.
    """
    for i, v in enumerate(data["clique"]):
        model.addConstr(x[v, i] == 1, name=f"clique_fix[{v}]")

def solve_coloring_mip(
    adj: np.ndarray,
    H: int,
    clique: tuple[int, ...] = (),
    time_limit: float = 60.0,
) -> tuple[int, np.ndarray] | None:
    """Assignment-model MIP for minimum vertex coloring with color budget H.

    Returns (number of colors, coloring array) or None if no solution was found.
    H should be a known upper bound, e.g. the DSATUR color count.
    """
    iu, jv = np.nonzero(np.triu(adj, k=1))
    data = {
        "n": adj.shape[0],
        "H": H,
        "edges": list(zip(iu.tolist(), jv.tolist())),
        "clique": clique,
    }
    model = gp.Model("vertex_coloring")
    model.Params.OutputFlag = 0
    model.Params.TimeLimit = time_limit
    model.Params.Symmetry = 2  # aggressive solver-side symmetry detection on top
    x = model.addVars(data["n"], H, vtype=GRB.BINARY, name="x")
    y = model.addVars(H, vtype=GRB.BINARY, name="y")
    add_assignment_constraints(model, x, data)
    add_conflict_constraints(model, x, y, data)
    add_color_ordering_constraints(model, y, data)
    add_clique_fixing_constraints(model, x, data)
    model.setObjective(gp.quicksum(y[c] for c in range(H)), GRB.MINIMIZE)
    model.optimize()
    ok = model.Status == GRB.OPTIMAL or (
        model.Status == GRB.TIME_LIMIT and model.SolCount > 0
    )
    if not ok:
        return None
    colors = np.array(
        [max(range(H), key=lambda c: x[v, c].X) for v in range(data["n"])], dtype=int
    )
    return int(round(model.ObjVal)), colors

if __name__ == "__main__":
    c5 = np.zeros((5, 5), dtype=bool)
    for u, v in [(0, 1), (1, 2), (2, 3), (3, 4), (4, 0)]:
        c5[u, v] = c5[v, u] = True
    result = solve_coloring_mip(c5, H=5, clique=(0, 1))
    print(result)
    # Expected: (3, coloring) — an odd cycle has chromatic number 3

CP-SAT is usually the stronger exact tool here: the disequality per edge propagates well, there are no big-M or linking artifacts, and clique fixing plugs in the same way. Minimize the maximum color index instead of counting used colors — with compact color use these coincide, and the max-objective propagates better.

import numpy as np
from ortools.sat.python import cp_model

def solve_coloring_cpsat(
    adj: np.ndarray,
    H: int,
    clique: tuple[int, ...] = (),
    time_limit: float = 30.0,
) -> tuple[int, np.ndarray] | None:
    """CP-SAT minimum coloring: integer color per vertex, disequality per edge.

    Returns (number of colors, coloring) or None. Fixing a clique to colors
    0..|Q|-1 breaks color symmetry and seeds the lower bound. For interval
    variables, channeling, and search strategies see the constraint-programming
    skill.
    """
    n = adj.shape[0]
    model = cp_model.CpModel()
    x = [model.NewIntVar(0, H - 1, f"x[{v}]") for v in range(n)]
    max_color = model.NewIntVar(0, H - 1, "max_color")
    iu, jv = np.nonzero(np.triu(adj, k=1))
    for u, v in zip(iu.tolist(), jv.tolist()):
        model.Add(x[u] != x[v])
    model.AddMaxEquality(max_color, x)
    for i, v in enumerate(clique):
        model.Add(x[v] == i)
    model.Minimize(max_color)
    solver = cp_model.CpSolver()
    solver.parameters.max_time_in_seconds = time_limit
    solver.parameters.num_workers = 8
    status = solver.Solve(model)
    if status not in (cp_model.OPTIMAL, cp_model.FEASIBLE):
        return None
    colors = np.array([solver.Value(x[v]) for v in range(n)], dtype=int)
    return int(solver.Value(max_color)) + 1, colors

if __name__ == "__main__":
    c5 = np.zeros((5, 5), dtype=bool)
    for u, v in [(0, 1), (1, 2), (2, 3), (3, 4), (4, 0)]:
        c5[u, v] = c5[v, u] = True
    print(solve_coloring_cpsat(c5, H=5, clique=(0, 1)))
    # Expected: (3, coloring) — matches the MIP result, typically in milliseconds

When you need the decision version (is the graph k-colorable for fixed k?), drop max_color and the objective, set the domains to 0..k-1, and call the solver once per k; CP-SAT's clause learning shines on these satisfiability-style runs.

Construction Heuristics: DSATUR and RLF

DSATUR (Brélaz 1979) colors the vertex with the highest saturation degree (number of distinct colors among its neighbors) first, breaking ties by degree. It is exact on bipartite graphs and is the default quick upper bound — also the standard source of the color budget H for the exact models and the starting point for tabucol. RLF (Leighton 1979) builds one maximal independent set (color class) at a time; it uses more time per color but typically fewer colors on dense graphs.

import numpy as np

def greedy_coloring(adj: np.ndarray, o

…

## Source & license

This open-source skill is cataloged on AgentStack and links to its original source — we do not rehost the code.

- **Author:** [hajibabaie](https://github.com/hajibabaie)
- **Source:** [hajibabaie/combinatorial-optimization-skills](https://github.com/hajibabaie/combinatorial-optimization-skills)
- **License:** MIT

Install and usage instructions live in the source repository linked above.