Funzioni di costo

In questa lezione impareremo a valutare una funzione di costo:

• Prima conosceremo le primitive Qiskit Runtime
• Definiremo una funzione di costo $C(\vec\theta)$. Si tratta di una funzione specifica al problema che ne definisce l'obiettivo che l'ottimizzatore deve minimizzare (o massimizzare)
• Definiremo una strategia di misura con le primitive Qiskit Runtime per ottimizzare la velocità rispetto alla precisione

 

Primitive

Tutti i sistemi fisici, siano essi classici o quantistici, possono esistere in stati diversi. Ad esempio, un'auto su una strada può avere una certa massa, posizione, velocità o accelerazione che ne caratterizzano lo stato. Analogamente, anche i sistemi quantistici possono avere configurazioni o stati diversi, ma si distinguono dai sistemi classici per il modo in cui gestiamo le misure e l'evoluzione degli stati. Ciò porta a proprietà uniche come la sovrapposizione e l'entanglement, esclusive della meccanica quantistica. Proprio come possiamo descrivere lo stato di un'auto usando proprietà fisiche come velocità o accelerazione, possiamo descrivere lo stato di un sistema quantistico tramite le osservabili, che sono oggetti matematici.

In meccanica quantistica, gli stati sono rappresentati da vettori colonna complessi normalizzati, o ket ($|\psi\rangle$), e le osservabili sono operatori lineari hermitiani ($\hat{H}=\hat{H}^{\dagger}$) che agiscono sui ket. Un autovettore ($|\lambda\rangle$) di un'osservabile è noto come autostato. Misurare un'osservabile per uno dei suoi autostati ($|\lambda\rangle$) restituirà il corrispondente autovalore ($\lambda$) come lettura.

Se ti stai chiedendo come misurare un sistema quantistico e cosa puoi misurare, Qiskit offre due primitive che possono aiutare:

• Sampler: Dato uno stato quantistico $|\psi\rangle$, questa primitiva ottiene la probabilità di ciascun possibile stato della base computazionale.
• Estimator: Data un'osservabile quantistica $\hat{H}$ e uno stato $|\psi\rangle$, questa primitiva calcola il valore atteso di $\hat{H}$.

La primitiva Sampler

La primitiva Sampler calcola la probabilità di ottenere ciascun possibile stato $|k\rangle$ dalla base computazionale, dato un circuito quantistico che prepara lo stato $|\psi\rangle$. Calcola

$$p_k = |\langle k | \psi \rangle|^2 \quad \forall k \in \mathbb{Z}_2^n \equiv \{0,1,\cdots,2^n-1\},$$

dove $n$ è il numero di qubit e $k$ la rappresentazione intera di qualsiasi possibile stringa binaria di output $\{0,1\}^n$ (ovvero, interi in base $2$).

Il Sampler di Qiskit Runtime esegue il circuito più volte su un dispositivo quantistico, effettuando misurazioni a ogni esecuzione, e ricostruisce la distribuzione di probabilità dalle stringhe di bit recuperate. Più esecuzioni (o shot) vengono effettuate, più accurati saranno i risultati, ma ciò richiede più tempo e risorse quantistiche.

Tuttavia, poiché il numero di possibili output cresce esponenzialmente con il numero di qubit $n$ (ovvero $2^n$), anche il numero di shot dovrà crescere esponenzialmente per catturare una distribuzione di probabilità densa. Pertanto, Sampler è efficiente solo per distribuzioni di probabilità sparse; in cui lo stato target $|\psi\rangle$ deve poter essere espresso come combinazione lineare degli stati della base computazionale, con il numero di termini che cresce al più polinomialmente con il numero di qubit:

$$|\psi\rangle = \sum^{\text{Poly}(n)}_k w_k |k\rangle.$$

Sampler può anche essere configurato per recuperare le probabilità da una sottosezione del circuito, rappresentando un sottoinsieme dei possibili stati totali.

La primitiva Estimator

