Fonctions de coût

Dans cette leçon, nous apprendrons à évaluer une fonction de coût :

• D'abord, nous découvrirons les primitives de Qiskit Runtime
• Nous définirons une fonction de coût $C(\vec\theta)$. C'est une fonction spécifique au problème qui définit l'objectif que l'optimiseur doit minimiser (ou maximiser)
• Nous définirons une stratégie de mesure avec les primitives de Qiskit Runtime pour équilibrer vitesse et précision

 

Primitives

Tous les systèmes physiques, qu'ils soient classiques ou quantiques, peuvent exister dans différents états. Par exemple, une voiture sur une route peut avoir une certaine masse, position, vitesse ou accélération qui caractérisent son état. De la même manière, les systèmes quantiques peuvent aussi avoir différentes configurations ou états, mais ils diffèrent des systèmes classiques dans la façon dont nous traitons les mesures et l'évolution de l'état. Cela donne lieu à des propriétés uniques comme la superposition et l'intrication, exclusives à la mécanique quantique. Tout comme nous pouvons décrire l'état d'une voiture en utilisant des propriétés physiques comme la vitesse ou l'accélération, nous pouvons aussi décrire l'état d'un système quantique en utilisant des observables, qui sont des objets mathématiques.

En mécanique quantique, les états sont représentés par des vecteurs colonnes complexes normalisés, ou kets ($|\psi\rangle$), et les observables sont des opérateurs linéaires hermitiens ($\hat{H}=\hat{H}^{\dagger}$) qui agissent sur les kets. Un vecteur propre ($|\lambda\rangle$) d'un observable est connu comme un état propre. Mesurer un observable pour l'un de ses états propres ($|\lambda\rangle$) nous donnera la valeur propre correspondante ($\lambda$) comme résultat.

Si tu te demandes comment mesurer un système quantique et ce que tu peux mesurer, Qiskit offre deux primitives qui peuvent t'aider :

• Sampler : étant donné un état quantique $|\psi\rangle$, cette primitive obtient la probabilité de chaque état possible de la base computationnelle.
• Estimator : étant donné un observable quantique $\hat{H}$ et un état $|\psi\rangle$, cette primitive calcule la valeur d'espérance de $\hat{H}$.

La primitive Sampler

La primitive Sampler calcule la probabilité d'obtenir chaque état possible $|k\rangle$ de la base computationnelle, étant donné un circuit quantique qui prépare l'état $|\psi\rangle$. Elle calcule

$$p_k = |\langle k | \psi \rangle|^2 \quad \forall k \in \mathbb{Z}_2^n \equiv \{0,1,\cdots,2^n-1\},$$

