Algorithmes quantiques : algorithmes quantiques variationnels

remarque

Takashi Imamichi (24 mai 2024)

Télécharge le PDF de la conférence originale. Remarque : certains extraits de code peuvent être devenus obsolètes puisqu'il s'agit d'images statiques.

Le temps de QPU approximatif pour exécuter cette expérience est de 9 minutes (testé sur un processeur Eagle).

(il se peut que ce notebook ne puisse pas s'exécuter dans le temps imparti sur l'Open Plan. Merci d'utiliser les ressources de calcul quantique avec discernement.)

1. Introduction

Ce tutoriel offre une vue d'ensemble d'un algorithme hybride quantique-classique, en se concentrant particulièrement sur le solveur variationnel quantique de valeurs propres (VQE) et l'algorithme d'optimisation approximative quantique (QAOA). L'objectif principal de ces algorithmes est de résoudre des problèmes d'optimisation en utilisant des circuits quantiques avec des portes quantiques paramétrées.

Malgré les avancées de l'informatique quantique, la présence de bruit dans les dispositifs quantiques actuels rend difficile l'extraction de résultats significatifs à partir de circuits quantiques profonds. Pour surmonter ce défi, VQE et QAOA adoptent une approche hybride quantique-classique, qui consiste à exécuter de manière itérative des circuits quantiques relativement courts grâce au calcul quantique, et à optimiser les paramètres des circuits quantiques paramétrés cibles grâce au calcul classique.

QAOA a le potentiel de fournir des solutions optimales aux problèmes cibles à l'échelle de l'utilité, grâce à l'application de diverses techniques de mitigation et de suppression d'erreurs. VQE a de nombreuses applications (comme la chimie quantique) pour lesquelles il est moins évolutif. Mais un certain nombre d'approches liées aux valeurs propres ont émergé pour compléter et enrichir VQE, notamment la diagonalisation de sous-espaces de Krylov et la diagonalisation quantique basée sur l'échantillonnage (SQD). Comprendre VQE constitue une première étape importante pour appréhender le large éventail d'algorithmes hybrides classique-quantique qui ont émergé.

Ce module décrit les concepts fondamentaux et l'implémentation de VQE et de QAOA. Les tutoriels suivants exploreront des sujets et techniques avancés pour faire monter en échelle ces algorithmes. Tu as besoin de la bibliothèque suivante dans ton environnement pour exécuter ce notebook. Si tu ne l'as pas encore installée, tu peux le faire en décommentant et en exécutant la cellule suivante.

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

# % pip install 'qiskit[visualization]' qiskit-ibm-runtime

2. Calcul de la valeur propre minimale d'un Hamiltonien simple

Nous allons commencer par appliquer VQE à un cas très simple, pour voir comment cela fonctionne. Nous allons calculer la valeur propre minimale de la matrice de Pauli $Z$ avec VQE. Nous allons commencer par importer quelques paquets généraux.

import numpy as np
from qiskit.circuit import ParameterVector, QuantumCircuit
from qiskit.primitives import StatevectorEstimator, StatevectorSampler
from qiskit.quantum_info import SparsePauliOp
from scipy.optimize import minimize

Nous définissons maintenant l'opérateur qui nous intéresse et l'affichons sous forme matricielle.

op = SparsePauliOp("Z")
op.to_matrix()

array([[ 1.+0.j,  0.+0.j],
       [ 0.+0.j, -1.+0.j]])

Il est facile d'obtenir les valeurs propres de manière classique, ce qui nous permet de vérifier notre travail. Cela peut devenir difficile à mesure que nous nous rapprochons de l'échelle de l'utilité. Ici, nous utilisons numpy.

# compute eigenvalues with numpy
result = np.linalg.eigh(op.to_matrix())
print("Eigenvalues:", result.eigenvalues)

Eigenvalues: [-1.  1.]

Pour obtenir les valeurs propres à l'aide d'un algorithme quantique variationnel, nous construisons un circuit avec des portes qui prennent des paramètres variationnels :

# define a variational form
param = ParameterVector("a", 3)
qc = QuantumCircuit(1, 1)
qc.u(param[0], param[1], param[2], 0)
qc_estimator = qc.copy()
qc.measure(0, 0)
qc.draw("mpl")

Si nous voulons estimer la valeur d'espérance d'un opérateur (comme $Z$), nous devons utiliser Estimator. Si nous voulons observer les états du système, nous utilisons Sampler.

sampler = StatevectorSampler()
estimator = StatevectorEstimator()

Nous pouvons calculer le nombre d'occurrences des chaînes de bits 0 et 1 avec des valeurs de paramètres aléatoires [1, 2, 3] à l'aide de Sampler.

# compute counts of bitstrings with random parameter values by Sampler
result = sampler.run([(qc, [1, 2, 3])]).result()
counts = result[0].data.c.get_counts()
counts

{'0': 783, '1': 241}

Nous savons que nous pouvons calculer la valeur d'espérance de Z par $\langle Z \rangle = p_0 - p_1$ avec les probabilités $\{0: p_0, 1: p_1\}$.

# compute the expectation value of Z based on the counts
(counts.get("0", 0) - counts.get("1", 0)) / sum(counts.values())

0.529296875

Ce circuit a fonctionné, mais les valeurs de paramètres choisies ne correspondaient pas à un état de très basse énergie (ou de faible valeur propre). La valeur propre obtenue est nettement supérieure au minimum. Le résultat est similaire avec Estimator.

Remarque : Estimator prend des circuits quantiques sans mesures.

result = estimator.run([(qc_estimator, op, [1, 2, 3])]).result()
result[0].data.evs

array(0.54030231)

Nous devrons parcourir l'espace des paramètres et trouver ceux qui donnent la valeur propre la plus basse. Nous créons une fonction qui reçoit les valeurs des paramètres de la forme variationnelle et renvoie la valeur d'espérance $\langle Z \rangle$.

# define a cost function to look for the minimum eigenvalue of Z
def cost(x):
    result = sampler.run([(qc, x)]).result()
    counts = result[0].data.c.get_counts()
    expval = (counts.get("0", 0) - counts.get("1", 0)) / sum(counts.values())
    # the following line shows the trajectory of the optimization
    print(expval, counts)
    return expval

Appliquons la fonction minimize de SciPy pour trouver la valeur propre minimale de Z.

# minimize the cost function with scipy's minimize
min_result = minimize(cost, [0, 0, 0], method="COBYLA", tol=1e-8)
min_result

1.0 {'0': 1024}
0.494140625 {'0': 765, '1': 259}
0.466796875 {'0': 751, '1': 273}
0.564453125 {'0': 801, '1': 223}
-0.4296875 {'1': 732, '0': 292}
-0.984375 {'1': 1016, '0': 8}
-0.8984375 {'1': 972, '0': 52}
-0.990234375 {'1': 1019, '0': 5}
-0.892578125 {'1': 969, '0': 55}
-0.986328125 {'1': 1017, '0': 7}
-0.861328125 {'1': 953, '0': 71}
-1.0 {'1': 1024}
-0.982421875 {'1': 1015, '0': 9}
-0.99609375 {'1': 1022, '0': 2}
-0.986328125 {'1': 1017, '0': 7}
-1.0 {'1': 1024}
-0.990234375 {'1': 1019, '0': 5}
-0.998046875 {'1': 1023, '0': 1}
-0.99609375 {'1': 1022, '0': 2}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-0.998046875 {'1': 1023, '0': 1}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-0.998046875 {'1': 1023, '0': 1}
-0.998046875 {'1': 1023, '0': 1}
-0.998046875 {'1': 1023, '0': 1}
-1.0 {'1': 1024}
-0.99609375 {'1': 1022, '0': 2}
-1.0 {'1': 1024}
-0.99609375 {'1': 1022, '0': 2}
-0.998046875 {'1': 1023, '0': 1}
-0.998046875 {'1': 1023, '0': 1}
-0.99609375 {'1': 1022, '0': 2}
-0.998046875 {'1': 1023, '0': 1}
-1.0 {'1': 1024}
-0.998046875 {'1': 1023, '0': 1}
-0.998046875 {'1': 1023, '0': 1}
-0.99609375 {'1': 1022, '0': 2}
-1.0 {'1': 1024}
-0.998046875 {'1': 1023, '0': 1}
-1.0 {'1': 1024}
-0.998046875 {'1': 1023, '0': 1}
-0.998046875 {'1': 1023, '0': 1}
-1.0 {'1': 1024}
-0.998046875 {'1': 1023, '0': 1}
-0.998046875 {'1': 1023, '0': 1}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-0.998046875 {'1': 1023, '0': 1}
-0.994140625 {'1': 1021, '0': 3}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}

