相位估计

简介

相位估计 (Phase estimation) 是量子傅里叶变换的一个重要应用,它是很多量子算法的基本步骤,其中包括Shor算法HHL算法

问题描述: 已知给定可作用于量子线路的酉矩阵\(U\)以及其本征向量\(\vert\psi\rangle\), 求其对应的本征值。更具体的,给定酉矩阵\(U\)和量子态\(\vert\psi\rangle\)且满足\(U\vert\psi\rangle=e^{2\pi i \theta}\vert\psi\rangle\),其中\(\theta\in[0,1)\),要求估计\(\theta\)的值。相位估计算法可以以不小于\(1-\epsilon\)的概率正确估计\(\theta\) 的前\(n\)比特。

说明:因为\(U\)是酉矩阵,其特征值的模必须为1。因此它的特征值可以写成\(e^{2\pi i \theta},\theta\in[0,1)\)的形式。

在实现相位估计算法时,我们需要假设存在能够提供以下操作的 黑盒 (即oracle):

  1. 该黑盒可以制备初态\(\vert\psi\rangle\)
  2. 该黑盒可以执行受控\(U^{2^j}\)操作。

输入:相位估计算法的输入由两个寄存器组成,

  1. \(n\)个量子比特组成第一个寄存器。该寄存器用于存储数值\(\left\lfloor 2^n\theta \right\rfloor\)。如果我们将\(\theta\)二进制展开: \(\theta=0.\theta_0\theta_1\cdots\),那么有 \(\left\lfloor 2^n\theta \right\rfloor = \theta_0\theta_1\cdots\theta_{n-1}\),即\(n\)量子比特存放\(\theta\)的前\(n\)比特。\(n\)的选择需要从两方面来考虑:我们希望能得到的估计精度,以及我们希望估计算法的成功概率。简而言之,需要的精度和成功概率越高,那么\(n\)也越大。
  2. \(m\)个比特组成第二个寄存器,该寄存器用于存放状态\(\vert\psi\rangle\)。在算法的整个过程中,该寄存器保持不变。

量子电路模型:相位估计的量子线路模型描述如下图所示(该图截取自Wikipedia)。

phase-estimation

以下我们一步步分析该算法的流程。我们使用\(\vert\Psi_t\rangle\)来表示系统经过\(t\)步演化之后的状态。系统的初始状态可以表示为:

\[\vert\Psi_0\rangle = \vert0\rangle_n\otimes\vert\psi\rangle.\]

第一步:创建叠加态。正如几乎所有的量子算法初始阶段所做的那样,我们将系统的第一个寄存器置于最大叠加态,该步骤可以用\(n\)的Hadamard操作完成。第一步之后,系统状态变成

\[\vert\Psi_1\rangle = \frac{1}{\sqrt{2n}}\sum_{t=0}^{2^n-1}\vert t \rangle \vert \psi \rangle.\]

这里,我们假设\(t\)是二进制展开表示的。

第二步:执行受控酉操作。通过计算,我们可以验证以下等式

\[U^{2j}\vert\psi\rangle = U^{2j-1}e^{2\pi i\theta}\vert\psi\rangle = \cdots = e^{2\pi i 2^j\theta}\vert\psi\rangle.\]

\(t\)的二进制展开为\(t = t_0\cdots t_{n-1}\),其中\(t_i\in\{0,1\}\) ,那么有

\[\operatorname{C-U}\left(\vert t \rangle\vert\psi\rangle\right) = \vert t \rangle\left(\prod_{i=0}^{n-1} U^{2^{i}\times t_i}\right)\vert\psi\rangle = \vert t \rangle U^{t}\vert\psi\rangle = e^{2\pi i t\theta}\vert t\rangle\vert\psi\rangle.\]

第二步之后,系统状态变成

\[\vert\Psi_2\rangle = \frac{1}{\sqrt{2^n}}\sum_{t=0}^{2^n-1}e^{2\pi i t\theta}\vert t \rangle\vert \psi\rangle.\]

易见,经过受控酉门演化之后,\(\theta\)被映射到系统态的概率幅上。

第三步:执行逆傅里叶变换。通过仔细分析\(\vert\Psi_2\rangle\),我们可以发现它可以通过在状态\(\vert\theta\rangle\vert\psi\rangle\)上执行傅里叶变换得到(在一个寄存器上执行)。 这也就是说,我们可以通过逆傅里叶变换将\(\vert\Psi_2\rangle\)映射到\(\vert\theta\rangle\vert\psi\rangle\),从而将\(\theta\)映射到第一个寄存器中。第三步之后,系统状态变成

\[\vert\Psi_3\rangle = \frac{1}{\sqrt{2^n}} \sum_{t=0}^{2^n-1} e^{2\pi i t\theta} \left( \frac{1}{\sqrt{2^n}} \sum_{x=0}^{2^n-1} e^{\frac{-2\pi i t x}{2^n}} \vert x\rangle \right) \vert\psi\rangle = \frac{1}{2^n} \sum_{t=0}^{2^n-1}\sum_{x=0}^{2^n-1} e^{-\frac{2\pi it}{2^n}\left(x - 2^n\theta\right)}\vert x\rangle\vert\psi\rangle.\]