où $n$ est le nombre de qubits, et $k$ la représentation entière de toute chaîne binaire de sortie possible $\{0,1\}^n$ (c'est-à-dire des entiers en base $2$).

Le Sampler de Qiskit Runtime exécute le circuit plusieurs fois sur un dispositif quantique, effectue des mesures à chaque exécution et reconstruit la distribution de probabilité à partir des chaînes de bits obtenues. Plus il y a d'exécutions (ou shots), plus les résultats seront précis, mais cela nécessite plus de temps et de ressources quantiques.

Cependant, étant donné que le nombre de sorties possibles croît exponentiellement avec le nombre de qubits $n$ (c'est-à-dire $2^n$), le nombre de shots devra également croître exponentiellement pour capturer une distribution de probabilité dense. Par conséquent, Sampler n'est efficace que pour des distributions de probabilité éparses ; où l'état cible $|\psi\rangle$ doit pouvoir s'exprimer comme une combinaison linéaire des états de la base computationnelle, avec le nombre de termes croissant au plus polynomialement avec le nombre de qubits :

$$|\psi\rangle = \sum^{\text{Poly}(n)}_k w_k |k\rangle.$$

Le Sampler peut également être configuré pour récupérer les probabilités d'une sous-section du circuit, ce qui représente un sous-ensemble du total des états possibles.

La primitive Estimator

La primitive Estimator calcule la valeur d'espérance d'un observable $\hat{H}$ pour un état quantique $|\psi\rangle$ ; où les probabilités de l'observable peuvent s'exprimer comme $p_\lambda = |\langle\lambda|\psi\rangle|^2$, $|\lambda\rangle$ étant les états propres de l'observable $\hat{H}$. La valeur d'espérance est alors définie comme la moyenne de tous les résultats possibles $\lambda$ (c'est-à-dire les valeurs propres de l'observable) d'une mesure de l'état $|\psi\rangle$, pondérée par les probabilités correspondantes :

$$\langle\hat{H}\rangle_\psi := \sum_\lambda p_\lambda \lambda = \langle \psi | \hat{H} | \psi \rangle$$

Cependant, calculer la valeur d'espérance d'un observable n'est pas toujours possible, car nous ne connaissons souvent pas sa base propre. L'Estimator de Qiskit Runtime utilise un processus algébrique complexe pour estimer la valeur d'espérance sur un dispositif quantique réel, en décomposant l'observable en une combinaison d'autres observables dont nous connaissons la base propre.

En termes plus simples, Estimator décompose tout observable qu'il ne sait pas mesurer en observables plus simples et mesurables appelés opérateurs de Pauli. Tout opérateur peut s'exprimer comme une combinaison de $4^n$ opérateurs de Pauli.

$$\hat{P}_k := \sigma_{k_{n-1}}\otimes \cdots \otimes \sigma_{k_0} \quad \forall k \in \mathbb{Z}_4^n \equiv \{0,1,\cdots,4^n-1\}, \\$$

tel que

$$\hat{H} = \sum^{4^n-1}_{k=0} w_k \hat{P}_k$$

où $n$ est le nombre de qubits, $k \equiv k_{n-1} \cdots k_0$ pour $k_l \in \mathbb{Z}_4 \equiv \{0, 1, 2, 3\}$ (c'est-à-dire des entiers en base $4$), et $(\sigma_0, \sigma_1, \sigma_2, \sigma_3) := (I, X, Y, Z)$.

Après avoir effectué cette décomposition, Estimator dérive un nouveau circuit $V_k|\psi\rangle$ pour chaque observable $\hat{P}_k$ (à partir du circuit original), afin de diagonaliser effectivement l'observable de Pauli dans la base computationnelle et de le mesurer. Nous pouvons facilement mesurer les observables de Pauli car nous connaissons $V_k$ à l'avance, ce qui n'est pas le cas en général pour d'autres observables.

Pour chaque $\hat{P}_{k}$, l'Estimator exécute le circuit correspondant sur un dispositif quantique plusieurs fois, mesure l'état de sortie dans la base computationnelle et calcule la probabilité $p_{kj}$ d'obtenir chaque sortie possible $j$. Il cherche ensuite la valeur propre $\lambda_{kj}$ de $P_k$ correspondant à chaque sortie $j$, multiplie par $w_k$ et additionne tous les résultats pour obtenir la valeur d'espérance de l'observable $\hat{H}$ pour l'état donné $|\psi\rangle$.

$$\langle\hat{H}\rangle_\psi = \sum_{k=0}^{4^n-1} w_k \sum_{j=0}^{2^n-1}p_{kj} \lambda_{kj},$$

Étant donné que calculer la valeur d'espérance de $4^n$ Paulis est peu pratique (c'est-à-dire que cela croît exponentiellement), Estimator ne peut être efficace que lorsqu'une grande quantité de $w_k$ sont nuls (c'est-à-dire une décomposition de Pauli éparse plutôt que dense). Formellement, nous disons que, pour que ce calcul soit efficacement soluble, le nombre de termes non nuls doit croître au plus polynomialement avec le nombre de qubits $n$ : $\hat{H} = \sum^{\text{Poly}(n)}_k w_k \hat{P}_k.$

Le lecteur peut noter l'hypothèse implicite que l'échantillonnage de probabilités doit également être efficace, comme expliqué pour Sampler, ce qui signifie

$$\langle\hat{H}\rangle_\psi = \sum_{k}^{\text{Poly}(n)} w_k \sum_{j}^{\text{Poly}(n)}p_{kj} \lambda_{kj}.$$

Exemple guidé pour calculer des valeurs d'espérance

Supposons l'état d'un qubit $|+\rangle := H|0\rangle = \frac{1}{\sqrt{2}}(|0\rangle + |1\rangle)$, et l'observable

$$\begin{aligned} \hat{H} & = \begin{pmatrix} -1 & 2 \\ 2 & 1 \\ \end{pmatrix}\\[1mm] & = 2X - Z \end{aligned}$$

avec la valeur d'espérance théorique suivante $\langle\hat{H}\rangle_+ = \langle+|\hat{H}|+\rangle = 2.$

Comme nous ne savons pas comment mesurer cet observable, nous ne pouvons pas calculer sa valeur d'espérance directement et nous devons le réexprimer comme $\langle\hat{H}\rangle_+ = 2\langle X \rangle_+ - \langle Z \rangle_+$. On peut montrer que cela donne le même résultat en notant que $\langle+|X|+\rangle = 1$ et $\langle+|Z|+\rangle = 0$.

Voyons comment calculer $\langle X \rangle_+$ et $\langle Z \rangle_+$ directement. Comme $X$ et $Z$ ne commutent pas (c'est-à-dire qu'ils ne partagent pas la même base propre), ils ne peuvent pas être mesurés simultanément, nous avons donc besoin des circuits auxiliaires :

# Added by doQumentation — required packages for this notebook
!pip install -q matplotlib numpy qiskit qiskit-aer qiskit-ibm-runtime rustworkx

from qiskit import QuantumCircuit
from qiskit.quantum_info import SparsePauliOp

# The following code will work for any other initial single-qubit state and observable
original_circuit = QuantumCircuit(1)
original_circuit.h(0)

H = SparsePauliOp(["X", "Z"], [2, -1])

aux_circuits = []
for pauli in H.paulis:
    aux_circ = original_circuit.copy()
    aux_circ.barrier()
    if str(pauli) == "X":
        aux_circ.h(0)
    elif str(pauli) == "Y":
        aux_circ.sdg(0)
        aux_circ.h(0)
    else:
        aux_circ.id(0)
    aux_circ.measure_all()
    aux_circuits.append(aux_circ)

original_circuit.draw("mpl")

# Auxiliary circuit for X
aux_circuits[0].draw("mpl")

# Auxiliary circuit for Z
aux_circuits[1].draw("mpl")

Nous pouvons maintenant effectuer le calcul manuellement en utilisant Sampler et vérifier les résultats avec Estimator :

from qiskit.primitives import StatevectorSampler, StatevectorEstimator
from qiskit.result import QuasiDistribution
import numpy as np

## SAMPLER
shots = 10000
sampler = StatevectorSampler()
job = sampler.run(aux_circuits, shots=shots)

# Run the sampler job and step through results
expvals = []
for index, pauli in enumerate(H.paulis):
    data_pub = job.result()[index].data
    bitstrings = data_pub.meas.get_bitstrings()
    counts = data_pub.meas.get_counts()
    quasi_dist = QuasiDistribution(
        {outcome: freq / shots for outcome, freq in counts.items()}
    )

    # Use the probabilities and known eigenvalues of Pauli operators to estimate the expectation value.
    val = 0

    if str(pauli) == "X":
        val += -1 * quasi_dist.get(1, 0)
        val += 1 * quasi_dist.get(0, 0)

    if str(pauli) == "Y":
        val += -1 * quasi_dist.get(1, 0)
        val += 1 * quasi_dist.get(0, 0)

    if str(pauli) == "Z":
        val += 1 * quasi_dist.get(0, 0)
        val += -1 * quasi_dist.get(1, 0)

    expvals.append(val)

# Print expectation values

print("Sampler results:")
for pauli, expval in zip(H.paulis, expvals):
    print(f"  >> Expected value of {str(pauli)}: {expval:.5f}")

total_expval = np.sum(H.coeffs * expvals).real
print(f"  >> Total expected value: {total_expval:.5f}")

# Use estimator for comparison
observables = [
    *H.paulis,
    H,
]  # Note: run for individual Paulis as well as full observable H

estimator = StatevectorEstimator()
job = estimator.run([(original_circuit, observables)])
estimator_expvals = job.result()[0].data.evs

# Print results
print("Estimator results:")
for obs, expval in zip(observables, estimator_expvals):
    if obs is not H:
        print(f"  >> Expected value of {str(obs)}: {expval:.5f}")
    else:
        print(f"  >> Total expected value: {expval:.5f}")

Sampler results:
  >> Expected value of X: 1.00000
  >> Expected value of Z: 0.00420
  >> Total expected value: 1.99580
Estimator results:
  >> Expected value of X: 1.00000
  >> Expected value of Z: 0.00000
  >> Total expected value: 2.00000

Rigueur mathématique (optionnel)

En exprimant $|\psi\rangle$ par rapport à la base des états propres de $\hat{H}$, $|\psi\rangle = \sum_\lambda a_\lambda |\lambda\rangle$, on déduit :

$$\begin{aligned} \langle \psi | \hat{H} | \psi \rangle & = \bigg(\sum_{\lambda'}a^*_{\lambda'} \langle \lambda'|\bigg) \hat{H} \bigg(\sum_{\lambda} a_\lambda | \lambda\rangle\bigg)\\[1mm] & = \sum_{\lambda}\sum_{\lambda'} a^*_{\lambda'}a_{\lambda} \langle \lambda'|\hat{H}| \lambda\rangle\\[1mm] & = \sum_{\lambda}\sum_{\lambda'} a^*_{\lambda'}a_{\lambda} \lambda \langle \lambda'| \lambda\rangle\\[1mm] & = \sum_{\lambda}\sum_{\lambda'} a^*_{\lambda'}a_{\lambda} \lambda \cdot \delta_{\lambda, \lambda'}\\[1mm] & = \sum_\lambda |a_\lambda|^2 \lambda\\[1mm] & = \sum_\lambda p_\lambda \lambda\\[1mm] \end{aligned}$$

Comme nous ne connaissons pas les valeurs propres ni les états propres de l'observable cible $\hat{H}$, nous devons d'abord considérer sa diagonalisation. Étant donné que $\hat{H}$ est hermitien, il existe une transformation unitaire $V$ tel que $\hat{H}=V^\dagger \Lambda V,$ où $\Lambda$ est la matrice diagonale des valeurs propres, de sorte que $\langle j | \Lambda | k \rangle = 0$ si $j\neq k$, et $\langle j | \Lambda | j \rangle = \lambda_j$.

Cela implique que la valeur d'espérance peut être réécrite comme :

$$\begin{aligned} \langle\psi|\hat{H}|\psi\rangle & = \langle\psi|V^\dagger \Lambda V|\psi\rangle\\[1mm] & = \langle\psi|V^\dagger \bigg(\sum_{j=0}^{2^n-1} |j\rangle \langle j|\bigg) \Lambda \bigg(\sum_{k=0}^{2^n-1} |k\rangle \langle k|\bigg) V|\psi\rangle\\[1mm] & = \sum_{j=0}^{2^n-1} \sum_{k=0}^{2^n-1}\langle\psi|V^\dagger |j\rangle \langle j| \Lambda |k\rangle \langle k| V|\psi\rangle\\[1mm] & = \sum_{j=0}^{2^n-1}\langle\psi|V^\dagger |j\rangle \langle j| \Lambda |j\rangle \langle j| V|\psi\rangle\\[1mm] & = \sum_{j=0}^{2^n-1}|\langle j| V|\psi\rangle|^2 \lambda_j\\[1mm] \end{aligned}$$

Étant donné que si un système est dans l'état $|\phi\rangle = V |\psi\rangle$ la probabilité de mesurer $| j\rangle$ est $p_j = |\langle j|\phi \rangle|^2$, la valeur d'espérance ci-dessus peut s'exprimer comme :

$$\langle\psi|\hat{H}|\psi\rangle = \sum_{j=0}^{2^n-1} p_j \lambda_j.$$

Il est très important de noter que les probabilités sont prises à partir de l'état $V |\psi\rangle$ plutôt que $|\psi\rangle$. C'est pourquoi la matrice $V$ est absolument nécessaire. Tu te demandes peut-être comment obtenir la matrice $V$ et les valeurs propres $\Lambda$. Si tu avais déjà les valeurs propres, il n'y aurait pas besoin d'utiliser un ordinateur quantique, puisque l'objectif des algorithmes variationnels est précisément de trouver ces valeurs propres de $\hat{H}$.

Heureusement, il existe une solution : toute matrice $2^n \times 2^n$ peut s'écrire comme une combinaison linéaire de $4^n$ produits tensoriels de $n$ matrices de Pauli et identités, toutes hermitiennes et unitaires avec $V$ et $\Lambda$ connus. C'est ce que l'Estimator de Runtime fait en interne en décomposant tout objet Operator en un SparsePauliOp.

Voici les opérateurs qui peuvent être utilisés :

$$\begin{array}{c|c|c|c} \text{Operator} & \sigma & V & \Lambda \\[1mm] \hline I & \sigma_0 = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} & V_0 = I & \Lambda_0 = I = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} \\[4mm] X & \sigma_1 = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} & V_1 = H =\frac{1}{\sqrt{2}} \begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix} & \Lambda_1 = \sigma_3 = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} \\[4mm] Y & \sigma_2 = \begin{pmatrix} 0 & -i \\ i & 0 \end{pmatrix} & V_2 = HS^\dagger =\frac{1}{\sqrt{2}} \begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix}\cdot \begin{pmatrix} 1 & 0 \\ 0 & -i \end{pmatrix} = \frac{1}{\sqrt{2}} \begin{pmatrix} 1 & -i \\ 1 & i \end{pmatrix}\quad & \Lambda_2 = \sigma_3 = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} \\[4mm] Z & \sigma_3 = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} & V_3 = I & \Lambda_3 = \sigma_3 = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} \end{array}$$

