Algoritmo di Shor

Stima di utilizzo: Tre secondi su un processore Eagle r3 (NOTA: questa è solo una stima. Il tempo di esecuzione potrebbe variare)

Risultati di apprendimento

Dopo aver seguito questo tutorial, gli utenti dovrebbero aver compreso:

• I fondamenti matematici dell'algoritmo di Shor per la fattorizzazione dei numeri interi
• Come eseguire un'istanza di esempio di questo algoritmo su un sistema hardware

Prerequisiti

Consigliamo agli utenti di acquisire familiarità con i seguenti argomenti prima di seguire questo tutorial:

• Fondamenti degli algoritmi quantistici.
• Stima delle fasi e fattorizzazione. In questo tutorial trattiamo alcuni di questi argomenti.

Sfondo

L'algoritmo di Shor, sviluppato da Peter Shor nel 1994, è un rivoluzionario algoritmo quantistico che consente di scomporre i numeri interi in tempo polinomiale. La sua importanza risiede nella capacità di scomporre in fattori grandi numeri interi in modo esponenzialmente più veloce rispetto a qualsiasi algoritmo classico conosciuto, mettendo a rischio la sicurezza di sistemi crittografici ampiamente utilizzati come l'RSA, che si basano proprio sulla difficoltà di scomporre in fattori numeri di grandi dimensioni. Se questo problema venisse risolto in modo efficiente su un computer quantistico sufficientemente potente, l'algoritmo di Shor potrebbe rivoluzionare settori quali la crittografia, la sicurezza informatica e la matematica computazionale, mettendo in evidenza il potere trasformativo dell'informatica quantistica.

Questa esercitazione si concentra sulla dimostrazione dell'algoritmo di Shor attraverso la fattorizzazione di 15 su un computer quantistico.

In primo luogo, definiamo il problema di ricerca dell'ordine e costruiamo i circuiti corrispondenti a partire dal protocollo di stima della fase quantistica. Successivamente, eseguiamo i circuiti di ricerca dell'ordine su hardware reale, utilizzando i circuiti a profondità più breve che possiamo transpilare. L'ultima sezione completa l'algoritmo di Shor collegando il problema della ricerca degli ordini alla fattorizzazione degli interi.

Concludiamo l'esercitazione con una discussione su altre dimostrazioni dell'algoritmo di Shor su hardware reale, concentrandoci sia sulle implementazioni generiche sia su quelle personalizzate per la fattorizzazione di numeri interi specifici, come 15 e 21.

Nota: questa esercitazione si concentra maggiormente sull'implementazione e sulla dimostrazione dei circuiti relativi all'algoritmo di Shor. Per una risorsa educativa approfondita sul materiale, consultare il corso Fundamentals of quantum algorithms del Dr. John Watrous e i documenti nella sezione Riferimenti.

Requisiti

Prima di iniziare questa esercitazione, assicuratevi di aver installato quanto segue:

• Qiskit SDK v2.0 o versioni successive, con supporto alla visualizzazione
• Qiskit Runtime v0.40 o più tardi (pip install qiskit-ibm-runtime)

Configura

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

Fase 1: mappare gli input classici su un problema quantistico

L'algoritmo di Shor per la fattorizzazione degli interi utilizza un problema intermedio noto come problema di ricerca dell'ordine. In questa sezione dimostriamo come risolvere il problema della ricerca dell'ordine utilizzando la stima quantistica della fase.

Problema di stima di fase

Nel problema della stima di fase, ci viene dato uno stato quantistico $\ket{\psi}$ di $n$ qubit, insieme a un circuito quantistico unitario che agisce su $n$ qubit. Ci è stato promesso che $\ket{\psi}$ è un autovettore della matrice unitaria $U$ che descrive l'azione del circuito e il nostro obiettivo è calcolare o approssimare l'autovalore $\lambda = e^{2 \pi i \theta}$ a cui corrisponde $\ket{\psi}$. In altre parole, il circuito deve fornire un'approssimazione del numero $\theta \in [0, 1)$ soddisfacente $U \ket{\psi}= e^{2 \pi i \theta} \ket{\psi}.$ L'obiettivo del circuito di stima della fase è quello di approssimare $\theta$ in $m$ bit. In termini matematici, vorremmo trovare $y$ tale che $\theta \approx y / 2^m$, dove $y \in {0, 1, 2, \dots, 2^{m-1}}$. L'immagine seguente mostra il circuito quantistico che stima $y$ in $m$ bit effettuando una misura su $m$ qubit.