不妨令\(2^n\theta = a + 2^n\delta\),其中\(a\)是距离\(2^n\theta\)最近的整数,差值\(2^n\delta\)满足\(0 \leq 2^n\delta \leq 1/2\)。 使用这个展开式,\(\vert\Psi_3\rangle\) 可以表示为

\[\vert\Psi_3\rangle = \frac{1}{2^n} \sum_{t=0}^{2^n-1}\sum{x=0}^{2^n-1} e^{-\frac{2\pi it}{2^n}\left(x - a - 2^n\delta\right)}\vert x \rangle\vert\psi\rangle = \frac{1}{2^n} \sum_{t=0}^{2^n-1}\sum_{x=0}^{2^n-1} e^{-\frac{2\pi it}{2^n}\left(x - a\right)}e^{2\pi it \delta} \vert x\rangle\vert\psi\rangle.\]

第四步:测量得到估计值。 我们在第一个寄存器上执行标准基测量,得到输出\(\vert a\rangle\)的概率为

\[ \begin{align}\begin{aligned} \operatorname{Pr}\left(\vert a \rangle\right) = \left\vert \langle a \vert\Psi_3\rangle \right\vert^2 = \frac{1}{2^{2n}}\left\vert\sum_{t=0}^{2^n-1}e^{2\pi it \delta} \right\vert^2 = \frac{1}{2^{2n}}\left\vert \frac{1- e^{2\pi i 2^n\delta}}{1 - 2^{2\pi i\delta}} \right\vert^2\\以下我们分析这个概率值:\end{aligned}\end{align} \]
  • 第一种情况: \(\delta = 0\). 在这种情况下,我们的估计算法是精确算法,因为\(\theta = a/2^n\)。测量后系统的状态变成\(\vert a \rangle\otimes\vert \psi \rangle\)

  • 第二种情况: \(\delta \neq 0\). 在这种情况下,我们必须求测量得到\(\vert a\rangle\)的概率下界\(\operatorname{Pr}\left(\vert a \rangle\right)\)。 因为\(2^n\delta\leq 1/2\),所以我们的估计算法概率下界为

    \[\operatorname{Pr}\left(\vert a \rangle\right) \geq \frac{4}{\pi^2} \approx 0.405\]

    这个下界说明我们可以以大的正确概率估计\(\theta\)的前\(n\)比特。通过增加\(n\)的数目,比如说额外再增加\(O(\log(1/\epsilon))\)比特, 我们可以将成功概率提升至\(1-\epsilon\)

HiQ实现

简介部分可以看出,为了实现相位估计,HiQ必须有能力实现以下操作:

  1. 控制酉门。在HiQ中,该操作可以通过元指令Control实现。
  2. 酉门的幂次方。在HiQ中,该操作可以通过元指令Loop实现。
  3. 逆傅里叶变换。在HiQ中,该操作可以通过内置操作QFT加上函数get_inverse(QFT)得到。

首先,我们写一个函数用于计算我们以不小于\(1-\epsilon\)的概率正确估计\(\theta\)\(n\)比特所需要的比特数目。该数目的数学表达式及其推导可见Nielsen的书\(^{[2]}\),Eq. (5.35):

def _bits_required_to_achieve_accuracy(n, epsilon):
    tmp = math.log2(2 + 1/(2*epsilon))
    tmp = math.ceil(tmp)
    n_updated = int(n + tmp)
    return n_updated

然后,我们实现简介部分的Step1-3,作为相位估计的核心模块。之所以不把Step 4作为一部分加入到核心模块中,是因为很多时候我们并不需要测量相位估计之后量子态,而是作为一个输入提供给其他函数。

def run_phase_estimation(eng, U, state, m, n):
    """
    用于估计酉变换的核心模块代码
    参数:
        eng (MainEngine): 用于运行相位估计算法的主引擎
        U (np.matrix): 需要被估计的酉门,以其矩阵形式描述
        state (qureg): 相位估计初始输入态|Psi_0> = |psi> \otimes |0>。
                和简介部分描述保持一致,该状态分为两个部分:
                 1. 前m个比特用于存储需要被估计的特征状态 |psi>
                 2. 后n个比特初始化为|0>,用于存储估计的相位
        m (int): 前m个比特用于存储需要被估计的特征状态 |psi>
        n (int): 后n个比特初始化为|0>,用于存储估计的相位
    """
    # 后n个比特相对于state的偏移值
    OFFSET = m
    # 第一步&第二步:设置为叠加态并且执行受控U操作
    for k in range(n):
        H | state[OFFSET+k]

        with Control(eng, state[OFFSET+k]):
            num = int(math.pow(2, k))
            with Loop(eng, num):
                U | state[:OFFSET]

    # 第三步:逆QFT操作
    # 第四步:交换量子比特
    for k in range(math.floor(n/2)):
        Swap | (state[OFFSET+k], state[OFFSET+n-k-1])
    # 第五步:执行最初的GFT操作
    get_inverse(QFT) | state[OFFSET:]

参考文献:

[2] Nielsen, M., & Chuang, I. (2010). Quantum Computation and Quantum Information: 10th Anniversary Edition.