一点五. 和 Leo 一起学量子算法: 相位估计

理论基础

相位估计 (Phase estimation) 是量子傅里叶变换 QFT 的一个重要应用, 它的重要性体现在它是很多量子算法的基础, 其中包括但不限于

数论中的求阶问题
因式分解 (可用于破解 RSA 加密)
HHL 算法求解线性系统　(由于求解线性系统的重要价值，它最近10年内该领域最重要的理论突破之一)

相位估计问题定义为

给定可作用于量子线路的幺正矩阵 U 以及其本征向量之一 $|b\rangle$ , 求其对应的本征值.

由于幺正矩阵的本征值一定是幺模的，于是该本征值可以被表示为 $\lambda = e^{2\pi i\phi}$ 。求本征值在这里等价于求相位 $\phi$ ，从算法名可以看出，接下来的算法实际求解的是相位 $\phi$ 。

具体算法的推导在经典教材Neilson 的Quantum Computation and Quantum Information这本书里面P221已经讲得非常到位，也可以看Wiki，这里不再赘述。结论而言，要实现如下线路

示意图中比特Q1-5为存放结果的寄存器，Q6-8为输入向量。U矩阵作为受控门电路作用在Q6-8上。这里需要用到第一章中的介绍的QFT线路，心急的同学照样可以戳notebook。

理论虽然简单，但在量子线路中真正使用相位估计有一些“坑”要注意，具体见文末附录。

代码实现

定义问题

为了演示，首先来定义问题中的矩阵U和其本征向量 $|b\rangle$ ，可以采用谱分解来构造 $U = VSV^\dagger$ ，其中V为任意幺正矩阵。S必须是幺模的对角矩阵以保证幺正性，它的对角元对应于U的本征值。

第一步先用QR分解构造任意 8 x 8 幺正矩阵V（对应于 3 个qubit）

# Compatibility of Julia 0.6 and Julia 0.7
# 下面的几行是为了保证Julia的版本兼容性
# Pkg.add("Compat") <- for import error, run this (如果import出错，运行此行)
using Compat

using Yao

"""generate a random unitary matrix"""
rand_unitary(N::Int) = qr(randn(ComplexF64, N, N))[1]

N = 3
V = rand_unitary(1<<N)

这里任意指定为任意幺正矩阵V中第3个向量为 $|b\rangle$ ，并同时指定它的相位。下面代码中的signs即为 S 矩阵的对角元。

######### Phase Estimation Problem ##########
phases = rand(1<<N)
ϕ = Int(0b111101)/(1<<6)
phases[3] = ϕ    # this 3 is randomly chosen
signs = exp.(2pi*im.*phases)
U = V*Diagonal(signs)*V'
b = V[:,3]

那么目标就是让 U 和 b 输入量子线路，获得相位 ϕ。如果你也是使用Atom编辑器，可以在下面的REPL看下具体的角度。同时，给出经典解法，其时间复杂度为 $O(4^n)$ ，n 为b 矢量对应的比特数。