Réécrivons donc $\hat{H}$ par rapport aux Paulis et identités :

$$\hat{H} = \sum_{k_{n-1}=0}^3... \sum_{k_0=0}^3 w_{k_{n-1}...k_0} \sigma_{k_{n-1}}\otimes ... \otimes \sigma_{k_0} = \sum_{k=0}^{4^n-1} w_k \hat{P}_k,$$

où $k = \sum_{l=0}^{n-1} 4^l k_l \equiv k_{n-1}...k_0$ pour $k_{n-1},...,k_0\in \{0,1,2,3\}$ (c'est-à-dire en base $4$), et $\hat{P}_{k} := \sigma_{k_{n-1}}\otimes ... \otimes \sigma_{k_0}$ :

$$\begin{aligned} \langle\psi|\hat{H}|\psi\rangle & = \sum_{k=0}^{4^n-1} w_k \sum_{j=0}^{2^n-1}|\langle j| V_k|\psi\rangle|^2 \langle j| \Lambda_k |j\rangle \\[1mm] & = \sum_{k=0}^{4^n-1} w_k \sum_{j=0}^{2^n-1}p_{kj} \lambda_{kj}, \\[1mm] \end{aligned}$$

où $V_k := V_{k_{n-1}}\otimes ... \otimes V_{k_0}$ et $\Lambda_k := \Lambda_{k_{n-1}}\otimes ... \otimes \Lambda_{k_0}$, tel que : $\hat{P_k}=V_k^\dagger \Lambda_k V_k.$

Fonctions de coût

En général, les fonctions de coût sont utilisées pour décrire l'objectif d'un problème et à quel point un état d'essai fonctionne bien par rapport à cet objectif. Cette définition peut s'appliquer à différents contextes, comme la chimie, l'apprentissage automatique, la finance, l'optimisation, entre autres.

Considérons un exemple simple : trouver l'état fondamental d'un système. Notre objectif est de minimiser la valeur d'espérance de l'observable qui représente l'énergie (Hamiltonien $\hat{\mathcal{H}}$) :

$$\min_{\vec\theta} \langle\psi(\vec\theta)|\hat{\mathcal{H}}|\psi(\vec\theta)\rangle$$

Nous pouvons utiliser l'Estimator pour évaluer la valeur d'espérance et la passer à un optimiseur pour qu'il la minimise. Si l'optimisation réussit, elle retournera un ensemble de valeurs de paramètres optimaux $\vec\theta^*$, à partir desquels nous pourrons construire l'état solution proposé $|\psi(\vec\theta^*)\rangle$ et calculer la valeur d'espérance observée comme $C(\vec\theta^*)$.

Note que nous ne pourrons minimiser la fonction de coût que pour l'ensemble limité d'états que nous considérons. Cela nous amène à deux possibilités distinctes :

• Notre ansatz ne définit pas l'état solution dans l'espace de recherche : dans ce cas, l'optimiseur ne trouvera jamais la solution et nous devrons expérimenter avec d'autres ansatzes capables de représenter notre espace de recherche avec plus de précision.
• L'optimiseur est incapable de trouver une solution valide : l'optimisation peut être définie de manière globale ou locale. Nous explorerons ce que cela signifie dans la section suivante.

En définitive, nous exécuterons un cycle d'optimisation classique tout en nous appuyant sur l'évaluation de la fonction de coût par un ordinateur quantique. De ce point de vue, on pourrait concevoir l'optimisation comme une tâche purement classique dans laquelle nous appelons un oracle quantique de boîte noire chaque fois que l'optimiseur a besoin d'évaluer la fonction de coût.

def cost_func_vqe(params, circuit, hamiltonian, estimator):
    """Return estimate of energy from estimator

    Parameters:
        params (ndarray): Array of ansatz parameters
        ansatz (QuantumCircuit): Parameterized ansatz circuit
        hamiltonian (SparsePauliOp): Operator representation of Hamiltonian
        estimator (Estimator): Estimator primitive instance

    Returns:
        float: Energy estimate
    """
    pub = (circuit, hamiltonian, params)
    cost = estimator.run([pub]).result()[0].data.evs
    return cost

from qiskit.circuit.library import TwoLocal

observable = SparsePauliOp.from_list([("XX", 1), ("YY", -3)])

reference_circuit = QuantumCircuit(2)
reference_circuit.x(0)

variational_form = TwoLocal(
    2,
    rotation_blocks=["rz", "ry"],
    entanglement_blocks="cx",
    entanglement="linear",
    reps=1,
)
ansatz = reference_circuit.compose(variational_form)

theta_list = (2 * np.pi * np.random.rand(1, 8)).tolist()
ansatz.decompose().draw("mpl")

D'abord, nous exécuterons ceci avec un simulateur : le StatevectorEstimator. C'est ce qui est recommandé pour le débogage, mais ensuite nous répéterons le calcul sur du matériel quantique réel. De plus en plus, les problèmes d'intérêt ne sont plus simulables classiquement sans des installations de supercalcul de pointe.

estimator = StatevectorEstimator()
cost = cost_func_vqe(theta_list, ansatz, observable, estimator)
print(cost)

[-0.58744589]

Nous allons maintenant exécuter sur un ordinateur quantique réel. Note les changements de syntaxe. Les étapes liées au pass_manager seront analysées plus en détail dans l'exemple suivant. Une étape particulièrement importante dans les algorithmes variationnels est l'utilisation d'une session de Qiskit Runtime. Ouvrir une session te permet d'exécuter plusieurs itérations d'un algorithme variationnel sans avoir à attendre dans une nouvelle file d'attente à chaque mise à jour des paramètres. C'est important si les temps d'attente sont longs et/ou si de nombreuses itérations sont nécessaires. Seuls les partenaires du réseau IBM Quantum® peuvent utiliser les sessions de Runtime. Si tu n'as pas accès aux sessions, tu peux réduire le nombre d'itérations que tu soumets à la fois et sauvegarder les paramètres les plus récents pour les utiliser dans des exécutions futures. Si tu soumets trop d'itérations ou que tu rencontres des temps d'attente trop longs, tu pourrais obtenir le code d'erreur 1217, qui fait référence à des délais longs entre les soumissions de jobs.

# Estimated usage: < 1 min. Benchmarked at 7 seconds on an Eagle processor
# Load necessary packages:

from qiskit_ibm_runtime import (
    QiskitRuntimeService,
    Session,
    EstimatorOptions,
    EstimatorV2 as Estimator,
)
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

# Select the least busy backend:

service = QiskitRuntimeService()
backend = service.least_busy(
    operational=True, min_num_qubits=ansatz.num_qubits, simulator=False
)
# Or get a specific backend:
# backend = service.backend("ibm_brisbane")

# Use a pass manager to transpile the circuit and observable for the specific backend being used:

pm = generate_preset_pass_manager(backend=backend, optimization_level=1)
isa_ansatz = pm.run(ansatz)
isa_observable = observable.apply_layout(layout=isa_ansatz.layout)

# Set estimator options
estimator_options = EstimatorOptions(resilience_level=1, default_shots=10_000)

# Open a Runtime session:

with Session(backend=backend) as session:
    estimator = Estimator(mode=session, options=estimator_options)
    cost = cost_func_vqe(theta_list, isa_ansatz, isa_observable, estimator)

session.close()
print(cost)

Note que les valeurs obtenues dans les deux calculs précédents sont très similaires. Les techniques pour améliorer les résultats seront analysées plus loin.

Exemple de correspondance avec des systèmes non physiques

Le problème de coupe maximale (Max-Cut) est un problème d'optimisation combinatoire qui consiste à diviser les sommets d'un graphe en deux ensembles disjoints de manière à maximiser le nombre d'arêtes entre les deux ensembles. Plus formellement, étant donné un graphe non dirigé $G=(V,E)$, où $V$ est l'ensemble des sommets et $E$ est l'ensemble des arêtes, le problème Max-Cut demande de diviser les sommets en deux sous-ensembles disjoints, $S$ et $T$, de sorte que le nombre d'arêtes ayant une extrémité dans $S$ et l'autre dans $T$ soit maximisé.

Nous pouvons appliquer Max-Cut pour résoudre différents problèmes, comme le regroupement (clustering), la conception de réseaux, les transitions de phase, entre autres. Commençons par créer un graphe du problème :

import rustworkx as rx
from rustworkx.visualization import mpl_draw

n = 4
G = rx.PyGraph()
G.add_nodes_from(range(n))
# The edge syntax is (start, end, weight)
edges = [(0, 1, 1.0), (0, 2, 1.0), (0, 3, 1.0), (1, 2, 1.0), (2, 3, 1.0)]
G.add_edges_from(edges)

mpl_draw(
    G, pos=rx.shell_layout(G), with_labels=True, edge_labels=str, node_color="#1192E8"
)

Ce problème peut s'exprimer comme un problème d'optimisation binaire. Pour chaque nœud $0 \leq i < n$, où $n$ est le nombre de nœuds du graphe (dans ce cas $n=4$), nous considérerons la variable binaire $x_i$. Cette variable aura la valeur $1$ si le nœud $i$ appartient à l'un des groupes que nous appellerons $1$, et $0$ s'il appartient à l'autre groupe, que nous appellerons $0$. Nous noterons également $w_{ij}$ (élément $(i,j)$ de la matrice d'adjacence $w$) le poids de l'arête allant du nœud $i$ au nœud $j$. Comme le graphe est non dirigé, $w_{ij}=w_{ji}$. Nous pouvons alors formuler notre problème comme la maximisation de la fonction de coût suivante :

