Funkcje kosztu

Podczas tej lekcji nauczymy się, jak obliczać funkcję kosztu:

• Najpierw poznamy prymitywy Qiskit Runtime
• Zdefiniujemy funkcję kosztu $C(\vec\theta)$. Jest to funkcja specyficzna dla problemu, która określa cel problemu do minimalizacji (lub maksymalizacji) przez optymalizator
• Zdefiniujemy strategię pomiaru z użyciem prymitywów Qiskit Runtime, aby zoptymalizować szybkość względem dokładności

 

Prymitywy

Wszystkie układy fizyczne, zarówno klasyczne, jak i kwantowe, mogą istnieć w różnych stanach. Na przykład samochód na drodze może mieć określoną masę, położenie, prędkość lub przyspieszenie, które charakteryzują jego stan. Podobnie układy kwantowe również mogą mieć różne konfiguracje lub stany, ale różnią się one od układów klasycznych tym, jak traktujemy pomiary i ewolucję stanu. Prowadzi to do unikalnych właściwości, takich jak superpozycja i splątanie, które są wyłączne dla mechaniki kwantowej. Podobnie jak możemy opisać stan samochodu przy użyciu właściwości fizycznych, takich jak prędkość czy przyspieszenie, tak samo możemy opisać stan układu kwantowego za pomocą obserwabli, które są obiektami matematycznymi.

W mechanice kwantowej stany są reprezentowane przez znormalizowane zespolone wektory kolumnowe, czyli kety ($|\psi\rangle$), a obserwable są hermitowskimi operatorami liniowymi ($\hat{H}=\hat{H}^{\dagger}$), które działają na ketach. Wektor własny ($|\lambda\rangle$) obserwabli jest nazywany stanem własnym. Pomiar obserwabli dla jednego z jej stanów własnych ($|\lambda\rangle$) da nam odpowiadającą wartość własną ($\lambda$) jako odczyt.

Jeśli zastanawiasz się, jak mierzyć układ kwantowy i co można mierzyć, Qiskit oferuje dwa prymitywy, które mogą pomóc:

• Sampler: Mając dany stan kwantowy $|\psi\rangle$, ten prymityw uzyskuje prawdopodobieństwo każdego możliwego stanu bazy obliczeniowej.
• Estimator: Mając daną obserwablę kwantową $\hat{H}$ i stan $|\psi\rangle$, ten prymityw oblicza wartość oczekiwaną $\hat{H}$.

Prymityw Sampler

Prymityw Sampler oblicza prawdopodobieństwo uzyskania każdego możliwego stanu $|k\rangle$ z bazy obliczeniowej, mając dany obwód kwantowy, który przygotowuje stan $|\psi\rangle$. Oblicza on

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

Gdzie $n$ to liczba kubitów, a $k$ to całkowita reprezentacja dowolnego możliwego wyjściowego ciągu binarnego $\{0,1\}^n$ (czyli liczby całkowite w bazie $2$).

Qiskit Runtime Sampler wykonuje obwód wielokrotnie na urządzeniu kwantowym, przeprowadzając pomiary przy każdym uruchomieniu, i rekonstruuje rozkład prawdopodobieństwa z odzyskanych ciągów bitów. Im więcej uruchomień (lub strzałów) wykonuje, tym dokładniejsze będą wyniki, ale wymaga to więcej czasu i zasobów kwantowych.

Jednakże, ponieważ liczba możliwych wyników rośnie wykładniczo wraz z liczbą kubitów $n$ (czyli $2^n$), liczba strzałów również będzie musiała rosnąć wykładniczo, aby uchwycić gęsty rozkład prawdopodobieństwa. Dlatego Sampler jest wydajny tylko dla rzadkich rozkładów prawdopodobieństwa; gdzie stan docelowy $|\psi\rangle$ musi być wyrażalny jako kombinacja liniowa stanów bazy obliczeniowej, przy czym liczba składników rośnie co najwyżej wielomianowo wraz z liczbą kubitów:

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

Sampler można również skonfigurować tak, aby pobierał prawdopodobieństwa z podsekcji obwodu, reprezentując podzbiór wszystkich możliwych stanów.

Prymityw Estimator

