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


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

我们将会做到

67 声望
5 粉丝
0 条评论
推荐阅读
理解 D-Wave 量子退火机的计时
混合求解器(Hybrid solvers):量子-传统混合的求解器。量子处理器(Quantum Processing Unit, QPU):计算单元量子机器指令(Quantum Machine Instruction,QMI):包括 Ising 模型与 QUBO 的参数和退火周期参...

osinovsky阅读 1.4k

Ubuntu20.04 从源代码编译安装 python3.10
Ubuntu 22.04 Release DateUbuntu 22.04 Jammy Jellyfish is scheduled for release on April 21, 2022If you’re ready to use Ubuntu 22.04 Jammy Jellyfish, you can either upgrade your current Ubuntu syste...

ponponon1阅读 3.9k

日常Python 代码片段整理
1、简单的 HTTP Web 服务器 {代码...} 2、单行循环List {代码...} 3、更新字典 {代码...} 4、拆分多行字符串 {代码...} 5、跟踪列表中元素的频率 {代码...} 6、不使用 Pandas 读取 CSV 文件 {代码...} 7、将列表...

墨城2阅读 288

Unicode 正则表达式(qbit)
前言本文根据《精通正则表达式》和 Unicode Regular Expressions 整理。本文的示例默认以 Python3 为实现语言,用到 Python3 的 re 模块或 regex 库。基本的 Unicode 属性分类 {代码...} 基本的 Unicode 子属性Le...

qbit阅读 4.3k

Python + Sqlalchemy 对数据库的批量插入或更新(Upsert)
由于不同数据库对这种 upsert 的实现机制不同,Sqlalchemy 也就不再试图做一致性的封装了,而是提供了各自的方言 API,具体到 Mysql,就是给 insert statement ,增加了 on_duplicate_key_update 方法。

songofhawk1阅读 1.9k评论 4

封面图
Go for 循环有时候真的很坑。。。
大家好,我是煎鱼。不知道有多少 Go 的面试题和泄露,都和 for 循环有关。今天我在周末认真一看,发现了 redefining for loop variable semantics 。著名的硬核大佬 Russ Cox 表示他一直在研究这个问题,并表示十...

煎鱼阅读 3.4k

打脸了兄弟们,Go1.20 arena 来了!
大家好,我是煎鱼。大概半年前,我写过一篇文章《Go 要违背初心吗?新提案:手动管理内存》。有兴趣了深入解的同学,可以再回顾一下。当时我们还想着 Go 团队应该不会接纳,至少不会那么快:懒得翻也可以看我再次...

煎鱼阅读 3.1k

我们将会做到

67 声望
5 粉丝
宣传栏