Algorytm Shora

Szacowane zużycie zasobów: trzy sekundy na procesorze Eagle r3 (UWAGA: to tylko szacunek. Twój czas wykonania może się różnić.)

Algorytm Shora, opracowany przez Petera Shora w 1994 roku, to przełomowy algorytm kwantowy do faktoryzacji liczb całkowitych w czasie wielomianowym. Jego znaczenie polega na zdolności do faktoryzacji dużych liczb całkowitych wykładniczo szybciej niż jakikolwiek znany algorytm klasyczny, co zagraża bezpieczeństwu powszechnie stosowanych systemów kryptograficznych, takich jak RSA, których bezpieczeństwo opiera się na trudności faktoryzacji dużych liczb. Poprzez efektywne rozwiązanie tego problemu na wystarczająco wydajnym komputerze kwantowym, algorytm Shora mógłby zrewolucjonizować dziedziny takie jak kryptografia, cyberbezpieczeństwo i matematyka obliczeniowa, podkreślając transformacyjną moc obliczeń kwantowych.

Ten samouczek skupia się na demonstracji algorytmu Shora poprzez faktoryzację liczby 15 na komputerze kwantowym.

Najpierw definiujemy problem wyznaczania rzędu i konstruujemy odpowiednie Circuit z protokołu kwantowej estymacji fazy. Następnie uruchamiamy Circuit wyznaczające rząd na prawdziwym sprzęcie, używając Circuit o najkrótszej głębokości, jakie jesteśmy w stanie skompilować. Ostatnia sekcja kończy algorytm Shora, łącząc problem wyznaczania rzędu z faktoryzacją liczb całkowitych.

Kończymy samouczek dyskusją na temat innych demonstracji algorytmu Shora na prawdziwym sprzęcie, skupiając się zarówno na implementacjach ogólnych, jak i tych dostosowanych do faktoryzacji konkretnych liczb, takich jak 15 i 21. Uwaga: Ten samouczek skupia się bardziej na implementacji i demonstracji Circuit dotyczących algorytmu Shora. Aby uzyskać szczegółowy materiał edukacyjny, zapoznaj się z kursem Podstawy algorytmów kwantowych autorstwa dr. Johna Watrouca oraz artykułami w sekcji Odnośniki.

Wymagania

Przed rozpoczęciem tego samouczka upewnij się, że masz zainstalowane następujące elementy:

• Qiskit SDK w wersji 2.0 lub nowszej, z obsługą wizualizacji
• Qiskit Runtime w wersji 0.40 lub nowszej (pip install qiskit-ibm-runtime)

Konfiguracja

# 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

Krok 1: Odwzorowanie klasycznych danych wejściowych na problem kwantowy

Tło

Algorytm Shora do faktoryzacji liczb całkowitych wykorzystuje pośredni problem znany jako problem wyznaczania rzędu. W tej sekcji pokazujemy, jak rozwiązać problem wyznaczania rzędu za pomocą kwantowej estymacji fazy.

Problem estymacji fazy

