Algorithme de Shor

Estimation de l'utilisation : Trois secondes sur un processeur Eagle r3 (NOTE : Il s'agit uniquement d'une estimation. Votre durée d'exécution peut varier.)

Résultats d'apprentissage

À l'issue de ce tutoriel, les utilisateurs devraient avoir compris :

• Les fondements mathématiques de l'algorithme de Shor pour la factorisation des nombres entiers
• Comment exécuter un exemple de cet algorithme sur du matériel informatique

Prérequis

Nous recommandons aux utilisateurs de se familiariser avec les sujets suivants avant de suivre ce tutoriel :

• Principes fondamentaux des algorithmes quantiques.
• Estimation de phase et factorisation. Nous abordons une partie de ces notions dans ce tutoriel.

Arrière-plan

L'algorithme de Shor, mis au point par Peter Shor en 1994, est un algorithme quantique révolutionnaire permettant de factoriser des nombres entiers en temps polynomial. Son importance réside dans sa capacité à factoriser de grands nombres entiers à une vitesse exponentiellement supérieure à celle de tout algorithme classique connu, ce qui menace la sécurité de systèmes cryptographiques largement utilisés, tels que RSA, qui reposent sur la difficulté de factoriser de grands nombres. En résolvant efficacement ce problème sur un ordinateur quantique suffisamment puissant, l'algorithme de Shor pourrait révolutionner des domaines tels que la cryptographie, la cybersécurité et les mathématiques computationnelles, soulignant ainsi le pouvoir transformateur de l'informatique quantique.

Ce tutoriel se concentre sur la démonstration de l'algorithme de Shor en factorisant 15 sur un ordinateur quantique.

Tout d'abord, nous définissons le problème de recherche d'ordre et construisons les circuits correspondants à partir du protocole d'estimation de phase quantique. Ensuite, nous exécutons les circuits de recherche d'ordre sur du matériel réel en utilisant les circuits les plus courts que nous pouvons transpiler. La dernière section complète l'algorithme de Shor en reliant le problème de la recherche d'ordre à la factorisation des nombres entiers.

Nous terminons le tutoriel par une discussion sur d'autres démonstrations de l'algorithme de Shor sur du matériel réel, en nous concentrant à la fois sur les implémentations génériques et celles adaptées à la factorisation d'entiers spécifiques tels que 15 et 21.

Note : Ce tutoriel se concentre davantage sur l'implémentation et la démonstration des circuits concernant l'algorithme de Shor. Pour une formation approfondie sur le sujet, veuillez vous référer au cours Fundamentals of quantum algorithms du Dr. John Watrous, ainsi qu'aux articles figurant dans la section Références.

Exigences

Avant de commencer ce tutoriel, assurez-vous que les éléments suivants sont installés :

• Qiskit SDK v2.0 ou plus tard, avec prise en charge de la visualisation
• Qiskit Runtime v0.40 ou plus tard (pip install qiskit-ibm-runtime)

Configuration

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

Étape 1 : Mettre en correspondance les entrées classiques avec un problème quantique

L'algorithme de Shor pour la factorisation des nombres entiers utilise un problème intermédiaire connu sous le nom de problème de recherche d'ordre. Dans cette section, nous démontrons comment résoudre le problème de recherche d'ordre à l'aide de l' estimation quantique de la phase.

Problème d'estimation de phase

Dans le problème de l'estimation de phase, on nous donne un état quantique $\ket{\psi}$ de $n$ qubits, ainsi qu'un circuit quantique unitaire qui agit sur $n$ qubits. On nous promet que $\ket{\psi}$ est un vecteur propre de la matrice unitaire $U$ qui décrit l'action du circuit, et notre objectif est de calculer ou d'approcher la valeur propre $\lambda = e^{2 \pi i \theta}$ à laquelle correspond $\ket{\psi}$. En d'autres termes, le circuit doit produire une approximation du nombre $\theta \in [0, 1)$ satisfaisante $U \ket{\psi}= e^{2 \pi i \theta} \ket{\psi}.$ L'objectif du circuit d'estimation de phase est d'obtenir une approximation de $\theta$ dans $m$ bits. Mathématiquement parlant, nous aimerions trouver $y$ tel que $\theta \approx y / 2^m$, où $y \in {0, 1, 2, \dots, 2^{m-1}}$. L'image suivante montre le circuit quantique qui estime $y$ en $m$ bits en effectuant une mesure sur $m$ qubits.

