比抛硬币游戏的蛮力算法更好

Question

我有一个问题，我觉得应该有一个众所周知的算法来解决它，它比蛮力更好，但我想不出一个，所以我在这里问。

问题如下：给定包含m概率的n排序（从低到高）列表，为每个列表选择一个索引，使得所选索引的总和小于比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 秒的时间来解决一个问题，并用三种不同的方法来解决它。考虑到您尝试该方法的秒数，您拥有每种方法解决问题的可能性的模型，但是如果您切换方法，您将失去之前尝试的方法的所有进度。您应该尝试哪些方法以及尝试多长时间？

Answer 1

最大化1 - math.prod([1 - p for p in chosen_probabilities])等于最小化math.prod([1 - p for p in chosen_probabilities])，也就是最小化这个目标的log，它是0-1个指标变量的线性函数，所以可以这样做整数规划公式。

我不能保证这会比蛮力好得多。问题是当p 接近于零时，math.log(1 - p) 很好地近似于-p。我的直觉是，对于非平凡的实例，它在质量上类似于使用整数规划来解决子集和，但并不是特别好。

如果您愿意满足于双准则近似方案（得到一个答案，使得所选索引的总和小于 m，那至少与总和小于 (1 - ε) 的最佳答案一样好) m) 然后你可以将概率四舍五入为 ε 的倍数，并使用动态规划得到一个算法，该算法在 n, m, 1/ε 的时间多项式中运行。

Answer 2

这是David Eisenstat's solution 的工作代码。

要理解实现，我认为先进行数学运算会有所帮助。

提醒一下，有n 列表，每个列表都有m 选项。 （在问题底部的激励示例中，每个列表代表一种解决问题的方法，并且给您m-1 秒来解决问题。每个列表都使得list[index] 给出解决问题的机会如果该方法运行index 秒，则使用该方法。）

我们将列表存储在一个名为d（在代码中命名为data）的矩阵中，其中矩阵中的每一行都是一个列表。 （因此，每一列代表一个索引，或者，如果遵循激励示例，则代表一个时间量。）

假设我们为列表i 选择索引j*，硬币正面朝上的概率计算为

我们希望将其最大化。

（为了解释这个等式背后的统计数据，我们计算 1 减去硬币不正面朝上的概率。硬币不正面朝上的概率是每次抛硬币不正面朝上的概率正面朝上。单次翻转不正面朝上的概率只是 1 减去正面朝上的概率。它正面朝上的概率就是我们选择的数字，d[i][j*]。因此，所有翻转落在尾巴上的总概率只是每个翻转落在尾巴上的概率的乘积。然后硬币落在正面的概率只是 1 减去所有翻转落在尾巴上的概率。）

正如大卫所指出的，这与最小化相同：

这与最小化相同：

相当于：

那么，既然这是线性求和，我们可以把它变成一个整数程序。

我们将最小化：

这让计算机通过允许它创建一个由 1 和 0 组成的 n 和 m 矩阵来选择索引，称为 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 索引的。）

谢谢大卫的建议！