La primitiva Estimator calcola il valore atteso di un'osservabile $\hat{H}$ per uno stato quantistico $|\psi\rangle$; dove le probabilità delle osservabili possono essere espresse come $p_\lambda = |\langle\lambda|\psi\rangle|^2$, essendo $|\lambda\rangle$ gli autostati dell'osservabile $\hat{H}$. Il valore atteso è quindi definito come la media di tutti i possibili risultati $\lambda$ (ovvero gli autovalori dell'osservabile) di una misurazione dello stato $|\psi\rangle$, pesata per le probabilità corrispondenti:

$$\langle\hat{H}\rangle_\psi := \sum_\lambda p_\lambda \lambda = \langle \psi | \hat{H} | \psi \rangle$$

Tuttavia, calcolare il valore atteso di un'osservabile non è sempre possibile, poiché spesso non ne conosciamo la base degli autostati. Il Estimator di Qiskit Runtime utilizza un complesso processo algebrico per stimare il valore atteso su un dispositivo quantistico reale, scomponendo l'osservabile in una combinazione di altre osservabili di cui conosciamo la base degli autostati.

In termini più semplici, Estimator scompone qualsiasi osservabile che non sa come misurare in osservabili più semplici e misurabili chiamati operatori di Pauli. Qualsiasi operatore può essere espresso come combinazione di $4^n$ operatori di Pauli.

$$\hat{P}_k := \sigma_{k_{n-1}}\otimes \cdots \otimes \sigma_{k_0} \quad \forall k \in \mathbb{Z}_4^n \equiv \{0,1,\cdots,4^n-1\}, \\$$

tale che

$$\hat{H} = \sum^{4^n-1}_{k=0} w_k \hat{P}_k$$

dove $n$ è il numero di qubit, $k \equiv k_{n-1} \cdots k_0$ per $k_l \in \mathbb{Z}_4 \equiv \{0, 1, 2, 3\}$ (ovvero, interi in base $4$), e $(\sigma_0, \sigma_1, \sigma_2, \sigma_3) := (I, X, Y, Z)$.

Dopo aver eseguito questa decomposizione, Estimator deriva un nuovo circuito $V_k|\psi\rangle$ per ciascuna osservabile $\hat{P}_k$ (dal circuito originale), per diagonalizzare efficacemente l'osservabile di Pauli nella base computazionale e misurarla. Possiamo misurare facilmente le osservabili di Pauli perché conosciamo $V_k$ in anticipo, il che non è il caso in generale per altre osservabili.

Per ciascun $\hat{P}_{k}$, Estimator esegue il circuito corrispondente su un dispositivo quantistico più volte, misura lo stato di output nella base computazionale e calcola la probabilità $p_{kj}$ di ottenere ciascun possibile output $j$. Cerca poi l'autovalore $\lambda_{kj}$ di $P_k$ corrispondente a ciascun output $j$, moltiplica per $w_k$ e somma tutti i risultati per ottenere il valore atteso dell'osservabile $\hat{H}$ per lo stato $|\psi\rangle$ dato.

$$\langle\hat{H}\rangle_\psi = \sum_{k=0}^{4^n-1} w_k \sum_{j=0}^{2^n-1}p_{kj} \lambda_{kj},$$

Poiché calcolare il valore atteso di $4^n$ Pauli è impraticabile (ovvero, cresce esponenzialmente), Estimator può essere efficiente solo quando un gran numero di $w_k$ è zero (ovvero, decomposizione di Pauli sparsa invece che densa). Formalmente diciamo che, affinché questo calcolo sia risolvibile in modo efficiente, il numero di termini non nulli deve crescere al più polinomialmente con il numero di qubit $n$: $\hat{H} = \sum^{\text{Poly}(n)}_k w_k \hat{P}_k.$

Il lettore potrebbe notare l'assunzione implicita che anche il campionamento delle probabilità deve essere efficiente come spiegato per Sampler, il che significa

$$\langle\hat{H}\rangle_\psi = \sum_{k}^{\text{Poly}(n)} w_k \sum_{j}^{\text{Poly}(n)}p_{kj} \lambda_{kj}.$$

Esempio guidato per calcolare i valori attesi

Supponiamo lo stato a singolo qubit $|+\rangle := H|0\rangle = \frac{1}{\sqrt{2}}(|0\rangle + |1\rangle)$, e l'osservabile

$$\begin{aligned} \hat{H} & = \begin{pmatrix} -1 & 2 \\ 2 & 1 \\ \end{pmatrix}\\[1mm] & = 2X - Z \end{aligned}$$

con il seguente valore atteso teorico $\langle\hat{H}\rangle_+ = \langle+|\hat{H}|+\rangle = 2.$

Poiché non sappiamo come misurare questa osservabile, non possiamo calcolarne il valore atteso direttamente e dobbiamo riscriverla come $\langle\hat{H}\rangle_+ = 2\langle X \rangle_+ - \langle Z \rangle_+$. Si può dimostrare che questo dà lo stesso risultato notando che $\langle+|X|+\rangle = 1$ e $\langle+|Z|+\rangle = 0$.

Vediamo come calcolare direttamente $\langle X \rangle_+$ e $\langle Z \rangle_+$. Poiché $X$ e $Z$ non commutano (ovvero non condividono la stessa base degli autostati), non possono essere misurati simultaneamente, pertanto servono i circuiti ausiliari:

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

from qiskit import QuantumCircuit
from qiskit.quantum_info import SparsePauliOp

# The following code will work for any other initial single-qubit state and observable
original_circuit = QuantumCircuit(1)
original_circuit.h(0)

H = SparsePauliOp(["X", "Z"], [2, -1])

aux_circuits = []
for pauli in H.paulis:
    aux_circ = original_circuit.copy()
    aux_circ.barrier()
    if str(pauli) == "X":
        aux_circ.h(0)
    elif str(pauli) == "Y":
        aux_circ.sdg(0)
        aux_circ.h(0)
    else:
        aux_circ.id(0)
    aux_circ.measure_all()
    aux_circuits.append(aux_circ)

original_circuit.draw("mpl")

# Auxiliary circuit for X
aux_circuits[0].draw("mpl")

# Auxiliary circuit for Z
aux_circuits[1].draw("mpl")

Possiamo ora eseguire il calcolo manualmente usando Sampler e verificare i risultati con Estimator:

from qiskit.primitives import StatevectorSampler, StatevectorEstimator
from qiskit.result import QuasiDistribution
import numpy as np

## SAMPLER
shots = 10000
sampler = StatevectorSampler()
job = sampler.run(aux_circuits, shots=shots)

# Run the sampler job and step through results
expvals = []
for index, pauli in enumerate(H.paulis):
    data_pub = job.result()[index].data
    bitstrings = data_pub.meas.get_bitstrings()
    counts = data_pub.meas.get_counts()
    quasi_dist = QuasiDistribution(
        {outcome: freq / shots for outcome, freq in counts.items()}
    )

    # Use the probabilities and known eigenvalues of Pauli operators to estimate the expectation value.
    val = 0

    if str(pauli) == "X":
        val += -1 * quasi_dist.get(1, 0)
        val += 1 * quasi_dist.get(0, 0)

    if str(pauli) == "Y":
        val += -1 * quasi_dist.get(1, 0)
        val += 1 * quasi_dist.get(0, 0)

    if str(pauli) == "Z":
        val += 1 * quasi_dist.get(0, 0)
        val += -1 * quasi_dist.get(1, 0)

    expvals.append(val)

# Print expectation values

print("Sampler results:")
for pauli, expval in zip(H.paulis, expvals):
    print(f"  >> Expected value of {str(pauli)}: {expval:.5f}")

total_expval = np.sum(H.coeffs * expvals).real
print(f"  >> Total expected value: {total_expval:.5f}")

# Use estimator for comparison
observables = [
    *H.paulis,
    H,
]  # Note: run for individual Paulis as well as full observable H

estimator = StatevectorEstimator()
job = estimator.run([(original_circuit, observables)])
estimator_expvals = job.result()[0].data.evs

# Print results
print("Estimator results:")
for obs, expval in zip(observables, estimator_expvals):
    if obs is not H:
        print(f"  >> Expected value of {str(obs)}: {expval:.5f}")
    else:
        print(f"  >> Total expected value: {expval:.5f}")

Sampler results:
  >> Expected value of X: 1.00000
  >> Expected value of Z: 0.00420
  >> Total expected value: 1.99580
Estimator results:
  >> Expected value of X: 1.00000
  >> Expected value of Z: 0.00000
  >> Total expected value: 2.00000

Rigore matematico (opzionale)

Esprimendo $|\psi\rangle$ rispetto alla base degli autostati di $\hat{H}$, $|\psi\rangle = \sum_\lambda a_\lambda |\lambda\rangle$, si ottiene:

$$\begin{aligned} \langle \psi | \hat{H} | \psi \rangle & = \bigg(\sum_{\lambda'}a^*_{\lambda'} \langle \lambda'|\bigg) \hat{H} \bigg(\sum_{\lambda} a_\lambda | \lambda\rangle\bigg)\\[1mm] & = \sum_{\lambda}\sum_{\lambda'} a^*_{\lambda'}a_{\lambda} \langle \lambda'|\hat{H}| \lambda\rangle\\[1mm] & = \sum_{\lambda}\sum_{\lambda'} a^*_{\lambda'}a_{\lambda} \lambda \langle \lambda'| \lambda\rangle\\[1mm] & = \sum_{\lambda}\sum_{\lambda'} a^*_{\lambda'}a_{\lambda} \lambda \cdot \delta_{\lambda, \lambda'}\\[1mm] & = \sum_\lambda |a_\lambda|^2 \lambda\\[1mm] & = \sum_\lambda p_\lambda \lambda\\[1mm] \end{aligned}$$

Poiché non conosciamo gli autovalori o gli autostati dell'osservabile target $\hat{H}$, dobbiamo prima considerarne la diagonalizzazione. Dato che $\hat{H}$ è hermitiana, esiste una trasformazione unitaria $V$ tale che $\hat{H}=V^\dagger \Lambda V,$ dove $\Lambda$ è la matrice diagonale degli autovalori, così $\langle j | \Lambda | k \rangle = 0$ se $j\neq k$, e $\langle j | \Lambda | j \rangle = \lambda_j$.

Ciò implica che il valore atteso può essere riscritto come:

$$\begin{aligned} \langle\psi|\hat{H}|\psi\rangle & = \langle\psi|V^\dagger \Lambda V|\psi\rangle\\[1mm] & = \langle\psi|V^\dagger \bigg(\sum_{j=0}^{2^n-1} |j\rangle \langle j|\bigg) \Lambda \bigg(\sum_{k=0}^{2^n-1} |k\rangle \langle k|\bigg) V|\psi\rangle\\[1mm] & = \sum_{j=0}^{2^n-1} \sum_{k=0}^{2^n-1}\langle\psi|V^\dagger |j\rangle \langle j| \Lambda |k\rangle \langle k| V|\psi\rangle\\[1mm] & = \sum_{j=0}^{2^n-1}\langle\psi|V^\dagger |j\rangle \langle j| \Lambda |j\rangle \langle j| V|\psi\rangle\\[1mm] & = \sum_{j=0}^{2^n-1}|\langle j| V|\psi\rangle|^2 \lambda_j\\[1mm] \end{aligned}$$

Dato che se un sistema si trova nello stato $|\phi\rangle = V |\psi\rangle$ la probabilità di misurare $| j\rangle$ è $p_j = |\langle j|\phi \rangle|^2$, il valore atteso sopra può essere espresso come:

$$\langle\psi|\hat{H}|\psi\rangle = \sum_{j=0}^{2^n-1} p_j \lambda_j.$$

È molto importante notare che le probabilità sono calcolate dallo stato $V |\psi\rangle$ invece che da $|\psi\rangle$. Ecco perché la matrice $V$ è assolutamente necessaria. Potresti chiederti come ottenere la matrice $V$ e gli autovalori $\Lambda$. Se avessi già gli autovalori, non ci sarebbe bisogno di usare un computer quantistico, poiché l'obiettivo degli algoritmi variazionali è proprio trovare questi autovalori di $\hat{H}$.

Fortunatamente, esiste un modo alternativo: qualsiasi matrice $2^n \times 2^n$ può essere scritta come combinazione lineare di $4^n$ prodotti tensoriali di $n$ matrici di Pauli e identità, tutte hermitiane e unitarie con $V$ e $\Lambda$ noti. Questo è ciò che fa internamente l'Estimator di Runtime, decomponendo qualsiasi oggetto Operator in un SparsePauliOp.

Ecco gli operatori che possono essere utilizzati:

$$\begin{array}{c|c|c|c} \text{Operator} & \sigma & V & \Lambda \\[1mm] \hline I & \sigma_0 = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} & V_0 = I & \Lambda_0 = I = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} \\[4mm] X & \sigma_1 = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} & V_1 = H =\frac{1}{\sqrt{2}} \begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix} & \Lambda_1 = \sigma_3 = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} \\[4mm] Y & \sigma_2 = \begin{pmatrix} 0 & -i \\ i & 0 \end{pmatrix} & V_2 = HS^\dagger =\frac{1}{\sqrt{2}} \begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix}\cdot \begin{pmatrix} 1 & 0 \\ 0 & -i \end{pmatrix} = \frac{1}{\sqrt{2}} \begin{pmatrix} 1 & -i \\ 1 & i \end{pmatrix}\quad & \Lambda_2 = \sigma_3 = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} \\[4mm] Z & \sigma_3 = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} & V_3 = I & \Lambda_3 = \sigma_3 = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} \end{array}$$

