Qiskit-汉化3.7 Shor's Algorithm

Dora Dōjo
1

Shor算法

Shor算法因在多项式时间内分解整数而闻名。由于最著名的经典算法需要超多项式时间来分解两个素数的乘积，因此广泛使用的密码系统RSA依赖于因子分解对于足够大的整数是不可能的。

本章我们将关注Shor算法的量子部分，它实际上解决了寻找周期的问题。由于因式分解问题可以转化为多项式时间内的周期发现问题，因此也可以使用一个高效的周期发现算法来高效地分解整数。现在我们可以证明，如果我们能够高效地计算 a^x\bmod N 的周期，那么我们也可以高效地因式分解。由于周期查找本身就是一个有价值的问题，我们将首先解决这个问题，然后在第5节讨论如何将其用于因子分解。

import matplotlib.pyplot as plt
import numpy as np
from qiskit import QuantumCircuit, Aer, transpile, assemble
from qiskit.visualization import plot_histogram
from math import gcd
from numpy.random import randint
import pandas as pd
from fractions import Fraction
print("Imports Successful")

1. 问题在于:周期查找

我们来看一下周期函数:

f(x) = a^x \bmod{N}
提醒：模和模运算（点击这里展开）

求模运算(缩写为mod)的意思是求一个数除以另一个数的余数。例如:

17 \bmod 5 = 2

因为 17 \div 5 = 3 ，余数 2 。(即 17 =(3\times 5)+ 2 )。在Python中，求模操作用%符号表示。这种行为在模运算中使用，其中数字在达到某个值(模)后“换行”。使用模运算，我们可以这样写:

17 = 2 \pmod 5

注意，这里 \pmod 5 适用于整个方程(因为它在括号中)，而不像上面的方程，它只适用于方程的左侧。

其中 a 和 N 是正整数， a 小于 N ，并且它们没有公约数。周期或顺序( r )是最小的(非零的)整数，满足以下条件:

a^r \bmod N = 1

我们可以在下面的图中看到这个函数的一个例子。请注意，点与点之间的线是为了帮助看到周期性，而不是表示x标记之间的中间值。

1

2.解决方案

Shor的解决方案是在酉算子上使用量子相位估计:

U|y\rangle \equiv |ay \bmod N \rangle

为了看看这有什么用，我们来算出U的特征态是什么样的。如果我们从这个状态开始 |1\rangle ，我们可以看到，U的每个连续应用将使寄存器的状态乘以 a \pmod N ，在 r 应用之后，我们将到达状态 |1\rangle 。例如， a = 3 和 N = 35 :

\begin{aligned} U|1\rangle &= |3\rangle & \\ U^2|1\rangle &= |9\rangle \\ U^3|1\rangle &= |27\rangle \\ & \vdots \\ U^{(r-1)}|1\rangle &= |12\rangle \\ U^r|1\rangle &= |1\rangle \end{aligned}

2

所以这个循环的叠加态( |u_0\rangle )就是 U 的特征态:

|u_0\rangle = \tfrac{1}{\sqrt{r}}\sum_{k=0}^{r-1}{|a^k \bmod N\rangle}
点击扩展：样例 $a=3, N=35$
\begin{aligned} |u_0\rangle &= \tfrac{1}{\sqrt{12}}(|1\rangle + |3\rangle + |9\rangle \dots + |4\rangle + |12\rangle) \\[10pt] U|u_0\rangle &= \tfrac{1}{\sqrt{12}}(U|1\rangle + U|3\rangle + U|9\rangle \dots + U|4\rangle + U|12\rangle) \\[10pt] &= \tfrac{1}{\sqrt{12}}(|3\rangle + |9\rangle + |27\rangle \dots + |12\rangle + |1\rangle) \\[10pt] &= |u_0\rangle \end{aligned}

这个特征态的特征值是1，这没什么意思。一个更有趣的本征态可能是每个计算基态的相位不同。具体来说，让我们看看 k -th状态的相位与 k 成比例的情况:

\begin{aligned} |u_1\rangle &= \tfrac{1}{\sqrt{r}}\sum_{k=0}^{r-1}{e^{-\tfrac{2\pi i k}{r}}|a^k \bmod N\rangle}\\[10pt] U|u_1\rangle &= e^{\tfrac{2\pi i}{r}}|u_1\rangle \end{aligned}
点击扩展：样例 $a=3, N=35$
\begin{aligned} |u_1\rangle &= \tfrac{1}{\sqrt{12}}(|1\rangle + e^{-\tfrac{2\pi i}{12}}|3\rangle + e^{-\tfrac{4\pi i}{12}}|9\rangle \dots + e^{-\tfrac{20\pi i}{12}}|4\rangle + e^{-\tfrac{22\pi i}{12}}|12\rangle) \\[10pt] U|u_1\rangle &= \tfrac{1}{\sqrt{12}}(|3\rangle + e^{-\tfrac{2\pi i}{12}}|9\rangle + e^{-\tfrac{4\pi i}{12}}|27\rangle \dots + e^{-\tfrac{20\pi i}{12}}|12\rangle + e^{-\tfrac{22\pi i}{12}}|1\rangle) \\[10pt] U|u_1\rangle &= e^{\tfrac{2\pi i}{12}}\cdot\tfrac{1}{\sqrt{12}}(e^{\tfrac{-2\pi i}{12}}|3\rangle + e^{-\tfrac{4\pi i}{12}}|9\rangle + e^{-\tfrac{6\pi i}{12}}|27\rangle \dots + e^{-\tfrac{22\pi i}{12}}|12\rangle + e^{-\tfrac{24\pi i}{12}}|1\rangle) \\[10pt] U|u_1\rangle &= e^{\tfrac{2\pi i}{12}}|u_1\rangle \end{aligned}

我们可以看到 r=12 出现在相位的分母中。

这是一个特别有趣的特征值，因为它包含 r 。事实上，必须包含 r 以确保 r 计算基状态之间的相位差是相等的。这不是这种行为的唯一特征态;为了进一步推广，我们可以将一个整数 s 乘以这个相位差，它将显示在我们的特征值中:

\begin{aligned} |u_s\rangle &= \tfrac{1}{\sqrt{r}}\sum_{k=0}^{r-1}{e^{-\tfrac{2\pi i s k}{r}}|a^k \bmod N\rangle}\\[10pt] U|u_s\rangle &= e^{\tfrac{2\pi i s}{r}}|u_s\rangle \end{aligned}
点击扩展：样例 a=3 , N=35
\begin{aligned} |u_s\rangle &= \tfrac{1}{\sqrt{12}}(|1\rangle + e^{-\tfrac{2\pi i s}{12}}|3\rangle + e^{-\tfrac{4\pi i s}{12}}|9\rangle \dots + e^{-\tfrac{20\pi i s}{12}}|4\rangle + e^{-\tfrac{22\pi i s}{12}}|12\rangle) \\[10pt] U|u_s\rangle &= \tfrac{1}{\sqrt{12}}(|3\rangle + e^{-\tfrac{2\pi i s}{12}}|9\rangle + e^{-\tfrac{4\pi i s}{12}}|27\rangle \dots + e^{-\tfrac{20\pi i s}{12}}|12\rangle + e^{-\tfrac{22\pi i s}{12}}|1\rangle) \\[10pt] U|u_s\rangle &= e^{\tfrac{2\pi i s}{12}}\cdot\tfrac{1}{\sqrt{12}}(e^{-\tfrac{2\pi i s}{12}}|3\rangle + e^{-\tfrac{4\pi i s}{12}}|9\rangle + e^{-\tfrac{6\pi i s}{12}}|27\rangle \dots + e^{-\tfrac{22\pi i s}{12}}|12\rangle + e^{-\tfrac{24\pi i s}{12}}|1\rangle) \\[10pt] U|u_s\rangle &= e^{\tfrac{2\pi i s}{12}}|u_s\rangle \end{aligned}

对于 s 的每个整数值，我们现在有一个唯一的特征态，其中

0 \leq s \leq r-1.

很方便，如果我们把所有的本征态加起来，不同的相位抵消了所有的计算基态，除了 |1\rangle

\tfrac{1}{\sqrt{r}}\sum_{s=0}^{r-1} |u_s\rangle = |1\rangle
点击扩展：样例 a=7, N=15

为此，我们将看一个较小的示例 a = 7 and N=15 . In this case r=4 :