Dans le circuit ci-dessus, les qubits supérieurs $m$ sont initiés dans l'état $\ket{0^m}$ et les qubits inférieurs $n$ sont initiés dans l'état $\ket{\psi}$, qui est promis à être un vecteur propre de $U$. Le premier ingrédient du circuit d'estimation de la phase sont les opérations unitaires contrôlées qui sont responsables de l'exécution d'un rebond de phase sur leur qubit de contrôle correspondant. Ces unités contrôlées sont exponentielles en fonction de la position du qubit de contrôle, du bit le moins significatif au bit le plus significatif. Comme $\ket{\psi}$ est un vecteur propre de $U$, l'état des qubits inférieurs de $n$ n'est pas affecté par cette opération, mais l'information de phase de la valeur propre se propage aux qubits supérieurs de $m$.

Il s'avère qu'après l'opération de rebond de phase via les unités contrôlées, tous les états possibles des qubits supérieurs $m$ sont orthonormés les uns par rapport aux autres pour chaque vecteur propre $\ket{\psi}$ de l'unité $U$. Par conséquent, ces états sont parfaitement distinguables et nous pouvons faire pivoter la base qu'ils forment vers la base de calcul pour effectuer une mesure. Une analyse mathématique montre que cette matrice de rotation correspond à la transformée de Fourier quantique inverse (QFT) dans l'espace de Hilbert $2^m$ -dimensionnel. L'intuition sous-jacente est que la structure périodique des opérateurs d'exponentiation modulaire est encodée dans l'état quantique, et que la QFT convertit cette périodicité en pics mesurables dans le domaine des fréquences.

Pour une compréhension plus approfondie de la raison pour laquelle le circuit QFT est utilisé dans l'algorithme de Shor, nous renvoyons le lecteur au cours Fundamentals of quantum algorithms (fondements des algorithmes quantiques ).

Nous sommes maintenant prêts à utiliser le circuit d'estimation de phase pour la recherche d'ordre.

Problème de recherche de commande

Pour définir le problème de la recherche d'ordre, nous commençons par quelques concepts de la théorie des nombres. Tout d'abord, pour tout entier positif donné $N$, définissez l'ensemble $\mathbb{Z}_N$ comme suit $\mathbb{Z}_N = \{0, 1, 2, \dots, N-1\}.$ Toutes les opérations arithmétiques dans $\mathbb{Z}_N$ sont effectuées modulo $N$. En particulier, tous les éléments $a \in \mathbb{Z}_n$ qui sont coprimes avec $N$ sont spéciaux et constituent $\mathbb{Z}^*_N$ sous la forme suivante $\mathbb{Z}^*_N = \{ a \in \mathbb{Z}_N : \mathrm{gcd}(a, N)=1 \}.$ Pour un élément $a \in \mathbb{Z}^*_N$, le plus petit entier positif $r$ tel que $a^r \equiv 1 \; (\mathrm{mod} \; N)$ est défini comme l' ordre de $a$ modulo $N$. Comme nous le verrons plus loin, trouver l'ordre d'un élément $a \in \mathbb{Z}^*_N$ nous permettra de factoriser $N$.

Pour construire le circuit de recherche d'ordre à partir du circuit d'estimation de phase, nous avons besoin de deux considérations. Premièrement, nous devons définir l'unité $U$ qui nous permettra de trouver l'ordre $r$, et deuxièmement, nous devons définir un vecteur propre $\ket{\psi}$ de $U$ pour préparer l'état initial du circuit d'estimation de phase.

Pour relier le problème de recherche d'ordre à l'estimation de phase, nous considérons l'opération définie sur un système dont les états classiques correspondent à $\mathbb{Z}_N$, que nous multiplions par un élément fixe $a \in \mathbb{Z}^*_N$. En particulier, nous définissons cet opérateur de multiplication $M_a$ tel que $M_a \ket{x} = \ket{ax \; (\mathrm{mod} \; N)}$ pour chaque $x \in \mathbb{Z}_N$. Notez qu'il est implicite que nous prenons le produit modulo $N$ à l'intérieur du ket du côté droit de l'équation. Une analyse mathématique montre que $M_a$ est un opérateur unitaire. En outre, il s'avère que $M_a$ possède des paires de vecteurs propres et de valeurs propres qui nous permettent de relier l'ordre $r$ de $a$ au problème de l'estimation de la phase. Plus précisément, pour tout choix de $j \in \{0, \dots, r-1\}$, nous avons que $\ket{\psi_j} = \frac{1}{\sqrt{r}} \sum^{r-1}_{k=0} \omega^{-jk}_{r} \ket{a^k}$ est un vecteur propre de $M_a$ dont la valeur propre correspondante est $\omega^{j}_{r}$, où $\omega^{j}_{r} = e^{2 \pi i \frac{j}{r}}.$

