Функції витрат

У цьому уроці ми навчимося обчислювати функцію витрат:

• Спочатку ознайомимося з примітивами Qiskit Runtime
• Визначимо функцію витрат $C(\vec\theta)$ — специфічну для задачі функцію, яка описує її мету та яку оптимізатор має мінімізувати (або максимізувати)
• Визначимо стратегію вимірювань за допомогою примітивів Qiskit Runtime для балансування між швидкістю та точністю

 

Примітиви

Усі фізичні системи — як класичні, так і квантові — можуть перебувати в різних станах. Наприклад, автомобіль на дорозі характеризується певною масою, положенням, швидкістю або прискоренням. Квантові системи так само можуть мати різні конфігурації або стани, але відрізняються від класичних тим, як ми маємо справу з вимірюваннями і еволюцією стану. Це призводить до унікальних властивостей — суперпозиції та заплутаності, — властивих виключно квантовій механіці. Подібно до того, як стан автомобіля можна описати фізичними величинами на кшталт швидкості чи прискорення, стан квантової системи можна описати за допомогою спостережуваних — математичних об'єктів.

У квантовій механіці стани представляють нормованими комплексними стовпчастими векторами, або кетами ($|\psi\rangle$), а спостережувані — ермітовими лінійними операторами ($\hat{H}=\hat{H}^{\dagger}$), що діють на кети. Власний вектор ($|\lambda\rangle$) спостережуваної називають власним станом. Вимірювання спостережуваної для одного з її власних станів ($|\lambda\rangle$) дасть відповідне власне значення ($\lambda$) як результат.

Якщо тебе цікавить, як вимірювати квантову систему і що саме можна вимірювати, Qiskit пропонує два примітиви:

• Sampler: для заданого квантового стану $|\psi\rangle$ цей примітив обчислює ймовірність кожного можливого стану обчислювального базису.
• Estimator: для заданої квантової спостережуваної $\hat{H}$ і стану $|\psi\rangle$ цей примітив обчислює очікуване значення $\hat{H}$.

Примітив Sampler

Примітив Sampler обчислює ймовірність отримання кожного можливого стану $|k\rangle$ обчислювального базису для квантового Circuit, що готує стан $|\psi\rangle$. Він обчислює

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

де $n$ — кількість Qubit, а $k$ — цілочисельне представлення будь-якого можливого вихідного бінарного рядка $\{0,1\}^n$ (тобто цілі числа за основою $2$).

Qiskit Runtime Sampler запускає Circuit кілька разів на квантовому пристрої, виконуючи вимірювання під час кожного запуску, та відновлює розподіл ймовірностей із отриманих бітових рядків. Чим більше запусків (або шотів) виконується, тим точнішими будуть результати, але це вимагає більше часу та квантових ресурсів.

Однак, оскільки кількість можливих виходів зростає експоненційно з кількістю Qubit $n$ (тобто $2^n$), кількість шотів також має зростати експоненційно, щоб охопити щільний розподіл ймовірностей. Тому Sampler ефективний лише для розріджених розподілів ймовірностей; при цьому цільовий стан $|\psi\rangle$ має виражатися як лінійна комбінація станів обчислювального базису, де кількість доданків зростає не швидше, ніж поліноміально від кількості Qubit:

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

Sampler також можна налаштувати так, щоб отримувати ймовірності з підсхеми Circuit, що відповідає підмножині всіх можливих станів.

Примітив Estimator

Примітив Estimator обчислює очікуване значення спостережуваної $\hat{H}$ для квантового стану $|\psi\rangle$; при цьому ймовірності спостережуваної можна виразити як $p_\lambda = |\langle\lambda|\psi\rangle|^2$, де $|\lambda\rangle$ — власні стани спостережуваної $\hat{H}$. Очікуване значення визначається як середнє всіх можливих результатів $\lambda$ (тобто власних значень спостережуваної) вимірювання стану $|\psi\rangle$, зважене відповідними ймовірностями:

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

Проте обчислення очікуваного значення спостережуваної не завжди можливе, оскільки ми часто не знаємо її власного базису. Qiskit Runtime Estimator використовує складний алгебраїчний процес для оцінки очікуваного значення на реальному квантовому пристрої, розкладаючи спостережувану на комбінацію інших спостережуваних, чий власний базис нам відомий.

