쇼어 알고리즘

사용 예상 시간: Eagle r3 프로세서에서 약 3초 (참고: 이는 예상치이며, 실제 실행 시간은 다를 수 있습니다.)

쇼어 알고리즘은 1994년 Peter Shor가 개발한 것으로, 정수를 다항식 시간 내에 인수분해하는 획기적인 양자 알고리즘입니다. 이 알고리즘의 중요성은 알려진 어떠한 고전 알고리즘보다 지수적으로 빠르게 큰 정수를 인수분해할 수 있다는 점에 있으며, 이는 큰 수의 인수분해가 어렵다는 점에 의존하는 RSA와 같은 광범위하게 사용되는 암호화 시스템의 보안을 위협합니다. 충분히 강력한 양자 컴퓨터에서 이 문제를 효율적으로 해결함으로써 쇼어 알고리즘은 암호학, 사이버 보안, 계산 수학 등의 분야를 혁신할 수 있으며, 이는 양자 계산의 변혁적인 힘을 잘 보여줍니다.

이 튜토리얼은 양자 컴퓨터에서 15를 인수분해함으로써 쇼어 알고리즘을 시연하는 데 초점을 맞춥니다.

먼저 위수 찾기 문제(order finding problem)를 정의하고 양자 위상 추정 프로토콜로부터 해당 Circuit을 구성합니다. 다음으로, 트랜스파일할 수 있는 가장 짧은 깊이의 Circuit을 사용하여 실제 하드웨어에서 위수 찾기 Circuit을 실행합니다. 마지막 섹션에서는 위수 찾기 문제를 정수 인수분해와 연결하여 쇼어 알고리즘을 완성합니다.

튜토리얼 마지막 부분에서는 실제 하드웨어에서의 다른 쇼어 알고리즘 시연 사례들을 살펴보며, 범용 구현과 15나 21과 같은 특정 정수의 인수분해에 맞춰진 구현 모두에 초점을 맞춥니다. 참고: 이 튜토리얼은 쇼어 알고리즘과 관련된 Circuit의 구현 및 시연에 더 집중합니다. 해당 내용에 대한 심층적인 교육 자료는 John Watrous 박사의 양자 알고리즘 기초 강좌와 참고문헌 섹션의 논문들을 참조해 주세요.

요구 사항

이 튜토리얼을 시작하기 전에 다음이 설치되어 있는지 확인하세요:

• 시각화 지원이 포함된 Qiskit SDK v2.0 이상
• Qiskit Runtime v0.40 이상 (pip install qiskit-ibm-runtime)

설정

# Added by doQumentation — required packages for this notebook
!pip install -q numpy pandas qiskit qiskit-ibm-runtime

import numpy as np
import pandas as pd
from fractions import Fraction
from math import floor, gcd, log

from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister
from qiskit.circuit.library import QFT, UnitaryGate
from qiskit.transpiler import CouplingMap, generate_preset_pass_manager
from qiskit.visualization import plot_histogram

from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_ibm_runtime import SamplerV2 as Sampler

1단계: 고전적 입력을 양자 문제로 변환하기

배경

쇼어의 정수 인수분해 알고리즘은 위수 찾기 문제라는 중간 문제를 활용합니다. 이 섹션에서는 양자 위상 추정을 사용하여 위수 찾기 문제를 푸는 방법을 시연합니다.

위상 추정 문제

위상 추정 문제에서는 $n$개의 Qubit으로 이루어진 양자 상태 $\ket{\psi}$와 $n$개의 Qubit에 작용하는 유니터리 양자 Circuit이 주어집니다. $\ket{\psi}$는 해당 Circuit의 작용을 설명하는 유니터리 행렬 $U$의 고유벡터임이 보장되며, 우리의 목표는 $\ket{\psi}$에 대응하는 고유값 $\lambda = e^{2 \pi i \theta}$를 계산하거나 근사하는 것입니다. 다시 말해, Circuit은 다음을 만족하는 $\theta \in [0, 1)$에 대한 근삿값을 출력해야 합니다. $U \ket{\psi}= e^{2 \pi i \theta} \ket{\psi}.$ 위상 추정 Circuit의 목표는 $m$비트로 $\theta$를 근사하는 것입니다. 수학적으로는 $\theta \approx y / 2^m$을 만족하는 $y$를 찾고자 하며, 여기서 $y \in {0, 1, 2, \dots, 2^{m-1}}$입니다. 다음 이미지는 $m$개의 Qubit에 대한 측정을 수행하여 $m$비트로 $y$를 추정하는 양자 Circuit을 보여줍니다. 위 Circuit에서 상위 $m$개의 Qubit은 $\ket{0^m}$ 상태로 초기화되고, 하위 $n$개의 Qubit은 $U$의 고유벡터임이 보장된 $\ket{\psi}$로 초기화됩니다. 위상 추정 Circuit의 첫 번째 구성 요소는 제어 Qubit으로 *위상 킥백(phase kickback)*을 수행하는 제어-유니터리 연산입니다. 이 제어 유니터리들은 최하위 비트부터 최상위 비트까지 제어 Qubit의 위치에 따라 지수적으로 적용됩니다. $\ket{\psi}$가 $U$의 고유벡터이므로, 하위 $n$개의 Qubit 상태는 이 연산에 의해 영향을 받지 않지만, 고유값의 위상 정보는 상위 $m$개의 Qubit으로 전파됩니다. 제어 유니터리를 통한 위상 킥백 연산 이후, 상위 $m$개의 Qubit의 모든 가능한 상태는 유니터리 $U$의 각 고유벡터 $\ket{\psi}$에 대해 서로 정규직교함이 밝혀져 있습니다. 따라서 이 상태들은 완벽하게 구별 가능하며, 측정을 위해 이들이 형성하는 기저를 계산 기저로 다시 회전시킬 수 있습니다. 수학적 분석에 따르면 이 회전 행렬은 $2^m$차원 힐베르트 공간에서의 역 양자 푸리에 변환(QFT)에 해당합니다. 이에 대한 직관은, 모듈러 지수 연산자의 주기적 구조가 양자 상태에 인코딩되고, QFT가 이 주기성을 주파수 도메인의 측정 가능한 피크로 변환한다는 것입니다.

