Shorův algoritmus

Odhadované využití: Tři sekundy na procesoru Eagle r3 (POZNÁMKA: Jedná se pouze o odhad. Skutečná doba běhu se může lišit.)

Shorův algoritmus, vyvinutý Peterem Shorem v roce 1994, je průlomový kvantový algoritmus pro faktorizaci celých čísel v polynomiálním čase. Jeho význam spočívá ve schopnosti faktorizovat velká celá čísla exponenciálně rychleji než jakýkoli známý klasický algoritmus, čímž ohrožuje bezpečnost široce používaných kryptografických systémů, jako je RSA, které se spoléhají na obtížnost faktorizace velkých čísel. Efektivním řešením tohoto problému na dostatečně výkonném kvantovém počítači by Shorův algoritmus mohl revolucionizovat oblasti jako kryptografie, kybernetická bezpečnost a výpočetní matematika, což podtrhuje transformativní sílu kvantové výpočetní techniky.

Tento tutoriál se zaměřuje na demonstraci Shorova algoritmu faktorizací čísla 15 na kvantovém počítači.

Nejprve definujeme problém hledání řádu a konstruujeme odpovídající Circuit z protokolu kvantového odhadování fáze. Dále spustíme Circuit pro hledání řádu na reálném hardwaru s využitím co nejkratších Circuit, které dokážeme transpilovat. Poslední část dokončuje Shorův algoritmus propojením problému hledání řádu s faktorizací celých čísel.

Tutoriál zakončíme diskusí o dalších demonstracích Shorova algoritmu na reálném hardwaru se zaměřením jak na obecné implementace, tak na ty přizpůsobené faktorizaci konkrétních čísel, jako jsou 15 a 21. Poznámka: Tento tutoriál se více zaměřuje na implementaci a demonstraci Circuit týkajících se Shorova algoritmu. Pro podrobný vzdělávací zdroj k tomuto tématu se prosím podívej na kurz Fundamentals of quantum algorithms od Dr. Johna Watrouse a na články v sekci Odkazy.

Požadavky

Před zahájením tohoto tutoriálu se ujisti, že máš nainstalováno:

• Qiskit SDK v2.0 nebo novější, s podporou vizualizace
• Qiskit Runtime v0.40 nebo novější (pip install qiskit-ibm-runtime)

Nastavení

# 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: Mapování klasických vstupů na kvantový problém

Pozadí

Shorův algoritmus pro faktorizaci celých čísel využívá pomocný problém označovaný jako hledání řádu. V této části ukážeme, jak řešit problém hledání řádu pomocí kvantového odhadování fáze.

Problém odhadování fáze

