Funcții de cost

În această lecție vom învăța cum să evaluăm o funcție de cost:

• Mai întâi, vom afla despre primitivele Qiskit Runtime
• Vom defini o funcție de cost $C(\vec\theta)$. Aceasta este o funcție specifică problemei care definește obiectivul problemei pentru ca optimizatorul să îl minimizeze (sau maximizeze)
• Vom defini o strategie de măsurare cu primitivele Qiskit Runtime pentru a optimiza viteza față de acuratețe

 

Primitive

Toate sistemele fizice, fie clasice, fie cuantice, pot exista în stări diferite. De exemplu, o mașină pe un drum poate avea o anumită masă, poziție, viteză sau accelerație care îi caracterizează starea. În mod similar, sistemele cuantice pot avea, de asemenea, configurații sau stări diferite, dar se deosebesc de sistemele clasice prin modul în care tratăm măsurătorile și evoluția stărilor. Aceasta conduce la proprietăți unice, cum ar fi superpunerea și entanglementul, exclusive mecanicii cuantice. Așa cum putem descrie starea unei mașini folosind proprietăți fizice precum viteza sau accelerația, putem descrie și starea unui sistem cuantic folosind observabile, care sunt obiecte matematice.

În mecanica cuantică, stările sunt reprezentate prin vectori coloană complecși normalizați, sau kets ($|\psi\rangle$), iar observabilele sunt operatori liniari Hermitian ($\hat{H}=\hat{H}^{\dagger}$) care acționează asupra kets. Un vector propriu ($|\lambda\rangle$) al unui observabil este cunoscut sub numele de stare proprie. Măsurarea unui observabil pentru una dintre stările sale proprii ($|\lambda\rangle$) ne va oferi valoarea proprie corespunzătoare ($\lambda$) ca rezultat.

Dacă te întrebi cum să măsori un sistem cuantic și ce poți măsura, Qiskit oferă două primitive care pot ajuta:

• Sampler: Dată o stare cuantică $|\psi\rangle$, această primitivă obține probabilitatea fiecărei stări posibile din baza computațională.
• Estimator: Dat un observabil cuantic $\hat{H}$ și o stare $|\psi\rangle$, această primitivă calculează valoarea așteptată a lui $\hat{H}$.

Primitiva Sampler

Primitiva Sampler calculează probabilitatea de a obține fiecare stare posibilă $|k\rangle$ din baza computațională, dată un circuit cuantic care pregătește starea $|\psi\rangle$. Aceasta calculează

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

Unde $n$ este numărul de qubiți, iar $k$ este reprezentarea întreagă a oricărui șir binar de ieșire posibil $\{0,1\}^n$ (adică, numere întregi în baza $2$).

Qiskit Runtime Sampler rulează circuitul de mai multe ori pe un dispozitiv cuantic, efectuând măsurători la fiecare rulare, și reconstruind distribuția de probabilitate din șirurile de biți recuperate. Cu cât mai multe rulări (sau shots) efectuează, cu atât rezultatele vor fi mai precise, dar aceasta necesită mai mult timp și resurse cuantice.

Cu toate acestea, deoarece numărul de ieșiri posibile crește exponențial cu numărul de qubiți $n$ (adică, $2^n$), numărul de shots va trebui să crească exponențial și el pentru a capta o distribuție de probabilitate densă. Prin urmare, Sampler este eficient doar pentru distribuții de probabilitate rare (sparse); unde starea țintă $|\psi\rangle$ trebuie să fie exprimabilă ca o combinație liniară a stărilor din baza computațională, cu numărul de termeni crescând cel mult polinomial cu numărul de qubiți:

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

Sampler poate fi, de asemenea, configurat pentru a recupera probabilitățile dintr-o subsecțiune a circuitului, reprezentând un subset din totalul stărilor posibile.

Primitiva Estimator

