Algoritmi cuantici: Algoritmi cuantici variaționali

notă

Takashi Imamichi (24 May 2024)

Descarcă pdf-ul al prelegerii originale. Rețineți că unele fragmente de cod ar putea deveni depreciate deoarece acestea sunt imagini statice.

Approximate QPU time to run this experiment is 9 minutes (tested on an Eagle processor).

(this notebook might not evaluate in the time allowed on the Open Plan. Please use quantum computing resources wisely.)

1. Introducere

Acest tutorial oferă o prezentare generală a unui algoritm hibrid cuantic-clasic, concentrându-se în special pe variational quantum eigensolver (VQE) și pe quantum approximate optimization algorithm (QAOA). Obiectivul principal al acestor algoritmi este să abordeze probleme de optimizare folosind Circuit-uri cuantice cu Gate-uri cuantice parametrizate.

În ciuda progreselor din domeniul calculului cuantic, prezența zgomotului în dispozitivele cuantice actuale face dificilă extragerea unor rezultate semnificative din Circuit-uri cuantice adânci. Pentru a depăși această provocare, VQE și QAOA adoptă o abordare hibridă cuantic-clasică, care implică executarea iterativă a unor Circuit-uri cuantice relativ scurte folosind calculul cuantic și optimizarea parametrilor Circuit-urilor cuantice parametrizate țintă folosind calculul clasic.

QAOA are potențialul de a furniza soluții optime pentru problemele țintă la scară de utilitate, datorită aplicării diverselor tehnici de atenuare și suprimare a erorilor. VQE are multe aplicații (precum chimia cuantică) în care este mai puțin scalabil. Totuși, au apărut o serie de abordări legate de valori proprii care vin să completeze și să augmenteze VQE, inclusiv diagonalizarea subspațiului Krylov și diagonalizarea cuantică bazată pe eșantionare (SQD). Înțelegerea VQE reprezintă un prim pas important în înțelegerea gamei largi de algoritmi hibrid cuantic-clasici care au apărut.

Acest modul descrie conceptele fundamentale și implementarea VQE și QAOA. Tutorialele ulterioare vor explora subiecte avansate și tehnici pentru scalarea acestor algoritmi. Ai nevoie de următoarea bibliotecă în mediul tău pentru a rula acest notebook. Dacă nu ai instalat-o încă, o poți instala decomentând și rulând celula de mai jos.

# 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. Calculul valorii proprii minime a unui Hamiltonian simplu

Vom începe prin a aplica VQE unui caz foarte simplu, pentru a vedea cum funcționează. Vom calcula valoarea proprie minimă a matricei Pauli $Z$ cu VQE. Vom începe prin importarea unor pachete generale.

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

Acum definim operatorul de interes și îl vizualizăm sub formă matriceală.

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

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

Este ușor să obținem valorile proprii în mod clasic, astfel că ne putem verifica munca. Acest lucru poate deveni dificil pe măsură ce scalăm spre utilitate. Aici folosim numpy.

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

Eigenvalues: [-1.  1.]

Pentru a obține valorile proprii folosind un algoritm cuantic variațional, construim un Circuit cu Gate-uri care primesc parametri variaționali:

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

Dacă vrei să estimezi valoarea așteptată a unui operator (precum $Z$), ar trebui să folosești Estimator. Dacă vrei să examinezi stările sistemului, folosești Sampler.

sampler = StatevectorSampler()
estimator = StatevectorEstimator()

Putem calcula numărătorile de șiruri de biți 0 și 1 cu valori aleatorii ale parametrilor [1, 2, 3] folosind 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}

Știm că putem calcula valoarea așteptată a lui Z prin $\langle Z \rangle = p_0 - p_1$ cu probabilitățile $\{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

Acest Circuit a funcționat, dar valorile parametrilor alese nu corespundeau unei stări cu energie (sau valoare proprie) foarte scăzută. Valoarea proprie obținută este considerabil mai mare decât minimul. Rezultatul este similar când se folosește Estimator.

Reține că Estimator primește Circuit-uri cuantice fără măsurători.

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

array(0.54030231)

Va trebui să căutăm prin parametri și să îi găsim pe cei care produc cea mai mică valoare proprie. Definim o funcție care primește valorile parametrilor formei variaționali și returnează valoarea așteptată $\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

Să aplicăm funcția minimize din SciPy pentru a găsi valoarea proprie minimă a lui 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 Exercițiu

Calculează valoarea proprie minimă a lui $Z \otimes Z$ cu 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

Soluțiile exercițiului

Definim operatorul de interes și îl vizualizăm sub formă matriceală.

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

Pentru a obține valorile proprii folosind un algoritm cuantic variațional, construim un Circuit cu Gate-uri care primesc parametri variaționali:

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

Dacă vrei să estimezi valoarea așteptată a unui operator (precum $Z \otimes Z$), ai folosi Estimator. Dacă vrei să examinezi stările sistemului, folosești 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

Acest Circuit a funcționat, dar valorile parametrilor alese nu corespundeau unei stări cu energie (sau valoare proprie) foarte scăzută. Valoarea proprie obținută este considerabil mai mare decât minimul. Rezultatul este similar când se folosește 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)