W problemie estymacji fazy dany jest stan kwantowy $\ket{\psi}$ złożony z $n$ Qubit, wraz z kwantowym Circuit działającym na $n$ Qubit. Zakłada się, że $\ket{\psi}$ jest wektorem własnym macierzy unitarnej $U$ opisującej działanie Circuit, a naszym celem jest obliczenie lub przybliżenie wartości własnej $\lambda = e^{2 \pi i \theta}$, do której $\ket{\psi}$ należy. Innymi słowy, Circuit powinien wypisać przybliżenie liczby $\theta \in [0, 1)$ spełniającej $U \ket{\psi}= e^{2 \pi i \theta} \ket{\psi}.$ Celem Circuit estymacji fazy jest przybliżenie $\theta$ na $m$ bitach. Matematycznie rzecz ujmując, chcemy znaleźć $y$ takie, że $\theta \approx y / 2^m$, gdzie $y \in {0, 1, 2, \dots, 2^{m-1}}$. Poniższy obraz przedstawia Circuit kwantowy szacujący $y$ na $m$ bitach poprzez pomiar na $m$ Qubit. W powyższym Circuit górne $m$ Qubit jest inicjowanych w stanie $\ket{0^m}$, a dolne $n$ Qubit — w stanie $\ket{\psi}$, który z założenia jest wektorem własnym $U$. Pierwszym składnikiem Circuit estymacji fazy są operacje sterowanej-unitarności odpowiedzialne za wykonanie phase kickback do odpowiadającego im Qubit kontrolnego. Te sterowane operacje unitarne są potęgowane zgodnie z pozycją Qubit kontrolnego, od najmniej znaczącego bitu do najbardziej znaczącego. Ponieważ $\ket{\psi}$ jest wektorem własnym $U$, stan dolnych $n$ Qubit nie jest przez tę operację zmieniony, ale informacja o fazie wartości własnej propaguje do górnych $m$ Qubit. Okazuje się, że po operacji phase kickback za pomocą sterowanych-unitarności, wszystkie możliwe stany górnych $m$ Qubit są wzajemnie ortonormalne dla każdego wektora własnego $\ket{\psi}$ unitarności $U$. Dlatego stany te są doskonale rozróżnialne i możemy obrócić bazę, którą tworzą, z powrotem do bazy obliczeniowej, aby dokonać pomiaru. Analiza matematyczna pokazuje, że ta macierz obrotu odpowiada odwrotnej kwantowej transformacie Fouriera (QFT) w przestrzeni Hilberta wymiaru $2^m$. Intuicja za tym stoi taka, że periodyczna struktura operatorów modularnego potęgowania jest zakodowana w stanie kwantowym, a QFT przekształca tę periodyczność na mierzalne piki w dziedzinie częstotliwości.

Aby głębiej zrozumieć, dlaczego Circuit QFT jest stosowany w algorytmie Shora, odsyłamy czytelnika do kursu Podstawy algorytmów kwantowych. Jesteśmy teraz gotowi, aby użyć Circuit estymacji fazy do wyznaczania rzędu.

Problem wyznaczania rzędu

Aby zdefiniować problem wyznaczania rzędu, zaczynamy od kilku pojęć z teorii liczb. Przede wszystkim, dla dowolnej danej dodatniej liczby całkowitej $N$, definiujemy zbiór $\mathbb{Z}_N$ jako $\mathbb{Z}_N = \{0, 1, 2, \dots, N-1\}.$ Wszystkie operacje arytmetyczne w $\mathbb{Z}_N$ wykonywane są modulo $N$. W szczególności, wszystkie elementy $a \in \mathbb{Z}_n$ będące wzajemnie pierwsze z $N$ są szczególne i tworzą $\mathbb{Z}^*_N$ jako $\mathbb{Z}^*_N = \{ a \in \mathbb{Z}_N : \mathrm{gcd}(a, N)=1 \}.$ Dla elementu $a \in \mathbb{Z}^*_N$, najmniejsza dodatnia liczba całkowita $r$ taka, że $a^r \equiv 1 \; (\mathrm{mod} \; N)$ jest definiowana jako rząd $a$ modulo $N$. Jak zobaczymy później, znalezienie rzędu $a \in \mathbb{Z}^*_N$ pozwoli nam faktoryzować $N$. Aby skonstruować Circuit wyznaczający rząd na podstawie Circuit estymacji fazy, potrzebujemy dwóch rzeczy. Po pierwsze, musimy zdefiniować unitarność $U$, która pozwoli nam znaleźć rząd $r$, a po drugie, musimy zdefiniować wektor własny $\ket{\psi}$ unitarności $U$, aby przygotować stan początkowy Circuit estymacji fazy.