Простіше кажучи, Estimator розкладає будь-яку спостережувану, яку він не знає як виміряти, на простіші вимірювані спостережувані — оператори Паулі. Будь-який оператор можна виразити як комбінацію $4^n$ операторів Паулі.

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

таким чином

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

де $n$ — кількість Qubit, $k \equiv k_{n-1} \cdots k_0$ при $k_l \in \mathbb{Z}_4 \equiv \{0, 1, 2, 3\}$ (тобто цілі числа за основою $4$), а $(\sigma_0, \sigma_1, \sigma_2, \sigma_3) := (I, X, Y, Z)$.

Після такого розкладання Estimator будує новий Circuit $V_k|\psi\rangle$ для кожної спостережуваної $\hat{P}_k$ (на основі вихідного Circuit), щоб ефективно діагоналізувати оператор Паулі в обчислювальному базисі і виміряти його. Ми можемо легко вимірювати оператори Паулі, оскільки знаємо $V_k$ заздалегідь — що, як правило, не виконується для довільних спостережуваних.

Для кожного $\hat{P}_{k}$ Estimator запускає відповідний Circuit на квантовому пристрої кілька разів, вимірює вихідний стан у обчислювальному базисі та обчислює ймовірність $p_{kj}$ отримання кожного можливого виходу $j$. Далі він знаходить власне значення $\lambda_{kj}$ оператора $P_k$, що відповідає кожному виходу $j$, множить на $w_k$ і підсумовує всі результати, щоб отримати очікуване значення спостережуваної $\hat{H}$ для стану $|\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},$$

Оскільки обчислення очікуваного значення для $4^n$ операторів Паулі є практично нездійсненним (тобто зростає експоненційно), Estimator може бути ефективним лише тоді, коли велика кількість $w_k$ дорівнює нулю (тобто використовується розріджений розклад за Паулі, а не щільний). Формально кажемо, що для того, щоб це обчислення було ефективно розв'язним, кількість ненульових доданків має зростати не швидше ніж поліноміально від кількості Qubit $n$: $\hat{H} = \sum^{\text{Poly}(n)}_k w_k \hat{P}_k.$

Уважний читач помітить неявне припущення про те, що вибірка ймовірностей також має бути ефективною — як пояснено для Sampler, — а отже

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

Покроковий приклад обчислення очікуваних значень

Припустимо, що маємо однокубітний стан $|+\rangle := H|0\rangle = \frac{1}{\sqrt{2}}(|0\rangle + |1\rangle)$ та спостережувану

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

з теоретичним очікуваним значенням $\langle\hat{H}\rangle_+ = \langle+|\hat{H}|+\rangle = 2.$

Оскільки ми не знаємо, як виміряти цю спостережувану безпосередньо, неможливо обчислити її очікуване значення прямим способом, і нам потрібно переписати його як $\langle\hat{H}\rangle_+ = 2\langle X \rangle_+ - \langle Z \rangle_+$. Можна показати, що це дає той самий результат, враховуючи, що $\langle+|X|+\rangle = 1$ та $\langle+|Z|+\rangle = 0$.

Подивимось, як обчислити $\langle X \rangle_+$ та $\langle Z \rangle_+$ безпосередньо. Оскільки $X$ і $Z$ не комутують (тобто не мають спільного власного базису), їх не можна вимірювати одночасно, тому потрібні допоміжні Circuit-и:

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

Тепер можемо виконати обчислення вручну за допомогою Sampler і перевірити результати через 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