Va trebui să căutăm prin parametri și să îi găsim pe cei care produc cea mai mică valoare proprie.

# 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

Am obținut o valoare proprie extrem de apropiată de minimul furnizat de 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. Optimizare cuantică cu tipare Qiskit

În acest tutorial vom afla despre tiparele Qiskit și optimizarea cuantică aproximativă. Un tipar Qiskit este un set intuitiv și repetabil de pași pentru implementarea unui flux de lucru de calcul cuantic: Vom aplica tiparele în contextul optimizării combinatoriale și vom arăta cum să rezolvăm problema Max-Cut folosind Quantum Approximate Optimization Algorithm (QAOA), o metodă iterativă hibridă (cuantică-clasică).

Reține că această parte despre QAOA se bazează pe „Partea 1: QAOA la scară mică" din tutorialul Quantum approximate optimization algorithm. Consultă tutorialul pentru a afla cum să scalezi.

3.1 Tiparul Qiskit (la scară mică) pentru optimizare

Această parte va folosi o problemă Max-Cut la scară mică pentru a ilustra pașii necesari rezolvării unei probleme de optimizare folosind un calculator cuantic. Problema MaxCut este o problemă de optimizare greu de rezolvat (mai precis, este o problemă NP-dificilă) cu o serie de aplicații diferite în clustering, știința rețelelor și fizica statistică. Acest tutorial consideră un graf de noduri conectate prin muchii și urmărește să împartă nodurile în două seturi prin „tăierea" muchiilor, astfel încât numărul de muchii tăiate să fie maximizat. Pentru a oferi un context înainte de maparea acestei probleme pe un algoritm cuantic, poți înțelege mai bine cum devine problema MaxCut o problemă de optimizare combinatorie clasică luând în considerare mai întâi minimizarea unei funcții $f(x)$

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

unde inputul $x$ este un vector ale cărui componente corespund fiecărui nod al unui graf. Apoi, constrângi fiecare dintre aceste componente să fie fie $0$, fie $1$ (care reprezintă includerea sau neincluderea în tăietură). Acest exemplu la scară mică folosește un graf cu $n=5$ noduri.

Poți scrie o funcție pentru o pereche de noduri $i,j$ care indică dacă muchia corespunzătoare $(i,j)$ se află în tăietură. De exemplu, funcția $x_i + x_j - 2 x_i x_j$ este 1 doar dacă unul dintre $x_i$ sau $x_j$ este 1 (ceea ce înseamnă că muchia se află în tăietură) și zero altfel. Problema maximizării muchiilor din tăietură poate fi formulată astfel:

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

care poate fi rescrisă ca o minimizare de forma

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

Minimul lui $f(x)$ în acest caz apare când numărul de muchii traversate de tăietură este maximal. După cum poți vedea, nu există nimic legat de calculul cuantic încă. Trebuie să reformulezi această problemă în ceva ce un calculator cuantic poate înțelege. Inițializează problema creând un graf cu $n=5$ noduri.

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 Pasul 1. Maparea inputurilor clasice la o problemă cuantică

Primul pas al tiparului este de a mapa problema clasică (graful) în Circuit-uri și operatori cuantici. Pentru a face acest lucru, există trei pași principali de urmat:

1. Utilizarea unei serii de reformulări matematice pentru a reprezenta această problemă folosind notația problemelor de Optimizare Binară Pătratică Fără Constrângeri (QUBO).
2. Rescrierea problemei de optimizare ca un Hamiltonian pentru care starea fundamentală corespunde soluției care minimizează funcția de cost.
3. Crearea unui Circuit cuantic care va pregăti starea fundamentală a acestui Hamiltonian printr-un proces similar cu recoacerea cuantică.

Notă: În metodologia QAOA, vrei în cele din urmă să ai un operator (Hamiltonian) care reprezintă funcția de cost a algoritmului nostru hibrid, precum și un Circuit parametrizat (ansatz) care reprezintă stări cuantice cu soluții candidate la problemă. Poți eșantiona din aceste stări candidate și apoi să le evaluezi folosind funcția de cost.

Graf → problemă de optimizare

Primul pas al mapării este o schimbare de notație. Următoarea formulă exprimă problema în notație QUBO:

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