Aby powiązać problem wyznaczania rzędu z estymacją fazy, rozważamy operację zdefiniowaną na układzie, którego stany klasyczne odpowiadają $\mathbb{Z}_N$, gdzie mnożymy przez stały element $a \in \mathbb{Z}^*_N$. W szczególności definiujemy ten operator mnożenia $M_a$ tak, że $M_a \ket{x} = \ket{ax \; (\mathrm{mod} \; N)}$ dla każdego $x \in \mathbb{Z}_N$. Należy zauważyć, że niejawnie przyjmujemy, że iloczyn modulo $N$ jest obliczany wewnątrz ket po prawej stronie równania. Analiza matematyczna pokazuje, że $M_a$ jest operatorem unitarnym. Co więcej, okazuje się, że $M_a$ ma pary wektor własny–wartość własna, które pozwalają nam powiązać rząd $r$ elementu $a$ z problemem estymacji fazy. Konkretnie, dla dowolnego wyboru $j \in \{0, \dots, r-1\}$ mamy, że $\ket{\psi_j} = \frac{1}{\sqrt{r}} \sum^{r-1}_{k=0} \omega^{-jk}_{r} \ket{a^k}$ jest wektorem własnym $M_a$, któremu odpowiada wartość własna $\omega^{j}_{r}$, gdzie $\omega^{j}_{r} = e^{2 \pi i \frac{j}{r}}.$ Przez obserwację widzimy, że wygodną parą wektor własny/wartość własna jest stan $\ket{\psi_1}$ z $\omega^{1}_{r} = e^{2 \pi i \frac{1}{r}}$. Dlatego gdybyśmy mogli znaleźć wektor własny $\ket{\psi_1}$, moglibyśmy oszacować fazę $\theta=1/r$ naszym Circuit kwantowym i tym samym uzyskać przybliżenie rzędu $r$. Jednak nie jest to łatwe i musimy rozważyć alternatywę.

Zastanówmy się, co dałby Circuit, gdybyśmy przygotowali stan obliczeniowy $\ket{1}$ jako stan początkowy. Nie jest to stan własny $M_a$, ale jest to równomierna superpozycja stanów własnych, które właśnie opisaliśmy. Innymi słowy, zachodzi następująca zależność: $\ket{1} = \frac{1}{\sqrt{r}} \sum^{r-1}_{k=0} \ket{\psi_k}$ Konsekwencją powyższego równania jest to, że jeśli ustawimy stan początkowy na $\ket{1}$, otrzymamy dokładnie taki sam wynik pomiaru, jak gdybyśmy wybrali $k \in \{ 0, \dots, r-1\}$ jednostajnie losowo i użyli $\ket{\psi_k}$ jako wektora własnego w Circuit estymacji fazy. Innymi słowy, pomiar górnych $m$ Qubit daje przybliżenie $y / 2^m$ wartości $k / r$, gdzie $k \in \{ 0, \dots, r-1\}$ jest wybrany jednostajnie losowo. Pozwala nam to poznać $r$ z dużą pewnością po kilku niezależnych uruchomieniach, co było naszym celem.

Operatory modularnego potęgowania

Do tej pory powiązaliśmy problem estymacji fazy z problemem wyznaczania rzędu, definiując $U = M_a$ i $\ket{\psi} = \ket{1}$ w naszym Circuit kwantowym. Dlatego ostatnim brakującym składnikiem jest znalezienie efektywnego sposobu definiowania potęg modularnych $M_a$ jako $M_a^k$ dla $k = 1, 2, 4, \dots, 2^{m-1}$. Aby wykonać to obliczenie, odkrywamy, że dla dowolnej potęgi $k$ możemy stworzyć Circuit dla $M_a^k$ nie przez $k$-krotną iterację Circuit dla $M_a$, lecz przez obliczenie $b = a^k \; \mathrm{mod} \; N$, a następnie użycie Circuit dla $M_b$. Ponieważ potrzebujemy tylko potęg będących samymi potęgami dwójki, możemy to zrobić klasycznie wydajnie za pomocą iteracyjnego potęgowania przez podnoszenie do kwadratu.

Krok 2: Optymalizacja problemu do wykonania na sprzęcie kwantowym

Konkretny przykład z $N = 15$ i $a=2$

Możemy tu zatrzymać się, aby omówić konkretny przykład i skonstruować Circuit wyznaczający rząd dla $N=15$. Należy zauważyć, że możliwe nietrywalne $a \in \mathbb{Z}_N^*$ dla $N=15$ to $a \in \{2, 4, 7, 8, 11, 13, 14 \}$. W tym przykładzie wybieramy $a=2$. Skonstruujemy operator $M_2$ oraz operatory modularnego potęgowania $M_2^k$. Działanie $M_2$ na stanach bazy obliczeniowej jest następujące. $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}$ Przez obserwację widzimy, że stany bazy są przestawiane, więc mamy macierz permutacji. Możemy skonstruować tę operację na czterech Qubit za pomocą bramek swap. Poniżej konstruujemy operacje $M_2$ i sterowanego-$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)

Gate działające na więcej niż dwóch Qubit zostaną dalej rozłożone na Gate dwuqubitowe.

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