Quindi riscriviamo $\hat{H}$ rispetto alle matrici di Pauli e alle identità:

$$\hat{H} = \sum_{k_{n-1}=0}^3... \sum_{k_0=0}^3 w_{k_{n-1}...k_0} \sigma_{k_{n-1}}\otimes ... \otimes \sigma_{k_0} = \sum_{k=0}^{4^n-1} w_k \hat{P}_k,$$

dove $k = \sum_{l=0}^{n-1} 4^l k_l \equiv k_{n-1}...k_0$ per $k_{n-1},...,k_0\in \{0,1,2,3\}$ (ovvero, base $4$), e $\hat{P}_{k} := \sigma_{k_{n-1}}\otimes ... \otimes \sigma_{k_0}$:

$$\begin{aligned} \langle\psi|\hat{H}|\psi\rangle & = \sum_{k=0}^{4^n-1} w_k \sum_{j=0}^{2^n-1}|\langle j| V_k|\psi\rangle|^2 \langle j| \Lambda_k |j\rangle \\[1mm] & = \sum_{k=0}^{4^n-1} w_k \sum_{j=0}^{2^n-1}p_{kj} \lambda_{kj}, \\[1mm] \end{aligned}$$

dove $V_k := V_{k_{n-1}}\otimes ... \otimes V_{k_0}$ e $\Lambda_k := \Lambda_{k_{n-1}}\otimes ... \otimes \Lambda_{k_0}$, tale che: $\hat{P_k}=V_k^\dagger \Lambda_k V_k.$

