0-1 背包问题在量子退火机上的简单实现

写在前面

平台、配置、接口都不是一成不变的，不要拘泥于某一版本或者某一 API，要是有更方便的呢？

本文主要来源自 D-Wave 官方文档，如有疏忽，请参考原文档。

（不要全文抄袭当作课程作业或者调研报告）

本文配置

• Windows 10
• Python==3.8.5
• dwave-ocean-sdk==3.1.1

如果使用 Leap 在线平台的求解器，一定要配置（前提是已经注册Leap leap），指令 dwave setup。你会看到

$ dwave setup

Optionally install non-open-source packages and configure your environment.

Do you want to select non-open-source packages to install (y/n)? [y]:

D-Wave Drivers
These drivers enable some automated performance-tuning features.
This package is available under the 'D-Wave EULA' license.
The terms of the license are available online: https://docs.ocean.dwavesys.com/eula
Install (y/n)? [y]:
Installing: D-Wave Drivers
Successfully installed D-Wave Drivers.

D-Wave Problem Inspector
This tool visualizes problems submitted to the quantum computer and the results returned.
This package is available under the 'D-Wave EULA' license.
The terms of the license are available online: https://docs.ocean.dwavesys.com/eula
Install (y/n)? [y]:
Installing: D-Wave Problem Inspector
Successfully installed D-Wave Problem Inspector.

Creating the D-Wave configuration file.
Configuration file not found; the default location is: /home/someone/.config/dwave/dwave.conf
Confirm configuration file path [/home/someone/.config/dwave/dwave.conf]:
Profile (create new) [prod]:
API endpoint URL [skip]:
Authentication token [skip]: ABC-1234567890abcdef1234567890abcdef
Default client class (qpu or sw) [qpu]:
Default solver [skip]:
Configuration saved.

跟着指示，选择是否同意，输入sapi，api token，就安装成功了。sapi 一般是 https://cloud.dwavesys.com/sapi/ ，token 在 Leap 页面就能找到。

基本概念

背包问题

假如给定一个背包，它的容量有限，只能装下总质量 \( W \) 的物品。现在有一些物品 \({x_i, i = 1\dots N}\)，对于每个物品 \(x_i\) 都有质量 \(w_i\) 以及价值 \(c_i\)。问怎样找到这样一组物品，使得它们可以被装入背包，同时总价值最大。

为了方便起见，我们设\( x_i \in [0, 1], i = 1,2, ... N \) 为变量。问题用数学建模表示为：

$$ Max~~~Cost(X) = \sum^{N}_{i = 1} c_i x_i \\ s.t. ~~~~~ Weight(X) = \sum^{N}_{i=1} w_i x_i \leq W $$

其中，\(Max\,\,Cost(x)\) 表示要最大化价值\(Cost(x)\)，\(s.t.~~Weight(X) \leq W\) 表示这个问题约束于质量限制。

为什么以0-1背包问题为例呢？因为“一切有界整型线性规划问题，都可以转化成等价的背包问题” （Transformation of integer programs to knapsack problems, Bradley, 197190005-7)）

问题模型

因为 D-Wave 做的是退火机，建模用的是 BQM，所以我们的问题必须建模成 Binary Quadratic Models (BQM)。从名字就可以得知，BQM 是：

• 布尔值变量（或者值域类似于{0, 1}只能取两个值的变量）
• 目标函数是 quadratic 的，也就是支持 \(x^2, x_i x_j, x\)

实际上还有个隐藏的条件是：

• 无约束

D-Wave 的库实际上还可以使用 Ising 模型和 QUBO（Quadratic Unconstrained Binary Optimization），都可以转化为 BQM，区别是表达方式。

$$ BQM: E(v) = \sum_{i=1}a_iv_i+\sum_{i<j} b_{i,j}v_iv_j+c, \\ where~~~~v_i \in \{−1,+1\}\ or\ \{0,1\} \\ QUBO: E(v) = \sum_{i\leq j} Q_{i,j}x_ix_j,~~~x_i \in \{0,1\} \\ Ising: E(v) = \sum_{i=1}h_is_i+\sum_{i<j}J_{i,j}s_is_j,~~~s_i \in \{−1,+1\}\\ $$

转化

为了在 Leap 上求解背包问题，还需要对不等式约束改写成等式、把整型变量改成布尔值变量。可能会有疑问，哪里来的整型变量？前文的建模一直在用 \(0/1\) 布尔二值变量啊？

不等式约束改写很简单，对于约束 \(W(X) \leq W\) 我们只需要使用一个变量 \(0 \leq y \leq W\)。这里上界自然是容纳的最大质量，注意下界并不是随便的一个\(0\)，在狭义的背包问题上，下界当然是什么都不装，自然是 \(0\)，而前面提到背包问题可以表示任何边界整型线性规划问题，下界自然也需要相应地考虑如何计算。