V problému odhadování fáze máme dán kvantový stav $\ket{\psi}$ o $n$ Qubitech spolu s unitárním kvantovým Circuit, který působí na $n$ Qubitech. Je zaručeno, že $\ket{\psi}$ je vlastním vektorem unitární matice $U$ popisující akci Circuit, a naším cílem je vypočítat nebo přiblížit vlastní hodnotu $\lambda = e^{2 \pi i \theta}$, které $\ket{\psi}$ odpovídá. Jinými slovy, Circuit by měl vydat aproximaci čísla $\theta \in [0, 1)$ splňující $U \ket{\psi}= e^{2 \pi i \theta} \ket{\psi}.$ Cílem Circuit pro odhadování fáze je přiblížit $\theta$ na $m$ bitů. Matematicky řečeno, chceme najít $y$ tak, aby $\theta \approx y / 2^m$, kde $y \in {0, 1, 2, \dots, 2^{m-1}}$. Následující obrázek ukazuje kvantový Circuit, který odhaduje $y$ na $m$ bitů provedením měření na $m$ Qubitech. Ve výše uvedeném Circuit je horních $m$ Qubitů inicializováno ve stavu $\ket{0^m}$ a spodních $n$ Qubitů ve stavu $\ket{\psi}$, který je zaručen jako vlastní vektor $U$. Prvním prvkem Circuit pro odhadování fáze jsou kontrolované unitární operace, které jsou zodpovědné za provádění phase kickback na příslušný řídicí Qubit. Tyto kontrolované unitární operace jsou umocňovány podle pozice řídicího Qubitu, od nejméně významného bitu po nejvýznamnější. Protože $\ket{\psi}$ je vlastním vektorem $U$, stav spodních $n$ Qubitů touto operací není ovlivněn, ale fázová informace vlastní hodnoty se propaguje na horních $m$ Qubitů. Ukázalo se, že po operaci phase kickback prostřednictvím kontrolovaných unitárních operací jsou všechny možné stavy horních $m$ Qubitů navzájem ortonormální pro každý vlastní vektor $\ket{\psi}$ unitárního $U$. Tyto stavy jsou tedy dokonale rozlišitelné a můžeme otočit bázi, kterou tvoří, zpět do výpočetní báze pro měření. Matematická analýza ukazuje, že tato rotační matice odpovídá inverzní kvantové Fourierově transformaci (QFT) v $2^m$-dimenzionálním Hilbertově prostoru. Intuice za tím je, že periodická struktura operátorů modulárního umocňování je zakódována v kvantovém stavu a QFT převádí tuto periodicitu na měřitelné vrcholy ve frekvenční doméně.

Pro hlubší pochopení toho, proč je Circuit QFT v Shorově algoritmu využíván, odkazujeme čtenáře na kurz Fundamentals of quantum algorithms. Nyní jsme připraveni použít Circuit pro odhadování fáze k hledání řádu.

Problém hledání řádu

Pro definování problému hledání řádu začneme s několika koncepty z teorie čísel. Nejprve pro libovolné kladné celé číslo $N$ definujeme množinu $\mathbb{Z}_N$ jako $\mathbb{Z}_N = \{0, 1, 2, \dots, N-1\}.$ Veškeré aritmetické operace v $\mathbb{Z}_N$ se provádějí modulo $N$. Zejména všechny prvky $a \in \mathbb{Z}_n$, které jsou nesoudělné s $N$, jsou speciální a tvoří $\mathbb{Z}^*_N$ jako $\mathbb{Z}^*_N = \{ a \in \mathbb{Z}_N : \mathrm{gcd}(a, N)=1 \}.$ Pro prvek $a \in \mathbb{Z}^*_N$ je nejmenší kladné celé číslo $r$ takové, že $a^r \equiv 1 \; (\mathrm{mod} \; N)$ definováno jako řád $a$ modulo $N$. Jak uvidíme později, nalezení řádu $a \in \mathbb{Z}^*_N$ nám umožní faktorizovat $N$. Pro sestavení Circuit pro hledání řádu z Circuit pro odhadování fáze potřebujeme dvě věci. Zaprvé musíme definovat unitární operátor $U$, který nám umožní najít řád $r$, a zadruhé musíme definovat vlastní vektor $\ket{\psi}$ operátoru $U$ pro přípravu počátečního stavu Circuit pro odhadování fáze.

