[PYTHON] What is Shor's algorithm? Is it true that a quantum computer can perform prime factorization in polynomial time? I looked it up!

Introduction

It takes $ O (\ sqrt {N}) $ to perform prime factorization on a classical computer, and even the most efficient algorithm currently being developed takes about an exponential time of digits. It seems that it will end up. However, if you try to factor the prime factors at the 300-digit level, it will take a long time to die. e? Prime factorization in polynomial time? I can do it! That is Shor's algorithm. In this article, the goal of this article is to implement the Shor's algorithm in qiskit and run it on a real quantum computer. This time we will factor 15 down. I think there are probably (or certainly) mistakes and typographical errors in this article, so if you find one, please kindly let me know ...

flow

To explain Shor's algorithm, various preparations are required, but first I would like to write the overall flow.

First, consider a natural number $ N $ that you want to factor into a natural number and a natural number $ x $ that are relatively prime (Euclidean algorithm can be used to determine whether they are relatively prime). Let's think about it in the world of modN below. The smallest r such that $ x ^ r = 1 $ is called an order. It is known that finding this order is effective for prime factorization. However, with a classical calculator, this calculation also takes an exponential time. The phase estimation algorithm solves this problem. The phase estimation algorithm approximates the phase of the eigenvalue (absolute value is 1) of the unitary matrix. Considering U such that $ U ^ r = I $ in the world of $ \ mod N $, r can be obtained from its eigenvalue. The Hadamard test and quantum Fourier transform described below are required as subroutines for this phase estimation algorithm.

Hadamard test

First, I will explain an algorithm called the Hadamard test.

The Hadamard test is realized by the following circuit ($ U $ is an arbitrary gate).

In this circuitq_0 = \left|0\right>, q_1 = \left|\varphi\right>Let's think about what happens when you pass it through. However\left|\varphi\right>IsUIs the eigenvector ofe^{i\lambda}will do. $ \left|0\right>\left|\varphi\right>\rightarrow \frac{1}{\sqrt{2}}(\left|0\right>\left|\varphi\right>+\left|1\right>\left|\varphi\right>)\rightarrow \frac{1}{\sqrt{2}}(\left|0\right>\left|\varphi\right>+e^{i\lambda}\left|1\right>\left|\varphi\right>)\rightarrow \left(\frac{1+e^{i\lambda}}{2}\left|0\right>+\frac{1-e^{i\lambda}}{2}\left|1\right>\right)\otimes\left|\varphi\right> $

control-The qubit after passing through the U gate\frac{1}{\sqrt{2}}(\left|0\right>+e^{i\lambda}\left|1\right>)\otimes \left|\varphi\right>It will be in the state. You can see that the phase of the eigenvalue of U comes out as the phase before (the first qubit). By further applying the Hadamard gate there, the phase of the eigenvalue of U can be obtained as the observation probability of the first qubit. It is strange that the observation probability changes as if the phase has changed from the second qubit, even though no operation has been performed on the control bit. However, if you want to know the phase of the eigenvalues in this way, you have to measure it many times. A phase estimation algorithm will come out as an algorithm to solve this problem.

Quantum Fourier transform

The quantum Fourier transform is an algorithm that performs the discrete Fourier transform. Assuming that the length of the sequence is $ 2 ^ n $, the discrete Fourier transform can be performed with $ O (n2 ^ n) $ in the fast Fourier transform that everyone loves in classical computers, but $ O (n ^ n ^) in the quantum Fourier transform. 2) It can be solved in $ n $ polynomial time called $!

First, let's recall the definition formula of the discrete Fourier transform. $ y_k = \frac{1}{\sqrt{2^n}} \sum_{j=0}^{2^n-1} x_j e^{i\frac{2\pi kj}{2^n}} $ Expressing this on qubit, $ |x\rangle = \sum_{j=0}^{2^n-1} x_j |j\rangle \rightarrow |y\rangle = \sum_{k = 0}^{2^n-1} y_k |k\rangle \\\ |j\rangle \rightarrow \frac{1}{\sqrt{2^n}} \sum_{k=0}^{2^n-1} e^{i\frac{2\pi jk}{2^n}}|k\rangle $ It will be.

In fact, this transformation is unitary, and the quantum computer can approximate any unitary transformation with a universal gate set, so this transformation is feasible on the quantum computer. However, I do not know what kind of gate to bite as it is, so I would like to transform it into a shape that is easier to use. If you expand k in binary and do your best, it will be as follows. $ |j\rangle \rightarrow\otimes_{l = 1}^{2^n} (|0\rangle+e^{\frac{2\pi i}{2^l}}|1\rangle) $ In this form, the Fourier transform is multiplied as the tensor product of each q-bit, so you can see that the Fourier transform can be realized by the Hadamard gate and the rotation gate (see the figure below).

