Быстрый алгоритм для поиска ровно k столбцов в двоичной матрице, так что сумма этих столбцов является 1-вектором

Question

Предположим, у меня есть двоичная матрица (M x N), в которой M и N могут быть большими. Я хочу найти ровно k столбцов (k относительно мало, скажем, меньше 10), чтобы сумма этих k столбцов была 1-вектором (все элементы равны 1). Достаточно одного решения. Есть ли для этого быстрый алгоритм?

Например, алгоритм, работающий с матрицей

1 0 0
1 0 0
1 1 0
0 1 1

с k = 2, должен возвращать столбцы 0 и 2, но не должен сообщать никаких решений, если k = 1 или k = 3.

Я пробовал два подхода:

1. Медленный комбинаторный подход, когда я пробую все (N выбирают k) комбинации и нахожу комбинацию, которая в сумме дает 1-вектор. Это выполняется за время O (N ^ k), что, очевидно, ужасно.
2. Рекурсивный подход, который быстрее, но все же выполняется за время O (N ^ k) в худшем случае. Код Python выглядит следующим образом:

import numpy as np

def recursiveFn(mat, col_used_bool, col_sum_to_date, cols_to_go):
    N = len(mat)
    if cols_to_go == 1:
        col_unused = 1 - col_sum_to_date
        if list(col_unused) in [list(i) for i in mat]:
            return (True, [col_unused])
        else:
            return (False, None)
    for col_id in range(N):
        if col_used_bool[col_id]:
            continue
        if 2 not in mat[col_id]+col_sum_to_date:
            col_used_bool[col_id] = True
            x = recursiveFn(mat, col_used_bool, mat[col_id]+col_sum_to_date, cols_to_go-1)
            col_used_bool[col_id] = False
            if x[0]:
                return (True, x[1] + [mat[col_id]])
    return (False, None)

exMat = [[1,1,1,0],[0,0,1,1],[0,0,0,1]] #input by colums
exMat = [np.asarray(i) for i in exMat]
k = 2
output = recursiveFn(mat = exMat, col_used_bool = [False for i in exMat], 
    col_sum_to_date = np.asarray([0 for i in exMat[0]]), cols_to_go = k)
print(output[1])
### prints this : [array([0, 0, 0, 1]), array([1, 1, 1, 0])]

Меня не устраивает ни один из этих подходов, и я чувствую, что существует более умный и быстрый алгоритм. Большое спасибо за вашу помощь. Это мой первый пост на StackOverflow, поэтому, пожалуйста, будьте осторожны со мной, если я где-то ошибся!

(Если интересно, я также задал тот же вопрос на Math Stack Exchange , но там меня меньше беспокоит эффективность алгоритмов c и больше интересуют математические методы.)

Answer 1

Моей первой попыткой было бы целочисленное программирование с использованием одного из доступных высокопроизводительных решателей (например, Cb c).

Предполагая некоторая разреженность в вашей матрице инцидентности, они будут очень эффективными и довольно общими (боковые ограничения / адаптации). Они также являются полными и могут доказать несостоятельность.

Простая формулировка будет выглядеть так:

Instance

c0 c1 c2
1  0  0  r0
1  0  0  r1
1  1  0  r2
0  1  1  r3

IP:

minimize(0)        # constant objective | pure feasibility problem

sum(c_i) = k       # target of columns chosen

r0 = 1 = c0        # r0 just showing the origin of the constraint; no real variable!
r1 = 1 = c0
r2 = 1 = c0 + c1
r3 = 1 = c1 + c2

c_i in {0, 1}      # all variables are binary

Можно было бы усилить эту формулировку дополнительными неравенствами, такими как неравенства клик ( конфликт-граф -> максимальные клики), но не уверен, помогает ли это. Хорошие решатели будут делать нечто подобное динамически генерировать разрезов .

Существует много теории. Одно ключевое слово может быть точное покрытие или все те проблемы упаковки / покрытия, которые очень похожи.

Простой пример кода:

import cvxpy as cp
import numpy as np

data = np.array([[1, 0, 0],
                 [1, 0, 0],
                 [1, 1, 0],
                 [0, 1, 1]])

def solve(k, data):
  c = cp.Variable(data.shape[1], boolean=True)

  con = [data * c == 1,
         cp.sum(c) == k,
         c >= 0,
         c <= 1]

  obj = cp.Minimize(0)
  
  problem = cp.Problem(obj, con)
  problem.solve(verbose=True, solver=cp.GLPK_MI)

  if(problem.status == 'optimal'):
    return np.where(np.isclose(c.value, 1.0) == True)[0]
  else:
    assert problem.status == 'infeasible'
    return None

print(solve(2, data))
print(solve(1, data))
print(solve(3, data))

# [0 2]
# None
# None

Примечания:

• Код использует cvxpy, который очень мощный, но не имеет некоторой расширенной поддержки целочисленного программирования.
  • Единственный простой в использовании некоммерческий решатель - GLPK, что очень хорошо, но обычно не может конкурировать с Cbc
  • Само алгебраическое c использование cvxpy вместе с некоторыми интерфейсными решениями приводит к необычным ограничениям переменных формулировке здесь

Answer 2

Как упоминалось в первом ответе, это проблема точного покрытия , которая является NP-сложной. Классический способ решения NP-сложной проблемы - это откат. Различные реализации могут дать совершенно разные результаты. здесь стоит протестировать алгоритм.

Однако из-за того, что нужно выбрать только небольшое количество k столбцов, я бы попробовал другой подход, т.е. классический алгоритм поиска с возвратом с логическим b[j] указывает, выбран ли столбец j, с двумя дополнительными приемами.

1. При добавлении столбца j к текущей сумме столбцов мы можем остановить процесс, как только появится " 2 ", нам не нужно ждать, пока будет вычислена окончательная сумма.

2. Вместо того, чтобы добавлять элементы столбца по одному, мы можем сгруппировать p элементов (соответствующих до p строк) каждого столбца в одно целое число, чтобы ускорить процесс суммирования столбцов. Для этого нужно выбрать базу. Небольшая база позволяет избежать слишком больших чисел (это важно для ограничения размера массива isValid [], см. Ниже).

База 2 невозможна: например, добавление (1 0) и (1 0) даст (0 1), которое по-прежнему является допустимым числом.

Поэтому я предлагаю использовать base 3 , что позволяет обнаруживать наличие ошибочной цифры «2» при суммировании. Например,

V(0 1 1 0) = 0*3**0 + 1*3**1 +1*3**2 + 0*3**3

На практике для анализа групп элементов «p» нам нужна логическая таблица размера «3 ** p», isValid[], которая позволит сразу определить, есть ли данное полученное целое число действительно. Эта таблица должна быть предварительно обработана во время фазы инициализации.

Мы знаем, что мы получили 1-вектор, когда все целые числа равны заданному c значению (3**p - 1)/2, учитывая, что последняя группа может иметь другой размер p' < p.

Из-за большого значения n можно проверить последний трюк:

Найдите допустимые решения для ряда строк n1 < n, а затем для каждого полученного решения-кандидата проверьте, действительно ли оно является решением для всех n строк.