Maximum Clique Problem, in depth: colouring bounds, bitset branch and bound, k-core pruning and when to stop

Maximal, maximum and the problems in disguise

Two words cause most of the confusion. A maximal clique cannot be extended by adding any one vertex. A maximum clique is the largest clique in the graph. Every maximum clique is maximal; the converse is false. In the example graph below, the triangle {5, 6, 7} is maximal, because no other vertex is adjacent to all three, yet the graph contains the 4-clique {1, 2, 3, 4}.

This matters because the best-known clique algorithm, Bron-Kerbosch, enumerates maximal cliques. You can find the maximum by enumerating all of them and keeping the biggest, and on sparse graphs with small degeneracy that works well. But a graph can have exponentially many maximal cliques (up to 3 to the power n/3, the Moon-Moser bound). A maximum-clique solver instead prunes every branch that provably cannot beat the best clique so far.

Two neighbouring problems are the same problem in disguise. A clique in G is an independent set in the complement graph, so maximum independent set is maximum clique on the complement. The complement of an independent set is a vertex cover, so minimum vertex cover is n minus the maximum independent set. A sparse graph has a dense complement, so pick the formulation your solver handles well. Vertex cover also has a 2-approximation, covered in vertex cover approximation, but that guarantee does not carry over to clique. Unless P = NP, clique cannot be approximated within a factor of n to the power 1 minus epsilon, for any epsilon above zero. In practice the useful choices are exact search and heuristics that come with no guarantee.

Branch and bound with a colouring bound

Branch and bound keeps three things: the current clique R, a candidate set P of vertices adjacent to everything in R, and the best clique found so far. Branching picks a vertex v from P and recurses on R plus v with candidates P intersected with the neighbours of v. After the branch returns, v is removed from P, so later branches never reconsider it. The bound asks one question: can R plus anything from P beat the best? If not, return.

The weakest bound is |R| + |P|, used by the early Carraghan-Pardalos style of solver. It is correct, because no clique can use more vertices than exist, but it is loose. The bound that changed practice, popularised by Tomita's MCQ family, is a greedy colouring of P. Vertices that share a colour are pairwise non-adjacent, so a clique can contain at most one vertex of each colour. If P can be coloured with k colours, then no clique inside P has more than k vertices, and |R| + k is a valid upper bound. Colouring is cheap and greedy colouring is not optimal, but the bound only has to be valid, not tight.

The colouring does double duty. Order P by colour class and branch from the highest colour downward. When you reach a vertex whose colour number c gives |R| + c no larger than the best size, every remaining vertex has colour c or less, so the whole rest of the loop can be cut with one comparison. Tomita's MCR and MCS refine the vertex order and recolouring, San Segundo's BBMC uses bitsets, and Ostergard's Cliquer instead reuses results for vertex suffixes as bounds.

A bitset implementation

The implementation stores each neighbourhood as a bitset. Python integers are arbitrary-precision, so 1 << v and & work for any n. In C or Rust you would use arrays of 64-bit words, which is what BBMC does. The candidate intersection P & adj[v] becomes a word-parallel AND, the main reason bitset solvers beat adjacency-list solvers on dense graphs by large factors.

def max_clique(n, edges):
    """Exact maximum clique. Vertices 0..n-1; Python ints serve as bitsets."""
    adj = [0] * n
    for u, v in edges:
        adj[u] |= 1 << v
        adj[v] |= 1 << u
    best = 0                                   # bitset of the best clique found

    def colour_order(P):
        """Greedy sequential colouring; vertices returned in non-decreasing colour."""
        order, colour = [], 0
        while P:
            colour += 1
            Q = P                              # vertices still eligible for this colour
            while Q:
                v = (Q & -Q).bit_length() - 1  # lowest set bit
                Q &= ~(1 << v) & ~adj[v]       # neighbours cannot share the colour
                P &= ~(1 << v)
                order.append((v, colour))
        return order

    def expand(R, size, P):
        nonlocal best
        for v, colour in reversed(colour_order(P)):
            if size + colour <= best.bit_count():
                return                         # every remaining vertex is bounded too
            NP = P & adj[v]
            if NP:
                expand(R | (1 << v), size + 1, NP)
            elif size + 1 > best.bit_count():
                best = R | (1 << v)            # R + v is maximal and a new record
            P &= ~(1 << v)

    expand(0, 0, (1 << n) - 1)
    return [v for v in range(n) if best >> v & 1]

Note int.bit_count(), which needs Python 3.10 or later. Two details are easy to get wrong. First, the prune is a return, not a continue: colours only decrease from there, so once one vertex is bounded, all the rest are too. Second, P must shrink after each branch, or the solver explores every clique once per vertex ordering. This code was cross-checked against brute force on 300 random graphs; do the same for any change, because a bound bug returns a valid but too-small clique rather than crashing.

Worked example

Take the eight-vertex graph in the figure, with edges 0-1, 0-2, 1-2, 1-3, 1-4, 2-3, 2-4, 3-4, 4-5, 5-6, 5-7, 6-7 and 0-7. At the root, P holds all eight vertices. Greedy colouring in index order gives colour 1 to {0, 3, 5}, colour 2 to {1, 6}, colour 3 to {2, 7} and colour 4 to {4}. The root bound is therefore 4, so the solver knows before it starts that no clique has more than four vertices.