쇼어 알고리즘에서 QFT Circuit이 사용되는 이유에 대한 더 깊은 이해를 위해서는 양자 알고리즘 기초 강좌를 참조해 주세요. 이제 위수 찾기에 위상 추정 Circuit을 사용할 준비가 되었습니다.

위수 찾기 문제

위수 찾기 문제를 정의하기 위해 몇 가지 수론 개념부터 시작합니다. 먼저, 임의의 양의 정수 $N$에 대해 집합 $\mathbb{Z}_N$을 다음과 같이 정의합니다. $\mathbb{Z}_N = \{0, 1, 2, \dots, N-1\}.$ $\mathbb{Z}_N$에서의 모든 산술 연산은 $N$을 법으로(modulo $N$) 수행됩니다. 특히, $N$과 서로소인 모든 원소 $a \in \mathbb{Z}_n$은 특별하며 다음과 같이 $\mathbb{Z}^*_N$을 구성합니다. $\mathbb{Z}^*_N = \{ a \in \mathbb{Z}_N : \mathrm{gcd}(a, N)=1 \}.$ 원소 $a \in \mathbb{Z}^*_N$에 대해 $a^r \equiv 1 \; (\mathrm{mod} \; N)$을 만족하는 가장 작은 양의 정수 $r$을 $a$의 $N$에 대한 *위수(order)*라고 정의합니다. 나중에 살펴보겠지만, $a \in \mathbb{Z}^*_N$의 위수를 찾으면 $N$을 인수분해할 수 있습니다. 위상 추정 Circuit으로부터 위수 찾기 Circuit을 구성하기 위해서는 두 가지를 고려해야 합니다. 첫째, 위수 $r$을 찾을 수 있게 해주는 유니터리 $U$를 정의해야 하고, 둘째, 위상 추정 Circuit의 초기 상태를 준비하기 위한 $U$의 고유벡터 $\ket{\psi}$를 정의해야 합니다.

위수 찾기 문제를 위상 추정과 연결하기 위해, 고전적 상태가 $\mathbb{Z}_N$에 대응하는 시스템에 정의된 연산을 고려합니다. 여기서 고정된 원소 $a \in \mathbb{Z}^*_N$을 곱합니다. 특히, 각 $x \in \mathbb{Z}_N$에 대해 $M_a \ket{x} = \ket{ax \; (\mathrm{mod} \; N)}$이 되도록 이 곱셈 연산자 $M_a$를 정의합니다. 방정식 우변의 ket 내부에서 $N$을 법으로 곱을 취한다는 것이 암묵적으로 포함되어 있음에 주의하세요. 수학적 분석에 따르면 $M_a$는 유니터리 연산자입니다. 더 나아가, $M_a$는 $a$의 위수 $r$을 위상 추정 문제와 연결할 수 있게 해주는 고유벡터 및 고유값 쌍을 가집니다. 구체적으로, 임의의 $j \in \{0, \dots, r-1\}$ 선택에 대해 $\ket{\psi_j} = \frac{1}{\sqrt{r}} \sum^{r-1}_{k=0} \omega^{-jk}_{r} \ket{a^k}$는 대응하는 고유값이 $\omega^{j}_{r}$인 $M_a$의 고유벡터입니다. 여기서 $\omega^{j}_{r} = e^{2 \pi i \frac{j}{r}}.$ 관찰에 의하면, 편리한 고유벡터/고유값 쌍은 $\omega^{1}_{r} = e^{2 \pi i \frac{1}{r}}$을 갖는 상태 $\ket{\psi_1}$임을 알 수 있습니다. 따라서 고유벡터 $\ket{\psi_1}$를 찾을 수 있다면, 양자 Circuit으로 위상 $\theta=1/r$을 추정하여 위수 $r$의 추정치를 얻을 수 있습니다. 하지만 이것이 쉽지 않으므로 대안을 고려해야 합니다.

