Algoritmi quantistici: algoritmi quantistici variazionali

nota

Takashi Imamichi (24 maggio 2024)

Scarica il PDF della lezione originale. Nota che alcuni frammenti di codice potrebbero risultare deprecati, poiché si tratta di immagini statiche.

Il tempo QPU approssimativo per eseguire questo esperimento è di 9 minuti (testato su un processore Eagle).

(questo notebook potrebbe non completare l'esecuzione nel tempo consentito dal piano Open. Si prega di utilizzare le risorse di quantum computing con giudizio.)

1. Introduzione

Questo tutorial fornisce una panoramica di un algoritmo ibrido quantistico-classico, con particolare attenzione al variational quantum eigensolver (VQE) e al quantum approximate optimization algorithm (QAOA). L'obiettivo principale di questi algoritmi è affrontare problemi di ottimizzazione impiegando circuiti quantistici con gate quantistici parametrizzati.

Nonostante i progressi nel quantum computing, la presenza di rumore nei dispositivi quantistici attuali rende difficile estrarre risultati significativi da circuiti quantistici profondi. Per superare questa sfida, VQE e QAOA adottano un approccio ibrido quantistico-classico, che consiste nell'eseguire iterativamente circuiti quantistici relativamente brevi tramite calcolo quantistico e nell'ottimizzare i parametri dei circuiti quantistici parametrizzati target tramite calcolo classico.

QAOA ha il potenziale di fornire soluzioni ottimali ai problemi target su scala utility, grazie all'applicazione di varie tecniche di error mitigation e suppression. VQE ha molte applicazioni (come la chimica quantistica) in cui è meno scalabile. Tuttavia, sono emersi diversi approcci basati sugli autovalori per completare e potenziare VQE, tra cui la diagonalizzazione del sottospazio di Krylov e la sampling-based quantum diagonalization (SQD). Comprendere VQE è un primo passo fondamentale per comprendere l'ampia gamma di algoritmi ibridi classico-quantistici emersi.

Questo modulo descrive i concetti fondamentali e l'implementazione di VQE e QAOA. Tutorial successivi esploreranno argomenti avanzati e tecniche per scalare questi algoritmi. Per eseguire questo notebook è necessaria la seguente libreria nel tuo ambiente. Se non l'hai ancora installata, puoi farlo decommentando ed eseguendo la cella seguente.

# 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. Calcolo dell'autovalore minimo di un semplice Hamiltoniano

Inizieremo applicando VQE a un caso molto semplice, per vedere come funziona. Calcoleremo l'autovalore minimo della matrice di Pauli $Z$ con VQE. Cominciamo importando alcuni pacchetti generali.

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

Definiamo ora l'operatore di interesse e visualizziamolo in forma matriciale.

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

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

È facile ottenere gli autovalori in modo classico, così possiamo verificare il nostro lavoro. Questo potrebbe diventare difficile man mano che ci avviciniamo alla scala utility. Qui usiamo numpy.

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

Eigenvalues: [-1.  1.]

Per ottenere gli autovalori tramite un algoritmo quantistico variazionale, costruiamo un circuito con gate che accettano parametri variazionali:

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

Se vogliamo stimare il valore di aspettazione di un operatore (come $Z$), dobbiamo usare Estimator. Se vogliamo osservare gli stati del sistema, usiamo Sampler.

sampler = StatevectorSampler()
estimator = StatevectorEstimator()

Possiamo calcolare i conteggi delle stringhe di bit 0 e 1 con valori di parametro casuali [1, 2, 3] usando 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}