message: Optimization terminated successfully.
 success: True
  status: 1
     fun: -1.0
       x: [ 3.182e+00  1.338e+00  1.664e-01]
    nfev: 63
   maxcv: 0.0

# check counts of bitstrings with the optimal parameters
result = sampler.run([(qc, min_result.x)]).result()
result[0].data.c.get_counts()

{'0': 1, '1': 1023}

2.1 Exercice

Calcule la valeur propre minimale de $Z \otimes Z$ avec VQE.

z2 = SparsePauliOp("ZZ")
print(z2)
print(z2.to_matrix())

SparsePauliOp(['ZZ'],
              coeffs=[1.+0.j])
[[ 1.+0.j  0.+0.j  0.+0.j  0.+0.j]
 [ 0.+0.j -1.+0.j  0.+0.j  0.+0.j]
 [ 0.+0.j  0.+0.j -1.+0.j  0.+0.j]
 [ 0.+0.j  0.+0.j  0.+0.j  1.+0.j]]

# compute eigenvalues with numpy

# define a variational form
# qc = ...

# compute counts of bitstrings with a random parameter values by Sampler
# result = sampler.run(...)
# result

# compute the expectation value of ZZ based on the counts

# verify the expectation value of ZZ with Estimator

# define a cost function to look for the minimum eigenvalue of ZZ
# def cost(x):
#    expval = ...
#    return expval