$$\begin{aligned} C(\vec{x}) & =\sum_{i,j=0}^n w_{ij} x_i(1-x_j)\\[1mm] & = \sum_{i,j=0}^n w_{ij} x_i - \sum_{i,j=0}^n w_{ij} x_ix_j\\[1mm] & = \sum_{i,j=0}^n w_{ij} x_i - \sum_{i=0}^n \sum_{j=0}^i 2w_{ij} x_ix_j \end{aligned}$$

Pour résoudre ce problème avec un ordinateur quantique, nous exprimerons la fonction de coût comme la valeur d'espérance d'un observable. Cependant, les observables que Qiskit prend en charge nativement sont composés d'opérateurs de Pauli, qui ont des valeurs propres $1$ et $-1$ au lieu de $0$ et $1$. C'est pourquoi nous effectuerons le changement de variable suivant :

Où $\vec{x}=(x_0,x_1,\cdots ,x_{n-1})$. Nous pouvons utiliser la matrice d'adjacence $w$ pour accéder facilement aux poids de toutes les arêtes. Cela nous servira pour obtenir notre fonction de coût :

$$z_i = 1-2x_i \rightarrow x_i = \frac{1-z_i}{2}$$

Cela implique que :

$$\begin{array}{lcl} x_i=0 & \rightarrow & z_i=1 \\ x_i=1 & \rightarrow & z_i=-1.\end{array}$$