Funzioni di costo

In generale, le funzioni di costo servono a descrivere l'obiettivo di un problema e quanto bene uno stato candidato si avvicina a tale obiettivo. Questa definizione si applica a vari esempi in chimica, machine learning, finanza, ottimizzazione e così via.

Consideriamo un semplice esempio: trovare lo stato fondamentale di un sistema. Il nostro obiettivo è minimizzare il valore atteso dell'osservabile che rappresenta l'energia (l'Hamiltoniano $\hat{\mathcal{H}}$):

$$\min_{\vec\theta} \langle\psi(\vec\theta)|\hat{\mathcal{H}}|\psi(\vec\theta)\rangle$$

Possiamo usare l'Estimator per calcolare il valore atteso e passare questo valore a un ottimizzatore affinché lo minimizzi. Se l'ottimizzazione ha successo, restituirà un insieme di valori ottimali dei parametri $\vec\theta^*$, dai quali saremo in grado di costruire lo stato soluzione proposto $|\psi(\vec\theta^*)\rangle$ e calcolare il valore atteso osservato come $C(\vec\theta^*)$.

Nota che saremo in grado di minimizzare la funzione di costo soltanto per l'insieme limitato di stati che stiamo considerando. Ciò ci porta a due possibilità distinte:

• Il nostro ansatz non definisce lo stato soluzione nello spazio di ricerca: In questo caso, l'ottimizzatore non troverà mai la soluzione e dobbiamo sperimentare altri ansatz che potrebbero rappresentare lo spazio di ricerca in modo più accurato.
• L'ottimizzatore non riesce a trovare questa soluzione valida: L'ottimizzazione può essere definita globalmente o localmente. Esploreremo il significato di questo nella sezione successiva.

In sostanza, eseguiremo un ciclo di ottimizzazione classica affidandoci però alla valutazione della funzione di costo su un computer quantistico. Da questa prospettiva, si potrebbe pensare all'ottimizzazione come a un processo puramente classico in cui chiamiamo un qualche oracolo quantistico a scatola nera ogni volta che l'ottimizzatore deve valutare la funzione di costo.

def cost_func_vqe(params, circuit, hamiltonian, estimator):
    """Return estimate of energy from estimator

    Parameters:
        params (ndarray): Array of ansatz parameters
        ansatz (QuantumCircuit): Parameterized ansatz circuit
        hamiltonian (SparsePauliOp): Operator representation of Hamiltonian
        estimator (Estimator): Estimator primitive instance

    Returns:
        float: Energy estimate
    """
    pub = (circuit, hamiltonian, params)
    cost = estimator.run([pub]).result()[0].data.evs
    return cost

from qiskit.circuit.library import TwoLocal

observable = SparsePauliOp.from_list([("XX", 1), ("YY", -3)])

reference_circuit = QuantumCircuit(2)
reference_circuit.x(0)

variational_form = TwoLocal(
    2,
    rotation_blocks=["rz", "ry"],
    entanglement_blocks="cx",
    entanglement="linear",
    reps=1,
)
ansatz = reference_circuit.compose(variational_form)

theta_list = (2 * np.pi * np.random.rand(1, 8)).tolist()
ansatz.decompose().draw("mpl")

Eseguiremo prima questo calcolo su un simulatore: lo StatevectorEstimator. Di solito è consigliabile per il debug, ma faremo seguire subito al run di debug un calcolo su hardware quantistico reale. Sempre più spesso, i problemi di interesse non sono più simulabili classicamente senza strutture di supercalcolo all'avanguardia.

estimator = StatevectorEstimator()
cost = cost_func_vqe(theta_list, ansatz, observable, estimator)
print(cost)

[-0.58744589]

Procederemo ora con l'esecuzione su un vero computer quantistico. Nota le differenze nella sintassi. I passaggi che coinvolgono il pass_manager saranno approfonditi nell'esempio successivo. Un passaggio di particolare importanza negli algoritmi variazionali è l'uso di una sessione Qiskit Runtime. Aprire una sessione consente di eseguire più iterazioni di un algoritmo variazionale senza dover attendere in una nuova coda ogni volta che i parametri vengono aggiornati. Questo è importante quando i tempi di coda sono lunghi e/o sono necessarie molte iterazioni. Solo i partner nella IBM Quantum® Network possono utilizzare le sessioni Runtime. Se non hai accesso alle sessioni, puoi ridurre il numero di iterazioni da inviare in un dato momento e salvare i parametri più recenti per usi futuri. Se invii troppe iterazioni o incontri tempi di coda eccessivamente lunghi, potresti ricevere il codice di errore 1217, che si riferisce a lunghi ritardi tra le sottomissioni di job.

# Estimated usage: < 1 min. Benchmarked at 7 seconds on an Eagle processor
# Load necessary packages:

from qiskit_ibm_runtime import (
    QiskitRuntimeService,
    Session,
    EstimatorOptions,
    EstimatorV2 as Estimator,
)
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

# Select the least busy backend:

service = QiskitRuntimeService()
backend = service.least_busy(
    operational=True, min_num_qubits=ansatz.num_qubits, simulator=False
)
# Or get a specific backend:
# backend = service.backend("ibm_brisbane")

# Use a pass manager to transpile the circuit and observable for the specific backend being used:

pm = generate_preset_pass_manager(backend=backend, optimization_level=1)
isa_ansatz = pm.run(ansatz)
isa_observable = observable.apply_layout(layout=isa_ansatz.layout)

# Set estimator options
estimator_options = EstimatorOptions(resilience_level=1, default_shots=10_000)

# Open a Runtime session:

with Session(backend=backend) as session:
    estimator = Estimator(mode=session, options=estimator_options)
    cost = cost_func_vqe(theta_list, isa_ansatz, isa_observable, estimator)

session.close()
print(cost)

Nota che i valori ottenuti dai due calcoli precedenti sono molto simili. Le tecniche per migliorare i risultati saranno discusse più avanti.

Esempio di mappatura su sistemi non fisici

Il problema del taglio massimo (Max-Cut) è un problema di ottimizzazione combinatoria che consiste nel dividere i vertici di un grafo in due insiemi disgiunti in modo da massimizzare il numero di archi tra i due insiemi. Più formalmente, dato un grafo non orientato $G=(V,E)$, dove $V$ è l'insieme dei vertici ed $E$ è l'insieme degli archi, il problema Max-Cut chiede di partizionare i vertici in due sottoinsiemi disgiunti, $S$ e $T$, in modo che il numero di archi con un estremo in $S$ e l'altro in $T$ sia massimizzato.

Possiamo applicare Max-Cut per risolvere vari problemi tra cui: clustering, progettazione di reti, transizioni di fase e così via. Inizieremo creando un grafo del problema:

import rustworkx as rx
from rustworkx.visualization import mpl_draw

n = 4
G = rx.PyGraph()
G.add_nodes_from(range(n))
# The edge syntax is (start, end, weight)
edges = [(0, 1, 1.0), (0, 2, 1.0), (0, 3, 1.0), (1, 2, 1.0), (2, 3, 1.0)]
G.add_edges_from(edges)

mpl_draw(
    G, pos=rx.shell_layout(G), with_labels=True, edge_labels=str, node_color="#1192E8"
)

Questo problema può essere espresso come un problema di ottimizzazione binaria. Per ogni nodo $0 \leq i < n$, dove $n$ è il numero di nodi del grafo (in questo caso $n=4$), considereremo la variabile binaria $x_i$. Questa variabile avrà valore $1$ se il nodo $i$ appartiene a uno dei gruppi che chiameremo $1$, e $0$ se appartiene all'altro gruppo, che chiameremo $0$. Indicheremo inoltre con $w_{ij}$ (elemento $(i,j)$ della matrice di adiacenza $w$) il peso dell'arco che va dal nodo $i$ al nodo $j$. Poiché il grafo non è orientato, $w_{ij}=w_{ji}$. Possiamo quindi formulare il nostro problema come la massimizzazione della seguente funzione di costo:

$$\begin{aligned} C(\vec{x}) & =\sum_{i,j=0}^n w_{ij} x_i(1-x_j)\\[1mm] & = \sum_{i,j=0}^n w_{ij} x_i - \sum_{i,j=0}^n w_{ij} x_ix_j\\[1mm] & = \sum_{i,j=0}^n w_{ij} x_i - \sum_{i=0}^n \sum_{j=0}^i 2w_{ij} x_ix_j \end{aligned}$$

Per risolvere questo problema con un computer quantistico, dobbiamo esprimere la funzione di costo come valore atteso di un osservabile. Tuttavia, gli osservabili che Qiskit ammette nativamente sono composti da operatori di Pauli, che hanno autovalori $1$ e $-1$ anziché $0$ e $1$. Per questo motivo effettueremo il seguente cambio di variabile:

Dove $\vec{x}=(x_0,x_1,\cdots ,x_{n-1})$. Possiamo usare la matrice di adiacenza $w$ per accedere comodamente ai pesi di tutti gli archi. Questo verrà usato per ottenere la nostra funzione di costo:

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

Questo implica che:

$$\begin{array}{lcl} x_i=0 & \rightarrow & z_i=1 \\ x_i=1 & \rightarrow & z_i=-1.\end{array}$$

Quindi la nuova funzione di costo che vogliamo massimizzare è:

$$\begin{aligned} C(\vec{z}) & = \sum_{i,j=0}^n w_{ij} \bigg(\frac{1-z_i}{2}\bigg)\bigg(1-\frac{1-z_j}{2}\bigg)\\[1mm] & = \sum_{i,j=0}^n \frac{w_{ij}}{4} - \sum_{i,j=0}^n \frac{w_{ij}}{4} z_iz_j\\[1mm] & = \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2} - \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2} z_iz_j \end{aligned}$$

