Shors Algorithmus

Geschätzter Rechenaufwand: Drei Sekunden auf einem Eagle-r3-Prozessor (HINWEIS: Dies ist nur eine Schätzung. Deine tatsächliche Laufzeit kann abweichen.)

Shors Algorithmus, der 1994 von Peter Shor entwickelt wurde, ist ein bahnbrechender Quantenalgorithmus zur Faktorisierung ganzer Zahlen in polynomialer Zeit. Seine Bedeutung liegt darin, dass er große ganze Zahlen exponentiell schneller faktorisieren kann als jeder bekannte klassische Algorithmus, was die Sicherheit weit verbreiteter kryptographischer Systeme wie RSA bedroht, die auf der Schwierigkeit der Faktorisierung großer Zahlen beruhen. Durch die effiziente Lösung dieses Problems auf einem hinreichend leistungsstarken Quantencomputer könnte Shors Algorithmus Bereiche wie Kryptographie, Cybersicherheit und rechnergestützte Mathematik revolutionieren und unterstreicht damit die transformative Kraft der Quantenberechnung.

Dieses Tutorial zeigt Shors Algorithmus, indem die Zahl 15 auf einem Quantencomputer faktorisiert wird.

Zunächst definieren wir das Ordnungsfindungsproblem und konstruieren zugehörige Circuits aus dem Protokoll der Quantenphasenschätzung. Anschließend führen wir die Ordnungsfindungs-Circuits auf echter Hardware aus und verwenden dabei die Circuit mit der geringsten Tiefe, die wir transpilieren können. Im letzten Abschnitt vervollständigen wir Shors Algorithmus, indem wir das Ordnungsfindungsproblem mit der ganzzahligen Faktorisierung verknüpfen.

Zum Abschluss des Tutorials diskutieren wir weitere Demonstrationen von Shors Algorithmus auf echter Hardware, mit Schwerpunkt sowohl auf generischen Implementierungen als auch auf solchen, die speziell auf die Faktorisierung bestimmter ganzer Zahlen wie 15 und 21 zugeschnitten sind. Hinweis: Dieses Tutorial konzentriert sich mehr auf die Implementierung und Demonstration der Circuits zu Shors Algorithmus. Für eine tiefergehende Erklärung des Lernstoffs empfehlen wir den Kurs Fundamentals of quantum algorithms von Dr. John Watrous sowie die Artikel im Abschnitt Referenzen.

Voraussetzungen

Stelle vor Beginn dieses Tutorials sicher, dass Folgendes installiert ist:

• Qiskit SDK v2.0 oder höher, mit Unterstützung für Visualisierung
• Qiskit Runtime v0.40 oder höher (pip install qiskit-ibm-runtime)

Einrichtung

# 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

Schritt 1: Klassische Eingaben auf ein Quantenproblem abbilden

Hintergrund

Shors Algorithmus zur ganzzahligen Faktorisierung nutzt ein vermittelndes Problem, das als Ordnungsfindungsproblem bekannt ist. In diesem Abschnitt zeigen wir, wie das Ordnungsfindungsproblem mithilfe der Quantenphasenschätzung gelöst werden kann.

Phasenschätzungsproblem