Pro propojení problému hledání řádu s odhadováním fáze uvažujeme operaci definovanou na systému, jehož klasické stavy odpovídají $\mathbb{Z}_N$, kde násobíme pevným prvkem $a \in \mathbb{Z}^*_N$. Konkrétně definujeme tento operátor násobení $M_a$ tak, že $M_a \ket{x} = \ket{ax \; (\mathrm{mod} \; N)}$ pro každé $x \in \mathbb{Z}_N$. Všimni si, že je implicitní, že bereme součin modulo $N$ uvnitř ket na pravé straně rovnice. Matematická analýza ukazuje, že $M_a$ je unitární operátor. Navíc se ukazuje, že $M_a$ má páry vlastních vektorů a vlastních hodnot, které nám umožňují propojit řád $r$ prvku $a$ s problémem odhadování fáze. Konkrétně, pro libovolnou volbu $j \in \{0, \dots, r-1\}$ platí, že $\ket{\psi_j} = \frac{1}{\sqrt{r}} \sum^{r-1}_{k=0} \omega^{-jk}_{r} \ket{a^k}$ je vlastním vektorem $M_a$, jehož odpovídající vlastní hodnota je $\omega^{j}_{r}$, kde $\omega^{j}_{r} = e^{2 \pi i \frac{j}{r}}.$ Pozorováním vidíme, že vhodným párem vlastní vektor/vlastní hodnota je stav $\ket{\psi_1}$ s $\omega^{1}_{r} = e^{2 \pi i \frac{1}{r}}$. Kdybychom tedy dokázali najít vlastní vektor $\ket{\psi_1}$, mohli bychom odhadnout fázi $\theta=1/r$ pomocí našeho kvantového Circuit a tím získat odhad řádu $r$. To však není snadné, a musíme proto zvážit alternativu.

Uvažujme, k jakému výsledku by Circuit dospěl, kdybychom jako počáteční stav připravili výpočetní stav $\ket{1}$. To není vlastní stav $M_a$, ale je to rovnoměrná superpozice vlastních stavů, které jsme právě popsali. Jinými slovy, platí následující vztah. $\ket{1} = \frac{1}{\sqrt{r}} \sum^{r-1}_{k=0} \ket{\psi_k}$ Z výše uvedené rovnice plyne, že pokud nastavíme počáteční stav na $\ket{1}$, dostaneme přesně stejný výsledek měření, jako kdybychom zvolili $k \in \{ 0, \dots, r-1\}$ rovnoměrně náhodně a použili $\ket{\psi_k}$ jako vlastní vektor v Circuit pro odhadování fáze. Jinými slovy, měření horních $m$ Qubitů dá aproximaci $y / 2^m$ hodnoty $k / r$, kde $k \in \{ 0, \dots, r-1\}$ je zvolen rovnoměrně náhodně. To nám po několika nezávislých spuštěních umožní zjistit $r$ s vysokou mírou jistoty, což byl náš cíl.

Operátory modulárního umocňování

Doposud jsme propojili problém odhadování fáze s problémem hledání řádu definováním $U = M_a$ a $\ket{\psi} = \ket{1}$ v našem kvantovém Circuit. Posledním zbývajícím prvkem je tedy najít efektivní způsob, jak definovat modulární mocniny $M_a$ jako $M_a^k$ pro $k = 1, 2, 4, \dots, 2^{m-1}$. Pro provedení tohoto výpočtu zjistíme, že pro libovolnou mocninu $k$ můžeme sestavit Circuit pro $M_a^k$ nikoli opakováním Circuit pro $M_a$ $k$-krát, ale místo toho výpočtem $b = a^k \; \mathrm{mod} \; N$ a následným použitím Circuit pro $M_b$. Protože potřebujeme pouze mocniny, které jsou samy mocninami 2, lze to klasicky efektivně provést pomocí iterativního umocňování.

Krok 2: Optimalizace problému pro spuštění na kvantovém hardwaru

Konkrétní příklad s $N = 15$ a $a=2$

Zde si můžeme udělat pauzu a probrat konkrétní příklad a sestavit Circuit pro hledání řádu pro $N=15$. Všimni si, že možné netriviální $a \in \mathbb{Z}_N^*$ pro $N=15$ jsou $a \in \{2, 4, 7, 8, 11, 13, 14 \}$. Pro tento příklad volíme $a=2$. Sestavíme operátor $M_2$ a operátory modulárního umocňování $M_2^k$. Akce $M_2$ na výpočetních bázových stavech je následující. $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}$ Pozorováním vidíme, že bázové stavy jsou přeuspořádány, takže máme permutační matici. Tuto operaci můžeme na čtyřech Qubitech sestavit pomocí Gate pro prohazování (swap). Níže sestavíme operace $M_2$ a kontrolovaný-$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 působící na více než dva Qubity budou dále dekomponovány na dvouQubitové Gate.

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