The loop starts with the highest colour, vertex 4. Its candidates are {1, 2, 3, 5}. Colouring them gives 1 and 5 colour 1, 2 colour 2 and 3 colour 3, so the bound is 1 + 3 = 4, which beats the empty incumbent. Branching on 3 leaves candidates {1, 2}, then branching on 2 leaves {1}, and adding 1 empties the candidate set with R = {1, 2, 3, 4}. That becomes the incumbent, of size 4. Back at the root, the next vertex is 7 with colour 3, and 0 + 3 is not greater than 4, so the loop returns. The search is finished, and the triangle {5, 6, 7} was never explored, because the colouring proved it could not win.

Bron-Kerbosch on the same graph reports all five maximal cliques, {1, 2, 3, 4}, {0, 1, 2}, {5, 6, 7}, {4, 5} and {0, 7}. That is fine here and hopeless on a dense graph with millions of them.

Preprocessing that does most of the work

On large sparse graphs, preprocessing does more work than the search does. Three steps carry most of the benefit.

1. A heuristic incumbent first. Process vertices in degeneracy order, greedily grow a clique from each one's later neighbours, and keep the best. It takes near-linear time and is often optimal or one short on real graphs. Every vertex the next step removes depends on this number.
2. k-core pruning. A vertex in a clique of size k has degree at least k - 1 inside the clique, so its core number is at least k - 1. With an incumbent of size b, any vertex whose core number plus one is no larger than b cannot belong to a strictly larger clique. Delete it, update degrees and repeat. In the example, vertices 0, 5, 6 and 7 have core number 2 and the rest have 3. A triangle incumbent removes 0, 5, 6 and 7; the 4-clique removes everything, proving optimality without search. On social and web graphs this step often removes nearly all vertices.
3. Vertex ordering. Number vertices so that the initial colouring is tight, typically by degeneracy order or decreasing degree. MCQ-style solvers are sensitive to this, and a bad order can cost orders of magnitude.

After pruning, search each remaining vertex's neighbourhood separately: find the largest clique containing v among v's later neighbours. Each subproblem is small and dense, where bitsets work best, and independent, so they spread across cores sharing one incumbent. A stale incumbent only slows pruning; it never makes it wrong.

The ILP formulation and when to use it

The textbook integer program has a 0-1 variable per vertex, maximises their sum, and requires x_u + x_v of at most 1 for every non-adjacent pair. It is weak: every x at one half is feasible, so the LP bound is n/2. The fix mirrors the colouring bound: one constraint per independent set, such as each colour class, allowing at most one of its vertices.

# pip install pulp   (any MILP solver works; this shows the tightened formulation)
import pulp

def clique_ilp(n, edges, colour_classes):
    E = {frozenset(e) for e in edges}
    m = pulp.LpProblem("max_clique", pulp.LpMaximize)
    x = [pulp.LpVariable(f"x{v}", cat="Binary") for v in range(n)]
    m += pulp.lpSum(x)
    for S in colour_classes:                    # independent sets: at most one vertex
        m += pulp.lpSum(x[v] for v in S) <= 1
    for u in range(n):                          # remaining non-edges, for correctness
        for v in range(u + 1, n):
            if frozenset((u, v)) not in E and not any(u in S and v in S for S in colour_classes):
                m += x[u] + x[v] <= 1
    m.solve(pulp.PULP_CBC_CMD(msg=False, timeLimit=60))
    return [v for v in range(n) if x[v].value() > 0.5]

Reach for the ILP when the clique is one part of a larger model: vertex weights, side constraints such as at most two members from one region, or a combined objective. For the pure problem, a dedicated combinatorial solver is usually much faster, though the ILP's dual bound tells you how far a timed-out incumbent might be from optimal.

Failure modes

Exact clique search fails in a few recognisable ways.

• Dense, mid-sized, structureless graphs. Random graphs with a few hundred vertices at density 0.9 or higher, and the DIMACS benchmark families built to be hard, can run for hours or longer. The colouring bound is loose when the clique number is far below the chromatic number, and no amount of tuning fixes that gap.
• Unbounded runtime in production. Always run with a deadline and report the incumbent plus the best remaining bound, the maximum of size plus colour over unexplored root branches, so callers know whether the answer is proven.
• Recursion depth. Depth equals clique size; a planted clique of thousands of vertices exceeds Python's default limit, so use an explicit stack.
• Silent wrong answers. A bound bug returns a smaller clique, not an error. Always verify the returned set is pairwise adjacent; that check catches invalid output but not suboptimal output, which is why the brute-force cross-check is part of the test suite.
• Wrong graph. Cliques in thresholded graphs, such as correlation above 0.8, are very sensitive to the threshold; report how the answer moves with it.

Trade-offs and graph classes

Graph class matters as much as algorithm. Chordal and interval graphs have polynomial-time maximum clique via a perfect elimination ordering, bipartite graphs have clique number at most 2, and planar graphs at most 4. The chromatic number is at least the clique number, so a large clique is a lower bound for graph colouring, and a good colouring is an upper bound for clique. If the problem arrives as constraints rather than a graph, a SAT or MaxSAT solver is a reasonable alternative, with clique-specific encodings in the literature.

What to do next

1. Write down whether you need the maximum, all maximal cliques, or a good clique fast. That one decision picks the algorithm.
2. Check the graph class first. Bipartite, chordal, interval or planar inputs have polynomial-time answers.
3. Compute a degeneracy-order heuristic incumbent and run k-core pruning before any exact search, and log how many vertices survive.
4. Implement or adopt a bitset colour-bound solver, and cross-check it against brute force on random graphs of up to 12 vertices.
5. Put a deadline on every run and return the incumbent with its upper bound, so callers know whether it is proven.
6. Switch to an ILP with independent-set cuts only when side constraints or weights require it.