于是这就有了整型变量\(y\)，最简单的二进制编码方法是 one-hot 的： \(y_0, y_1, \dots y_W\)，当 \(y_k = 1\) 而 \(y_{\neq k} = 0\) 指示 \(y = k\)。哈密顿算符（Hamiltonian）之于退火机，如同损失函数之于机器学习。这里列出到目前为止，0-1背包问题的 Hamiltonian：

$$ H = H_A + H_B \\ H_A = A(1 - \sum^{W}_{k=1} y_k)^2 + A(\sum^{W}_{k=1} k y_k - \sum_{i} w_ix_i)^2 \\ H_B = -B\sum_{i} c_ix_i $$

其中，\(H_A\) 扮演了约束的角色。因为退火机寻找最低能量态，所以 \(H\) 越小说明优化得越好。\(H_A\) 的前一项指示 one-hot 约束，因为总质量每次只可能是一个值。第二项指示了总质量应该等于所有物品质量之和。最优时，\(H_A = 0\)， 约束被满足。同理 \(H_B\) 是总价值的负数，因为我们希望它越大越好。对于常量参数 \(A\) 与 \(B\) 的选择，应该满足：\(0 < B~max(c_i) < A\) （Lucas建议）。

显然 one-hot 编码在这个通用量子计算机还没有的世代是不现实的，因为它一共使用了 \(N + W\) 个量子位，当 \(W\) 非常大的时候是对 Qbit 的浪费。可以用自然的二进制编码，只需要将 \(k y_k\) 改写成 \(2^k y_k\)，同时去掉 \(y_i\) 只可以有一个取 \(y_i = 1\) 的约束.我们记 \(M = ⌊logW⌋\)，用于约束的变量便成了 \({y_i, i = 0\dots M-1}\)。对于一般的 \(2^M \leq W \leq 2^{M+1}\)，需要多一位\(y_{M} = W - (2^M - 1)\)。我们记 \(y_i\) 对应的值

$$ k_i = 2^i,~i \in [0, M-1] \\ k_{M} = W - (2^M - 1) \\ $$

最终：

$$ H_A = A(\sum_{i} k_i y_i - \sum_{i} w_ix_i)^2 \\ H_B = -B\sum_{i} c_ix_i $$

设 \(A = max(c_i), B=1\) 展开：

$$ H = A \ H'\\ H' = \sum_{i=1}^{N} (w_i^2 - c_i) x_i + 2\sum_{i=1}^{N}\sum_{j=i+1}^{N} w_i w_j x_i x_j \\ - 2\sum_{i=1}^{N}\sum_{j=0}^{M}w_i k_j x_i y_j \\ + \sum_{j=0}^{M}k_i^2 y_i + 2 \sum_{i=0}^{M}\sum_{j=i+1}^{M} k_i k_j y_i y_j $$

很容易得知 \(x^2 = x, y^2 = y\)，上式已作替换。

实现

# passed on 2020.12.11
from dwave.system import LeapHybridSampler
from math import log2, floor
import dimod

# 0-1 背包问题原型
W = 70
N = 7
c = [35, 85, 30, 50, 70, 80, 55]
w = [12, 27, 11, 17, 20, 10, 15]

# 其他变量
A = max(c)
M = floor(log2(W))
k = [2**i for i in range(M)] + [W + 1 - 2**M]

# BQM
bqm = dimod.AdjVectorBQM(dimod.Vartype.BINARY)

# x 项
for i in range(N):
    bqm.set_linear('x' + str(i), A * (w[i]**2) - c[i])

# x-x 项
for i in range(N):
    for j in range(i + 1, N):
        key = ('x' + str(i), 'x' + str(j))
        bqm.quadratic[key] = 2 * A * w[i] * w[j]

# x-y 项
for i in range(N):
    for j in range(M + 1):
        key = ('x' + str(i), 'y' + str(j))
        bqm.quadratic[key] = -2 * A * w[i] * k[j]

# y 项
for i in range(M + 1):
    bqm.set_linear('y' + str(i), A * (k[i]**2))

# y-y 项
for i in range(M + 1):
    for j in range(i + 1, M + 1):
        key = ('y' + str(i), 'y' + str(j))
        bqm.quadratic[key] = 2 * A * k[i] * k[j]

# 求解
sampler = LeapHybridSampler()
sampleset = sampler.sample(bqm)
sample = sampleset.first.sample
energy = sampleset.first.energy

# 哪些物品被选中了
selected = []
for varname, value in sample.items():
    if value and varname.startswith('x'): # x*
        selected.append(int(varname[1:]))
selected = sorted(selected)

weight_sum = 0
cost_sum = 0
for i in selected:
    weight_sum += w[i]
    cost_sum += c[i]

print('能量：', energy)
print('被选中的物品：', selected)
print('总质量：', weight_sum)
print('总价值：', cost_sum)

得到的结果：

能量： -270.0
被选中的物品： [0, 2, 4, 5, 6]
总质量： 68
总价值： 270

参考

量子信息学请参考 John Preskill CS219，基础概念可以读一下第一章讲义。

问题建模请参考神文 Ising formulations of many NP problems, Andrew Lucas, 2014