Primitiva Estimator calculează valoarea de așteptare a unui observabil $\hat{H}$ pentru o stare cuantică $|\psi\rangle$; unde probabilitățile observabilului pot fi exprimate ca $p_\lambda = |\langle\lambda|\psi\rangle|^2$, fiind $|\lambda\rangle$ stările proprii ale observabilului $\hat{H}$. Valoarea de așteptare este definită ca media tuturor rezultatelor posibile $\lambda$ (adică, valorile proprii ale observabilului) ale unei măsurători a stării $|\psi\rangle$, ponderate de probabilitățile corespunzătoare:

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

Cu toate acestea, calcularea valorii de așteptare a unui observabil nu este întotdeauna posibilă, deoarece adesea nu cunoaștem eigenbasis-ul său. Qiskit Runtime Estimator utilizează un proces algebric complex pentru a estima valoarea de așteptare pe un dispozitiv cuantic real, descompunând observabilul într-o combinație de alți observabili al căror eigenbasis îl cunoaștem.

Mai simplu, Estimator descompune orice observabil pe care nu știe cum să îl măsoare în observabile mai simple, măsurabile, numite operatori Pauli. Orice operator poate fi exprimat ca o combinație de $4^n$ operatori 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\}, \\$$

astfel că

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

unde $n$ este numărul de qubiți, $k \equiv k_{n-1} \cdots k_0$ pentru $k_l \in \mathbb{Z}_4 \equiv \{0, 1, 2, 3\}$ (adică, numere întregi în baza $4$), și $(\sigma_0, \sigma_1, \sigma_2, \sigma_3) := (I, X, Y, Z)$.

După efectuarea acestei descompuneri, Estimator derivă un nou circuit $V_k|\psi\rangle$ pentru fiecare observabil $\hat{P}_k$ (din circuitul original), pentru a diagonaliza efectiv observabilul Pauli în baza computațională și a-l măsura. Putem măsura cu ușurință observabilele Pauli deoarece cunoaștem $V_k$ în avans, ceea ce nu este cazul în general pentru alți observabili.

Pentru fiecare $\hat{P}_{k}$, Estimator rulează circuitul corespunzător pe un dispozitiv cuantic de mai multe ori, măsoară starea de ieșire în baza computațională și calculează probabilitatea $p_{kj}$ de a obține fiecare ieșire posibilă $j$. Apoi caută valoarea proprie $\lambda_{kj}$ a lui $P_k$ corespunzătoare fiecărei ieșiri $j$, înmulțește cu $w_k$ și adaugă toate rezultatele împreună pentru a obține valoarea așteptată a observabilului $\hat{H}$ pentru starea dată $|\psi\rangle$.

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

Deoarece calcularea valorii de așteptare a $4^n$ operatori Pauli este impractică (adică, creștere exponențială), Estimator poate fi eficient doar atunci când un număr mare de $w_k$ sunt zero (adică, descompunere Pauli rară în loc de densă). Formal spunem că, pentru ca această calculare să fie eficient rezolvabilă, numărul de termeni nenuli trebuie să crească cel mult polinomial cu numărul de qubiți $n$: $\hat{H} = \sum^{\text{Poly}(n)}_k w_k \hat{P}_k.$

Cititorul poate observa ipoteza implicită că eșantionarea probabilităților trebuie să fie, de asemenea, eficientă, așa cum s-a explicat pentru Sampler, ceea ce înseamnă

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

Exemplu ghidat pentru calculul valorilor de așteptare

Să presupunem starea cu un singur qubit $|+\rangle := H|0\rangle = \frac{1}{\sqrt{2}}(|0\rangle + |1\rangle)$, și observabilul

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

cu următoarea valoare de așteptare teoretică $\langle\hat{H}\rangle_+ = \langle+|\hat{H}|+\rangle = 2.$

Deoarece nu știm cum să măsurăm acest observabil, nu putem calcula direct valoarea sa de așteptare și trebuie să îl reexprimăm ca $\langle\hat{H}\rangle_+ = 2\langle X \rangle_+ - \langle Z \rangle_+$. Se poate arăta că evaluează același rezultat observând că $\langle+|X|+\rangle = 1$ și $\langle+|Z|+\rangle = 0$.

