Algoritmi quantistici: algoritmi quantistici variazionali

1. Introduzione

Questo tutorial fornisce una panoramica di un algoritmo ibrido quantistico-classico, concentrandosi in particolare sul variational quantum eigensolver (VQE) e sull'algoritmo di ottimizzazione approssimativa quantistica (QAOA). L'obiettivo primario di questi algoritmi è quello di affrontare problemi di ottimizzazione utilizzando circuiti quantistici con porte quantistiche parametrizzate.

Nonostante i progressi dell'informatica quantistica, la presenza di rumore negli attuali dispositivi quantistici rende difficile estrarre risultati significativi dai circuiti quantistici profondi. Per superare questa sfida, VQE e QAOA adottano un approccio ibrido quantistico-classico, che prevede l'esecuzione iterativa di circuiti quantistici relativamente brevi utilizzando la computazione quantistica e l'ottimizzazione dei parametri dei circuiti quantistici parametrizzati di destinazione utilizzando la computazione classica.

Il QAOA ha il potenziale per fornire soluzioni ottimali ai problemi target su scala utility, grazie all'applicazione di varie tecniche di mitigazione e soppressione degli errori. Il VQE ha molte applicazioni (come la chimica quantistica) in cui è meno scalabile. Tuttavia, sono emersi diversi approcci legati agli autovalori per completare e aumentare la VQE, tra cui la diagonalizzazione del sottospazio di Krylov e la diagonalizzazione quantistica basata sul campionamento (SQD). La comprensione del VQE è un primo passo importante per capire l'ampia gamma di algoritmi ibridi classico-quantistici che sono emersi.

Questo modulo descrive i concetti fondamentali e l'implementazione di VQE e QAOA. Ulteriori esercitazioni esploreranno argomenti e tecniche avanzate per scalare questi algoritmi.

Per eseguire questo notebook è necessario disporre della seguente libreria nel proprio ambiente. Se non l'avete ancora installato, potete farlo togliendo il commento ed eseguendo la seguente cella.

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

2. Calcolo dell'autovalore minimo di un hamiltoniano semplice

Inizieremo applicando la VQE a un caso molto semplice, per vedere come funziona. Calcoleremo l'autovalore minimo della matrice di Pauli $Z$ con VQE. Inizieremo 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 lo vediamo in forma matriciale.

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

Output:

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

È facile ottenere gli autovalori in modo classico, quindi possiamo verificare il nostro lavoro. Questo potrebbe diventare difficile man mano che si procede verso l'utilità. Qui utilizziamo numpy.

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

Output:

Eigenvalues: [-1.  1.]

Per ottenere gli autovalori utilizzando un algoritmo quantistico variazionale, costruiamo un circuito con porte 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")

Output:

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

sampler = StatevectorSampler()
estimator = StatevectorEstimator()

Possiamo calcolare i conteggi delle bitstring 0 e 1 con i valori dei parametri casuali [1, 2, 3] utilizzando 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}

Sappiamo che possiamo calcolare il valore di aspettativa di Z da $\langle Z \rangle = p_0 - p_1$ con le 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())

Output:

0.529296875

Questo circuito funzionava, ma i valori dei parametri scelti non corrispondevano a uno stato a bassissima energia (o a bassi autovalori). L'autovalore ottenuto è di gran lunga superiore al minimo. Il risultato è simile quando si utilizza lo stimatore.

Si noti che Estimator prende i circuiti quantistici senza misurazioni.

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

Output:

array(0.54030231)

Dovremo cercare tra i parametri e trovare quelli che producono l'autovalore più basso. Creiamo una funzione che riceve i valori dei parametri della forma variazionale e restituisce il valore dell'aspettativa $\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 SciPy's minimize 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

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 Esercizio fisico

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

Soluzioni dell'esercizio

Definiamo l'operatore di interesse e lo visualizziamo in forma matriciale.

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

Per ottenere gli autovalori utilizzando un algoritmo quantistico variazionale, costruiamo un circuito con porte 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")

