Nákladové funkce

Během této lekce se naučíme, jak vyhodnotit nákladovou funkci:

• Nejprve se seznámíme s primitivy Qiskit Runtime
• Definujeme nákladovou funkci $C(\vec\theta)$. Jde o funkci specifickou pro daný problém, která definuje cíl problému, jenž má optimalizátor minimalizovat (nebo maximalizovat)
• Definování strategie měření pomocí primitiv Qiskit Runtime pro optimalizaci rychlosti vůči přesnosti

 

Primitiva

Všechny fyzikální systémy, ať už klasické nebo kvantové, mohou existovat v různých stavech. Například auto na silnici může mít určitou hmotnost, polohu, rychlost nebo zrychlení, které charakterizují jeho stav. Podobně i kvantové systémy mohou mít různé konfigurace nebo stavy, ale od klasických systémů se liší tím, jak zacházíme s měřeními a s vývojem stavů. To vede k jedinečným vlastnostem, jako je superpozice a provázanost, které jsou výhradně záležitostí kvantové mechaniky. Stejně jako můžeme popsat stav auta pomocí fyzikálních veličin, jako je rychlost nebo zrychlení, můžeme popsat i stav kvantového systému pomocí pozorovatelných veličin, což jsou matematické objekty.

V kvantové mechanice jsou stavy reprezentovány normalizovanými komplexními sloupcovými vektory neboli kety ($|\psi\rangle$) a pozorovatelné veličiny jsou hermitovské lineární operátory ($\hat{H}=\hat{H}^{\dagger}$), které na tyto kety působí. Vlastní vektor ($|\lambda\rangle$) pozorovatelné veličiny se nazývá vlastní stav. Měřením pozorovatelné veličiny na některém z jejích vlastních stavů ($|\lambda\rangle$) získáme jako výstup odpovídající vlastní hodnotu ($\lambda$).

Pokud se ptáš, jak změřit kvantový systém a co vlastně můžeš měřit, Qiskit ti nabízí dvě primitiva, která mohou pomoci:

• Sampler: Pro daný kvantový stav $|\psi\rangle$ toto primitivum získává pravděpodobnost každého možného stavu výpočetní báze.
• Estimator: Pro danou kvantovou pozorovatelnou veličinu $\hat{H}$ a stav $|\psi\rangle$ toto primitivum počítá střední hodnotu $\hat{H}$.

Primitivum Sampler

Primitivum Sampler počítá pravděpodobnost získání každého možného stavu $|k\rangle$ z výpočetní báze pro daný kvantový circuit, který připravuje stav $|\psi\rangle$. Počítá

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

kde $n$ je počet qubitů a $k$ celočíselná reprezentace libovolného možného výstupního binárního řetězce $\{0,1\}^n$ (tedy celá čísla v základu $2$).

Qiskit Runtime Sampler spouští circuit na kvantovém zařízení opakovaně, při každém běhu provádí měření a z získaných bitových řetězců rekonstruuje rozdělení pravděpodobností. Čím více běhů (neboli shotů) provede, tím přesnější budou výsledky, to ale vyžaduje více času a kvantových prostředků.

Protože však počet možných výstupů roste exponenciálně s počtem qubitů $n$ (tedy $2^n$), bude muset i počet shotů růst exponenciálně, aby zachytil husté rozdělení pravděpodobnosti. Proto je Sampler efektivní pouze pro řídká rozdělení pravděpodobnosti; kde cílový stav $|\psi\rangle$ musí být vyjádřitelný jako lineární kombinace stavů výpočetní báze, přičemž počet členů roste nejvýše polynomiálně s počtem qubitů:

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

Sampler lze také nakonfigurovat tak, aby získával pravděpodobnosti z podčásti circuitu, reprezentující podmnožinu všech možných stavů.

Primitivum Estimator