Prymityw Estimator oblicza wartość oczekiwaną obserwabli $\hat{H}$ dla stanu kwantowego $|\psi\rangle$; gdzie prawdopodobieństwa obserwabli można wyrazić jako $p_\lambda = |\langle\lambda|\psi\rangle|^2$, przy czym $|\lambda\rangle$ to stany własne obserwabli $\hat{H}$. Wartość oczekiwana jest zatem zdefiniowana jako średnia wszystkich możliwych wyników $\lambda$ (czyli wartości własnych obserwabli) pomiaru stanu $|\psi\rangle$, ważona odpowiednimi prawdopodobieństwami:

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

Jednak obliczenie wartości oczekiwanej obserwabli nie zawsze jest możliwe, ponieważ często nie znamy jej bazy własnej. Qiskit Runtime Estimator wykorzystuje złożony proces algebraiczny do oszacowania wartości oczekiwanej na rzeczywistym urządzeniu kwantowym, rozkładając obserwablę na kombinację innych obserwabli, których bazę własną znamy.

Mówiąc prościej, Estimator rozkłada każdą obserwablę, której nie wie, jak zmierzyć, na prostsze, mierzalne obserwable zwane operatorami Pauli. Każdy operator można wyrazić jako kombinację $4^n$ operatorów 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\}, \\$$

tak, że

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

gdzie $n$ to liczba kubitów, $k \equiv k_{n-1} \cdots k_0$ dla $k_l \in \mathbb{Z}_4 \equiv \{0, 1, 2, 3\}$ (czyli liczby całkowite w bazie $4$) oraz $(\sigma_0, \sigma_1, \sigma_2, \sigma_3) := (I, X, Y, Z)$.

Po wykonaniu tej dekompozycji, Estimator wyprowadza nowy obwód $V_k|\psi\rangle$ dla każdej obserwabli $\hat{P}_k$ (z oryginalnego obwodu), aby skutecznie zdiagonalizować obserwablę Pauli w bazie obliczeniowej i zmierzyć ją. Możemy łatwo mierzyć obserwable Pauli, ponieważ znamy $V_k$ z góry, co nie jest generalnie prawdziwe dla innych obserwabli.

Dla każdego $\hat{P}_{k}$, Estimator uruchamia odpowiedni obwód na urządzeniu kwantowym wielokrotnie, mierzy stan wyjściowy w bazie obliczeniowej i oblicza prawdopodobieństwo $p_{kj}$ uzyskania każdego możliwego wyjścia $j$. Następnie szuka wartości własnej $\lambda_{kj}$ $P_k$ odpowiadającej każdemu wyjściu $j$, mnoży przez $w_k$ i sumuje wszystkie wyniki, aby otrzymać wartość oczekiwaną obserwabli $\hat{H}$ dla danego stanu $|\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},$$

Ponieważ obliczanie wartości oczekiwanej $4^n$ operatorów Pauli jest niepraktyczne (czyli rośnie wykładniczo), Estimator może być wydajny tylko wtedy, gdy duża liczba $w_k$ wynosi zero (czyli rzadka dekompozycja Pauli zamiast gęstej). Formalnie mówimy, że aby to obliczenie było efektywnie rozwiązywalne, liczba niezerowych składników musi rosnąć co najwyżej wielomianowo wraz z liczbą kubitów $n$: $\hat{H} = \sum^{\text{Poly}(n)}_k w_k \hat{P}_k.$

Czytelnik może zauważyć niejawne założenie, że próbkowanie prawdopodobieństwa również musi być wydajne, jak wyjaśniono dla Sampler, co oznacza

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

Przewodnik na przykładzie obliczania wartości oczekiwanych

Załóżmy stan jednokubitowy $|+\rangle := H|0\rangle = \frac{1}{\sqrt{2}}(|0\rangle + |1\rangle)$ oraz obserwablę

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

o następującej teoretycznej wartości oczekiwanej $\langle\hat{H}\rangle_+ = \langle+|\hat{H}|+\rangle = 2.$

Ponieważ nie wiemy, jak zmierzyć tę obserwablę, nie możemy obliczyć jej wartości oczekiwanej bezpośrednio i musimy wyrazić ją ponownie jako $\langle\hat{H}\rangle_+ = 2\langle X \rangle_+ - \langle Z \rangle_+$. Można pokazać, że daje ona ten sam wynik, zauważając, że $\langle+|X|+\rangle = 1$ oraz $\langle+|Z|+\rangle = 0$.