Să vedem cum să calculăm $\langle X \rangle_+$ și $\langle Z \rangle_+$ direct. Deoarece $X$ și $Z$ nu comută (adică, nu împart același eigenbasis), nu pot fi măsurate simultan, prin urmare avem nevoie de circuitele auxiliare:

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

Putem acum efectua calculul manual folosind Sampler și verifica rezultatele cu 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

Rigoare matematică (opțional)

Exprimând $|\psi\rangle$ în raport cu baza stărilor proprii ale lui $\hat{H}$, $|\psi\rangle = \sum_\lambda a_\lambda |\lambda\rangle$, rezultă:

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

Deoarece nu cunoaștem valorile proprii sau stările proprii ale observabilului țintă $\hat{H}$, mai întâi trebuie să luăm în considerare diagonalizarea sa. Dat că $\hat{H}$ este Hermitian, există o transformare unitară $V$ astfel că $\hat{H}=V^\dagger \Lambda V,$ unde $\Lambda$ este matricea diagonală a valorilor proprii, astfel că $\langle j | \Lambda | k \rangle = 0$ dacă $j\neq k$, și $\langle j | \Lambda | j \rangle = \lambda_j$.

Aceasta implică că valoarea de așteptare poate fi rescrisă ca:

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

Dat că dacă un sistem se află în starea $|\phi\rangle = V |\psi\rangle$, probabilitatea de a măsura $| j\rangle$ este $p_j = |\langle j|\phi \rangle|^2$, valoarea de așteptare de mai sus poate fi exprimată ca:

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

Este foarte important de reținut că probabilitățile sunt luate din starea $V |\psi\rangle$ în loc de $|\psi\rangle$. De aceea matricea $V$ este absolut necesară. S-ar putea să te întrebi cum să obții matricea $V$ și valorile proprii $\Lambda$. Dacă ai deja valorile proprii, atunci nu ar mai fi nevoie să folosești un calculator cuantic, deoarece scopul algoritmilor variaționali este tocmai găsirea acestor valori proprii ale lui $\hat{H}$.

Din fericire, există o soluție: orice matrice $2^n \times 2^n$ poate fi scrisă ca o combinație liniară de $4^n$ produse tensoriale de $n$ matrici Pauli și identități, toate fiind atât Hermitian, cât și unitare, cu $V$ și $\Lambda$ cunoscute. Aceasta este ceea ce face Estimator din Runtime intern, descompunând orice obiect Operator într-un SparsePauliOp.

Iată operatorii care pot fi folosiți:

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

Deci să rescriem $\hat{H}$ în raport cu operatorii Pauli și identități:

$$\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,$$

unde $k = \sum_{l=0}^{n-1} 4^l k_l \equiv k_{n-1}...k_0$ pentru $k_{n-1},...,k_0\in \{0,1,2,3\}$ (adică, baza $4$), și $\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}$$

unde $V_k := V_{k_{n-1}}\otimes ... \otimes V_{k_0}$ și $\Lambda_k := \Lambda_{k_{n-1}}\otimes ... \otimes \Lambda_{k_0}$, astfel că: $\hat{P_k}=V_k^\dagger \Lambda_k V_k.$

Funcții de cost

În general, funcțiile de cost sunt folosite pentru a descrie obiectivul unei probleme și cât de bine performează o stare de testare față de acel obiectiv. Această definiție poate fi aplicată la diverse exemple din chimie, învățare automată, finanțe, optimizare și altele.

Să considerăm un exemplu simplu de găsire a stării fundamentale a unui sistem. Obiectivul nostru este să minimizăm valoarea de așteptare a observabilului care reprezintă energia (Hamiltonian $\hat{\mathcal{H}}$):

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

Putem folosi Estimator pentru a evalua valoarea de așteptare și a transmite această valoare unui optimizator pentru a o minimiza. Dacă optimizarea reușește, va returna un set de valori optime ale parametrilor $\vec\theta^*$, din care vom putea construi starea soluție propusă $|\psi(\vec\theta^*)\rangle$ și calcula valoarea de așteptare observată ca $C(\vec\theta^*)$.

