Багатопродуктні формули для зменшення похибки Троттера

Очікуване використання: чотири хвилини на процесорі Heron r2 (ПРИМІТКА: Це лише оцінка. Твій час виконання може відрізнятися.)

Результати навчання

• Як мультипродуктні формули (MPF) зменшують похибку Троттера при моделюванні гамільтоніана, комбінуючи очікувані значення з кількох неглибоких схем

• Коли MPF корисні порівняно зі стандартними продуктними формулами, а коли вони не є правильним інструментом

• Як обчислити статичні та динамічні коефіцієнти MPF за допомогою пакета qiskit_addon_mpf

• Як виконати робочий процес MPF від початку до кінця на апаратному забезпеченні IBM Quantum®, включно з транспіляцією, зменшенням похибок та постобробкою

Передумови

• Методи компіляції для схем моделювання гамільтоніана — знайомить з продуктними схемами Троттера в Qiskit.

• Продуктні формули в Qiskit, зокрема класи синтезу SuzukiTrotter та LieTrotter.

• Примітиви Qiskit та інтерфейс Estimator.

Передумови

Що таке мультипродуктні формули?

При моделюванні квантових систем на квантовому комп'ютері центральним завданням є наближення оператора часової еволюції $e^{-iHt}$ для гамільтоніана $H$. Стандартний підхід використовує продуктні формули (PF), також відомі як розкладання Троттера-Судзукі. Вони розкладають $H = \sum_{a=1}^d F_a$ на члени, окремі унітарні оператори яких $e^{-iF_a t}$ ефективно реалізувати, а потім наближають повну еволюцію як впорядкований добуток цих простіших унітарних операторів.

Продуктна формула першого порядку (Лі-Троттера) має вигляд:

$$S_1(t) := \prod_{a=1}^d e^{-i F_a t},$$

яка вносить квадратичну похибку: $S_1(t) = e^{-iHt} + \mathcal{O}(t^2)$. Симетричні формули вищого порядку $S_{2\chi}(t)$, де $\chi$ позначає порядок симетричної продуктної формули (див. посилання [1]), збігаються швидше як $e^{-iHt} + \mathcal{O}(t^{2\chi+1})$, але за рахунок глибших схем на крок.

Щоб зменшити похибку при фіксованому порядку $\chi$, зазвичай розбивають загальний час еволюції $t$ на $k$ менших кроків Троттера. Кожен крок наближає $e^{-iHt/k}$ продуктною формулою, і кроки об'єднуються:

$$e^{-iHt} \approx \left[S_{2\chi}(t/k)\right]^k.$$

Для симетричної формули порядку $2\chi$ залишкова похибка Троттера тоді масштабується як $\mathcal{O}\!\left(t^{2\chi+1} / k^{2\chi}\right)$. Тож збільшення $k$ швидко пригнічує похибку Троттера — але це також лінійно поглиблює схему, а на зашумленому апаратному забезпеченні це означає більше накопиченого шуму вентилів. Це протиріччя між похибкою Троттера (сприяє більшому $k$) і апаратним шумом (сприяє меншому $k$) — це саме те, для вирішення чого призначені мультипродуктні формули. Зауваж, що MPF стосуються комбінування результатів з різних виборів $k$ при фіксованому порядку $\chi$ — вони не змінюють порядок базової продуктної формули.

Мультипродуктні формули (MPF) [1] будують зважену лінійну комбінацію очікуваних значень, отриманих з кількох неглибших схем Троттера, кожна з яких використовує різну кількість кроків Троттера $k_1, k_2, \ldots, k_r$ (набір з $r$ значень кроків):

$$\langle A \rangle_{\text{MPF}}(t) = \sum_{j=1}^r x_j \, \langle A \rangle_{k_j}(t),$$

де $\langle A \rangle_{k_j}(t)$ — очікуване значення спостережуваної величини $A$ у момент часу $t$, оцінене зі схеми Троттера з $k_j$ кроками, а коефіцієнти $\{x_j\}_{j=1}^r$ вибираються так, щоб провідні члени похибки Троттера в комбінації скасовувалися. Ми повернемося до цього виразу на Кроці 4, де оцінимо його явно, щоб об'єднати наші результати Троттера. Ключовий практичний момент полягає в тому, що найглибша схема в MPF потребує лише $k_{\max}$ кроків, що значно менше за єдине $k$, яке було б потрібне для досягнення тієї ж ефективної похибки Троттера безпосередньо. Неглибші схеми роблять підхід MPF краще пристосованим до зашумленого апаратного забезпечення.

Як визначаються коефіцієнти?

Існує два сімейства коефіцієнтів MPF:

Статичні коефіцієнти не залежать від гамільтоніана, початкового стану та часу еволюції. Вони знаходяться шляхом розв'язання лінійної системи $Ax = b$, яка забезпечує скасування провідних членів похибки Троттера. Для набору кроків Троттера $\{k_j\}_{j=1}^r$, що використовуються з симетричною продуктною формулою порядку $2\chi$, розкладання похибки Троттера в оберненому степені $k_j$ приводить до рівнянь-обмежень вигляду:

$$\sum_{j=1}^r x_j = 1, \quad \sum_{j=1}^r \frac{x_j}{k_j^{\eta_n}} = 0 \quad (n = 0, \ldots, r-2),$$

де цілі показники степеня $\{\eta_n\}$ — це порядки послідовних членів похибки Троттера для обраної продуктної формули. Для симетричної PF порядку $2\chi$ провідна похибка в $\left[S_{2\chi}(t/k)\right]^k$ масштабується як $1/k^{2\chi}$, з подальшими поправками при $1/k^{2\chi+2}, 1/k^{2\chi+4}, \ldots$ — тож показники степеня $\eta_n = 2\chi + 2n$. Для несиметричних PF внесок дають як непарні, так і парні степені, і $\eta_n = 2\chi + n$. Див. посилання [1] для повного виведення. Перше рівняння в системі вище забезпечує незміщеність (MPF відтворює точне очікуване значення в границі $k_j \to \infty$), а решта $r-1$ рівнянь послідовно скасовують перші $r-1$ членів похибки Троттера. Коли отримана $L_1$-норма $\|x\|_1$ занадто велика (що посилює шум вибірки), замість цього можна розв'язати наближену оптимізацію, яка обмежує $\|x\|_1$, мінімізуючи $\|Ax - b\|$.

Динамічні коефіцієнти [2], [3] додатково залежать від гамільтоніана, початкового стану та часу еволюції $t$. Вони мінімізують відстань за нормою Фробеніуса між справжнім еволюціонованим у часі станом і наближенням MPF:

$$\|\rho(t) - \mu^D(t)\|_F^2 = 1 + \sum_{i,j} M_{ij}(t)\, x_i(t)\, x_j(t) - 2\sum_i L_i(t)\, x_i(t),$$