Sappiamo che possiamo calcolare il valore di aspettazione di Z tramite $\langle Z \rangle = p_0 - p_1$ con probabilità $\{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

Questo circuito ha funzionato, ma i valori dei parametri scelti non corrispondevano a uno stato a bassa energia (o basso autovalore). L'autovalore ottenuto è notevolmente più alto del minimo. Il risultato è simile quando si usa estimator.

Nota che Estimator accetta circuiti quantistici senza misurazioni.

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

array(0.54030231)

Dovremo cercare tra i parametri e trovare quelli che producono l'autovalore più basso. Definiamo una funzione che riceva i valori dei parametri della forma variazionale e restituisca il valore di aspettazione $\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

Applichiamo la funzione minimize di SciPy per trovare l'autovalore minimo di 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 Esercizio

Calcola l'autovalore minimo di $Z \otimes Z$ con 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

Soluzioni dell'esercizio

Definiamo l'operatore di interesse e visualizziamolo in forma matriciale.

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

Per ottenere gli autovalori tramite un algoritmo quantistico variazionale, costruiamo un circuito con gate che accettano parametri variazionali:

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

Se vogliamo stimare il valore di aspettazione di un operatore (come $Z \otimes Z$), useremmo Estimator. Se vogliamo osservare gli stati del sistema, usiamo 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

Questo circuito ha funzionato, ma i valori dei parametri scelti non corrispondevano a uno stato a bassa energia (o basso autovalore). L'autovalore ottenuto è notevolmente più alto del minimo. Il risultato è simile quando si usa 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)

Dovremo cercare tra i parametri e trovare quelli che producono l'autovalore più basso.

# 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

Abbiamo ottenuto un autovalore estremamente vicino al minimo fornito da 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. Ottimizzazione quantistica con i pattern Qiskit

In questo how-to impareremo i pattern Qiskit e l'ottimizzazione approssimata quantistica. Un pattern Qiskit è un insieme intuitivo e ripetibile di passi per implementare un flusso di lavoro nel quantum computing: Applicheremo i pattern al contesto dell'ottimizzazione combinatoria e mostreremo come risolvere il problema Max-Cut usando il Quantum Approximate Optimization Algorithm (QAOA), un metodo iterativo ibrido (quantistico-classico).

Nota che questa parte su QAOA è basata su "Part 1: Small-scale QAOA" del tutorial Quantum approximate optimization algorithm. Consulta il tutorial per imparare come scalarlo.

3.1 Pattern Qiskit (su piccola scala) per l'ottimizzazione

Questa parte utilizzerà un problema Max-Cut su piccola scala per illustrare i passi necessari per risolvere un problema di ottimizzazione con un computer quantistico. Il problema Max-Cut è un problema di ottimizzazione difficile da risolvere (più precisamente, è un problema NP-hard) con diverse applicazioni nel clustering, nella scienza delle reti e nella fisica statistica. Questo tutorial considera un grafo di nodi collegati da archi e mira a partizionare i nodi in due insiemi "tagliando" gli archi, in modo da massimizzare il numero di archi tagliati. Per fornire un contesto prima di mappare questo problema su un algoritmo quantistico, puoi capire meglio come il problema Max-Cut diventi un problema di ottimizzazione combinatoria classica considerando prima la minimizzazione di una funzione $f(x)$

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

dove l'input $x$ è un vettore le cui componenti corrispondono a ciascun nodo del grafo. Poi si vincola ciascuna di queste componenti a essere $0$ o $1$ (che rappresentano l'essere inclusi o meno nel taglio). Questo esempio su piccola scala usa un grafo con $n=5$ nodi.

Puoi scrivere una funzione di una coppia di nodi $i,j$ che indica se l'arco corrispondente $(i,j)$ è nel taglio. Per esempio, la funzione $x_i + x_j - 2 x_i x_j$ vale 1 solo se uno tra $x_i$ e $x_j$ è 1 (cioè l'arco è nel taglio) e zero altrimenti. Il problema di massimizzare gli archi nel taglio può essere formulato come

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

che può essere riscritto come una minimizzazione nella forma

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

Il minimo di $f(x)$ in questo caso si ottiene quando il numero di archi attraversati dal taglio è massimo. Come si può notare, non c'è ancora nulla che riguardi il quantum computing. È necessario riformulare il problema in qualcosa che un computer quantistico possa comprendere. Inizializza il problema creando un grafo con $n=5$ nodi.

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 Passo 1. Mappare gli input classici in un problema quantistico

Il primo passo del pattern consiste nel mappare il problema classico (il grafo) in circuit e operatori quantistici. Per fare ciò, ci sono tre passi principali:

1. Utilizzare una serie di riformulazioni matematiche per rappresentare il problema con la notazione dei problemi di Ottimizzazione Binaria Quadratica Non Vincolata (QUBO).
2. Riscrivere il problema di ottimizzazione come un Hamiltoniano per cui lo stato fondamentale corrisponde alla soluzione che minimizza la funzione di costo.
3. Creare un circuito quantistico che preparerà lo stato fondamentale di questo Hamiltoniano tramite un processo simile al quantum annealing.

Nota: Nella metodologia QAOA, si vuole avere un operatore (Hamiltoniano) che rappresenti la funzione di costo dell'algoritmo ibrido, nonché un circuito parametrizzato (Ansatz) che rappresenti gli stati quantistici con le soluzioni candidate al problema. È possibile campionare questi stati candidati e poi valutarli usando la funzione di costo.

Grafo → problema di ottimizzazione

Il primo passo della mappatura è un cambio di notazione. Di seguito viene espresso il problema in notazione QUBO:

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

