python - 在二进制矩阵中精确找到 k 列的快速算法,使得这些列的总和是 1 向量
问题描述
假设我有一个 (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,则不应报告任何解。
我尝试了两种方法:
- 我尝试所有(N 选择 k)组合并找到总和为 1 向量的组合的慢速组合方法。这在 O(N^k) 时间内运行,这显然是可怕的。
- 一种递归方法,速度更快,但仍然在 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 上问过同样的问题,但我不太关心算法效率,更关心数学技术。)
解决方案
我的第一次尝试是使用一个可用的高性能求解器(例如Cbc)进行整数编程。
假设您的关联矩阵中有一些稀疏性,这些将非常有效并且非常普遍(侧面约束/适应)。它们也是完整的,并且可能能够证明不可行。
一个简单的公式如下所示:
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
- cvxpy 的非常代数的使用以及一些接口决策导致了不寻常的变量边界作为这里的约束公式
- 唯一易于使用的非商业求解器是