де $M_{ij}(t) = \mathrm{Tr}[\rho_{k_i}(t)\,\rho_{k_j}(t)]$ — матриця Грама перекриттів між еволюціонованими за Троттером станами для різної кількості кроків $k_i, k_j$, а $L_i(t) = \mathrm{Tr}[\rho(t)\,\rho_{k_i}(t)]$ вимірює перекриття з (наближеним) точним станом. У цьому туторіалі ці величини обчислюються ефективно за допомогою методів тензорних мереж, зокрема бекендів на основі TeNPy у qiskit_addon_mpf.

Коли використовувати MPF

MPF найбільш корисні, коли:

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

• Тобі потрібні точні очікувані значення, а не повна підготовка стану. MPF працюють на рівні очікуваних значень — вони комбінують класичні числа, а не квантові стани. Тому вони ідеально підходять для оцінки спостережуваних величин при використанні примітиву Estimator.

• Ти комбінуєш помірну кількість значень кроків Троттера. Зазвичай комбінування $r = 3$–$5$ різних значень кроків $k_j$ достатньо для скасування кількох провідних членів похибки Троттера, зберігаючи при цьому $\|x\|_1$ керованим.

Коли MPF можуть не допомогти

• Дуже короткі часи еволюції. Коли $t$ достатньо мале, щоб одна формула Троттера низького порядку вже була точною, накладні витрати на запуск кількох схем є непотрібними.

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

• Значення кроків Троттера, що порушують режим збіжності. Виведення статичних коефіцієнтів розкладає кожен окремий $\left[S_{2\chi}(t/k_j)\right]^{k_j}$ у ряд за $t/k_j$; це розкладання добре збігається лише коли $t/k_{\min} \lesssim 1$. Якщо $k_{\min}$ обрано занадто малим для заданого $t$, найнеглибша схема далеко виходить за межі пертурбативного режиму, скасовані MPF члени похибки вищого порядку стають великими, і скасування може вимагати великих коефіцієнтів. $L_1$-норма $\|x\|_1$ є практичним діагностичним показником: коли $\|x\|_1 \gg 1$, накладні витрати на вибірку $\propto \|x\|_1^2$ можуть переважити зменшення похибки Троттера. Див. посібник з вибору кроків Троттера для деталей.

Що охоплює цей туторіал

Цей туторіал проводить через повний робочий процес MPF від початку до кінця у два етапи. Спочатку приклад малого симулятора (10-кубітний ланцюг Гейзенберга) демонструє, як налаштувати задачу, обчислити статичні та динамічні коефіцієнти MPF та порівняти отримані очікувані значення з точною діагоналізацією. Потім приклад великого апаратного забезпечення (50-кубітний ланцюг XXZ) показує, як транспілювати, виконувати на апаратному забезпеченні IBM Quantum зі зменшенням похибок та постобробляти результати за допомогою коефіцієнтів MPF. Протягом усього туторіалу ми використовуємо пакет qiskit_addon_mpf разом зі стандартними інструментами Qiskit.

Вимоги

Перед початком цього посібника переконайся, що у тебе встановлено наступне:

• Qiskit SDK v2.0 або пізніша версія, з підтримкою візуалізації

• Qiskit Runtime v0.22 або пізніша версія (pip install qiskit-ibm-runtime)

• Симулятор Qiskit Aer (pip install qiskit-aer)

• Доповнення MPF Qiskit з бекендом TeNPy (pip install "qiskit-addon-mpf[tenpy]")

• Утиліти доповнень Qiskit (pip install qiskit-addon-utils)

• SciPy (pip install scipy)

Налаштування

Нижче ми збираємо всі імпорти пакетів, що використовуються протягом цього туторіалу, в одній комірці. Ми також визначаємо транспіляторний прохід CollectAndCollapse, який об'єднує сусідні обертання rxx і ryy в один XXPlusYYGate. Цей прохід застосовується як під час побудови схеми у Кроці 1 (щоб зберегти низьку кількість вентилів), так і опосередковано, коли ми витягуємо шарувату структуру для динамічного MPF у Кроці 4 (TeNPy очікує двокубітні вентилі, а не пари необ'єднаних обертань).

# Added by doQumentation — installs the packages this notebook needs if they are missing
import importlib.util

_needed = {"matplotlib": "matplotlib", "numpy": "numpy", "qiskit": "qiskit", "qiskit_addon_mpf": "qiskit-addon-mpf", "qiskit_addon_utils": "qiskit-addon-utils", "qiskit_aer": "qiskit-aer", "qiskit_ibm_runtime": "qiskit-ibm-runtime", "scipy": "scipy", "tenpy": "physics-tenpy"}
_missing = [pip for module, pip in _needed.items()
            if importlib.util.find_spec(module) is None]
# One at a time, so a package that fails to install does not block the others
for _pip in _missing:
    %pip install -q {_pip}
if not _missing:
    print("\u2713 All required packages are installed")

import warnings

import numpy as np
import matplotlib.pyplot as plt
from functools import partial
from copy import deepcopy

from qiskit import QuantumCircuit
from qiskit.quantum_info import Pauli, SparsePauliOp, Statevector
from qiskit.synthesis import SuzukiTrotter
from qiskit.transpiler import CouplingMap, PassManager
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
from qiskit.circuit.library import XXPlusYYGate
from qiskit.transpiler.passes.optimization.collect_and_collapse import (
    CollectAndCollapse,
    collect_using_filter_function,
    collapse_to_operation,
)

from qiskit_aer import AerSimulator
from qiskit_ibm_runtime import EstimatorV2 as Estimator, QiskitRuntimeService

from qiskit_addon_utils.problem_generators import (
    generate_xyz_hamiltonian,
    generate_time_evolution_circuit,
)
from qiskit_addon_utils.slicing import slice_by_depth
from qiskit_addon_mpf.static import setup_static_lse
from qiskit_addon_mpf.dynamic import setup_dynamic_lse
from qiskit_addon_mpf.costs import (
    setup_exact_problem,
    setup_sum_of_squares_problem,
    setup_frobenius_problem,
)
from qiskit_addon_mpf.backends.tenpy_layers import (
    LayerModel,
    LayerwiseEvolver,
)
from qiskit_addon_mpf.backends.tenpy_tebd import MPOState, MPS_neel_state

from scipy.linalg import expm

# Suppress TeNPy's `unit_cell_width` future-API warning. The default
# (`unit_cell_width=len(sites)`) is correct for Chain lattices, which is what
# `CouplingMap.from_line(...)` produces here, so the warning is informational.
warnings.filterwarnings(
    "ignore",
    message=r".*unit_cell_width.*",
    category=UserWarning,
)

# --- Helper: collect XX + YY rotations into a single gate ---
def filter_function(node):
    return node.op.name in {"rxx", "ryy"}

collect_function = partial(
    collect_using_filter_function,
    filter_function=filter_function,
    split_blocks=True,
    min_block_size=1,
)

def collapse_to_xx_plus_yy(block):
    param = 0.0
    for node in block.data:
        param += node.operation.params[0]
    return XXPlusYYGate(param)

collapse_function = partial(
    collapse_to_operation,
    collapse_function=collapse_to_xx_plus_yy,
)

pm = PassManager()
pm.append(CollectAndCollapse(collect_function, collapse_function))

Приклад малого симулятора