Primitivum Estimator počítá střední hodnotu pozorovatelné veličiny $\hat{H}$ pro kvantový stav $|\psi\rangle$; kde pravděpodobnosti pozorovatelné veličiny lze vyjádřit jako $p_\lambda = |\langle\lambda|\psi\rangle|^2$, přičemž $|\lambda\rangle$ jsou vlastní stavy pozorovatelné veličiny $\hat{H}$. Střední hodnota je pak definována jako průměr všech možných výsledků $\lambda$ (tedy vlastních hodnot pozorovatelné veličiny) měření stavu $|\psi\rangle$, vážený odpovídajícími pravděpodobnostmi:

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

Výpočet střední hodnoty pozorovatelné veličiny však není vždy možný, protože často neznáme její vlastní bázi. Qiskit Runtime Estimator využívá složitý algebraický postup k odhadu střední hodnoty na reálném kvantovém zařízení tak, že pozorovatelnou veličinu rozloží na kombinaci jiných pozorovatelných, jejichž vlastní bázi známe.

Zjednodušeně řečeno, Estimator rozkládá libovolnou pozorovatelnou veličinu, kterou neumí změřit, na jednodušší, měřitelné pozorovatelné zvané Pauliho operátory. Libovolný operátor lze vyjádřit jako kombinaci $4^n$ Pauliho operátorů.

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

tak, že

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

kde $n$ je počet qubitů, $k \equiv k_{n-1} \cdots k_0$ pro $k_l \in \mathbb{Z}_4 \equiv \{0, 1, 2, 3\}$ (tedy celá čísla v základu $4$) a $(\sigma_0, \sigma_1, \sigma_2, \sigma_3) := (I, X, Y, Z)$.

Po provedení tohoto rozkladu Estimator odvodí nový circuit $V_k|\psi\rangle$ pro každou pozorovatelnou veličinu $\hat{P}_k$ (z původního circuitu), aby efektivně diagonalizoval Pauliho pozorovatelnou ve výpočetní bázi a změřil ji. Pauliho pozorovatelné můžeme snadno změřit, protože $V_k$ známe předem, což obecně pro jiné pozorovatelné neplatí.

Pro každé $\hat{P}_{k}$ Estimator spustí odpovídající circuit na kvantovém zařízení několikrát, změří výstupní stav ve výpočetní bázi a vypočítá pravděpodobnost $p_{kj}$ získání každého možného výstupu $j$. Poté najde vlastní hodnotu $\lambda_{kj}$ operátoru $P_k$ odpovídající každému výstupu $j$, vynásobí ji $w_k$ a všechny výsledky sečte, čímž získá střední hodnotu pozorovatelné veličiny $\hat{H}$ pro daný stav $|\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},$$

Protože výpočet střední hodnoty $4^n$ Pauliho operátorů je nepraktický (tedy roste exponenciálně), Estimator může být efektivní pouze tehdy, když je velké množství $w_k$ nulových (tedy řídký Pauliho rozklad namísto hustého). Formálně říkáme, že aby tento výpočet byl efektivně řešitelný, musí počet nenulových členů růst nejvýše polynomiálně s počtem qubitů $n$: $\hat{H} = \sum^{\text{Poly}(n)}_k w_k \hat{P}_k.$

Čtenář si může povšimnout implicitního předpokladu, že pravděpodobnostní vzorkování musí být také efektivní, jak bylo vysvětleno pro Sampler, což znamená

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

Vedený příklad výpočtu středních hodnot

Předpokládejme jednoqubitový stav $|+\rangle := H|0\rangle = \frac{1}{\sqrt{2}}(|0\rangle + |1\rangle)$ a pozorovatelnou

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

s následující teoretickou střední hodnotou $\langle\hat{H}\rangle_+ = \langle+|\hat{H}|+\rangle = 2.$

Protože neumíme tuto pozorovatelnou změřit přímo, nemůžeme její střední hodnotu spočítat přímo a musíme ji přepsat jako $\langle\hat{H}\rangle_+ = 2\langle X \rangle_+ - \langle Z \rangle_+$. Lze ukázat, že se vyhodnotí na stejný výsledek, vezmeme-li v úvahu, že $\langle+|X|+\rangle = 1$ a $\langle+|Z|+\rangle = 0$.