In fact, Shor's algorithm uses the inverse quantum Fourier transform, which is the inverse transform, so this is also shown below. This just hung the gate from the opposite side (but the rotation is reverse).

By the way, you need to be careful about the reading direction of the qubit. Here, we will read in the order of q0q1q2q3. That is, if $ q0 = 1, q1 = q2 = q3 = 0 $, then $ x = 1000 (2) = 8 $. For the time being, H is the same as the Hadamard gate and U1 is the same as the Rz gate, and only $ | 1 \ rangle $ is the gate that has the phase of $ e ^ {i \ theta} $. For more information, click here (https://quantum-computing.ibm.com/support/guides/gate-overview?section=5d00d964853ef8003c6d6820#)

Phase estimation algorithm

The phase algorithm is an improved version of the Hadamard test.

Let the eigenvalue of the unitary matrix $ U $ of size $ 2 ^ n $ be $ e ^ {2 \ pi i \ lambda} $, and the corresponding eigenvector be $ | \ varphi \ rangle $. If $ \ lambda $ can be represented in binary as a decimal of $ m , $ |y\rangle = \frac{1}{\sqrt{2^n}}\sum_{k = 0}^{2^m-1}e^{2\pi ij\lambda}|j\rangle $$ If you can create the state, you can find $ \ lambda $ by inverse discrete Fourier expansion (compare with the definition formula in the previous chapter).

M for phase estimation algorithm+Uses n qubits for input|0\rangle|\varphi\rangleI will put in. First, by applying the Hadamard gate to the first m qubits|0\rangleFrom|2^m-1\rangleMake a superposition state up to. nextiofk桁目が1ならば後ろof|\varphi\rangleToU^kTo apply. Then, on the same principle as the Hadamard test|i\rangle\rightarrow e^{2\pi ik\lambda}|i\rangleToなります。これをk = 1, 2, \cdots mで繰り返すことで、上位mビットTo上of状態を実現することができます！（jを二進数で考えてみてください）回路図は下ofようToなります。回路図を見た方がわかりやすいかもしれません。

Application to prime factorization

Why is it possible to factor the order r once it is found? First, let's assume that r is an even number (it seems that if you take it randomly, it will be an even number with a good probability). Then, the following transformations can be made. $ x^r = 1 \mod N\\ (x^{r/2}-1)(x^{r/2}+1) = 0\mod N $ In other words(x^{r/2}-1)When(x^{r/2}+1)のどちらかはNWhen非自明な公約数を持つこWhenになります。どちらもがN　の倍数Whenなってしまう確率は低いです。公約数はユークリッドの互助法によって、古典計算機で高速に計算できます。そこでU^r = II want to make a matrix that is easy,U|i\rangle = |i\cdot x \mod N\rangleWhenなるように定義すれば良いです。ただしi\geq Nについては何もしないこWhenにします。実はこれもユニタリ行列であるこWhenが示せます。これに位相推定アルゴリズムを用いれば位数rが十分な確率で求まります。一つ問題があるのは、位相推定アルゴリズムでは、固有ベクトルを利用していましたが、今回はわかりません。しかし、入力Whenしては\left|1\right>Is sufficient. Because, if expanded with eigenvectors, the value obtained by the phase estimation algorithm is one of the eigenvalues, and any of thems/r (s = 0, 1, \dots r-1)Whenいう形で表せるため、連分数アルゴリズムによってrを高速に求められるからです(s = 0Ignore if)。連分数アルゴリズムは小数を連分数で近似するこWhenでもっWhenもらしい分数Whenしての形を推定します。詳しいこWhenはcontinued fraction algorithmでググりましょう。一応実装は後に載せてあります。

Implementation in Qiskit

Well, here is the production. Let's run the algorithm we saw above in a simulator and factor $ N = 15 $ by $ x = 4 $. This time, we will use a quantum calculation library called qiskit. For more information on qiskit, see Official Documentation. All execution is done on jupyter notebook.

The quantum computer part is as follows. The implementation of U mentioned in the previous chapter, but since it is difficult to implement general modulo operation, we have taken advantage of the fact that $ N = 15, x = 4 $ ($ {U} ^ {2 ^). For i} $, I've also used the fact that the order is 2, so I've omitted it ...), but this time it's okay because it just gets a feel for it.