# minimize the cost function with scipy's minimize
# min_result = minimize(cost, [...], method="COBYLA", tol=1e-8)
# min_result

# check counts of bitstrings with the optimal parameter values
# result = sampler.run(qc, min_result.x).result()
# result

Solutions de l'exercice

Nous définissons l'opérateur qui nous intéresse et l'affichons sous forme matricielle.

z2 = SparsePauliOp("ZZ")
print(z2)
print(z2.to_matrix())

SparsePauliOp(['ZZ'],
              coeffs=[1.+0.j])
[[ 1.+0.j  0.+0.j  0.+0.j  0.+0.j]
 [ 0.+0.j -1.+0.j  0.+0.j  0.+0.j]
 [ 0.+0.j  0.+0.j -1.+0.j  0.+0.j]
 [ 0.+0.j  0.+0.j  0.+0.j  1.+0.j]]

Pour obtenir les valeurs propres à l'aide d'un algorithme quantique variationnel, nous construisons un circuit avec des portes qui prennent des paramètres variationnels :

# define a variational form
param = ParameterVector("a", 6)
qc = QuantumCircuit(2, 2)
qc.u(param[0], param[1], param[2], 0)
qc.u(param[3], param[4], param[5], 1)
qc_estimator = qc.copy()
qc.measure([0, 1], [0, 1])
qc.draw("mpl")

Si nous voulons estimer la valeur d'espérance d'un opérateur (comme $Z \otimes Z$), nous devrions utiliser Estimator. Si nous voulons observer les états du système, nous utilisons Sampler.

sampler = StatevectorSampler()
estimator = StatevectorEstimator()

# compute counts of bitstrings with random parameter values by Sampler
result = sampler.run([(qc, [1, 2, 3, 4, 5, 6])]).result()
counts = result[0].data.c.get_counts()
counts

{'10': 661, '11': 203, '01': 47, '00': 113}

# compute the expectation value of ZZ based on the counts
(
    counts.get("00", 0)
    - counts.get("01", 0)
    - counts.get("10", 0)
    + counts.get("11", 0)
) / sum(counts.values())

-0.3828125

Ce circuit a fonctionné, mais les valeurs de paramètres choisies ne correspondaient pas à un état de très basse énergie (ou de faible valeur propre). La valeur propre obtenue est nettement supérieure au minimum. Le résultat est similaire avec Estimator.

# verify the expectation value of ZZ with Estimator
result = estimator.run([(qc_estimator, z2, [1, 2, 3, 4, 5, 6])]).result()
result[0].data.evs

array(-0.35316516)

Nous devrons parcourir l'espace des paramètres et trouver ceux qui donnent la valeur propre la plus basse.

# define a cost function to look for the minimum eigenvalue of ZZ
def cost(x):
    result = sampler.run([(qc, x)]).result()
    counts = result[0].data.c.get_counts()
    expval = (
        counts.get("00", 0)
        - counts.get("01", 0)
        - counts.get("10", 0)
        + counts.get("11", 0)
    ) / sum(counts.values())
    print(expval, counts)
    return expval

# minimize the cost function with scipy's minimize
min_result = minimize(cost, [0, 0, 0, 0, 0, 0], method="COBYLA", tol=1e-8)
min_result

1.0 {'00': 1024}
0.578125 {'00': 808, '01': 216}
0.5234375 {'00': 780, '01': 244}
0.548828125 {'00': 793, '01': 231}
0.3515625 {'00': 637, '10': 164, '11': 55, '01': 168}
0.3359375 {'00': 638, '11': 46, '10': 174, '01': 166}
0.283203125 {'00': 602, '10': 181, '01': 186, '11': 55}
-0.087890625 {'01': 414, '00': 184, '10': 143, '11': 283}
0.236328125 {'10': 27, '11': 623, '01': 364, '00': 10}
-0.0625 {'11': 261, '01': 403, '00': 219, '10': 141}
0.248046875 {'01': 366, '11': 628, '00': 11, '10': 19}
-0.0625 {'10': 145, '11': 254, '01': 399, '00': 226}
0.228515625 {'01': 373, '11': 609, '00': 20, '10': 22}
0.0546875 {'11': 376, '10': 273, '01': 211, '00': 164}
-0.447265625 {'01': 731, '10': 10, '11': 267, '00': 16}
-0.71484375 {'01': 871, '11': 99, '00': 47, '10': 7}
-0.46484375 {'01': 741, '00': 253, '10': 9, '11': 21}
-0.87890625 {'01': 962, '00': 39, '11': 23}
-0.640625 {'00': 176, '01': 837, '11': 8, '10': 3}
-0.88671875 {'01': 966, '00': 41, '11': 17}
-0.994140625 {'01': 1021, '11': 3}
-0.91796875 {'01': 982, '11': 35, '00': 7}
-0.994140625 {'01': 1021, '11': 2, '00': 1}
-0.939453125 {'01': 993, '00': 31}
-0.990234375 {'01': 1019, '11': 5}
-0.90234375 {'01': 974, '00': 21, '11': 29}
-0.98046875 {'01': 1014, '11': 10}
-0.994140625 {'01': 1021, '00': 3}
-0.990234375 {'01': 1019, '11': 4, '00': 1}
-0.98828125 {'01': 1018, '11': 6}
-0.990234375 {'01': 1019, '11': 4, '00': 1}
-0.994140625 {'01': 1021, '11': 2, '00': 1}
-0.99609375 {'01': 1022, '11': 2}
-0.998046875 {'01': 1023, '00': 1}
-0.99609375 {'01': 1022, '00': 2}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '00': 1}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '00': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '00': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '00': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '00': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-0.99609375 {'01': 1022, '00': 1, '11': 1}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '00': 1}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-0.99609375 {'01': 1022, '11': 1, '00': 1}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '00': 1}
-0.994140625 {'01': 1021, '00': 3}
-0.998046875 {'01': 1023, '00': 1}
-0.99609375 {'01': 1022, '11': 2}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '00': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '11': 1}
-0.99609375 {'01': 1022, '11': 2}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}