Beim Phasenschätzungsproblem ist uns ein Quantenzustand $\ket{\psi}$ aus $n$ Qubits gegeben, zusammen mit einem unitären Quantum Circuit, der auf $n$ Qubits wirkt. Es wird vorausgesetzt, dass $\ket{\psi}$ ein Eigenvektor der unitären Matrix $U$ ist, die die Wirkung des Circuits beschreibt, und unser Ziel ist es, den Eigenwert $\lambda = e^{2 \pi i \theta}$ zu berechnen oder anzunähern, dem $\ket{\psi}$ entspricht. Mit anderen Worten: Der Circuit soll eine Näherung an die Zahl $\theta \in [0, 1)$ ausgeben, die $U \ket{\psi}= e^{2 \pi i \theta} \ket{\psi}$ erfüllt. Das Ziel des Phasenschätzungs-Circuits ist es, $\theta$ auf $m$ Bits anzunähern. Mathematisch gesprochen möchten wir $y$ so finden, dass $\theta \approx y / 2^m$ gilt, wobei $y \in {0, 1, 2, \dots, 2^{m-1}}$. Das folgende Bild zeigt den Quantum Circuit, der $y$ auf $m$ Bits schätzt, indem eine Messung an $m$ Qubits durchgeführt wird. Im obigen Circuit werden die oberen $m$ Qubits im Zustand $\ket{0^m}$ initialisiert und die unteren $n$ Qubits in $\ket{\psi}$, das bekanntlich ein Eigenvektor von $U$ ist. Die erste Zutat des Phasenschätzungs-Circuits sind die kontrollierten unitären Operationen, die dafür verantwortlich sind, einen Phase Kickback auf ihr jeweiliges Kontroll-Qubit durchzuführen. Diese kontrollierten Unitären werden entsprechend der Position des Kontroll-Qubits potenziert, angefangen beim niedrigstwertigen Bit bis hin zum höchstwertigen Bit. Da $\ket{\psi}$ ein Eigenvektor von $U$ ist, wird der Zustand der unteren $n$ Qubits durch diese Operation nicht beeinflusst, aber die Phaseninformation des Eigenwertes pflanzt sich auf die oberen $m$ Qubits fort. Es stellt sich heraus, dass nach der Phase-Kickback-Operation über kontrollierte Unitäre alle möglichen Zustände der oberen $m$ Qubits für jeden Eigenvektor $\ket{\psi}$ des Unitären $U$ orthonormal zueinander sind. Daher sind diese Zustände perfekt unterscheidbar, und wir können die Basis, die sie bilden, zurück in die Rechenbasis drehen, um eine Messung vorzunehmen. Eine mathematische Analyse zeigt, dass diese Rotationsmatrix der inversen Quanten-Fourier-Transformation (QFT) im $2^m$-dimensionalen Hilbertraum entspricht. Die Intuition dahinter ist, dass die periodische Struktur der modularen Exponentiationsoperatoren im Quantenzustand kodiert ist und die QFT diese Periodizität in messbare Peaks im Frequenzbereich umwandelt.

Für ein tiefergehendes Verständnis, warum der QFT-Circuit in Shors Algorithmus eingesetzt wird, verweisen wir auf den Kurs Fundamentals of quantum algorithms. Wir sind nun bereit, den Phasenschätzungs-Circuit für die Ordnungsfindung zu verwenden.

Ordnungsfindungsproblem

Um das Ordnungsfindungsproblem zu definieren, beginnen wir mit einigen zahlentheoretischen Konzepten. Zunächst definieren wir für jede gegebene positive ganze Zahl $N$ die Menge $\mathbb{Z}_N$ als $\mathbb{Z}_N = \{0, 1, 2, \dots, N-1\}.$ Alle arithmetischen Operationen in $\mathbb{Z}_N$ werden modulo $N$ durchgeführt. Insbesondere sind alle Elemente $a \in \mathbb{Z}_n$, die teilerfremd zu $N$ sind, besonders und bilden $\mathbb{Z}^*_N$ als $\mathbb{Z}^*_N = \{ a \in \mathbb{Z}_N : \mathrm{gcd}(a, N)=1 \}.$ Für ein Element $a \in \mathbb{Z}^*_N$ wird die kleinste positive ganze Zahl $r$ mit $a^r \equiv 1 \; (\mathrm{mod} \; N)$ als die Ordnung von $a$ modulo $N$ definiert. Wie wir später sehen werden, ermöglicht uns das Finden der Ordnung eines $a \in \mathbb{Z}^*_N$ die Faktorisierung von $N$. Um den Ordnungsfindungs-Circuit aus dem Phasenschätzungs-Circuit zu konstruieren, benötigen wir zwei Überlegungen. Erstens müssen wir die Unitäre $U$ definieren, die uns das Finden der Ordnung $r$ ermöglicht, und zweitens müssen wir einen Eigenvektor $\ket{\psi}$ von $U$ definieren, um den Anfangszustand des Phasenschätzungs-Circuits vorzubereiten.