Output:

Se vogliamo stimare il valore di aspettativa di un operatore (come $Z \otimes Z$ ), useremo Estimator. Se vogliamo esaminare 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

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

Questo circuito funzionava, ma i valori dei parametri scelti non corrispondevano a uno stato a bassissima energia (o a bassi autovalori). L'autovalore ottenuto è di gran lunga superiore al minimo. Il risultato è simile quando si utilizza lo stimatore.

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

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

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

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

Output:

{'01': 1024}

3. Ottimizzazione quantistica con i modelli Qiskit

In questa guida impareremo a conoscere i modelli Qiskit e l'ottimizzazione approssimativa quantistica. Un modello Qiskit è un insieme intuitivo e ripetibile di passaggi per l'implementazione di un flusso di lavoro di calcolo quantistico:

Applicheremo questi modelli al contesto dell 'ottimizzazione combinatoria e mostreremo come risolvere il problema del taglio massimo utilizzando l 'algoritmo di ottimizzazione approssimativa quantistica (QAOA), un metodo iterativo ibrido (quantistico-classico).

Si noti che questa parte di QAOA si basa sulla "Parte 1: QAOA su piccola scala" del tutorial sull' algoritmo di ottimizzazione approssimativa quantistica. Vedere il tutorial per imparare a ridimensionarlo.

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

In questa sezione useremo un problema di "max-cut" su piccola scala per illustrare i passaggi necessari per risolvere un problema di ottimizzazione utilizzando un computer quantistico.

Il problema del max-cut è un problema di ottimizzazione di difficile risoluzione (più precisamente, è un problema NP-difficile) che trova numerose applicazioni nel clustering, nella scienza delle reti e nella fisica statistica. Questo tutorial prende in esame un grafo costituito da nodi collegati da spigoli e mira a suddividere i nodi in due insiemi "tagliando" gli spigoli, in modo tale da massimizzare il numero di spigoli tagliati.

Per contestualizzare il problema prima di tradurlo in un algoritmo quantistico, è possibile comprendere meglio come il problema del taglio massimo si trasformi in un problema di ottimizzazione combinatoria classica considerando innanzitutto la minimizzazione di una funzione $f(x)$

dove l'ingresso $x$ è un vettore le cui componenti corrispondono a ciascun nodo di un grafo. Quindi, vincoliamo ciascuno di questi componenti a essere $0$ o $1$ (che rappresentano l'inclusione o l'esclusione dal taglio). Questo esempio su piccola scala utilizza un grafo con $n=5$ nodi.

Si potrebbe scrivere una funzione di una coppia di nodi $i,j$ che indichi se il bordo corrispondente $(i,j)$ è nel taglio. Per esempio, la funzione $x_i + x_j - 2 x_i x_j$ è 1 solo se uno dei due $x_i$ o $x_j$ sono 1 (il che significa che il bordo è nel taglio) e zero altrimenti. Il problema della massimizzazione degli spigoli nel taglio può essere formulato come

che può essere riscritta come una minimizzazione della forma

Il minimo di $f(x)$ in questo caso si ha quando il numero di bordi attraversati dal taglio è massimo. Come si può vedere, non c'è ancora nulla che riguardi l'informatica quantistica. È necessario riformulare questo problema in qualcosa che un computer quantistico possa capire.

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

Output:

3.2 Passaggio 1. Mappare gli input classici su un problema quantistico

Il primo passo del modello consiste nel mappare il problema classico (grafo) in circuiti e operatori quantistici. A tal fine, sono tre le fasi principali da seguire:

1. Utilizzare una serie di riformulazioni matematiche per rappresentare questo problema utilizzando la notazione dei problemi di ottimizzazione binaria non vincolata quadratica (QUBO).
2. Riscrivere il problema di ottimizzazione come un'hamiltoniana per la quale lo stato fondamentale corrisponde alla soluzione che minimizza la funzione di costo.
3. Creare un circuito quantistico che prepari lo stato fondamentale di questa hamiltoniana attraverso un processo simile alla ricottura quantistica.

