和 Leo 一起学量子算法: 二. Grover-Search 和量子推断 (Quantum Inference)

今天给大家讲一个量子计算和机器学习相联系的例子，量子推断，概念很强大，内容很朴实。

推断是机器学习中非常重要的概念，主要目的是获得条件概率 p(x|y)，意思是在发生 y 事件的前提下，x发生的概率，这在推荐系统，图片修复等领域有重要的应用。而真正把机器学习中的Bayesian推断和量子推断结合是在文章 Quantum inference on Bayesian networks 中，这个工作里面介绍的两种方案都是 Grover Search 的直系近亲。两者能联系起来的桥梁是量子力学波函数的概率论解释，或者说Born’s rule，基于该解释的概率模型在 @sym cheng 的文章 Information Perspective to Probabilistic Modeling: Boltzmann Machines versus Born Machines中被称为 Bron Machine。而今天我们的主要任务有两个

实现 Grover Search 基本算法
用量子推断还原一张破碎的图片

老司机可以直接看notebook

理论基础

Grover Search解决的问题是

假定有一个 oracle 线路存在，它选择性的翻转一个（或几个）构型的符号，求这些被翻转的构型？

理论推导可以在Neilson 的Quantum Computation and Quantum Information这本书里面P248找到。简而言之实现上面的示意图中的线路，基本思路是：

A 虚线框内的初始态制备
Grover迭代直到特定的次数

这个算法也有一些坑要注意，见文后附录。

代码实现

Grover Search 基本线路

1) Householder 线路

Grover Search的基本单元是Grover迭代, 也就是虚线框 C 内的操作，其精髓是构造Householder 算符 $R(|ψ⟩)=1-2|ψ⟩⟨ψ|$ , 即 B 虚线框内的内容。假设 $|ψ⟩=A|0⟩$ ，那么Householder 算符可以表达为 $R(|ψ⟩)=A(1-2|0⟩⟨0|)A^†$ 。

using Yao
using Yao.Blocks
using Compat
using Compat.Test
using StatsBase

"""
A way to construct oracle, e.g. inference_oracle([1,2,-3,5]) will
invert the sign when a qubit configuration matches: 1=>1, 2=>1, 3=>0, 5=>1.
"""
function inference_oracle(locs::Vector{Int})
    control(locs[1:end-1], abs(locs[end]) => (locs[end]>0 ? Z : chain(phase(π), Z)))
end

function reflectblock(A::MatrixBlock{N}) where N
    chain(N, A |> adjoint, inference_oracle(-collect(1:N)), A)
end

上述代码中，inference_oracle(-collect(1:N))实现了操作 $1-2|0⟩⟨0|$ ，即当且仅当比特处在全 0 的状态翻转相位。

一般在 Grover Search 中, 我们选择均匀叠加的态作为初态, 它可以用 Hadamard 门作用于每个 qubit 上面来构造, 它对应的虚线框 A 内的幺正操作. 同时，它对应于B 框内的 Householder 操作如下

nbit = 12
A = repeat(nbit, H)
ref = reflectblock(A)

打印出这个模块你会看到

julia> ref
Total: 12, DataType: Complex{Float64}
chain
├─ repeat on (1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12)
│  └─ H gate
├─ control(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11)
│  └─ (12,)=>chain
│     ├─ Global Phase Gate:3.141592653589793
│     └─ Z gate
└─ repeat on (1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12)
   └─ H gate

随后，不妨看看这样构造的 Householder 线路是否符合我们的预期～

@testset "test reflect" begin
    reg = rand_state(nbit)
    ref_vec = apply!(zero_state(nbit), A) |> statevec
    v0 = reg |> statevec
 @test -2*(ref_vec'*v0)*ref_vec + v0 ≈ apply!(copy(reg), ref) |> statevec
end

你应该会看到如下输出

Test Summary: | Pass  Total
test reflect  |    1      1
Base.Test.DefaultTestSet("test reflect", Any[], 1, false)

2). 构造 oracle 线路

有了 Householder 线路，就已经成功了一半啦！那么赶紧写个 oracle 测试下吧。以12个qubit为例，构造一个oracle 让波函数 100-105 位置上有翻转如下

# first, construct the oracle with desired state in the range 100-105.
oracle!(reg::DefaultRegister) = (reg.state[100:105,:]*=-1; reg)