Par observation, nous voyons qu'une paire de vecteurs propres/valeurs propres pratique est l'état $\ket{\psi_1}$ avec $\omega^{1}_{r} = e^{2 \pi i \frac{1}{r}}$. Par conséquent, si nous pouvions trouver le vecteur propre $\ket{\psi_1}$, nous pourrions estimer la phase $\theta=1/r$ avec notre circuit quantique et donc obtenir une estimation de l'ordre $r$. Cependant, ce n'est pas facile à faire, et nous devons envisager une alternative.

Examinons le résultat du circuit si nous préparons l'état de calcul $\ket{1}$ comme état initial. Il ne s'agit pas d'un état propre de $M_a$, mais de la superposition uniforme des états propres que nous venons de décrire. En d'autres termes, la relation suivante s'applique. $\ket{1} = \frac{1}{\sqrt{r}} \sum^{r-1}_{k=0} \ket{\psi_k}$ L'implication de l'équation ci-dessus est que si nous fixons l'état initial à $\ket{1}$, nous obtiendrons exactement le même résultat de mesure que si nous avions choisi $k \in \{ 0, \dots, r-1\}$ uniformément au hasard et utilisé $\ket{\psi_k}$ comme vecteur propre dans le circuit d'estimation de la phase. En d'autres termes, une mesure des $m$ qubits supérieurs donne une approximation $y / 2^m$ de la valeur $k / r$ où $k \in \{ 0, \dots, r-1\}$ est choisi uniformément au hasard. Cela nous permet d'apprendre $r$ avec un haut degré de confiance après plusieurs essais indépendants, ce qui était notre objectif.

Opérateurs d'exponentiation modulaires

Jusqu'à présent, nous avons lié le problème de l'estimation de la phase au problème de la recherche d'ordre en définissant $U = M_a$ et $\ket{\psi} = \ket{1}$ dans notre circuit quantique. Par conséquent, le dernier ingrédient restant consiste à trouver un moyen efficace de définir les exponentielles modulaires de $M_a$ comme $M_a^k$ pour $k = 1, 2, 4, \dots, 2^{m-1}$. Pour effectuer ce calcul, nous constatons que pour toute puissance $k$ que nous choisissons, nous pouvons créer un circuit pour $M_a^k$ non pas en itérant $k$ fois le circuit pour $M_a$, mais plutôt en calculant $b = a^k \; \mathrm{mod} \; N$ et en utilisant ensuite le circuit pour $M_b$. Comme nous n'avons besoin que des puissances qui sont elles-mêmes des puissances de 2, nous pouvons le faire de manière classiquement efficace en utilisant la quadrature itérative.

Étape 2 : Optimiser le problème pour l'exécution sur du matériel quantique

Exemple concret avec $N = 15$ et $a=2$

Nous pouvons nous arrêter ici pour discuter d'un exemple spécifique et construire le circuit de recherche d'ordre pour $N=15$. Notez que les $a \in \mathbb{Z}_N^*$ non triviaux possibles pour $N=15$ sont $a \in \{2, 4, 7, 8, 11, 13, 14 \}$. Pour cet exemple, nous choisissons $a=2$. Nous construirons l'opérateur $M_2$ et les opérateurs d'exponentiation modulaire $M_2^k$.

L'action de $M_2$ sur les états de base de calcul est la suivante. $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}$ Par observation, nous pouvons voir que les états de base sont mélangés, nous avons donc une matrice de permutation. Nous pouvons réaliser cette opération sur quatre qubits à l'aide de portes de permutation. Ci-dessous, nous construisons les opérations $M_2$ et $M_2$ contrôlées.

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:

Les portes agissant sur plus de deux qubits seront décomposées en portes à deux qubits.

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

Output:

Nous devons maintenant construire les opérateurs d'exponentiation modulaire. Pour obtenir une précision suffisante dans l'estimation de la phase, nous utiliserons huit qubits pour la mesure de l'estimation. Par conséquent, nous devons construire $M_b$ avec $b = a^{2^k} \; (\mathrm{mod} \; N)$ pour chaque $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]