Zobaczmy, jak bezpośrednio obliczyć $\langle X \rangle_+$ oraz $\langle Z \rangle_+$. Ponieważ $X$ i $Z$ nie komutują (czyli nie dzielą tej samej bazy własnej), nie można ich mierzyć jednocześnie, dlatego potrzebujemy obwodów pomocniczych:

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

Możemy teraz ręcznie przeprowadzić obliczenia za pomocą Sampler i sprawdzić wyniki na 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

Rygor matematyczny (opcjonalnie)

Wyrażając $|\psi\rangle$ względem bazy stanów własnych operatora $\hat{H}$, $|\psi\rangle = \sum_\lambda a_\lambda |\lambda\rangle$, otrzymujemy:

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

Ponieważ nie znamy wartości własnych ani stanów własnych docelowej obserwabli $\hat{H}$, musimy najpierw rozważyć jej diagonalizację. Biorąc pod uwagę, że $\hat{H}$ jest hermitowski, istnieje przekształcenie unitarne $V$ takie, że $\hat{H}=V^\dagger \Lambda V,$ gdzie $\Lambda$ jest diagonalną macierzą wartości własnych, zatem $\langle j | \Lambda | k \rangle = 0$ jeśli $j\neq k$ oraz $\langle j | \Lambda | j \rangle = \lambda_j$.

Oznacza to, że wartość oczekiwaną można przepisać 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}$$

Biorąc pod uwagę, że jeśli układ znajduje się w stanie $|\phi\rangle = V |\psi\rangle$, to prawdopodobieństwo pomiaru $| j\rangle$ wynosi $p_j = |\langle j|\phi \rangle|^2$, powyższą wartość oczekiwaną można wyrazić jako:

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

Bardzo ważne jest, aby zauważyć, że prawdopodobieństwa są brane ze stanu $V |\psi\rangle$, a nie $|\psi\rangle$. Dlatego macierz $V$ jest absolutnie niezbędna. Być może zastanawiasz się, jak uzyskać macierz $V$ i wartości własne $\Lambda$. Gdybyś już miał wartości własne, nie byłoby potrzeby używania komputera kwantowego, ponieważ celem algorytmów wariacyjnych jest znalezienie tych wartości własnych operatora $\hat{H}$.

Na szczęście istnieje sposób na obejście tego: dowolna macierz $2^n \times 2^n$ może być zapisana jako kombinacja liniowa $4^n$ iloczynów tensorowych $n$ macierzy Pauliego i identyczności, z których wszystkie są zarówno hermitowskie, jak i unitarne, ze znanymi $V$ i $\Lambda$. Tak właśnie działa wewnętrznie Estimator Runtime, dekomponując dowolny obiekt Operator na SparsePauliOp.

Oto operatory, których można użyć:

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

Przepiszmy zatem $\hat{H}$ względem macierzy Pauliego i identyczności:

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

gdzie $k = \sum_{l=0}^{n-1} 4^l k_l \equiv k_{n-1}...k_0$ dla $k_{n-1},...,k_0\in \{0,1,2,3\}$ (czyli w bazie $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}$$

gdzie $V_k := V_{k_{n-1}}\otimes ... \otimes V_{k_0}$ oraz $\Lambda_k := \Lambda_{k_{n-1}}\otimes ... \otimes \Lambda_{k_0}$, takie że: $\hat{P_k}=V_k^\dagger \Lambda_k V_k.$

Funkcje kosztu

Ogólnie rzecz biorąc, funkcje kosztu służą do opisania celu problemu oraz tego, jak dobrze stan próbny spełnia ten cel. Tę definicję można stosować w różnych przykładach z chemii, uczenia maszynowego, finansów, optymalizacji itd.

Rozważmy prosty przykład znajdowania stanu podstawowego układu. Naszym celem jest zminimalizowanie wartości oczekiwanej obserwabli reprezentującej energię (Hamiltonian $\hat{\mathcal{H}}$):

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

Możemy użyć Estimator do wyznaczenia wartości oczekiwanej i przekazać tę wartość do optymalizatora w celu minimalizacji. Jeśli optymalizacja zakończy się sukcesem, zwróci ona zestaw optymalnych wartości parametrów $\vec\theta^*$, z których będziemy w stanie skonstruować proponowany stan rozwiązania $|\psi(\vec\theta^*)\rangle$ i obliczyć obserwowaną wartość oczekiwaną jako $C(\vec\theta^*)$.

