Atténuation des erreurs de lecture pour la primitive Sampler avec M3

Estimation d'utilisation : moins d'une minute sur un processeur Heron r2 (REMARQUE : il s'agit uniquement d'une estimation. Votre temps d'exécution peut varier.)

Contexte

Contrairement à la primitive Estimator, la primitive Sampler ne dispose pas d'un support intégré pour l'atténuation des erreurs. Plusieurs des méthodes prises en charge par l'Estimator sont spécifiquement conçues pour les valeurs d'espérance, et ne sont donc pas applicables à la primitive Sampler. Une exception est l'atténuation des erreurs de lecture, qui est une méthode très efficace également applicable à la primitive Sampler.

L'addon Qiskit M3 implémente une méthode efficace pour l'atténuation des erreurs de lecture. Ce tutoriel explique comment utiliser l'addon Qiskit M3 pour atténuer les erreurs de lecture de la primitive Sampler.

Qu'est-ce qu'une erreur de lecture ?

Immédiatement avant la mesure, l'état d'un registre de qubits est décrit par une superposition d'états de base computationnels, ou par une matrice densité. La mesure du registre de qubits dans un registre de bits classiques se déroule alors en deux étapes. Premièrement, la mesure quantique proprement dite est effectuée. Cela signifie que l'état du registre de qubits est projeté sur un unique état de base caractérisé par une chaîne de $1$ et de $0$. La seconde étape consiste à lire la chaîne de bits caractérisant cet état de base et à l'écrire dans la mémoire d'un ordinateur classique. Nous appelons cette étape la lecture (readout). Il s'avère que la seconde étape (lecture) engendre plus d'erreurs que la première étape (projection sur les états de base). Cela s'explique aisément si vous vous rappelez que la lecture nécessite de détecter un état quantique microscopique et de l'amplifier au niveau macroscopique. Un résonateur de lecture est couplé au qubit (transmon), subissant ainsi un très faible décalage de fréquence. Une impulsion micro-onde est alors réfléchie par le résonateur, subissant à son tour de petites modifications de ses caractéristiques. L'impulsion réfléchie est ensuite amplifiée et analysée. Il s'agit d'un processus délicat sujet à de nombreuses erreurs.

Le point important est que, bien que la mesure quantique et la lecture soient toutes deux sujettes à des erreurs, cette dernière engendre l'erreur dominante, appelée erreur de lecture, qui est l'objet de ce tutoriel.

Contexte théorique

Si la chaîne de bits échantillonnée (stockée dans la mémoire classique) diffère de la chaîne de bits caractérisant l'état quantique projeté, nous disons qu'une erreur de lecture s'est produite. Ces erreurs sont observées comme étant aléatoires et non corrélées d'un échantillon à l'autre. Il s'est avéré utile de modéliser l'erreur de lecture comme un canal classique bruité. C'est-à-dire que pour chaque paire de chaînes de bits $i$ et $j$, il existe une probabilité fixe que la vraie valeur $j$ soit incorrectement lue comme $i$.

Plus précisément, pour chaque paire de chaînes de bits $(i, j)$, il existe une probabilité (conditionnelle) ${M}_{i,j}$ que $i$ soit lu, étant donné que la vraie valeur est $j.$ C'est-à-dire,

$${M}_{i,j} = \Pr(\text{readout value is } i | \text{true value is } j) \text{ for } i,j \in (0,...,2^n - 1), \tag{1}$$

où $n$ est le nombre de bits dans le registre de lecture. Pour être concret, nous supposons que $i$ est un entier décimal dont la représentation binaire est la chaîne de bits qui étiquette les états de base computationnels. Nous appelons la matrice $2^n \times 2^n$ ${M}$ la matrice d'assignation. Pour une valeur vraie $j$ fixée, la somme de la probabilité sur tous les résultats bruités $i$ doit donner $1$. C'est-à-dire

$$\sum_{i=0}^{2^n - 1} {M}_{i,j} = 1 \text{ for all } j$$

Une matrice sans entrées négatives qui satisfait (1) est dite stochastique à gauche. Une matrice stochastique à gauche est aussi appelée stochastique par colonnes car chacune de ses colonnes somme à $1$. Nous déterminons expérimentalement les valeurs approximatives de chaque élément ${M}_{i,j}$ en préparant de manière répétée chaque état de base $|j \rangle$ puis en calculant les fréquences d'occurrence des chaînes de bits échantillonnées.

Si une expérience consiste à estimer une distribution de probabilité sur les chaînes de bits de sortie par échantillonnage répété, alors nous pouvons utiliser ${M}$ pour atténuer l'erreur de lecture au niveau de la distribution. La première étape consiste à répéter un circuit fixe d'intérêt de nombreuses fois, créant un histogramme des chaînes de bits échantillonnées. L'histogramme normalisé est la distribution de probabilité mesurée sur les $2^n$ chaînes de bits possibles, que nous notons ${\tilde{p}} \in \mathbb{R}^{2^n}$. La probabilité (estimée) ${{\tilde{p}}}_i$ d'échantillonner la chaîne de bits $i$ est égale à la somme sur toutes les vraies chaînes de bits $j$, chacune pondérée par la probabilité qu'elle soit confondue avec $i$. Cette relation sous forme matricielle est

$${\tilde{p}} = {M} {\vec{p}}, \tag{2},$$

où ${\vec{p}}$ est la distribution vraie. En d'autres termes, l'erreur de lecture a pour effet de multiplier la distribution idéale sur les chaînes de bits ${\vec{p}}$ par la matrice d'assignation ${M}$ pour produire la distribution observée ${\tilde{p}}$. Nous avons mesuré ${\tilde{p}}$ et ${M}$, mais n'avons pas d'accès direct à ${\vec{p}}$. En principe, nous obtiendrons la vraie distribution de chaînes de bits pour notre circuit en résolvant numériquement l'équation (2) pour ${\vec{p}}$.

Avant de poursuivre, il convient de noter quelques caractéristiques importantes de cette approche naïve.

• En pratique, l'équation (2) n'est pas résolue en inversant ${M}$. Les routines d'algèbre linéaire des bibliothèques logicielles emploient des méthodes plus stables, précises et efficaces.
• Lors de l'estimation de ${M}$, nous avons supposé que seules des erreurs de lecture se produisaient. En particulier, nous supposons qu'il n'y avait pas d'erreurs de préparation d'état ni de mesure quantique — ou du moins qu'elles étaient atténuées par ailleurs. Dans la mesure où cette hypothèse est valide, ${M}$ représente réellement uniquement l'erreur de lecture. Mais lorsque nous utilisons ${M}$ pour corriger une distribution mesurée sur les chaînes de bits, nous ne faisons pas une telle hypothèse. En fait, nous nous attendons à ce qu'un circuit intéressant introduise du bruit, par exemple des erreurs de portes. La « vraie » distribution inclut toujours les effets de toutes les erreurs qui ne sont pas atténuées par ailleurs.

Cette méthode, bien qu'utile dans certaines circonstances, souffre de quelques limitations.

Les ressources en espace et en temps nécessaires pour estimer ${M}$ croissent exponentiellement en $n$ :

• L'estimation de ${M}$ et ${\tilde{p}}$ est sujette à une erreur statistique due à l'échantillonnage fini. Ce bruit peut être rendu aussi petit que souhaité au prix de plus de shots (jusqu'à l'échelle de temps de dérive des paramètres matériels qui entraînent des erreurs systématiques dans ${M}$). Cependant, si aucune hypothèse n'est faite sur les chaînes de bits observées lors de l'atténuation, le nombre de shots requis pour estimer ${M}$ croît au moins exponentiellement en $n$.
• ${M}$ est une matrice $2^n \times 2^n$. Quand $n>10$, la quantité de mémoire requise pour stocker ${M}$ est supérieure à la mémoire disponible sur un ordinateur portable puissant.

Les limitations supplémentaires sont :

• La distribution récupérée ${\vec{p}}$ peut avoir une ou plusieurs probabilités négatives (tout en sommant à un). Une solution consiste à minimiser $||{M} {\vec{p}} - {\tilde{p}}||^2$ sous la contrainte que chaque entrée de ${\vec{p}}$ soit non négative. Cependant, le temps d'exécution d'une telle méthode est de plusieurs ordres de grandeur supérieur à la résolution directe de l'équation (2).
• Cette procédure d'atténuation fonctionne au niveau d'une distribution de probabilité sur les chaînes de bits. En particulier, elle ne peut pas corriger une erreur dans une chaîne de bits observée individuellement.

Addon Qiskit M3 : passage à l'échelle pour des chaînes de bits plus longues

La résolution de l'équation (2) à l'aide de routines standard d'algèbre linéaire numérique est limitée à des chaînes de bits de longueur maximale d'environ 10 bits. M3, cependant, peut traiter des chaînes de bits beaucoup plus longues. Deux propriétés clés de M3 rendent cela possible :

• Les corrélations d'ordre trois et supérieur dans les erreurs de lecture entre collections de bits sont supposées négligeables et sont ignorées. En principe, au prix de plus de shots, on pourrait estimer également les corrélations d'ordre supérieur.
• Plutôt que de construire ${M}$ explicitement, nous utilisons une matrice effective beaucoup plus petite qui enregistre les probabilités uniquement pour les chaînes de bits collectées lors de la construction de ${\tilde{p}}$.

À un niveau élevé, la procédure fonctionne comme suit.

Premièrement, nous construisons des blocs élémentaires à partir desquels nous pouvons construire une description simplifiée et effective de ${M}$. Ensuite, nous exécutons le circuit d'intérêt de manière répétée et collectons les chaînes de bits que nous utilisons pour construire à la fois ${\tilde{p}}$ et, à l'aide des blocs élémentaires, une matrice ${M}$ effective.

Plus précisément,

• Les matrices d'assignation à un seul qubit sont estimées pour chaque qubit. Pour ce faire, nous préparons de manière répétée le registre de qubits dans l'état tout-zéro $|0 ... 0 \rangle$ puis dans l'état tout-un $|1 ... 1 \rangle$, et nous enregistrons la probabilité pour chaque qubit d'être lu incorrectement.

• Les corrélations d'ordre trois et supérieur sont supposées négligeables et sont ignorées.

  Au lieu de cela, nous construisons un nombre $n$ de matrices d'assignation à un seul qubit de taille $2 \times 2$, et un nombre $n(n-1)/2$ de matrices d'assignation à deux qubits de taille $4 \times 4$. Ces matrices d'assignation à un et deux qubits sont stockées pour une utilisation ultérieure.

• Après avoir échantillonné un circuit de manière répétée pour construire ${\tilde{p}}$, nous construisons une approximation effective de ${M}$ en utilisant uniquement les chaînes de bits échantillonnées lors de la construction de ${\tilde{p}}$. Cette matrice effective est construite à l'aide des matrices à un et deux qubits décrites dans le point précédent. La dimension linéaire de cette matrice est au plus de l'ordre du nombre de shots utilisés dans la construction de ${\tilde{p}}$, ce qui est bien inférieur à la dimension $2^n$ de la matrice d'assignation complète ${M}$.

Pour les détails techniques sur M3, vous pouvez consulter Scalable Mitigation of Measurement Errors on Quantum Computers.

Application de M3 à un algorithme quantique

Nous appliquerons l'atténuation de lecture de M3 au problème du décalage caché. Le problème du décalage caché, et les problèmes étroitement liés tels que le problème du sous-groupe caché, ont été initialement conçus dans un contexte tolérant aux fautes (plus précisément, avant qu'il ne soit prouvé que les QPU tolérants aux fautes étaient possibles !). Mais ils sont également étudiés avec les processeurs disponibles. Un exemple d'accélération algorithmique exponentielle obtenue pour une variante du problème du décalage caché sur des QPU IBM® de 127 qubits peut être trouvé dans cet article (version arXiv).