Alors la nouvelle fonction de coût que nous voulons maximiser est :

$$\begin{aligned} C(\vec{z}) & = \sum_{i,j=0}^n w_{ij} \bigg(\frac{1-z_i}{2}\bigg)\bigg(1-\frac{1-z_j}{2}\bigg)\\[1mm] & = \sum_{i,j=0}^n \frac{w_{ij}}{4} - \sum_{i,j=0}^n \frac{w_{ij}}{4} z_iz_j\\[1mm] & = \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2} - \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2} z_iz_j \end{aligned}$$

De plus, la tendance naturelle d'un ordinateur quantique est de trouver des minima (généralement l'énergie la plus basse) plutôt que des maxima, donc au lieu de maximiser $C(\vec{z})$ nous allons minimiser :

$$-C(\vec{z}) = \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2} z_iz_j - \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2}$$

Maintenant que nous avons une fonction de coût à minimiser dont les variables peuvent prendre les valeurs $-1$ et $1$, nous pouvons établir l'analogie suivante avec la Pauli $Z$ :

$$z_i \equiv Z_i = \overbrace{I}^{n-1}\otimes ... \otimes \overbrace{Z}^{i} \otimes ... \otimes \overbrace{I}^{0}$$

En d'autres termes, la variable $z_i$ sera équivalente à une porte $Z$ agissant sur le qubit $i$. De plus :

$$Z_i|x_{n-1}\cdots x_0\rangle = z_i|x_{n-1}\cdots x_0\rangle \rightarrow \langle x_{n-1}\cdots x_0 |Z_i|x_{n-1}\cdots x_0\rangle = z_i$$

Alors l'observable que nous allons considérer est :