Inoltre, la tendenza naturale di un computer quantistico è quella di trovare minimi (di solito l'energia più bassa) anziché massimi, quindi invece di massimizzare $C(\vec{z})$ minimizzeremo:

$$-C(\vec{z}) = \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2} z_iz_j - \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2}$$

Ora che abbiamo una funzione di costo da minimizzare le cui variabili possono assumere i valori $-1$ e $1$, possiamo fare la seguente analogia con il Pauli $Z$:

$$z_i \equiv Z_i = \overbrace{I}^{n-1}\otimes ... \otimes \overbrace{Z}^{i} \otimes ... \otimes \overbrace{I}^{0}$$

In altre parole, la variabile $z_i$ sarà equivalente a un gate $Z$ che agisce sul qubit $i$. Inoltre:

$$Z_i|x_{n-1}\cdots x_0\rangle = z_i|x_{n-1}\cdots x_0\rangle \rightarrow \langle x_{n-1}\cdots x_0 |Z_i|x_{n-1}\cdots x_0\rangle = z_i$$

L'osservabile che considereremo è quindi:

$$\hat{H} = \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2} Z_iZ_j$$

al quale dovremo aggiungere successivamente il termine indipendente:

$$\texttt{offset} = - \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2}$$

L'operatore è una combinazione lineare di termini con operatori Z sui nodi collegati da un arco (si ricorda che il qubit 0 è il più a destra): $IIZZ + IZIZ + IZZI + ZIIZ + ZZII$. Una volta costruito l'operatore, l'ansatz per l'algoritmo QAOA può essere facilmente costruito usando il circuito QAOAAnsatz dalla libreria di circuiti Qiskit.