Um das Ordnungsfindungsproblem mit der Phasenschätzung zu verknüpfen, betrachten wir die Operation, die auf einem System definiert ist, dessen klassische Zustände $\mathbb{Z}_N$ entsprechen, wobei wir mit einem festen Element $a \in \mathbb{Z}^*_N$ multiplizieren. Wir definieren diesen Multiplikationsoperator $M_a$ so, dass $M_a \ket{x} = \ket{ax \; (\mathrm{mod} \; N)}$ für jedes $x \in \mathbb{Z}_N$ gilt. Es wird implizit vorausgesetzt, dass wir das Produkt modulo $N$ innerhalb des Kets auf der rechten Seite der Gleichung nehmen. Eine mathematische Analyse zeigt, dass $M_a$ ein unitärer Operator ist. Außerdem stellt sich heraus, dass $M_a$ Eigenvektor-Eigenwert-Paare hat, die es uns ermöglichen, die Ordnung $r$ von $a$ mit dem Phasenschätzungsproblem zu verknüpfen. Für jede Wahl von $j \in \{0, \dots, r-1\}$ gilt konkret, dass $\ket{\psi_j} = \frac{1}{\sqrt{r}} \sum^{r-1}_{k=0} \omega^{-jk}_{r} \ket{a^k}$ ein Eigenvektor von $M_a$ ist, dessen zugehöriger Eigenwert $\omega^{j}_{r}$ ist, wobei $\omega^{j}_{r} = e^{2 \pi i \frac{j}{r}}.$ Durch Beobachtung sehen wir, dass ein praktisches Eigenvektor-Eigenwert-Paar der Zustand $\ket{\psi_1}$ mit $\omega^{1}_{r} = e^{2 \pi i \frac{1}{r}}$ ist. Wenn wir also den Eigenvektor $\ket{\psi_1}$ finden könnten, könnten wir die Phase $\theta=1/r$ mit unserem Quantum Circuit schätzen und damit eine Schätzung der Ordnung $r$ erhalten. Das ist jedoch nicht einfach, und wir müssen eine Alternative in Betracht ziehen.

Überlegen wir, was der Circuit ergeben würde, wenn wir den Rechenzustand $\ket{1}$ als Anfangszustand vorbereiten. Dies ist kein Eigenzustand von $M_a$, aber es ist die gleichmäßige Superposition der Eigenzustände, die wir gerade beschrieben haben. Mit anderen Worten gilt folgende Beziehung: $\ket{1} = \frac{1}{\sqrt{r}} \sum^{r-1}_{k=0} \ket{\psi_k}$ Die Bedeutung der obigen Gleichung ist, dass wir, wenn wir den Anfangszustand auf $\ket{1}$ setzen, genau das gleiche Messergebnis erhalten, als hätten wir $k \in \{ 0, \dots, r-1\}$ gleichmäßig zufällig gewählt und $\ket{\psi_k}$ als Eigenvektor im Phasenschätzungs-Circuit verwendet. Mit anderen Worten liefert eine Messung der oberen $m$ Qubits eine Näherung $y / 2^m$ an den Wert $k / r$, wobei $k \in \{ 0, \dots, r-1\}$ gleichmäßig zufällig gewählt wird. Dies ermöglicht es uns, $r$ nach mehreren unabhängigen Läufen mit hoher Konfidenz zu bestimmen, was unser Ziel war.

Modulare Exponentiationsoperatoren

