比抛硬币游戏的蛮力算法更好
Better than brute force algorithms for a coin-flipping game
我有一个问题，我觉得应该有一个比蛮力更好的著名算法来解决它，但我想不出一个，所以我在这里问。
问题如下：给定n个排序（从低到高）包含m概率的列表，为每个列表选择一个索引这样所选索引的总和小于 m。然后，对于每个列表，我们抛一枚硬币，硬币正面朝上的机会等于该列表所选索引处的概率。最大限度地提高硬币正面朝上至少一次的机会。
是否有解决此问题的算法比仅靠蛮力更好？
这个问题看起来与背包问题最相似，只是背包中物品的价值不仅仅是背包中物品的总和。 （写在 Python 中，而不是 sum(p for p in chosen_probabilities)，而是 1 - math.prod([1 - p for p in chosen_probabilities])）而且，考虑到背包中已有的物品，您可以添加的物品有限制。例如，如果特定列表的 index = 3 项目已经在背包中，则不允许为同一列表添加带有 index = 2 的项目（因为您只能为每个列表选择一个索引列表）。因此，根据背包中已有的物品，某些物品可以和不能添加到背包中。
线性优化不起作用，因为列表中的值不是线性增加的，最终硬币概率与所选概率不是线性关系，我们的约束是总和 正如 David 指出的那样，线性优化 将 如果您使用二进制变量来挑选索引和处理非线性的对数。
编辑：
我发现解释这个问题背后的动机有助于理解它。想象一下，你有 10 秒的时间来解决一个问题，并用三种不同的方法来解决它。考虑到您尝试该方法的秒数，您拥有每种方法解决问题的可能性的模型，但如果您切换方法，您将失去之前尝试的方法的所有进展。您应该尝试哪些方法以及多长时间？
最大化1 - math.prod([1 - p for p in chosen_probabilities])相当于最小化math.prod([1 - p for p in chosen_probabilities])，相当于最小化这个objective的对数，是0-1指标变量的线性函数，所以你可以用这种方式做一个整数规划公式。
我不能保证这会比蛮力好得多。问题是当 p 接近于零时，math.log(1 - p) 很接近于 -p。我的直觉是，对于非平凡的实例，它在性质上类似于使用整数规划来求解子集和，但效果不是特别好。
如果您愿意接受双标准近似方案（得到一个答案，使所选索引的总和小于 m，这至少与总和小于 (1 − ε) m) 那么您可以将概率四舍五入为 ε 的倍数，并使用动态规划得到一个在 n, m, 1/ε 的时间多项式中运行的算法。
这是 David Eisenstat's solution 的工作代码。
要理解实现，我认为先完成数学运算会有所帮助。
提醒一下，有 n 个列表，每个列表都有 m 个选项。 （在问题底部的激励示例中，每个列表代表一种解决问题的方法，你有 m-1 秒来解决问题。每个列表都是 list[index] 给机会如果该方法 运行 持续 index 秒，则可以使用该方法解决问题。）
我们让列表存储在一个名为d（代码中名为data）的矩阵中，矩阵中的每一行都是一个列表。 （因此每一列代表一个索引，或者，如果遵循激励示例，则代表一段时间。）
假定我们为列表 i 选择索引 j*，则硬币正面朝上的概率计算为
)
我们希望将其最大化。
（为了解释这个等式背后的统计数据，我们计算 1 减去硬币没有正面朝上的概率。硬币没有正面朝上的概率是每次抛硬币没有正面朝上的概率不会正面朝上。单次翻转没有正面朝上的概率只是 1 减去正面朝上的概率。而它确实正面朝上的概率是我们选择的数字，d[i][j*]。因此，所有抛掷落在反面的总概率只是每个抛掷落在反面的概率的乘积。然后硬币落在正面的概率只是 1 减去所有抛掷落在反面的概率尾巴。）
正如大卫指出的那样，这与最小化相同：
)
这与最小化相同：
&space;%5Cright))
相当于：
)
那么，因为这是线性和，我们可以把它变成一个整数规划。
我们将最小化：
&space;x_%7Bi,&space;j%7D)
这允许计算机通过允许它创建一个 n by m 1 和 0 的矩阵来选择索引，称为 x，其中 1 挑选特定索引。然后我们将定义规则，以便它不会挑选出无效的索引集。
第一个规则是你必须为每个列表选择一个索引：
第二条规则是您必须遵守所选索引总和必须等于或小于 m 的约束：
就是这样！然后我们可以告诉计算机根据这些规则最小化该总和。它会吐出一个 x 矩阵，每行有一个 1 来告诉我们它为该行的列表选择了哪个索引。
在代码中（使用激励示例），实现为：
'''
Requirements:
cvxopt==1.2.6
cvxpy==1.1.10
ecos==2.0.7.post1
numpy==1.20.1
osqp==0.6.2.post0
qdldl==0.1.5.post0
scipy==1.6.1
scs==2.1.2
'''
import math
import cvxpy as cp
import numpy as np
# number of methods
n = 3
# if you have 10 seconds, there are 11 options for each method (0 seconds, 1 second, ..., 10 seconds)
m = 11
# method A has 30% chance of working if run for at least 3 seconds
# equivalent to [0, 0, 0, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3]
A_list = [0, 0, 0] + [0.3] * (m - 3)
# method B has 30% chance of working if run for at least 3 seconds
# equivalent to [0, 0, 0, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3]
B_list = [0, 0, 0] + [0.3] * (m - 3)
# method C has 40% chance of working if run for 4 seconds, 30% otherwise
# equivalent to [0.3, 0.3, 0.3, 0.3, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4]
C_list = [0.3, 0.3, 0.3, 0.3] + [0.4] * (m - 4)
data = [A_list, B_list, C_list]
# do the logarithm
log_data = []
for row in data:
log_row = []
for col in row:
# deal with domain exception
if col == 1:
new_col = float('-inf')
else:
new_col = math.log(1 - col)
log_row.append(new_col)
log_data.append(log_row)
log_data = np.array(log_data)
x = cp.Variable((n, m), boolean=True)
objective = cp.Minimize(cp.sum(cp.multiply(log_data, x)))
# the current solver doesn't work with equalities, so each equality must be split into two inequalities.
# see https://github.com/cvxgrp/cvxpy/issues/1112
one_choice_per_method_constraint = [cp.sum(x[i]) <= 1 for i in range(n)] + [cp.sum(x[i]) >= 1 for i in range(n)]
# constrain the solution to not use more time than is allowed
# note that the time allowed is (m - 1), not m, because time is 1-indexed and the lists are 0-indexed
js = np.tile(np.array(list(range(m))), (n, 1))
time_constraint = [cp.sum(cp.multiply(js, x)) <= m - 1, cp.sum(cp.multiply(js, x)) >= m - 1]
constraints = one_choice_per_method_constraint + time_constraint
prob = cp.Problem(objective, constraints)
result = prob.solve()
def compute_probability(data, choices):
# compute 1 - ((1 - p1) * (1 - p2) * ...)
return 1 - np.prod(np.add(1, -np.multiply(data, choices)))
print("Choices:")
print(x.value)
'''
Choices:
[[0. 0. 0. 1. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 1. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 1. 0. 0. 0. 0. 0. 0.]]
'''
print("Chance of success:")
print(compute_probability(data, x.value))
'''
Chance of success:
0.7060000000000001
'''
我们有它！计算机已经正确判定运行宁方法A 3秒，方法B 3秒，方法C 4秒是最优的。 （请记住，x 矩阵是 0 索引的，而时间是 1 索引的。）
谢谢大卫的建议！
我有一个问题，我觉得应该有一个比蛮力更好的著名算法来解决它，但我想不出一个，所以我在这里问。
问题如下：给定n个排序（从低到高）包含m概率的列表，为每个列表选择一个索引这样所选索引的总和小于 m。然后，对于每个列表，我们抛一枚硬币，硬币正面朝上的机会等于该列表所选索引处的概率。最大限度地提高硬币正面朝上至少一次的机会。
是否有解决此问题的算法比仅靠蛮力更好？
这个问题看起来与背包问题最相似，只是背包中物品的价值不仅仅是背包中物品的总和。 （写在 Python 中，而不是 sum(p for p in chosen_probabilities)，而是 1 - math.prod([1 - p for p in chosen_probabilities])）而且，考虑到背包中已有的物品，您可以添加的物品有限制。例如，如果特定列表的 index = 3 项目已经在背包中，则不允许为同一列表添加带有 index = 2 的项目（因为您只能为每个列表选择一个索引列表）。因此，根据背包中已有的物品，某些物品可以和不能添加到背包中。
线性优化不起作用，因为列表中的值不是线性增加的，最终硬币概率与所选概率不是线性关系，我们的约束是总和 正如 David 指出的那样，线性优化 将 如果您使用二进制变量来挑选索引和处理非线性的对数。
编辑：
我发现解释这个问题背后的动机有助于理解它。想象一下，你有 10 秒的时间来解决一个问题，并用三种不同的方法来解决它。考虑到您尝试该方法的秒数，您拥有每种方法解决问题的可能性的模型，但如果您切换方法，您将失去之前尝试的方法的所有进展。您应该尝试哪些方法以及多长时间？
最大化1 - math.prod([1 - p for p in chosen_probabilities])相当于最小化math.prod([1 - p for p in chosen_probabilities])，相当于最小化这个objective的对数，是0-1指标变量的线性函数，所以你可以用这种方式做一个整数规划公式。
我不能保证这会比蛮力好得多。问题是当 p 接近于零时，math.log(1 - p) 很接近于 -p。我的直觉是，对于非平凡的实例，它在性质上类似于使用整数规划来求解子集和，但效果不是特别好。
如果您愿意接受双标准近似方案（得到一个答案，使所选索引的总和小于 m，这至少与总和小于 (1 − ε) m) 那么您可以将概率四舍五入为 ε 的倍数，并使用动态规划得到一个在 n, m, 1/ε 的时间多项式中运行的算法。
这是 David Eisenstat's solution 的工作代码。
要理解实现，我认为先完成数学运算会有所帮助。
提醒一下，有 n 个列表，每个列表都有 m 个选项。 （在问题底部的激励示例中，每个列表代表一种解决问题的方法，你有 m-1 秒来解决问题。每个列表都是 list[index] 给机会如果该方法 运行 持续 index 秒，则可以使用该方法解决问题。）
我们让列表存储在一个名为d（代码中名为data）的矩阵中，矩阵中的每一行都是一个列表。 （因此每一列代表一个索引，或者，如果遵循激励示例，则代表一段时间。）
假定我们为列表 i 选择索引 j*，则硬币正面朝上的概率计算为
我们希望将其最大化。
（为了解释这个等式背后的统计数据，我们计算 1 减去硬币没有正面朝上的概率。硬币没有正面朝上的概率是每次抛硬币没有正面朝上的概率不会正面朝上。单次翻转没有正面朝上的概率只是 1 减去正面朝上的概率。而它确实正面朝上的概率是我们选择的数字，d[i][j*]。因此，所有抛掷落在反面的总概率只是每个抛掷落在反面的概率的乘积。然后硬币落在正面的概率只是 1 减去所有抛掷落在反面的概率尾巴。）
正如大卫指出的那样，这与最小化相同：
这与最小化相同：
相当于：
那么，因为这是线性和，我们可以把它变成一个整数规划。
我们将最小化：
这允许计算机通过允许它创建一个 n by m 1 和 0 的矩阵来选择索引，称为 x，其中 1 挑选特定索引。然后我们将定义规则，以便它不会挑选出无效的索引集。
第一个规则是你必须为每个列表选择一个索引：
第二条规则是您必须遵守所选索引总和必须等于或小于 m 的约束：
就是这样！然后我们可以告诉计算机根据这些规则最小化该总和。它会吐出一个 x 矩阵，每行有一个 1 来告诉我们它为该行的列表选择了哪个索引。
在代码中（使用激励示例），实现为：
'''
Requirements:
cvxopt==1.2.6
cvxpy==1.1.10
ecos==2.0.7.post1
numpy==1.20.1
osqp==0.6.2.post0
qdldl==0.1.5.post0
scipy==1.6.1
scs==2.1.2
'''
import math
import cvxpy as cp
import numpy as np
# number of methods
n = 3
# if you have 10 seconds, there are 11 options for each method (0 seconds, 1 second, ..., 10 seconds)
m = 11
# method A has 30% chance of working if run for at least 3 seconds
# equivalent to [0, 0, 0, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3]
A_list = [0, 0, 0] + [0.3] * (m - 3)
# method B has 30% chance of working if run for at least 3 seconds
# equivalent to [0, 0, 0, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3]
B_list = [0, 0, 0] + [0.3] * (m - 3)
# method C has 40% chance of working if run for 4 seconds, 30% otherwise
# equivalent to [0.3, 0.3, 0.3, 0.3, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4]
C_list = [0.3, 0.3, 0.3, 0.3] + [0.4] * (m - 4)
data = [A_list, B_list, C_list]
# do the logarithm
log_data = []
for row in data:
log_row = []
for col in row:
# deal with domain exception
if col == 1:
new_col = float('-inf')
else:
new_col = math.log(1 - col)
log_row.append(new_col)
log_data.append(log_row)
log_data = np.array(log_data)
x = cp.Variable((n, m), boolean=True)
objective = cp.Minimize(cp.sum(cp.multiply(log_data, x)))
# the current solver doesn't work with equalities, so each equality must be split into two inequalities.
# see https://github.com/cvxgrp/cvxpy/issues/1112
one_choice_per_method_constraint = [cp.sum(x[i]) <= 1 for i in range(n)] + [cp.sum(x[i]) >= 1 for i in range(n)]
# constrain the solution to not use more time than is allowed
# note that the time allowed is (m - 1), not m, because time is 1-indexed and the lists are 0-indexed
js = np.tile(np.array(list(range(m))), (n, 1))
time_constraint = [cp.sum(cp.multiply(js, x)) <= m - 1, cp.sum(cp.multiply(js, x)) >= m - 1]
constraints = one_choice_per_method_constraint + time_constraint
prob = cp.Problem(objective, constraints)
result = prob.solve()
def compute_probability(data, choices):
# compute 1 - ((1 - p1) * (1 - p2) * ...)
return 1 - np.prod(np.add(1, -np.multiply(data, choices)))
print("Choices:")
print(x.value)
'''
Choices:
[[0. 0. 0. 1. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 1. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 1. 0. 0. 0. 0. 0. 0.]]
'''
print("Chance of success:")
print(compute_probability(data, x.value))
'''
Chance of success:
0.7060000000000001
'''
我们有它！计算机已经正确判定运行宁方法A 3秒，方法B 3秒，方法C 4秒是最优的。 （请记住，x 矩阵是 0 索引的，而时间是 1 索引的。）
谢谢大卫的建议！