Zauważ, że będziemy w stanie minimalizować funkcję kosztu jedynie dla ograniczonego zbioru stanów, które rozważamy. Prowadzi to do dwóch odrębnych możliwości:

• Nasz ansatz nie definiuje stanu rozwiązania w całej przestrzeni poszukiwań: Jeśli tak jest, nasz optymalizator nigdy nie znajdzie rozwiązania i musimy eksperymentować z innymi ansatzami, które mogą dokładniej reprezentować naszą przestrzeń poszukiwań.
• Nasz optymalizator nie jest w stanie znaleźć tego prawidłowego rozwiązania: Optymalizacja może być zdefiniowana globalnie i lokalnie. Zgłębimy znaczenie tych pojęć w dalszej części.

Podsumowując, będziemy wykonywać klasyczną pętlę optymalizacyjną, ale polegając na ewaluacji funkcji kosztu na komputerze kwantowym. Z tej perspektywy można myśleć o optymalizacji jako o czysto klasycznym przedsięwzięciu, w którym wywołujemy pewną wyrocznię kwantową typu czarna skrzynka za każdym razem, gdy optymalizator musi obliczyć funkcję kosztu.

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

Najpierw wykonamy to za pomocą symulatora: StatevectorEstimator. Jest to zazwyczaj zalecane do debugowania, ale zaraz po przebiegu debugowania przeprowadzimy obliczenia na prawdziwym sprzęcie kwantowym. Coraz częściej interesujące problemy nie są już klasycznie symulowalne bez najnowocześniejszych obiektów superkomputerowych.

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

[-0.58744589]

Teraz przejdziemy do uruchomienia na prawdziwym komputerze kwantowym. Zwróć uwagę na zmiany w składni. Kroki dotyczące pass_manager zostaną omówione dalej w następnym przykładzie. Jednym ze szczególnie ważnych kroków w algorytmach wariacyjnych jest użycie sesji Qiskit Runtime. Rozpoczęcie sesji umożliwia uruchamianie wielu iteracji algorytmu wariacyjnego bez oczekiwania w nowej kolejce za każdym razem, gdy parametry są aktualizowane. Jest to istotne, jeśli czasy oczekiwania w kolejce są długie i/lub potrzeba wielu iteracji. Tylko partnerzy sieci IBM Quantum® Network mogą korzystać z sesji Runtime. Jeśli nie masz dostępu do sesji, możesz zmniejszyć liczbę iteracji wysyłanych w danym momencie i zapisać najnowsze parametry do wykorzystania w przyszłych przebiegach. Jeśli wyślesz zbyt wiele iteracji lub napotkasz zbyt długie czasy oczekiwania w kolejce, możesz napotkać kod błędu 1217, który odnosi się do długich opóźnień między wysyłaniem zadań.

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

Zauważ, że wartości uzyskane z dwóch powyższych obliczeń są bardzo podobne. Techniki ulepszania wyników zostaną omówione w dalszej części.

Przykład odwzorowania na układy niefizyczne

Problem maksymalnego cięcia (Max-Cut) to problem optymalizacji kombinatorycznej, który polega na podziale wierzchołków grafu na dwa rozłączne zbiory w taki sposób, aby liczba krawędzi między tymi dwoma zbiorami była maksymalna. Bardziej formalnie, mając nieskierowany graf $G=(V,E)$, gdzie $V$ jest zbiorem wierzchołków, a $E$ jest zbiorem krawędzi, problem Max-Cut polega na podziale wierzchołków na dwa rozłączne podzbiory $S$ i $T$ w taki sposób, aby liczba krawędzi, których jeden koniec znajduje się w $S$, a drugi w $T$, była maksymalna.

Możemy zastosować Max-Cut do rozwiązania różnych problemów, w tym: klastrowania, projektowania sieci, przejść fazowych itp. Zaczniemy od utworzenia grafu problemu:

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

Problem ten można wyrazić jako problem optymalizacji binarnej. Dla każdego węzła $0 \leq i < n$, gdzie $n$ jest liczbą węzłów grafu (w tym przypadku $n=4$), rozważymy zmienną binarną $x_i$. Zmienna ta będzie miała wartość $1$, jeśli węzeł $i$ należy do jednej z grup, którą oznaczymy $1$, oraz $0$, jeśli należy do drugiej grupy, którą oznaczymy jako $0$. Oznaczymy również jako $w_{ij}$ (element $(i,j)$ macierzy sąsiedztwa $w$) wagę krawędzi biegnącej od węzła $i$ do węzła $j$. Ponieważ graf jest nieskierowany, $w_{ij}=w_{ji}$. Wówczas nasz problem możemy sformułować jako maksymalizację następującej funkcji kosztu:

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