julia> ϕ
0.953125
# classical solution, notice log is a periodic function of 2π*im
julia> log(b'*U*b)/(2π*im) + 1
0.953125 - 1.546102970155876e-17im

量子线路解相位估计

定义两个寄存器与U操作

reg1初始化为0态，用于放置最后结果，易推算使用6个qubit结果最高精度是1/64。reg2用于存放矢量 b。工场函数matrixgate可定义任意矩阵门，用它来构造U矩阵对应的幺正门UG 。

M = 6
reg1 = zero_state(M)
reg2 = register(b)
UG = matrixgate(U)

2. 作用A虚线框内的Hadamard门

作用一个门到寄存器可以通过apply!函数，Julia里面会改变输入参数的函数按照惯例会在最后加!，这样用起来比较安全。

apply!(reg1, repeat(M, H))

工厂函数repeat定义了一个在多个比特上做同样操作的门，可以像下面这样

julia> repeat(5, H)
Total: 5, DataType: Complex{Float64}
repeat on (1, 2, 3, 4, 5)
└─ H gate


julia> repeat(5, H, (1,2,3))
Total: 5, DataType: Complex{Float64}
repeat on (1, 2, 3)
└─ H gate

3. 作用 B 虚线框内的控制线路

# construct a control circuit
reg = join(reg1, reg2)
control_circuit = chain(M+N)
for i = 1:M
    push!(control_circuit, control(M+N, (i,), (M+1:M+N...,)=>UG))
    if i != M
        UG = matrixgate(mat(UG) * mat(UG))
    end
end
apply!(reg, control_circuit)

首先，用join函数把reg1, reg2组合成一个新的大寄存器，控制线路是个ChainBlock（API仿照了Vector，所以支持了push!，索引等操作），其组成单元ControlBlock的受控比特位置是(M+1:M+N...,)，等价于tuple (M+1, M+2, ..., M+N)。让我们在REPL里面看下这个线路

julia> control_circuit
Total: 9, DataType: Complex{Float64}
chain
├─ control(1)
│  └─ (7, 8, 9)=>GeneralMatrixGate(2^3 × 2^3; Array{Complex{Float64},2})
├─ control(2)
│  └─ (7, 8, 9)=>GeneralMatrixGate(2^3 × 2^3; Array{Complex{Float64},2})
├─ control(3)
│  └─ (7, 8, 9)=>GeneralMatrixGate(2^3 × 2^3; Array{Complex{Float64},2})
├─ control(4)
│  └─ (7, 8, 9)=>GeneralMatrixGate(2^3 × 2^3; Array{Complex{Float64},2})
├─ control(5)
│  └─ (7, 8, 9)=>GeneralMatrixGate(2^3 × 2^3; Array{Complex{Float64},2})
└─ control(6)
   └─ (7, 8, 9)=>GeneralMatrixGate(2^3 × 2^3; Array{Complex{Float64},2})

4. 把第一章中定义的IQFT线路作用在前M个qubit上

focus!(reg, 1:M)
apply!(reg, QFT(M) |> adjoint)

这里focus!函数可以修改寄存器的active qubit ，active qubit 意思是可操作的比特，这样随后作用的门就只看到第1 – M个qubit，以保证和QFT门的大小一致。其实，focus!不仅可以分割“系统”和“环境”，有时候还可以改变量子比特的编号。

julia> focus!(rand_state(9), [1, 3, 9, 5, 6, 2])
DefaultRegister{1, Complex{Float64}}
    active qubits: 6/9

5. 测量结果

res = measure(reg, 10)

measure函数第二个参数是测量次数，这里的测量不是物理测量，它不改变寄存器。物理测是measure!，因为它有副作用，会改变输入的寄存器，对应量子态塌缩，从而没法选择测量次数。

julia> res
10-element Array{Int64,1}:
 47
 47
 47
 47
 47
 47
 47
 47
 47
 47
julia> bin(47, M)
"101111"

看下结果，一致，都是47，如何解读呢？不如把它变成2进制看看，发现结果还是不对，眼尖的同学发现结果大小端反了呀（还记得上次说的QFT作用完会调换大小端吗？）。不要慌，我们的Yao.Intrinsics里面定义了很多二进制操作，其中breflect可以胜任这个翻转比特的需求。

using Yao.Intrinsics
res = breflect(M, res[1])/(1<<M)
println("Phase is 2π * $(res), the exact value is 2π * $ϕ")

最后你会看到输出结果

Phase is 2π * 0.953125, the exact value is 2π * 0.953125

看到这行输出的同学，恭喜你，Get了相位估计！你的下一步计划是什么呢？用Shor算法去破译金库密码？还是用HHL做量子机器学习？

那么最后，肯定有人疑惑为什么不像QFT那样把相位估计做成一个MatrixBlock呢？答案是做不到，原因见附录。

第一章（包括这个一点五）的全部内容可以在这个notebook找到，Enjoy the coding ～

附录

相位估计算法有哪些坑 (Caveats)？

“坑” 是一个量子算法的重要组成部分，是一个量子算法所依赖的重要前提假设。

如果涉及的算法用到了相位估计却没有好好的考虑这些坑, 那就有可能变成纸上谈兵. 对于相位估计, 它主要有以下两个坑

对于整数 k, 必须要有一种方法可以高效的计算量子门 $C-U^{2^k}$
保证初态是 U 的本征态.

在求阶，因式分解和HHL算法的文献里，都可以看到对这些约束的考虑。

为何相位估计不是MatrixBlock？

如果相位估计是理想的，那么它并不会改变reg2里面的态，于是这个操作对应于单位矩阵。

然而现实的往往更加复杂，下面的情况可能导致它影响输入态

输入态不严格处于U的本征态
辅助比特过少而不能完美表达 ϕ

归根结底还是因为有了测量的存在，无法简单的把输出简单的定义为矩阵操作的结果。

虽然它不是不是MatrixBlock，Yao.jl任然有办法把它集成到线路框架中去，相位估计可以作为函数被FunctionBlock打包加入到其它线路，此时它的属性就是AbstractBlock，是比MatrixBlock更低阶的公民。

后记

QFT作为常用算法，已经被我们加入了标准库的Yao.Zoo里面。

可以用如下方式调用

using Yao.Zoo
qft = QFTBlock{5}()

其效果虽然等价于原来的QFT，但由于模拟的时候可以当“黑箱”，里面用了经典FFT加相应的变换来得到输出来保证高效，即保证正确性的同时并不是老老实实的计算了那么多门。

由于 Yao.jl 的目标之一是白箱，你可以通过如下方式得到一个黑箱线路的白箱版本：

white_qft = openbox(qft)

这样，你就可以得到一个忠实的QFT线路，并且看到里面的内容啦～

不妨试试用QFTBlock替代之前写的QFT吧！

来源：知乎 www.zhihu.com

作者：Leo

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