Математична строгість (необов'язково)

Виражаючи $|\psi\rangle$ через базис власних станів $\hat{H}$, а саме $|\psi\rangle = \sum_\lambda a_\lambda |\lambda\rangle$, отримуємо:

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

Оскільки ми не знаємо власних значень чи власних станів цільової спостережуваної $\hat{H}$, спочатку потрібно розглянути її діагоналізацію. З огляду на те, що $\hat{H}$ є ермітовою, існує унітарне перетворення $V$ таке, що $\hat{H}=V^\dagger \Lambda V,$ де $\Lambda$ — діагональна матриця власних значень: $\langle j | \Lambda | k \rangle = 0$ при $j\neq k$ та $\langle j | \Lambda | j \rangle = \lambda_j$.

Це означає, що очікуване значення можна переписати як:

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

Враховуючи, що якщо система перебуває в стані $|\phi\rangle = V |\psi\rangle$, ймовірність виміряти $| j\rangle$ дорівнює $p_j = |\langle j|\phi \rangle|^2$, наведене очікуване значення можна виразити як:

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

Важливо зазначити, що ймовірності беруться зі стану $V |\psi\rangle$, а не зі стану $|\psi\rangle$. Саме тому матриця $V$ є абсолютно необхідною. Можливо, тебе цікавить, як отримати матрицю $V$ і власні значення $\Lambda$. Якби власні значення вже були відомі, квантовий комп'ютер взагалі не знадобився б — адже мета варіаційних алгоритмів і полягає у пошуку цих власних значень $\hat{H}$.

На щастя, існує вихід: будь-яку матрицю розміру $2^n \times 2^n$ можна записати як лінійну комбінацію $4^n$ тензорних добутків із $n$ матриць Паулі та одиничних матриць, кожен з яких є одночасно ермітовим і унітарним із відомими $V$ і $\Lambda$. Саме це Runtime Estimator робить внутрішньо, розкладаючи будь-який об'єкт Operator у SparsePauliOp.

Ось оператори, які можна використовувати:

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

Отже, перепишемо $\hat{H}$ через матриці Паулі та одиничні матриці:

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

де $k = \sum_{l=0}^{n-1} 4^l k_l \equiv k_{n-1}...k_0$ при $k_{n-1},...,k_0\in \{0,1,2,3\}$ (тобто за основою $4$), і $\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}$$

де $V_k := V_{k_{n-1}}\otimes ... \otimes V_{k_0}$ і $\Lambda_k := \Lambda_{k_{n-1}}\otimes ... \otimes \Lambda_{k_0}$, причому: $\hat{P_k}=V_k^\dagger \Lambda_k V_k.$

Функції витрат

Загалом функції витрат використовуються для опису мети задачі та оцінки того, наскільки добре пробний стан відповідає цій меті. Це визначення застосовне в різних галузях — хімії, машинному навчанні, фінансах, оптимізації тощо.

Розглянемо простий приклад — пошук основного стану системи. Наша мета — мінімізувати очікуване значення спостережуваної, що представляє енергію (Гамільтоніан $\hat{\mathcal{H}}$):

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

Ми можемо використовувати Estimator для обчислення очікуваного значення та передавати його оптимізатору для мінімізації. Якщо оптимізація пройде успішно, вона поверне набір оптимальних значень параметрів $\vec\theta^*$, на основі яких ми побудуємо запропонований розв'язок $|\psi(\vec\theta^*)\rangle$ і обчислимо спостережуване очікуване значення $C(\vec\theta^*)$.

Зверни увагу: ми зможемо мінімізувати функцію витрат лише в межах тих станів, які ми розглядаємо. Це призводить до двох окремих можливостей:

• Наш ансатц не охоплює розв'язуючий стан у просторі пошуку: у цьому разі оптимізатор ніколи не знайде розв'язок, і нам потрібно спробувати інші ансатци, здатні точніше представити простір пошуку.
• Оптимізатор не може знайти цей допустимий розв'язок: оптимізація може бути глобальною або локальною. Що це означає, ми розглянемо далі.

Загалом ми виконуємо класичний цикл оптимізації, але покладаємось на квантовий комп'ютер для обчислення функції витрат. З цієї точки зору оптимізацію можна сприймати як суто класичну задачу, де щоразу, коли оптимізатору потрібно обчислити функцію витрат, він звертається до певного чорноящикового квантового оракула.

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

Спочатку виконаємо це за допомогою симулятора: StatevectorEstimator. Зазвичай це рекомендується для налагодження, але одразу після налагодження ми перейдемо до розрахунку на реальному квантовому залізі. Дедалі більше задач, що становлять практичний інтерес, вже неможливо класично моделювати без надпотужних суперкомп'ютерів.

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

[-0.58744589]