Podívejme se, jak spočítat $\langle X \rangle_+$ a $\langle Z \rangle_+$ přímo. Protože $X$ a $Z$ nekomutují (tj. nesdílejí stejnou vlastní bázi), nemohou být měřeny současně, a proto potřebujeme pomocné obvody:

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

Nyní můžeme výpočet provést ručně pomocí Sampleru a zkontrolovat výsledky pomocí Estimatoru:

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

Matematická přísnost (volitelné)

Vyjádříme-li $|\psi\rangle$ v bázi vlastních stavů $\hat{H}$, $|\psi\rangle = \sum_\lambda a_\lambda |\lambda\rangle$, dostaneme:

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

Protože neznáme vlastní hodnoty ani vlastní stavy cílové pozorovatelné $\hat{H}$, musíme nejprve uvažovat její diagonalizaci. Vzhledem k tomu, že $\hat{H}$ je hermitovská, existuje unitární transformace $V$ taková, že $\hat{H}=V^\dagger \Lambda V,$ kde $\Lambda$ je diagonální matice vlastních hodnot, tedy $\langle j | \Lambda | k \rangle = 0$ pro $j\neq k$ a $\langle j | \Lambda | j \rangle = \lambda_j$.

Z toho plyne, že střední hodnotu lze přepsat jako:

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

Pokud je systém ve stavu $|\phi\rangle = V |\psi\rangle$, je pravděpodobnost naměření $| j\rangle$ rovna $p_j = |\langle j|\phi \rangle|^2$, a výše uvedenou střední hodnotu lze tedy vyjádřit jako:

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

Je velmi důležité poznamenat, že pravděpodobnosti jsou brány ze stavu $V |\psi\rangle$, nikoliv ze stavu $|\psi\rangle$. Proto je matice $V$ naprosto nezbytná. Možná tě napadá, jak matici $V$ a vlastní hodnoty $\Lambda$ získat. Kdybys vlastní hodnoty již měl/a k dispozici, nebylo by potřeba používat kvantový počítač, protože cílem variačních algoritmů je právě tyto vlastní hodnoty $\hat{H}$ nalézt.

Naštěstí existuje způsob, jak tuto překážku obejít: každou matici $2^n \times 2^n$ lze zapsat jako lineární kombinaci $4^n$ tenzorových součinů $n$ Pauliho matic a jednotkových matic, přičemž všechny jsou hermitovské i unitární se známými $V$ a $\Lambda$. To je právě to, co Runtime Estimator dělá interně – rozkládá libovolný objekt Operator na SparsePauliOp.

Zde jsou Operátory, které lze použít:

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

Přepišme tedy $\hat{H}$ pomocí Pauliho matic a jednotkových matic:

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

kde $k = \sum_{l=0}^{n-1} 4^l k_l \equiv k_{n-1}...k_0$ pro $k_{n-1},...,k_0\in \{0,1,2,3\}$ (tj. v základu $4$) a $\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}$$

kde $V_k := V_{k_{n-1}}\otimes ... \otimes V_{k_0}$ a $\Lambda_k := \Lambda_{k_{n-1}}\otimes ... \otimes \Lambda_{k_0}$, takže: $\hat{P_k}=V_k^\dagger \Lambda_k V_k.$

Cenové funkce

Cenové funkce se obecně používají k popisu cíle problému a k tomu, jak dobře si zkušební stav vede vzhledem k tomuto cíli. Tuto definici lze aplikovat na různé příklady v chemii, strojovém učení, financích, optimalizaci a tak dále.

Uvažujme jednoduchý příklad hledání základního stavu systému. Naším cílem je minimalizovat střední hodnotu pozorovatelné veličiny představující energii (Hamiltonián $\hat{\mathcal{H}}$):

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