계산 기저 상태 $\ket{1}$을 초기 상태로 준비하면 Circuit이 어떤 결과를 낼지 생각해 보겠습니다. 이것은 $M_a$의 고유 상태가 아니지만, 앞서 설명한 고유 상태들의 균일한 중첩입니다. 다시 말해, 다음 관계가 성립합니다. $\ket{1} = \frac{1}{\sqrt{r}} \sum^{r-1}_{k=0} \ket{\psi_k}$ 위 방정식의 의미는, 초기 상태를 $\ket{1}$로 설정하면 $k \in \{ 0, \dots, r-1\}$을 균일하게 임의로 선택하여 위상 추정 Circuit에서 고유벡터로 $\ket{\psi_k}$를 사용한 것과 정확히 동일한 측정 결과를 얻는다는 것입니다. 즉, 상위 $m$개의 Qubit에 대한 측정은 균일하게 임의로 선택된 $k \in \{ 0, \dots, r-1\}$에서 값 $k / r$에 대한 근삿값 $y / 2^m$을 산출합니다. 이를 통해 여러 번의 독립적인 실행 후 높은 신뢰도로 $r$을 알아낼 수 있으며, 이것이 우리의 목표였습니다.

모듈러 지수 연산자

지금까지 우리는 양자 Circuit에서 $U = M_a$와 $\ket{\psi} = \ket{1}$을 정의함으로써 위상 추정 문제를 위수 찾기 문제와 연결했습니다. 따라서 마지막으로 남은 것은 $k = 1, 2, 4, \dots, 2^{m-1}$에 대해 $M_a^k$로서 $M_a$의 모듈러 지수를 효율적으로 정의하는 방법을 찾는 것입니다. 이 계산을 수행하기 위해, 선택한 임의의 거듭제곱 $k$에 대해 $M_a$ Circuit을 $k$번 반복하는 대신 $b = a^k \; \mathrm{mod} \; N$을 계산한 다음 $M_b$ Circuit을 사용하여 $M_a^k$에 대한 Circuit을 생성할 수 있음을 알 수 있습니다. 2의 거듭제곱에 해당하는 멱수만 필요하므로, 반복 제곱법을 사용하여 이를 고전적으로 효율적으로 수행할 수 있습니다.

2단계: 양자 하드웨어 실행을 위한 문제 최적화

$N = 15$, $a=2$를 이용한 구체적 예시

여기서 잠시 구체적인 예시를 살펴보고 $N=15$에 대한 위수 찾기 Circuit을 구성해 보겠습니다. $N=15$에 대해 가능한 비자명 원소 $a \in \mathbb{Z}_N^*$는 $a \in \{2, 4, 7, 8, 11, 13, 14 \}$임에 유의하세요. 이 예시에서는 $a=2$를 선택합니다. $M_2$ 연산자와 모듈러 지수 연산자 $M_2^k$를 구성하겠습니다. $M_2$의 계산 기저 상태에 대한 작용은 다음과 같습니다. $M_2 \ket{0} = \ket{0} \quad M_2 \ket{5} = \ket{10} \quad M_2 \ket{10} = \ket{5}$ $M_2 \ket{1} = \ket{2} \quad M_2 \ket{6} = \ket{12} \quad M_2 \ket{11} = \ket{7}$ $M_2 \ket{2} = \ket{4} \quad M_2 \ket{7} = \ket{14} \quad M_2 \ket{12} = \ket{9}$ $M_2 \ket{3} = \ket{6} \quad M_2 \ket{8} = \ket{1} \quad M_2 \ket{13} = \ket{11}$ $M_2 \ket{4} = \ket{8} \quad M_2 \ket{9} = \ket{3} \quad M_2 \ket{14} = \ket{13}$ 관찰에 의하면 기저 상태들이 뒤섞이므로 순열 행렬임을 알 수 있습니다. 이 연산을 스왑 Gate를 사용하여 4개의 Qubit으로 구성할 수 있습니다. 아래에서 $M_2$와 제어-$M_2$ 연산을 구성합니다.

def M2mod15():
    """
    M2 (mod 15)
    """
    b = 2
    U = QuantumCircuit(4)

    U.swap(2, 3)
    U.swap(1, 2)
    U.swap(0, 1)

    U = U.to_gate()
    U.name = f"M_{b}"

    return U

# Get the M2 operator
M2 = M2mod15()

# Add it to a circuit and plot
circ = QuantumCircuit(4)
circ.compose(M2, inplace=True)
circ.decompose(reps=2).draw(output="mpl", fold=-1)

def controlled_M2mod15():
    """
    Controlled M2 (mod 15)
    """
    b = 2
    U = QuantumCircuit(4)

    U.swap(2, 3)
    U.swap(1, 2)
    U.swap(0, 1)

    U = U.to_gate()
    U.name = f"M_{b}"
    c_U = U.control()

    return c_U

# Get the controlled-M2 operator
controlled_M2 = controlled_M2mod15()

# Add it to a circuit and plot
circ = QuantumCircuit(5)
circ.compose(controlled_M2, inplace=True)
circ.decompose(reps=1).draw(output="mpl", fold=-1)