Nyní musíme sestavit operátory modulárního umocňování. Pro dostatečnou přesnost při odhadování fáze použijeme osm Qubitů pro měření odhadu. Proto musíme sestavit $M_b$ s $b = a^{2^k} \; (\mathrm{mod} \; N)$ pro každé $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 vidíme ze seznamu hodnot $b$, kromě $M_2$, který jsme dříve sestavili, musíme také sestavit $M_4$ a $M_1$. Všimni si, že $M_1$ působí triviálně na výpočetní bázové stavy, takže je to jednoduše operátor identity.

$M_4$ působí na výpočetní bázové stavy takto. $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}$

Tuto permutaci lze tedy sestavit následující operací prohazování.

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 působící na více než dva Qubity budou dále dekomponovány na dvouqubitové Gate.

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

Viděli jsme, že operátory $M_b$ pro dané $b \in \mathbb{Z}^*_N$ jsou permutační operace. Díky relativně malé velikosti permutačního problému, který máme zde, protože $N=15$ vyžaduje pouze čtyři Qubity, jsme dokázali tyto operace přímo syntetizovat pomocí Gate SWAP vizuální inspekcí. Obecně to nemusí být škálovatelný přístup. Místo toho bychom možná museli permutační matici sestavit explicitně a použít třídu UnitaryGate z Qiskit a transpilační metody pro syntézu této permutační matice. To však může vést k výrazně hlubším Circuit. Příklad následuje.

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

Porovnejme tyto počty s hloubkou zkompilovaného Circuit z naší ruční implementace 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 vidíme, přístup s permutační maticí vedl k výrazně hlubšímu Circuit i pro jediný Gate $M_2$ ve srovnání s naší ruční implementací. Proto budeme pokračovat s naší předchozí implementací operací $M_b$. Nyní jsme připraveni sestavit úplný Circuit pro hledání řádu pomocí dříve definovaných kontrolovaných operátorů modulárního umocňování. V následujícím kódu také importujeme Circuit QFT z knihovny Qiskit Circuit library, který využívá Hadamardovy Gate na každém Qubitu, sérii kontrolovaných-U1 (nebo Z, v závislosti na fázi) Gate a vrstvu 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)

Všimni si, že jsme vynechali kontrolované operace modulárního umocňování ze zbývajících řídicích Qubitů, protože $M_1$ je operátor identity. Všimni si, že později v tomto tutoriálu spustíme tento Circuit na Backend ibm_marrakesh. K tomu transpilujeme Circuit pro tento konkrétní Backend a reportujeme hloubku Circuit a počty 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: Spuštění pomocí Qiskit primitiv

Nejprve si probereme, co bychom teoreticky získali, kdybychom tento Circuit spustili na ideálním simulátoru. Níže máme sadu výsledků simulace výše uvedeného Circuitu s 1024 snímky. Jak vidíme, dostaneme přibližně rovnoměrné rozdělení přes čtyři bitové řetězce na řídicích Qubitech.

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

plot_histogram(counts)

Měřením řídicích Qubitů získáme osmibitový odhad fáze operátoru $M_a$. Tuto binární reprezentaci můžeme převést na desetinné číslo a najít naměřenou fázi. Jak vidíme z výše uvedeného histogramu, byly naměřeny čtyři různé bitové řetězce a každý z nich odpovídá hodnotě fáze následovně.

# 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

Připomeňme, že každá naměřená fáze odpovídá $\theta = k / r$, kde $k$ je vzorkováno rovnoměrně náhodně z $\{0, 1, \dots, r-1 \}$. Proto můžeme použít algoritmus řetězových zlomků k pokusu o nalezení $k$ a řádu $r$. Python má tuto funkci vestavěnou. Pomocí modulu fractions můžeme převést desetinné číslo na objekt Fraction, například:

Fraction(0.666)