message: Optimization terminated successfully.
 success: True
  status: 1
     fun: -0.998046875
       x: [ 3.167e+00  6.940e-01  1.033e+00 -2.894e-02  8.933e-01
            1.885e+00]
    nfev: 128
   maxcv: 0.0

message: Optimization terminated successfully.
 success: True
  status: 1
     fun: -0.99609375
       x: [ 3.098e+00 -5.402e-01  1.091e+00 -1.004e-02  3.615e-01
            6.913e-01]
    nfev: 115
   maxcv: 0.0

Nous avons obtenu une valeur propre extrêmement proche du minimum donné par numpy.

# check counts of bitstrings with the optimal parameters
result = sampler.run([(qc, min_result.x)]).result()
result[0].data.c.get_counts()

{'01': 1024}

3. Optimisation quantique avec les patterns Qiskit

Dans ce guide pratique, nous allons découvrir les patterns Qiskit et l'optimisation approximative quantique. Un pattern Qiskit est un ensemble intuitif et reproductible d'étapes permettant de mettre en œuvre un flux de travail de calcul quantique : Nous allons appliquer ces patterns au contexte de l'optimisation combinatoire et montrer comment résoudre le problème de coupe maximale (max-cut) à l'aide de l'algorithme d'optimisation approximative quantique (QAOA), une méthode itérative hybride (quantique-classique).

Remarque : cette partie sur le QAOA s'appuie sur la « Partie 1 : QAOA à petite échelle » du tutoriel Algorithme d'optimisation approximative quantique. Consulte ce tutoriel pour apprendre à le faire monter en échelle.

3.1 Pattern Qiskit pour l'optimisation (à petite échelle)

Cette partie utilise un problème de coupe maximale à petite échelle pour illustrer les étapes nécessaires à la résolution d'un problème d'optimisation à l'aide d'un ordinateur quantique. Le problème de coupe maximale est un problème d'optimisation difficile à résoudre (plus précisément, un problème NP-difficile) qui a de nombreuses applications en regroupement (clustering), en science des réseaux et en physique statistique. Ce tutoriel considère un graphe de nœuds reliés par des arêtes et vise à répartir les nœuds en deux ensembles en « coupant » des arêtes, de sorte que le nombre d'arêtes coupées soit maximal. Pour donner un peu de contexte avant de transposer ce problème en algorithme quantique, tu peux mieux comprendre comment le problème de coupe maximale devient un problème d'optimisation combinatoire classique en considérant d'abord la minimisation d'une fonction $f(x)$

$$\min_{x\in \{0, 1\}^n}f(x),$$

où l'entrée $x$ est un vecteur dont les composantes correspondent à chaque nœud d'un graphe. Ensuite, chacune de ces composantes est contrainte à être $0$ ou $1$ (ce qui représente l'inclusion ou non dans la coupe). Cet exemple à petite échelle utilise un graphe de $n=5$ nœuds.

Tu pourrais écrire une fonction pour une paire de nœuds $i,j$ qui indique si l'arête correspondante $(i,j)$ fait partie de la coupe. Par exemple, la fonction $x_i + x_j - 2 x_i x_j$ vaut 1 uniquement si $x_i$ ou $x_j$ (mais pas les deux) vaut 1 (ce qui signifie que l'arête fait partie de la coupe), et zéro sinon. Le problème consistant à maximiser les arêtes dans la coupe peut être formulé comme

$$\max_{x\in \{0, 1\}^n} \sum_{(i,j)} x_i + x_j - 2 x_i x_j,$$

qui peut se réécrire comme une minimisation de la forme

$$\min_{x\in \{0, 1\}^n} \sum_{(i,j)} 2 x_i x_j - x_i - x_j.$$