Aby rozwiązać ten problem za pomocą komputera kwantowego, wyrazimy funkcję kosztu jako wartość oczekiwaną pewnej obserwabli. Jednakże obserwable, które Qiskit natywnie obsługuje, składają się z operatorów Pauli, które mają wartości własne $1$ i $-1$ zamiast $0$ i $1$. Dlatego dokonamy następującej zamiany zmiennych:

Gdzie $\vec{x}=(x_0,x_1,\cdots ,x_{n-1})$. Możemy wykorzystać macierz sąsiedztwa $w$, aby wygodnie uzyskać dostęp do wag wszystkich krawędzi. Zostanie to użyte do uzyskania naszej funkcji kosztu:

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

Oznacza to, że:

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

Zatem nowa funkcja kosztu, którą chcemy zmaksymalizować, to:

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

Co więcej, naturalną tendencją komputera kwantowego jest znajdowanie minimów (zwykle najniższej energii) zamiast maksimów, więc zamiast maksymalizować $C(\vec{z})$, będziemy minimalizować:

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

Teraz, gdy mamy funkcję kosztu do minimalizacji, której zmienne mogą przyjmować wartości $-1$ i $1$, możemy wykonać następującą analogię z Pauli $Z$:

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

Innymi słowy, zmienna $z_i$ będzie odpowiadać bramce $Z$ działającej na kubit $i$. Ponadto:

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

Wtedy obserwablą, którą rozważymy, jest:

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

do której będziemy musieli później dodać niezależny człon:

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

Operator jest liniową kombinacją członów z operatorami Z na węzłach połączonych krawędzią (przypomnijmy, że kubit 0 znajduje się najbardziej na prawo): $IIZZ + IZIZ + IZZI + ZIIZ + ZZII$. Po skonstruowaniu operatora ansatz dla algorytmu QAOA można łatwo zbudować, używając obwodu QAOAAnsatz z biblioteki obwodów 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

Ponieważ Runtime Estimator bezpośrednio przyjmuje Hamiltonian i sparametryzowany ansatz oraz zwraca niezbędną energię, funkcja kosztu dla instancji QAOA jest dość prosta:

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

Wrócimy do tego przykładu w sekcji Zastosowania, aby zbadać, jak wykorzystać optymalizator do iterowania po przestrzeni przeszukiwania. Ogólnie rzecz biorąc, obejmuje to:

• Wykorzystanie optymalizatora do znalezienia optymalnych parametrów
• Powiązanie optymalnych parametrów z ansatzem w celu znalezienia wartości własnych
• Przetłumaczenie wartości własnych na naszą definicję problemu

Strategia pomiaru: szybkość kontra dokładność

Jak wspomniano, używamy zaszumionego komputera kwantowego jako czarnej skrzynki wyroczni, gdzie szum może sprawić, że uzyskiwane wartości będą niedeterministyczne, prowadząc do losowych fluktuacji, które z kolei zaszkodzą — lub nawet całkowicie uniemożliwią — zbieżność pewnych optymalizatorów do proponowanego rozwiązania. Jest to ogólny problem, który musimy rozwiązywać, stopniowo badając użyteczność kwantową i zmierzając w kierunku przewagi kwantowej:

Możemy wykorzystać opcje tłumienia błędów i łagodzenia błędów prymitywów Qiskit Runtime, aby poradzić sobie z szumem i zmaksymalizować użyteczność dzisiejszych komputerów kwantowych.

Tłumienie błędów

Tłumienie błędów odnosi się do technik stosowanych do optymalizacji i transformacji obwodu podczas kompilacji w celu zminimalizowania błędów. Jest to podstawowa technika obsługi błędów, która zwykle skutkuje pewnym klasycznym narzutem przetwarzania wstępnego do całkowitego czasu wykonania. Narzut ten obejmuje transpilację obwodów do uruchomienia na sprzęcie kwantowym poprzez:

• Wyrażenie obwodu za pomocą natywnych bramek dostępnych w systemie kwantowym
• Odwzorowanie wirtualnych kubitów na fizyczne kubity
• Dodawanie operacji SWAP na podstawie wymagań łączności
• Optymalizację bramek 1Q i 2Q
• Dodawanie dynamicznego rozprzęgania do bezczynnych kubitów, aby zapobiec efektom dekoherencji.

Prymitywy umożliwiają korzystanie z technik tłumienia błędów poprzez ustawienie opcji optimization_level i wybór zaawansowanych opcji transpilacji. W późniejszym kursie zagłębimy się w różne metody konstrukcji obwodów w celu poprawy wyników, ale w większości przypadków zalecamy ustawienie optimization_level=3.

Zwizualizujemy wartość zwiększenia optymalizacji w procesie transpilacji, patrząc na przykładowy obwód z prostym idealnym zachowaniem.

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

Powyższy obwód może dawać sinusoidalne wartości oczekiwane danej obserwabli, pod warunkiem że wstawimy fazy obejmujące odpowiedni przedział, taki jak $[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]

Możemy użyć symulatora, aby pokazać użyteczność zoptymalizowanej transpilacji. Poniżej powrócimy do korzystania z prawdziwego sprzętu, aby zademonstrować użyteczność łagodzenia błędów. Użyjemy QiskitRuntimeService, aby uzyskać prawdziwy backend (w tym przypadku ibm_brisbane), oraz AerSimulator, aby symulować ten backend, w tym jego zachowanie związane z szumem.

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)

Możemy teraz użyć menedżera przebiegów, aby transpilować obwód do "architektury zestawu instrukcji" lub ISA backendu. Jest to nowe wymaganie w Qiskit Runtime: wszystkie obwody przesyłane do backendu muszą być zgodne z ograniczeniami celu backendu, co oznacza, że muszą być napisane w kategoriach ISA backendu — czyli zestawu instrukcji, które urządzenie może zrozumieć i wykonać. Te ograniczenia celu są zdefiniowane przez takie czynniki, jak natywne bramki bazowe urządzenia, jego łączność kubitów oraz — gdy ma to znaczenie — specyfikacje impulsów i innych czasów instrukcji.

Zauważ, że w obecnym przypadku zrobimy to dwukrotnie: raz z optimization_level = 0, a raz ustawione na 3. Za każdym razem użyjemy prymitywu Estimator do oszacowania wartości oczekiwanych obserwabli przy różnych wartościach fazy.

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

Na koniec możemy wykreślić wyniki i widzimy, że precyzja obliczeń była dość dobra nawet bez optymalizacji, ale zdecydowanie poprawiła się po zwiększeniu optymalizacji do poziomu 3. Zauważ, że w głębszych, bardziej skomplikowanych obwodach różnica między poziomami optymalizacji 0 i 3 będzie prawdopodobnie bardziej znacząca. Jest to bardzo prosty obwód użyty jako zabawkowy 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()

Łagodzenie błędów

Łagodzenie błędów odnosi się do technik, które pozwalają użytkownikom redukować błędy obwodu poprzez modelowanie szumu urządzenia w czasie wykonania. Zazwyczaj skutkuje to kwantowym narzutem przetwarzania wstępnego związanym z trenowaniem modelu oraz klasycznym narzutem przetwarzania końcowego w celu złagodzenia błędów w surowych wynikach za pomocą wygenerowanego modelu.

Opcja resilience_level prymitywu Qiskit Runtime określa poziom odporności budowanej przeciwko błędom. Wyższe poziomy generują dokładniejsze wyniki kosztem dłuższych czasów przetwarzania z powodu narzutu próbkowania kwantowego. Poziomy odporności mogą być używane do konfiguracji kompromisu między kosztem a dokładnością podczas stosowania łagodzenia błędów w zapytaniach prymitywów.

Wdrażając jakąkolwiek technikę łagodzenia błędów, oczekujemy, że obciążenie w naszych wynikach zostanie zmniejszone w stosunku do poprzedniego, niezłagodzonego obciążenia. W niektórych przypadkach obciążenie może nawet zniknąć. Wiąże się to jednak z kosztem. Gdy zmniejszamy obciążenie w naszych szacowanych wielkościach, zmienność statystyczna wzrasta (czyli wariancja), co możemy uwzględnić, dalej zwiększając liczbę pomiarów na obwód w naszym procesie próbkowania. Wprowadzi to narzut poza tym potrzebnym do zmniejszenia obciążenia, więc nie jest to robione domyślnie. Możemy łatwo włączyć to zachowanie, dostosowując liczbę pomiarów na obwód w options.executions.shots, jak pokazano w przykładzie poniżej.