2개 이상의 Qubit에 작용하는 Gate는 2-Qubit Gate로 추가 분해됩니다.

circ.decompose(reps=2).draw(output="mpl", fold=-1)

이제 모듈러 지수 연산자를 구성해야 합니다. 위상 추정에서 충분한 정밀도를 얻기 위해 추정 측정에 8개의 Qubit을 사용합니다. 따라서 각 $k = 0, 1, \dots, 7$에 대해 $b = a^{2^k} \; (\mathrm{mod} \; N)$을 갖는 $M_b$를 구성해야 합니다.

def a2kmodN(a, k, N):
    """Compute a^{2^k} (mod N) by repeated squaring"""
    for _ in range(k):
        a = int(np.mod(a**2, N))
    return a

k_list = range(8)
b_list = [a2kmodN(2, k, 15) for k in k_list]

print(b_list)

[2, 4, 1, 1, 1, 1, 1, 1]

$b$ 값 목록에서 볼 수 있듯이, 이전에 구성한 $M_2$ 외에도 $M_4$와 $M_1$을 구성해야 합니다. $M_1$은 계산 기저 상태에 자명하게 작용하므로 단순히 항등 연산자임에 유의하세요.

$M_4$의 계산 기저 상태에 대한 작용은 다음과 같습니다. $M_4 \ket{0} = \ket{0} \quad M_4 \ket{5} = \ket{5} \quad M_4 \ket{10} = \ket{10}$ $M_4 \ket{1} = \ket{4} \quad M_4 \ket{6} = \ket{9} \quad M_4 \ket{11} = \ket{14}$ $M_4 \ket{2} = \ket{8} \quad M_4 \ket{7} = \ket{13} \quad M_4 \ket{12} = \ket{3}$ $M_4 \ket{3} = \ket{12} \quad M_4 \ket{8} = \ket{2} \quad M_4 \ket{13} = \ket{7}$ $M_4 \ket{4} = \ket{1} \quad M_4 \ket{9} = \ket{6} \quad M_4 \ket{14} = \ket{11}$

따라서 이 순열은 다음 스왑 연산으로 구성할 수 있습니다.

def M4mod15():
    """
    M4 (mod 15)
    """
    b = 4
    U = QuantumCircuit(4)

    U.swap(1, 3)
    U.swap(0, 2)

    U = U.to_gate()
    U.name = f"M_{b}"

    return U

# Get the M4 operator
M4 = M4mod15()

# Add it to a circuit and plot
circ = QuantumCircuit(4)
circ.compose(M4, inplace=True)
circ.decompose(reps=2).draw(output="mpl", fold=-1)

def controlled_M4mod15():
    """
    Controlled M4 (mod 15)
    """
    b = 4
    U = QuantumCircuit(4)

    U.swap(1, 3)
    U.swap(0, 2)

    U = U.to_gate()
    U.name = f"M_{b}"
    c_U = U.control()

    return c_U

# Get the controlled-M4 operator
controlled_M4 = controlled_M4mod15()

# Add it to a circuit and plot
circ = QuantumCircuit(5)
circ.compose(controlled_M4, inplace=True)
circ.decompose(reps=1).draw(output="mpl", fold=-1)

두 큐비트 이상에 작용하는 Gate는 추가로 2-큐비트 Gate로 분해됩니다.

circ.decompose(reps=2).draw(output="mpl", fold=-1)

앞서 살펴본 것처럼, 주어진 $b \in \mathbb{Z}^*_N$에 대한 $M_b$ 연산자는 순열(permutation) 연산입니다. 여기서 다루는 순열 문제의 크기가 비교적 작기 때문에($N=15$는 4개의 큐비트만 필요합니다), 직관적으로 SWAP Gate를 이용해 이 연산을 직접 합성할 수 있었습니다. 하지만 일반적으로 이 방법은 확장 가능한 접근법이 아닐 수 있습니다. 대신 순열 행렬을 명시적으로 구성하고, Qiskit의 UnitaryGate 클래스와 트랜스파일 메서드를 사용해 순열 행렬을 합성해야 할 수도 있습니다. 그러나 이 방법은 Circuit 깊이가 상당히 깊어질 수 있습니다. 아래에 예시가 나와 있습니다.

def mod_mult_gate(b, N):
    """
    Modular multiplication gate from permutation matrix.
    """
    if gcd(b, N) > 1:
        print(f"Error: gcd({b},{N}) > 1")
    else:
        n = floor(log(N - 1, 2)) + 1
        U = np.full((2**n, 2**n), 0)
        for x in range(N):
            U[b * x % N][x] = 1
        for x in range(N, 2**n):
            U[x][x] = 1
        G = UnitaryGate(U)
        G.name = f"M_{b}"
        return G

# Let's build M2 using the permutation matrix definition
M2_other = mod_mult_gate(2, 15)

# Add it to a circuit
circ = QuantumCircuit(4)
circ.compose(M2_other, inplace=True)
circ = circ.decompose()