Le minimum de $f(x)$ correspond dans ce cas au nombre maximal d'arêtes traversées par la coupe. Comme tu peux le voir, rien ici ne relève encore du calcul quantique. Il faut reformuler ce problème sous une forme qu'un ordinateur quantique puisse comprendre. Initialise ton problème en créant un graphe de $n=5$ nœuds.

import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import rustworkx as rx
from rustworkx.visualization import mpl_draw

n = 5

graph = rx.PyGraph()
graph.add_nodes_from(range(1, n + 1))
edge_list = [
    (0, 1, 1.0),
    (0, 2, 1.0),
    (1, 2, 1.0),
    (1, 3, 1.0),
    (2, 4, 1.0),
    (3, 4, 1.0),
]
graph.add_edges_from(edge_list)
pos = rx.spring_layout(graph, seed=2)
mpl_draw(graph, node_size=600, pos=pos, with_labels=True, labels=str)

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

La première étape du pattern consiste à transposer le problème classique (le graphe) en circuits et opérateurs quantiques. Pour cela, il y a trois étapes principales :

1. Utiliser une série de reformulations mathématiques pour représenter ce problème à l'aide de la notation des problèmes d'optimisation binaire quadratique sans contrainte (QUBO).
2. Réécrire le problème d'optimisation sous forme d'un Hamiltonien dont l'état fondamental correspond à la solution qui minimise la fonction de coût.
3. Créer un circuit quantique qui préparera l'état fondamental de cet Hamiltonien par un processus similaire au recuit quantique.

Remarque : dans la méthodologie QAOA, l'objectif final est de disposer d'un opérateur (Hamiltonien) qui représente la fonction de coût de notre algorithme hybride, ainsi que d'un circuit paramétré (ansatz) qui représente des états quantiques correspondant à des solutions candidates du problème. Tu peux échantillonner ces états candidats, puis les évaluer à l'aide de la fonction de coût.

Graphe → problème d'optimisation

La première étape de la transposition consiste en un changement de notation. Ce qui suit exprime le problème en notation QUBO :

$$\min_{x\in \{0, 1\}^n}x^T Q x,$$

où $Q$ est une matrice $n\times n$ de nombres réels, $n$ correspond au nombre de nœuds de ton graphe, $x$ est le vecteur de variables binaires introduit ci-dessus, et $x^T$ indique la transposée du vecteur $x$.

Problem name: maxcut

Minimize
  2*x_1*x_2 + 2*x_1*x_3 + 2*x_2*x_3 + 2*x_2*x_4 + 2*x_3*x_5 + 2*x_4*x_5 - 2*x_1
  - 3*x_2 - 3*x_3 - 2*x_4 - 2*x_5

Subject to
  No constraints

  Binary variables (5)
    x_1 x_2 x_3 x_4 x_5

Problème d'optimisation → Hamiltonien

Tu peux ensuite reformuler le problème QUBO sous forme d'un Hamiltonien (ici, une matrice qui représente l'énergie d'un système) :

$$H_C=\sum_{ij}Q_{ij}Z_iZ_j + \sum_i b_iZ_i.$$

Étapes de reformulation du problème QAOA vers l'Hamiltonien

Pour montrer comment le problème QAOA peut être réécrit de cette manière, remplace d'abord les variables binaires $x_i$ par un nouvel ensemble de variables $z_i\in\{-1, 1\}$ via

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

Tu peux voir ici que si $x_i$ vaut $0$, alors $z_i$ doit valoir $1$. Lorsque les $x_i$ sont remplacés par les $z_i$ dans le problème d'optimisation ($x^TQx$), on peut obtenir une formulation équivalente.

$$x^TQx=\sum_{ij}Q_{ij}x_ix_j \\ =\frac{1}{4}\sum_{ij}Q_{ij}(1-z_i)(1-z_j) \\=\frac{1}{4}\sum_{ij}Q_{ij}z_iz_j-\frac{1}{4}\sum_{ij}(Q_{ij}+Q_{ji})z_i + \frac{n^2}{4}.$$

Si l'on définit maintenant $b_i=-\sum_{j}(Q_{ij}+Q_{ji})$, et que l'on retire le préfacteur ainsi que le terme constant $n^2$, on arrive aux deux formulations équivalentes du même problème d'optimisation.

$$min_{x\in\{0,1\}^n} x^TQx\Longleftrightarrow \min_{z\in\{-1,1\}^n}z^TQz + b^Tz$$

Ici, $b$ dépend de $Q$. Remarque que pour obtenir $z^TQz + b^Tz$, nous avons abandonné le facteur 1/4 ainsi qu'un décalage constant de $n^2$, qui ne jouent aucun rôle dans l'optimisation.

Maintenant, pour obtenir une formulation quantique du problème, on promeut les variables $z_i$ en une matrice de Pauli $Z$, telle qu'une matrice $2\times 2$ de la forme

$$Z_i = \begin{pmatrix}1 & 0 \\ 0 & -1\end{pmatrix}.$$

Lorsque tu substitues ces matrices dans le problème d'optimisation ci-dessus, tu obtiens l'Hamiltonien suivant