Bisher haben wir das Phasenschätzungsproblem mit dem Ordnungsfindungsproblem verknüpft, indem wir $U = M_a$ und $\ket{\psi} = \ket{1}$ in unserem Quantum Circuit definiert haben. Daher besteht die letzte verbleibende Zutat darin, einen effizienten Weg zu finden, modulare Potenzen von $M_a$ als $M_a^k$ für $k = 1, 2, 4, \dots, 2^{m-1}$ zu definieren. Zur Durchführung dieser Berechnung stellen wir fest, dass wir für jede gewählte Potenz $k$ einen Circuit für $M_a^k$ erstellen können – nicht indem wir den Circuit für $M_a$ $k$-mal iterieren, sondern indem wir $b = a^k \; \mathrm{mod} \; N$ berechnen und dann den Circuit für $M_b$ verwenden. Da wir nur Potenzen benötigen, die selbst Zweierpotenzen sind, können wir dies klassisch effizient durch iteriertes Quadrieren bewerkstelligen.

Schritt 2: Problem für die Ausführung auf Quantenhardware optimieren

Konkretes Beispiel mit $N = 15$ und $a=2$

Hier können wir innehalten, um ein konkretes Beispiel zu besprechen und den Ordnungsfindungs-Circuit für $N=15$ zu konstruieren. Beachte, dass mögliche nichttriviale $a \in \mathbb{Z}_N^*$ für $N=15$ die Werte $a \in \{2, 4, 7, 8, 11, 13, 14 \}$ sind. Für dieses Beispiel wählen wir $a=2$. Wir konstruieren den Operator $M_2$ und die modularen Exponentiationsoperatoren $M_2^k$. Die Wirkung von $M_2$ auf die Rechenbasis-Zustände ist wie folgt. $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}$ Durch Inspektion erkennen wir, dass die Basiszustände umsortiert werden, sodass wir eine Permutationsmatrix haben. Wir können diese Operation auf vier Qubits mit SWAP-Gates konstruieren. Im Folgenden konstruieren wir die Operationen $M_2$ und kontrolliertes-$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)

Gates, die auf mehr als zwei Qubits wirken, werden weiter in Zwei-Qubit-Gates zerlegt.

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

Nun müssen wir die modularen Exponentiationsoperatoren konstruieren. Um eine ausreichende Genauigkeit bei der Phasenschätzung zu erreichen, verwenden wir acht Qubits für die Schätzmessung. Daher müssen wir $M_b$ mit $b = a^{2^k} \; (\mathrm{mod} \; N)$ für jedes $k = 0, 1, \dots, 7$ konstruieren.

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]

Wie wir aus der Liste der $b$-Werte erkennen können, müssen wir neben $M_2$, das wir zuvor konstruiert haben, auch $M_4$ und $M_1$ bauen. Beachte, dass $M_1$ trivial auf die Rechenbasis-Zustände wirkt und daher einfach der Identitätsoperator ist.

$M_4$ wirkt auf die Rechenbasis-Zustände wie folgt. $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}$

Daher kann diese Permutation mit der folgenden SWAP-Operation konstruiert werden.

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)

Gates, die auf mehr als zwei Qubits wirken, werden weiter in Zwei-Qubit-Gates zerlegt.

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

Wir haben gesehen, dass $M_b$-Operatoren für ein gegebenes $b \in \mathbb{Z}^*_N$ Permutationsoperationen sind. Aufgrund der vergleichsweise kleinen Größe des Permutationsproblems, das wir hier haben – da $N=15$ nur vier Qubits erfordert –, konnten wir diese Operationen direkt durch Inspektion mit SWAP-Gates synthetisieren. Im Allgemeinen ist dies möglicherweise kein skalierbarer Ansatz. Stattdessen müssten wir die Permutationsmatrix explizit konstruieren und Qiskits UnitaryGate-Klasse sowie Transpilierungsmethoden verwenden, um diese Permutationsmatrix zu synthetisieren. Dies kann jedoch zu deutlich tieferen Circuits führen. Ein Beispiel folgt.

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

Vergleichen wir diese Zählwerte mit der kompilierten Circuit-Tiefe unserer manuellen Implementierung des $M_2$-Gates.

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