Dans ce qui suit, toute l'arithmétique est booléenne. C'est-à-dire que pour $a, b \in \mathbb{Z}_2 = \{0, 1\}$, l'addition $a + b$ est la fonction OU exclusif (XOR) logique. De plus, la multiplication $a \times b$ (ou $a b$) est la fonction ET (AND) logique. Pour $x, y \in \{0, 1\}^n$, $x + y$ est défini par application bit à bit du XOR. Le produit scalaire $\cdot: {\mathbb{Z}_2^n} \rightarrow \mathbb{Z}_2$ est défini par $x \cdot y = \sum_i x_i y_i$.

Opérateur de Hadamard et transformée de Fourier

Dans l'implémentation d'algorithmes quantiques, il est très courant d'utiliser l'opérateur de Hadamard comme transformée de Fourier. Les états de base computationnels sont parfois appelés états classiques. Ils sont en correspondance biunivoque avec les chaînes de bits classiques. L'opérateur de Hadamard à $n$ qubits sur les états classiques peut être vu comme une transformée de Fourier sur l'hypercube booléen :

$$H^{\otimes n} = \frac{1}{\sqrt{2^n}} \sum_{x,y \in {\mathbb{Z}_2^n}} (-1)^{x \cdot y} {|{y}\rangle}{\langle{x}|}.$$

