量子相位估计（QPE）

量子相位估计在很多量子算法中都发挥着重要作用，而且是 Shor 算法的重要一环，需要解决的问题大致如下：

给定一个酉矩阵 $U$ 和它的特征向量 $|u\rangle$，我们想要找到相应的特征值，$U|u\rangle = e^{2\pi i\varphi} |u\rangle$ ( 相位 $\varphi$ 未知 )。由于 $U$ 是一个酉矩阵， 我们可以保证 $0\leq\varphi<1,.$

算法部分

为了阐明相位估计算法的运行原理，考虑一个例子。假设 $$ U = \begin{bmatrix} 1 &0 \\ 0 & e^{2i\pi/5} \end{bmatrix},, $$ 它其中一个特征向量是 $|u\rangle = |1\rangle$，因此我们想要估计的相位是 $\varphi = 1/5$，$U|1\rangle = e^{2\pi i\varphi} |1\rangle$。一般来说，QPE 可以分为四个部分：

• 预处理

  QPE 需要两个量子存储器，第一个量子存储器被称为计数存储器，包含 $n$ 个量子比特且全部初始化为 $|0\rangle$。我们将会以 2 进制的形式估计未知的相位 $\varphi$。$n$ 的大小决定了最后获得结果的精度。第二个量子存储器含有 $m$ 个量子比特，需要被初始化为特征向量 $|u\rangle$。在我们的例子中 $m=1，n=3$。由于第一个计数量子存储器只有三个量子比特，精度不会太高。这两个量子存储器的初始态应该是 $$ |\psi_0\rangle = |0\rangle^{\otimes n} \otimes |u\rangle =|000\rangle \otimes |1\rangle,. $$

  我们可以在模拟器内轻松初始化：

  qubit_num = 4
  shots = 1000
  phase = 2 * np.pi / 5
  
  
  def main():
    """
    main
    """
    # Create environment
    env = QEnv()
    # Choose backend Baidu Local Quantum Simulator-Sim2
    env.backend(BackendName.LocalBaiduSim2)
  
    # Initialize qubits
    q = env.Q.createList(qubit_num)
  
    # Prepare eigenstate |1> = X|0> on the ancilla qubit
    X(q[3])

• 制造叠加态 $H$

  接下来我们在第一个储存器的所有量子比特上施加阿达玛门。由于它们的初始值为 $|0\rangle$，因此制作出一个叠加态后，我们得到的是 $$ |\psi_1\rangle = (H|0\rangle)^{\otimes n} \otimes |u\rangle = \frac{1}{\sqrt{2^n}}(|0\rangle + |1\rangle)^{\otimes n} \otimes |u\rangle $$

  在我们的例子中，$n=3$ 并且 $|u\rangle = |1\rangle$。在模拟器上，你可以按如下的方法使用阿达玛门

      # Superposition
      H(q[0])
      H(q[1])
      H(q[2])