Wie wir sehen können, führte der Permutationsmatrix-Ansatz selbst für ein einzelnes $M_2$-Gate zu einem deutlich tieferen Circuit als unsere manuelle Implementierung. Daher fahren wir mit unserer früheren Implementierung der $M_b$-Operationen fort. Nun sind wir bereit, den vollständigen Ordnungsfindungs-Circuit mit unseren zuvor definierten kontrollierten modularen Exponentiationsoperatoren zu konstruieren. Im folgenden Code importieren wir auch den QFT-Circuit aus der Qiskit Circuit Library, der Hadamard-Gates auf jedem Qubit, eine Reihe von kontrollierten U1- (oder Z-, je nach Phase) Gates und eine Schicht von SWAP-Gates verwendet.

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

Beachte, dass wir die kontrollierten modularen Exponentiationsoperationen von den verbleibenden Kontroll-Qubits weggelassen haben, weil $M_1$ der Identitätsoperator ist. Beachte, dass wir diesen Circuit später in diesem Tutorial auf dem Backend ibm_marrakesh ausführen werden. Dazu transpilieren wir den Circuit für dieses spezifische Backend und berichten über die Circuit-Tiefe und Gate-Anzahlen.

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

Schritt 3: Ausführung mit Qiskit-Primitiven

Zunächst besprechen wir, was wir theoretisch erhalten würden, wenn wir diesen Circuit auf einem idealen Simulator ausführen. Unten haben wir eine Reihe von Simulationsergebnissen des obigen Circuits mit 1024 Shots. Wie wir sehen können, erhalten wir eine annähernd gleichmäßige Verteilung über vier Bitstrings der Kontroll-Qubits.

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

plot_histogram(counts)

Durch das Messen der Kontroll-Qubits erhalten wir eine Acht-Bit-Phasenschätzung des $M_a$-Operators. Wir können diese Binärdarstellung in Dezimalzahlen umwandeln, um die gemessene Phase zu finden. Wie wir aus dem obigen Histogramm sehen können, wurden vier verschiedene Bitstrings gemessen, und jeder von ihnen entspricht einem Phasenwert wie folgt.

# 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

Erinnere dich: Jede gemessene Phase entspricht $\theta = k / r$, wobei $k$ gleichmäßig zufällig aus $\{0, 1, \dots, r-1 \}$ gezogen wird. Daher können wir den Kettenbruchalgorithmus verwenden, um $k$ und die Ordnung $r$ zu finden. Python hat diese Funktionalität bereits eingebaut. Wir können das fractions-Modul verwenden, um einen Float in ein Fraction-Objekt umzuwandeln, zum Beispiel:

Fraction(0.666)

Fraction(5998794703657501, 9007199254740992)

Da dies Brüche liefert, die das Ergebnis exakt zurückgeben (in diesem Fall 0.6660000...), kann dies zu unhandlichen Ergebnissen wie dem obigen führen. Wir können die Methode .limit_denominator() verwenden, um den Bruch zu erhalten, der unserem Float am nächsten kommt, mit einem Nenner unterhalb eines bestimmten Werts:

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

Fraction(2, 3)

Das ist viel übersichtlicher. Die Ordnung (r) muss kleiner als N sein, also setzen wir den maximalen Nenner auf 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

Wir können sehen, dass zwei der gemessenen Eigenwerte uns das richtige Ergebnis geliefert haben: $r=4$, und wir sehen, dass Shors Ordnungsfindungsalgorithmus eine Möglichkeit des Scheiterns hat. Diese schlechten Ergebnisse entstehen, weil $k = 0$ ist oder weil $k$ und $r$ nicht teilerfremd sind – und anstelle von $r$ erhalten wir dann einen Teiler von $r$. Die einfachste Lösung dafür ist, das Experiment einfach zu wiederholen, bis wir ein befriedigendes Ergebnis für $r$ erhalten. Bisher haben wir das Ordnungsfindungsproblem für $N=15$ mit $a=2$ mithilfe des Phasenschätzungs-Circuits auf einem Simulator implementiert. Der letzte Schritt von Shors Algorithmus besteht darin, das Ordnungsfindungsproblem mit dem ganzzahligen Faktorisierungsproblem zu verknüpfen. Dieser letzte Teil des Algorithmus ist rein klassisch und kann auf einem klassischen Computer gelöst werden, nachdem die Phasenmessungen von einem Quantencomputer erhalten wurden. Daher verschieben wir den letzten Teil des Algorithmus, bis wir demonstriert haben, wie wir den Ordnungsfindungs-Circuit auf echter Hardware ausführen können.