Considérons un état ${|{s}\rangle}$ correspondant à une chaîne de bits fixe $s$. En appliquant $H^{\otimes n}$, et en utilisant ${\langle {x}|{s}\rangle} = \delta_{x,s}$, nous voyons que la transformée de Fourier de ${|{s}\rangle}$ peut s'écrire

$$H^{\otimes n} {|{s}\rangle} = \frac{1}{\sqrt{2^n}} \sum_{y \in {\mathbb{Z}_2^n}} (-1)^{s \cdot y} {|{y}\rangle}.$$

L'opérateur de Hadamard est son propre inverse, c'est-à-dire $H^{\otimes n} H^{\otimes n} = (H H)^{\otimes n} = I^{\otimes n}$. Ainsi, la transformée de Fourier inverse est le même opérateur, $H^{\otimes n}$. Explicitement, nous avons,

$${|{s}\rangle} = H^{\otimes n} H^{\otimes n} {|{s}\rangle} = H^{\otimes n} \frac{1}{\sqrt{2^n}} \sum_{y \in {\mathbb{Z}_2^n}} (-1)^{s \cdot y} {|{y}\rangle}.$$

Le problème du décalage caché

Nous considérons un exemple simple d'un problème de décalage caché. Le problème consiste à identifier un décalage constant dans l'entrée d'une fonction. La fonction que nous considérons est le produit scalaire. C'est le membre le plus simple d'une grande classe de fonctions qui admettent une accélération quantique pour le problème du décalage caché via des techniques similaires à celles présentées ci-dessous.

Soient $x,y \in {\mathbb{Z}_2^m}$ des chaînes de bits de longueur $m$. Nous définissons ${f}: {\mathbb{Z}_2^m} \times {\mathbb{Z}_2^m} \rightarrow \{-1,1\}$ par

$${f}(x, y) = (-1)^{x \cdot y}.$$

Soient $a,b \in {\mathbb{Z}_2^m}$ des chaînes de bits fixées de longueur $m$. Nous définissons de plus $g: {\mathbb{Z}_2^m} \times {\mathbb{Z}_2^m} \rightarrow \{-1,1\}$ par

$$g(x, y) = {f}(x+a, y+b) = (-1)^{(x+a) \cdot (y+b)},$$

où $a$ et $b$ sont des paramètres (cachés). On nous donne deux boîtes noires, l'une implémentant $f$, et l'autre $g$. Nous supposons que nous savons qu'elles calculent les fonctions définies ci-dessus, sauf que nous ne connaissons ni $a$ ni $b$. Le jeu consiste à déterminer les chaînes de bits cachées (décalages) $a$ et $b$ en effectuant des requêtes à $f$ et $g$. Il est clair que si nous jouons le jeu de manière classique, nous avons besoin de $O(2m)$ requêtes pour déterminer $a$ et $b$. Par exemple, nous pouvons interroger $g$ avec toutes les paires de chaînes telles que l'un des éléments de la paire est entièrement nul, et l'autre élément a exactement un bit mis à $1$. À chaque requête, nous apprenons un élément de $a$ ou de $b$. Cependant, nous verrons que, si les boîtes noires sont implémentées comme des circuits quantiques, nous pouvons déterminer $a$ et $b$ avec une seule requête à chacune de $f$ et $g$.