Observă că vom putea minimiza funcția de cost doar pentru setul limitat de stări pe care le luăm în considerare. Aceasta ne conduce la două posibilități separate:

• ansatz-ul nostru nu definește starea soluție în tot spațiul de căutare: Dacă acesta este cazul, optimizatorul nu va găsi niciodată soluția și trebuie să experimentăm cu alte ansatz-uri care ar putea reprezenta spațiul nostru de căutare mai precis.
• Optimizatorul nostru este incapabil să găsească această soluție validă: Optimizarea poate fi definită global și local. Vom explora ce înseamnă aceasta în secțiunea ulterioară.

În general, vom efectua o buclă de optimizare clasică, bazându-ne pe evaluarea funcției de cost pe un calculator cuantic. Din această perspectivă, s-ar putea considera optimizarea ca o sarcină pur clasică în care apelăm un oracol cuantic de tip cutie neagră de fiecare dată când optimizatorul trebuie să evalueze funcția de cost.

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

Vom efectua mai întâi acest lucru folosind un simulator: StatevectorEstimator. Aceasta este de obicei recomandabilă pentru depanare, dar vom urma imediat rularea de depanare cu un calcul pe hardware cuantic real. Din ce în ce mai mult, problemele de interes nu mai sunt simulabile clasic fără facilități de supercalculare de ultimă generație.

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

[-0.58744589]

Vom continua acum cu rularea pe un calculator cuantic real. Observă modificările de sintaxă. Pașii care implică pass_manager vor fi discutați mai departe în exemplul următor. Un pas de importanță deosebită în algoritmii variaționali este utilizarea unei sesiuni Qiskit Runtime. Deschiderea unei sesiuni îți permite să rulezi mai multe iterații ale unui algoritm variațional fără a aștepta într-o coadă nouă de fiecare dată când parametrii sunt actualizați. Aceasta este importantă dacă timpii de coadă sunt lungi și/sau sunt necesare multe iterații. Numai partenerii din IBM Quantum® Network pot folosi sesiunile Runtime. Dacă nu ai acces la sesiuni, poți reduce numărul de iterații pe care le trimiți la un moment dat și salva cei mai recenți parametri pentru utilizare în rulări viitoare. Dacă trimiți prea multe iterații sau întâlnești timpi de coadă prea lungi, poți întâlni codul de eroare 1217, care se referă la întârzieri mari între trimiteri de joburi.

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

Observă că valorile obținute din cele două calculări de mai sus sunt foarte similare. Tehnicile pentru îmbunătățirea rezultatelor vor fi discutate mai departe mai jos.

Exemplu de mapare la sisteme ne-fizice

Problema tăieturii maxime (Max-Cut) este o problemă de optimizare combinatorică care implică împărțirea vârfurilor unui graf în două seturi disjuncte astfel încât numărul de muchii dintre cele două seturi să fie maximizat. Mai formal, dat un graf neorientat $G=(V,E)$, unde $V$ este setul de vârfuri și $E$ este setul de muchii, problema Max-Cut cere partiționarea vârfurilor în două subseturi disjuncte, $S$ și $T$, astfel încât numărul de muchii cu un capăt în $S$ și celălalt în $T$ să fie maximizat.

Putem aplica Max-Cut pentru a rezolva diverse probleme, inclusiv: clasificare, proiectarea rețelelor, tranziții de fază și altele. Vom începe prin crearea unui graf al problemei:

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

Această problemă poate fi exprimată ca o problemă de optimizare binară. Pentru fiecare nod $0 \leq i < n$, unde $n$ este numărul de noduri ale grafului (în acest caz $n=4$), vom considera variabila binară $x_i$. Această variabilă va avea valoarea $1$ dacă nodul $i$ se află în unul din grupuri pe care îl vom eticheta $1$ și $0$ dacă se află în celălalt grup, pe care îl vom eticheta ca $0$. De asemenea, vom nota cu $w_{ij}$ (elementul $(i,j)$ al matricei de adiacență $w$) greutatea muchiei care merge de la nodul $i$ la nodul $j$. Deoarece graful este neorientat, $w_{ij}=w_{ji}$. Putem formula astfel problema noastră ca maximizarea următoarei funcții de cost:

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