Pomocí Estimator můžeme vyhodnotit střední hodnotu a předat tuto hodnotu optimalizátoru k minimalizaci. Pokud je optimalizace úspěšná, vrátí sadu optimálních hodnot parametrů $\vec\theta^*$, ze kterých budeme moci sestavit navrhovaný stavový vektor $|\psi(\vec\theta^*)\rangle$ a vypočítat pozorovanou střední hodnotu jako $C(\vec\theta^*)$.

Všimni si, že cenovou funkci budeme moci minimalizovat pouze pro omezený soubor stavů, které zvažujeme. To nás vede ke dvěma různým možnostem:

• Náš ansatz nedefinuje řešení v celém prohledávaném prostoru: V takovém případě optimalizátor řešení nikdy nenajde a musíme experimentovat s jinými ansatzy, které by mohly lépe reprezentovat náš prohledávaný prostor.
• Náš optimalizátor není schopen toto platné řešení nalézt: Optimalizace může být definována globálně i lokálně. Co to znamená, prozkoumáme v pozdější části.

Celkově vzato budeme provádět klasickou optimalizační smyčku, přičemž se budeme spoléhat na vyhodnocení cenové funkce kvantovým počítačem. Z tohoto pohledu by bylo možné na optimalizaci nahlížet jako na čistě klasické úsilí, při němž pokaždé, když optimalizátor potřebuje vyhodnotit cenovou funkci, voláme jakýsi kvantový orákulum černé skříňky.

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

Nejprve to provedeme pomocí simulátoru: StatevectorEstimator. To se obvykle doporučuje pro ladění, ale ihned po ladicím spuštění provedeme výpočet na skutečném kvantovém hardware. Problémy, které nás zajímají, jsou stále méně klasicky simulovatelné bez nejmodernějších superpočítačových zařízení.

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

[-0.58744589]

Nyní přistoupíme ke spuštění na skutečném kvantovém počítači. Všimni si změn v syntaxi. Kroky zahrnující pass_manager budou podrobněji probrány v dalším příkladu. Jedním obzvláště důležitým krokem ve variačních algoritmech je použití Session Qiskit Runtime. Zahájení session ti umožňuje spouštět více iterací variačního algoritmu, aniž bys pokaždé čekal/a v nové frontě, když se aktualizují parametry. To je důležité, pokud jsou časy front dlouhé nebo je potřeba mnoho iterací. Runtime sessions mohou používat pouze partneři v síti IBM Quantum® Network. Pokud k sessions nemáš přístup, můžeš snížit počet iterací, které odešleš najednou, a uložit nejnovější parametry pro použití v budoucích spuštěních. Pokud odešleš příliš mnoho iterací nebo narazíš na příliš dlouhé časy front, může se objevit chybový kód 1217, který označuje dlouhé prodlevy mezi odesláním úloh.

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

Všimni si, že hodnoty získané z obou výpočtů výše jsou velmi podobné. Techniky pro zlepšení výsledků budou dále probrány níže.

Příklad mapování na nefyzikální systémy

Problém maximálního řezu (Max-Cut) je kombinatorický optimalizační problém, který spočívá v rozdělení vrcholů grafu do dvou disjunktních množin tak, aby byl maximalizován počet hran mezi těmito dvěma množinami. Formálněji: pro neorientovaný graf $G=(V,E)$, kde $V$ je množina vrcholů a $E$ je množina hran, se problém Max-Cut ptá, jak rozdělit vrcholy do dvou disjunktních podmnožin $S$ a $T$ tak, aby byl maximalizován počet hran, jejichž jeden koncový bod leží v $S$ a druhý v $T$.

Max-Cut lze použít k řešení různých problémů, včetně: shlukování, návrhu sítí, fázových přechodů a dalších. Začneme vytvořením grafu problému:

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

