Chromatic Polynomial, in depth: why colourings are counted by a polynomial, deletion-contraction by hand, two tested exact algorithms and their limits

Definition and why it is a polynomial

Let G have n vertices and m edges, with no loops. A proper k-colouring assigns each vertex one of k colours so that the two ends of every edge differ. Colours are labelled, so swapping red and blue gives a different colouring.

Some graphs can be counted directly. With no edges every vertex is free: k^n. For the complete graph K_n the first vertex has k choices, the next k - 1, and so on: k(k-1)(k-2)...(k-n+1). For a tree, colour a root (k choices) and then each further vertex only has to differ from its parent: k(k-1)^(n-1). The cycle C_n gives (k-1)^n + (-1)^n (k-1).

Why is the count always a polynomial? Group colourings by the partition of vertices into colour classes. Each class must be an independent set. If a_j is the number of ways to split the vertices into exactly j non-empty independent sets, then assigning distinct colours to the j classes can be done in k(k-1)...(k-j+1) ways, so

P(G, k) = sum over j of  a_j * k(k-1)(k-2)...(k-j+1)

Each falling factorial is a polynomial in k, so P is too. This form also hands us an algorithm.

Deletion-contraction

The classical tool is deletion-contraction. Take any edge e = uv. Colourings of G - e (the graph with e removed) fall into two groups: those where u and v differ, which are exactly the proper colourings of G, and those where u and v share a colour, which correspond one-to-one with colourings of G / e, the graph with u and v merged into a single vertex. So

P(G, k) = P(G - e, k) - P(G / e, k)

Both graphs on the right have one edge fewer, so the recursion ends at edgeless graphs, where P = k^n. When contracting, parallel edges merge into one; a loop would mean the graph has no proper colouring, but it cannot arise when contracting a simple graph along an existing edge.

Worked example: the diamond graph

Work through the diamond: a 4-cycle 0-1-2-3 with a chord 0-2.

1. C4 first. Pick edge 3-0. Deleting it leaves the path 0-1-2-3, a tree: k(k-1)^3. Contracting it merges 3 into 0 and gives a triangle: k(k-1)(k-2). So P(C4) = k(k-1)^3 - k(k-1)(k-2) = k(k-1)(k^2 - 3k + 3) = k^4 - 4k^3 + 6k^2 - 3k. The cycle formula agrees: (k-1)^4 + (k-1) expands to the same thing.
2. Diamond. Pick the chord 0-2. Deleting it leaves C4. Contracting it merges 0 and 2; the edges 0-1 and 2-1 become one edge, as do 0-3 and 2-3, leaving a path on three vertices: k(k-1)^2.
3. Subtract. P(diamond) = k(k-1)(k^2 - 3k + 3) - k(k-1)^2 = k(k-1)(k^2 - 4k + 4) = k(k-1)(k-2)^2 = k^4 - 5k^3 + 8k^2 - 4k.
4. Evaluate. P(2) = 0, so the diamond is not 2-colourable (it has triangles). P(3) = 3 x 2 x 1 x 1 = 6: choose colours for the triangle 0-1-2 in 6 ways, and vertex 3 must take the colour of 1, the only one that differs from both 0 and 2. P(4) = 48.

The factorised form has a reason. The diamond is chordal, and for chordal graphs P(G, k) is the product of (k - d_i), where d_i counts each vertex's already-placed neighbours along a perfect elimination ordering: here 0, 1, 2 and 2 earlier neighbours. The chordal graphs article shows how to find such an order in linear time, which makes the polynomial of any chordal graph cheap to compute.

Coefficients and other properties

Whitney showed in 1932 that the coefficients count edge subsets: P(G, k) is the sum over all edge subsets A of (-1)^|A| k^c(A), where c(A) is the number of connected components of the spanning subgraph with edges A. Several facts follow, and each is a free test for your code:

• Degree n, leading coefficient 1, constant term 0 (for n > 0).
• The coefficient of k^(n-1) is -m, minus the number of edges.
• The coefficient of k^(n-2) is C(m, 2) minus the number of triangles. C4: 6. Diamond: 10 - 2 = 8. K4: 15 - 4 = 11.
• Coefficients alternate in sign, and the lowest non-zero power of k equals the number of connected components.
• Disjoint union multiplies: P(G union H) = P(G) P(H). Gluing two graphs along a shared complete subgraph K_r gives P(G) P(H) / P(K_r).
• Stanley (1973): (-1)^n P(G, -1) counts the acyclic orientations of G. For C4 this is 14, which is 2^4 orientations minus the 2 directed cycles.
• Real roots: none are negative, none lie in (0, 1), and Jackson proved in 1993 that none lie in (1, 32/27]. Integer roots 0, 1, ..., chi - 1 are always present, where chi is the chromatic number.

These properties also explain the link to physics. P(G, q) is the zero-temperature partition function of the antiferromagnetic q-state Potts model, and the chromatic polynomial is a specialisation of the Tutte polynomial, which unifies it with spanning-tree counts like those in the matrix-tree theorem article.

Two exact algorithms in code

Two exact algorithms are worth having. Polynomials are stored as integer coefficient lists, lowest degree first. The first is deletion-contraction with memoisation on the remaining vertex and edge sets:

from functools import lru_cache

