Set Cover in Python: From Brute Force to Optimal MILP with a Sparse Constraint Matrix

Learn how to solve the set cover problem in Python using SciPy's MILP and a sparse matrix model. Avoid brute force and get optimal, scalable covers fast.

Covering a finite universe with the fewest subsets looks deceptively simple, until you try to brute-force your way through. The task is to select the minimum number of keys from a dict of sets so that their union equals the full universe; if no such selection exists, the result should be empty. For small inputs, enumerating all combinations works, but as the number of keys approaches the hundreds, that approach quickly becomes impractical.

Problem setup and a brute-force baseline

Consider a mapping whose values are subsets of a total set. The goal is to pick the smallest number of keys whose sets cover the universe.

import copy

pool = {'a': {1,2,8}, 'b': {3,1,2,6}, 'c': {0,4,1,2}, 'd': {9}, 'e': {2,5},
        'f': {4,8}, 'g': {0,9}, 'h': {7,2,3}, 'i': {5,6,3}, 'j': {4,6,8}}

universe = set(range(10))

picks_ok = []

def walk(current, remaining):
    for idx in range(len(remaining)):
        next_pick = copy.copy(current)
        next_pick.append(remaining[idx])
        rest_labels = remaining[idx+1:]
        if {elem for name in next_pick for elem in pool[name]} == universe:
            picks_ok.append(next_pick)
        walk(next_pick, rest_labels)

start_pick = []
labels = sorted(pool.keys())
walk(start_pick, labels)
answer = sorted(picks_ok, key=lambda r: len(r))
print(answer[0])
# ['a', 'c', 'd', 'h', 'i']

This exhaustive recursion tries every combination of keys, collects those whose union equals the universe, and returns the smallest. It works, but its cost grows explosively with the number of keys.

What is really happening

The task is an instance of the set cover problem. The objective is to minimize how many subsets are chosen while ensuring that every element of the universe is covered by at least one chosen subset. This problem is known to be computationally intensive; consequently, searching all combinations does not scale.

A practical approach: MILP with a sparse cover matrix

For inputs of this size, a compact and fast way to solve the problem is to formulate it as a Mixed-Integer Linear Program (MILP). Each key becomes a binary decision variable. The objective minimizes the number of variables set to one. Constraints enforce that each element of the universe is covered by at least one selected set. Building these constraints as a sparse matrix keeps it efficient.

import time

import numpy as np
from scipy.optimize import milp, Bounds, LinearConstraint
from scipy.sparse import coo_array

groups = {
    'a': {1,2,8}, 'b': {3,1,2,6}, 'c': {0,4,1,2}, 'd': {9}, 'e': {2,5},
    'f': {4,8}, 'g': {0,9}, 'h': {7,2,3}, 'i': {5,6,3}, 'j': {4,6,8},
}

u_size = 1 + max(val for subset in groups.values() for val in subset)
k_count = len(groups)

objective = np.ones(k_count)
var_bounds = Bounds(lb=0, ub=1)

nz = np.array(
    [
        (row, col)
        for row in range(u_size)
        for col, subset in enumerate(groups.values())
        if row in subset
    ],
    dtype=np.int32,
).T

covering = LinearConstraint(
    A=coo_array((
        np.ones(nz.shape[1]),
        nz,
    )).tocsc(),
    lb=1,
)

t0 = time.perf_counter()
sol = milp(c=objective, integrality=1, bounds=var_bounds, constraints=covering)
t1 = time.perf_counter()
assert sol.success

print(sol.message)
print(f'Cost: {sol.fun} in {(t1 - t0)*1e3:.2f} ms')
print(
    'Assigned:',
    ' '.join(
        name
        for name, flag in zip(groups.keys(), sol.x)
        if flag
    ),
)

Optimization terminated successfully. (HiGHS Status 7: Optimal)
Cost: 5.0 in 2.58 ms
Assigned: c g h i j

The result is a minimum-size cover for the provided mapping.

Scaling to around a hundred keys

The same idea applies when the dictionary grows. The following snippet generates a case with one hundred keys and solves it in the same way.

import time

import numpy as np
from scipy.optimize import milp, Bounds, LinearConstraint
from scipy.sparse import coo_array

rng = np.random.default_rng(seed=0)
u_size = 10
k_count = 100

# Random instance: each key maps to 1..4 unique elements in [0, u_size)
sets_rand = {
    i: rng.choice(a=u_size, size=rng.integers(low=1, high=4, endpoint=True),
                  replace=False, shuffle=False)
    for i in range(k_count)
}

objective = np.ones(k_count)
var_bounds = Bounds(lb=0, ub=1)

nz = np.array(
    [
        (row, col)
        for row in range(u_size)
        for col, subset in enumerate(sets_rand.values())
        if row in subset
    ],
    dtype=np.int32,
).T

covering = LinearConstraint(
    A=coo_array((np.ones(nz.shape[1]), nz)).tocsc(),
    lb=1,
)

t0 = time.perf_counter()
sol = milp(c=objective, integrality=1, bounds=var_bounds, constraints=covering)
t1 = time.perf_counter()

print(sol.message)
print(f'Cost: {sol.fun} in {(t1 - t0)*1e3:.2f} ms')
print(
    'Assigned:',
    ' '.join(
        str(k)
        for k, flag in zip(sets_rand.keys(), sol.x)
        if flag > 0.5
    ),
)
for (k, subset), flag in zip(sets_rand.items(), sol.x):
    if flag > 0.5:
        print(f'{k}: {subset}')

Optimization terminated successfully. (HiGHS Status 7: Optimal)
Cost: 3.0 in 18.49 ms
Assigned: 0 82 99
0: [4 7 2 3]
82: [6 0 5 8]
99: [1 5 4 9]

On typical hardware, inputs around twelve hundred nonzero entries in the constraint matrix are handled in about a second.

Why this matters

Even though the set cover problem is computationally intensive, a well-posed MILP with a sparse representation can solve realistic instances quickly. Instead of exploring all combinations, you delegate the search to an optimizer that exploits structure, pruning, and efficient linear algebra under the hood. The result is reproducible, optimal, and robust for the problem sizes shown.

Takeaways

If you need the smallest collection of keys that covers a universe, treat it as set cover and express it as a MILP with binary variables, a simple “cover each element at least once” constraint set, and a sum-minimizing objective. For the illustrated sizes, building a sparse constraint matrix and calling a solver yields optimal answers within milliseconds to seconds. When you only need any one optimal combination, the MILP solution naturally provides exactly that.