Крок 1: Відображення класичних вхідних даних на квантову задачу

Ми починаємо з 10-кубітної моделі Гейзенберга на лінії, використовуючи стан Нееля $\vert 0101\ldots01 \rangle$ як початковий стан. Гамільтоніан має вигляд:

$$\hat{\mathcal{H}}_{\text{Heis}} = J \sum_{i=1}^{L-1} \left(X_i X_{i+1} + Y_i Y_{i+1} + Z_i Z_{i+1}\right),$$

де $J$ — сила зв'язку між найближчими сусідами. Ми вимірюємо ZZ-корелятор $Z_{L/2-1} Z_{L/2}$ на парі кубітів у середині ланцюга та використовуємо кроки Троттера $k_j = [1, 2, 4]$ з продуктною формулою другого порядку.

L = 10

# Generate coupling map and Hamiltonian
coupling_map = CouplingMap.from_line(L, bidirectional=False)

hamiltonian = generate_xyz_hamiltonian(
    coupling_map,
    coupling_constants=(1.0, 1.0, 1.0),
    ext_magnetic_field=(0.0, 0.0, 0.0),
)
print(hamiltonian)

SparsePauliOp(['IIIIIIIXXI', 'IIIIIIIYYI', 'IIIIIIIZZI', 'IIIIIXXIII', 'IIIIIYYIII', 'IIIIIZZIII', 'IIIXXIIIII', 'IIIYYIIIII', 'IIIZZIIIII', 'IXXIIIIIII', 'IYYIIIIIII', 'IZZIIIIIII', 'IIIIIIIIXX', 'IIIIIIIIYY', 'IIIIIIIIZZ', 'IIIIIIXXII', 'IIIIIIYYII', 'IIIIIIZZII', 'IIIIXXIIII', 'IIIIYYIIII', 'IIIIZZIIII', 'IIXXIIIIII', 'IIYYIIIIII', 'IIZZIIIIII', 'XXIIIIIIII', 'YYIIIIIIII', 'ZZIIIIIIII'],
              coeffs=[1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j,
 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j,
 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j])