Comme nous pouvons le voir dans la liste des valeurs de $b$, en plus de $M_2$ que nous avons construit précédemment, nous devons également construire $M_4$ et $M_1$. Notez que $M_1$ agit trivialement sur les états de la base de calcul, il s'agit donc simplement de l'opérateur d'identité.

$M_4$ agit sur les états de base de calcul comme suit. $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}$

Par conséquent, cette permutation peut être construite avec l'opération de permutation suivante.

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:

Les portes agissant sur plus de deux qubits seront décomposées en portes à deux qubits.

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

Output:

Nous avons vu que les opérateurs $M_b$ pour un $b \in \mathbb{Z}^*_N$ donné sont des opérations de permutation. En raison de la taille relativement petite du problème de permutation que nous avons ici, puisque $N=15$ ne nécessite que quatre qubits, nous avons pu synthétiser ces opérations directement avec des portes SWAP par inspection. En général, il ne s'agit pas d'une approche évolutive. Au lieu de cela, nous pourrions avoir besoin de construire la matrice de permutation explicitement, et d'utiliser la classe UnitaryGate et les méthodes de transpilation de Qiskit pour synthétiser cette matrice de permutation. Toutefois, cela peut entraîner des circuits beaucoup plus profonds. En voici un exemple.

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

Comparons ces chiffres avec la profondeur du circuit compilé de notre implémentation manuelle de la porte $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})

Comme nous pouvons le constater, l'approche de la matrice de permutation a donné lieu à un circuit nettement plus profond, même pour une seule porte $M_2$, par rapport à notre mise en œuvre manuelle. Par conséquent, nous poursuivrons la mise en œuvre des opérations $M_b$.

Nous sommes maintenant prêts à construire le circuit de recherche d'ordre complet en utilisant les opérateurs d'exponentiation modulaire contrôlée que nous avons définis précédemment. Dans le code suivant, nous importons également le circuit QFT de la bibliothèque Qiskit Circuit, qui utilise des portes de Hadamard sur chaque qubit, une série de portes controlled-U1 (ou Z, selon la phase) et une couche de portes de permutation.

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

Notez que nous avons omis les opérations d'exponentiation modulaire contrôlée pour les qubits de contrôle restants, car $M_1$ est l'opérateur d'identité.

Notez que plus loin dans ce tutoriel, nous exécuterons ce circuit sur le backend ibm_marrakesh . Pour ce faire, nous transposons le circuit en fonction de ce backend spécifique et indiquons la profondeur du circuit et le nombre de portes.

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

Étape 3 : Exécutez à l'aide d' Qiskit primitives

Tout d'abord, nous discutons de ce que nous obtiendrions théoriquement si nous faisions fonctionner ce circuit sur un simulateur idéal. Nous présentons ci-dessous une série de résultats de simulation du circuit ci-dessus en utilisant 1024 plans. Comme nous pouvons le constater, nous obtenons une distribution approximativement uniforme sur quatre chaînes de bits pour les qubits de contrôle.

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

plot_histogram(counts)

Output:

En mesurant les qubits de contrôle, nous obtenons une estimation de phase sur huit bits de l'opérateur $M_a$. Nous pouvons convertir cette représentation binaire en décimale pour trouver la phase mesurée. Comme le montre l'histogramme ci-dessus, quatre chaînes de bits différentes ont été mesurées et chacune d'entre elles correspond à une valeur de phase comme suit.

# 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

Rappelons que la phase mesurée quelconque correspond à $\theta = k / r$ où $k$ est échantillonné uniformément au hasard à partir de $\{0, 1, \dots, r-1 \}$. Par conséquent, nous pouvons utiliser l'algorithme des fractions continues pour tenter de trouver $k$ et l'ordre $r$. Python a cette fonctionnalité intégrée. Nous pouvons utiliser le module fractions pour transformer un flotteur en un objet Fraction , par exemple :

Fraction(0.666)

Output:

Fraction(5998794703657501, 9007199254740992)

Parce que cela donne des fractions qui retournent exactement le résultat (dans ce cas, 0.6660000...), cela peut donner des résultats bizarres comme celui ci-dessus. Nous pouvons utiliser la méthode .limit_denominator() pour obtenir la fraction qui ressemble le plus à notre flotteur, avec un dénominateur inférieur à une certaine valeur :

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