Tento problém lze vyjádřit jako binární optimalizační problém. Pro každý uzel $0 \leq i < n$, kde $n$ je počet uzlů grafu (v tomto případě $n=4$), budeme uvažovat binární proměnnou $x_i$. Tato proměnná bude mít hodnotu $1$, pokud uzel $i$ patří do skupiny, kterou označíme jako $1$, a hodnotu $0$, pokud patří do druhé skupiny, kterou označíme jako $0$. Budeme také označovat jako $w_{ij}$ (prvek $(i,j)$ matice sousednosti $w$) váhu hrany vedoucí z uzlu $i$ do uzlu $j$. Protože je graf neorientovaný, platí $w_{ij}=w_{ji}$. Náš problém pak lze formulovat jako maximalizaci následující účelové funkce:

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

Abychom mohli tento problém řešit na kvantovém počítači, vyjádříme účelovou funkci jako střední hodnotu pozorovatelné veličiny. Pozorovatelné veličiny, které Qiskit nativně podporuje, se skládají z Pauliho operátorů s vlastními hodnotami $1$ a $-1$ místo $0$ a $1$. Proto provedeme následující změnu proměnných:

Kde $\vec{x}=(x_0,x_1,\cdots ,x_{n-1})$. Matici sousednosti $w$ lze pohodlně využít pro přístup k váhám všech hran. Ta bude použita k získání naší účelové funkce:

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

Z toho plyne:

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

Nová účelová funkce, kterou chceme maximalizovat, je tedy:

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

Kvantový počítač navíc přirozeně hledá minima (zpravidla nejnižší energii) namísto maxim, takže místo maximalizace $C(\vec{z})$ budeme minimalizovat:

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

Nyní, když máme účelovou funkci k minimalizaci, jejíž proměnné mohou nabývat hodnot $-1$ a $1$, můžeme provést následující analogii s Pauliho operátorem $Z$:

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

Jinými slovy, proměnná $z_i$ bude ekvivalentní Gate $Z$ působícímu na Qubit $i$. Navíc:

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

Pozorovatelná veličina, kterou budeme uvažovat, je tedy:

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

ke které bude třeba následně přičíst konstantní člen:

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

Operátor je lineární kombinací členů s operátory Z na uzlech spojených hranou (připomeň, že 0. Qubit je nejdále vpravo): $IIZZ + IZIZ + IZZI + ZIIZ + ZZII$. Jakmile je operátor sestaven, ansatz pro algoritmus QAOA lze snadno vytvořit pomocí Circuit QAOAAnsatz z knihovny Qiskit circuit library.

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

Protože Estimator z Runtime přímo přijímá hamiltonián a parametrizovaný ansatz a vrací potřebnou energii, je účelová funkce pro instanci QAOA poměrně jednoduchá:

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

K tomuto příkladu se vrátíme v části Aplikace, kde prozkoumáme, jak využít optimalizátor k procházení prohledávaného prostoru. Obecně to zahrnuje:

• Využití optimalizátoru k nalezení optimálních parametrů
• Navázání optimálních parametrů na ansatz za účelem nalezení vlastních hodnot
• Převod vlastních hodnot na definici našeho problému

Strategie měření: rychlost versus přesnost

Jak bylo zmíněno, používáme zašuměný kvantový počítač jako black-box orákulum, kde šum může způsobit, že získané hodnoty budou nedeterministické, což vede k náhodným výkyvům, jež zase mohou uškodit — nebo dokonce zcela zabránit — konvergenci určitých optimalizátorů k navrhovanému řešení. Jde o obecný problém, který musíme řešit, jak postupně zkoumáme kvantovou užitečnost a pokračujeme směrem ke kvantové výhodě:

Ke zvládnutí šumu a maximalizaci užitečnosti dnešních kvantových počítačů můžeme využít možnosti potlačení chyb a zmírnění chyb primitiv Qiskit Runtime.

Potlačení chyb

Potlačení chyb označuje techniky používané k optimalizaci a transformaci Circuit během kompilace s cílem minimalizovat chyby. Jde o základní techniku zpracování chyb, která obvykle přináší určitou režii klasického předzpracování do celkové doby běhu. Tato režie zahrnuje transpilaci Circuit pro spuštění na kvantovém hardwaru prostřednictvím:

• Vyjádření Circuit pomocí nativních Gate dostupných na kvantovém systému
• Mapování virtuálních Qubit na fyzické Qubit
• Přidání SWAP operací na základě požadavků na konektivitu
• Optimalizace 1Q a 2Q Gate
• Přidání dynamického oddělování (dynamical decoupling) k nečinným Qubit, aby se předešlo efektům dekoherence.

Primitiva umožňují využívat techniky potlačení chyb nastavením volby optimization_level a výběrem pokročilých možností transpilace. V pozdějším kurzu se budeme podrobněji věnovat různým metodám konstrukce Circuit pro zlepšení výsledků, ale pro většinu případů doporučujeme nastavit optimization_level=3.

Hodnotu zvyšující se optimalizace v procesu transpilace vizualizujeme na příkladu Circuit s jednoduchou ideální chování.

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

Výše uvedená Circuit může poskytovat sinusoidální střední hodnoty daného pozorovatelné veličiny za předpokladu, že vložíme fáze pokrývající vhodný interval, například $[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]

Pomocí simulátoru můžeme ukázat užitečnost optimalizované transpilace. Níže se vrátíme k používání skutečného hardwaru, abychom demonstrovali užitečnost zmírnění chyb. Pomocí QiskitRuntimeService získáme skutečný Backend (v tomto případě ibm_brisbane) a použijeme AerSimulator k simulaci tohoto Backend včetně jeho šumového chování.

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)

Nyní můžeme použít pass manager k transpilaci Circuit do „instrukční sady architektury" (ISA) Backend. Toto je nový požadavek v Qiskit Runtime: všechny Circuit odeslané do Backend musí odpovídat omezením cíle Backend, to znamená, že musí být zapsány ve smyslu ISA Backend — tedy sady instrukcí, které zařízení dokáže pochopit a vykonat. Tato cílová omezení jsou definována faktory jako nativní základní Gate zařízení, konektivita jeho Qubit a — pokud je to relevantní — specifikace pulzů a dalších instrukčních časování.

Všimni si, že v tomto případě to provedeme dvakrát: jednou s optimization_level = 0 a jednou s hodnotou nastavenou na 3. Pokaždé použijeme primitiv Estimator k odhadnutí středních hodnot pozorovatelné veličiny při různých hodnotách fáze.

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

Nakonec můžeme zobrazit výsledky a vidíme, že přesnost výpočtu byla poměrně dobrá i bez optimalizace, ale rozhodně se zlepšila zvýšením optimalizace na úroveň 3. Všimni si, že u hlubších a složitějších Circuit bude rozdíl mezi úrovněmi optimalizace 0 a 3 pravděpodobně výraznější. Toto je velmi jednoduchá Circuit použitá jako hračkový model.

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

Zmírnění chyb

Zmírnění chyb označuje techniky, které umožňují uživatelům snižovat chyby v Circuit tím, že modelují šum zařízení v době provádění. Typicky to vede k overhead při kvantovém předzpracování spojeném s trénováním modelu a k overhead při klasickém postprocesingu, který slouží ke zmírnění chyb v surových výsledcích pomocí vygenerovaného modelu.

Možnost resilience_level u Qiskit Runtime primitivů určuje míru odolnosti vůči chybám. Vyšší úrovně generují přesnější výsledky za cenu delší doby zpracování kvůli overhead při kvantovém vzorkování. Úrovně odolnosti lze použít ke konfiguraci kompromisu mezi náklady a přesností při aplikaci zmírnění chyb na dotaz na primitiva.

