Как решить задачу покрытия множеств с помощью MILP в SciPy

Показываем, как свести задачу покрытия множеств к MILP и решить её в Python с SciPy: бинарные переменные, разреженная матрица ограничений и масштабирование.

Покрыть конечный универсум минимальным числом подмножеств кажется обманчиво простой задачей, пока не попробуешь решить её грубой силой. Нужно выбрать минимальное количество ключей из словаря множеств так, чтобы их объединение совпадало с полным универсумом; если такого набора нет, результат должен быть пустым. На небольших входах перебор всех комбинаций срабатывает, но по мере того как число ключей приближается к сотням, этот подход быстро становится непрактичным.

Постановка задачи и базовый перебор

Рассмотрим отображение, значения которого — подмножества общего множества. Цель — выбрать наименьшее число ключей, чьи множества покрывают универсум.

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']

Эта исчерпывающая рекурсия перебирает все комбинации ключей, отбирает те, чьё объединение равно универсуму, и возвращает минимальную. Это работает, но стоимость растёт взрывным образом с увеличением числа ключей.

Что на самом деле происходит

Эта задача — частный случай проблемы покрытия множеств. Требуется минимизировать количество выбранных подмножеств при условии, что каждый элемент универсума покрыт хотя бы одним из них. Известно, что задача вычислительно трудная; следовательно, поиск по всем комбинациям не масштабируется.

Практический подход: MILP с разреженной матрицей покрытия

Для входов такого масштаба компактный и быстрый способ — сформулировать её как смешанно-целочисленную линейную задачу (MILP). Каждый ключ превращается в бинарную переменную. Целевая функция минимизирует число переменных, равных единице. Ограничения требуют, чтобы каждый элемент универсума был покрыт хотя бы одним выбранным множеством. Формирование этих ограничений в виде разреженной матрицы сохраняет эффективность.

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

Результат — покрытие минимального размера для заданного отображения.

Масштабирование до примерно сотни ключей

Тот же подход применим, когда словарь растёт. Следующий фрагмент генерирует пример со ста ключами и решает его тем же способом.

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]

На типичном оборудовании входы с примерно 1200 ненулевыми элементами в матрице ограничений обрабатываются примерно за секунду.

Почему это важно

Хотя задача покрытия множеств вычислительно трудна, корректно поставленная MILP с разреженным представлением позволяет быстро решать реалистичные экземпляры. Вместо перебора всех комбинаций вы делегируете поиск оптимизатору, который эксплуатирует структуру задачи, отсечение и эффективную линейную алгебру «под капотом». Результат воспроизводим, оптимален и надёжен для показанных размеров.

Выводы

Если вам нужна минимальная коллекция ключей, покрывающая универсум, рассматривайте её как set cover и выражайте в виде MILP с бинарными переменными, набором простых ограничений «покрыть каждый элемент хотя бы раз» и суммирующей минимизируемой целевой функцией. Для продемонстрированных размеров формирование разреженной матрицы ограничений и запуск решателя дают оптимальные ответы в пределах от миллисекунд до секунд. Когда требуется любая одна оптимальная комбинация, решение MILP естественным образом именно её и предоставляет.