$$\hat{H} = \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2} Z_iZ_j$$

auquel nous devrons ajouter le terme indépendant ensuite :

$$\texttt{offset} = - \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2}$$

L'opérateur est une combinaison linéaire de termes avec des opérateurs Z sur les nœuds connectés par une arête (rappelle-toi que le qubit 0 est le plus à droite) : $IIZZ + IZIZ + IZZI + ZIIZ + ZZII$. Une fois l'opérateur construit, l'ansatz pour l'algorithme QAOA peut être facilement construit en utilisant le circuit QAOAAnsatz de la bibliothèque de circuits de Qiskit.

from qiskit.circuit.library import QAOAAnsatz
from qiskit.quantum_info import SparsePauliOp

hamiltonian = SparsePauliOp.from_list(
    [("IIZZ", 1), ("IZIZ", 1), ("IZZI", 1), ("ZIIZ", 1), ("ZZII", 1)]
)

ansatz = QAOAAnsatz(hamiltonian, reps=2)
# Draw
ansatz.decompose(reps=3).draw("mpl")

# Sum the weights, and divide by 2

offset = -sum(edge[2] for edge in edges) / 2
print(f"""Offset: {offset}""")

Offset: -2.5

Avec l'Estimator de Runtime prenant directement un Hamiltonien et un ansatz paramétré, et retournant l'énergie nécessaire, la fonction de coût pour une instance de QAOA est assez simple :

def cost_func(params, ansatz, hamiltonian, estimator):
    """Return estimate of energy from estimator

    Parameters:
        params (ndarray): Array of ansatz parameters
        ansatz (QuantumCircuit): Parameterized ansatz circuit
        hamiltonian (SparsePauliOp): Operator representation of Hamiltonian
        estimator (Estimator): Estimator primitive instance

    Returns:
        float: Energy estimate
    """
    pub = (ansatz, hamiltonian, params)
    cost = estimator.run([pub]).result()[0].data.evs
    #    cost = estimator.run(ansatz, hamiltonian, parameter_values=params).result().values[0]
    return cost

import numpy as np

x0 = 2 * np.pi * np.random.rand(ansatz.num_parameters)

estimator = StatevectorEstimator()
cost = cost_func_vqe(x0, ansatz, hamiltonian, estimator)
print(cost)

1.473098768180865

# Estimated usage: < 1 min, benchmarked at 6 seconds on ibm_osaka, 5-23-24
# Load some necessary packages:

from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_ibm_runtime import Session, EstimatorV2 as Estimator

# Select the least busy backend:

backend = service.least_busy(
    operational=True, min_num_qubits=ansatz.num_qubits, simulator=False
)

# Or get a specific backend:
# backend = service.backend("ibm_brisbane")

# Use a pass manager to transpile the circuit and observable for the specific backend being used:

pm = generate_preset_pass_manager(backend=backend, optimization_level=1)
isa_ansatz = pm.run(ansatz)
isa_hamiltonian = hamiltonian.apply_layout(layout=isa_ansatz.layout)

# Set estimator options
estimator_options = EstimatorOptions(resilience_level=1, default_shots=10_000)

# Open a Runtime session:

with Session(backend=backend) as session:
    estimator = Estimator(mode=session, options=estimator_options)
    cost = cost_func_vqe(x0, isa_ansatz, isa_hamiltonian, estimator)

# Close session after done
session.close()
print(cost)

1.1120776913677988

Nous reviendrons à cet exemple dans la section Applications pour explorer comment utiliser un optimiseur pour parcourir l'espace de recherche. En termes généraux, cela inclut :

• Utiliser un optimiseur pour trouver les paramètres optimaux
• Lier les paramètres optimaux à l'ansatz pour trouver les valeurs propres
• Traduire les valeurs propres en notre définition du problème

Stratégie de mesure : vitesse versus précision

Comme nous l'avons mentionné, nous utilisons un ordinateur quantique bruité comme un oracle de boîte noire, où le bruit peut rendre les valeurs récupérées non déterministes, ce qui génère des fluctuations aléatoires qui, à leur tour, rendront difficile — voire empêcheront complètement — la convergence de certains optimiseurs vers une solution proposée. C'est un problème général que nous devons aborder à mesure que nous explorons progressivement l'utilité quantique et avançons vers l'avantage quantique :

Nous pouvons utiliser les options de suppression et d'atténuation des erreurs des primitives de Qiskit Runtime pour faire face au bruit et maximiser l'utilité des ordinateurs quantiques actuels.

Suppression des erreurs

La suppression des erreurs fait référence aux techniques utilisées pour optimiser et transformer un circuit pendant la compilation afin de minimiser les erreurs. Il s'agit d'une technique basique de gestion des erreurs qui introduit normalement un certain surcoût de prétraitement classique dans le temps d'exécution total. Ce surcoût inclut la transpilation des circuits pour les exécuter sur du matériel quantique en :

• Exprimant le circuit en utilisant les portes natives disponibles sur le système quantique
• Mappant les qubits virtuels sur les qubits physiques
• Ajoutant des opérations SWAP selon les exigences de connectivité
• Optimisant les portes à 1 et 2 qubits
• Ajoutant du découplage dynamique aux qubits inactifs pour prévenir les effets de la décohérence.

Les primitives permettent d'utiliser des techniques de suppression des erreurs en configurant l'option optimization_level et en sélectionnant des options avancées de transpilation. Dans un cours ultérieur, nous approfondirons différentes méthodes de construction de circuits pour améliorer les résultats, mais dans la plupart des cas nous recommandons d'utiliser optimization_level=3.

Nous allons visualiser l'intérêt d'augmenter l'optimisation dans le processus de transpilation en examinant un circuit d'exemple avec un comportement idéal simple.

from qiskit.circuit import Parameter, QuantumCircuit
from qiskit.quantum_info import SparsePauliOp

theta = Parameter("theta")

qc = QuantumCircuit(2)
qc.x(1)
qc.h(0)
qc.cp(theta, 0, 1)
qc.h(0)
observables = SparsePauliOp.from_list([("ZZ", 1)])

qc.draw("mpl")