Dans le contexte de la complexité algorithmique, une boîte noire est appelée un oracle. En plus d'être opaque, un oracle a la propriété de consommer l'entrée et de produire la sortie instantanément, n'ajoutant rien au budget de complexité de l'algorithme dans lequel il est intégré. En fait, dans le cas présent, les oracles implémentant $f$ et $g$ se révéleront être efficaces.

Circuits quantiques pour $f$ et $g$

Nous avons besoin des ingrédients suivants pour implémenter $f$ et $g$ comme circuits quantiques.

Pour des états classiques à un seul qubit ${|{x_1}\rangle}, {|{y_1}\rangle}$, avec $x_1,y_1 \in \mathbb{Z}_2$, la porte $Z$ contrôlée ${CZ}$ peut s'écrire

$${CZ} {|{x_1}\rangle}{|{y_1}\rangle}{x_1} = (-1)^{x_1 y_1} {|{x_1}\rangle}{x_1}{|{y_1}\rangle}.$$

Nous opérerons avec $m$ portes CZ, une sur $(x_1, y_1)$, une sur $(x_2, y_2)$, et ainsi de suite jusqu'à $(x_m, y_m)$. Nous appelons cet opérateur ${CZ}_{x,y}$.

$U_f = {CZ}_{x,y}$ est une version quantique de ${f} = {f}(x,y)$ :

$$U_f {|{x}\rangle}{|{y}\rangle} = {CZ}_{x,y} {|{x}\rangle}{|{y}\rangle} = (-1)^{x \cdot y} {|{x}\rangle}{|{y}\rangle}.$$

Nous devons également implémenter un décalage de chaîne de bits. Nous notons l'opérateur sur le registre $x$ $X^{a_1}\cdots X^{a_m}$ par $X_a$ et de même sur le registre $y$ $X_b = X^{b_1}\cdots X^{b_m}$. Ces opérateurs appliquent $X$ partout où un bit individuel est $1$, et l'identité $I$ partout où il est $0$. Nous avons alors

$$X_a X_b {|{x}\rangle}{|{y}\rangle} = {|{x+a}\rangle}{|{y+b}\rangle}.$$

La seconde boîte noire $g$ est implémentée par l'unitaire $U_g$, donné par

$$U_g = X_aX_b {CZ}_{x,y} X_aX_b.$$

Pour le voir, nous appliquons les opérateurs de droite à gauche sur l'état ${|{x}\rangle}{|{y}\rangle}$. Premièrement

$$X_a X_b {|{x}\rangle}{|{y}\rangle} = {|{x+a}\rangle}{|{y+b}\rangle}.$$

Puis,

$${CZ}_{x,y} {|{x+a}\rangle}{|{y+b}\rangle} = (-1)^{(x+a)\cdot (y+b)} {|{x+a}\rangle}{|{y+b}\rangle}.$$

Enfin,

$$X^a X^b (-1)^{(x+a)\cdot (y+b)} {|{x+a}\rangle}{|{y+b}\rangle} = (-1)^{(x+a)\cdot (y+b)} {|{x}\rangle}{|{y}\rangle},$$

ce qui est bien la version quantique de $f(x+a, y+b)$.

L'algorithme du décalage caché

Maintenant nous assemblons les pièces pour résoudre le problème du décalage caché. Nous commençons par appliquer des portes de Hadamard aux registres initialisés à l'état tout-zéro.

$$H^{\otimes 2m} = H^{\otimes m} \otimes H^{\otimes m} {{|{0}\rangle}^{\otimes m}}{{|{0}\rangle}^{\otimes m}} = \frac{1}{\sqrt{2^{2m}}} \sum_{x, y \in {\mathbb{Z}_2^m}} (-1)^{x \cdot y} {|{x}\rangle}{|{y}\rangle}.$$

Ensuite, nous interrogeons l'oracle $g$ pour obtenir

$$U_g H^{\otimes 2m} {{|{0}\rangle}^{\otimes m}}{{|{0}\rangle}^{\otimes m}} = \frac{1}{\sqrt{2^{2m}}} \sum_{x, y \in {\mathbb{Z}_2^m}} (-1)^{(x+a) \cdot (y+b)} {|{x}\rangle}{|{y}\rangle}$$ $$\approx \frac{1}{\sqrt{2^{2m}}} \sum_{x, y \in {\mathbb{Z}_2^m}} (-1)^{x \cdot y + x \cdot b + y \cdot a} {|{x}\rangle}{|{y}\rangle}.$$

Dans la dernière ligne, nous avons omis le facteur de phase globale constant $(-1)^{a \cdot b}$, et notons l'égalité à une phase près par $\approx$. Ensuite, l'application de l'oracle $f$ introduit un autre facteur $(-1)^{x \cdot y}$, annulant celui déjà présent. Nous avons alors :