# Transpile the circuit and get the depth
coupling_map = CouplingMap.from_line(4)
pm = generate_preset_pass_manager(coupling_map=coupling_map)
transpiled_circ = pm.run(circ)

print(f"qubits: {circ.num_qubits}")
print(
    f"2q-depth: {transpiled_circ.depth(lambda x: x.operation.num_qubits==2)}"
)
print(f"2q-size: {transpiled_circ.size(lambda x: x.operation.num_qubits==2)}")
print(f"Operator counts: {transpiled_circ.count_ops()}")
transpiled_circ.decompose().draw(
    output="mpl", fold=-1, style="clifford", idle_wires=False
)

qubits: 4
2q-depth: 94
2q-size: 96
Operator counts: OrderedDict({'cx': 45, 'swap': 32, 'u': 24, 'u1': 7, 'u3': 4, 'unitary': 3, 'circuit-335': 1, 'circuit-338': 1, 'circuit-341': 1, 'circuit-344': 1, 'circuit-347': 1, 'circuit-350': 1, 'circuit-353': 1, 'circuit-356': 1, 'circuit-359': 1, 'circuit-362': 1, 'circuit-365': 1, 'circuit-368': 1, 'circuit-371': 1, 'circuit-374': 1, 'circuit-377': 1, 'circuit-380': 1})

이 카운트를 $M_2$ Gate의 수동 구현 방식으로 컴파일된 Circuit 깊이와 비교해 보겠습니다.

# Get the M2 operator from our manual construction
M2 = M2mod15()

# Add it to a circuit
circ = QuantumCircuit(4)
circ.compose(M2, inplace=True)
circ = circ.decompose(reps=3)

# Transpile the circuit and get the depth
coupling_map = CouplingMap.from_line(4)
pm = generate_preset_pass_manager(coupling_map=coupling_map)
transpiled_circ = pm.run(circ)

print(f"qubits: {circ.num_qubits}")
print(
    f"2q-depth: {transpiled_circ.depth(lambda x: x.operation.num_qubits==2)}"
)
print(f"2q-size: {transpiled_circ.size(lambda x: x.operation.num_qubits==2)}")
print(f"Operator counts: {transpiled_circ.count_ops()}")
transpiled_circ.draw(
    output="mpl", fold=-1, style="clifford", idle_wires=False
)

qubits: 4
2q-depth: 9
2q-size: 9
Operator counts: OrderedDict({'cx': 9})

보면 알 수 있듯이, 순열 행렬 방식은 수동 구현 방식에 비해 단일 $M_2$ Gate만으로도 Circuit 깊이가 상당히 깊어졌습니다. 따라서 이전에 구현한 $M_b$ 연산 방식을 계속 사용하겠습니다. 이제 앞서 정의한 제어형 모듈러 지수(controlled modular exponentiation) 연산자를 사용해 전체 위수 탐색 Circuit을 구성할 준비가 되었습니다. 아래 코드에서는 Qiskit Circuit 라이브러리에서 QFT Circuit도 가져옵니다. QFT Circuit은 각 Qubit에 Hadamard Gate를 적용하고, 일련의 제어형 U1(또는 위상에 따라 Z) Gate와 SWAP Gate 레이어를 사용합니다.

# Order finding problem for N = 15 with a = 2
N = 15
a = 2

# Number of qubits
num_target = floor(log(N - 1, 2)) + 1  # for modular exponentiation operators
num_control = 2 * num_target  # for enough precision of estimation

# List of M_b operators in order
k_list = range(num_control)
b_list = [a2kmodN(2, k, 15) for k in k_list]

# Initialize the circuit
control = QuantumRegister(num_control, name="C")
target = QuantumRegister(num_target, name="T")
output = ClassicalRegister(num_control, name="out")
circuit = QuantumCircuit(control, target, output)

# Initialize the target register to the state |1>
circuit.x(num_control)

# Add the Hadamard gates and controlled versions of the
# multiplication gates
for k, qubit in enumerate(control):
    circuit.h(k)
    b = b_list[k]
    if b == 2:
        circuit.compose(
            M2mod15().control(), qubits=[qubit] + list(target), inplace=True
        )
    elif b == 4:
        circuit.compose(
            M4mod15().control(), qubits=[qubit] + list(target), inplace=True
        )
    else:
        continue  # M1 is the identity operator

# Apply the inverse QFT to the control register
circuit.compose(QFT(num_control, inverse=True), qubits=control, inplace=True)

# Measure the control register
circuit.measure(control, output)

circuit.draw("mpl", fold=-1)

나머지 제어 Qubit에서의 제어형 모듈러 지수 연산은 $M_1$이 항등(identity) 연산자이므로 생략하였습니다. 이 튜토리얼의 후반부에서는 이 Circuit을 ibm_marrakesh Backend에서 실행할 것입니다. 이를 위해 해당 특정 Backend에 맞게 Circuit을 트랜스파일하고, Circuit 깊이 및 Gate 카운트를 보고합니다.

service = QiskitRuntimeService()
backend = service.backend("ibm_marrakesh")
pm = generate_preset_pass_manager(optimization_level=2, backend=backend)

transpiled_circuit = pm.run(circuit)