Fraction(5998794703657501, 9007199254740992)

Protože toto vrací zlomky, které vrací výsledek přesně (v tomto případě 0.6660000...), může to dávat nepříjemné výsledky jako ten výše. Metodu .limit_denominator() můžeme použít k získání zlomku, který nejlépe odpovídá našemu desetinnému číslu, se jmenovatelem pod určitou hodnotou:

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

Fraction(2, 3)

To je mnohem lepší. Řád (r) musí být menší než N, takže nastavíme maximální jmenovatel 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

Vidíme, že dvě z naměřených vlastních hodnot nám poskytly správný výsledek: $r=4$, a vidíme, že Shorův algoritmus pro hledání řádu má šanci selhat. Tyto špatné výsledky jsou způsobeny tím, že $k = 0$, nebo tím, že $k$ a $r$ nejsou nesoudělná čísla — a místo $r$ dostaneme faktor $r$. Nejjednodušším řešením je prostě opakovat experiment, dokud nedostaneme uspokojivý výsledek pro $r$. Doposud jsme implementovali problém hledání řádu pro $N=15$ s $a=2$ pomocí Circuitu odhadu fáze na simulátoru. Posledním krokem Shorova algoritmu bude propojit problém hledání řádu s problémem faktorizace celých čísel. Tato poslední část algoritmu je čistě klasická a lze ji vyřešit na klasickém počítači po získání fázových měření z kvantového počítače. Proto odkládáme poslední část algoritmu na dobu po tom, co ukážeme, jak můžeme spustit Circuit hledání řádu na skutečném hardwaru.

Spuštění na hardwaru

Nyní můžeme spustit Circuit hledání řádu, který jsme dříve transpilovali pro ibm_marrakesh. Zde se obracíme k dynamickému oddělování (DD) pro potlačení chyb a k twirlingování Gate pro účely zmírnění chyb. DD zahrnuje aplikaci sekvencí přesně načasovaných řídicích pulzů na kvantové zařízení, čímž efektivně průměruje nežádoucí interakce s prostředím a dekoherenci. Twirlingování Gate naproti tomu náhodně transformuje konkrétní kvantové Gate tak, aby převedlo koherentní chyby na Pauliho chyby, které se hromadí lineárně, nikoli kvadraticky. Obě techniky jsou často kombinovány ke zvýšení koherence a věrnosti kvantových výpočtů.

# 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

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 vidíme, získali jsme stejné bitové řetězce s nejvyššími počty. Protože kvantový hardware obsahuje šum, dochází k úniku do jiných bitových řetězců, které můžeme statisticky odfiltrovat.

# 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: Následné zpracování a vrácení výsledku v požadovaném klasickém formátu

Faktorizace celých čísel

Doposud jsme diskutovali, jak můžeme implementovat problém hledání řádu pomocí Circuitu odhadu fáze. Nyní propojíme problém hledání řádu s faktorizací celých čísel, čímž dokončíme Shorův algoritmus. Všimni si, že tato část algoritmu je klasická. Nyní to ukážeme na našem příkladu $N = 15$ a $a = 2$. Připomeňme, že naměřená fáze je $k / r$, kde $a^r \; (\textrm{mod} \; N) = 1$ a $k$ je náhodné celé číslo mezi $0$ a $r - 1$. Z této rovnice máme $(a^r - 1) \; (\textrm{mod} \; N) = 0,$ což znamená, že $N$ musí dělit $a^r-1$. Pokud je $r$ také sudé, pak můžeme psát $a^r -1 = (a^{r/2}-1)(a^{r/2}+1).$ Pokud $r$ není sudé, nemůžeme pokračovat a musíme to zkusit znovu s jinou hodnotou $a$; jinak existuje vysoká pravděpodobnost, že největší společný dělitel $N$ a buď $a^{r/2}-1$, nebo $a^{r/2}+1$ je vlastní faktor $N$.