unde $Q$ este o matrice $n\times n$ de numere reale, $n$ corespunde numărului de noduri din graful tău, $x$ este vectorul variabilelor binare introdus mai sus, iar $x^T$ indică transpusa vectorului $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

Problemă de optimizare → Hamiltonian

Poți reformula apoi problema QUBO ca un Hamiltonian (aici, o matrice care reprezintă energia unui sistem):

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

Pași de reformulare de la problema QAOA la Hamiltonian

Pentru a demonstra cum poate fi rescrisă problema QAOA în acest mod, înlocuiești mai întâi variabilele binare $x_i$ cu un nou set de variabile $z_i\in\{-1, 1\}$ prin

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

Poți vedea că dacă $x_i$ este $0$, atunci $z_i$ trebuie să fie $1$. Când $x_i$-urile sunt substituite cu $z_i$-urile în problema de optimizare ($x^TQx$), se poate obține o formulare echivalentă.

$$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}.$$

Acum dacă definim $b_i=-\sum_{j}(Q_{ij}+Q_{ji})$, eliminăm prefactorul și termenul constant $n^2$, ajungem la cele două formulări echivalente ale aceleiași probleme de optimizare.

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

Aici, $b$ depinde de $Q$. Reține că pentru a obține $z^TQz + b^Tz$ am eliminat factorul 1/4 și un offset constant de $n^2$ care nu joacă niciun rol în optimizare.

Acum, pentru a obține o formulare cuantică a problemei, promovăm variabilele $z_i$ la o matrice Pauli $Z$, cum ar fi o matrice $2\times 2$ de forma

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

Când substitui aceste matrice în problema de optimizare de mai sus, obții următorul Hamiltonian

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

Reține de asemenea că matricele $Z$ sunt integrate în spațiul computațional al calculatorului cuantic, adică un spațiu Hilbert de dimensiune $2^n\times 2^n$. Prin urmare, trebuie să înțelegi termeni precum $Z_iZ_j$ ca produsul tensorial $Z_i\otimes Z_j$ integrat în spațiul Hilbert $2^n\times 2^n$. De exemplu, într-o problemă cu cinci variabile de decizie, termenul $Z_1Z_3$ înseamnă $I\otimes Z_3\otimes I\otimes Z_1\otimes I$ unde $I$ este matricea identitate $2\times 2$.

Acest Hamiltonian se numește Hamiltonianul funcției de cost. Are proprietatea că starea sa fundamentală corespunde soluției care minimizează funcția de cost $f(x)$. Prin urmare, pentru a rezolva problema de optimizare, trebuie acum să pregătești starea fundamentală a lui $H_C$ (sau o stare cu un overlap ridicat cu aceasta) pe calculatorul cuantic. Apoi, eșantionând din această stare va produce, cu o probabilitate ridicată, soluția pentru $\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

Hamiltonian → Circuit cuantic

Hamiltonianul $H_C$ conține definiția cuantică a problemei tale. Acum poți crea un Circuit cuantic care va ajuta la eșantionarea unor soluții bune din calculatorul cuantic. QAOA este inspirat de recoacerea cuantică și aplică straturi alternante de operatori în Circuit-ul cuantic.

Ideea generală este să pornești din starea fundamentală a unui sistem cunoscut, $H^{\otimes n}|0\rangle$ de mai sus, și apoi să directionezi sistemul spre starea fundamentală a operatorului de cost care te interesează. Aceasta se realizează aplicând operatorii $\exp\{-i\gamma_k H_C\}$ și $\exp\{-i\beta_k H_m\}$ cu unghiurile $\gamma_1,...,\gamma_p$ și $\beta_1,...,\beta_p~$.

Circuit-ul cuantic pe care îl generezi este parametrizat de $\gamma_i$ și $\beta_i$, astfel că poți testa diferite valori ale lui $\gamma_i$ și $\beta_i$ și să eșantionezi din starea rezultată. În acest caz vom încerca un exemplu cu 1 strat QAOA care conține doi parametri: $\gamma_1$ și $\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 Pasul 2. Optimizarea Circuit-urilor pentru execuția pe hardware cuantic

Circuit-ul de mai sus conține o serie de abstracții utile pentru a te gândi la algoritmi cuantici, dar care nu pot fi rulate pe hardware. Pentru a putea rula pe un QPU, Circuit-ul trebuie să treacă printr-o serie de operații care alcătuiesc pasul de transpilare sau optimizare a Circuit-ului din tipar.

Biblioteca Qiskit oferă o serie de pase de transpilare care acoperă o gamă largă de transformări ale Circuit-urilor. Trebuie să te asiguri că Circuit-ul tău este optimizat pentru scopul tău.