W tym kursie zbadamy te modele łagodzenia błędów na wysokim poziomie, aby zilustrować łagodzenie błędów, które prymitywy Qiskit Runtime mogą wykonać bez konieczności znajomości pełnych szczegółów implementacji.

Twirled readout error extinction (T-REx)

Twirled readout error extinction (T-REx) wykorzystuje technikę znaną jako twirling Pauliego w celu zmniejszenia szumu wprowadzanego podczas procesu pomiaru kwantowego. Technika ta nie zakłada konkretnej postaci szumu, co czyni ją bardzo ogólną i skuteczną.

Ogólny przepływ pracy:

1. Pozyskaj dane dla stanu zerowego z losowymi odwróceniami bitów (Pauli X przed pomiarem)
2. Pozyskaj dane dla pożądanego (zaszumionego) stanu z losowymi odwróceniami bitów (Pauli X przed pomiarem)
3. Oblicz specjalną funkcję dla każdego zestawu danych i podziel.  

Możemy ustawić to za pomocą options.resilience_level = 1, co pokazano w przykładzie poniżej.

Ekstrapolacja do zerowego szumu

Ekstrapolacja do zerowego szumu (ZNE, zero noise extrapolation) działa poprzez wstępne wzmocnienie szumu w obwodzie przygotowującym żądany stan kwantowy, uzyskanie pomiarów dla kilku różnych poziomów szumu, a następnie wykorzystanie tych pomiarów do wnioskowania o wyniku pozbawionym szumu.

Ogólny przebieg pracy:

1. Wzmocnij szum obwodu dla kilku współczynników szumu
2. Uruchom każdy obwód ze wzmocnionym szumem
3. Ekstrapoluj z powrotem do granicy zerowego szumu

 

Możemy to ustawić za pomocą options.resilience_level = 2. Możemy dodatkowo zoptymalizować to, eksplorując różne noise_factors, noise_amplifiers i extrapolators, ale wykracza to poza zakres tego kursu. Zachęcamy do eksperymentowania z tymi opcjami opisanymi tutaj.

Każda metoda wiąże się z własnym narzutem: kompromisem między liczbą wymaganych obliczeń kwantowych (czas) a dokładnością naszych wyników:

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

Korzystanie z opcji mitigacji i tłumienia w Qiskit Runtime

Oto jak obliczyć wartość oczekiwaną przy jednoczesnym użyciu mitigacji i tłumienia błędów w Qiskit Runtime. Możemy wykorzystać dokładnie ten sam obwód i obserwablę co poprzednio, ale tym razem utrzymując poziom optymalizacji na poziomie 2, a teraz dostrajając odporność lub stosowane techniki mitigacji błędów. Ten proces mitigacji błędów odbywa się wielokrotnie w całej pętli optymalizacji.

Tę część wykonujemy na rzeczywistym sprzęcie, ponieważ mitigacja błędów nie jest dostępna w symulatorach.

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

Tak jak poprzednio, możemy wykreślić otrzymane wartości oczekiwane jako funkcję kąta fazowego dla trzech zastosowanych poziomów mitigacji błędów. Z dużym trudem można dostrzec, że mitigacja błędów nieznacznie poprawia wyniki. Ponownie, efekt ten jest znacznie bardziej wyraźny w głębszych, bardziej skomplikowanych obwodach.

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

Podsumowanie

Dzięki tej lekcji nauczyłeś się, jak stworzyć funkcję kosztu:

• Tworzenie funkcji kosztu
• Jak wykorzystać prymitywy Qiskit Runtime do mitigacji i tłumienia szumu
• Jak zdefiniować strategię pomiaru w celu optymalizacji kompromisu między szybkością a dokładnością

Oto nasze wysokopoziomowe obciążenie wariacyjne:

Nasza funkcja kosztu jest uruchamiana podczas każdej iteracji pętli optymalizacyjnej. Następna lekcja będzie poświęcona temu, jak klasyczny optymalizator wykorzystuje naszą ewaluację funkcji kosztu do wybierania nowych parametrów.

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