$$U_f U_g H^{\otimes 2m} {{|{0}\rangle}^{\otimes m}}{{|{0}\rangle}^{\otimes m}} \approx \frac{1}{\sqrt{2^{2m}}} \sum_{x, y \in {\mathbb{Z}_2^m}} (-1)^{x \cdot b + y \cdot a} {|{x}\rangle}{|{y}\rangle}.$$

L'étape finale consiste à appliquer la transformée de Fourier inverse, $H^{\otimes 2m} = H^{\otimes m} \otimes H^{\otimes m}$, donnant

$$H^{\otimes 2m} U_f U_g H^{\otimes 2m} {{|{0}\rangle}^{\otimes m}}{{|{0}\rangle}^{\otimes m}} \approx {|{b}\rangle}{|{a}\rangle}.$$

Le circuit est terminé. En l'absence de bruit, l'échantillonnage des registres quantiques renverra les chaînes de bits $b, a$ avec une probabilité de $1$.

Le produit scalaire booléen est un exemple de ce qu'on appelle les fonctions courbes (bent functions). Nous ne définirons pas les fonctions courbes ici, mais nous notons simplement qu'elles « résistent au maximum aux attaques qui cherchent à exploiter une dépendance des sorties par rapport à un sous-espace linéaire des entrées. » Cette citation provient de l'article Quantum algorithms for highly non-linear Boolean functions, qui propose des algorithmes de décalage caché efficaces pour plusieurs classes de fonctions courbes. L'algorithme de ce tutoriel apparaît dans la Section 3.1 de cet article.

Dans le cas plus général, le circuit pour trouver un décalage caché $s \in \mathbb{Z}^n$ est

$$H^{\otimes n} U_{\tilde{f}} H^{\otimes n} U_g H^{\otimes n} {|{0}\rangle}^{\otimes n} = {|{s}\rangle}.$$

Dans le cas général, $f$ et $g$ sont des fonctions d'une seule variable. Notre exemple du produit scalaire a cette forme si nous posons $f(x, y) \to f(z)$, avec $z$ égal à la concaténation de $x$ et $y$, et $s$ égal à la concaténation de $a$ et $b$. Le cas général nécessite exactement deux oracles : un oracle pour $g$ et un pour $\tilde{f}$, où ce dernier est une fonction connue sous le nom de duale de la fonction courbe $f$. La fonction produit scalaire possède la propriété d'auto-dualité $\tilde{f}=f$.

Dans notre circuit pour le décalage caché sur le produit scalaire, nous avons omis la couche intermédiaire de portes de Hadamard qui apparaît dans le circuit du cas général. Alors que dans le cas général cette couche est nécessaire, nous avons gagné un peu de profondeur en l'omettant, au prix d'un peu de post-traitement car la sortie est ${|{b}\rangle}{|{a}\rangle}$ au lieu du ${|{a}\rangle}{|{b}\rangle}$ souhaité.

Prérequis

Avant de commencer ce tutoriel, assurez-vous d'avoir installé les éléments suivants :

• Qiskit SDK v2.1 ou ultérieur, avec le support de visualisation
• Qiskit Runtime v0.41 ou ultérieur (pip install qiskit-ibm-runtime)
• Addon Qiskit M3 v3.0 (pip install mthree)

Configuration

# Added by doQumentation — installs packages not in the Binder environment
%pip install -q mthree

from collections.abc import Iterator, Sequence
from random import Random
from qiskit.circuit import (
    CircuitInstruction,
    QuantumCircuit,
    QuantumRegister,
    Qubit,
)
from qiskit.circuit.library import CZGate, HGate, XGate
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
from qiskit_ibm_runtime import QiskitRuntimeService
import timeit
import matplotlib.pyplot as plt
from qiskit_ibm_runtime import SamplerV2 as Sampler
import mthree

Étape 1 : Transposer les entrées classiques en un problème quantique

Tout d'abord, nous écrivons les fonctions pour implémenter le problème du décalage caché sous forme de QuantumCircuit.

def apply_hadamards(qubits: Sequence[Qubit]) -> Iterator[CircuitInstruction]:
    """Apply a Hadamard gate to every qubit."""
    for q in qubits:
        yield CircuitInstruction(HGate(), [q], [])

def apply_shift(
    qubits: Sequence[Qubit], shift: int
) -> Iterator[CircuitInstruction]:
    """Apply X gates where the bits of the shift are equal to 1."""
    for i, q in zip(range(shift.bit_length()), qubits):
        if shift >> i & 1:
            yield CircuitInstruction(XGate(), [q], [])

def oracle_f(qubits: Sequence[Qubit]) -> Iterator[CircuitInstruction]:
    """Apply the f oracle."""
    for i in range(0, len(qubits) - 1, 2):
        yield CircuitInstruction(CZGate(), [qubits[i], qubits[i + 1]])

def oracle_g(
    qubits: Sequence[Qubit], shift: int
) -> Iterator[CircuitInstruction]:
    """Apply the g oracle."""
    yield from apply_shift(qubits, shift)
    yield from oracle_f(qubits)
    yield from apply_shift(qubits, shift)

def determine_hidden_shift(
    qubits: Sequence[Qubit], shift: int
) -> Iterator[CircuitInstruction]:
    """Determine the hidden shift."""
    yield from apply_hadamards(qubits)
    yield from oracle_g(qubits, shift)
    # We omit this layer in exchange for post processing
    # yield from apply_hadamards(qubits)
    yield from oracle_f(qubits)
    yield from apply_hadamards(qubits)