Nel circuito di cui sopra, i qubit superiori $m$ sono iniziati nello stato $\ket{0^m}$ e i qubit inferiori $n$ sono iniziati nello stato $\ket{\psi}$, che è promesso essere un autovettore di $U$. Il primo ingrediente del circuito di stima della fase sono le operazioni controllate-unitarie che sono responsabili dell'esecuzione di un contraccolpo di fase al loro corrispondente qubit di controllo. Queste unità controllate vengono esponenziate in base alla posizione del qubit di controllo, che va dal bit meno significativo al bit più significativo. Poiché $\ket{\psi}$ è un autovalore di $U$, lo stato dei qubit inferiori di $n$ non è influenzato da questa operazione, ma l'informazione di fase dell'autovalore si propaga ai qubit superiori di $m$.

Si scopre che, dopo l'operazione di contraccolpo di fase tramite le unità controllate, tutti i possibili stati dei qubit top $m$ sono ortonormali tra loro per ogni autovettore $\ket{\psi}$ dell'unità $U$. Pertanto, questi stati sono perfettamente distinguibili e possiamo ruotare la base che formano verso la base computazionale per effettuare una misura. Un'analisi matematica mostra che questa matrice di rotazione corrisponde alla trasformata quantistica di Fourier (QFT) inversa in $2^m$ spazio di Hilbert. L'intuizione è che la struttura periodica degli operatori di esponenziazione modulare è codificata nello stato quantistico e la QFT converte questa periodicità in picchi misurabili nel dominio della frequenza.

Per una comprensione più approfondita del motivo per cui il circuito QFT è impiegato nell'algoritmo di Shor, rimandiamo il lettore al corso Fondamenti di algoritmi quantistici.

Siamo ora pronti a utilizzare il circuito di stima della fase per la ricerca dell'ordine.

Problema con l'ordine

Per definire il problema della ricerca di ordini, iniziamo con alcuni concetti di teoria dei numeri. In primo luogo, per ogni dato intero positivo $N$, definire l'insieme $\mathbb{Z}_N$ come $\mathbb{Z}_N = \{0, 1, 2, \dots, N-1\}.$ Tutte le operazioni aritmetiche in $\mathbb{Z}_N$ sono eseguite modulo $N$. In particolare, tutti gli elementi $a \in \mathbb{Z}_n$ che sono coprimi con $N$ sono speciali e costituiscono $\mathbb{Z}^*_N$ come $\mathbb{Z}^*_N = \{ a \in \mathbb{Z}_N : \mathrm{gcd}(a, N)=1 \}.$ Per un elemento $a \in \mathbb{Z}^*_N$, il più piccolo intero positivo $r$ tale che $a^r \equiv 1 \; (\mathrm{mod} \; N)$ sia definito come l' ordine di $a$ modulo $N$. Come vedremo in seguito, trovare l'ordine di un $a \in \mathbb{Z}^*_N$ ci permetterà di fattorizzare $N$.

Per costruire il circuito di ricerca dell'ordine dal circuito di stima della fase, sono necessarie due considerazioni. In primo luogo, dobbiamo definire l'unità $U$ che ci permetterà di trovare l'ordine $r$, e in secondo luogo, dobbiamo definire un autovettore $\ket{\psi}$ di $U$ per preparare lo stato iniziale del circuito di stima della fase.

Per collegare il problema dell'order finding alla stima di fase, consideriamo l'operazione definita su un sistema i cui stati classici corrispondono a $\mathbb{Z}_N$, dove moltiplichiamo per un elemento fisso $a \in \mathbb{Z}^*_N$. In particolare, definiamo questo operatore di moltiplicazione $M_a$ tale che $M_a \ket{x} = \ket{ax \; (\mathrm{mod} \; N)}$ per ogni $x \in \mathbb{Z}_N$. Si noti che è implicito che stiamo prendendo il prodotto modulo $N$ all'interno del ket sul lato destro dell'equazione. Un'analisi matematica dimostra che $M_a$ è un operatore unitario. Inoltre, si scopre che $M_a$ ha coppie di autovettori e autovalori che ci permettono di collegare l'ordine $r$ di $a$ al problema della stima della fase. In particolare, per qualsiasi scelta di $j \in \{0, \dots, r-1\}$, si ha che $\ket{\psi_j} = \frac{1}{\sqrt{r}} \sum^{r-1}_{k=0} \omega^{-jk}_{r} \ket{a^k}$ è un autovalore di $M_a$ il cui corrispondente autovalore è $\omega^{j}_{r}$, dove $\omega^{j}_{r} = e^{2 \pi i \frac{j}{r}}.$