Teraz musimy skonstruować operatory modularnego potęgowania. Aby uzyskać wystarczającą precyzję w estymacji fazy, użyjemy ośmiu Qubit do pomiaru estymacji. Dlatego musimy skonstruować $M_b$ z $b = a^{2^k} \; (\mathrm{mod} \; N)$ dla każdego $k = 0, 1, \dots, 7$.

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]

Jak widzimy z listy wartości $b$, oprócz $M_2$, który wcześniej skonstruowaliśmy, musimy również zbudować $M_4$ i $M_1$. Należy zauważyć, że $M_1$ działa trywialnie na stanach bazy obliczeniowej, więc jest po prostu operatorem tożsamości.

$M_4$ działa na stanach bazy obliczeniowej w następujący sposób. $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}$

Zatem tę permutację można skonstruować za pomocą następującej operacji swap.

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 działające na więcej niż dwóch Qubit zostaną dalej rozłożone na Gate dwuqubitowe.

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

Widzieliśmy, że operatory $M_b$ dla danego $b \in \mathbb{Z}^*_N$ są operacjami permutacji. Ze względu na stosunkowo mały rozmiar problemu permutacji, który tutaj rozważamy — ponieważ $N=15$ wymaga tylko czterech Qubit — byliśmy w stanie zsyntetyzować te operacje bezpośrednio za pomocą bramek SWAP przez inspekcję. Ogólnie rzecz biorąc, może to nie być skalowalne podejście. Zamiast tego może być konieczne jawne skonstruowanie macierzy permutacji i użycie klasy UnitaryGate Qiskit oraz metod transpilacji do syntezy tej macierzy permutacji. Może to jednak prowadzić do znacznie głębszych Circuit. Przykład poniżej.

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})

Porównajmy te liczby z głębokością skompilowanego Circuit naszej ręcznej implementacji Gate $M_2$.

# 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})

Jak widzimy, podejście z macierzą permutacji dało znacznie głębszy Circuit nawet dla pojedynczego Gate $M_2$ w porównaniu z naszą ręczną implementacją. Dlatego będziemy kontynuować z naszą poprzednią implementacją operacji $M_b$. Jesteśmy teraz gotowi skonstruować pełny Circuit wyznaczający rząd przy użyciu wcześniej zdefiniowanych sterowanych operatorów modularnego potęgowania. W poniższym kodzie importujemy również Circuit QFT z biblioteki Circuit Qiskit, który używa bramek Hadamarda na każdym Qubit, serii sterowanych Gate U1 (lub Z, w zależności od fazy) i warstwy bramek swap.

# 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)

Należy zauważyć, że pominęliśmy sterowane operacje modularnego potęgowania z pozostałych Qubit kontrolnych, ponieważ $M_1$ jest operatorem tożsamości. Należy zauważyć, że w dalszej części tego samouczka uruchomimy ten Circuit na Backend ibm_marrakesh. Aby to zrobić, kompilujemy Circuit zgodnie z tym konkretnym Backend i raportujemy głębokość Circuit oraz liczbę 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})

Krok 3: Wykonanie z użyciem prymitywów Qiskit

Najpierw omówimy, co teoretycznie uzyskalibyśmy, gdybyśmy uruchomili ten Circuit na idealnym symulatorze. Poniżej przedstawiamy zestaw wyników symulacji powyższego Circuit przy użyciu 1024 próbek. Jak widać, otrzymujemy w przybliżeniu jednostajny rozkład na czterech ciągach bitów dla kubitów sterujących.

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

plot_histogram(counts)

Mierząc kubity sterujące, uzyskujemy ośmiobitową estymację fazy operatora $M_a$. Możemy przekształcić tę reprezentację binarną na dziesiętną, aby znaleźć zmierzoną fazę. Jak widać z powyższego histogramu, zmierzono cztery różne ciągi bitów, a każdy z nich odpowiada wartości fazy w następujący sposób.

# 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

Przypomnijmy, że każda zmierzona faza odpowiada $\theta = k / r$, gdzie $k$ jest losowane jednostajnie z $\{0, 1, \dots, r-1 \}$. Możemy zatem użyć algorytmu ułamków łańcuchowych, aby spróbować znaleźć $k$ i rząd $r$. Python ma tę funkcjonalność wbudowaną. Możemy użyć modułu fractions, aby zamienić liczbę zmiennoprzecinkową na obiekt Fraction, na przykład:

Fraction(0.666)

Fraction(5998794703657501, 9007199254740992)

Ponieważ daje to ułamki, które zwracają dokładny wynik (w tym przypadku 0.6660000...), może to dawać nieprzyjemne rezultaty, jak powyżej. Możemy użyć metody .limit_denominator(), aby uzyskać ułamek najbardziej zbliżony do naszej liczby zmiennoprzecinkowej, z mianownikiem poniżej pewnej wartości:

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

Fraction(2, 3)

To jest o wiele ładniejszy wynik. Rząd (r) musi być mniejszy od N, więc ustawimy maksymalny mianownik na 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

Widzimy, że dwie z zmierzonych wartości własnych dały nam poprawny wynik: $r=4$, a algorytm Shora do wyznaczania rzędu może zawieść. Te złe wyniki wynikają z tego, że $k = 0$ lub że $k$ i $r$ nie są wzajemnie pierwsze — i zamiast $r$ otrzymujemy czynnik $r$. Najprostszym rozwiązaniem jest po prostu powtarzanie eksperymentu, aż uzyskamy satysfakcjonujący wynik dla $r$. Do tej pory zaimplementowaliśmy problem wyznaczania rzędu dla $N=15$ z $a=2$, używając Circuit estymacji fazy na symulatorze. Ostatnim krokiem algorytmu Shora będzie powiązanie problemu wyznaczania rzędu z problemem faktoryzacji liczb całkowitych. Ta ostatnia część algorytmu jest czysto klasyczna i może być rozwiązana na klasycznym komputerze po uzyskaniu pomiarów fazy z komputera kwantowego. Dlatego odłożymy tę ostatnią część algorytmu na po demonstracji uruchomienia Circuit wyznaczania rzędu na rzeczywistym sprzęcie.

Uruchomienia na sprzęcie

Teraz możemy uruchomić Circuit wyznaczania rzędu, który wcześniej skompilowaliśmy dla ibm_marrakesh. Sięgamy tu po dynamiczne odsprzęganie (DD) do tłumienia błędów oraz twirling bramek do mitygacji błędów. DD polega na stosowaniu sekwencji precyzyjnie zsynchronizowanych impulsów sterujących do urządzenia kwantowego, efektywnie uśredniając niepożądane interakcje ze środowiskiem i dekoherencję. Twirling bramek natomiast losowo przekształca określone bramki kwantowe, zamieniając błędy koherentne w błędy Pauliego, które kumulują się liniowo, a nie kwadratowo. Obie techniki są często łączone, aby poprawić koherencję i wierność obliczeń kwantowych.

# 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))

Jak widać, uzyskaliśmy te same ciągi bitów z najwyższymi zliczeniami. Ponieważ sprzęt kwantowy ma szumy, następuje pewien wyciek do innych ciągów bitów, który możemy statystycznie odfiltrować.

# 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}

Krok 4: Post-processing i zwrócenie wyniku w pożądanym formacie klasycznym

Faktoryzacja liczb całkowitych

Do tej pory omówiliśmy, jak możemy zaimplementować problem wyznaczania rzędu przy użyciu Circuit estymacji fazy. Teraz łączymy problem wyznaczania rzędu z faktoryzacją liczb całkowitych, co uzupełnia algorytm Shora. Należy zauważyć, że ta część algorytmu jest klasyczna. Demonstrujemy to teraz na naszym przykładzie $N = 15$ i $a = 2$. Przypomnijmy, że zmierzona faza to $k / r$, gdzie $a^r \; (\textrm{mod} \; N) = 1$ a $k$ jest losową liczbą całkowitą między $0$ a $r - 1$. Z tego równania mamy $(a^r - 1) \; (\textrm{mod} \; N) = 0,$ co oznacza, że $N$ musi dzielić $a^r-1$. Jeśli $r$ jest parzyste, możemy napisać $a^r -1 = (a^{r/2}-1)(a^{r/2}+1).$ Jeśli $r$ nie jest parzyste, nie możemy iść dalej i musimy spróbować ponownie z inną wartością $a$; w przeciwnym razie istnieje duże prawdopodobieństwo, że największy wspólny dzielnik $N$ i $a^{r/2}-1$ lub $a^{r/2}+1$ jest właściwym czynnikiem $N$.