from qiskit.circuit.library import QAOAAnsatz
from qiskit.quantum_info import SparsePauliOp

hamiltonian = SparsePauliOp.from_list(
    [("IIZZ", 1), ("IZIZ", 1), ("IZZI", 1), ("ZIIZ", 1), ("ZZII", 1)]
)

ansatz = QAOAAnsatz(hamiltonian, reps=2)
# Draw
ansatz.decompose(reps=3).draw("mpl")

# Sum the weights, and divide by 2

offset = -sum(edge[2] for edge in edges) / 2
print(f"""Offset: {offset}""")

Offset: -2.5

Con l'Estimator Runtime che riceve direttamente un Hamiltoniano e un ansatz parametrizzato restituendo l'energia necessaria, la funzione di costo per un'istanza QAOA è piuttosto semplice:

def cost_func(params, ansatz, hamiltonian, estimator):
    """Return estimate of energy from estimator

    Parameters:
        params (ndarray): Array of ansatz parameters
        ansatz (QuantumCircuit): Parameterized ansatz circuit
        hamiltonian (SparsePauliOp): Operator representation of Hamiltonian
        estimator (Estimator): Estimator primitive instance

    Returns:
        float: Energy estimate
    """
    pub = (ansatz, hamiltonian, params)
    cost = estimator.run([pub]).result()[0].data.evs
    #    cost = estimator.run(ansatz, hamiltonian, parameter_values=params).result().values[0]
    return cost

import numpy as np

x0 = 2 * np.pi * np.random.rand(ansatz.num_parameters)

estimator = StatevectorEstimator()
cost = cost_func_vqe(x0, ansatz, hamiltonian, estimator)
print(cost)

1.473098768180865

# Estimated usage: < 1 min, benchmarked at 6 seconds on ibm_osaka, 5-23-24
# Load some necessary packages:

from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_ibm_runtime import Session, EstimatorV2 as Estimator

# Select the least busy backend:

backend = service.least_busy(
    operational=True, min_num_qubits=ansatz.num_qubits, simulator=False
)

# Or get a specific backend:
# backend = service.backend("ibm_brisbane")

# Use a pass manager to transpile the circuit and observable for the specific backend being used:

pm = generate_preset_pass_manager(backend=backend, optimization_level=1)
isa_ansatz = pm.run(ansatz)
isa_hamiltonian = hamiltonian.apply_layout(layout=isa_ansatz.layout)

# Set estimator options
estimator_options = EstimatorOptions(resilience_level=1, default_shots=10_000)

# Open a Runtime session:

with Session(backend=backend) as session:
    estimator = Estimator(mode=session, options=estimator_options)
    cost = cost_func_vqe(x0, isa_ansatz, isa_hamiltonian, estimator)

# Close session after done
session.close()
print(cost)

1.1120776913677988

Ritorneremo su questo esempio nella sezione Applicazioni per esplorare come sfruttare un ottimizzatore per iterare nello spazio di ricerca. In linea generale, questo include:

• Sfruttare un ottimizzatore per trovare i parametri ottimali
• Associare i parametri ottimali all'ansatz per trovare gli autovalori
• Tradurre gli autovalori nella nostra definizione del problema

Strategia di misura: velocità versus precisione

Come accennato, stiamo usando un computer quantistico rumoroso come oracolo a scatola nera, dove il rumore può rendere i valori recuperati non deterministici, portando a fluttuazioni casuali che, a loro volta, ostacoleranno — o addirittura impediranno completamente — la convergenza di certi ottimizzatori verso una soluzione proposta. Questo è un problema generale che dobbiamo affrontare man mano che esploriamo in modo incrementale l'utilità quantistica e progrediamo verso il vantaggio quantistico:

Possiamo usare le opzioni di soppressione e mitigazione degli errori delle Primitive Qiskit Runtime per affrontare il rumore e massimizzare l'utilità dei computer quantistici odierni.

Soppressione degli errori

La soppressione degli errori si riferisce alle tecniche usate per ottimizzare e trasformare un circuito durante la compilazione al fine di minimizzare gli errori. È una tecnica di base per la gestione degli errori che di solito comporta un certo overhead di pre-elaborazione classica sul tempo di esecuzione complessivo. L'overhead include la trascompilazione dei circuiti per farli girare su hardware quantistico tramite:

• Esprimere il circuito usando i gate nativi disponibili su un sistema quantistico
• Mappare i qubit virtuali sui qubit fisici
• Aggiungere SWAP in base ai requisiti di connettività
• Ottimizzare i gate a 1 e 2 qubit
• Aggiungere il decoupling dinamico ai qubit inattivi per prevenire gli effetti della decoerenza.

Le Primitive consentono di usare tecniche di soppressione degli errori impostando l'opzione optimization_level e selezionando opzioni di trascompilazione avanzate. In un corso successivo, approfondiremo diversi metodi di costruzione dei circuiti per migliorare i risultati, ma per la maggior parte dei casi raccomandiamo di impostare optimization_level=3.

Visualizzeremo il valore dell'aumento dell'ottimizzazione nel processo di trascompilazione esaminando un circuito di esempio con un comportamento ideale semplice.

from qiskit.circuit import Parameter, QuantumCircuit
from qiskit.quantum_info import SparsePauliOp

theta = Parameter("theta")

qc = QuantumCircuit(2)
qc.x(1)
qc.h(0)
qc.cp(theta, 0, 1)
qc.h(0)
observables = SparsePauliOp.from_list([("ZZ", 1)])