def padd(p, q, sign=1):
    n = max(len(p), len(q))
    p = p + [0] * (n - len(p)); q = q + [0] * (n - len(q))
    out = [a + sign * b for a, b in zip(p, q)]
    while len(out) > 1 and out[-1] == 0:
        out.pop()
    return out

def chromatic_dc(n, edges):
    @lru_cache(maxsize=None)
    def go(verts, es):
        if not es:
            return tuple([0] * len(verts) + [1])           # k^|V|
        e = min(es, key=lambda f: tuple(sorted(f)))
        u, v = sorted(e)
        deleted = es - {e}
        merged = set()
        for f in deleted:                                    # contract v into u
            a, b = tuple(f)
            merged.add(frozenset((u if a == v else a, u if b == v else b)))
        a = list(go(verts, frozenset(deleted)))
        b = list(go(verts - {v}, frozenset(merged)))
        return tuple(padd(a, b, -1))
    return list(go(frozenset(range(n)), frozenset(frozenset(e) for e in edges)))

Without memoisation the call tree satisfies roughly T(n, m) = T(n, m-1) + T(n-1, m-1), which grows like Fibonacci numbers in n + m. The second algorithm follows the partition formula. It counts a_j for every subset of vertices with a dynamic program over bitmasks, anchoring each block on the lowest vertex so each partition is counted once, in O(3^n n) time:

def chromatic_subset_dp(n, edges):
    adj = [0] * n
    for u, v in edges:
        adj[u] |= 1 << v; adj[v] |= 1 << u
    indep = [True] * (1 << n)
    for S in range(1, 1 << n):
        low = (S & -S).bit_length() - 1
        rest = S & (S - 1)
        indep[S] = indep[rest] and not (adj[low] & rest)
    a = [[0] * (n + 1) for _ in range(1 << n)]   # a[S][j]: S split into j blocks
    a[0][0] = 1
    for S in range(1, 1 << n):
        low = S & -S
        rest = S ^ low
        T = rest
        while True:                               # blocks that contain low
            B = T | low
            if indep[B]:
                for j in range(1, n + 1):
                    a[S][j] += a[S ^ B][j - 1]
            if T == 0:
                break
            T = (T - 1) & rest
    poly = [0]
    for j in range(1, n + 1):
        if a[-1][j]:
            ff = [1]
            for i in range(j):                    # falling factorial in k
                ff = [0] + ff
                for t in range(len(ff) - 1):
                    ff[t] -= i * ff[t + 1]
            poly = padd(poly, [a[-1][j] * x for x in ff])
    return poly

Both functions agree on the cycle, diamond, K4, C5 and a disconnected test graph, and both match brute-force enumeration of all k^n assignments for k up to 4. On the Petersen graph (10 vertices, 15 edges) the subset DP returns P(3) = 120, which brute force over all 3^10 assignments confirms. Always keep a brute-force counter beside code like this; it is three lines and catches nearly every bug on graphs of up to about ten vertices.

The subset DP is the same style as the Held-Karp dynamic program in the TSP article. Björklund, Husfeldt and Koivisto (2009) improved it with inclusion-exclusion and fast subset transforms to O*(2^n) time, meaning 2^n times a polynomial factor.

Complexity and what scales

No efficient general algorithm is expected. Deciding whether P(G, 3) > 0 is the NP-complete 3-colouring problem, and computing P(G, k) for any fixed integer k of at least 3 is #P-hard. The easy evaluations are at k = 0 and 1 (zero once there is an edge) and k = 2, which is 2^c for a bipartite graph with c components and 0 otherwise. Practical choices by graph size:

Two engineering details matter. First, split the graph before computing: handle each connected component and each block separately and multiply, using the gluing rule. Second, prefer edges whose removal disconnects or simplifies the graph when choosing e in deletion-contraction; a good choice can shrink the tree dramatically.

Failure modes

• Floating-point evaluation. Coefficients alternate in sign and grow quickly (the Petersen polynomial already has 4305 in it), so evaluating the expanded form in floating point cancels badly at large n. Use exact integers, or evaluate from the falling-factorial form, whose terms are all non-negative for integer k at least n.
• Multigraph input. Parallel edges do not change the polynomial, but a self-loop makes it zero. Normalise input before computing, and reject loops explicitly.
• Confusing the questions. The chromatic number needs only the smallest k with P(G, k) > 0; computing the full polynomial for that is wasteful. Use a colouring solver instead.
• Memory blow-up. The subset DP stores (n + 1) 2^n counts; at n = 25 that is over 800 million big integers. Check sizes before running.
• Unlabelled versus labelled colours. P counts labelled colourings. To count partitions into at most k colour classes, use the a_j values instead.

What to do next

1. Compute P(C5) by hand with deletion-contraction and confirm (k-1)^5 - (k-1).
2. Run both functions above on a few graphs from your own work and compare them against a brute-force counter for small k.
3. Add the coefficient checks (leading 1, -m, C(m, 2) minus triangles, alternating signs) as assertions in your test suite.
4. For a scheduling or assignment problem, compute P(G, k) / k^n to see what fraction of random assignments are conflict-free, and use it to judge whether random restarts can work.
5. If your graphs are chordal or tree-like, implement the elimination-order product before reaching for exponential algorithms.