Nota: nella metodologia QAOA, in ultima analisi, si vuole avere un operatore (hamiltoniano ) che rappresenti la funzione di costo del nostro 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 utilizzando la funzione di costo.

Grafico → problema di ottimizzazione

Il primo passo della mappatura è un cambiamento di notazione, di seguito viene espresso il problema in notazione QUBO:

dove $Q$ è una matrice $n\times n$ di numeri reali, $n$ corrisponde al numero di nodi del grafo, $x$ è il vettore di variabili binarie introdotto sopra e $x^T$ indica la trasposizione 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' hamiltoniana (in questo caso, una matrice che rappresenta l'energia di un sistema):

Fasi di riformulazione dal problema QAOA all'Hamiltoniana

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

Qui si può notare che se $x_i$ è $0$, allora $z_i$ deve essere $1$. Sostituendo le $x_i$ con le $z_i$ nel problema di ottimizzazione ( $x^TQx$ ), si ottiene una formulazione equivalente.

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

Qui, $b$ dipende da $Q$. Si noti che per ottenere $z^TQz + b^Tz$ abbiamo eliminato il fattore 1/4 e un offset costante di $n^2$ che non hanno alcun 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

Sostituendo queste matrici nel problema di ottimizzazione di cui sopra, si ottiene la seguente hamiltoniana

Ricordiamo inoltre che le matrici $Z$ sono incorporate nello spazio computazionale del computer quantistico, cioè uno spazio di Hilbert di dimensioni $2^n\times 2^n$. Pertanto, si devono intendere termini come $Z_iZ_j$ come il prodotto tensoriale $Z_i\otimes Z_j$ incorporato nello spazio di Hilbert $2^n\times 2^n$. Ad esempio, in un problema con cinque variabili decisionali, il termine $Z_1Z_3$ è inteso come $I\otimes Z_3\otimes I\otimes Z_1\otimes I$ dove $I$ è la matrice identità $2\times 2$.

Questa hamiltoniana è chiamata funzione di costo hamiltoniana. Ha la proprietà che il suo stato fondamentale corrisponde alla soluzione di minimizza la funzione di costo $f(x)$. Pertanto, per risolvere il problema di ottimizzazione è necessario preparare lo stato fondamentale di $H_C$ (o uno stato con un'elevata sovrapposizione con esso) sul computer quantistico. Quindi, il campionamento da questo stato produrrà, con un'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)

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

Circuito quantistico hamiltoniano

L'Hamiltonian $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 si ispira alla ricottura quantistica e applica strati alternati di operatori nel circuito quantistico.

L'idea generale è quella di partire dallo stato fondamentale di un sistema noto, $H^{\otimes n}|0\rangle$, e poi indirizzare il sistema verso lo stato fondamentale dell'operatore di costo a cui si è interessati. Ciò avviene applicando gli operatori $\exp\{-i\gamma_k H_C\}$ e $\exp\{-i\beta_k H_m\}$ con gli angoli $\gamma_1,...,\gamma_p$ e $\beta_1,...,\beta_p~$.

Il circuito quantistico generato è parametrizzato da $\gamma_i$ e $\beta_i$, quindi è possibile provare diversi valori di $\gamma_i$ e $\beta_i$ e campionare lo stato risultante.

In questo caso proveremo 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")

Output:

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

Output:

circuit.parameters

Output:

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

3.3 Fase 2. Ottimizzare i circuiti per l'esecuzione dell'hardware quantistico

Il circuito di cui sopra contiene una serie di astrazioni utili per pensare agli algoritmi quantistici, ma non è possibile eseguirli sull'hardware. Per poter essere eseguito su una QPU, il circuito deve essere sottoposto a una serie di operazioni che costituiscono la fase di transpilazione o ottimizzazione del circuito del modello.