def run_hidden_shift_circuit(n_qubits, rng):
    hidden_shift = rng.getrandbits(n_qubits)

    qubits = QuantumRegister(n_qubits, name="q")
    circuit = QuantumCircuit.from_instructions(
        determine_hidden_shift(qubits, hidden_shift), qubits=qubits
    )
    circuit.measure_all()
    # Format the hidden shift as a string.
    hidden_shift_string = format(hidden_shift, f"0{n_qubits}b")
    return (circuit, hidden_shift, hidden_shift_string)

def display_circuit(circuit):
    return circuit.remove_final_measurements(inplace=False).draw(
        "mpl", idle_wires=False, scale=0.5, fold=-1
    )

Nous commencerons par un petit exemple :

n_qubits = 6
random_seed = 12345
rng = Random(random_seed)
circuit, hidden_shift, hidden_shift_string = run_hidden_shift_circuit(
    n_qubits, rng
)

print(f"Hidden shift string {hidden_shift_string}")

display_circuit(circuit)

Hidden shift string 011010

Étape 2 : Optimiser les circuits pour l'exécution sur du matériel quantique

job_tags = [
    f"shift {hidden_shift_string}",
    f"n_qubits {n_qubits}",
    f"seed = {random_seed}",
]
job_tags

['shift 011010', 'n_qubits 6', 'seed = 12345']

# Uncomment this to run the circuits on a quantum computer on IBMCloud.
service = QiskitRuntimeService()
backend = service.least_busy(
    operational=True, simulator=False, min_num_qubits=100
)

# from qiskit_ibm_runtime.fake_provider import FakeMelbourneV2
# backend = FakeMelbourneV2()
# backend.refresh(service)

print(f"Using backend {backend.name}")

def get_isa_circuit(circuit, backend):
    pass_manager = generate_preset_pass_manager(
        optimization_level=3, backend=backend, seed_transpiler=1234
    )
    isa_circuit = pass_manager.run(circuit)
    return isa_circuit

isa_circuit = get_isa_circuit(circuit, backend)
display_circuit(isa_circuit)

Using backend ibm_kingston

Étape 3 : Exécuter les circuits à l'aide des primitives Qiskit

# submit job for solving the hidden shift problem using the Sampler primitive
NUM_SHOTS = 50_000

def run_sampler(backend, isa_circuit, num_shots):
    sampler = Sampler(mode=backend)
    sampler.options.environment.job_tags
    pubs = [(isa_circuit, None, NUM_SHOTS)]
    job = sampler.run(pubs)
    return job

def setup_mthree_mitigation(isa_circuit, backend):
    # retrieve the final qubit mapping so mthree knows which qubits to calibrate
    qubit_mapping = mthree.utils.final_measurement_mapping(isa_circuit)

    # submit jobs for readout error calibration
    mit = mthree.M3Mitigation(backend)
    mit.cals_from_system(qubit_mapping, rep_delay=None)

    return mit, qubit_mapping

job = run_sampler(backend, isa_circuit, NUM_SHOTS)
mit, qubit_mapping = setup_mthree_mitigation(isa_circuit, backend)

Étape 4 : Post-traiter et renvoyer les résultats au format classique

Dans la discussion théorique ci-dessus, nous avons déterminé que pour l'entrée $ab$, nous attendons la sortie $ba$. Une complication supplémentaire est que, afin d'avoir un circuit plus simple (avant transpilation), nous avons inséré les portes CZ requises entre des paires de qubits voisins. Cela revient à entrelacer les chaînes de bits $a$ et $b$ sous la forme $a1 b1 a2 b2 \ldots$. La chaîne de sortie $ba$ sera entrelacée de manière similaire : $b1 a1 b2 a2 \ldots$. La fonction unscramble ci-dessous transforme la chaîne de sortie de $b1 a1 b2 a2 \ldots$ en $a1 b1 a2 b2 \ldots$ afin que les chaînes d'entrée et de sortie puissent être comparées directement.

# retrieve bitstring counts
def get_bitstring_counts(job):
    result = job.result()
    pub_result = result[0]
    counts = pub_result.data.meas.get_counts()
    return counts, pub_result

counts, pub_result = get_bitstring_counts(job)

La distance de Hamming entre deux chaînes de bits est le nombre d'indices auxquels les bits diffèrent.

def hamming_distance(s1, s2):
    weight = 0
    for c1, c2 in zip(s1, s2):
        (c1, c2) = (int(c1), int(c2))
        if (c1 == 1 and c2 == 1) or (c1 == 0 and c2 == 0):
            weight += 1

    return weight

# Replace string of form a1b1a2b2... with b1a1b2a1...
# That is, reverse order of successive pairs of bits.
def unscramble(bitstring):
    ps = [bitstring[i : i + 2][::-1] for i in range(0, len(bitstring), 2)]
    return "".join(ps)