Hardware-Läufe

Nun können wir den Ordnungsfindungs-Circuit ausführen, den wir zuvor für ibm_marrakesh transpiliert haben. Hier greifen wir auf dynamisches Entkoppeln (DD) zur Fehlerunterdrückung und auf Gate Twirling zur Fehlerminderung zurück. DD beinhaltet das Anwenden von Sequenzen präzise getakteter Steuerimpulse auf ein Quantengerät, wodurch unerwünschte Umweltwechselwirkungen und Dekohärenz effektiv herausgemittelt werden. Gate Twirling hingegen randomisiert bestimmte Quantum Gates, um kohärente Fehler in Pauli-Fehler umzuwandeln, die sich linear statt quadratisch akkumulieren. Beide Techniken werden oft kombiniert, um die Kohärenz und Fidelity von Quantenberechnungen zu verbessern.

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

Wie wir sehen können, haben wir dieselben Bitstrings mit den höchsten Zählwerten erhalten. Da Quantenhardware Rauschen aufweist, gibt es einen gewissen Leckstrom zu anderen Bitstrings, den wir statistisch herausfiltern können.

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

Schritt 4: Nachverarbeitung und Rückgabe des Ergebnisses im gewünschten klassischen Format

Ganzzahlige Faktorisierung

Bisher haben wir besprochen, wie wir das Ordnungsfindungsproblem mithilfe eines Phasenschätzungs-Circuits implementieren können. Nun verknüpfen wir das Ordnungsfindungsproblem mit der ganzzahligen Faktorisierung, womit wir Shors Algorithmus vervollständigen. Beachte, dass dieser Teil des Algorithmus klassisch ist. Wir demonstrieren dies nun anhand unseres Beispiels mit $N = 15$ und $a = 2$. Erinnere dich daran, dass die gemessene Phase $k / r$ ist, wobei $a^r \; (\textrm{mod} \; N) = 1$ und $k$ eine zufällige ganze Zahl zwischen $0$ und $r - 1$ ist. Aus dieser Gleichung ergibt sich $(a^r - 1) \; (\textrm{mod} \; N) = 0,$ was bedeutet, dass $N$ durch $a^r-1$ teilbar sein muss. Wenn $r$ auch gerade ist, können wir $a^r -1 = (a^{r/2}-1)(a^{r/2}+1)$ schreiben. Wenn $r$ nicht gerade ist, können wir nicht weiter vorgehen und müssen es mit einem anderen Wert für $a$ erneut versuchen; andernfalls besteht eine hohe Wahrscheinlichkeit, dass der größte gemeinsame Teiler von $N$ und entweder $a^{r/2}-1$ oder $a^{r/2}+1$ ein echter Faktor von $N$ ist.

Da einige Läufe des Algorithmus statistisch scheitern werden, wiederholen wir diesen Algorithmus, bis mindestens ein Faktor von $N$ gefunden wurde. Die folgende Zelle wiederholt den Algorithmus, bis mindestens ein Faktor von $N=15$ gefunden wurde. Wir verwenden die Ergebnisse des obigen Hardware-Laufs, um in jeder Iteration die Phase und den entsprechenden Faktor zu schätzen.

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

Diskussion

Verwandte Arbeiten

In diesem Abschnitt diskutieren wir andere bedeutende Arbeiten, die Shors Algorithmus auf echter Hardware demonstriert haben.