Dall'osservazione, vediamo che una conveniente coppia autovettore/valore è lo stato $\ket{\psi_1}$ con $\omega^{1}_{r} = e^{2 \pi i \frac{1}{r}}$. Pertanto, se riuscissimo a trovare l'autovalore $\ket{\psi_1}$, potremmo stimare la fase $\theta=1/r$ con il nostro circuito quantistico e quindi ottenere una stima dell'ordine $r$. Tuttavia, non è facile farlo e dobbiamo considerare un'alternativa.

Consideriamo il risultato del circuito se prepariamo lo stato computazionale $\ket{1}$ come stato iniziale. Questo non è un autostato di $M_a$, ma è la sovrapposizione uniforme degli autostati appena descritti. In altre parole, vale la seguente relazione. $\ket{1} = \frac{1}{\sqrt{r}} \sum^{r-1}_{k=0} \ket{\psi_k}$ L'implicazione dell'equazione precedente è che se impostiamo lo stato iniziale su $\ket{1}$, otterremo esattamente lo stesso risultato di misura che avremmo ottenuto se avessimo scelto $k \in \{ 0, \dots, r-1\}$ in modo uniforme e casuale e avessimo usato $\ket{\psi_k}$ come autovettore nel circuito di stima della fase. In altre parole, una misura dei primi $m$ qubit produce un'approssimazione $y / 2^m$ al valore $k / r$ dove $k \in \{ 0, \dots, r-1\}$ è scelto uniformemente a caso. Questo ci permette di imparare $r$ con un alto grado di confidenza dopo diverse esecuzioni indipendenti, che era il nostro obiettivo.

Operatori di esponenziazione modulari

Finora abbiamo collegato il problema della stima della fase al problema della ricerca dell'ordine definendo $U = M_a$ e $\ket{\psi} = \ket{1}$ nel nostro circuito quantistico. Pertanto, l'ultimo ingrediente rimasto è trovare un modo efficiente per definire gli esponenziali modulari di $M_a$ come $M_a^k$ per $k = 1, 2, 4, \dots, 2^{m-1}$. Per eseguire questo calcolo, scopriamo che per qualsiasi potenza $k$ scelta, possiamo creare un circuito per $M_a^k$ non iterando $k$ per il circuito per $M_a$, ma calcolando $b = a^k \; \mathrm{mod} \; N$ e poi usando il circuito per $M_b$. Poiché ci servono solo le potenze che sono potenze di 2, possiamo eseguire questa operazione in modo classicamente efficiente utilizzando la quadratura iterativa.

Fase 2: Ottimizzazione del problema per l'esecuzione su hardware quantistico

Esempio specifico con $N = 15$ e $a=2$

Possiamo fermarci qui per discutere un esempio specifico e costruire il circuito di ricerca dell'ordine per $N=15$. Si noti che i possibili $a \in \mathbb{Z}_N^*$ non banali per $N=15$ sono $a \in \{2, 4, 7, 8, 11, 13, 14 \}$. Per questo esempio, scegliamo $a=2$. Costruiremo l'operatore $M_2$ e gli operatori di esponenziazione modulare $M_2^k$.

L'azione di $M_2$ sugli stati base computazionali è la seguente. $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}$ Osservando, si nota che gli stati base sono mescolati, quindi si ha una matrice di permutazione. Possiamo costruire questa operazione su quattro qubit con porte di scambio. Di seguito, si costruiscono le operazioni $M_2$ e $M_2$ controllate.

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)

Output:

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)

Output:

Le porte che agiscono su più di due qubit saranno ulteriormente scomposte in porte a due qubit.

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

Output:

Ora dobbiamo costruire gli operatori di esponenziazione modulare. Per ottenere una precisione sufficiente nella stima della fase, utilizzeremo otto qubit per la misura della stima. Pertanto, è necessario costruire $M_b$ con $b = a^{2^k} \; (\mathrm{mod} \; N)$ per ogni $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)

Output:

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

Come possiamo vedere dall'elenco dei valori di $b$, oltre a $M_2$ che abbiamo costruito in precedenza, dobbiamo costruire anche $M_4$ e $M_1$. Si noti che $M_1$ agisce banalmente sugli stati della base computazionale, quindi è semplicemente l'operatore di identità.

$M_4$ agisce sugli stati base computazionali come segue. $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}$