$$H_C=\sum_{ij}Q_{ij}Z_iZ_j + \sum_i b_iZ_i.$$

Rappelle-toi également que les matrices $Z$ sont plongées dans l'espace computationnel de l'ordinateur quantique, c'est-à-dire un espace de Hilbert de taille $2^n\times 2^n$. Par conséquent, tu dois comprendre des termes tels que $Z_iZ_j$ comme le produit tensoriel $Z_i\otimes Z_j$ plongé dans l'espace de Hilbert $2^n\times 2^n$. Par exemple, dans un problème comportant cinq variables de décision, le terme $Z_1Z_3$ signifie $I\otimes Z_3\otimes I\otimes Z_1\otimes I$, où $I$ est la matrice identité $2\times 2$.

Cet Hamiltonien est appelé l'Hamiltonien de la fonction de coût. Il a la propriété que son état fondamental correspond à la solution qui minimise la fonction de coût $f(x)$. Par conséquent, pour résoudre ton problème d'optimisation, tu dois maintenant préparer l'état fondamental de $H_C$ (ou un état ayant un fort recouvrement avec celui-ci) sur l'ordinateur quantique. Ensuite, l'échantillonnage à partir de cet état donnera, avec une forte probabilité, la solution de $\min~f(x)$.

def build_max_cut_operator(graph: rx.PyGraph) -> tuple[SparsePauliOp, float]:
    sp_list = []
    constant = 0
    for s, t in graph.edge_list():
        w = graph.get_edge_data(s, t)
        sp_list.append(("ZZ", [s, t], w / 2))
        constant -= 1 / 2
    return SparsePauliOp.from_sparse_list(
        sp_list, num_qubits=graph.num_nodes()
    ), constant

cost_hamiltonian, constant = build_max_cut_operator(graph)
print("Cost Function Hamiltonian:", cost_hamiltonian)
print("Constant:", constant)

Cost Function Hamiltonian: SparsePauliOp(['IIIZZ', 'IIZIZ', 'IIZZI', 'IZIZI', 'ZIZII', 'ZZIII'],
              coeffs=[0.5+0.j, 0.5+0.j, 0.5+0.j, 0.5+0.j, 0.5+0.j, 0.5+0.j])
Constant: -3.0

Hamiltonien → circuit quantique

L'Hamiltonien $H_C$ contient la définition quantique de ton problème. Tu peux maintenant créer un circuit quantique qui aidera à échantillonner de bonnes solutions à partir de l'ordinateur quantique. Le QAOA s'inspire du recuit quantique et applique des couches alternées d'opérateurs dans le circuit quantique.

L'idée générale est de partir de l'état fondamental d'un système connu, $H^{\otimes n}|0\rangle$ ci-dessus, puis de piloter le système vers l'état fondamental de l'opérateur de coût qui nous intéresse. Cela se fait en appliquant les opérateurs $\exp\{-i\gamma_k H_C\}$ et $\exp\{-i\beta_k H_m\}$ avec les angles $\gamma_1,...,\gamma_p$ et $\beta_1,...,\beta_p~$.

Le circuit quantique que tu génères est paramétré par $\gamma_i$ et $\beta_i$, de sorte que tu peux essayer différentes valeurs de $\gamma_i$ et $\beta_i$ et échantillonner l'état résultant. Dans ce cas, nous allons essayer un exemple avec 1 couche QAOA qui contient deux paramètres : $\gamma_1$ et $\beta_1$.

from qiskit.circuit.library import QAOAAnsatz

circuit = QAOAAnsatz(cost_operator=cost_hamiltonian, reps=1)
circuit.measure_all()
circuit.draw("mpl")

circuit.decompose(reps=3).draw("mpl", fold=-1)

circuit.parameters

ParameterView([ParameterVectorElement(β[0]), ParameterVectorElement(γ[0])])

3.3 Étape 2. Optimiser les circuits pour une exécution sur du matériel quantique

Le circuit ci-dessus contient une série d'abstractions utiles pour raisonner sur les algorithmes quantiques, mais qu'il n'est pas possible d'exécuter tel quel sur le matériel. Pour pouvoir s'exécuter sur une QPU, le circuit doit subir une série d'opérations qui constituent l'étape de transpilation ou d'optimisation de circuit du pattern.

La bibliothèque Qiskit propose une série de passes de transpilation couvrant un large éventail de transformations de circuits. Tu dois t'assurer que ton circuit est optimisé pour ton objectif.

La transpilation peut impliquer plusieurs étapes, telles que :

• Le mappage initial des qubits du circuit (comme les variables de décision) vers les qubits physiques du dispositif.
• Le déroulement (unrolling) des instructions du circuit quantique vers les instructions natives que le backend comprend.
• Le routage des qubits du circuit qui interagissent, vers des qubits physiques adjacents les uns aux autres.
• La suppression d'erreurs en ajoutant des portes à un seul qubit pour supprimer le bruit grâce au découplage dynamique.

Plus d'informations sur la transpilation sont disponibles dans notre documentation.