dove $Q$ è una matrice $n\times n$ di numeri reali, $n$ corrisponde al numero di nodi nel grafo, $x$ è il vettore di variabili binarie introdotto sopra, e $x^T$ indica la trasposta del vettore $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

Problema di ottimizzazione → Hamiltoniano

È quindi possibile riformulare il problema QUBO come un Hamiltoniano (qui, una matrice che rappresenta l'energia di un sistema):

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

Passi di riformulazione dal problema QAOA all'Hamiltoniano

Per dimostrare come il problema QAOA possa essere riscritto in questo modo, si sostituiscono prima le variabili binarie $x_i$ con un nuovo insieme di variabili $z_i\in\{-1, 1\}$ tramite

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

Qui si può vedere che se $x_i$ è $0$, allora $z_i$ deve essere $1$. Quando le $x_i$ vengono sostituite dalle $z_i$ nel problema di ottimizzazione ($x^TQx$), si ottiene una formulazione equivalente.

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

Ora, se definiamo $b_i=-\sum_{j}(Q_{ij}+Q_{ji})$, rimuoviamo il prefattore e il termine costante $n^2$, arriviamo alle due formulazioni equivalenti dello stesso problema di ottimizzazione.

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

Qui, $b$ dipende da $Q$. Nota che per ottenere $z^TQz + b^Tz$ abbiamo eliminato il fattore 1/4 e un offset costante di $n^2$ che non hanno ruolo nell'ottimizzazione.

Ora, per ottenere una formulazione quantistica del problema, si promuovono le variabili $z_i$ a una matrice di Pauli $Z$, come una matrice $2\times 2$ della forma

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

Sostituendo queste matrici nel problema di ottimizzazione sopra, si ottiene il seguente Hamiltoniano

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

Ricorda anche che le matrici $Z$ sono incorporate nello spazio computazionale del computer quantistico, ovvero uno spazio di Hilbert di dimensione $2^n\times 2^n$. Pertanto, termini come $Z_iZ_j$ vanno intesi come il prodotto tensoriale $Z_i\otimes Z_j$ incorporato nello spazio di Hilbert $2^n\times 2^n$. Per esempio, in un problema con cinque variabili decisionali, il termine $Z_1Z_3$ va inteso come $I\otimes Z_3\otimes I\otimes Z_1\otimes I$ dove $I$ è la matrice identità $2\times 2$.

Questo Hamiltoniano è chiamato Hamiltoniano della funzione di costo. Ha la proprietà che il suo stato fondamentale corrisponde alla soluzione che minimizza la funzione di costo $f(x)$. Pertanto, per risolvere il tuo problema di ottimizzazione è ora necessario preparare lo stato fondamentale di $H_C$ (o uno stato con un'alta sovrapposizione con esso) sul computer quantistico. Campionando da questo stato si otterrà, con alta probabilità, la soluzione di $\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

Hamiltoniano → circuito quantistico

L'Hamiltoniano $H_C$ contiene la definizione quantistica del problema. Ora è possibile creare un circuito quantistico che aiuterà a campionare buone soluzioni dal computer quantistico. Il QAOA è ispirato al quantum annealing e applica livelli alternati di operatori nel circuito quantistico.

L'idea generale è quella di partire dallo stato fondamentale di un sistema noto, $H^{\otimes n}|0\rangle$ sopra, e poi guidare il sistema verso lo stato fondamentale dell'operatore di costo di interesse. Questo viene fatto applicando gli operatori $\exp\{-i\gamma_k H_C\}$ e $\exp\{-i\beta_k H_m\}$ con angoli $\gamma_1,...,\gamma_p$ e $\beta_1,...,\beta_p~$.

Il circuito quantistico generato è parametrizzato da $\gamma_i$ e $\beta_i$, così si possono provare diversi valori di $\gamma_i$ e $\beta_i$ e campionare dallo stato risultante. In questo caso utilizzeremo un esempio con 1 livello QAOA che contiene due parametri: $\gamma_1$ e $\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 Passo 2. Ottimizzare i circuit per l'esecuzione su hardware quantistico

Il circuito sopra contiene una serie di astrazioni utili per ragionare sugli algoritmi quantistici, ma non eseguibili sull'hardware. Per poter essere eseguito su una QPU, il circuito deve subire una serie di operazioni che costituiscono il passo di transpilazione o ottimizzazione del circuit del pattern.

La libreria Qiskit offre una serie di pass di transpilazione che coprono un'ampia gamma di trasformazioni dei circuit. È necessario assicurarsi che il circuito sia ottimizzato per il proprio scopo.

La transpilazione può comportare diversi passi, ad esempio:

• Mappatura iniziale dei qubit nel circuito (come le variabili decisionali) ai qubit fisici sul dispositivo.
• Srotolamento (unrolling) delle istruzioni nel circuito quantistico alle istruzioni native dell'hardware che il Backend comprende.
• Routing dei qubit nel circuito che interagiscono verso qubit fisici adiacenti tra loro.
• Soppressione degli errori aggiungendo gate a qubit singolo per sopprimere il rumore con il dynamical decoupling.

Ulteriori informazioni sulla transpilazione sono disponibili nella nostra documentazione.

Il codice seguente trasforma e ottimizza il circuito astratto in un formato pronto per l'esecuzione su uno dei dispositivi accessibili tramite cloud usando il servizio Qiskit IBM® Runtime.

Nota che puoi testare i tuoi programmi localmente con la "modalità di test locale" prima di inviarli ai computer quantistici reali. Ulteriori informazioni sulla modalità di test locale sono disponibili nella documentazione.

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 Passo 3. Eseguire usando le primitive Qiskit

Nel flusso di lavoro QAOA, i parametri QAOA ottimali vengono trovati in un ciclo di ottimizzazione iterativo, che esegue una serie di valutazioni del circuito e usa un ottimizzatore classico per trovare i parametri ottimali $\beta_k$ e $\gamma_k$. Questo ciclo di esecuzione viene svolto tramite i seguenti passi:

1. Definire i parametri iniziali
2. Istanziare una nuova Session contenente il ciclo di ottimizzazione e la primitiva usata per campionare il circuito
3. Una volta trovato un insieme ottimale di parametri, eseguire il circuito un'ultima volta per ottenere una distribuzione finale che verrà utilizzata nel passo di post-elaborazione.

Definire il circuito con i parametri iniziali

Iniziamo con parametri scelti arbitrariamente.

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

Definire il Backend e la primitiva di esecuzione

Usa le primitive Qiskit Runtime per interagire con i Backend IBM®. Le due primitive sono Sampler ed Estimator, e la scelta della primitiva dipende dal tipo di misura che si vuole eseguire sul computer quantistico. Per la minimizzazione di $H_C$, si utilizza l'Estimator poiché la misura della funzione di costo è semplicemente il valore atteso di $\langle H_C \rangle$.

Esecuzione

Le primitive offrono diverse modalità di esecuzione per pianificare i carichi di lavoro sui dispositivi quantistici, e un flusso di lavoro QAOA viene eseguito iterativamente in una session. Puoi collegare la funzione di costo basata sul Sampler alla routine di minimizzazione di SciPy per trovare i parametri ottimali.

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'ottimizzatore è riuscito a ridurre il costo e trovare parametri migliori per il circuito.

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

Una volta trovati i parametri ottimali per il circuito, è possibile assegnare questi parametri e campionare la distribuzione finale ottenuta con i parametri ottimizzati. Qui è dove va utilizzata la primitiva Sampler, poiché è la distribuzione di probabilità delle misurazioni delle bitstringhe che corrisponde al taglio ottimale del grafo.

Nota: Ciò significa preparare uno stato quantistico $\psi$ nel computer e poi misurarlo. Una misura farà collassare lo stato in un singolo stato della base computazionale - per esempio, 010101110000... - che corrisponde a una soluzione candidata $x$ al nostro problema di ottimizzazione iniziale ($\max f(x)$ o $\min f(x)$ a seconda del compito).

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 Passo 4. Post-elaborazione, restituzione del risultato in formato classico

Il passo di post-elaborazione interpreta l'output del campionamento per restituire una soluzione al problema originale. In questo caso, si è interessati alla bitstringha con la probabilità più alta, poiché questa determina il taglio ottimale. Le simmetrie del problema consentono quattro possibili soluzioni, e il processo di campionamento restituirà una di esse con una probabilità leggermente più alta, ma si può vedere nella distribuzione riportata di seguito che quattro delle bitstringhe sono distintamente più probabili delle altre.

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

Visualizzare il taglio migliore

Dalla bitstringha ottimale, è possibile visualizzare questo taglio sul grafo originale.

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)

E calcolare il valore del taglio. La soluzione non è ottimale a causa del rumore (il valore del taglio della soluzione ottimale è 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

Questo conclude il tutorial QAOA su piccola scala. Imparerai come adattare QAOA su scala utility in "Part 2: scale it up!" del tutorial Quantum approximate optimization algorithm.

# Check Qiskit version
import qiskit

qiskit.__version__

'2.0.2'