Output:

Fraction(2, 3)

C'est beaucoup plus beau. L'ordre (r) doit être inférieur à N, nous fixerons donc le dénominateur maximal à 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

Nous pouvons constater que deux des valeurs propres mesurées nous ont fourni le bon résultat : $r=4$ et nous pouvons voir que l'algorithme de Shor pour la recherche d'ordre a une chance d'échouer. Ces mauvais résultats sont dus au fait que $k = 0$, ou que $k$ et $r$ ne sont pas coprimes - et qu'au lieu de $r$, on nous donne un facteur de $r$. La solution la plus simple consiste à répéter l'expérience jusqu'à ce que nous obtenions un résultat satisfaisant pour $r$.

Jusqu'à présent, nous avons mis en œuvre le problème de recherche d'ordre pour $N=15$ avec $a=2$ en utilisant le circuit d'estimation de phase sur un simulateur. La dernière étape de l'algorithme de Shor consistera à relier le problème de la recherche d'ordre au problème de la factorisation des nombres entiers. Cette dernière partie de l'algorithme est purement classique et peut être résolue sur un ordinateur classique après que les mesures de phase ont été obtenues à partir d'un ordinateur quantique. Par conséquent, nous reportons la dernière partie de l'algorithme jusqu'à ce que nous ayons démontré comment nous pouvons faire fonctionner le circuit de recherche d'ordre sur du matériel réel.

Exécution du matériel

Nous pouvons maintenant exécuter le circuit de recherche d'ordre que nous avons transposé précédemment pour ibm_marrakesh. Nous nous tournons ici vers le découplage dynamique (DD) pour supprimer les erreurs, et vers le tourbillonnement des portes pour atténuer les erreurs. Le DD consiste à appliquer à un dispositif quantique des séquences d'impulsions de contrôle chronométrées avec précision, ce qui permet d'éliminer les interactions environnementales indésirables et la décohérence. Le tournoiement de portes, quant à lui, randomise des portes quantiques spécifiques afin de transformer les erreurs cohérentes en erreurs de Pauli, qui s'accumulent de manière linéaire plutôt que quadratique. Ces deux techniques sont souvent combinées pour améliorer la cohérence et la fidélité des calculs quantiques.

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

Comme nous pouvons le constater, nous obtenons les mêmes chaînes de bits avec les nombres les plus élevés. Le matériel quantique étant bruyant, il y a des fuites vers d'autres chaînes de bits, que nous pouvons filtrer statistiquement.

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

Étape 4 : Post-traitement et restitution du résultat dans le format classique souhaité

Factorisation des nombres entiers

Jusqu'à présent, nous avons examiné comment nous pouvons mettre en œuvre le problème de recherche d'ordre à l'aide d'un circuit d'estimation de phase. Maintenant, nous relions le problème de la recherche d'ordre à la factorisation des nombres entiers, ce qui complète l'algorithme de Shor. Notez que cette partie de l'algorithme est classique.

Nous allons maintenant le démontrer en utilisant notre exemple de $N = 15$ et $a = 2$. Rappelons que la phase que nous avons mesurée est $k / r$, où $a^r \; (\textrm{mod} \; N) = 1$ et $k$ sont des nombres entiers aléatoires compris entre $0$ et $r - 1$. Cette équation nous donne $(a^r - 1) \; (\textrm{mod} \; N) = 0,$, ce qui signifie que $N$ doit diviser $a^r-1$. Si $r$ est également pair, nous pouvons écrire $a^r -1 = (a^{r/2}-1)(a^{r/2}+1).$ Si $r$ n'est pas pair, nous ne pouvons pas aller plus loin et devons réessayer avec une valeur différente pour $a$; sinon, il y a une forte probabilité que le plus grand commun diviseur de $N$ et soit $a^{r/2}-1$, soit $a^{r/2}+1$ soit un facteur propre de $N$.

Étant donné que certaines exécutions de l'algorithme échoueront statistiquement, nous répéterons cet algorithme jusqu'à ce qu'au moins un facteur de $N$ soit trouvé.

La cellule ci-dessous répète l'algorithme jusqu'à ce qu'au moins un facteur de $N=15$ soit trouvé. Nous utiliserons les résultats de l'exécution matérielle ci-dessus pour deviner la phase et le facteur correspondant à chaque itération.

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

La discussion

Travaux connexes

Dans cette section, nous examinons d'autres travaux importants qui ont démontré l'algorithme de Shor sur du matériel réel.