Le code suivant transforme et optimise le circuit abstrait dans un format prêt à l'exécution sur l'un des dispositifs accessibles via le cloud, en utilisant le service Qiskit IBM® Runtime.

Remarque : tu peux tester tes programmes localement grâce au « mode de test local » avant de les envoyer à de véritables ordinateurs quantiques. Plus d'informations sur le mode de test local sont disponibles dans la documentation.

from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

# Use a quantum device
service = QiskitRuntimeService()
backend = service.least_busy(min_num_qubits=127)
# backend = service.backend("ibm_kingston")

# You can test your programs locally with a fake backend (local testing mode)
# backend = FakeBrisbane()

print(backend)

# Create pass manager for transpilation
pm = generate_preset_pass_manager(optimization_level=3, backend=backend)

candidate_circuit = pm.run(circuit)
candidate_circuit.draw("mpl", fold=False, idle_wires=False)

service = QiskitRuntimeService(channel="ibm_quantum_platform")

<IBMBackend('ibm_strasbourg')>

3.4 Étape 3. Exécuter à l'aide des primitives Qiskit

Dans le flux de travail QAOA, les paramètres QAOA optimaux sont trouvés dans une boucle d'optimisation itérative, qui exécute une série d'évaluations de circuits et utilise un optimiseur classique pour trouver les paramètres $\beta_k$ et $\gamma_k$ optimaux. Cette boucle d'exécution se déroule selon les étapes suivantes :

1. Définir les paramètres initiaux
2. Instancier une nouvelle Session contenant la boucle d'optimisation et la primitive utilisée pour échantillonner le circuit
3. Une fois un ensemble optimal de paramètres trouvé, exécuter le circuit une dernière fois pour obtenir une distribution finale qui sera utilisée à l'étape de post-traitement.

Définir le circuit avec les paramètres initiaux

Nous commençons avec des paramètres choisis arbitrairement.

initial_gamma = np.pi
initial_beta = np.pi / 2
init_params = [initial_gamma, initial_beta]

Définir le backend et la primitive d'exécution

Utilise les primitives Qiskit Runtime pour interagir avec les backends IBM®. Les deux primitives sont Sampler et Estimator, et le choix de la primitive dépend du type de mesure que tu souhaites effectuer sur l'ordinateur quantique. Pour la minimisation de $H_C$, utilise Estimator, car la mesure de la fonction de coût est simplement la valeur d'espérance de $\langle H_C \rangle$.

Exécution

Les primitives offrent une variété de modes d'exécution pour planifier des charges de travail sur des dispositifs quantiques, et un flux de travail QAOA s'exécute de manière itérative dans une session. Tu peux brancher la fonction de coût basée sur Sampler dans la routine de minimisation de SciPy pour trouver les paramètres optimaux.

def cost_func_estimator(params, ansatz, hamiltonian, estimator):
    # transform the observable defined on virtual qubits to
    # an observable defined on all physical qubits
    isa_hamiltonian = hamiltonian.apply_layout(ansatz.layout)

    pub = (ansatz, isa_hamiltonian, params)
    job = estimator.run([pub])

    results = job.result()[0]
    cost = results.data.evs

    objective_func_vals.append(cost)

    return cost

from qiskit_ibm_runtime import Session, EstimatorV2
from scipy.optimize import minimize

objective_func_vals = []  # Global variable
with Session(backend=backend) as session:
    # If using qiskit-ibm-runtime<0.24.0, change `mode=` to `session=`
    estimator = EstimatorV2(mode=session)
    estimator.options.default_shots = 1000

    # Set simple error suppression/mitigation options
    estimator.options.dynamical_decoupling.enable = True
    estimator.options.dynamical_decoupling.sequence_type = "XY4"
    estimator.options.twirling.enable_gates = True
    estimator.options.twirling.num_randomizations = "auto"

    result = minimize(
        cost_func_estimator,
        init_params,
        args=(candidate_circuit, cost_hamiltonian, estimator),
        method="COBYLA",
        tol=1e-2,
    )
    print(result)

message: Optimization terminated successfully.
 success: True
  status: 1
     fun: -0.6557925874481715
       x: [ 2.873e+00  9.414e-01]
    nfev: 21
   maxcv: 0.0

L'optimiseur a pu réduire le coût et trouver de meilleurs paramètres pour le circuit.

plt.figure(figsize=(12, 6))
plt.plot(objective_func_vals)
plt.xlabel("Iteration")
plt.ylabel("Cost")
plt.show()

Une fois que tu as trouvé les paramètres optimaux du circuit, tu peux les assigner et échantillonner la distribution finale obtenue avec les paramètres optimisés. C'est ici que la primitive Sampler doit être utilisée, puisque c'est la distribution de probabilité des mesures de chaînes de bits qui correspond à la coupe optimale du graphe.

Remarque : cela implique de préparer un état quantique $\psi$ dans l'ordinateur, puis de le mesurer. Une mesure fait s'effondrer l'état en un seul état de base computationnelle — par exemple, 010101110000... — qui correspond à une solution candidate $x$ de notre problème d'optimisation initial ($\max f(x)$ ou $\min f(x)$ selon la tâche).