• 量子黑箱 $\hat{O}$

  QPE 中的量子黑箱是一些受控酉矩阵，计数储存器的第 $j$ 个量子比特是控制比特，而受控酉矩阵 $U^{2^j}$ 则是作用在第二个量子储存器整体。令 $j = 2$，由于第二个储存器始终保持在特征能量态 $|u\rangle$ 上，我们得到 $$ U^{4}|u\rangle = U^{3}U|u\rangle = e^{2\pi i\varphi} U^{3}|u\rangle = e^{2\pi i(2^1\varphi)}U^{2}|u\rangle = \cdots = e^{2\pi i (2^2\varphi)} |u\rangle $$

  可以看出这里有明显的递归关系， $$ U^{2^j}|u\rangle = U^{2^j-1}U|u\rangle = U^{2^j-1}e^{2\pi i\varphi} |u\rangle = \cdots = e^{2 \pi i (2^j \varphi)} |u\rangle, \quad j\in [0, n-1] $$

  然后我们需要搞明白这些受控酉矩阵到底都做了些什么。考虑第一储存器中的一般量子态 $$ |\phi\rangle = \alpha|0\rangle + \beta|1\rangle $$ 施加受控 U 门在整个量子系统上，可以得到 $$ CU|\phi\rangle \otimes |u\rangle = \alpha|0\rangle \otimes |u\rangle + \beta \cdot U|1\rangle \otimes |u\rangle = \alpha|0\rangle \otimes |u\rangle + \beta \cdot e^{2 \pi i\varphi}|1\rangle \otimes |u\rangle $$

  这个技巧通常被称为相位反冲 (phase kickback)。注意到它的作用就是，通过 U 在第二储存器的作用，给控制比特增加了一个未知的相对相位 $\varphi$。绝大多数情况下，我们的直觉是控制比特应该不受影响的这里就比较特殊。随后我们在量子态 $|\psi_1\rangle$ 上施加 $\hat{O}$ (Control-$U^{2^j}$) $$ \hat{O} |\psi_1\rangle = \frac{1}{\sqrt{2^n}}(|0\rangle + e^{2\pi i (2^{n-1}\varphi)}|1\rangle) \otimes\cdots \otimes (|0\rangle + e^{2\pi i (2^{0}\varphi)}|1\rangle) \otimes |u\rangle $$ 在我们 3 量子比特的例子中， $$ |\psi_2\rangle = \hat{O} |\psi_1\rangle =\frac{1}{\sqrt{8}}(|0\rangle + e^{2\pi i (4\varphi)}|1\rangle) \otimes(|0\rangle + e^{2\pi i (2\varphi)}|1\rangle) \otimes (|0\rangle + e^{2\pi i (1\varphi)}|1\rangle) \otimes |1\rangle $$ 目前为止，我们构造的量子电路是

  

  好像还是有些不清楚，这些 $U^{2^j}$ 到底是什么东西？其实就是重复施加 $2^j$ 次受控 U 门

  

  在我们的模拟器中，受控 U 门可以由一个三参数受控旋转门生成。其中广义旋转门的定义是： $$ R(\theta, \phi, \varphi) = \begin{bmatrix} \cos(\frac{\theta}{2}) & -e^{i\varphi}\sin(\frac{\theta}{2})\\ e^{i\phi}\sin(\frac{\theta}{2}) & e^{i(\phi+\varphi)} \cos(\frac{\theta}{2}) \end{bmatrix} $$

  当这三个参数取一些特定值的时候，我们就可以得到想要的酉矩阵。这里可以令 $\varphi' = 2\pi/5$。 $$ R(0, 0, \varphi') = \begin{bmatrix} 1 &0\\ 0 & e^{i\varphi'} \end{bmatrix} $$ 代码如下

      # Control-U gates
      CU(0, 0, phase)(q[0], q[3])
  
      CU(0, 0, phase)(q[1], q[3])
      CU(0, 0, phase)(q[1], q[3])
  
      CU(0, 0, phase)(q[2], q[3])
      CU(0, 0, phase)(q[2], q[3])
      CU(0, 0, phase)(q[2], q[3])
      CU(0, 0, phase)(q[2], q[3])

• 量子傅里叶变换 (QFT)

  我们已经很接近最终答案了。不过，这是在傅里叶基向量的意义下。我们需要回归到经典的计算基去做测量，然后就可以读出 3 比特的近似相位值 $\hat{\varphi}$。关于量子傅里叶变换的详细讲解会单独开一个文档。介绍 QFT 之前，我们需要引入记号 $\varphi = 0.\varphi_0\cdots \varphi_{n-1}$ 来代表二进制分数 $\varphi = \varphi_0/2+\cdots + \varphi_{n-1}/2^n$，其中每一个 $\varphi_i \in {0,1}$。在我们的例子中，由于只选取了 $n=3$ 个量子比特来估计相位，最好的精度只到达到以下两者： $$ \varphi_0 =0, \varphi_1 =0, \varphi_2 = 1 \quad \Rightarrow \quad \hat{\varphi}{-} = \frac{1}{2}\varphi_0 + \frac{1}{4}\varphi_1 + \frac{1}{8}\varphi_2 = \frac{1}{8} $$ 还有 $$ \varphi_0 =0, \varphi_1 = 1, \varphi_2 = 0 \quad \Rightarrow \quad \hat{\varphi}{+} = \frac{1}{2}\varphi_0 + \frac{1}{4}\varphi_1 + \frac{1}{8}\varphi_2 = \frac{1}{4} $$ 我们期望测量后概率振幅集中在这个范围内。如果我们想增加精度，或者说减小这个范围，那增加第一储存器中的量子比特数量 $n$ 即可。回到之前得到的 $|\psi_2 \rangle$ 可以用引入的新记号表示为： $$ |\psi_2\rangle =\frac{1}{\sqrt{8}}(|0\rangle + e^{2\pi i (0.\varphi_2)}|1\rangle) \otimes(|0\rangle + e^{2\pi i (0.\varphi_1\varphi_2)}|1\rangle) \otimes (|0\rangle + e^{2\pi i (0.\varphi_0\varphi_1\varphi_2)}|1\rangle) \otimes |1\rangle $$ 为什么能这么做呢？不妨考察， $$ 2^1 \varphi = 2^1* (0.\varphi_0\varphi_1\varphi_2) = 2*( \frac{1}{2}\varphi_0 + \frac{1}{4}\varphi_1 + \frac{1}{8}\varphi_2) = \varphi_0 + \frac{1}{2}\varphi_1 + \frac{1}{4}\varphi_2 = 0.\varphi_1\varphi_2 $$ 同样的， $$ 2^2 \varphi = 2^2* (0.\varphi_0\varphi_1\varphi_2) = 4*( \frac{1}{2}\varphi_0 + \frac{1}{4}\varphi_1 + \frac{1}{8}\varphi_2) = 2\varphi_0 + \varphi_1 + \frac{1}{2}\varphi_2 = 0.\varphi_2 $$ 推广可以得到 $$ 2^j \varphi = 2^j* (0.\varphi_0\cdots\varphi_{n-1}) = 2^j*( \frac{1}{2}\varphi_0 + \cdots + \frac{1}{2^n}\varphi_{n-1}) = 2^{j-1}\varphi_0 + \cdots+ \frac{1}{2}\varphi_j + \cdots = 0.\varphi_j\cdots\varphi_{n-1} $$ 现在我们可以去理解量子傅里叶变换能做些什么了，它把一个计算基中的元素映射到傅里叶基向量 $$ QFT_n |\varphi_0, \cdots, \varphi_{n-1}\rangle = \frac{1}{\sqrt{2^n}} (|0\rangle + e^{2\pi i (0.\varphi_{n-1})}|1\rangle) \otimes \cdots \otimes (|0\rangle + e^{2\pi i (0.\varphi_0\cdots\varphi_{n-1})}|1\rangle) $$ 这就是我们之前看到的 $|\psi_2\rangle$ 的形式。既然我们需要回到计算基中去测量，不如在 $|\psi_2\rangle$ 施加 3 量子比特傅里叶逆变换。 $$ QFT^{-1}3\otimes I |\psi_2 \rangle = |\varphi_0, \varphi_1, \varphi_2\rangle \otimes |u\rangle $$ 这个变换的矩阵表示是 $$ QFT^{-1}{3} = \frac{1}{\sqrt{8}} \begin{bmatrix} 1 &1 &1 &1 &1 &1 &1 &1\\ 1 &\sqrt{-i} &-i &-\sqrt{i} &-1 &-\sqrt{-i} &i &\sqrt{i}\\ 1 &-i &-1 &i &1 &-i &-1 &i\\ 1 &-\sqrt{i} &i &\sqrt{-i} &-1 &i &-i &-\sqrt{-i}\\ 1 &-1 &1 &-1 &1 &-1 &1 &-1\\ 1 &-\sqrt{-i} &i &\sqrt{i} &-1 &\sqrt{-i} &-i &-\sqrt{i}\\ 1 &i &-1 &-i &1 &i &-1 &-i\\ 1 &\sqrt{i} &i &-\sqrt{-i} &-1 &-\sqrt{i} &-i &\sqrt{-i} \end{bmatrix} $$ 电路表示是

  

  注意到如下两者量子门的等价性， $$ R(0, 0, -\pi/4) = \begin{bmatrix} 1 &0\\ 0 & e^{-i\frac{\pi}{4}} \end{bmatrix} =T^\dagger,. $$ 还有 $$ R(0, 0, -\pi/2) = \begin{bmatrix} 1 &0\\ 0 & e^{-i\frac{\pi}{2}} \end{bmatrix} = \begin{bmatrix} 1 &0\\ 0 & -i \end{bmatrix} =S^\dagger,. $$ 有了这些之后，傅里叶逆变换就可以按如下方法在模拟器中制得

      # 3-qubit inverse QFT
      SWAP(q[0], q[2])
      H(q[0])
      CU(0, 0, -np.pi / 2)(q[0], q[1])
      H(q[1])
      CU(0, 0, -np.pi / 4)(q[0], q[2])
      CU(0, 0, -np.pi / 2)(q[1], q[2])
      H(q[2])

  • 测量

    费尽千辛万苦，我们终于可以读取出精心获得的成果 $\hat{\varphi} = 0.\varphi_0\varphi_1\varphi_2$ 结合上面所有步骤，我们得到测量结果:

    

    我们可以看到概率最高的态是 $|\varphi_0, \varphi_1, \varphi_2\rangle = |010\rangle$，正好对应着相对相位的估计值 $$ \hat{\varphi}_{3-qubit} = 0.\varphi_0\varphi_1\varphi_2 = \frac{1}{2}\varphi_0 + \frac{1}{4}\varphi_1 + \frac{1}{8}\varphi_2 = 0.25 $$

    最终答案有 $25%$ 的误差。但这是三个量子比特的计数寄存器能得到的最好结果了。如果我们用一个四量子比特的储存器，那么精度会更高 (误差按照 $1/2^n$ 的规模减小)。

代码示范

下面是完整代码，你可以修改数值做一些不同的尝试

import numpy as np

import sys
sys.path.append('../../..')  # "from QCompute import *" requires this
from QCompute import *

matchSdkVersion('Python 3.3.3')

qubit_num = 4
shots = 1000
phase = 2 * np.pi / 5


def main():
    """
    main
    """
    # Create environment
    env = QEnv()
    # Choose backend Baidu Local Quantum Simulator-Sim2
    env.backend(BackendName.LocalBaiduSim2)

    # Initialize qubits
    q = env.Q.createList(qubit_num)

    # Prepare eigenstate |1> = X|0> on the ancilla qubit
    X(q[3])

    # Superposition
    H(q[0])
    H(q[1])
    H(q[2])

    # Control-U gates
    CU(0, 0, phase)(q[0], q[3])

    CU(0, 0, phase)(q[1], q[3])
    CU(0, 0, phase)(q[1], q[3])

    CU(0, 0, phase)(q[2], q[3])
    CU(0, 0, phase)(q[2], q[3])
    CU(0, 0, phase)(q[2], q[3])
    CU(0, 0, phase)(q[2], q[3])

    # 3-qubit inverse QFT
    SWAP(q[0], q[2])
    H(q[0])
    CU(0, 0, -np.pi / 2)(q[0], q[1])
    H(q[1])
    CU(0, 0, -np.pi / 4)(q[0], q[2])
    CU(0, 0, -np.pi / 2)(q[1], q[2])
    H(q[2])

    # Measurement result
    MeasureZ(*env.Q.toListPair())
    taskResult = env.commit(shots, fetchMeasure=True)
    return taskResult['counts']


if __name__ == '__main__':
    main()

讨论

可以看到，相位估计算法就是为了寻找一个给定酉矩阵的特征值的。这就是它为什么有时候会被称为特征值估计。你可以试着更改相位的数值为 $\varphi = 2\pi/3$，然后重新跑一下程序。看看你自己能不能解释得到的结果？

参考文献

Nielsen, Michael A. & Isaac L. Chuang (2001). Quantum computation and quantum information (Repr. ed.). Cambridge [u.a.]: Cambridge Univ. Press. ISBN 978-0521635035.