Pentru a rezolva această problemă cu un calculator cuantic, vom exprima funcția de cost ca valoarea de așteptare a unui observabil. Cu toate acestea, observabilele pe care Qiskit le acceptă nativ constau din operatori Pauli, care au valorile proprii $1$ și $-1$ în loc de $0$ și $1$. De aceea vom face următoarea schimbare de variabilă:

Unde $\vec{x}=(x_0,x_1,\cdots ,x_{n-1})$. Putem folosi matricea de adiacență $w$ pentru a accesa convenabil ponderile tuturor muchiilor. Aceasta va fi folosită pentru a obține funcția noastră de cost:

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

Aceasta implică că:

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

Deci noua funcție de cost pe care dorim să o maximizăm este:

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

Mai mult, tendința naturală a unui calculator cuantic este de a găsi minime (de obicei cea mai joasă energie) în loc de maxime, deci în loc să maximizăm $C(\vec{z})$ vom minimiza:

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

Acum că avem o funcție de cost de minimizat ale cărei variabile pot lua valorile $-1$ și $1$, putem face următoarea analogie cu Pauli $Z$:

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

Cu alte cuvinte, variabila $z_i$ va fi echivalentă cu un Gate $Z$ care acționează pe qubit $i$. Mai mult:

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

Deci observabilul pe care îl vom lua în considerare este:

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

la care va trebui să adăugăm termenul independent ulterior:

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

Operatorul este o combinație liniară de termeni cu operatori Z pe nodurile conectate printr-o muchie (reamintește că qubit-ul 0 se află cel mai în dreapta): $IIZZ + IZIZ + IZZI + ZIIZ + ZZII$. Odată ce operatorul este construit, ansatz-ul pentru algoritmul QAOA poate fi ușor construit folosind circuitul QAOAAnsatz din biblioteca de circuite 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

Cu Estimator din Runtime care primește direct un Hamiltonian și un ansatz parametrizat, și returnează energia necesară, funcția de cost pentru o instanță QAOA este destul de simplă:

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

Vom reveni la acest exemplu în secțiunea Aplicații pentru a explora cum să utilizăm un optimizator pentru a parcurge spațiul de căutare. În general, aceasta include:

• Utilizarea unui optimizator pentru a găsi parametrii optimi
• Legarea parametrilor optimi la ansatz pentru a găsi valorile proprii
• Traducerea valorilor proprii la definiția problemei noastre

Strategia de măsurare: viteză versus acuratețe

Așa cum am menționat, folosim un calculator cuantic zgomotos ca oracol de tip cutie neagră, unde zgomotul poate face valorile recuperate nedeterministe, conducând la fluctuații aleatorii care, la rândul lor, vor afecta — sau chiar împiedica complet — convergența anumitor optimizatori la o soluție propusă. Aceasta este o problemă generală pe care trebuie să o abordăm pe măsură ce explorăm incremental utilitatea cuantică și progresăm spre avantajul cuantic:

Putem folosi opțiunile de suprimare a erorilor și de atenuare a erorilor din Qiskit Runtime Primitives pentru a aborda zgomotul și a maximiza utilitatea calculatoarelor cuantice de astăzi.

Suprimarea erorilor

Suprimarea erorilor se referă la tehnicile utilizate pentru a optimiza și transforma un circuit în timpul compilării pentru a minimiza erorile. Aceasta este o tehnică de bază de gestionare a erorilor care de obicei rezultă în ceva overhead de preprocesare clasică pentru durata totală de execuție. Overhead-ul include transpunerea circuitelor pentru a rula pe hardware cuantic prin:

• Exprimarea circuitului folosind porțile native disponibile pe un sistem cuantic
• Maparea qubiților virtuali la qubiți fizici
• Adăugarea de SWAP-uri pe baza cerințelor de conectivitate
• Optimizarea porților 1Q și 2Q
• Adăugarea de decuplare dinamică la qubiții inactivi pentru a preveni efectele decoherenței.

Primitivele permit utilizarea tehnicilor de suprimare a erorilor prin setarea opțiunii optimization_level și selectarea opțiunilor avansate de transpilare. Într-un curs ulterior, vom aprofunda diferite metode de construcție a circuitelor pentru a îmbunătăți rezultatele, dar pentru majoritatea cazurilor, recomandăm setarea optimization_level=3.

Vom vizualiza valoarea creșterii optimizării în procesul de transpilare prin examinarea unui exemplu de circuit cu un comportament ideal simplu.

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

Circuitul de mai sus poate produce valori de așteptare sinusoidale ale observabilului dat, cu condiția să inserăm faze care acoperă un interval adecvat, cum ar fi $[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]

Putem folosi un simulator pentru a demonstra utilitatea unei transpilări optimizate. Vom reveni mai jos la utilizarea hardware-ului real pentru a demonstra utilitatea atenuării erorilor. Vom folosi QiskitRuntimeService pentru a obține un Backend real (în acest caz, ibm_brisbane), și AerSimulator pentru a simula acel Backend, inclusiv comportamentul său zgomotos.

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)

Putem acum folosi un pass manager pentru a transpune circuitul în „arhitectura setului de instrucțiuni" sau ISA a Backend-ului. Aceasta este o nouă cerință în Qiskit Runtime: toate circuitele trimise unui Backend trebuie să respecte constrângerile țintei Backend-ului, adică trebuie să fie scrise în termenii ISA a Backend-ului — adică setul de instrucțiuni pe care dispozitivul le poate înțelege și executa. Aceste constrângeri ale țintei sunt definite de factori precum porțile native ale dispozitivului, conectivitatea qubiților și — atunci când este relevant — specificațiile de sincronizare a pulsurilor și altor instrucțiuni.

Observă că în prezentul caz, vom face aceasta de două ori: o dată cu optimization_level = 0, și o dată cu acesta setat la 3. De fiecare dată vom folosi primitiva Estimator pentru a estima valorile de așteptare ale observabilului la diferite valori ale fazei.

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

În final, putem reprezenta grafic rezultatele și observăm că precizia calculului a fost destul de bună chiar și fără optimizare, dar s-a îmbunătățit definitiv prin creșterea optimizării la nivelul 3. Observă că în circuite mai profunde și mai complexe, diferența dintre nivelurile de optimizare 0 și 3 este probabil să fie mai semnificativă. Acesta este un circuit foarte simplu folosit ca model de jucărie.

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

Atenuarea erorilor

Atenuarea erorilor se referă la tehnicile care permit utilizatorilor să reducă erorile din circuit prin modelarea zgomotului dispozitivului la momentul execuției. De obicei, aceasta rezultă în overhead de preprocesare cuantică legat de antrenarea modelului și overhead de postprocesare clasică pentru a atenua erorile din rezultatele brute folosind modelul generat.

Opțiunea resilience_level a primitivei Qiskit Runtime specifică cantitatea de reziliență construită împotriva erorilor. Nivelurile mai ridicate generează rezultate mai precise cu prețul unor timpi de procesare mai lungi din cauza overhead-ului de eșantionare cuantică. Nivelurile de reziliență pot fi folosite pentru a configura compromisul dintre cost și acuratețe atunci când se aplică atenuarea erorilor la interogarea ta de primitivă.

Atunci când implementăm orice tehnică de atenuare a erorilor, ne așteptăm ca deviația din rezultatele noastre să fie redusă față de deviația anterioară, neatenuată. În unele cazuri, deviația poate chiar dispărea. Cu toate acestea, aceasta vine cu un cost. Pe măsură ce reducem deviația din cantitățile noastre estimate, variabilitatea statistică va crește (adică varianța), pe care o putem compensa prin creșterea ulterioară a numărului de shots per circuit în procesul nostru de eșantionare. Aceasta va introduce overhead dincolo de cel necesar pentru reducerea deviației, deci nu se face implicit. Putem opta ușor pentru acest comportament ajustând numărul de shots per circuit în options.executions.shots, după cum se arată în exemplul de mai jos.