Le circuit ci-dessus peut produire des valeurs d'espérance sinusoïdales de l'observable donné, à condition que nous insérions des phases couvrant un intervalle approprié, comme $[0,2\pi]$.

## Setup phases
import numpy as np

phases = np.linspace(0, 2 * np.pi, 50)

# phases need to be expressed as a list of lists in order to work
individual_phases = [[phase] for phase in phases]

Nous pouvons utiliser un simulateur pour démontrer l'utilité d'une transpilation optimisée. Plus loin, nous utiliserons à nouveau du matériel réel pour montrer l'utilité de l'atténuation des erreurs. Nous utiliserons QiskitRuntimeService pour obtenir un backend réel (dans ce cas, ibm_brisbane) et AerSimulator pour simuler ce backend, y compris son comportement de bruit.

from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_aer import AerSimulator

# get a real backend from the runtime service
service = QiskitRuntimeService()
backend = service.backend("ibm_brisbane")

# generate a simulator that mimics the real quantum system with the latest calibration results
backend_sim = AerSimulator.from_backend(backend)

Nous pouvons maintenant utiliser un pass manager pour transpiler le circuit vers le « jeu d'instructions d'architecture » ou ISA du backend. C'est une nouvelle exigence dans Qiskit Runtime : tous les circuits envoyés à un backend doivent se conformer aux contraintes de la cible du backend, ce qui signifie qu'ils doivent être écrits en termes de l'ISA du backend, c'est-à-dire le jeu d'instructions que le dispositif peut comprendre et exécuter. Ces contraintes de la cible sont définies par des facteurs comme les portes de base natives du dispositif, la connectivité entre qubits et, le cas échéant, les spécifications de pulses et autres durées d'instructions.

Note que dans ce cas, nous le ferons deux fois : une fois avec optimization_level = 0 et une autre avec le niveau fixé à 3. Dans chaque cas, nous utiliserons la primitive Estimator pour estimer les valeurs d'espérance de l'observable pour différentes valeurs de phase.

# Import estimator and specify that we are using the simulated backend:

from qiskit_ibm_runtime import EstimatorV2 as Estimator

estimator = Estimator(mode=backend_sim)

circuit = qc

# Use a pass manager to transpile the circuit and observable for the backend being simulated.
# Start with no optimization:

from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

pm = generate_preset_pass_manager(backend=backend_sim, optimization_level=0)
isa_circuit = pm.run(circuit)
isa_observables = observables.apply_layout(layout=isa_circuit.layout)

noisy_exp_values = []
pub = (isa_circuit, isa_observables, [individual_phases])
cost = estimator.run([pub]).result()[0].data.evs
noisy_exp_values = cost[0]

# Repeat above steps, but now with optimization = 3:

exp_values_with_opt_es = []
pm = generate_preset_pass_manager(backend=backend_sim, optimization_level=3)
isa_circuit = pm.run(circuit)
isa_observables = observables.apply_layout(layout=isa_circuit.layout)

pub = (isa_circuit, isa_observables, [individual_phases])
cost = estimator.run([pub]).result()[0].data.evs
exp_values_with_opt_es = cost[0]

Enfin, nous pouvons tracer les résultats. Nous voyons que la précision du calcul était assez bonne même sans optimisation, mais s'est clairement améliorée en passant le niveau d'optimisation à 3. Note que dans des circuits plus profonds et complexes, la différence entre les niveaux d'optimisation 0 et 3 sera probablement plus significative. Ceci est un circuit très simple qui sert de modèle jouet.

import matplotlib.pyplot as plt

plt.plot(phases, noisy_exp_values, "o", label="opt=0")
plt.plot(phases, exp_values_with_opt_es, "o", label="opt=3")
plt.plot(phases, 2 * np.sin(phases / 2) ** 2 - 1, label="ideal")
plt.ylabel("Expectation")
plt.legend()
plt.show()

Atténuation des erreurs

L'atténuation des erreurs fait référence à des techniques qui permettent de réduire les erreurs d'un circuit en modélisant le bruit du dispositif au moment de l'exécution. En général, cela implique un surcoût de prétraitement quantique lié à l'entraînement du modèle et un surcoût de post-traitement classique pour atténuer les erreurs dans les résultats bruts à partir du modèle généré.

L'option resilience_level des primitives de Qiskit Runtime spécifie le niveau de résilience face aux erreurs que l'on souhaite construire. Les niveaux plus élevés génèrent des résultats plus précis au prix de temps de traitement plus longs, dus au surcoût d'échantillonnage quantique plus important. Les niveaux de résilience peuvent être utilisés pour configurer l'équilibre entre coût et précision lors de l'application de l'atténuation des erreurs à ta requête de primitives.