\begin{aligned} \tfrac{1}{2}(\quad|u_0\rangle &= \tfrac{1}{2}(|1\rangle \hphantom{e^{-\tfrac{2\pi i}{12}}}+ |7\rangle \hphantom{e^{-\tfrac{12\pi i}{12}}} + |4\rangle \hphantom{e^{-\tfrac{12\pi i}{12}}} + |13\rangle)\dots \\[10pt] + |u_1\rangle &= \tfrac{1}{2}(|1\rangle + e^{-\tfrac{2\pi i}{4}}|7\rangle + e^{-\tfrac{\hphantom{1}4\pi i}{4}}|4\rangle + e^{-\tfrac{\hphantom{1}6\pi i}{4}}|13\rangle)\dots \\[10pt] + |u_2\rangle &= \tfrac{1}{2}(|1\rangle + e^{-\tfrac{4\pi i}{4}}|7\rangle + e^{-\tfrac{\hphantom{1}8\pi i}{4}}|4\rangle + e^{-\tfrac{12\pi i}{4}}|13\rangle)\dots \\[10pt] + |u_3\rangle &= \tfrac{1}{2}(|1\rangle + e^{-\tfrac{6\pi i}{4}}|7\rangle + e^{-\tfrac{12\pi i}{4}}|4\rangle + e^{-\tfrac{18\pi i}{4}}|13\rangle)\quad) = |1\rangle \\[10pt] \end{aligned}

由于计算基态 |1\rangle 是这些本征态的叠加，这意味着如果我们使用状态 |1\rangle 对 U 进行QPE，我们将测量一个相位:

\phi = \frac{s}{r}

其中 s 是 0 和 r-1 之间的随机整数。最后，我们在 \phi 上使用连分式算法来求 r 。电路图看起来像这样(注意，这张图使用了Qiskit的量子比特排序约定):
3

接下来，我们将使用Qiskit的模拟器演示Shor算法。在本演示中，我们将不作解释地提供 U 的电路，但在第4节中，我们将讨论如何有效地构造 U^{2^j} 的电路。

3.使用Qiskit实现

在这个例子中，我们将解决 a=7 和 N=15 的周期查找问题。我们提供 U 的电路，其中:

U|y\rangle = |ay\bmod 15\rangle

没有解释。为了创建 U^x ，我们只需重复电路 x 次。在下一节中，我们将讨论高效创建这些电路的一般方法。函数c_amod15对受控U门进行a次，重复power次。

def c_amod15(a, power):
    """Controlled multiplication by a mod 15"""
    if a not in [2,4,7,8,11,13]:
        raise ValueError("'a' must be 2,4,7,8,11 or 13")
    U = QuantumCircuit(4)        
    for iteration in range(power):
        if a in [2,13]:
            U.swap(2,3)
            U.swap(1,2)
            U.swap(0,1)
        if a in [7,8]:
            U.swap(0,1)
            U.swap(1,2)
            U.swap(2,3)
        if a in [4, 11]:
            U.swap(1,3)
            U.swap(0,2)
        if a in [7,11,13]:
            for q in range(4):
                U.x(q)
    U = U.to_gate()
    U.name = "%i^%i mod 15" % (a, power)
    c_U = U.control()
    return c_U

我们将使用8个计数量子比特:

# Specify variables
n_count = 8  # number of counting qubits
a = 7

我们还导入了用于量子傅里叶变换的电路(你可以在量子傅里叶变换一章中阅读有关量子傅里叶变换的更多信息):