optimized_circuit = candidate_circuit.assign_parameters(result.x)
optimized_circuit.draw("mpl", fold=False, idle_wires=False)

from qiskit_ibm_runtime import SamplerV2

# If using qiskit-ibm-runtime<0.24.0, change `mode=` to `backend=`
sampler = SamplerV2(mode=backend)

# Set simple error suppression/mitigation options
sampler.options.dynamical_decoupling.enable = True
sampler.options.dynamical_decoupling.sequence_type = "XY4"
sampler.options.twirling.enable_gates = True
sampler.options.twirling.num_randomizations = "auto"

pub = (optimized_circuit,)
job = sampler.run([pub], shots=int(1e4))
counts_int = job.result()[0].data.meas.get_int_counts()
counts_bin = job.result()[0].data.meas.get_counts()
shots = sum(counts_int.values())
final_distribution_int = {key: val / shots for key, val in counts_int.items()}
final_distribution_bin = {key: val / shots for key, val in counts_bin.items()}
print(final_distribution_int)

{12: 0.0652, 31: 0.0089, 4: 0.0085, 13: 0.0731, 26: 0.0256, 28: 0.0246, 17: 0.0405, 25: 0.0591, 20: 0.031, 15: 0.0221, 8: 0.017, 21: 0.0371, 14: 0.0461, 16: 0.0229, 19: 0.0723, 23: 0.0199, 22: 0.0478, 18: 0.0708, 24: 0.0165, 6: 0.0525, 7: 0.0155, 5: 0.0245, 3: 0.0231, 29: 0.0121, 30: 0.0062, 10: 0.0363, 1: 0.0097, 9: 0.042, 27: 0.0094, 11: 0.0349, 0: 0.0129, 2: 0.0119}

3.5 Étape 4. Post-traitement, retour du résultat au format classique

L'étape de post-traitement interprète la sortie de l'échantillonnage pour renvoyer une solution à ton problème d'origine. Dans ce cas, tu t'intéresses à la chaîne de bits la plus probable, car c'est elle qui détermine la coupe optimale. Les symétries du problème permettent quatre solutions possibles, et le processus d'échantillonnage en renverra une avec une probabilité légèrement plus élevée, mais tu peux voir dans la distribution tracée ci-dessous que quatre des chaînes de bits sont nettement plus probables que le reste.

# auxiliary functions to sample most likely bitstring
def to_bitstring(integer, num_bits):
    result = np.binary_repr(integer, width=num_bits)
    return [int(digit) for digit in result]

keys = list(final_distribution_int.keys())
values = list(final_distribution_int.values())
most_likely = keys[np.argmax(np.abs(values))]
most_likely_bitstring = to_bitstring(most_likely, len(graph))
most_likely_bitstring.reverse()

print("Result bitstring:", most_likely_bitstring)

Result bitstring: [1, 0, 1, 1, 0]

import matplotlib.pyplot as plt

matplotlib.rcParams.update({"font.size": 10})
final_bits = final_distribution_bin
values = np.abs(list(final_bits.values()))
top_4_values = sorted(values, reverse=True)[:4]
positions = []
for value in top_4_values:
    positions.append(np.where(values == value)[0])
fig = plt.figure(figsize=(11, 6))
ax = fig.add_subplot(1, 1, 1)
plt.xticks(rotation=45)
plt.title("Result Distribution")
plt.xlabel("Bitstrings (reversed)")
plt.ylabel("Probability")
ax.bar(list(final_bits.keys()), list(final_bits.values()), color="tab:grey")
for p in positions:
    ax.get_children()[p[0].item()].set_color("tab:purple")
plt.show()

Visualiser la meilleure coupe

À partir de la chaîne de bits optimale, tu peux ensuite visualiser cette coupe sur le graphe d'origine.

colors = ["tab:grey" if i == 0 else "tab:purple" for i in most_likely_bitstring]
mpl_draw(graph, node_size=600, pos=pos, with_labels=True, labels=str, node_color=colors)

Et calculons la valeur de la coupe. La solution n'est pas optimale à cause du bruit (la valeur de coupe de la solution optimale est 5).

from typing import Sequence

def evaluate_sample(x: Sequence[int], graph: rx.PyGraph) -> float:
    assert len(x) == len(
        list(graph.nodes())
    ), "The length of x must coincide with the number of nodes in the graph."
    return sum(
        x[u] * (1 - x[v]) + x[v] * (1 - x[u]) for u, v in list(graph.edge_list())
    )

cut_value = evaluate_sample(most_likely_bitstring, graph)
print("The value of the cut is:", cut_value)

The value of the cut is: 5

Ceci conclut le tutoriel sur le QAOA à petite échelle. Tu apprendras à adapter le QAOA à l'échelle de l'utilité dans la « Partie 2 : passons à l'échelle ! » du tutoriel Algorithme d'optimisation approximative quantique.

# Check Qiskit version
import qiskit

qiskit.__version__

'2.0.2'