qc.draw("mpl")

Il circuito qui sopra può produrre valori attesi sinusoidali dell'osservabile dato, a condizione di inserire fasi che coprano un intervallo appropriato, come $[0,2\pi]$.

## Setup phases
import numpy as np

phases = np.linspace(0, 2 * np.pi, 50)

# phases need to be expressed as a list of lists in order to work
individual_phases = [[phase] for phase in phases]

Possiamo usare un simulatore per mostrare l'utilità di una trascompilazione ottimizzata. Torneremo più avanti all'uso di hardware reale per dimostrare l'utilità della mitigazione degli errori. Useremo QiskitRuntimeService per ottenere un backend reale (in questo caso, ibm_brisbane) e AerSimulator per simulare quel backend, incluso il suo comportamento rumoroso.

from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_aer import AerSimulator

# get a real backend from the runtime service
service = QiskitRuntimeService()
backend = service.backend("ibm_brisbane")

# generate a simulator that mimics the real quantum system with the latest calibration results
backend_sim = AerSimulator.from_backend(backend)

Possiamo ora usare un pass manager per trascompilare il circuito nell'"architettura del set di istruzioni" (ISA) del backend. Questo è un nuovo requisito in Qiskit Runtime: tutti i circuiti inviati a un backend devono essere conformi ai vincoli del target del backend, il che significa che devono essere scritti in termini dell'ISA del backend — ovvero l'insieme di istruzioni che il dispositivo è in grado di comprendere ed eseguire. Questi vincoli del target sono definiti da fattori come i gate base nativi del dispositivo, la connettività dei suoi qubit e — quando rilevante — le specifiche di temporizzazione degli impulsi e delle altre istruzioni.

Nota che nel caso presente lo faremo due volte: una con optimization_level = 0 e una con il valore impostato a 3. Ogni volta useremo la primitiva Estimator per stimare i valori attesi dell'osservabile a diversi valori di fase.

# Import estimator and specify that we are using the simulated backend:

from qiskit_ibm_runtime import EstimatorV2 as Estimator

estimator = Estimator(mode=backend_sim)

circuit = qc

# Use a pass manager to transpile the circuit and observable for the backend being simulated.
# Start with no optimization:

from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

pm = generate_preset_pass_manager(backend=backend_sim, optimization_level=0)
isa_circuit = pm.run(circuit)
isa_observables = observables.apply_layout(layout=isa_circuit.layout)

noisy_exp_values = []
pub = (isa_circuit, isa_observables, [individual_phases])
cost = estimator.run([pub]).result()[0].data.evs
noisy_exp_values = cost[0]

# Repeat above steps, but now with optimization = 3:

exp_values_with_opt_es = []
pm = generate_preset_pass_manager(backend=backend_sim, optimization_level=3)
isa_circuit = pm.run(circuit)
isa_observables = observables.apply_layout(layout=isa_circuit.layout)

pub = (isa_circuit, isa_observables, [individual_phases])
cost = estimator.run([pub]).result()[0].data.evs
exp_values_with_opt_es = cost[0]

Infine, possiamo tracciare i risultati e osserviamo che la precisione del calcolo era abbastanza buona anche senza ottimizzazione, ma è definitivamente migliorata aumentando l'ottimizzazione al livello 3. Nota che in circuiti più profondi e complessi, la differenza tra i livelli di ottimizzazione 0 e 3 è probabilmente più significativa. Questo è un circuito molto semplice usato come modello giocattolo.

import matplotlib.pyplot as plt

plt.plot(phases, noisy_exp_values, "o", label="opt=0")
plt.plot(phases, exp_values_with_opt_es, "o", label="opt=3")
plt.plot(phases, 2 * np.sin(phases / 2) ** 2 - 1, label="ideal")
plt.ylabel("Expectation")
plt.legend()
plt.show()

Mitigazione degli errori

La mitigazione degli errori si riferisce alle tecniche che consentono agli utenti di ridurre gli errori dei circuiti modellando il rumore del dispositivo al momento dell'esecuzione. In genere, questo comporta un overhead di pre-elaborazione quantistica legato all'addestramento del modello e un overhead di post-elaborazione classica per mitigare gli errori nei risultati grezzi usando il modello generato.

L'opzione resilience_level della primitiva Qiskit Runtime specifica il grado di resistenza agli errori da costruire. Livelli più alti producono risultati più accurati al costo di tempi di elaborazione più lunghi dovuti all'overhead di campionamento quantistico. I livelli di resilienza possono essere usati per configurare il compromesso tra costo e accuratezza quando si applica la mitigazione degli errori alla query sulla primitiva.

Quando si implementa una qualsiasi tecnica di mitigazione degli errori, ci aspettiamo che il bias nei nostri risultati si riduca rispetto al bias precedente non mitigato. In alcuni casi, il bias può addirittura scomparire. Tuttavia, questo ha un costo. Man mano che riduciamo il bias nelle nostre quantità stimate, la variabilità statistica aumenterà (ovvero la varianza), di cui possiamo tenere conto aumentando ulteriormente il numero di shot per circuito nel processo di campionamento. Questo introduce un overhead aggiuntivo rispetto a quanto necessario per ridurre il bias, quindi non viene eseguito per impostazione predefinita. Possiamo abilitare facilmente questo comportamento regolando il numero di shot per circuito in options.executions.shots, come mostrato nell'esempio seguente.

Per questo corso, esploreremo questi modelli di mitigazione degli errori ad alto livello per illustrare la mitigazione che le primitive Qiskit Runtime possono eseguire senza richiedere i dettagli di implementazione completi.

Twirled readout error extinction (T-REx)

La Twirled readout error extinction (T-REx) usa una tecnica nota come Pauli twirling per ridurre il rumore introdotto durante il processo di misurazione quantistica. Questa tecnica non assume nessuna forma specifica di rumore, il che la rende molto generale ed efficace.

Flusso di lavoro complessivo:

1. Acquisire dati per lo stato zero con bit flip casuali (Pauli X prima della misurazione)
2. Acquisire dati per lo stato desiderato (rumoroso) con bit flip casuali (Pauli X prima della misurazione)
3. Calcolare la funzione speciale per ogni insieme di dati e dividere.

 

Possiamo impostarlo con options.resilience_level = 1, come dimostrato nell'esempio seguente.

Estrapolazione a rumore zero

L'estrapolazione a rumore zero (ZNE) funziona amplificando prima il rumore nel circuito che prepara lo stato quantistico desiderato, ottenendo misurazioni per diversi livelli di rumore, e usando quelle misurazioni per inferire il risultato privo di rumore.

Flusso di lavoro complessivo:

1. Amplificare il rumore del circuito per diversi fattori di rumore
2. Eseguire ogni circuito con rumore amplificato
3. Estrapolare fino al limite di rumore zero

 

Possiamo impostarlo con options.resilience_level = 2. Possiamo ottimizzare ulteriormente esplorando una varietà di noise_factors, noise_amplifiers ed extrapolators, ma questo va al di là dello scopo di questo corso. Ti incoraggiamo a sperimentare con queste opzioni come descritto qui.

Ogni metodo comporta il proprio overhead associato: un compromesso tra il numero di calcoli quantistici necessari (tempo) e l'accuratezza dei nostri risultati:

$$\begin{array}{c|c|c|c} \text{Methods} & R=1 \text{, T-REx} & R=2 \text{, ZNE} \\[1mm] \hline \text{Assumptions} & \text{None} & \text{Ability to scale noise} \\[1mm] \text{Qubit overhead} & 1 & 1 \\[1mm] \text{Sampling overhead} & 2 & N_{\text{noise-factors}} \\[1mm] \text{Bias} & 0 & \mathcal{O}(\lambda^{N_{\text{noise-factors}}}) \\[1mm] \end{array}$$

Uso delle opzioni di mitigazione e soppressione di Qiskit Runtime

Ecco come calcolare un valore atteso usando la mitigazione e la soppressione degli errori in Qiskit Runtime. Possiamo usare esattamente lo stesso circuito e osservabile di prima, ma questa volta mantenendo il livello di ottimizzazione fisso a 2 e regolando invece la resilienza ovvero le tecniche di mitigazione degli errori utilizzate. Questo processo di mitigazione degli errori avviene più volte durante un ciclo di ottimizzazione.

Eseguiamo questa parte su hardware reale, poiché la mitigazione degli errori non è disponibile sui simulatori.

# Estimated usage: 8 minutes, benchmarked on an Eagle processor, 5-23-24

from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_ibm_runtime import (
    Session,
    EstimatorOptions,
    EstimatorV2 as Estimator,
)

# We select the least busy backend

# Select the least busy backend
# backend = service.least_busy(
#    operational=True, min_num_qubits=ansatz.num_qubits, simulator=False
# )

# Or use a specific backend
backend = service.backend("ibm_brisbane")

# Initialize some variables to save the results from different runs:

exp_values_with_em0_es = []
exp_values_with_em1_es = []
exp_values_with_em2_es = []

# Use a pass manager to optimize the circuit and observables for the backend chosen:

pm = generate_preset_pass_manager(backend=backend, optimization_level=2)
isa_circuit = pm.run(circuit)
isa_observables = observables.apply_layout(layout=isa_circuit.layout)

# Open a session and run with no error mitigation:

estimator_options = EstimatorOptions(resilience_level=0, default_shots=10_000)

with Session(backend=backend) as session:
    estimator = Estimator(mode=session, options=estimator_options)

    pub = (isa_circuit, isa_observables, [individual_phases])
    cost = estimator.run([pub]).result()[0].data.evs

session.close()

exp_values_with_em0_es = cost[0]

# Open a session and run with resilience = 1:

estimator_options = EstimatorOptions(resilience_level=1, default_shots=10_000)

with Session(backend=backend) as session:
    estimator = Estimator(mode=session, options=estimator_options)

    pub = (isa_circuit, isa_observables, [individual_phases])
    cost = estimator.run([pub]).result()[0].data.evs

session.close()

exp_values_with_em1_es = cost[0]

# Open a session and run with resilience = 2:

estimator_options = EstimatorOptions(resilience_level=2, default_shots=10_000)

with Session(backend=backend) as session:
    estimator = Estimator(mode=session, options=estimator_options)

    pub = (isa_circuit, isa_observables, [individual_phases])
    cost = estimator.run([pub]).result()[0].data.evs

session.close()

exp_values_with_em2_es = cost[0]

Come prima, possiamo tracciare i valori attesi risultanti in funzione dell'angolo di fase per i tre livelli di mitigazione degli errori utilizzati. Con grande difficoltà, si può osservare che la mitigazione degli errori migliora leggermente i risultati. Anche qui, questo effetto è molto più pronunciato in circuiti più profondi e complessi.

import matplotlib.pyplot as plt

plt.plot(phases, exp_values_with_em0_es, "o", label="unmitigated")
plt.plot(phases, exp_values_with_em1_es, "o", label="resil = 1")
plt.plot(phases, exp_values_with_em2_es, "o", label="resil = 2")
plt.plot(phases, 2 * np.sin(phases / 2) ** 2 - 1, label="ideal")
plt.ylabel("Expectation")
plt.legend()
plt.show()

Riepilogo

Con questa lezione hai imparato come creare una funzione di costo:

• Creare una funzione di costo
• Come sfruttare le primitive Qiskit Runtime per mitigare e sopprimere il rumore
• Come definire una strategia di misura per ottimizzare velocità vs precisione

Ecco il nostro workload variazionale ad alto livello:

La nostra funzione di costo viene eseguita a ogni iterazione del ciclo di ottimizzazione. La prossima lezione esplorerà come l'ottimizzatore classico usa la valutazione della nostra funzione di costo per selezionare nuovi parametri.

import qiskit
import qiskit_ibm_runtime

print(qiskit.version.get_version_info())
print(qiskit_ibm_runtime.version.get_version_info())

1.1.0
0.23.0