Lors de l'implémentation de toute technique d'atténuation des erreurs, nous nous attendons à ce que le biais de nos résultats soit réduit par rapport au biais précédent non atténué. Dans certains cas, le biais peut même disparaître. Cependant, cela a un coût. À mesure que nous réduisons le biais dans nos quantités estimées, la variabilité statistique augmente (c'est-à-dire la variance), ce que nous pouvons compenser en augmentant encore le nombre de shots par circuit dans notre processus d'échantillonnage. Cela introduit un surcoût supplémentaire au-delà de celui nécessaire pour réduire le biais, c'est pourquoi cela n'est pas appliqué par défaut. Nous pouvons facilement activer ce comportement en ajustant le nombre de shots par circuit dans options.executions.shots, comme montré dans l'exemple ci-dessous.

Dans ce cours, nous explorerons ces modèles d'atténuation des erreurs à un niveau élevé pour illustrer l'atténuation que peuvent réaliser les primitives de Qiskit Runtime, sans avoir besoin de connaître tous les détails d'implémentation.

Extinction des erreurs de lecture par twirling (T-REx)

L'extinction des erreurs de lecture par twirling (T-REx) utilise une technique connue sous le nom de Pauli twirling pour réduire le bruit introduit pendant le processus de mesure quantique. Cette technique ne suppose aucune forme spécifique de bruit, ce qui la rend très générale et efficace.

Flux de travail général :

1. Acquérir des données pour l'état zéro avec des bits aléatoirement inversés (Pauli X avant la mesure)
2. Acquérir des données pour l'état désiré (bruité) avec des bits aléatoirement inversés (Pauli X avant la mesure)
3. Calculer la fonction spéciale pour chaque ensemble de données et diviser.

 

Nous pouvons configurer cela avec options.resilience_level = 1, comme montré dans l'exemple ci-dessous.

Extrapolation à bruit zéro

L'extrapolation à bruit zéro (ZNE, pour Zero Noise Extrapolation) fonctionne en amplifiant d'abord le bruit dans le circuit qui prépare l'état quantique désiré, en obtenant des mesures pour plusieurs niveaux de bruit différents et en utilisant ces mesures pour inférer le résultat sans bruit.

Flux de travail général :

1. Amplifier le bruit du circuit pour plusieurs facteurs de bruit
2. Exécuter chaque circuit avec le bruit amplifié
3. Extrapoler jusqu'à la limite de bruit zéro

 

Nous pouvons configurer cela avec options.resilience_level = 2. Nous pouvons l'optimiser davantage en explorant différentes valeurs de noise_factors, noise_amplifiers et extrapolators, mais cela dépasse le cadre de ce cours. Nous t'encourageons à expérimenter avec ces options décrites ici.

Chaque méthode entraîne son propre surcoût associé : un équilibre entre le nombre de calculs quantiques nécessaires (temps) et la précision des résultats :

$$\begin{array}{c|c|c|c} \text{Methods} & R=1 \text{, T-REx} & R=2 \text{, ZNE} \\[1mm] \hline \text{Assumptions} & \text{None} & \text{Ability to scale noise} \\[1mm] \text{Qubit overhead} & 1 & 1 \\[1mm] \text{Sampling overhead} & 2 & N_{\text{noise-factors}} \\[1mm] \text{Bias} & 0 & \mathcal{O}(\lambda^{N_{\text{noise-factors}}}) \\[1mm] \end{array}$$

Utilisation des options d'atténuation et de suppression de Qiskit Runtime

Ci-dessous, nous montrons comment calculer une valeur d'espérance en utilisant l'atténuation et la suppression des erreurs dans Qiskit Runtime. Nous pouvons exploiter exactement le même circuit et observable qu'avant, mais cette fois en maintenant le niveau d'optimisation fixé à 2 et en ajustant la résilience ou les techniques d'atténuation des erreurs utilisées. Ce processus d'atténuation des erreurs se produit plusieurs fois tout au long de la boucle d'optimisation.

Nous effectuons cette partie sur du matériel réel, car l'atténuation des erreurs n'est pas disponible sur les simulateurs.

# Estimated usage: 8 minutes, benchmarked on an Eagle processor, 5-23-24

from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_ibm_runtime import (
    Session,
    EstimatorOptions,
    EstimatorV2 as Estimator,
)

# We select the least busy backend

# Select the least busy backend
# backend = service.least_busy(
#    operational=True, min_num_qubits=ansatz.num_qubits, simulator=False
# )

# Or use a specific backend
backend = service.backend("ibm_brisbane")

# Initialize some variables to save the results from different runs:

exp_values_with_em0_es = []
exp_values_with_em1_es = []
exp_values_with_em2_es = []

# Use a pass manager to optimize the circuit and observables for the backend chosen:

pm = generate_preset_pass_manager(backend=backend, optimization_level=2)
isa_circuit = pm.run(circuit)
isa_observables = observables.apply_layout(layout=isa_circuit.layout)

# Open a session and run with no error mitigation:

estimator_options = EstimatorOptions(resilience_level=0, default_shots=10_000)

with Session(backend=backend) as session:
    estimator = Estimator(mode=session, options=estimator_options)

    pub = (isa_circuit, isa_observables, [individual_phases])
    cost = estimator.run([pub]).result()[0].data.evs

session.close()

exp_values_with_em0_es = cost[0]

# Open a session and run with resilience = 1:

estimator_options = EstimatorOptions(resilience_level=1, default_shots=10_000)

with Session(backend=backend) as session:
    estimator = Estimator(mode=session, options=estimator_options)

    pub = (isa_circuit, isa_observables, [individual_phases])
    cost = estimator.run([pub]).result()[0].data.evs

session.close()

exp_values_with_em1_es = cost[0]

# Open a session and run with resilience = 2:

estimator_options = EstimatorOptions(resilience_level=2, default_shots=10_000)

with Session(backend=backend) as session:
    estimator = Estimator(mode=session, options=estimator_options)

    pub = (isa_circuit, isa_observables, [individual_phases])
    cost = estimator.run([pub]).result()[0].data.evs

session.close()

exp_values_with_em2_es = cost[0]

Comme précédemment, nous pouvons tracer les valeurs d'espérance résultantes en fonction de l'angle de phase pour les trois niveaux d'atténuation des erreurs utilisés. Avec une certaine difficulté, on peut apprécier que l'atténuation des erreurs améliore légèrement les résultats. Encore une fois, cet effet est beaucoup plus prononcé dans des circuits plus profonds et complexes.

import matplotlib.pyplot as plt

plt.plot(phases, exp_values_with_em0_es, "o", label="unmitigated")
plt.plot(phases, exp_values_with_em1_es, "o", label="resil = 1")
plt.plot(phases, exp_values_with_em2_es, "o", label="resil = 2")
plt.plot(phases, 2 * np.sin(phases / 2) ** 2 - 1, label="ideal")
plt.ylabel("Expectation")
plt.legend()
plt.show()

Résumé

Avec cette leçon, tu as appris à créer une fonction de coût :

• Créer une fonction de coût
• Comment exploiter les primitives de Qiskit Runtime pour atténuer et supprimer le bruit
• Comment définir une stratégie de mesure pour optimiser vitesse versus précision

Voici notre flux de travail variationnel de haut niveau :

Notre fonction de coût s'exécute à chaque itération de la boucle d'optimisation. La prochaine leçon explorera comment l'optimiseur classique utilise l'évaluation de notre fonction de coût pour sélectionner de nouveaux paramètres.

import qiskit
import qiskit_ibm_runtime

print(qiskit.version.get_version_info())
print(qiskit_ibm_runtime.version.get_version_info())

1.1.0
0.23.0