# Observable: ZZ on the middle pair of qubits
observable = SparsePauliOp.from_sparse_list(
    [("ZZ", (L // 2 - 1, L // 2), 1.0)], num_qubits=L
)
print(observable)

SparsePauliOp(['IIIIZZIIII'],
              coeffs=[1.+0.j])

# MPF parameters
mpf_trotter_steps = [1, 2, 4]
order = 2
symmetric = False

trotter_times = np.arange(0.5, 1.55, 0.1)
exact_evolution_times = np.arange(trotter_times[0], 1.55, 0.05)

Побудова схем Троттера

Ми створюємо схеми, що реалізують наближені часові еволюції Троттера для кожної точки часу та кожної кількості кроків Троттера. Прохід CollectAndCollapse, визначений у розділі Налаштування, збирає обертання XX і YY в окремі вентилі XX+YY, щоб підготуватися до ефективнішого моделювання тензорних мереж пізніше.

# Initial Neel state preparation
initial_state_circ = QuantumCircuit(L)
initial_state_circ.x([i for i in range(L) if i % 2 != 0])

all_circs = []
for total_time in trotter_times:
    mpf_trotter_circs = [
        generate_time_evolution_circuit(
            hamiltonian,
            time=total_time,
            synthesis=SuzukiTrotter(reps=num_steps, order=order),
        )
        for num_steps in mpf_trotter_steps
    ]

    mpf_trotter_circs = pm.run(
        mpf_trotter_circs
    )  # Collect XX and YY into XX + YY

    mpf_circuits = [
        initial_state_circ.compose(circuit) for circuit in mpf_trotter_circs
    ]
    all_circs.append(mpf_circuits)

mpf_circuits[-1].draw("mpl", fold=-1)

Крок 2: Оптимізація задачі для виконання на квантовому обладнанні

Для прикладу малого масштабу ми орієнтуємось на симулятор Aer. Перед тим, як схеми будуть готові до виконання, відбуваються два перетворення:

1. Збирання вентилів на рівні моделювання гамільтоніана. У комірці Налаштування ми побудували прохід CollectAndCollapse, який об'єднує сусідні обертання rxx і ryy в один XXPlusYYGate. Ми вже застосували цей прохід, коли будували схеми Троттера у Кроці 1 (виклик pm.run(...)). Це зменшує кількість двокубітних вентилів і створює структуру, яка більш придатна для моделювання тензорних мереж для обчислення динамічних коефіцієнтів пізніше.

2. Зниження до ISA симулятора. Нижче ми запускаємо готовий менеджер проходів Qiskit з optimization_level=3, щоб знизити кожну схему Троттера до архітектури набору інструкцій (ISA) симулятора.

aer_sim = AerSimulator()
pm_sim = generate_preset_pass_manager(backend=aer_sim, optimization_level=3)

isa_circs_all_times = [
    pm_sim.run([deepcopy(c) for c in mpf_circuits])
    for mpf_circuits in all_circs
]

Крок 3: Виконання з використанням примітивів Qiskit

Для прикладу малого масштабу ми запускаємо трансльовані під ISA кола Троттера через примітив EstimatorV2 на основі Aer. Це дає нам безшумове еталонне значення для кожної пари $(k_j, t)$ — це значення $\langle A \rangle_{k_j}(t)$, які MPF об'єднає на Кроці 4. Ми проходимо по діапазону часів еволюції, щоб потім побудувати повну часову криву для кожної окремої формули добутку та для MPF.

estimator = Estimator(mode=aer_sim)

mpf_expvals_all_times, mpf_stds_all_times = [], []
for isa_circuits in isa_circs_all_times:
    result = estimator.run(
        [(circuit, observable) for circuit in isa_circuits], precision=0.005
    ).result()
    mpf_expvals_all_times.append([res.data.evs for res in result])
    mpf_stds_all_times.append([res.data.stds for res in result])

Крок 4: Постобробка та повернення результату у бажаному класичному форматі

Крок 4 — це той крок, де MPF власне будується. Хоча коефіцієнти $x_j$ обчислюються саме тут (і для динамічного варіанту це обчислення може бути ресурсомістким), концептуально вони є класичним рецептом об'єднання квантових вимірювань з Кроку 3 в одне скориговане очікуване значення — тому ми розглядаємо весь процес обчислення коефіцієнтів та об'єднання як постобробку.

Щоб оцінити, наскільки добре MPF відтворює справжню динаміку, ми спочатку обчислюємо точні еволюціонована у часі очікувані значення шляхом безпосереднього експоненціювання гамільтоніана. Це можливо лише тому, що $L = 10$; у прикладі великого масштабу на апаратному забезпеченні нижче нам доведеться натомість покладатися на оцінки за допомогою тензорних мереж.

exact_expvals = []
for t in exact_evolution_times:
    exp_H = expm(-1j * t * hamiltonian.to_matrix())
    initial_state = Statevector(initial_state_circ).data
    time_evolved_state = exp_H @ initial_state

    exact_obs = (
        time_evolved_state.conj()
        @ observable.to_matrix()
        @ time_evolved_state
    ).real
    exact_expvals.append(exact_obs)

Статичні коефіцієнти MPF

Статичні MPF використовують коефіцієнти $x_j$, які не залежать від часу еволюції, гамільтоніана та початкового стану. Ми складаємо лінійну систему $Ax = b$, описану в розділі «Довідкова інформація», і розв'язуємо її для коефіцієнтів. Матриця $A$ визначається кількістю кроків Троттера $k_j$, порядком $\chi$ формули добутку та тим, чи є формула симетричною (що визначає показники степенів $\eta_n$).

Для нашого прикладу малого масштабу ми використовуємо $k_j = [1, 2, 4]$ з несиметричною формулою Судзукі-Троттера порядку $2\chi=2$ (тобто $\chi=1$ та $\eta_n = 2 + n$, що дає $\eta_0 = 2,\, \eta_1 = 3$). Система набуває вигляду:

$$A = \begin{bmatrix} 1 & 1 & 1\\ 1 & \frac{1}{2^2} & \frac{1}{4^2} \\ 1 & \frac{1}{2^3} & \frac{1}{4^3} \\ \end{bmatrix}, \quad b = \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix}.$$

Перший рядок забезпечує незміщеність ($\sum_j x_j = 1$); другий та третій рядки скасовують провідні члени похибки Троттера $1/k^2$ та наступного порядку $1/k^3$ відповідно.

Налаштування LSE

Ми використовуємо setup_static_lse з qiskit_addon_mpf.static, щоб скласти матрицю $A$ та вектор правої частини $b$, описані вище. Матриця $A$ залежить не лише від $k_j$, а й від нашого вибору формули добутку — зокрема, від її порядку $\chi$ та того, чи є вона симетричною. Прапорець symmetric визначає схему показників степенів $\eta_n$ (симетричні формули дають лише члени похибки Троттера парного порядку; див. посилання [1]). Зауваж, що, як показано в посиланні [2], встановлення symmetric=True не є строго необхідним, навіть коли базова формула добутку є симетричною — несиметрична LSE залишається дійсною (вона накладає додаткові непотрібні обмеження).

Для нашого прикладу ми вже встановили order = 2 та symmetric = False на Кроці 1.

lse = setup_static_lse(mpf_trotter_steps, order=order, symmetric=symmetric)

Перевір побудовану матрицю $A$ та вектор $b$, щоб переконатися, що вони відповідають системі, записаній вище.

lse.A

array([[1.      , 1.      , 1.      ],
       [1.      , 0.25    , 0.0625  ],
       [1.      , 0.125   , 0.015625]])

lse.b

array([1., 0., 0.])

Маючи LSE, ми розв'язуємо для статичних коефіцієнтів $x_j$ за допомогою lse.solve() (це прямий розв'язок $x = A^{-1}b$).

mpf_coeffs = lse.solve()
print(
    f"The static coefficients associated with the ansatze are: {mpf_coeffs}"
)

The static coefficients associated with the ansatze are: [ 0.04761905 -0.57142857  1.52380952]

Оптимізація для $x$ з використанням точної моделі

Як альтернативу обчисленню $x = A^{-1}b$, ти можеш використовувати setup_exact_model для побудови екземпляра cvxpy.Problem, який використовує LSE як обмеження і чиє оптимальне рішення дасть $x$.

model_exact, coeffs_exact = setup_exact_problem(lse)
model_exact.solve()
print(coeffs_exact.value)

[ 0.04761905 -0.57142857  1.52380952]

print(
    "L1 norm of the exact coefficients:",
    np.linalg.norm(coeffs_exact.value, ord=1),
)

L1 norm of the exact coefficients: 2.1428571428556378

Оптимізація для $x$ з використанням наближеної моделі

Може статися, що $L_1$-норма для обраного набору значень $k_j$ вважається занадто високою. Якщо це так, і ти не можеш вибрати інший набір значень $k_j$, ти можеш використовувати наближений розв'язок, що обмежує $L_1$-норму до обраного порогу, мінімізуючи $\|Ax - b\|$. Ознайомся з посібником Як використовувати наближену модель.

model_approx, coeffs_approx = setup_sum_of_squares_problem(
    lse, max_l1_norm=1.5
)
model_approx.solve()
print(coeffs_approx.value)
print(
    "L1 norm of the approximate coefficients:",
    np.linalg.norm(coeffs_approx.value, ord=1),
)

[-1.10294118e-03 -2.48897059e-01  1.25000000e+00]
L1 norm of the approximate coefficients: 1.5

Динамічні коефіцієнти MPF

Статичний MPF скасовує члени похибки Троттера у спосіб, незалежний від гамільтоніана та стану, тому він не обов'язково дає найменшу можливу похибку апроксимації для заданого гамільтоніана та початкового стану. Динамічний MPF (посилання [2], [3]) натомість знаходить залежні від часу коефіцієнти $x_i(t)$, які мінімізують відстань за нормою Фробеніуса $\|\rho(t) - \mu^D(t)\|_F^2$ у кожен момент часу $t$. Як показано в розділі «Довідкова інформація», це потребує матриці перекриття $M_{ij}(t)$ між еволюціонованими Троттером станами та перекриття $L_i(t)$ з точним станом — обидва з яких ми оцінюємо за допомогою бекендів тензорних мереж (TeNPy) у qiskit_addon_mpf.

Щоб скласти динамічну LSE, нам потрібні три складові:

1. Фабрика наближеного еволюціонера, яку додон запускатиме для кожного $k_j$, щоб отримати $\rho_{k_j}(t)$ у вигляді MPS/MPO. Ми будуємо її зі шаруватої структури кола Троттера порядку $2$ (один шар на slice_by_depth), обгорнутого як LayerwiseEvolver з параметрами усічення TeNPy.

2. Фабрика точного еволюціонера, яка дає високоточне еталонне значення $\rho(t)$. Ми використовуємо коло Судзукі-Троттера четвертого порядку з малим кроком за часом (dt=0.1, order=4) як замінник точної еволюції.

3. Фабрика тотожності та початковий стан MPS, які задають початкові умови для симуляції TeNPy.

Комірка нижче будує фабрику наближеного еволюціонера.

# Create approximate time-evolution circuits
single_2nd_order_circ = generate_time_evolution_circuit(
    hamiltonian, time=1.0, synthesis=SuzukiTrotter(reps=1, order=order)
)
single_2nd_order_circ = pm.run(single_2nd_order_circ)  # collect XX and YY

# Find layers in the circuit
layers = slice_by_depth(single_2nd_order_circ, max_slice_depth=1)

# Create tensor network models
models = [
    LayerModel.from_quantum_circuit(layer, conserve="Sz") for layer in layers
]

# Create the time-evolution object
approx_factory = partial(
    LayerwiseEvolver,
    layers=models,
    options={
        "preserve_norm": False,
        "trunc_params": {
            "chi_max": 64,
            "svd_min": 1e-8,
            "trunc_cut": None,
        },
        "max_delta_t": 2,
    },
)

попередження

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

Ми апроксимуємо точний стан, еволюціонований у часі, за допомогою формули Судзукі-Троттера четвертого порядку, використовуючи малий часовий крок dt=0.1. Параметри усічення TeNPy можуть впливати на точність, тому важливо дослідити діапазон значень.

single_4th_order_circ = generate_time_evolution_circuit(
    hamiltonian, time=1.0, synthesis=SuzukiTrotter(reps=1, order=4)
)
single_4th_order_circ = pm.run(single_4th_order_circ)
exact_model_layers = [
    LayerModel.from_quantum_circuit(layer, conserve="Sz")
    for layer in slice_by_depth(single_4th_order_circ, max_slice_depth=1)
]

exact_factory = partial(
    LayerwiseEvolver,
    layers=exact_model_layers,
    dt=0.1,
    options={
        "preserve_norm": False,
        "trunc_params": {
            "chi_max": 64,
            "svd_min": 1e-8,
            "trunc_cut": None,
        },
        "max_delta_t": 2,
    },
)

Нарешті, ми визначаємо identity_factory, яка дає початковий стан MPO, і готуємо початковий стан Нееля у вигляді MPS, що відповідає ґратці, використовуваній у шаруватій моделі Троттера.

def identity_factory():
    return MPOState.initialize_from_lattice(models[0].lat, conserve=True)

mps_initial_state = MPS_neel_state(models[0].lat)

Маючи фабрики на місці, ми тепер обчислюємо динамічні коефіцієнти в кожен момент часу еволюції. Для кожного $t$ setup_dynamic_lse будує відповідні матриці перекриття за допомогою TeNPy, а setup_frobenius_problem повертає cvxpy.Problem, що мінімізує вартість за нормою Фробеніуса. Розв'язувач повертає коефіцієнти $x_j(t)$, підібрані для цього часу; ми збираємо їх у mpf_dynamic_coeffs_list. Якщо розв'язувач не спрацьовує для певного $t$, ми повертаємося до нульових коефіцієнтів, щоб цикл продовжувався.

mpf_dynamic_coeffs_list = []
for t in trotter_times:
    print(f"Computing dynamic coefficients for time={t}")
    lse = setup_dynamic_lse(
        mpf_trotter_steps,
        t,
        identity_factory,
        exact_factory,
        approx_factory,
        mps_initial_state,
    )
    problem, coeffs = setup_frobenius_problem(lse)
    try:
        problem.solve()
        mpf_dynamic_coeffs_list.append(coeffs.value)
    except Exception as error:
        mpf_dynamic_coeffs_list.append(np.zeros(len(mpf_trotter_steps)))
        print(error, "Calculation Failed for time", t)
    print("")

Computing dynamic coefficients for time=0.5

Computing dynamic coefficients for time=0.6

Computing dynamic coefficients for time=0.7

Computing dynamic coefficients for time=0.7999999999999999

Computing dynamic coefficients for time=0.8999999999999999

Computing dynamic coefficients for time=0.9999999999999999

Computing dynamic coefficients for time=1.0999999999999999

Computing dynamic coefficients for time=1.1999999999999997

Computing dynamic coefficients for time=1.2999999999999998

Computing dynamic coefficients for time=1.4

Computing dynamic coefficients for time=1.4999999999999998

Поєднання очікуваних значень Троттера з коефіцієнтами MPF

Тепер ми обчислюємо $\langle A \rangle_{\text{MPF}}(t) = \sum_j x_j \, \langle A \rangle_{k_j}(t)$ для кожного набору коефіцієнтів (статичний точний, статичний наближений та динамічний), поширюємо стандартні похибки на кожне коло та наносимо отриману часову серію проти кривої точної діагоналізації.

sym = {1: "^", 2: "s", 4: "p"}
# Get expectation values at all times for each Trotter step
for k, step in enumerate(mpf_trotter_steps):
    trotter_curve, trotter_curve_error = [], []
    for trotter_expvals, trotter_stds in zip(
        mpf_expvals_all_times, mpf_stds_all_times
    ):
        trotter_curve.append(trotter_expvals[k])
        trotter_curve_error.append(trotter_stds[k])

    plt.errorbar(
        trotter_times,
        trotter_curve,
        yerr=trotter_curve_error,
        alpha=0.5,
        markersize=4,
        marker=sym[step],
        color="grey",
        label=f"{mpf_trotter_steps[k]} Trotter steps",
    )

# Get expectation values at all times for the static MPF with exact coeffs
exact_mpf_curve, exact_mpf_curve_error = [], []
for trotter_expvals, trotter_stds in zip(
    mpf_expvals_all_times, mpf_stds_all_times
):
    mpf_std = np.sqrt(
        sum(
            [
                (coeff**2) * (std**2)
                for coeff, std in zip(coeffs_exact.value, trotter_stds)
            ]
        )
    )
    exact_mpf_curve_error.append(mpf_std)
    exact_mpf_curve.append(trotter_expvals @ coeffs_exact.value)

plt.errorbar(
    trotter_times,
    exact_mpf_curve,
    yerr=exact_mpf_curve_error,
    markersize=4,
    marker="o",
    label="Static MPF - Exact",
    color="purple",
)

# Get expectation values at all times for the static MPF with approximate coeffs
approx_mpf_curve, approx_mpf_curve_error = [], []
for trotter_expvals, trotter_stds in zip(
    mpf_expvals_all_times, mpf_stds_all_times
):
    mpf_std = np.sqrt(
        sum(
            [
                (coeff**2) * (std**2)
                for coeff, std in zip(coeffs_approx.value, trotter_stds)
            ]
        )
    )
    approx_mpf_curve_error.append(mpf_std)
    approx_mpf_curve.append(trotter_expvals @ coeffs_approx.value)

plt.errorbar(
    trotter_times,
    approx_mpf_curve,
    yerr=approx_mpf_curve_error,
    markersize=4,
    marker="o",
    label="Static MPF - Approx",
    color="orange",
)

# Get expectation values at all times for the dynamic MPF
dynamic_mpf_curve, dynamic_mpf_curve_error = [], []
for trotter_expvals, trotter_stds, dynamic_coeffs in zip(
    mpf_expvals_all_times, mpf_stds_all_times, mpf_dynamic_coeffs_list
):
    mpf_std = np.sqrt(
        sum(
            [
                (coeff**2) * (std**2)
                for coeff, std in zip(dynamic_coeffs, trotter_stds)
            ]
        )
    )
    dynamic_mpf_curve_error.append(mpf_std)
    dynamic_mpf_curve.append(trotter_expvals @ dynamic_coeffs)

plt.errorbar(
    trotter_times,
    dynamic_mpf_curve,
    yerr=dynamic_mpf_curve_error,
    markersize=4,
    marker="o",
    label="Dynamic MPF",
    color="pink",
)

# Exact expectation values
plt.plot(
    exact_evolution_times,
    exact_expvals,
    color="red",
    linestyle="--",
    label="Exact time-evolution",
)

plt.title(f"$\\langle Z_{{{L//2-1}}} Z_{{{L//2}}} \\rangle$ vs time")
plt.xlabel("Time")
plt.ylabel("Expectation Value")
plt.legend(loc="upper center", bbox_to_anchor=(0.5, -0.2), ncol=2)
plt.grid(alpha=0.1)
plt.tight_layout()
plt.show()

Графік вище ілюструє взаємодію між похибкою Троттера та похибкою вибірки.

• Похибка Троттера. Окремі формули добутку (сірі маркери) все більше відхиляються від точної кривої зі зростанням часу. Коло $k=1$ має найбільше відхилення і є найдрібнішим, але воно вже перебуває в режимі, де $t/k \gtrsim 1$, тому провідний член похибки $1/k^{2}$ великий. Комбінації MPF (кольорові маркери) скасовують кілька з цих провідних членів похибки Троттера, тому вони відстежують точну криву значно точніше, ніж будь-яке окреме коло $k_j$. Залишковий розрив відображає члени похибки Троттера вищого порядку, які MPF не скасовує: статичний MPF порядку $2$, $r=3$ усуває лише перші два порядки похибки, і при великих $t/k_{\min}$ незвільнений хвіст врешті-решт домінує — тому MPF не гарантує, що дуже дрібні кола залишаються точними в будь-який момент часу.

• Похибка вибірки. Ширші межі похибки на кривих MPF є прямим наслідком лінійної комбінації: поширення незалежних стандартних похибок $\sigma_{k_j}$ для кожного кола дає загальну дисперсію $\sigma_{\text{MPF}}^2 = \sum_j x_j^2 \, \sigma_{k_j}^2$. Тому чим більша $\|x\|_2$ (а на практиці $\|x\|_1$, яку ми контролюємо), тим більше запусків (shots) потрібно для досягнення заданої цільової невизначеності. Це компроміс, що лежить в основі опції наближеного розв'язувача в розділі «Довідкова інформація»: ми обмежуємо $\|x\|_1$, щоб утримати ці накладні витрати у прийнятних межах. Важливо, що, на відміну від похибки Троттера, похибка вибірки зменшується як $1/\sqrt{N_{\text{shots}}}$, тому її завжди можна зменшити, витративши більше запусків.

У прикладі великого масштабу на апаратному забезпеченні нижче апаратний шум входить як додаткове джерело похибки для кожного $\langle A \rangle_{k_j}$, яке подібним чином посилюється коефіцієнтами MPF. Ми побачимо, як пом'якшення похибок взаємодіє з MPF у цьому розділі.

Великомасштабний приклад на апаратному забезпеченні

У цьому розділі ми масштабуємо задачу за межі того, що можна точно симулювати. Ми відтворюємо деякі результати, показані в посиланні [3], використовуючи ланцюг XXZ на 50 кубітів у момент часу $t = 3$. Ми дотримуємося того ж чотирикрокового робочого процесу, що й у прикладі малого масштабу, тепер орієнтуючись на реальне квантове апаратне забезпечення з пом'якшенням похибок. Як і в шаблоні, кожен крок позначено безпосередньо в коді, і один крок може охоплювати кілька комірок, коли варто розглянути його проміжні результати. Відповідність дзеркалить приклад малого масштабу: визначення гамільтоніана, вибір параметрів Троттера, обчислення коефіцієнтів MPF (статичних і динамічних) та побудова кіл. Ключові відмінності такі:

• Гамільтоніан XXZ на 50 вузлах з випадковими зв'язками, взятими з $\mathcal{U}(0.5, 1.5)$ (посилання [3]).

• Симетрична формула Троттера другого порядку з $k_j = [3, 4, 6]$ (тобто $\chi=1$, symmetric=True).

• Один фіксований час еволюції $t = 3$. При $k_{\min}=3$ це дає $t/k_{\min}=1$, зберігаючи дрібні складові в режимі збіжності Троттера, де модель провідної похибки, на яку покладається MPF, є дійсною.

• Додатковий запуск для порівняння з одним колом із $k = 10$ кроками Троттера, використаний як базовий рівень. Ми обрали $k = 10$, оскільки його дволкубітна глибина на апаратному забезпеченні є більшою за найглибшу складову MPF ($k_{\max}=6$) плюс накладні витрати на запуск кількох кіл MPF — достатньо глибокою, щоб бути обмеженою шумом, що є режимом, у якому очікується, що комбінація MPF перевершить базовий рівень з одним колом. Це порівняння «одне глибоке коло» з комбінацією MPF, а не коло, що цілиться в ефективну похибку Троттера MPF (що вимагало б значно більше кроків).

Зауваж, що, хоча ми все ще на Кроці 1 тут (відображення та побудова кола), ми також попередньо обчислюємо динамічні коефіцієнти поряд зі статичними в цій комірці. Динамічні коефіцієнти залежать від $H$ та $t$, але не від квантових вимірювань, тому їх можна обчислити будь-коли до Кроку 4. Ми робимо це зараз, щоб зберегти все налаштування, специфічне для MPF, в одному місці.

# -------------------------Step 1-------------------------
L = 50
coupling_map = CouplingMap.from_line(L, bidirectional=False)

# XXZ Hamiltonian with random couplings (Ref. [3])
np.random.seed(0)
even_edges = list(coupling_map.get_edges())[::2]
odd_edges = list(coupling_map.get_edges())[1::2]

Js = np.random.uniform(0.5, 1.5, size=L)
hamiltonian = SparsePauliOp(Pauli("I" * L))
for i, edge in enumerate(even_edges + odd_edges):
    hamiltonian += SparsePauliOp.from_sparse_list(
        [
            ("XX", (edge), 2 * Js[i]),
            ("YY", (edge), 2 * Js[i]),
            ("ZZ", (edge), 4 * Js[i]),
        ],
        num_qubits=L,
    )

observable = SparsePauliOp.from_sparse_list(
    [("ZZ", (L // 2 - 1, L // 2), 1.0)], num_qubits=L
)

total_time = 3
mpf_trotter_steps = [3, 4, 6]
order = 2
symmetric = True

# Static coefficients
lse = setup_static_lse(mpf_trotter_steps, order=order, symmetric=symmetric)
mpf_coeffs = lse.solve()
print(f"Static coefficients: {mpf_coeffs}")
print(f"L1 norm: {np.linalg.norm(mpf_coeffs, ord=1)}")

model_approx, coeffs_approx = setup_sum_of_squares_problem(
    lse, max_l1_norm=2.0
)
model_approx.solve()
print(f"Approximate coefficients: {coeffs_approx.value}")
print(f"L1 norm (approx): {np.linalg.norm(coeffs_approx.value, ord=1)}")

# -------------------------Dynamic coefficients-------------------------
single_2nd_order_circ = generate_time_evolution_circuit(
    hamiltonian, time=1.0, synthesis=SuzukiTrotter(reps=1, order=order)
)
single_2nd_order_circ = pm.run(single_2nd_order_circ)

layers = slice_by_depth(single_2nd_order_circ, max_slice_depth=1)
models = [
    LayerModel.from_quantum_circuit(layer, conserve="Sz") for layer in layers
]

approx_factory = partial(
    LayerwiseEvolver,
    layers=models,
    options={
        "preserve_norm": False,
        "trunc_params": {"chi_max": 64, "svd_min": 1e-8, "trunc_cut": None},
        "max_delta_t": 4,
    },
)

single_4th_order_circ = generate_time_evolution_circuit(
    hamiltonian, time=1.0, synthesis=SuzukiTrotter(reps=1, order=4)
)
single_4th_order_circ = pm.run(single_4th_order_circ)
exact_model_layers = [
    LayerModel.from_quantum_circuit(layer, conserve="Sz")
    for layer in slice_by_depth(single_4th_order_circ, max_slice_depth=1)
]

exact_factory = partial(
    LayerwiseEvolver,
    layers=exact_model_layers,
    dt=0.1,
    options={
        "preserve_norm": False,
        "trunc_params": {"chi_max": 64, "svd_min": 1e-8, "trunc_cut": None},
        "max_delta_t": 3,
    },
)

def identity_factory():
    return MPOState.initialize_from_lattice(models[0].lat, conserve=True)

mps_initial_state = MPS_neel_state(models[0].lat)

print(f"Computing dynamic coefficients for time={total_time}")
lse_dyn = setup_dynamic_lse(
    mpf_trotter_steps,
    total_time,
    identity_factory,
    exact_factory,
    approx_factory,
    mps_initial_state,
)
problem, coeffs_dyn = setup_frobenius_problem(lse_dyn)
try:
    problem.solve()
    mpf_dynamic_coeffs = coeffs_dyn.value
except Exception as error:
    mpf_dynamic_coeffs = np.zeros(len(mpf_trotter_steps))
    print(error, "Calculation Failed")

# -------------------------Step 1 (cont): Build circuits-------------------------
mpf_circuits = []
for k in mpf_trotter_steps:
    circuit = QuantumCircuit(L)
    circuit.x([i for i in range(L) if i % 2])
    trotter_circ = generate_time_evolution_circuit(
        hamiltonian,
        synthesis=SuzukiTrotter(reps=k, order=order),
        time=total_time,
    )
    circuit.compose(trotter_circ, qubits=range(L), inplace=True)
    mpf_circuits.append(circuit)

# Baseline "single deep circuit" comparison run with k=10 Trotter steps.
# Its two-qubit depth is deeper than the deepest MPF constituent (k_max=6) plus
# the overhead of running multiple circuits, pushing it into the noise-limited
# regime where MPF is expected to outperform. It does NOT target the MPF's effective
# Trotter error (which would require many more steps).
comp_circuit = QuantumCircuit(L)
comp_circuit.x([i for i in range(L) if i % 2])
trotter_circ = generate_time_evolution_circuit(
    hamiltonian,
    synthesis=SuzukiTrotter(reps=10, order=order),
    time=total_time,
)
comp_circuit.compose(trotter_circ, qubits=range(L), inplace=True)
mpf_circuits.append(comp_circuit)

Static coefficients: [ 0.42857143 -1.82857143  2.4       ]
L1 norm: 4.65714285714286
Approximate coefficients: [-0.4942491   0.40206845  1.09218065]
L1 norm (approx): 1.9884981979026675
Computing dynamic coefficients for time=3

Тепер ми оптимізуємо кола для обраного бекенду. Ми використовуємо готовий менеджер проходів Qiskit з optimization_level=3, який автоматично обирає добрий набір фізичних кубітів і маршрутизує кожне коло на топологію пристрою.

# -------------------------Step 2-------------------------
service = QiskitRuntimeService()
# backend = service.least_busy(operational=True, simulator=False, min_num_qubits=L)
backend = service.backend("ibm_fez")
print(backend)

transpiler = generate_preset_pass_manager(
    optimization_level=3, backend=backend
)
transpiled_circuits = [transpiler.run(circ) for circ in mpf_circuits]

isa_observables = [
    observable.apply_layout(circ.layout) for circ in transpiled_circuits
]

<IBMBackend('ibm_fez')>

Запуск глибших кіл на реальному апаратному забезпеченні вимагає агресивного пом'якшення похибок. Ми вмикаємо динамічне розчеплення, твірлінг вентилів і вимірювань, пом'якшення похибок вимірювання та екстраполяцію до нульового шуму (ZNE). Зауваж, що фактори шуму ZNE, які ми використовуємо тут (1, 1.2, 1.4), менші, ніж у сценарії з дрібним колом, оскільки глибші складові MPF вже близькі до порогу шуму, і великі посилення шуму виштовхнули б їх за межу, де екстраполяція ZNE є надійною.

Ми надсилаємо всі чотири кола (три складові MPF при $k_j = [3, 4, 6]$ плюс базовий рівень $k = 10$) в одному завданні Estimator.

# -------------------------Step 3-------------------------
estimator = Estimator(mode=backend)
estimator.options.default_shots = 30000

# Error suppression/mitigation
estimator.options.dynamical_decoupling.enable = True
estimator.options.twirling.enable_gates = True
estimator.options.twirling.enable_measure = True
estimator.options.twirling.num_randomizations = "auto"
estimator.options.twirling.strategy = "active-accum"
estimator.options.resilience.measure_mitigation = True
estimator.options.experimental.execution_path = "gen3-turbo"

estimator.options.resilience.zne_mitigation = True
estimator.options.resilience.zne.noise_factors = (1, 1.2, 1.4)
estimator.options.resilience.zne.extrapolator = "linear"

estimator.options.environment.job_tags = ["TUT_MPF"]

job_50 = estimator.run(
    [
        (circ, observable)
        for circ, observable in zip(transpiled_circuits, isa_observables)
    ]
)

Ми отримуємо очікувані значення та стандартні відхилення для кожного кола з результату завдання, а потім об'єднуємо їх з кожним набором коефіцієнтів MPF точно так само, як у прикладі малого масштабу: $\langle A \rangle_{\text{MPF}} = \sum_j x_j \, \langle A \rangle_{k_j}$, з поширеною дисперсією $\sigma^2 = \sum_j x_j^2 \sigma_{k_j}^2$.

# -------------------------Step 4-------------------------
result = job_50.result()
evs = [res.data.evs for res in result]
std = [res.data.stds for res in result]

print(evs)
print(std)

[array(-0.07916195), array(-0.04479681), array(-0.2560756), array(-0.06045848)]
[array(0.04605538), array(0.10056336), array(0.14426151), array(0.04059092)]

exact_mpf_std = np.sqrt(
    sum([(coeff**2) * (std**2) for coeff, std in zip(mpf_coeffs, std[:3])])
)
print(
    "Exact static MPF expectation value: ",
    evs[:3] @ mpf_coeffs,
    "+-",
    exact_mpf_std,
)
approx_mpf_std = np.sqrt(
    sum(
        [
            (coeff**2) * (std**2)
            for coeff, std in zip(coeffs_approx.value, std[:3])
        ]
    )
)
print(
    "Approximate static MPF expectation value: ",
    evs[:3] @ coeffs_approx.value,
    "+-",
    approx_mpf_std,
)
dynamic_mpf_std = np.sqrt(
    sum(
        [
            (coeff**2) * (std**2)
            for coeff, std in zip(mpf_dynamic_coeffs, std[:3])
        ]
    )
)
print(
    "Dynamic MPF expectation value: ",
    evs[:3] @ mpf_dynamic_coeffs,
    "+-",
    dynamic_mpf_std,
)

Exact static MPF expectation value:  -0.5665938395816946 +- 0.3925273058119915
Approximate static MPF expectation value:  -0.25856647611537903 +- 0.164249927266166
Dynamic MPF expectation value:  -0.12667812062949296 +- 0.06059471006973169

sym = {3: "^", 4: "s", 6: "p"}
for k, step in enumerate(mpf_trotter_steps):
    plt.errorbar(
        k,
        evs[k],
        yerr=std[k],
        alpha=0.5,
        markersize=4,
        marker=sym[step],
        color="grey",
        label=f"{mpf_trotter_steps[k]} Trotter steps",
    )

plt.errorbar(
    3,
    evs[-1],
    yerr=std[-1],
    alpha=0.5,
    markersize=8,
    marker="x",
    color="blue",
    label="10 Trotter steps",
)

plt.errorbar(
    4,
    evs[:3] @ mpf_coeffs,
    yerr=exact_mpf_std,
    markersize=4,
    marker="o",
    color="purple",
    label="Static MPF",
)

plt.errorbar(
    5,
    evs[:3] @ coeffs_approx.value,
    yerr=approx_mpf_std,
    markersize=4,
    marker="o",
    color="orange",
    label="Approximate static MPF",
)

plt.errorbar(
    6,
    evs[:3] @ mpf_dynamic_coeffs,
    yerr=dynamic_mpf_std,
    markersize=4,
    marker="o",
    color="pink",
    label="Dynamic MPF",
)

exact_obs = -0.24384471447172074  # Calculated via Tensor Network calculation
plt.axhline(
    y=exact_obs, linestyle="--", color="red", label="Exact time-evolution"
)

plt.title(
    f"$\\langle Z_{{{L//2-1}}} Z_{{{L//2}}} \\rangle$ at time {total_time} for the different methods"
)
plt.xlabel("Method")
plt.ylabel("Expectation Value")
plt.legend(loc="upper center", bbox_to_anchor=(0.5, -0.2), ncol=2)
plt.grid(alpha=0.1)
plt.tight_layout()
plt.show()

Кілька спостережень про наведені вище результати на апаратному забезпеченні:

• Заглиблення не безкоштовне на апаратному забезпеченні. Базові рівні з одним колом розповідають про це прямо: коло $k = 6$ по суті точне ($-0.256$ проти еталонного значення $-0.244$), тоді як глибший базовий рівень $k = 10$ гірший ($-0.061$, відхилення на $\sim 0.18$), а не кращий. Коли похибка Троттера вже мала, додавання кроків здебільшого поглиблює коло і накопичує більше шуму вентилів та декогеренції. Це саме той режим, для якого побудовано MPF: досягти точності глибокого кола, використовуючи лише дрібні складові.

• MPF з малою нормою перевершує глибоке окреме коло. Наближений статичний MPF (обмежений на $\|x\|_1 \approx 2$) дає $-0.259$, у межах $\sim 0.015$ від еталонного значення і значно ближче, ніж базовий рівень $k = 10$. Динамічний MPF ($-0.127$) також помітно перевершує цей базовий рівень. Обидва об'єднують лише дрібні кола $k_j = [3, 4, 6]$, проте відновлюють відповідь, яку глибоке окреме коло не змогло отримати.

• Норма коефіцієнтів важливіша за математичну оптимальність. Точний статичний MPF має $\|x\|_1 = 4.66$ і є найгіршим оцінювачем з усіх ($-0.567$, відхилення понад $0.3$): велика норма коефіцієнтів посилює залишковий шум вентилів, декогеренцію та похибку ZNE для кожного $\langle A \rangle_{k_j}$ приблизно на той самий фактор, переважаючи скасування похибки Троттера, яке вона забезпечує. Обмеження норми (наближений статичний розв'язувач, $\|x\|_1 \approx 2$) усуває це переважання і дає найкращу оцінку — навіть попри те, що його коефіцієнти вже не скасовують провідну похибку Троттера точно.

• Окремі дрібні кола все ще можуть бути конкурентоспроможними. Одинока складова $k = 6$ ($-0.256$) сама по собі тут по суті точна — у цьому запуску вона навіть трохи ближча, ніж наближений статичний MPF. Проблема в тому, що заздалегідь невідомо, яке саме $k$ потрапляє у «солодку точку» «збіжне, але ще не обмежене шумом», а безпечний на вигляд вибір просто поглибити ($k = 10$), щоб гарантувати збіжність Троттера, є саме тим, що зазнає невдачі. MPF дає принципову комбінацію дрібних кіл, яка не вимагає вгадування правильної глибини.

Практичний висновок полягає в тому, що на апаратному забезпеченні MPF слід поєднувати із сильним пом'якшенням похибок для кожного окремого $\langle A \rangle_{k_j}$, $L_1$-норму коефіцієнтів слід тримати помірною (використовуй наближений розв'язувач або динамічний MPF), а кроки Троттера $k_j$ слід обирати так, щоб $t/k_{\min} \lesssim 1$ — тут $k_{\min} = 3$ при $t = 3$ дає $t/k_{\min} = 1$, зберігаючи складові в межах збіжного режиму, де модель провідної похибки, на яку покладається статичний MPF, є дійсною. З цими виборами MPF з малою нормою тут відповідають збіжному окремому колу, тоді як наївний базовий рівень «просто поглиблюй» — ні, відновлюючи перевагу глибини над точністю, показану в посиланні [3]. Зауваж також, що окремі запуски шумні — при іншому надсиланні того самого завдання (або на іншому бекенді) точний порядок може змінитися; надійні тенденції такі: MPF з малим $\|x\|_1$ показують добрі результати, точний статичний MPF з великим $\|x\|_1$ посилюється апаратним шумом, а надто глибоке окреме коло обмежене шумом.

Наступні кроки

Рекомендації

Якщо ця робота видалася тобі цікавою, тебе можуть зацікавити такі матеріали:

• Як обрати кроки Троттера для MPF — практичні поради щодо вибору значень $k_j$, щоб уникнути нестабільностей
• Як використовувати наближену модель — налаштування обмеження $L_1$-норми та опцій розв'язувача для наближеного статичного MPF
• Довідник API qiskit-addon-mpf — повна документація для статичних, динамічних та модулів бекендів

Посилання

[1] Vazquez, A. C., Egger, D. J., Ochsner, D., & Woerner, S. Well-conditioned multi-product formulas for hardware-friendly Hamiltonian simulation. Quantum, 7, 1067 (2023)

[2] Zhuk, S., Robertson, N. F., & Bravyi, S. Trotter error bounds and dynamic multi-product formulas for Hamiltonian simulation. Physical Review Research, 6(3), 033309 (2024)

[3] Robertson, N. F., et al. Tensor network enhanced dynamic multiproduct formulas. arXiv:2407.17405 (2024)