# transform it into a function block, so it can be put inside a `Sequential`.
fb_oracle = FunctionBlock{:Oracle}(reg->oracle!(reg))

"""
ratio of components in a wavefunction that flip sign under oracle.
"""
function prob_match_oracle(psi::DefaultRegister, oracle)
    fliped_reg = apply!(register(ones(Complex128, 1<<nqubits(psi))), oracle)
    match_mask = fliped_reg |> statevec |> real .< 0
    norm(statevec(psi)[match_mask])^2
end

这里prob_match_oracle函数衡量了搜索到 oracle 中翻转位置的概率。这里还用到了FunctionBlock，可以把作用在Register上的 inplace 函数定义为一种 AbstractBlock从而可以被塞入量子线路。

3). Grover Search 整体线路

随后定义初态，并计算需要迭代的次数

# uniform state as initial state
psi0 = apply!(zero_state(nbit), A)

# the number of grover steps that can make it reach first maximum overlap.
num_grover_step(prob::Real) = Int(round(pi/4/sqrt(prob)))-1
niter = num_grover_step(prob_match_oracle(psi0, fb_oracle))

# construct the whole circuit
gb = sequence(sequence(fb_oracle, ref) for i = 1:niter)

这里的构造线路用了 sequence 函数, 它和 chain 不同在于, 它可以放非 MatrixBlock, 在此做出区分以方便不同方法的派发. 因为要迭代接近 20 次，所以这里不画出壮观的线路。为了能看到搜索的过程，我们分步骤运行线路如下

for (i, blk) in enumerate(gb)
    apply!(psi0, blk)
    overlap = prob_match_oracle(psi0, fb_oracle)
    println("step $i, overlap = $overlap")
end

运行结果如下

step 1, overlap = 0.01313214562833301
step 2, overlap = 0.03619369756224663
step 3, overlap = 0.07010978618384153
step 4, overlap = 0.11408666758254363
step 5, overlap = 0.16709514342697387
step 6, overlap = 0.22789464746611884
step 7, overlap = 0.29506227870937746
step 8, overlap = 0.3670261018170988
step 9, overlap = 0.44210193536698855
step 10, overlap = 0.518532767034414
step 11, overlap = 0.5945298732465308
step 12, overlap = 0.6683146809800716
step 13, overlap = 0.7381603920041205
step 14, overlap = 0.8024323954287373
step 15, overlap = 0.8596265227777784
step 16, overlap = 0.9084042502960318
step 17, overlap = 0.947624024645163
step 18, overlap = 0.9763679788679578
step 19, overlap = 0.9939634133826709
step 20, overlap = 0.9999985392841676

我们可以看到Grover Search算法仅通过20步就在维度为 4096 的Hilbert空间中找到了那几个翻转比特位，因为它的算法复杂度仅仅为 $O(\sqrt{N/M})$ ，其中 N 和 M 分别为全Hilbert空间大小和需要搜索的Hilbert空间的大小。这里要注意的是，如果继续迭代，这个overlap的数值反而会减小哦，尝试下你会发现它会呈现周期性。

这样的算法可以用于检索无序的数据库，更加酷的是，它还可以从联合概率分布 p(x, y) 中得到条件概率分布 p(x | y)，即量子推断。

量子推断

1) 数据准备

我们有三张 5 x 3 的图片

using Yao.Intrinsics

x1 = [0 1 0; 0 1 0; 0 1 0; 0 1 0; 0 1 0]
x2 = [1 1 1; 0 0 1; 1 1 1; 1 0 0; 1 1 1]
x0 = [1 1 1; 1 0 1; 1 0 1; 1 0 1; 1 1 1]

nbit = 15
v = zeros(1<<nbit)

# they occur with different probabilities.
for (x, p) in [(x0, 0.7), (x1, 0.29), (x2,0.01)]
    v[(x |> vec |> BitArray |> packbits)+1] = sqrt(p)
end

这里把三张图片以不同概率编码到了波函数 v 中，如果把这三个图片画出来，你会看到它们其实对应 0,1,2 三个数字

编码于量子波函数的三张图片, 代表 0, 1, 2 三个数字，几率分别为0.7，0.29和0.01

2) inference 线路构造

这里，搜索的oracle的含义变成了前提条件: 在1, 2, 3, 4, 5 比特上分别观测得到 10111. 改 oracle 等价于如下图所示的破损图片的补全问题, 或者说通过右侧已知条件, 推断蓝色区域的未知部分.