print(
    f"2q-depth: {transpiled_circuit.depth(lambda x: x.operation.num_qubits==2)}"
)
print(
    f"2q-size: {transpiled_circuit.size(lambda x: x.operation.num_qubits==2)}"
)
print(f"Operator counts: {transpiled_circuit.count_ops()}")
transpiled_circuit.draw(
    output="mpl", fold=-1, style="clifford", idle_wires=False
)

2q-depth: 187
2q-size: 260
Operator counts: OrderedDict({'sx': 521, 'rz': 354, 'cz': 260, 'measure': 8, 'x': 4})

3단계: Qiskit 프리미티브를 사용한 실행

먼저, 이 Circuit을 이상적인 시뮬레이터에서 실행했을 때 이론적으로 얻을 수 있는 결과에 대해 논의합니다. 아래에는 1024 shots으로 위 Circuit을 시뮬레이션한 결과가 있습니다. 보면 알 수 있듯이, 제어 Qubit에 대한 네 개의 비트열에 걸쳐 거의 균일한 분포를 얻게 됩니다.

# Obtained from the simulator
counts = {"00000000": 264, "01000000": 268, "10000000": 249, "11000000": 243}

plot_histogram(counts)

제어 Qubit을 측정함으로써 $M_a$ 연산자의 8비트 위상 추정값을 얻습니다. 이 이진 표현을 십진수로 변환하여 측정된 위상을 구할 수 있습니다. 위의 히스토그램에서 볼 수 있듯이 네 개의 서로 다른 비트열이 측정되었으며, 각각은 다음과 같이 위상 값에 대응합니다.

# Rows to be displayed in table
rows = []
# Corresponding phase of each bitstring
measured_phases = []

for output in counts:
    decimal = int(output, 2)  # Convert bitstring to decimal
    phase = decimal / (2**num_control)  # 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 ** num_control} = {phase:.2f}",
        ]
    )

# Print the rows in a table
headers = ["Register Output", "Phase"]
df = pd.DataFrame(rows, columns=headers)
print(df)

Register Output           Phase
0  00000000(bin) =   0(dec)    0/256 = 0.00
1  01000000(bin) =  64(dec)   64/256 = 0.25
2  10000000(bin) = 128(dec)  128/256 = 0.50
3  11000000(bin) = 192(dec)  192/256 = 0.75

측정된 위상은 $\theta = k / r$ ($k$는 $\{0, 1, \dots, r-1 \}$에서 균일하게 무작위로 샘플링됨)에 해당한다는 점을 기억하세요. 따라서 연분수 알고리즘을 사용하여 $k$와 위수 $r$을 구할 수 있습니다. Python에는 이 기능이 내장되어 있습니다. fractions 모듈을 사용하여 부동소수점 수를 Fraction 객체로 변환할 수 있습니다. 예를 들어:

Fraction(0.666)

Fraction(5998794703657501, 9007199254740992)

이 방법은 결과를 정확하게 반환하는 분수를 제공하기 때문에(이 경우 0.6660000...), 위와 같이 복잡한 결과가 나올 수 있습니다. .limit_denominator() 메서드를 사용하면 분모가 특정 값 이하이면서 부동소수점 수에 가장 근접한 분수를 구할 수 있습니다:

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

Fraction(2, 3)

훨씬 깔끔합니다. 위수(r)는 N보다 작아야 하므로, 최대 분모를 15로 설정합니다:

# Rows to be displayed in a table
rows = []

for phase in measured_phases:
    frac = Fraction(phase).limit_denominator(15)
    rows.append(
        [phase, f"{frac.numerator}/{frac.denominator}", frac.denominator]
    )

# Print the rows in a table
headers = ["Phase", "Fraction", "Guess for r"]
df = pd.DataFrame(rows, columns=headers)
print(df)

Phase Fraction  Guess for r
0   0.00      0/1            1
1   0.25      1/4            4
2   0.50      1/2            2
3   0.75      3/4            4

측정된 고유값 중 두 개가 올바른 결과 $r=4$를 제공했음을 알 수 있으며, Shor 알고리즘의 위수 찾기가 실패할 가능성이 있다는 것도 확인할 수 있습니다. 좋지 않은 결과가 나오는 이유는 $k = 0$이거나, $k$와 $r$이 서로소가 아니기 때문입니다. 이 경우 $r$ 대신 $r$의 약수가 주어집니다. 가장 간단한 해결책은 만족스러운 $r$ 값을 얻을 때까지 실험을 반복하는 것입니다.

지금까지 시뮬레이터에서 위상 추정 Circuit을 사용하여 $a=2$로 $N=15$의 위수 찾기 문제를 구현했습니다. Shor 알고리즘의 마지막 단계는 위수 찾기 문제를 정수 인수분해 문제와 연결하는 것입니다. 알고리즘의 이 마지막 부분은 순수하게 고전적이며, 양자 컴퓨터에서 위상 측정값을 얻은 후 고전 컴퓨터에서 수행할 수 있습니다. 따라서, 위수 찾기 Circuit을 실제 하드웨어에서 실행하는 방법을 보여준 후에 알고리즘의 마지막 부분을 처리하겠습니다.