Die wegweisende Arbeit [3] von IBM® demonstrierte Shors Algorithmus zum ersten Mal und faktorisierte die Zahl 15 in ihre Primfaktoren 3 und 5 mit einem Sieben-Qubit-Kernspinresonanz-Quantencomputer (NMR). Ein weiteres Experiment [4] faktorisierte 15 mithilfe photonischer Qubits. Durch die Verwendung eines einzelnen, mehrfach recycelten Qubits und die Kodierung des Arbeitsregisters in höherdimensionale Zustände reduzierten die Forscher die benötigte Anzahl von Qubits auf ein Drittel des Standardprotokolls und nutzten einen Zwei-Photonen-kompilierten Algorithmus. Eine bedeutende Veröffentlichung zur Demonstration von Shors Algorithmus ist [5], die Kitaevs iterative Phasenschätzung [8] verwendet, um den Qubit-Bedarf des Algorithmus zu reduzieren. Die Autoren verwendeten sieben Kontroll-Qubits und vier Cache-Qubits sowie die Implementierung modularer Multiplizierer. Diese Implementierung erfordert jedoch Mid-Circuit-Messungen mit Feed-Forward-Operationen und Qubit-Recycling mit Reset-Operationen. Diese Demonstration wurde auf einem Ionenfallen-Quantencomputer durchgeführt.

Neuere Arbeiten [6] konzentrierten sich auf die Faktorisierung von 15, 21 und 35 auf IBM Quantum®-Hardware. Ähnlich wie frühere Arbeiten verwendeten die Forscher eine kompilierte Version des Algorithmus, die eine semi-klassische Quanten-Fourier-Transformation nach Kitaev einsetzt, um die Anzahl der physikalischen Qubits und Gates zu minimieren. Eine jüngste Arbeit [7] führte ebenfalls eine Proof-of-Concept-Demonstration zur Faktorisierung der ganzen Zahl 21 durch. Diese Demonstration verwendete ebenfalls eine kompilierte Version der Quantenphasenschätzungsroutine und baute auf der früheren Demonstration von [4] auf. Die Autoren gingen über diese Arbeit hinaus, indem sie eine Konfiguration aus approximativen Toffoli-Gates mit Restphasenverschiebungen verwendeten. Der Algorithmus wurde auf IBM-Quantenprozessoren mit nur fünf Qubits implementiert, und das Vorhandensein von Verschränkung zwischen den Kontroll- und Register-Qubits wurde erfolgreich verifiziert.

Skalierung des Algorithmus

Wir merken an, dass RSA-Verschlüsselung typischerweise Schlüsselgrößen in der Größenordnung von 2048 bis 4096 Bit umfasst. Der Versuch, eine 2048-Bit-Zahl mit Shors Algorithmus zu faktorisieren, würde zu einem Quantum Circuit mit Millionen von Qubits führen, einschließlich des Fehlerkorrektur-Overheads, und einer Circuit-Tiefe in der Größenordnung einer Milliarde, was die Grenzen aktueller Quantenhardware übersteigt. Daher wird Shors Algorithmus entweder optimierte Circuit-Konstruktionsmethoden oder robuste Quantenfehlerkorrektur erfordern, um für das Brechen moderner kryptographischer Systeme praktisch anwendbar zu sein. Wir verweisen auf [9] für eine ausführlichere Diskussion der Ressourcenschätzung für Shors Algorithmus.

Aufgabe

Herzlichen Glückwunsch zum Abschluss des Tutorials! Jetzt ist ein guter Zeitpunkt, dein Verständnis zu testen. Kannst du den Circuit zur Faktorisierung von 21 konstruieren? Du kannst ein $a$ deiner eigenen Wahl auswählen. Du musst die Bit-Genauigkeit des Algorithmus festlegen, um die Anzahl der Qubits zu wählen, und du musst die modularen Exponentiationsoperatoren $M_a$ entwerfen. Wir ermutigen dich, dies selbst auszuprobieren und dann die Methoden in Abb. 9 von [6] und Abb. 2 von [7] nachzulesen.

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

    # Your code here
    pass

Referenzen

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.