Pertanto, questa permutazione può essere costruita con la seguente operazione di scambio.

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)

Output:

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)

Output:

Le porte che agiscono su più di due qubit saranno ulteriormente scomposte in porte a due qubit.

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

Output:

Abbiamo visto che gli operatori di $M_b$ per un dato $b \in \mathbb{Z}^*_N$ sono operazioni di permutazione. A causa delle dimensioni relativamente ridotte del problema di permutazione che abbiamo qui, dato che $N=15$ richiede solo quattro qubit, siamo stati in grado di sintetizzare queste operazioni direttamente con le porte di SWAP mediante un'ispezione. In generale, questo potrebbe non essere un approccio scalabile. Invece, potrebbe essere necessario costruire esplicitamente la matrice di permutazione e utilizzare la classe UnitaryGate di Qiskit e i metodi di transpilazione per sintetizzare questa matrice di permutazione. Tuttavia, questo può portare a circuiti molto più profondi. Segue un esempio.

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
)

Output:

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

Confrontiamo questi conteggi con la profondità del circuito compilato della nostra implementazione manuale della porta $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
)

Output:

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

Come si può notare, l'approccio della matrice di permutazione ha prodotto un circuito significativamente più profondo anche per un singolo gate $M_2$ rispetto alla nostra implementazione manuale. Pertanto, continueremo con la nostra precedente implementazione delle operazioni di $M_b$.

Ora siamo pronti a costruire il circuito di ricerca dell'ordine completo utilizzando gli operatori di esponenziazione modulare controllata definiti in precedenza. Nel codice che segue, importiamo anche il circuito QFT dalla libreria Qiskit Circuit, che utilizza porte Hadamard su ogni qubit, una serie di porte controlled-U1 (o Z, a seconda della fase) e uno strato di porte 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)

Output:

Si noti che abbiamo omesso le operazioni di esponenziazione modulare controllata dai restanti qubit di controllo perché $M_1$ è l'operatore di identità.

Si noti che più avanti in questa esercitazione si eseguirà questo circuito sul backend ibm_marrakesh . A tal fine, transpiliamo il circuito in base a questo backend specifico e riportiamo la profondità del circuito e il numero di porte.

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
)

Output:

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

Passaggio 3: eseguire utilizzando Qiskit primitives

In primo luogo, discutiamo ciò che si otterrebbe teoricamente se si eseguisse questo circuito su un simulatore ideale. Di seguito sono riportati i risultati della simulazione del circuito di cui sopra con 1024 scatti. Come si vede, si ottiene una distribuzione approssimativamente uniforme su quattro stringhe di bit sui qubit di controllo.

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

plot_histogram(counts)

Output:

Misurando i qubit di controllo, otteniamo una stima di fase a otto bit dell'operatore $M_a$. Possiamo convertire questa rappresentazione binaria in decimale per trovare la fase misurata. Come si può vedere dall'istogramma sopra riportato, sono state misurate quattro diverse stringhe di bit, ognuna delle quali corrisponde a un valore di fase come segue.

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

Output:

            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

Ricordiamo che la fase misurata corrisponde a $\theta = k / r$ dove $k$ è campionata in modo uniformemente casuale da $\{0, 1, \dots, r-1 \}$. Pertanto, possiamo usare l'algoritmo delle frazioni continue per cercare di trovare $k$ e l'ordine $r$. Python ha questa funzionalità incorporata. Possiamo usare il modulo fractions per trasformare un galleggiante in un oggetto Fraction , ad esempio:

Fraction(0.666)

Output:

Fraction(5998794703657501, 9007199254740992)

Poiché fornisce frazioni che restituiscono esattamente il risultato (in questo caso, 0.6660000...), può dare risultati strani come quello sopra. Possiamo utilizzare il metodo .limit_denominator() per ottenere la frazione che più si avvicina al nostro galleggiante, con un denominatore inferiore a un certo valore:

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

Output:

Fraction(2, 3)

Questo è molto più bello. L'ordine (r) deve essere inferiore a N, quindi il denominatore massimo sarà 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)

Output:

   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

Possiamo notare che due degli autovalori misurati ci hanno fornito il risultato corretto: $r=4$, e possiamo vedere che l'algoritmo di Shor per la ricerca dell'ordine ha la possibilità di fallire. Questi cattivi risultati sono dovuti al fatto che $k = 0$, o perché $k$ e $r$ non sono coprimari - e invece di $r$, ci viene dato un fattore di $r$. La soluzione più semplice è semplicemente ripetere l'esperimento finché non si ottiene un risultato soddisfacente per $r$.

