Algorithmes quantiques : algorithmes quantiques variationnels

1. Introduction

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

Malgré les progrès 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 relever ce défi, VQE et QAOA adoptent une approche hybride quantique-classique, qui implique l'exécution itérative de circuits quantiques relativement courts à l'aide de l'informatique quantique et l'optimisation des paramètres des circuits quantiques paramétrés cibles à l'aide de l'informatique classique.

La QAOA a le potentiel de fournir des solutions optimales aux problèmes ciblés à l'échelle d'un service public, grâce à l'application de diverses techniques d'atténuation et de suppression des erreurs. La VQE a de nombreuses applications (comme la chimie quantique) dans lesquelles elle est moins évolutive. Mais un certain nombre d'approches liées aux valeurs propres sont apparues pour compléter et accroître l'EQV, notamment la diagonalisation du sous-espace de Krylov et la diagonalisation quantique basée sur l'échantillonnage (SQD). La compréhension de la VQE est une première étape importante dans la compréhension du large éventail d'algorithmes hybrides classiques-quantiques qui ont vu le jour.

Ce module décrit les concepts fondamentaux et la mise en œuvre de VQE et de QAOA. D'autres tutoriels exploreront des sujets avancés et des techniques pour augmenter la taille de ces algorithmes.

Vous avez besoin de la bibliothèque suivante dans votre environnement pour exécuter ce cahier. Si vous ne l'avez pas encore installé, vous pouvez le faire en décommentant et en exécutant la cellule suivante.

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

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

Nous commencerons par appliquer l'EQV à un cas très simple, afin de voir comment il fonctionne. Nous calculerons la valeur propre minimale de la matrice de Pauli $Z$ avec VQE. Nous commencerons 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 allons maintenant définir l'opérateur qui nous intéresse et le représenter sous forme de matrice.

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

Output:

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 pourrait s'avérer difficile au fur et à mesure que nous nous rapprochons de l'utilité. Nous utilisons ici numpy.

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

Output:

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

Output:

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

sampler = StatevectorSampler()
estimator = StatevectorEstimator()

Nous pouvons calculer le nombre de 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

Output:

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

Nous savons que nous pouvons calculer la valeur espérée 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())

Output:

0.529296875

Ce circuit a fonctionné, mais les valeurs des paramètres choisis ne correspondaient pas à un état de très faible énergie (ou de faible valeur propre). La valeur propre obtenue est nettement supérieure à la valeur minimale. Le résultat est similaire lorsque l'on utilise l'estimateur.

Notez que l'Estimateur prend des circuits quantiques sans mesures.

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

Output:

array(0.54030231)

Nous devrons rechercher parmi les paramètres ceux qui produisent la valeur propre la plus faible. Nous créons une fonction qui reçoit les valeurs des paramètres de la forme variationnelle et renvoie la valeur de l'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 SciPy's minimize 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

Output:

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

Output:

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

2.1 Exercice

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

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

Output:

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 le représentons sous forme de matrice.

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

Output:

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

Output:

Si nous voulons estimer la valeur de l'espérance d'un opérateur (comme $Z \otimes Z$ ), nous utiliserons Estimator. Si nous voulons examiner les états du système, nous utilisons l'échantillonneur.

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

Output:

{'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())

Output:

-0.3828125

Ce circuit a fonctionné, mais les valeurs des paramètres choisis ne correspondaient pas à un état de très faible énergie (ou de faible valeur propre). La valeur propre obtenue est nettement supérieure à la valeur minimale. Le résultat est similaire lorsque l'on utilise l'estimateur.

# 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

Output:

array(-0.35316516)

Nous devrons rechercher parmi les paramètres ceux qui produisent la valeur propre la plus faible.

# 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

Output:

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

Output:

{'01': 1024}

3. Optimisation quantique avec les modèles Qiskit

Dans ce guide pratique, nous apprendrons à connaître les motifs Qiskit et l'optimisation approximative quantique. Un modèle Qiskit est un ensemble d'étapes intuitives et reproductibles pour la mise en œuvre d'un flux de travail d'informatique quantique :

Nous appliquerons ces modèles au domaine de l'optimisation combinatoire et montrerons comment résoudre le problème du coupure maximale à l'aide de l 'algorithme d'optimisation approximative quantique (QAOA), une méthode itérative hybride (quantique-classique).

Notez que cette partie de QAOA est basée sur la "Partie 1 : QAOA à petite échelle" du tutoriel sur l'algorithme d'optimisation approximative quantique. Voir le tutoriel pour savoir comment l'agrandir.

3.1 Modèle Qiskit (à petite échelle) pour l'optimisation

Cette partie s'appuiera sur un problème de « max-cut » à 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 du « max-cut » est un problème d'optimisation difficile à résoudre (plus précisément, il s'agit d'un problème NP-difficile) qui trouve de nombreuses applications dans le regroupement de données, la science des réseaux et la physique statistique. Ce tutoriel porte sur un graphe composé de nœuds reliés par des arêtes et vise à partitionner ces nœuds en deux ensembles en « coupant » des arêtes, de manière à maximiser le nombre d'arêtes coupées.

Pour replacer ce problème dans son contexte avant de le transposer en algorithme quantique, vous comprendrez mieux comment le problème du coupé maximal se transforme en un problème d'optimisation combinatoire classique en considérant d'abord la minimisation d'une fonction $f(x)$

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

Vous pouvez écrire une fonction pour une paire de nœuds $i,j$ qui indique si l'arête correspondante $(i,j)$ est dans la coupe. Par exemple, la fonction $x_i + x_j - 2 x_i x_j$ n'est égale à 1 que si l'une des valeurs $x_i$ ou $x_j$ est égale à 1 (ce qui signifie que le bord est dans la coupe) et à zéro dans le cas contraire. Le problème de la maximisation des arêtes dans la coupe peut être formulé comme suit

qui peut être réécrite comme une minimisation de la forme

Dans ce cas, le minimum de $f(x)$ est atteint lorsque le nombre d'arêtes traversées par la coupe est maximal. Comme vous pouvez le constater, il n'y a encore rien en rapport avec l'informatique quantique. Vous devez reformuler ce problème de manière à ce qu'un ordinateur quantique puisse le comprendre.

Initialisez votre problème en créant un graphe avec $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)

Output:

3.2 Étape 1. Mapper les entrées classiques à un problème quantique

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

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 la forme d'un hamiltonien pour lequel l'état fondamental correspond à la solution qui minimise la fonction de coût.
3. Créez 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, vous souhaitez 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 les états quantiques avec des solutions candidates au problème. Vous pouvez prélever un échantillon de ces états candidats, puis les évaluer à l'aide de la fonction de coût.

Graphique → problème d'optimisation

La première étape de la mise en correspondance est un changement de notation. Le problème est exprimé ci-dessous en notation QUBO :

où $Q$ est une matrice $n\times n$ de nombres réels, $n$ correspond au nombre de nœuds dans votre 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

Vous pouvez alors reformuler le problème QUBO sous la forme d'un hamiltonien (ici, une matrice qui représente l'énergie d'un système) :

Étapes de reformulation du problème QAOA à l'hamiltonien

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

On voit ici que si $x_i$ est $0$, alors $z_i$ doit être $1$. Lorsque les $x_i$ sont remplacés par les $z_i$ dans le problème d'optimisation ( $x^TQx$ ), on obtient une formulation équivalente.

Si nous définissons $b_i=-\sum_{j}(Q_{ij}+Q_{ji})$, supprimons le préfacteur et le terme constant $n^2$, nous obtenons deux formulations équivalentes du même problème d'optimisation.

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

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

En remplaçant ces matrices dans le problème d'optimisation ci-dessus, on obtient l'hamiltonien suivant

Rappelez-vous également que les matrices $Z$ sont intégrées dans l'espace de calcul de l'ordinateur quantique, c'est-à-dire un espace de Hilbert de taille $2^n\times 2^n$. Par conséquent, vous devez comprendre des termes tels que $Z_iZ_j$ comme le produit tensoriel $Z_i\otimes Z_j$ intégré 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$ est compris comme signifiant $I\otimes Z_3\otimes I\otimes Z_1\otimes I$ où $I$ est la matrice d'identité $2\times 2$.

Cet hamiltonien est appelé fonction de coût Hamiltonien. Il a la propriété que son état fondamental correspond à la solution que minimise la fonction de coût $f(x)$. Par conséquent, pour résoudre votre problème d'optimisation, vous devez maintenant préparer l'état fondamental de $H_C$ (ou un état ayant un fort recouvrement avec lui) sur l'ordinateur quantique. L'échantillonnage de cet état permet alors, avec une forte probabilité, d'obtenir 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)

Output:

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

Circuit quantique hamiltonien

L'hamiltonien $H_C$ contient la définition quantique de votre problème. Vous pouvez maintenant créer un circuit quantique qui permettra d' échantillonner les bonnes solutions 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 d'orienter le système vers l'état fondamental de l'opérateur de coût qui vous intéresse. Pour ce faire, les opérateurs $\exp\{-i\gamma_k H_C\}$ et $\exp\{-i\beta_k H_m\}$ sont appliqués avec les angles $\gamma_1,...,\gamma_p$ et $\beta_1,...,\beta_p~$.