def qft_dagger(n):
    """n-qubit QFTdagger the first n qubits in circ"""
    qc = QuantumCircuit(n)
    # Don't forget the Swaps!
    for qubit in range(n//2):
        qc.swap(qubit, n-qubit-1)
    for j in range(n):
        for m in range(j):
            qc.cp(-np.pi/float(2**(j-m)), m, j)
        qc.h(j)
    qc.name = "QFT†"
    return qc

有了这些构件，我们可以轻松地构建Shor算法的电路:

# Create QuantumCircuit with n_count counting qubits
# plus 4 qubits for U to act on
qc = QuantumCircuit(n_count + 4, n_count)

# Initialize counting qubits
# in state |+>
for q in range(n_count):
    qc.h(q)
    
# And auxiliary register in state |1>
qc.x(n_count)

# Do controlled-U operations
for q in range(n_count):
    qc.append(c_amod15(a, 2**q), 
             [q] + [i+n_count for i in range(4)])

# Do inverse-QFT
qc.append(qft_dagger(n_count), range(n_count))

# Measure circuit
qc.measure(range(n_count), range(n_count))
qc.draw(fold=-1)  # -1 means 'do not fold'

3
3855×322 7.3 KB

让我们看看我们测量的结果:

aer_sim = Aer.get_backend('aer_simulator')
t_qc = transpile(qc, aer_sim)
qobj = assemble(t_qc)
results = aer_sim.run(qobj).result()
counts = results.get_counts()
plot_histogram(counts)

4

由于我们有8个量子比特，这些结果对应于测量的相位:

rows, measured_phases = [], []
for output in counts:
    decimal = int(output, 2)  # Convert (base 2) string to decimal
    phase = decimal/(2**n_count)  # Find corresponding eigenvalue
    measured_phases.append(phase)
    # Add these values to the rows in our table:
    rows.append([f"{output}(bin) = {decimal:>3}(dec)", 
                 f"{decimal}/{2**n_count} = {phase:.2f}"])
# Print the rows in a table
headers=["Register Output", "Phase"]
df = pd.DataFrame(rows, columns=headers)
print(df)

现在我们可以使用连分式算法来尝试找到 s 和 r 。Python内置了这个功能:我们可以使用fractions模块将float对象转换为Fraction对象，例如:

Fraction(0.666)

因为这给出了精确返回结果的分数(在本例中是0.6660000…)，这可能会给出像上面这样粗糙的结果。我们可以使用.limit_denominator()方法来得到与浮点数最相似的分数，而且分母要小于某个特定值:

# Get fraction that most closely resembles 0.666
# with denominator < 15
Fraction(0.666).limit_denominator(15)

好多了!阶数(r)必须小于N，因此我们将最大分母设置为15:

rows = []
for phase in measured_phases:
    frac = Fraction(phase).limit_denominator(15)
    rows.append([phase, f"{frac.numerator}/{frac.denominator}", frac.denominator])
# Print as a table
headers=["Phase", "Fraction", "Guess for r"]
df = pd.DataFrame(rows, columns=headers)
print(df)

我们可以看到，两个特征值的测量结果是正确的: r=4 ，并且我们可以看到Shor算法失败的几率。这些不好的结果是因为 s=0 ，或者是因为 s 和 r 不是互质的，我们得到的因子不是 r 。最简单的解决方法是重复这个实验，直到我们得到一个满意的 r 的结果。

快速练习

  • 修改上面的电路，使 $a$的值为 2,8,11 和 13 。你得到了什么结果?为什么?

4.模幂运算

您可能已经注意到，通过重复 j 来重复 U 门来创建 U^{2^j} 门的方法不会导致多项式时间算法。我们想要一种创建运算符的方法:

U^{2^j}|y\rangle = |a^{2^j}y \bmod N \rangle

它与 j 呈多项式增长。幸运的是,计算:

a^{2^j} \bmod N

高效是可能的。经典计算机可以使用一种称为重复平方的算法来计算指数。在我们的例子中，由于我们只处理 2^j 的指数形式，重复平方算法变得非常简单:

def a2jmodN(a, j, N):
    """Compute a^{2^j} (mod N) by repeated squaring"""
    for i in range(j):
        a = np.mod(a**2, N)
    return a

a2jmodN(7, 2049, 53)

如果在Python中可以实现高效的算法，那么我们可以在量子计算机上使用相同的算法。不幸的是，尽管可以用 j 进行多项式扩展，模幂电路却并不简单，是Shor算法的瓶颈。对初学者友好的实现可以在参考[1]中找到。

5.从周期查找中提取因子

并不是所有的保理问题都是困难的;我们可以立即发现一个偶数，并知道它的一个因数是2。事实上，对于难以因式分解的数，有一些特定的选择标准，但基本思想是选择两个大素数的乘积。

一般的因式分解算法首先会检查是否有因式分解整数的快捷方式(该整数是偶数吗?在使用Shor算法查找最坏情况下的句点之前，计算 N = a^b ?)由于我们的目标是算法的量子部分，我们将直接跳到N是两个素数之积的情况。

Example:Factoring 15
为了看一个分解少量量子比特的例子，我们将分解15，我们都知道它是不太大的质数3和5的乘积。

N=15

第一步是在 1 和 N-1 之间选择一个随机数 a :

np.random.seed(1) # This is to make sure we get reproduceable results
a = randint(2, 15)
print(a)

接下来，我们快速检查它是不是 N 的一个非平凡因子:

from math import gcd # greatest common divisor
gcd(a, N)

太好了。接下来，我们对a=7和N=15进行Shor求阶算法。记住，我们测量的相位为 s/r ，其中:

a^r \bmod N = 1

其中 s 是0到 r-1 之间的随机整数。

def qpe_amod15(a):
    n_count = 8
    qc = QuantumCircuit(4+n_count, n_count)
    for q in range(n_count):
        qc.h(q)     # Initialize counting qubits in state |+>
    qc.x(3+n_count) # And auxiliary register in state |1>
    for q in range(n_count): # Do controlled-U operations
        qc.append(c_amod15(a, 2**q), 
                 [q] + [i+n_count for i in range(4)])
    qc.append(qft_dagger(n_count), range(n_count)) # Do inverse-QFT
    qc.measure(range(n_count), range(n_count))
    # Simulate Results
    aer_sim = Aer.get_backend('aer_simulator')
    # Setting memory=True below allows us to see a list of each sequential reading
    t_qc = transpile(qc, aer_sim)
    qobj = assemble(t_qc, shots=1)
    result = aer_sim.run(qobj, memory=True).result()
    readings = result.get_memory()
    print("Register Reading: " + readings[0])
    phase = int(readings[0],2)/(2**n_count)
    print("Corresponding Phase: %f" % phase)
    return phase

从这个阶段，我们可以很容易地找到 r :

phase = qpe_amod15(a) # Phase = s/r
Fraction(phase).limit_denominator(15) # Denominator should (hopefully!) tell us r

frac = Fraction(phase).limit_denominator(15)
s, r = frac.numerator, frac.denominator
print(r)

现在我们有 r ，我们可以用它找到 N 的因子:

a^r \bmod N = 1

紧接着：

(a^r - 1) \bmod N = 0

也就是说 N 必须除以 a^r-1 。如果$r$也是偶数，那么我们可以这样写:

a^r -1 = (a^{r/2}-1)(a^{r/2}+1)

(如果 r 不是偶数，我们不能进一步，必须对 a 使用不同的值再次尝试)。那么 N 和 a^{r/2}-1 或 a^{r/2}+1 的最大公约数很可能是 N [2]的真因子:

guesses = [gcd(a**(r//2)-1, N), gcd(a**(r//2)+1, N)]
print(guesses)

下面的单元格重复该算法，直到找到至少一个15的因子。您应该尝试重新运行单元格几次，看看它的行为。

a = 7
factor_found = False
attempt = 0
while not factor_found:
    attempt += 1
    print("\nAttempt %i:" % attempt)
    phase = qpe_amod15(a) # Phase = s/r
    frac = Fraction(phase).limit_denominator(N) # Denominator should (hopefully!) tell us r
    r = frac.denominator
    print("Result: r = %i" % r)
    if phase != 0:
        # Guesses for factors are gcd(x^{r/2} ±1 , 15)
        guesses = [gcd(a**(r//2)-1, N), gcd(a**(r//2)+1, N)]
        print("Guessed Factors: %i and %i" % (guesses[0], guesses[1]))
        for guess in guesses:
            if guess not in [1,N] and (N % guess) == 0: # Check to see if guess is a factor
                print("*** Non-trivial factor found: %i ***" % guess)
                factor_found = True

6.参考文献

  1. Stephane Beauregard, Circuit for Shor’s algorithm using 2n+3 qubits, arXiv:quant-ph/0205095

  2. M. Nielsen and I. Chuang, Quantum Computation and Quantum Information, Cambridge Series on Information and the Natural Sciences (Cambridge University Press, Cambridge, 2000). (Page 633)

import qiskit.tools.jupyter
%qiskit_version_table

版本信息

Qiskit Software Version
qiskit-terra 0.22.0
qiskit-aer 0.11.0
qiskit-ibmq-provider 0.19.2
qiskit 0.39.0
qiskit-nature 0.4.1
qiskit-finance 0.3.2
qiskit-optimization 0.4.0
qiskit-machine-learning 0.4.0
System information
Python version 3.8.13
Python compiler Clang 13.1.6 (clang-1316.0.21.2.5)
Python build default, Aug 29 2022 05:17:23
OS Darwin
CPUs 8
Memory (Gb) 32.0
Thu Nov 24 11:55:19 2022 GMT