Тепер перейдемо до запуску на реальному квантовому комп'ютері. Зверни увагу на зміни синтаксису. Кроки, пов'язані з pass_manager, будуть детально розглянуті в наступному прикладі. Особливо важливим кроком у варіаційних алгоритмах є використання сесії Qiskit Runtime. Відкриття сесії дозволяє виконувати кілька ітерацій варіаційного алгоритму без очікування в новій черзі щоразу, коли оновлюються параметри. Це важливо, якщо черги довгі та/або потрібна велика кількість ітерацій. Лише партнери мережі IBM Quantum® Network можуть використовувати сесії Runtime. Якщо у тебе немає доступу до сесій, можна зменшити кількість ітерацій, що подаються за один раз, та зберігати найновіші параметри для подальших запусків. Якщо подати надто багато ітерацій або стикнутися з надто довгими чергами, може виникнути код помилки 1217, що вказує на тривалі затримки між поданнями завдань.

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

Зверни увагу, що значення, отримані з двох наведених вище розрахунків, дуже схожі. Техніки покращення результатів будуть розглянуті нижче.

Приклад відображення на нефізичні системи

Задача максимального розрізу (Max-Cut) — це задача комбінаторної оптимізації, яка полягає в поділі вершин графа на дві непересічні множини таким чином, щоб кількість ребер між цими двома множинами була максимальною. Формально: дано неорієнтований граф $G=(V,E)$, де $V$ — множина вершин, а $E$ — множина ребер; задача Max-Cut полягає в розбитті вершин на дві непересічні підмножини $S$ і $T$ так, щоб кількість ребер з одним кінцем у $S$ та іншим у $T$ була максимальною.

Max-Cut можна застосовувати для розв'язання різних задач, зокрема: кластеризації, проєктування мереж, фазових переходів тощо. Почнемо зі створення графа задачі:

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

Цю задачу можна виразити як задачу бінарної оптимізації. Для кожного вузла $0 \leq i < n$, де $n$ — кількість вузлів графа (у цьому випадку $n=4$), введемо бінарну змінну $x_i$. Ця змінна матиме значення $1$, якщо вузол $i$ входить до групи, яку ми позначимо як $1$, і $0$, якщо він входить до іншої групи, позначеної як $0$. Також позначимо через $w_{ij}$ (елемент $(i,j)$ матриці суміжності $w$) вагу ребра між вузлами $i$ та $j$. Оскільки граф неорієнтований, $w_{ij}=w_{ji}$. Тоді задачу можна сформулювати як максимізацію такої функції вартості:

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

Щоб розв'язати цю задачу на квантовому комп'ютері, ми виразимо функцію вартості як математичне сподівання спостережуваної величини. Однак спостережувані, що їх Qiskit підтримує нативно, складаються з операторів Паулі, власні значення яких дорівнюють $1$ та $-1$, а не $0$ та $1$. Тому зробимо таку заміну змінних:

Де $\vec{x}=(x_0,x_1,\cdots ,x_{n-1})$. Для зручного доступу до ваг усіх ребер скористаємося матрицею суміжності $w$. Це дозволить отримати нашу функцію вартості:

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

Це означає, що:

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

Отже, нова функція вартості, яку потрібно максимізувати, матиме вигляд:

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

Крім того, квантовий комп'ютер за своєю природою прагне знаходити мінімуми (зазвичай найнижчу енергію), а не максимуми, тому замість максимізації $C(\vec{z})$ будемо мінімізувати:

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

Тепер, коли ми маємо функцію вартості для мінімізації, змінні якої можуть приймати значення $-1$ і $1$, можна провести аналогію з оператором Паулі $Z$:

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

Іншими словами, змінна $z_i$ буде еквівалентна Gate $Z$, що діє на Qubit $i$. Крім того:

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

Тоді спостережувана, яку ми розглядатимемо, матиме вигляд:

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

до якої після цього необхідно додати незалежний доданок:

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

Оператор є лінійною комбінацією доданків із операторами Z на вузлах, з'єднаних ребром (нагадаємо, що нульовий Qubit розташований крайнім праворуч): $IIZZ + IZIZ + IZZI + ZIIZ + ZZII$. Після побудови оператора ансац для алгоритму QAOA легко сконструювати за допомогою Circuit QAOAAnsatz з бібліотеки схем 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

Оскільки Runtime Estimator безпосередньо приймає гамільтоніан і параметризований ансац та повертає необхідну енергію, функція вартості для екземпляра QAOA є досить простою:

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

До цього прикладу ми повернемося в розділі «Застосування», щоб розглянути, як використовувати оптимізатор для перебору простору пошуку. Загалом це включає:

• Використання оптимізатора для знаходження оптимальних параметрів
• Підстановку оптимальних параметрів в ансац для знаходження власних значень
• Перетворення власних значень відповідно до постановки нашої задачі

Стратегія вимірювання: швидкість проти точності

Як уже зазначалося, ми використовуємо зашумлений квантовий комп'ютер як оракул «чорного ящика», де шум може робити отримані значення недетермінованими, призводячи до випадкових флуктуацій, які, у свою чергу, заважають — або навіть повністю перешкоджають — збіжності певних оптимізаторів до запропонованого розв'язку. Це загальна проблема, яку необхідно вирішувати в міру поступового дослідження квантової корисності та просування до квантової переваги:

Для усунення шумів і максимального використання можливостей сучасних квантових комп'ютерів можна скористатися опціями придушення та пом'якшення помилок у Qiskit Runtime Primitives.

Придушення помилок

Придушення помилок — це методи оптимізації та перетворення Circuit під час компіляції з метою мінімізації помилок. Це базова техніка обробки помилок, яка зазвичай призводить до певних класичних витрат на попередню обробку overhead у загальному часі виконання. Ці витрати включають транспіляцію схем для запуску на квантовому залізі шляхом:

• Вираження Circuit через нативні Gate, доступні на квантовій системі
• Відображення віртуальних Qubit на фізичні Qubit
• Додавання операцій SWAP відповідно до вимог зв'язності
• Оптимізації Gate для одного та двох Qubit
• Додавання динамічного розв'язання до простоюючих Qubit для запобігання ефектам декогеренції

Примітиви дозволяють застосовувати техніки придушення помилок шляхом встановлення опції optimization_level та вибору розширених параметрів транспіляції. У наступному курсі ми детально розглянемо різні методи побудови Circuit для покращення результатів, але в більшості випадків рекомендуємо встановлювати optimization_level=3.

Ми наочно продемонструємо цінність підвищення рівня оптимізації в процесі транспіляції на прикладі Circuit з простою ідеальною поведінкою.

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

Наведена вище Circuit може давати синусоїдальні математичні сподівання вказаної спостережуваної за умови підстановки фаз у відповідному діапазоні, наприклад $[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]

Для демонстрації корисності оптимізованої транспіляції скористаємося симулятором. Нижче ми повернемося до використання реального залізка для демонстрації ефективності пом'якшення помилок. Ми використаємо QiskitRuntimeService, щоб отримати реальний Backend (у цьому випадку ibm_brisbane), і AerSimulator для його симуляції, включно з поведінкою шумів.

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)

Тепер можна скористатися менеджером проходів для транспіляції Circuit у «архітектуру набору інструкцій» (ISA) Backend. Це нова вимога в Qiskit Runtime: усі Circuit, що подаються на Backend, мають відповідати обмеженням цільового пристрою, тобто бути написані в термінах ISA Backend — набору інструкцій, які пристрій може зрозуміти та виконати. Ці цільові обмеження визначаються такими факторами, як нативні базові Gate пристрою, його зв'язність Qubit, а також — коли це актуально — специфікації синхронізації імпульсів та інших інструкцій.

Зверни увагу, що в цьому випадку ми зробимо це двічі: один раз із optimization_level = 0, і один раз зі значенням 3. Кожного разу ми використовуватимемо примітив Estimator для оцінки математичних сподівань спостережуваної при різних значеннях фази.

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

Нарешті, ми можемо побудувати графіки результатів і побачимо, що точність обчислень була досить хорошою навіть без оптимізації, але помітно покращилась після підвищення рівня оптимізації до 3. Зверни увагу, що для глибших та складніших Circuit різниця між рівнями оптимізації 0 та 3, швидше за все, буде значно більш суттєвою. Ця Circuit є дуже простою і використовується лише як навчальна модель.

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

Пом'якшення помилок

Пом'якшення помилок — це методи, які дозволяють користувачам зменшити помилки Circuit шляхом моделювання шумів пристрою під час виконання. Зазвичай це призводить до витрат на квантову попередню обробку, пов'язаних із навчанням моделі, та на класичну постобробку для пом'якшення помилок у необроблених результатах за допомогою побудованої моделі.