Le circuit quantique que vous générez est paramétré par $\gamma_i$ et $\beta_i$, de sorte que vous pouvez 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 une couche de 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")

Output:

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

Output:

circuit.parameters

Output:

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

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

Le circuit ci-dessus contient une série d'abstractions utiles pour réfléchir aux algorithmes quantiques, mais impossibles à exécuter sur le matériel. Pour pouvoir fonctionner sur une QPU, le circuit doit subir une série d'opérations qui constituent l'étape de transpilation ou d' optimisation du circuit du modèle.

La bibliothèque Qiskit offre une série de passes de transpilation qui répondent à un large éventail de transformations de circuits. Vous devez vous assurer que votre circuit est optimisé pour votre objectif.

La transpilation peut comporter plusieurs étapes, telles que

• Mappage initial des qubits du circuit (tels que les variables de décision) aux qubits physiques de l'appareil.
• Déroulement des instructions dans le circuit quantique vers les instructions natives du matériel que le backend comprend.
• Routage de tous les qubits du circuit qui interagissent vers des qubits physiques adjacents.
• Suppression des erreurs par l'ajout de portes à qubit unique pour supprimer le bruit avec découplage dynamique.

De plus amples informations sur la transpilation sont disponibles dans notre documentation.

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

Notez que vous pouvez tester vos programmes localement grâce au "mode de test local" avant de les envoyer à de véritables ordinateurs quantiques. De plus amples 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)

Output:

  service = QiskitRuntimeService(channel="ibm_quantum_platform")

<IBMBackend('ibm_strasbourg')>

3.4 Étape 3. Exécuter à l'aide des primitives « IBM Quantum »

Dans le flux de travail du QAOA, les paramètres optimaux du QAOA 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 optimaux $\beta_k$ et $\gamma_k$. Cette boucle d'exécution est exécutée en suivant les étapes suivantes :

1. Définir les paramètres initiaux
2. Instanciation d'un nouveau site Session contenant la boucle d'optimisation et la primitive utilisée pour échantillonner le circuit
3. Une fois qu'un ensemble optimal de paramètres a été trouvé, exécuter le circuit une dernière fois pour obtenir une distribution finale qui sera utilisée dans 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

Utilisez les primitives de la bibliothèque « IBM Quantum » pour interagir avec les backends de « IBM® ». Ces deux primitives sont « Sampler » et « Estimator »; le choix de la primitive dépend du type de mesure que vous souhaitez effectuer sur l'ordinateur quantique. Pour minimiser la fonction de coût « $H_C$ », utilisez l'estimateur « Estimator », car la mesure de cette fonction de coût correspond simplement à l'espérance de « $\langle H_C \rangle$ ».

Exécuter

Les primitives offrent une variété de modes d'exécution pour planifier les charges de travail sur les dispositifs quantiques, et un flux de travail QAOA s'exécute de manière itérative dans une session.

Vous pouvez introduire la fonction de coût basée sur l'échantillonneur dans la routine de minimisation 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)

Output:

 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 permis de réduire le coût et de 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()

Output:

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

Remarque : il s'agit de préparer un état quantique $\psi$ dans l'ordinateur et de le mesurer. Une mesure réduira l'état en un seul état de base de calcul - par exemple, 010101110000... - qui correspond à une solution candidate $x$ à notre problème d'optimisation initial ( $\max f(x)$ ou $\min f(x)$ en fonction de la tâche).

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

Output:

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)

Output:

{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, renvoyer le résultat dans un format classique

L'étape de post-traitement interprète le résultat de l'échantillonnage afin de fournir une solution au problème initial. Dans ce cas, vous vous intéressez à la chaîne de bits ayant la probabilité la plus élevée, car elle détermine le découpage optimal. Les symétries du problème permettent quatre solutions possibles, et le processus d'échantillonnage renverra l'une d'entre elles avec une probabilité légèrement supérieure, mais vous pouvez voir dans la distribution graphique ci-dessous que quatre des chaînes de bits sont nettement plus probables que les autres.

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

Output:

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

Output:

Visualiser la meilleure coupe

À partir de la chaîne de bits optimale, vous pouvez ensuite visualiser cette coupe sur le graphique original.

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)

Output:

Et calculer la valeur de la coupe. La solution n'est pas optimale en raison du bruit (la valeur de coupure de la solution optimale est de 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)

Output:

The value of the cut is: 5

Ceci conclut le tutoriel sur l'AQAO à petite échelle. Vous apprendrez comment adapter le QAOA à l'échelle d'un service public dans la "Partie 2 : passer à l'échelle supérieure" du tutoriel sur l 'algorithme d'optimisation approximative quantique.

# Check Qiskit version
import qiskit

qiskit.__version__

Output:

'2.0.2'