Transpilarea poate implica mai mulți pași, cum ar fi:

• Maparea inițială a qubiților din Circuit (cum ar fi variabilele de decizie) pe qubiți fizici de pe dispozitiv.
• Desfășurarea instrucțiunilor din Circuit-ul cuantic în instrucțiunile native hardware pe care le înțelege Backend-ul.
• Rutarea oricăror qubiți din Circuit care interacționează la qubiți fizici adiacenți unii altora.
• Suprimarea erorilor prin adăugarea de Gate-uri pe un singur qubit pentru a suprima zgomotul prin decuplare dinamică.

Mai multe informații despre transpilare sunt disponibile în documentația noastră.

Codul următor transformă și optimizează Circuit-ul abstract într-un format gata de executat pe unul dintre dispozitivele accesibile prin cloud folosind serviciul Qiskit IBM® Runtime.

Reține că poți testa programele local prin „modul de testare locală" înainte de a le trimite la calculatoare cuantice reale. Mai multe informații despre modul de testare locală sunt disponibile în documentație.

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 Pasul 3. Execuție folosind primitivele Qiskit

În fluxul de lucru QAOA, parametrii optimi QAOA sunt găsiți într-o buclă de optimizare iterativă, care rulează o serie de evaluări ale Circuit-ului și folosește un optimizator clasic pentru a găsi parametrii optimi $\beta_k$ și $\gamma_k$. Această buclă de execuție este realizată prin următorii pași:

1. Definirea parametrilor inițiali
2. Instanțierea unei noi Session care conține bucla de optimizare și primitiva folosită pentru a eșantiona Circuit-ul
3. Odată ce un set optim de parametri este găsit, se execută Circuit-ul o dată în plus pentru a obține o distribuție finală care va fi folosită în pasul de post-procesare.

Definirea Circuit-ului cu parametri inițiali

Începem cu parametri aleși arbitrar.

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

Definirea Backend-ului și primitivei de execuție

Folosește primitivele Qiskit Runtime pentru a interacționa cu Backend-urile IBM®. Cele două primitive sunt Sampler și Estimator, iar alegerea primitivei depinde de tipul de măsurătoare pe care vrei să o rulezi pe calculatorul cuantic. Pentru minimizarea lui $H_C$, folosește Estimator deoarece măsurarea funcției de cost este pur și simplu valoarea așteptată a lui $\langle H_C \rangle$.

Rulare

Primitivele oferă o varietate de moduri de execuție pentru a programa lucrările pe dispozitivele cuantice, iar un flux de lucru QAOA rulează iterativ într-o sesiune. Poți conecta funcția de cost bazată pe Sampler la rutina de minimizare SciPy pentru a găsi parametrii optimi.

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

Optimizatorul a reușit să reducă costul și să găsească parametri mai buni pentru Circuit.

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

Odată ce ai găsit parametrii optimi pentru Circuit, poți atribui acești parametri și eșantiona distribuția finală obținută cu parametrii optimizați. Aici este momentul în care primitiva Sampler ar trebui folosită, deoarece distribuția de probabilitate a măsurătorilor de șiruri de biți corespunde tăieturii optime a grafului.

Notă: Aceasta înseamnă pregătirea unei stări cuantice $\psi$ în calculator și apoi măsurarea sa. O măsurătoare va prăbuși starea într-o singură stare de bază computațională — de exemplu, 010101110000... — care corespunde unei soluții candidate $x$ la problema noastră inițială de optimizare ($\max f(x)$ sau $\min f(x)$ în funcție de sarcină).

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 Pasul 4. Post-procesare, returnarea rezultatului în format clasic

Pasul de post-procesare interpretează outputul eșantionării pentru a returna o soluție la problema ta originală. În acest caz, te interesează șirul de biți cu cea mai mare probabilitate, deoarece acesta determină tăietura optimă. Simetriile din problemă permit patru soluții posibile, iar procesul de eșantionare va returna una dintre ele cu o probabilitate ușor mai mare, dar poți vedea în distribuția reprezentată mai jos că patru dintre șirurile de biți sunt considerabil mai probabile decât restul.

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

Vizualizarea celei mai bune tăieturi

Din șirul de biți optim, poți vizualiza această tăietură pe graful 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)

Și se calculează valoarea tăieturii. Soluția nu este optimă din cauza zgomotului (valoarea tăieturii soluției optime este 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

Acesta încheie tutorialul QAOA la scară mică. Vei afla cum să adaptezi QAOA la scară de utilitate în „Partea 2: scalează!" din tutorialul Quantum approximate optimization algorithm.

# Check Qiskit version
import qiskit

qiskit.__version__

'2.0.2'