Pentru acest curs, vom explora aceste modele de atenuare a erorilor la un nivel înalt pentru a ilustra atenuarea erorilor pe care primitivele Qiskit Runtime le pot efectua fără a necesita detalii complete de implementare.

Extinctia erorii de citire prin twirling (T-REx)

Extinctia erorii de citire prin twirling (T-REx) folosește o tehnică cunoscută sub numele de Pauli twirling pentru a reduce zgomotul introdus în timpul procesului de măsurare cuantică. Această tehnică nu presupune nicio formă specifică de zgomot, ceea ce o face foarte generală și eficientă.

Fluxul de lucru general:

1. Achiziționarea datelor pentru starea zero cu flip-uri de biți aleatorii (Pauli X înainte de măsurare)
2. Achiziționarea datelor pentru starea dorită (zgomotoasă) cu flip-uri de biți aleatorii (Pauli X înainte de măsurare)
3. Calcularea funcției speciale pentru fiecare set de date și împărțirea.

 

Putem seta aceasta cu options.resilience_level = 1, demonstrat în exemplul de mai jos.

Extrapolarea la zero zgomot

Extrapolarea la zero zgomot (ZNE) funcționează prin amplificarea mai întâi a zgomotului din circuitul care pregătește starea cuantică dorită, obținerea măsurătorilor pentru mai multe niveluri diferite de zgomot, și utilizarea acelor măsurători pentru a deduce rezultatul fără zgomot.

Fluxul de lucru general:

1. Amplificarea zgomotului circuitului pentru mai mulți factori de zgomot
2. Rularea fiecărui circuit cu zgomot amplificat
3. Extrapolarea înapoi la limita de zero zgomot

 

Putem seta aceasta cu options.resilience_level = 2. Putem optimiza aceasta în continuare prin explorarea unei varietăți de noise_factors, noise_amplifiers și extrapolators, dar aceasta este în afara domeniului acestui curs. Te încurajăm să experimentezi cu aceste opțiuni descrise aici.

Fiecare metodă vine cu propriul overhead asociat: un compromis între numărul de calcule cuantice necesare (timp) și acuratețea rezultatelor noastre:

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

Utilizarea opțiunilor de atenuare și suprimare din Qiskit Runtime

Iată cum să calculezi o valoare de așteptare folosind atenuarea și suprimarea erorilor în Qiskit Runtime. Putem folosi exact același circuit și observabil ca înainte, dar de această dată păstrând nivelul de optimizare fix la nivelul 2 și reglând acum reziliența sau tehnica(ile) de atenuare a erorilor utilizate. Acest proces de atenuare a erorilor are loc de mai multe ori pe parcursul unei bucle de optimizare.

Efectuăm această parte pe hardware real, deoarece atenuarea erorilor nu este disponibilă pe simulatoare.

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

Ca și înainte, putem reprezenta grafic valorile de așteptare rezultate în funcție de unghiul de fază pentru cele trei niveluri de atenuare a erorilor utilizate. Cu mare dificultate, se poate observa că atenuarea erorilor îmbunătățește ușor rezultatele. Din nou, acest efect este mult mai pronunțat în circuite mai profunde și mai complexe.

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

Rezumat

Cu această lecție, ai învățat cum să creezi o funcție de cost:

• Crearea unei funcții de cost
• Cum să utilizezi primitivele Qiskit Runtime pentru a atenua și suprima zgomotul
• Cum să definești o strategie de măsurare pentru a optimiza viteza față de acuratețe

Iată sarcina de lucru variaționistă la nivel înalt:

Funcția noastră de cost rulează la fiecare iterație a buclei de optimizare. Lecția următoare va explora modul în care optimizatorul clasic utilizează evaluarea funcției noastre de cost pentru a selecta parametri noi.

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