Ponieważ niektóre uruchomienia algorytmu statystycznie zawiodą, będziemy powtarzać ten algorytm, aż zostanie znaleziony co najmniej jeden czynnik $N$. Poniższa komórka powtarza algorytm, aż zostanie znaleziony co najmniej jeden czynnik $N=15$. Użyjemy wyników powyższego uruchomienia na sprzęcie, aby w każdej iteracji zgadnąć fazę i odpowiadający jej czynnik.

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 ***

Dyskusja

Powiązane prace

W tej sekcji omawiamy inne przełomowe prace, które zademonstrowały algorytm Shora na rzeczywistym sprzęcie.

Przełomowa praca [3] firmy IBM® zademonstrowała algorytm Shora po raz pierwszy, faktoryzując liczbę 15 na jej pierwsze czynniki 3 i 5 przy użyciu siedmiokubitowego komputera kwantowego z magnetycznym rezonansem jądrowym (NMR). Inny eksperyment [4] sfaktoryzował 15 przy użyciu kubitów fotonicznych. Stosując jeden Qubit wielokrotnie używany i kodując rejestr roboczy w stanach wyższego wymiaru, badacze zredukowali wymaganą liczbę kubitów do jednej trzeciej standardowego protokołu, korzystając ze skompilowanego algorytmu dwufotonowego. Znaczącą pracą w demonstracji algorytmu Shora jest [5], która używa iteracyjnej estymacji fazy Kitajewa [8] do zmniejszenia wymagań qubitowych algorytmu. Autorzy użyli siedmiu kubitów sterujących i czterech kubitów podręcznych, wraz z implementacją modularnych mnożników. Ta implementacja wymaga jednak śródukładowych pomiarów z operacjami feed-forward oraz recyklingu kubitów z operacjami resetowania. Demonstracja ta została przeprowadzona na komputerze kwantowym z pułapką jonową.

Nowsze prace [6] skupiły się na faktoryzacji 15, 21 i 35 na sprzęcie IBM Quantum®. Podobnie jak poprzednie prace, badacze użyli skompilowanej wersji algorytmu, która zastosowała półklasyczną kwantową transformatę Fouriera zaproponowaną przez Kitajewa, aby zminimalizować liczbę fizycznych kubitów i bramek. Najnowsza praca [7] przeprowadziła również demonstrację koncepcji faktoryzacji liczby całkowitej 21. Ta demonstracja również obejmowała użycie skompilowanej wersji procedury kwantowej estymacji fazy i opierała się na poprzedniej demonstracji z pracy [4]. Autorzy poszli dalej, używając konfiguracji przybliżonych bramek Toffoli z resztkowymi przesunięciami fazy. Algorytm zaimplementowano na procesorach kwantowych IBM używając jedynie pięciu kubitów, i pomyślnie zweryfikowano obecność splątania między kubitami sterującymi a kubitami rejestrów.

Skalowalność algorytmu

Należy zauważyć, że szyfrowanie RSA zazwyczaj wiąże się z rozmiarami kluczy rzędu 2048 do 4096 bitów. Próba faktoryzacji 2048-bitowej liczby algorytmem Shora spowoduje Circuit kwantowy z milionami kubitów, włączając narzut korekcji błędów, i głębokością Circuit rzędu miliarda, co przekracza możliwości obecnego sprzętu kwantowego. Dlatego algorytm Shora będzie wymagał zoptymalizowanych metod konstruowania Circuit lub solidnej kwantowej korekcji błędów, aby być praktycznie wykonalnym w łamaniu nowoczesnych systemów kryptograficznych. Odsyłamy do [9] po bardziej szczegółową dyskusję na temat estymacji zasobów dla algorytmu Shora.

Zadanie

Gratulacje z okazji ukończenia tutoriala! To dobry moment, aby sprawdzić swoje rozumienie. Czy możesz spróbować skonstruować Circuit do faktoryzacji 21? Możesz wybrać $a$ według własnego uznania. Będziesz musiał zdecydować o dokładności bitowej algorytmu, aby wybrać liczbę kubitów, oraz zaprojektować operatory modularnego potęgowania $M_a$. Zachęcamy do samodzielnej próby, a następnie zapoznania się z metodologiami przedstawionymi na rys. 9 w [6] i rys. 2 w [7].

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

    # Your code here
    pass

Odnośniki

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.