Při implementaci jakékoli techniky zmírnění chyb očekáváme, že bias v našich výsledcích bude snížen oproti předchozímu, nemitigovanému biasu. V některých případech může bias dokonce zmizet. To však přichází za cenu. Jak snižujeme bias v odhadovaných veličinách, statistická variabilita se zvýší (tedy rozptyl), což můžeme kompenzovat dalším zvýšením počtu shotů na Circuit v našem procesu vzorkování. To způsobí overhead nad rámec toho, který je potřeba ke snížení biasu, takže se to ve výchozím nastavení neprovádí. Toto chování můžeme snadno zapnout úpravou počtu shotů na Circuit v options.executions.shots, jak je ukázáno v níže uvedeném příkladu.

V tomto kurzu budeme zkoumat tyto modely zmírnění chyb na vysoké úrovni, abychom ilustrovali zmírnění chyb, které mohou Qiskit Runtime primitivy provádět, aniž by bylo nutné znát veškeré detaily implementace.

Twirled readout error extinction (T-REx)

Twirled readout error extinction (T-REx) využívá techniku zvanou Pauli twirling ke snížení šumu zavedeného při procesu kvantového měření. Tato technika nepředpokládá žádnou konkrétní formu šumu, což ji činí velmi obecnou a účinnou.

Celkový postup:

1. Získej data pro nulový stav s náhodným překlápěním bitů (Pauli X před měřením)
2. Získej data pro požadovaný (zašuměný) stav s náhodným překlápěním bitů (Pauli X před měřením)
3. Vypočítej speciální funkci pro každou datovou sadu a vyděl.

 

Toto nastavení lze provést pomocí options.resilience_level = 1, jak je ukázáno v níže uvedeném příkladu.

Zero noise extrapolation

Zero noise extrapolation (ZNE) funguje tak, že nejprve zesiluje šum v Circuit, která připravuje požadovaný kvantový stav, získává měření pro několik různých úrovní šumu a pomocí těchto měření odvozuje výsledek bez šumu.

Celkový postup:

1. Zesil šum Circuit pro několik šumových faktorů
2. Spusť každý Circuit se zesíleným šumem
3. Extrapoluj zpět na limit nulového šumu

 

Toto nastavení lze provést pomocí options.resilience_level = 2. Můžeme to dále optimalizovat průzkumem různých noise_factors, noise_amplifiers a extrapolators, ale to je mimo rozsah tohoto kurzu. Doporučujeme ti experimentovat s těmito možnostmi, jak jsou popsány zde.

Každá metoda přináší vlastní overhead: kompromis mezi počtem potřebných kvantových výpočtů (čas) a přesností našich výsledků:

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

Použití možností zmírnění a potlačení chyb v Qiskit Runtime

Zde je ukázka, jak vypočítat střední hodnotu při použití zmírnění a potlačení chyb v Qiskit Runtime. Můžeme použít přesně stejný Circuit a pozorovatelnou jako dříve, tentokrát však s pevně nastavenou úrovní optimalizace na úrovni 2 a s laděním odolnosti neboli použité techniky (technik) zmírnění chyb. Tento proces zmírnění chyb probíhá vícekrát v průběhu optimalizační smyčky.

Tuto část provádíme na reálném hardwaru, protože zmírnění chyb není k dispozici na simulátorech.

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

Stejně jako dříve můžeme vykreslit výsledné střední hodnoty jako funkci fázového úhlu pro tři použité úrovně zmírnění chyb. S velkým úsilím lze vidět, že zmírnění chyb výsledky mírně zlepšuje. Tento efekt je opět mnohem výraznější u hlubších a složitějších Circuitů.

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

Shrnutí

V této lekci ses naučil(a), jak vytvořit cost funkci:

• Vytvoření cost funkce
• Jak využít primitiva Qiskit Runtime ke zmírnění a potlačení šumu
• Jak definovat strategii měření pro optimalizaci rychlosti versus přesnosti

Zde je náš přehled variačního pracovního postupu:

Naše cost funkce se spouští při každé iteraci optimalizační smyčky. Následující lekce prozkoumá, jak klasický optimalizátor využívá vyhodnocení naší cost funkce k výběru nových parametrů.

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