Параметр resilience_level примітива Qiskit Runtime визначає ступінь стійкості до помилок. Вищі рівні дають точніші результати ціною збільшення часу обробки через накладні витрати на квантову вибірку. Рівні стійкості можна використовувати для налаштування балансу між вартістю та точністю при застосуванні пом'якшення помилок до запиту примітива.

При застосуванні будь-якої техніки пом'якшення помилок ми очікуємо, що зміщення у наших результатах зменшиться порівняно з попередніми, непом'якшеними результатами. У деяких випадках зміщення може навіть зникнути повністю. Проте за це доводиться платити: зменшуючи зміщення в оцінюваних величинах, ми збільшуємо статистичну варіативність (тобто дисперсію), яку можна компенсувати, збільшивши кількість вимірювань на Circuit у процесі вибірки. Це вносить додаткові витрати понад ті, що необхідні для зменшення зміщення, тому за замовчуванням цього не робиться. Ми можемо легко ввімкнути таку поведінку, відрегулювавши кількість вимірювань на Circuit через options.executions.shots, як показано в прикладі нижче.

У цьому курсі ми розглянемо моделі пом'якшення помилок на високому рівні, щоб проілюструвати можливості пом'якшення помилок у примітивах Qiskit Runtime без заглиблення в деталі реалізації.

Згасання помилок зчитування за допомогою скручування (T-REx)

Згасання помилок зчитування за допомогою скручування (T-REx) використовує техніку, відому як скручування Паулі, для зменшення шуму, що виникає під час квантового вимірювання. Ця техніка не передбачає жодної конкретної форми шуму, що робить її дуже загальною та ефективною.

Загальний робочий процес:

1. Отримати дані для нульового стану з рандомізованими бітовими інверсіями (Pauli X перед вимірюванням)
2. Отримати дані для бажаного (зашумленого) стану з рандомізованими бітовими інверсіями (Pauli X перед вимірюванням)
3. Обчислити спеціальну функцію для кожного набору даних і поділити результати.

 

Це можна задати через options.resilience_level = 1, що показано в прикладі нижче.

Екстраполяція нульового шуму

Екстраполяція нульового шуму (ZNE) працює шляхом спочатку посилення шуму в Circuit, що готує бажаний квантовий стан, отримання вимірювань для кількох різних рівнів шуму, а потім використання цих вимірювань для висновку про беззашумний результат.

Загальний робочий процес:

1. Посилити шум Circuit для кількох коефіцієнтів шуму
2. Запустити кожну схему з посиленим шумом
3. Екстраполювати назад до межі нульового шуму

 

Це можна задати через options.resilience_level = 2. Можна додатково оптимізувати, досліджуючи різні noise_factors, noise_amplifiers та extrapolators, але це виходить за рамки цього курсу. Пропонуємо тобі поекспериментувати з цими параметрами, описаними тут.

Кожен метод пов'язаний із власними накладними витратами — компромісом між кількістю необхідних квантових обчислень (часом) та точністю результатів:

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

Використання опцій пом'якшення та придушення помилок у Qiskit Runtime

Ось як обчислити математичне сподівання з використанням пом'якшення та придушення помилок у Qiskit Runtime. Ми можемо використати ту саму Circuit і ту саму спостережувану, що й раніше, але цього разу зафіксуємо рівень оптимізації на значенні 2 і натомість будемо змінювати рівень стійкості або техніку(-и) пом'якшення помилок. Цей процес пом'якшення помилок відбувається кілька разів впродовж циклу оптимізації.

Цю частину ми виконуємо на реальному залізі, оскільки пом'якшення помилок недоступне на симуляторах.

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

Як і раніше, ми можемо побудувати графіки отриманих математичних сподівань як функцію кута фази для трьох рівнів пом'якшення помилок. Можна побачити, що пом'якшення помилок дещо покращує результати. Знову ж таки, цей ефект значно помітніший у глибших та складніших 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()

Підсумок

У цьому уроці ти навчився(-лась) створювати функцію вартості:

• Створювати функцію вартості
• Використовувати примітиви Qiskit Runtime для пом'якшення та придушення шумів
• Визначати стратегію вимірювання для оптимізації балансу між швидкістю та точністю

Ось наше варіаційне навантаження на високому рівні:

Наша функція вартості виконується на кожній ітерації циклу оптимізації. У наступному уроці ми розглянемо, як класичний оптимізатор використовує обчислене значення функції вартості для вибору нових параметрів.

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