Protože některé běhy algoritmu statisticky selžou, budeme tento algoritmus opakovat, dokud nenajdeme alespoň jeden faktor $N$. Buňka níže opakuje algoritmus, dokud není nalezen alespoň jeden faktor $N=15$. Použijeme výsledky hardwarového běhu výše k odhadu fáze a odpovídajícího faktoru v každé iteraci.

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

Diskuse

Související práce

V této části diskutujeme o dalších milníkových pracích, které demonstrovaly Shorův algoritmus na skutečném hardwaru.

Přelomová práce [3] od IBM® demonstrovala Shorův algoritmus poprvé, přičemž faktorizovala číslo 15 na jeho prvočinitele 3 a 5 pomocí sedmiqubitového kvantového počítače na bázi nukleární magnetické rezonance (NMR). Jiný experiment [4] faktorizoval 15 pomocí fotonových Qubitů. Výzkumníci snížili požadovaný počet Qubitů na třetinu standardního protokolu opakovaným využíváním jediného Qubitu a kódováním pracovního registru do stavů vyšší dimenze, přičemž použili dvoufotonový kompilovaný algoritmus. Významnou prací v demonstraci Shorova algoritmu je [5], která používá Kitaevovu iterativní techniku odhadu fáze [8] ke snížení požadavků algoritmu na Qubity. Autoři použili sedm řídicích Qubitů a čtyři cache Qubity spolu s implementací modulárních multiplikátorů. Tato implementace však vyžaduje měření uprostřed Circuitu s operacemi dopředné vazby a recyklaci Qubitů s operacemi resetování. Tato demonstrace byla provedena na iontovém kvantovém počítači.

Novější práce [6] se zaměřila na faktorizaci čísel 15, 21 a 35 na hardwaru IBM Quantum®. Podobně jako předchozí práce použili výzkumníci kompilovanou verzi algoritmu, která využívá poloklasickou kvantovou Fourierovu transformaci navrženou Kitaevem k minimalizaci počtu fyzických Qubitů a Gate. Nejnovější práce [7] také provedla demonstraci konceptu faktorizace celého čísla 21. Tato demonstrace rovněž zahrnovala použití kompilované verze rutiny kvantového odhadu fáze a navázala na předchozí demonstraci [4]. Autoři tuto práci rozšířili použitím konfigurace aproximovaných Toffoliho Gate se zbytkovými fázovými posuny. Algoritmus byl implementován na kvantových procesorech IBM s použitím pouhých pěti Qubitů a přítomnost provázanosti mezi řídicími a registrovými Qubity byla úspěšně ověřena.

Škálování algoritmu

Poznamenáváme, že šifrování RSA obvykle zahrnuje délky klíčů v řádu 2048 až 4096 bitů. Pokus o faktorizaci 2048bitového čísla pomocí Shorova algoritmu povede ke kvantovému Circuitu s miliony Qubitů, včetně režijních nákladů na opravu chyb, a hloubkou Circuitu v řádu miliardy, což je za hranicemi možností současného kvantového hardwaru. Proto bude Shorův algoritmus vyžadovat buď optimalizované metody konstrukce Circuitu, nebo robustní kvantovou opravu chyb, aby byl prakticky využitelný pro prolomení moderních kryptografických systémů. Odkazujeme tě na [9] pro podrobnější diskusi o odhadech zdrojů pro Shorův algoritmus.

Výzva

Gratulujeme k dokončení tutoriálu! Nyní je skvělý čas otestovat své porozumění. Zkus sestavit Circuit pro faktorizaci čísla 21? Můžeš si zvolit vlastní hodnotu $a$. Budeš muset rozhodnout o bitové přesnosti algoritmu, abys zvolil počet Qubitů, a budeš muset navrhnout operátory modulárního umocňování $M_a$. Doporučujeme ti to zkusit sám a poté si přečíst o metodologiích znázorněných na Obr. 9 v [6] a Obr. 2 v [7].

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

    # Your code here
    pass

Reference

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.