하드웨어 실행

이제 ibm_marrakesh에 맞게 트랜스파일한 위수 찾기 Circuit을 실행할 수 있습니다. 여기서는 오류 억제를 위해 동적 디커플링(DD)을, 오류 완화를 위해 게이트 트월링을 활용합니다. DD는 정밀하게 타이밍된 제어 펄스 시퀀스를 양자 장치에 인가하여 원하지 않는 환경과의 상호작용 및 디코히어런스를 효과적으로 평균화합니다. 반면 게이트 트월링은 특정 양자 Gate를 무작위화하여 결맞음 오류를 2차가 아닌 1차적으로 누적되는 파울리 오류로 변환합니다. 두 기법을 결합하면 양자 계산의 결맞음과 충실도를 더욱 향상시킬 수 있습니다.

# Sampler primitive to obtain the probability distribution
sampler = Sampler(backend)

# Turn on dynamical decoupling with sequence XpXm
sampler.options.dynamical_decoupling.enable = True
sampler.options.dynamical_decoupling.sequence_type = "XpXm"
# Enable gate twirling
sampler.options.twirling.enable_gates = True

# Assign tags before executing
sampler.options.environment.job_tags = ["TUT_SA"]

pub = transpiled_circuit
job = sampler.run([pub], shots=1024)

result = job.result()[0]
counts = result.data["out"].get_counts()

plot_histogram(counts, figsize=(35, 5))

보면 알 수 있듯이, 가장 높은 카운트를 가진 동일한 비트열들을 얻었습니다. 양자 하드웨어에는 노이즈가 있으므로 다른 비트열로 약간의 누설이 발생하며, 이는 통계적으로 필터링할 수 있습니다.

# Dictionary of bitstrings and their counts to keep
counts_keep = {}
# Threshold to filter
threshold = np.max(list(counts.values())) / 2

for key, value in counts.items():
    if value > threshold:
        counts_keep[key] = value

print(counts_keep)

{'00000000': 58, '01000000': 41, '11000000': 42, '10000000': 40}

4단계: 후처리 및 원하는 고전 형식으로 결과 반환

정수 인수분해

지금까지 위상 추정 Circuit을 사용하여 위수 찾기 문제를 구현하는 방법을 논의했습니다. 이제 위수 찾기 문제를 정수 인수분해와 연결하여 Shor 알고리즘을 완성합니다. 이 부분의 알고리즘은 고전적임을 참고하세요.