Le travail fondamental [3] de IBM® a démontré l'algorithme de Shor pour la première fois, en factorisant le nombre 15 en ses facteurs premiers 3 et 5 à l'aide d'un ordinateur quantique à résonance magnétique nucléaire (RMN) à sept qubits. Une autre expérience [4] a factorisé 15 en utilisant des qubits photoniques. En utilisant un seul qubit recyclé plusieurs fois et en codant le registre de travail dans des états de dimension supérieure, les chercheurs ont réduit le nombre de qubits requis à un tiers de celui du protocole standard, en utilisant un algorithme compilé à deux photons. Un article important dans la démonstration de l'algorithme de Shor est [5], qui utilise la technique d'estimation de phase itérative de Kitaev [8] pour réduire le nombre de qubits requis par l'algorithme. Les auteurs ont utilisé sept qubits de contrôle et quatre qubits de cache, ainsi que la mise en œuvre de multiplicateurs modulaires. Cette mise en œuvre nécessite toutefois des mesures en milieu de circuit avec des opérations de feed-forward et un recyclage des qubits avec des opérations de reset. Cette démonstration a été réalisée sur un ordinateur quantique à piège à ions.

Des travaux plus récents [6] ont porté sur la factorisation de 15, 21 et 35 sur le matériel Quantum® de IBM. Comme pour les travaux précédents, les chercheurs ont utilisé une version compilée de l'algorithme qui utilise une transformée de Fourier quantique semi-classique telle que proposée par Kitaev pour minimiser le nombre de qubits et de portes physiques. Un travail plus récent [7] a également effectué une démonstration de principe pour la factorisation de l'entier 21. Cette démonstration impliquait également l'utilisation d'une version compilée de la routine d'estimation de la phase quantique, et s'appuyait sur la démonstration précédente de [4]. Les auteurs sont allés plus loin en utilisant une configuration de portes de Toffoli approximatives avec des déphasages résiduels. L'algorithme a été mis en œuvre sur les processeurs quantiques IBM en utilisant seulement cinq qubits, et la présence d'intrication entre les qubits de contrôle et de registre a été vérifiée avec succès.

Mise à l'échelle de l'algorithme

Nous notons que le cryptage RSA implique généralement des tailles de clés de l'ordre de 2048 à 4096 bits. Tenter de factoriser un nombre de 2048 bits avec l'algorithme de Shor aboutira à un circuit quantique avec des millions de qubits, y compris la correction d'erreur et une profondeur de circuit de l'ordre du milliard, ce qui dépasse les limites d'exécution du matériel quantique actuel. Par conséquent, l'algorithme de Shor nécessitera soit des méthodes de construction de circuits optimisées, soit une correction d'erreur quantique robuste pour être viable dans la pratique afin de casser les systèmes cryptographiques modernes. Nous vous renvoyons à [9] pour une discussion plus détaillée sur l'estimation des ressources pour l'algorithme de Shor.

Le défi

Félicitations pour avoir terminé le tutoriel! C'est le moment idéal pour tester votre compréhension. Pourriez-vous essayer de construire le circuit pour factoriser 21? Vous pouvez sélectionner un site $a$ de votre choix. Vous devrez décider de la précision des bits de l'algorithme pour choisir le nombre de qubits, et vous devrez concevoir les opérateurs d'exponentiation modulaire $M_a$. Nous vous encourageons à essayer par vous-même, puis à lire les méthodologies présentées dans la Fig. 9 de [6] et la Fig. 2 de [7].

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

    # Your code here
    pass

Références

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 Cours « Fondements des algorithmes quantiques » dispensé par le 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, et Muir Kumph. " Étude expérimentale de l'algorithme de factorisation de Shor à l'aide de l'expérience IBM Q. " Physical Review A 100.1 (2019) : 012305.
7. Skosana, Unathi, et Mark Tame. " Démonstration de l'algorithme de factorisation de Shor pour N=21 sur les processeurs quantiques IBM " Rapports scientifiques 11.1 (2021) : 16599.
8. Kitaev, A. Yu. " Mesures quantiques et problème du stabilisateur abélien " arXiv preprint quant-ph/9511026 (1995).
9. Gidney, Craig, et Martin Ekerå. " Comment factoriser des entiers RSA de 2048 bits en 8 heures en utilisant 20 millions de qubits bruyants Quantum 5 (2021) : 433.