def find_hidden_shift_bitstring(counts, hidden_shift_string):
    # convert counts to probabilities
    probs = {
        unscramble(bitstring): count / NUM_SHOTS
        for bitstring, count in counts.items()
    }

    # Retrieve the most probable bitstring.
    most_probable = max(probs, key=lambda x: probs[x])

    print(f"Expected hidden shift string: {hidden_shift_string}")
    if most_probable == hidden_shift_string:
        print("Most probable bitstring matches hidden shift 😊.")
    else:
        print("Most probable bitstring didn't match hidden shift ☹️.")
    print("Top 10 bitstrings and their probabilities:")
    display(
        {
            k: (v, hamming_distance(hidden_shift_string, k))
            for k, v in sorted(
                probs.items(), key=lambda x: x[1], reverse=True
            )[:10]
        }
    )

    return probs, most_probable

probs, most_probable = find_hidden_shift_bitstring(
    counts, hidden_shift_string
)

Expected hidden shift string: 011010
Most probable bitstring matches hidden shift 😊.
Top 10 bitstrings and their probabilities:

{'011010': (0.9743, 6),
 '001010': (0.00812, 5),
 '010010': (0.0063, 5),
 '011000': (0.00554, 5),
 '011011': (0.00492, 5),
 '011110': (0.00044, 5),
 '001000': (0.00012, 4),
 '010000': (8e-05, 4),
 '001011': (6e-05, 4),
 '000010': (6e-05, 4)}

Enregistrons la probabilité de la chaîne de bits la plus probable avant d'appliquer l'atténuation des erreurs de lecture avec M3.

max_probability_before_M3 = probs[most_probable]
max_probability_before_M3

0.9743

Maintenant, nous appliquons la correction de lecture apprise par M3 aux comptages. La fonction apply_corrections renvoie une distribution de quasi-probabilité. Il s'agit d'une liste d'objets float dont la somme est $1$. Cependant, certaines valeurs peuvent être négatives.

def perform_mitigation(mit, counts, qubit_mapping):
    # mitigate readout error
    quasis = mit.apply_correction(counts, qubit_mapping)

    # print results
    most_probable_after_m3 = unscramble(max(quasis, key=lambda x: quasis[x]))

    is_hidden_shift_identified = most_probable_after_m3 == hidden_shift_string
    if is_hidden_shift_identified:
        print("Most probable bitstring matches hidden shift 😊.")
    else:
        print("Most probable bitstring didn't match hidden shift ☹️.")
    print("Top 10 bitstrings and their quasi-probabilities:")
    topten = {
        unscramble(k): f"{v:.2e}"
        for k, v in sorted(quasis.items(), key=lambda x: x[1], reverse=True)[
            :10
        ]
    }
    max_probability_after_M3 = float(topten[most_probable_after_m3])
    display(topten)

    return max_probability_after_M3, is_hidden_shift_identified

print(f"Expected hidden shift string: {hidden_shift_string}")
max_probability_after_M3, is_hidden_shift_identified = perform_mitigation(
    mit, counts, qubit_mapping
)

Expected hidden shift string: 011010
Most probable bitstring matches hidden shift 😊.
Top 10 bitstrings and their quasi-probabilities:

{'011010': '1.01e+00',
 '001010': '8.75e-04',
 '001000': '7.38e-05',
 '010000': '4.51e-05',
 '111000': '2.18e-05',
 '001011': '1.74e-05',
 '000010': '6.42e-06',
 '011001': '-7.18e-06',
 '011000': '-4.53e-04',
 '010010': '-1.28e-03'}

Comparer l'identification de la chaîne de bits du décalage caché avant et après l'application de la correction M3

def compare_before_and_after_M3(
    max_probability_before_M3,
    max_probability_after_M3,
    is_hidden_shift_identified,
):
    is_probability_improved = (
        max_probability_after_M3 > max_probability_before_M3
    )
    print(f"Most probable probability before M3: {max_probability_before_M3}")
    print(f"Most probable probability after M3: {max_probability_after_M3}")
    if is_hidden_shift_identified and is_probability_improved:
        print("Readout error mitigation effective! 😊")
    else:
        print("Readout error mitigation not effective. ☹️")

compare_before_and_after_M3(
    max_probability_before_M3,
    max_probability_after_M3,
    is_hidden_shift_identified,
)

Most probable probability before M3: 0.9743
Most probable probability after M3: 1.01
Readout error mitigation effective! 😊

Tracer l'évolution du temps CPU requis par M3 en fonction du nombre de shots

# Collect samples for numbers of shots varying from 5000 to 25000.
shots_range = range(5000, NUM_SHOTS + 1, 2500)
times = []
for shots in shots_range:
    print(f"Applying M3 correction to {shots} shots...")
    t0 = timeit.default_timer()
    _ = mit.apply_correction(
        pub_result.data.meas.slice_shots(range(shots)).get_counts(),
        qubit_mapping,
    )
    t1 = timeit.default_timer()
    print(f"\tDone in {t1 - t0} seconds.")
    times.append(t1 - t0)

fig, ax = plt.subplots()
ax.plot(shots_range, times, "o--")
ax.set_xlabel("Shots")
ax.set_ylabel("Time (s)")
ax.set_title("Time to apply M3 correction")

Applying M3 correction to 5000 shots...
	Done in 0.003321983851492405 seconds.
Applying M3 correction to 7500 shots...
	Done in 0.004425413906574249 seconds.
Applying M3 correction to 10000 shots...
	Done in 0.006366567220538855 seconds.
Applying M3 correction to 12500 shots...
	Done in 0.0071477219462394714 seconds.
Applying M3 correction to 15000 shots...
	Done in 0.00860048783943057 seconds.
Applying M3 correction to 17500 shots...
	Done in 0.010026784148067236 seconds.