$N = 15$, $a = 2$ 예시를 사용하여 이를 설명합니다. 측정된 위상이 $k / r$ (여기서 $a^r \; (\textrm{mod} \; N) = 1$ 이고 $k$는 $0$과 $r - 1$ 사이의 임의 정수)임을 기억하세요. 이 방정식에서 $(a^r - 1) \; (\textrm{mod} \; 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$의 진인수일 높은 확률이 있습니다.

알고리즘의 일부 실행은 통계적으로 실패하므로, $N$의 인수를 하나 이상 찾을 때까지 이 알고리즘을 반복합니다.

아래 셀은 $N=15$의 인수를 하나 이상 찾을 때까지 알고리즘을 반복합니다. 각 반복에서 위상 및 그에 대응하는 인수를 추측하기 위해 위의 하드웨어 실행 결과를 사용합니다.

a = 2
N = 15

FACTOR_FOUND = False
num_attempt = 0

while not FACTOR_FOUND:
    print(f"\nATTEMPT {num_attempt}:")
    # Here, we get the bitstring by iterating over outcomes
    # of a previous hardware run with multiple shots.
    # Instead, we can also perform a single-shot measurement
    # here in the loop.
    bitstring = list(counts_keep.keys())[num_attempt]
    num_attempt += 1
    # Find the phase from measurement
    decimal = int(bitstring, 2)
    phase = decimal / (2**num_control)  # phase = k / r
    print(f"Phase: theta = {phase}")

    # Guess the order from phase
    frac = Fraction(phase).limit_denominator(N)
    r = frac.denominator  # order = r
    print(f"Order of {a} modulo {N} estimated as: r = {r}")

    if phase != 0:
        # Guesses for factors are gcd(a^{r / 2} ± 1, 15)
        if r % 2 == 0:
            x = pow(a, r // 2, N) - 1
            d = gcd(x, N)
            if d > 1:
                FACTOR_FOUND = True
                print(f"*** Non-trivial factor found: {x} ***")

ATTEMPT 0:
Phase: theta = 0.0
Order of 2 modulo 15 estimated as: r = 1

ATTEMPT 1:
Phase: theta = 0.25
Order of 2 modulo 15 estimated as: r = 4
*** Non-trivial factor found: 3 ***

토론

관련 연구

이 섹션에서는 실제 하드웨어에서 Shor 알고리즘을 시연한 다른 주요 연구들을 논의합니다.

IBM®의 선구적인 연구 [3]는 7 Qubit 핵자기공명(NMR) 양자 컴퓨터를 사용하여 숫자 15를 소인수 3과 5로 처음으로 인수분해함으로써 Shor 알고리즘을 최초로 시연했습니다. 또 다른 실험 [4]은 광자 Qubit을 사용하여 15를 인수분해했습니다. 단일 Qubit을 여러 번 재활용하고 작업 레지스터를 고차원 상태로 인코딩함으로써, 연구자들은 표준 프로토콜에 비해 필요한 Qubit 수를 3분의 1로 줄였으며, 2-광자 컴파일된 알고리즘을 활용했습니다. Shor 알고리즘 시연에서 중요한 논문 [5]는 Kitaev의 반복 위상 추정 [8] 기법을 사용하여 알고리즘의 Qubit 요구 사항을 줄였습니다. 저자들은 7개의 제어 Qubit과 4개의 캐시 Qubit을 사용했으며, 모듈식 곱셈기 구현도 포함되었습니다. 그러나 이 구현은 순방향 연산을 통한 중간 Circuit 측정과 초기화 연산을 통한 Qubit 재활용이 필요합니다. 이 시연은 이온 트랩 양자 컴퓨터에서 수행되었습니다.

최근 연구 [6]는 IBM Quantum® 하드웨어에서 15, 21, 35를 인수분해하는 데 집중했습니다. 이전 연구와 유사하게, 연구자들은 물리적 Qubit과 Gate 수를 최소화하기 위해 Kitaev가 제안한 반고전적 양자 푸리에 변환을 사용하는 알고리즘의 컴파일된 버전을 활용했습니다. 가장 최근의 연구 [7]도 정수 21을 인수분해하는 개념 증명 시연을 수행했습니다. 이 시연에서도 컴파일된 버전의 양자 위상 추정 루틴을 사용했으며, [4]의 이전 시연을 발전시켰습니다. 저자들은 잔류 위상 이동이 있는 근사 Toffoli Gate 구성을 사용함으로써 이 연구를 발전시켰습니다. 알고리즘은 단 5개의 Qubit만을 사용하여 IBM 양자 프로세서에서 구현되었으며, 제어 Qubit과 레지스터 Qubit 사이의 얽힘이 성공적으로 검증되었습니다.

알고리즘의 확장성

RSA 암호화는 일반적으로 2048~4096비트 규모의 키 크기를 사용합니다. Shor 알고리즘으로 2048비트 수를 인수분해하려면 오류 수정 오버헤드를 포함하여 수백만 개의 Qubit과 수십억 수준의 Circuit 깊이가 필요한 양자 Circuit이 요구되며, 이는 현재 양자 하드웨어의 실행 한계를 훨씬 초과합니다. 따라서 Shor 알고리즘이 현대 암호 시스템을 실질적으로 해독하려면 최적화된 Circuit 구성 방법이나 강력한 양자 오류 수정이 필요합니다. Shor 알고리즘의 자원 추정에 대한 자세한 논의는 [9]를 참고하세요.

도전 과제

튜토리얼을 완료하신 것을 축하합니다! 이제 이해도를 테스트해 볼 좋은 기회입니다. 21을 인수분해하는 Circuit을 구성해 보실 수 있을까요? $a$ 값은 직접 선택하시면 됩니다. Qubit 수를 결정하기 위해 알고리즘의 비트 정밀도를 정해야 하고, 모듈식 지수 연산자 $M_a$를 설계해야 합니다. 직접 시도해 보신 후, [6]의 Fig. 9와 [7]의 Fig. 2에 나와 있는 방법론을 읽어보실 것을 권장합니다.

def M_a_mod21():
    """
    M_a (mod 21)
    """

    # Your code here
    pass

참고문헌

1. Shor, Peter W. "Polynomial-time algorithms for prime factorization and discrete logarithms on a quantum computer." SIAM review 41.2 (1999): 303-332.
2. IBM Quantum Fundamentals of Quantum Algorithms course by Dr. John Watrous.
3. Vandersypen, Lieven MK, et al. "Experimental realization of Shor's quantum factoring algorithm using nuclear magnetic resonance." Nature 414.6866 (2001): 883-887.
4. Martin-Lopez, Enrique, et al. "Experimental realization of Shor's quantum factoring algorithm using qubit recycling." Nature photonics 6.11 (2012): 773-776.
5. Monz, Thomas, et al. "Realization of a scalable Shor algorithm." Science 351.6277 (2016): 1068-1070.
6. Amico, Mirko, Zain H. Saleem, and Muir Kumph. "Experimental study of Shor's factoring algorithm using the IBM Q Experience." Physical Review A 100.1 (2019): 012305.
7. Skosana, Unathi, and Mark Tame. "Demonstration of Shor's factoring algorithm for N=21 on IBM quantum processors." Scientific reports 11.1 (2021): 16599.
8. Kitaev, A. Yu. "Quantum measurements and the Abelian stabilizer problem." arXiv preprint quant-ph/9511026 (1995).
9. Gidney, Craig, and Martin Ekerå. "How to factor 2048 bit RSA integers in 8 hours using 20 million noisy qubits." Quantum 5 (2021): 433.