这是一张破坏的图片, 蓝色部分未知, 已知的是图中褐色和部分为1, 浅褐色的部分为 0.

这里，由于篇幅，这里不展示如何用幺正运算构造Householder线路（翻转线路），而是借用Yao.Blocks里面的基本模块ReflectBlock，它可以用reflect工厂函数构造。

构造 oracle 的时候要注意, 由于 Yao.jl 的小端在后, 所以之前 inference 给的前提条件中指定的1-5号 bit 等价于图中1-5红色标注的位置. 通过计算可以得到迭代次数 niter 应该为 7

rb = reflect(copy(v))
psi0 = register(v)

# we want to find the digits with the first 5 qubits [1, 0, 1, 1, 1].
evidense = [1, -2, 3, 4, 5]
oracle_infer = inference_oracle(evidense)(nbit)

niter = num_grover_step(prob_match_oracle(psi0, oracle_infer))
gb_infer = chain(nbit, chain(oracle_infer, rb) for i = 1:niter)

看下这个线路

julia> gb_infer
Total: 15, DataType: Complex{Float64}
chain
├─ chain
│  ├─ control(1, 2, 3, 4, 5)
│  │  └─ (6,)=>Z gate
│  └─ ReflectBlock(N = 15)
├─ chain
│  ├─ control(1, 2, 3, 4, 5)
│  │  └─ (6,)=>Z gate
│  └─ ReflectBlock(N = 15)
├─ chain
│  ├─ control(1, 2, 3, 4, 5)
│  │  └─ (6,)=>Z gate
│  └─ ReflectBlock(N = 15)
├─ chain
│  ├─ control(1, 2, 3, 4, 5)
│  │  └─ (6,)=>Z gate
│  └─ ReflectBlock(N = 15)
├─ chain
│  ├─ control(1, 2, 3, 4, 5)
│  │  └─ (6,)=>Z gate
│  └─ ReflectBlock(N = 15)
├─ chain
│  ├─ control(1, 2, 3, 4, 5)
│  │  └─ (6,)=>Z gate
│  └─ ReflectBlock(N = 15)
└─ chain
   ├─ control(1, 2, 3, 4, 5)
   │  └─ (6,)=>Z gate
   └─ ReflectBlock(N = 15)

还和之前一样迭代搜索

for (i, blk) in enumerate(gb_infer)
    apply!(psi0, blk)
    p_target = prob_match_oracle(psi0, oracle_infer)
    println("step $(i-1), overlap^2 = $p_target")
end

这段代码的输出为

step 1, overlap = 0.08761600000000003
step 2, overlap = 0.23055362560000003
step 3, overlap = 0.4161715569049601
step 4, overlap = 0.6150679135961745
step 5, overlap = 0.7957375127737549
step 6, overlap = 0.9295622899279725
step 7, overlap = 0.9953444003575993

我们又成功的找到了那个极大点!

3) 结果分析

最后我们可以通过下面的方式处理输出, 把它 reshape 成 5×3 的图片

pl = psi0 |> probs
config = findn(pl.>0.5)[] - 1 |> bitarray(nbit)
res = reshape(config, 5,3)

可以看到输出为

5×3 BitArray{2}:
  true   true   true
 false  false   true
  true   true   true
  true  false  false
  true   true   true

把它可视化后发现，它是2！

通过量子推断，我们得到1,2,3,4,5位置分别为1,0,1,1,1的图片为“2”

显然, 这是正确的推断! 这个技术有一个非常实际的应用场景, 那就是作为量子机器学习中, 量子推荐系统的一部分.

看到这里的同学恭喜你，你Get了最前沿的量子推断！是时候打开notebook看下啦～

附录

Grover Search的坑

通常的Grover迭代往往无法提前知道要迭代的次数, 好在不知道次数这个问题可以通过一些方式避免，比如 Rejection Sampling。
必须知道输入波函数的构造方式。

后记

如何用 Grover-Search 做 inference 其实是吉林大学优秀后辈吴宇峰最先给我看了 Bars & Stripes 这个数据库的 inference 实现后, 令我眼前一亮, 觉得可以写一个类似的 Tutorial, 所以他才是这个章节幕后的英雄 ~

来源：知乎 www.zhihu.com

作者：Leo

【知乎日报】千万用户的选择，做朋友圈里的新鲜事分享大牛。
点击下载