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

osinovsky

写在前面

平台、配置、接口都不是一成不变的,不要拘泥于某一版本或者某一 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

阅读 1.3k

Osinovsky的学习笔记
一个计算机专业在读学生的笔记,内容可能是专业课以及其他阅读材料。

我们将会做到

67 声望
4 粉丝
0 条评论

我们将会做到

67 声望
4 粉丝
文章目录
宣传栏