Finora abbiamo implementato il problema della ricerca dell'ordine per $N=15$ con $a=2$ utilizzando il circuito di stima della fase su un simulatore. L'ultimo passo dell'algoritmo di Shor consiste nel mettere in relazione il problema della ricerca dell'ordine con il problema della fattorizzazione degli interi. Quest'ultima parte dell'algoritmo è puramente classica e può essere risolta su un computer classico dopo che le misure di fase sono state ottenute da un computer quantistico. Pertanto, rimandiamo l'ultima parte dell'algoritmo a quando avremo dimostrato come eseguire il circuito di ricerca degli ordini su un hardware reale.

Esecuzioni hardware

Ora possiamo eseguire il circuito di ricerca dell'ordine che abbiamo precedentemente tradotto per ibm_marrakesh. In questo caso ci rivolgiamo al disaccoppiamento dinamico (DD) per la soppressione degli errori e al gate twirling per la mitigazione degli errori. Il DD comporta l'applicazione di sequenze di impulsi di controllo precisamente temporizzati a un dispositivo quantistico, eliminando in modo efficace le interazioni ambientali indesiderate e la decoerenza. Il gate twirling, invece, randomizza specifiche porte quantistiche per trasformare gli errori coerenti in errori di Pauli, che si accumulano linearmente anziché quadraticamente. Entrambe le tecniche sono spesso combinate per migliorare la coerenza e la fedeltà delle computazioni quantistiche.

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

Output:

Come possiamo vedere, abbiamo ottenuto le stesse bitstring con i conteggi più alti. Poiché l'hardware quantistico ha un rumore, c'è una certa dispersione in altre stringhe di bit, che possiamo filtrare statisticamente.

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

Output:

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

Fase 4: Post-elaborazione e restituzione del risultato nel formato classico desiderato

Fattorizzazione dei numeri interi

Finora abbiamo discusso come implementare il problema della ricerca dell'ordine utilizzando un circuito di stima della fase. Ora colleghiamo il problema della ricerca dell'ordine alla fattorizzazione dei numeri interi, completando così l'algoritmo di Shor. Si noti che questa parte dell'algoritmo è classica.

Lo dimostriamo ora con l'esempio di $N = 15$ e $a = 2$. Ricordiamo che la fase che abbiamo misurato è $k / r$, dove $a^r \; (\textrm{mod} \; N) = 1$ e $k$ è un intero casuale tra $0$ e $r - 1$. Da questa equazione, abbiamo $(a^r - 1) \; (\textrm{mod} \; N) = 0,$ che significa che $N$ deve dividere $a^r-1$. Se anche $r$ è pari, allora possiamo scrivere $a^r -1 = (a^{r/2}-1)(a^{r/2}+1).$ Se $r$ non è pari, non possiamo andare oltre e dobbiamo riprovare con un valore diverso per $a$; altrimenti, c'è un'alta probabilità che il massimo comun divisore di $N$ e $a^{r/2}-1$, o $a^{r/2}+1$ sia un fattore proprio di $N$.

Poiché alcune esecuzioni dell'algoritmo falliranno statisticamente, ripeteremo l'algoritmo finché non sarà trovato almeno un fattore di $N$.

La cella sottostante ripete l'algoritmo finché non viene trovato almeno un fattore di $N=15$. Utilizzeremo i risultati dell'esecuzione hardware di cui sopra per indovinare la fase e il fattore corrispondente in ogni iterazione.

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

Output:

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

Discussione

Lavoro correlato

In questa sezione, discutiamo altre pietre miliari che hanno dimostrato l'algoritmo di Shor su hardware reale.

Il lavoro fondamentale [3] di IBM® ha dimostrato per la prima volta l'algoritmo di Shor, fattorizzando il numero 15 nei suoi fattori primi 3 e 5 utilizzando un computer quantistico a risonanza magnetica nucleare (NMR) a sette qubit. Un altro esperimento [4] ha fatto la fattorizzazione di 15 utilizzando qubit fotonici. Impiegando un singolo qubit riciclato più volte e codificando il registro di lavoro in stati a più alta dimensione, i ricercatori hanno ridotto il numero di qubit richiesto a un terzo di quello del protocollo standard, utilizzando un algoritmo compilato a due fotoni. Un documento significativo nella dimostrazione dell'algoritmo di Shor è [5], che utilizza la tecnica di stima iterativa della fase di Kitaev [8] per ridurre il requisito di qubit dell'algoritmo. Gli autori hanno utilizzato sette qubit di controllo e quattro qubit di cache, insieme all'implementazione di moltiplicatori modulari. Questa implementazione, tuttavia, richiede misure a metà circuito con operazioni di feed-forward e riciclo dei qubit con operazioni di reset. Questa dimostrazione è stata effettuata su un computer quantistico a trappola ionica.