from qiskit import *
from qiskit.providers.ibmq import least_busy
from qiskit.visualization import plot_histogram
from qiskit.tools.monitor import job_monitor
import math

q = QuantumRegister(8, "q")
c = ClassicalRegister(4, "c")
circuit = QuantumCircuit(q, c)
circuit.x(7)
#Hadamard transform
for i in range(4):
    circuit.h(i)
#U4^1 gate
circuit.cswap(3, 4, 6)
circuit.cswap(3, 5, 7)
circuit.barrier()
#U4^2 gate
#U4^4 gate
#U4^8 gate
circuit.barrier()
#inverse quantum Fourier transform
circuit.swap(0, 3)
circuit.swap(1, 2)
circuit.h(3)
circuit.cu1(-np.pi/2, 3, 2)
circuit.h(2)
circuit.barrier()
circuit.cu1(-np.pi/4, 3, 1)
circuit.cu1(-np.pi/2, 2, 1)
circuit.h(1)
circuit.barrier()
circuit.cu1(-np.pi/8, 3, 0)
circuit.cu1(-np.pi/4, 2, 0)
circuit.cu1(-np.pi/2, 1, 0)
circuit.h(0)
circuit.barrier()
#measure
for i in range(4):
    circuit.measure(q[i], c[3-i])#Binary x=I try to be q1q2q3q4
circuit.draw(output='mpl')

Now, the rest is the classical calculation part. I implemented continued fraction expansion as follows (it's quite suspicious because I just wrote it properly and tried a little sample ...). I try to end when the relative error is less than or equal to eps.

eps = 0.01
def continued_fractions_algorithm(x):
    res = [int(x)]
    x-=int(x)
    if x/(res[0]+0.1)>eps:
        a = continued_fractions_algorithm(1/x)
        res+=a
    return res

Now, let's estimate the amount of calculation (as a quantum computer). Assuming that the natural number you want to factor into is $ n $ bit, the discrete Fourier transform part is $ O (n ^ 2) $ and the phase estimation algorithm part is $ O (n ^ 3) $. The classic part doesn't become a bottleneck, so I was able to solve it in polynomial time of $ O (n ^ 3) $ in total! After that, if you execute the following code,

def shor_algorithm(use_simulator):
    if use_simulator:
        backend = BasicAer.get_backend('qasm_simulator')
    else:
        backend = IBMQ.get_provider().get_backend('ibmq_16_melbourne')
    flag = False
    job = execute(circuit, backend, shots = N)
    job_monitor(job)
    result = job.result()
    counts = result.get_counts(circuit)
    measures = np.array(list(map(lambda x:int(x, 2), counts.keys())), dtype = np.int64)
    probs = np.array(list(counts.values()))/N
    for i in range(5):
        output = np.random.choice(measures, p = probs)
        a = continued_fractions_algorithm(output/2**4)
        r , s =1,  a[-1]
        for i in range(len(a)-1)[::-1]:
            r, s = s, a[i]*s+r 
        if r % 15 == 0 or s == 0:
            continue
        d = math.gcd(15, 4**(r-1)-1)
        if d != 1:
            flag = True
            break
    plot_histogram(result.get_counts())
    if flag:
        print('{0} devides 15 with r = {1}'.format(d, r))
    else:
        print('the algorithm failed')    
    return result
 
%matplotlib inline
use_simulator = True
result = shor_algorithm(use_simulator)
plot_histogram(result.get_counts())

And succeeded in factoring!

Calculation on the actual machine

You can operate IBM Q, a quantum computer that can be used for free from qiskit. The official page is here. You need to create an account for the operation.

First, execute the following code.

from qiskit import IBMQ
my_token = ""
IBMQ.save_account(my_token)
provider = IBMQ.load_account()

In my_token, enter the token obtained from My Account on the official page.

Since I am using 8qubit this time, I will use ibmq_16_melbourne. It will take some time, but if you run it with use_simulator = False in the above code,

What?

Summary

How was it? I tried factorization with a quantum computer this time, but as we saw in the previous chapter, there is a considerable error in the actual machine. Of course, this time I'm trying something simple that anyone can use, but it seems that it will take time for quantum computers to be put into practical use if even factorization of 15 is difficult.

reference

Quantum Native Dojo https://dojo.qulacs.org/ja/latest/ Kenjiro Miyano, Akira Furusawa, Introduction to Quantum Computers, Nihon Hyoronsha 2016 Michael A. Nilsen,Isaac L. Chuang Quantum Computation and Quantum Information 10th Anniversary Edition 2010