Applying M3 correction to 20000 shots...
	Done in 0.011459112167358398 seconds.
Applying M3 correction to 22500 shots...
	Done in 0.012727141845971346 seconds.
Applying M3 correction to 25000 shots...
	Done in 0.01406092382967472 seconds.
Applying M3 correction to 27500 shots...
	Done in 0.01546052098274231 seconds.
Applying M3 correction to 30000 shots...
	Done in 0.016769016161561012 seconds.
Applying M3 correction to 32500 shots...
	Done in 0.019537431187927723 seconds.
Applying M3 correction to 35000 shots...
	Done in 0.019739801064133644 seconds.
Applying M3 correction to 37500 shots...
	Done in 0.021093040239065886 seconds.
Applying M3 correction to 40000 shots...
	Done in 0.022840639110654593 seconds.
Applying M3 correction to 42500 shots...
	Done in 0.023974396288394928 seconds.
Applying M3 correction to 45000 shots...
	Done in 0.026412792038172483 seconds.
Applying M3 correction to 47500 shots...
	Done in 0.026364430785179138 seconds.
Applying M3 correction to 50000 shots...
	Done in 0.02820305060595274 seconds.

Text(0.5, 1.0, 'Time to apply M3 correction')

Interprétation du graphique

Le graphique ci-dessus montre que le temps requis pour appliquer la correction M3 croît linéairement avec le nombre de shots.

Passage à l'échelle

n_qubits = 80
rng = Random(12345)
circuit, hidden_shift, hidden_shift_string = run_hidden_shift_circuit(
    n_qubits, rng
)

print(f"Hidden shift string {hidden_shift_string}")

Hidden shift string 00000010100110101011101110010001010000110011101001101010101001111001100110000111

isa_circuit = get_isa_circuit(circuit, backend)

job = run_sampler(backend, isa_circuit, NUM_SHOTS)
mit, qubit_mapping = setup_mthree_mitigation(isa_circuit, backend)

counts, pub_result = get_bitstring_counts(job)

probs, most_probable = find_hidden_shift_bitstring(
    counts, hidden_shift_string
)

Expected hidden shift string: 00000010100110101011101110010001010000110011101001101010101001111001100110000111
Most probable bitstring matches hidden shift 😊.
Top 10 bitstrings and their probabilities:

{'00000010100110101011101110010001010000110011101001101010101001111001100110000111': (0.50402,
  80),
 '00000010100110101011101110010001010000110011100001101010101001111001100110000111': (0.0396,
  79),
 '00000010100110101011101110010001010000110011101001101010101001111001100100000111': (0.0323,
  79),
 '00000010100110101011101110010001010000110011101001101010101001101001100110000111': (0.01936,
  79),
 '00000010100110101011101110010011010000110011101001101010101001111001100110000111': (0.01432,
  79),
 '00000010100110101011101110010001010000110011101001101010101001011001100110000111': (0.0101,
  79),
 '00000010100110101011101110010001010000110011101001101010101001110001100110000111': (0.00924,
  79),
 '00000010100110101011101110010001010000010011101001101010101001111001100110000111': (0.00908,
  79),
 '00000010100110101011101110010001010000110011101001101010101001111001100110000111': (0.00888,
  79),
 '00000010100110101011101110010001010000110011101001100010101001111001100110000111': (0.0082,
  79)}

Nous voyons que la bonne chaîne de bits du décalage caché est trouvée. De plus, les neuf chaînes de bits suivantes les plus probables ne diffèrent que d'une seule position. Enregistrons la probabilité la plus élevée :

max_probability_before_M3 = probs[most_probable]
max_probability_before_M3

0.50402

print(f"Expected hidden shift string: {hidden_shift_string}")
max_probability_after_M3, is_hidden_shift_identified = perform_mitigation(
    mit, counts, qubit_mapping
)

Expected hidden shift string: 00000010100110101011101110010001010000110011101001101010101001111001100110000111
Most probable bitstring matches hidden shift 😊.
Top 10 bitstrings and their quasi-probabilities:

{'00000010100110101011101110010001010000110011101001101010101001111001100110000111': '9.85e-01',
 '00000010100110101011101110010001010000110011100001101010101001111001100110000111': '6.84e-03',
 '00000010100110101011100110010001010000110011101001101010101001111001100110000111': '3.87e-03',
 '00000010100110101011101110010011010000110011101001101010101001111001100110000111': '3.42e-03',
 '00000010100110101011101110010001010000110011101001101010101001111001100100000111': '3.30e-03',
 '00000010100110101011101110010001010000110011101001101010101001110001100110000111': '3.28e-03',
 '00000010100010101011101110010001010000110011101001101010101001111001100110000111': '2.62e-03',
 '00000010100110101011101110010001010000110011101001101010101001101001100110000111': '2.43e-03',
 '00000010100110101011101110010000010000110011101001101010101001111001100110000111': '1.73e-03',
 '00000010100110101011101110010001010000110011101001101010101001111001000110000111': '1.63e-03'}

compare_before_and_after_M3(
    max_probability_before_M3,
    max_probability_after_M3,
    is_hidden_shift_identified,
)

Most probable probability before M3: 0.54348
Most probable probability after M3: 0.99
Readout error mitigation effective! 😊

Les résultats montrent que l'erreur de lecture était la source d'erreur dominante et que l'atténuation M3 a été efficace.