Un lavoro più recente [6] si è concentrato sulla fattorizzazione di 15, 21 e 35 sull'hardware IBM Quantum®. Analogamente a quanto fatto in precedenza, i ricercatori hanno utilizzato una versione compilata dell'algoritmo che impiega una trasformata di Fourier quantistica semiclassica, come proposto da Kitaev, per ridurre al minimo il numero di qubit e porte fisiche. Un lavoro più recente [7] ha anche eseguito una dimostrazione proof-of-concept per la fattorizzazione del numero intero 21. Questa dimostrazione prevedeva anche l'uso di una versione compilata della routine di stima della fase quantistica e si basava sulla precedente dimostrazione di [4]. Gli autori sono andati oltre questo lavoro utilizzando una configurazione di porte Toffoli approssimate con sfasamenti residui. L'algoritmo è stato implementato sui processori quantistici IBM utilizzando solo cinque qubit e la presenza di entanglement tra i qubit di controllo e di registro è stata verificata con successo.

Scalabilità dell'algoritmo

Si noti che la crittografia RSA comporta in genere chiavi di dimensioni dell'ordine di 2048-4096 bit. Il tentativo di fattorizzare un numero di 2048 bit con l'algoritmo di Shor comporterà un circuito quantistico con milioni di qubit, compreso l'overhead di correzione degli errori e una profondità del circuito dell'ordine del miliardo, che è al di là dei limiti di esecuzione dell'attuale hardware quantistico. Pertanto, l'algoritmo di Shor richiederà metodi ottimizzati di costruzione dei circuiti o una robusta correzione quantistica degli errori per essere praticamente praticabile per la violazione dei moderni sistemi crittografici. Rimandiamo a [9] per una discussione più dettagliata sulla stima delle risorse per l'algoritmo di Shor.

Sfida

Congratulazioni per aver terminato il tutorial! Questo è un ottimo momento per verificare la vostra comprensione. Si può provare a costruire il circuito per la fattorizzazione di 21? È possibile selezionare un sito $a$ di propria scelta. Dovrete decidere la precisione in bit dell'algoritmo per scegliere il numero di qubit e dovrete progettare gli operatori di esponenziazione modulare $M_a$. Vi invitiamo a provare voi stessi e a leggere le metodologie illustrate nella Fig. 9 di [6] e nella Fig. 2 di [7].

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

    # Your code here
    pass

Riferimenti

1. Shor, Peter W. "Algoritmi in tempo polinomiale per la fattorizzazione dei primi e i logaritmi discreti su un computer quantistico " Rassegna SIAM 41.2 (1999): 303-332.
2. IBM Quantum Corso "Fondamenti degli algoritmi quantistici" tenuto dal dott. John Watrous.
3. Vandersypen, Lieven MK, et al. "Realizzazione sperimentale dell'algoritmo di fattorizzazione quantistica di Shor mediante risonanza magnetica nucleare " Nature 414.6866 (2001): 883-887.
4. Martin-Lopez, Enrique, et al. "Realizzazione sperimentale dell'algoritmo di fattorizzazione quantistica di Shor utilizzando il riciclo dei qubit " Nature photonics 6.11 (2012): 773-776.
5. Monz, Thomas, et al. "Realizzazione di un algoritmo Shor scalabile " Science 351.6277 (2016): 1068-1070.
6. Amico, Mirko, Zain H. Saleem e Muir Kumph. "Studio sperimentale dell'algoritmo di fattorizzazione di Shor utilizzando l'esperienza di IBM Q " Physical Review A 100.1 (2019): 012305.
7. Skosana, Unathi e Mark Tame. "Dimostrazione dell'algoritmo di fattorizzazione di Shor per N=21 su processori quantistici IBM " Rapporti scientifici 11.1 (2021): 16599.
8. Kitaev, A. Yu. "Misure quantistiche e problema dello stabilizzatore abeliano " arXiv preprint quant-ph/9511026 (1995).
9. Gidney, Craig e Martin Ekerå. "Come fattorizzare interi RSA a 2048 bit in 8 ore usando 20 milioni di qubit rumorosi " Quantum 5 (2021): 433.