La libreria Qiskit offre una serie di passaggi di trasposizione che soddisfano un'ampia gamma di trasformazioni circuitali. È necessario assicurarsi che il circuito sia ottimizzato per il proprio scopo.

La trasposizione potrebbe comportare diverse fasi, quali:

• Mappatura iniziale dei qubit del circuito (come le variabili decisionali) ai qubit fisici del dispositivo.
• Srotolamento delle istruzioni del circuito quantistico alle istruzioni native dell'hardware che il backend comprende.
• Instradamento dei qubit del circuito che interagiscono verso qubit fisici adiacenti.
• Soppressione degli errori mediante l'aggiunta di porte a singolo bit per sopprimere il rumore con disaccoppiamento dinamico.

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 il cloud utilizzando il servizio Qiskit IBM® Runtime.

Si noti che è possibile testare i 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)

Output:

  service = QiskitRuntimeService(channel="ibm_quantum_platform")

<IBMBackend('ibm_strasbourg')>

3.4 Fase 3. Eseguire utilizzando le primitive di IBM Quantum

Nel flusso di lavoro QAOA, i parametri ottimali di QAOA vengono trovati in un ciclo di ottimizzazione iterativa, che esegue una serie di valutazioni del circuito e utilizza un ottimizzatore classico per trovare i parametri ottimali $\beta_k$ e $\gamma_k$. Questo ciclo di esecuzione viene eseguito attraverso i seguenti passaggi:

1. Definire i parametri iniziali
2. Istanziare un nuovo Session contenente il loop di ottimizzazione e la primitiva utilizzata 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 nella fase di post-processing.

Definire il circuito con i parametri iniziali

Si parte da parametri scelti in modo arbitrario.

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

Definire il backend e la primitiva di esecuzione

Utilizzare le primitive " IBM Quantum " per interagire con i backend " IBM® ". Le due primitive sono Sampler ed Estimator, e la scelta della primitiva dipende dal tipo di misurazione che si desidera eseguire sul computer quantistico. Per minimizzare $H_C$, utilizzare Estimator poiché il valore della funzione di costo corrisponde semplicemente al valore atteso di $\langle H_C \rangle$.

Esegui

Le primitive offrono una varietà di modalità di esecuzione per programmare i carichi di lavoro sui dispositivi quantistici e un flusso di lavoro QAOA viene eseguito iterativamente in una sessione.

È possibile inserire la funzione di costo basata sul campionatore nella 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)

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

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

Output:

Una volta trovati i parametri ottimali per il circuito, è possibile assegnarli e campionare la distribuzione finale ottenuta con i parametri ottimizzati. È qui che si dovrebbe usare la primitiva Sampler, poiché è la distribuzione di probabilità delle misure delle stringhe di bit che corrisponde al taglio ottimale del grafo.

Nota: ciò significa preparare uno stato quantico $\psi$ nel computer e poi misurarlo. Una misurazione farà collassare lo stato in un singolo stato 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)

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 Fase 4. Post-elaborazione, restituisce il risultato in formato classico

La fase di post-elaborazione interpreta il risultato del campionamento per fornire una soluzione al problema originale. In questo caso, si è interessati alla stringa di bit con la probabilità più alta, poiché questa determina il taglio ottimale. Le simmetrie del problema consentono quattro possibili soluzioni e il processo di campionamento ne restituirà una con una probabilità leggermente superiore, ma si può notare nella distribuzione tracciata qui sotto che quattro delle stringhe di bit 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)

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:

Visualizza il taglio migliore

Dalla stringa di bit ottimale, è possibile visualizzare questo taglio sul grafico 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)

Output:

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

Output:

The value of the cut is: 5

Si conclude così l'esercitazione su piccola scala di QAOA. Imparerete come adattare QAOA a livello di utility in "Parte 2: scalare!" del tutorial sull' algoritmo di ottimizzazione approssimativa quantistica.

# Check Qiskit version
import qiskit

qiskit.__version__

Output:

'2.0.2'