Перейти до основного вмісту

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

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

Передумови

Цей посібник демонструє, як використовувати багатопродуктну формулу (MPF) для досягнення меншої похибки Троттера для нашого спостережуваного порівняно з тією, яка виникає при найглибшому колі Троттера, яке ми фактично виконаємо. MPF зменшують похибку Троттера гамільтонової динаміки через зважену комбінацію кількох виконань кола. Розглянемо задачу знаходження очікувань спостережуваних для квантового стану ρ(t)=eiHtρ(0)eiHt\rho(t)=e^{-i H t} \rho(0) e^{i H t} з гамільтоніаном HH. Можна використовувати продуктні формули (PF) для апроксимації еволюції в часі eiHte^{-i H t}, виконавши наступне:

  • Записати гамільтоніан HH як H=a=1dFa,H=\sum_{a=1}^d F_a, де FaF_a є ермітовими операторами такими, що кожен відповідний унітарний оператор може бути ефективно реалізований на квантовому пристрої.
  • Апроксимувати члени FaF_a, які не комутують один з одним.

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

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

яка має квадратичний член похибки S1(t)=eiHt+O(t2)S_1(t)=e^{-i H t}+\mathcal{O}\left(t^{2}\right). Можна також використовувати продуктні формули вищих порядків (формули Лі-Троттера-Сузукі), які збігаються швидше і визначаються рекурсивно як:

S2(t):=a=1deiFat/2a=1deiFat/2S_2(t):=\prod_{a=1}^d e^{-i F_a t/2}\prod_{a=1}^d e^{-i F_a t/2}

S2χ(t):=S2χ2(sχt)2S2χ2((14sχ)t)S2χ2(sχt)2,S_{2 \chi}(t):= S_{2 \chi -2}(s_{\chi}t)^2 S_{2 \chi -2}((1-4s_{\chi})t)S_{2 \chi -2}(s_{\chi}t)^2,

де χ\chi є порядком симетричної PF і sp=(441/(2p1))1s_p = \left( 4 - 4^{1/(2p-1)} \right)^{-1}. Для довгих еволюцій у часі можна розділити часовий інтервал tt на kk інтервалів, які називаються кроками Троттера, тривалістю t/kt/k і апроксимувати еволюцію в часі на кожному інтервалі продуктною формулою порядку χ\chi SχS_{\chi}. Таким чином, PF порядку χ\chi для оператора еволюції в часі за kk кроків Троттера має вигляд:

Sχk(t)=[Sχ(tk)]k=eiHt+O(t(tk)χ) S_{\chi}^{k}(t) = \left[ S_{\chi} \left( \frac{t}{k} \right)\right]^k = e^{-i H t}+O\left(t \left( \frac{t}{k} \right)^{\chi} \right)

де член похибки зменшується зі збільшенням кількості кроків Троттера kk і порядку χ\chi PF.

Для заданого цілого числа k1k \geq 1 і продуктної формули Sχ(t)S_{\chi}(t), апроксимований стан з еволюцією в часі ρk(t)\rho_k(t) може бути отриманий з ρ0\rho_0 шляхом застосування kk ітерацій продуктної формули Sχ(tk)S_{\chi}\left(\frac{t}{k}\right).

ρk(t)=Sχ(tk)kρ0Sχ(tk)k\rho_k(t)=S_{\chi}\left(\frac{t}{k}\right)^k \rho_0 S_{\chi}\left(\frac{t}{k}\right)^{-k}

ρk(t)\rho_k(t) є апроксимацією для ρ(t)\rho(t) з похибкою апроксимації Троттера ||ρk(t)ρ(t)\rho_k(t)-\rho(t) ||. Якщо ми розглянемо лінійну комбінацію апроксимацій Троттера ρ(t)\rho(t):

μ(t)=jlxjρjkj(tkj)+деяка залишкова похибка Троттера,\mu(t) = \sum_{j}^{l} x_j \rho^{k_j}_{j}\left(\frac{t}{k_j}\right) + \text{деяка залишкова похибка Троттера} \, ,

де xjx_j є нашими ваговими коефіцієнтами, ρjkj\rho^{k_j}_j є матрицею густини, що відповідає чистому стану, отриманому шляхом еволюції початкового стану з продуктною формулою, SχkjS^{k_j}_{\chi}, що включає kjk_j кроків Троттера, і j1,...,lj \in {1, ..., l} індексує кількість PF, що складають MPF. Всі члени в μ(t)\mu(t) використовують ту саму продуктну формулу Sχ(t)S_{\chi}(t) як свою базу. Мета полягає в поліпшенні ||ρk(t)ρ(t)\rho_k(t)-\rho(t) \| шляхом знаходження μ(t)\mu(t) з ще меншим μ(t)ρ(t)\|\mu(t)-\rho(t)\|.

  • μ(t)\mu(t) не обов'язково має бути фізичним станом, оскільки xix_i не обов'язково мають бути додатними. Мета тут полягає в мінімізації похибки в очікуваному значенні спостережуваних, а не в знаходженні фізичної заміни для ρ(t)\rho(t).
  • kjk_j визначає як глибину кола, так і рівень апроксимації Троттера. Менші значення kjk_j призводять до коротших кіл, які зазнають менших похибок кола, але будуть менш точною апроксимацією бажаного стану.

Ключовим тут є те, що залишкова похибка Троттера, дана μ(t)\mu(t), є меншою, ніж похибка Троттера, яку можна було б отримати, просто використовуючи найбільше значення kjk_j.

Ви можете розглядати корисність цього з двох точок зору:

  1. Для фіксованого бюджету кроків Троттера, які Ви можете виконати, Ви можете отримати результати з похибкою Троттера, яка є меншою в загальному.
  2. Маючи деяку цільову кількість кроків Троттера, яка є занадто великою для виконання, Ви можете використовувати MPF для знаходження набору кіл меншої глибини для виконання, що призведе до подібної похибки Троттера.

Вимоги

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

  • Qiskit SDK v1.0 або пізніша версія, з підтримкою візуалізації
  • Qiskit Runtime v0.22 або пізніша версія (pip install qiskit-ibm-runtime)
  • Доповнення MPF Qiskit (pip install qiskit_addon_mpf)
  • Утиліти доповнень Qiskit (pip install qiskit_addon_utils)
  • Бібліотека Quimb (pip install quimb)
  • Бібліотека Qiskit Quimb (pip install qiskit-quimb)
  • Numpy v0.21 для сумісності між пакетами (pip install numpy==0.21)

Частина I. Приклад малого масштабу

Дослідження стабільності MPF

Немає очевидних обмежень щодо вибору кількості кроків Троттера kjk_j, які складають стан MPF μ(t)\mu(t). Однак їх потрібно обирати обережно, щоб уникнути нестабільностей в отриманих очікуваних значеннях, обчислених з μ(t)\mu(t). Хорошим загальним правилом є встановлення найменшого кроку Троттера kmink_{\text{min}} так, щоб t/kmin<1t/k_{\text{min}} \lt 1. Якщо Ви хочете дізнатися більше про це та про те, як вибрати інші значення kjk_j, зверніться до посібника Як вибрати кроки Троттера для MPF.

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

# Added by doQumentation — required packages for this notebook
!pip install -q matplotlib numpy qiskit qiskit-addon-mpf qiskit-addon-utils qiskit-aer qiskit-ibm-runtime rustworkx scipy
import numpy as np

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)

Для цього прикладу ми будемо використовувати стан Ніля як початковий стан Neel=0101...01\vert \text{Neel} \rangle = \vert 0101...01 \rangle та модель Гейзенберга на лінії з 10 вузлів для гамільтоніана, що керує еволюцією в часі

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

де JJ є силою зв'язку для ребер найближчих сусідів.

from qiskit.transpiler import CouplingMap
from rustworkx.visualization import graphviz_draw
from qiskit_addon_utils.problem_generators import generate_xyz_hamiltonian
import numpy as np

L = 10

# Generate some coupling map to use for this example
coupling_map = CouplingMap.from_line(L, bidirectional=False)
graphviz_draw(coupling_map.graph, method="circo")

# Get a qubit operator describing the Heisenberg field model
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])

Спостережуваним, яке ми будемо вимірювати, є намагніченість на парі кубітів в середині ланцюга.

from qiskit.quantum_info import SparsePauliOp

observable = SparsePauliOp.from_sparse_list(
[("ZZ", (L // 2 - 1, L // 2), 1.0)], num_qubits=L
)
print(observable)
SparsePauliOp(['IIIIZZIIII'],
coeffs=[1.+0.j])

Ми визначаємо прохід транспілятора для збору обертань XX та YY в колі як єдиний гейт XX+YY. Це дозволить нам використовувати властивості збереження спіну TeNPy під час обчислення MPO, значно прискорюючи обчислення.

from qiskit.circuit.library import XXPlusYYGate
from qiskit.transpiler import PassManager
from qiskit.transpiler.passes.optimization.collect_and_collapse import (
CollectAndCollapse,
collect_using_filter_function,
collapse_to_operation,
)
from functools import partial

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

Потім ми створюємо кола, що реалізують апроксимовані еволюції в часі Троттера.

from qiskit.synthesis import SuzukiTrotter
from qiskit_addon_utils.problem_generators import (
generate_time_evolution_circuit,
)
from qiskit import QuantumCircuit

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

Output of the previous code cell

Далі ми обчислюємо очікувані значення з еволюцією в часі з кіл Троттера.

from copy import deepcopy
from qiskit_aer import AerSimulator
from qiskit_ibm_runtime import EstimatorV2 as Estimator
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

aer_sim = AerSimulator()
estimator = Estimator(mode=aer_sim)

mpf_expvals_all_times, mpf_stds_all_times = [], []
for t, mpf_circuits in zip(trotter_times, all_circs):
mpf_expvals = []
circuits = [deepcopy(circuit) for circuit in mpf_circuits]
pm_sim = generate_preset_pass_manager(
backend=aer_sim, optimization_level=3
)
isa_circuits = pm_sim.run(circuits)
result = estimator.run(
[(circuit, observable) for circuit in isa_circuits], precision=0.005
).result()
mpf_expvals = [res.data.evs for res in result]
mpf_stds = [res.data.stds for res in result]
mpf_expvals_all_times.append(mpf_expvals)
mpf_stds_all_times.append(mpf_stds)

Ми також обчислюємо точні очікувані значення для порівняння.

from scipy.linalg import expm
from qiskit.quantum_info import Statevector

exact_expvals = []
for t in exact_evolution_times:
# Exact expectation values
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 є такими, де значення xjx_j не залежать від часу еволюції, tt. Розглянемо PF порядку χ=1\chi = 1 з kjk_j кроками Троттера, це можна записати як:

S1kj(tkj)=eiHt+n=1Antn+1kjnS_1^{k_j}\left( \frac{t}{k_j} \right)=e^{-i H t}+ \sum_{n=1}^{\infty} A_n \frac{t^{n+1}}{k_j^n}

де AnA_n є матрицями, які залежать від комутаторів членів FaF_a в розкладі гамільтоніана. Важливо зазначити, що самі AnA_n не залежать від часу та кількості кроків Троттера kjk_j. Тому можливо скасувати члени похибки нижчого порядку, що вносять внесок у μ(t)\mu(t), з ретельним вибором ваг xjx_j лінійної комбінації. Щоб скасувати похибку Троттера для перших l1l-1 членів (вони дадуть найбільші внески, оскільки вони відповідають меншій кількості кроків Троттера) у виразі для μ(t)\mu(t), коефіцієнти xjx_j повинні задовольняти наступні рівняння:

j=1lxj=1\sum_{j=1}^l x_j = 1 j=1l1xjkjn=0\sum_{j=1}^{l-1} \frac{x_j}{k_j^{n}} = 0

з n=0,...l2n=0, ... l-2. Перше рівняння гарантує, що немає зміщення в побудованому стані μ(t)\mu(t), тоді як друге рівняння забезпечує скасування похибок Троттера. Для PF вищого порядку друге рівняння стає j=1l1xjkjη=0\sum_{j=1}^{l-1} \frac{x_j}{k_j^{\eta}} = 0, де η=χ+2n\eta = \chi + 2n для симетричних PF і η=χ+n\eta = \chi + n в іншому випадку, з n=0,...,l2n=0, ..., l-2. Результуюча похибка (Посилання [1],[2]) тоді становить

ϵ=O(tl+1k1l). \epsilon = \mathcal{O} \left( \frac{t^{l+1}}{k_1^l} \right).

Визначення статичних коефіцієнтів MPF для заданого набору значень kjk_j зводиться до розв'язання лінійної системи рівнянь, визначеної двома рівняннями вище для змінних xjx_j: Ax=bAx=b. Де xx є нашими коефіцієнтами інтересу, AA є матрицею, яка залежить від kjk_j і типу PF, яку ми використовуємо (SS), а bb є вектором обмежень. Конкретно:

A0,j=1A_{0,j} = 1 Ai>0,j=kj(χ+s(i1))A_{i>0,j} = k_{j}^{-(\chi + s(i-1))} b0=1b_0 = 1 bi>0=0b_{i>0} = 0

де χ\chi є order, ss дорівнює 22, якщо symmetric має значення True, і 11 в іншому випадку, kjk_{j} є trotter_steps, а xx є змінними для розв'язання. Індекси ii і jj починаються з 00. Ми також можемо візуалізувати це в матричній формі:

A=[A0,0A0,1A0,2...A1,0A1,1A1,2...A2,0A2,1A2,2...............]=[111...k0(χ+s(11))k1(χ+s(11))k2(χ+s(11))...k0(χ+s(21))k1(χ+s(21))k2(χ+s(21))...............]A = \begin{bmatrix} A_{0,0} & A_{0,1} & A_{0,2} & ... \\ A_{1,0} & A_{1,1} & A_{1,2} & ... \\ A_{2,0} & A_{2,1} & A_{2,2} & ... \\ ... & ... & ... & ... \end{bmatrix} = \begin{bmatrix} 1 & 1 & 1 & ... \\ k_{0}^{-(\chi + s(1-1))} & k_{1}^{-(\chi + s(1-1))} & k_{2}^{-(\chi + s(1-1))} & ... \\ k_{0}^{-(\chi + s(2-1))} & k_{1}^{-(\chi + s(2-1))} & k_{2}^{-(\chi + s(2-1))} & ... \\ ... & ... & ... & ... \end{bmatrix}

і

b=[b0b1b2...]=[100...]b = \begin{bmatrix} b_{0} \\ b_{1} \\ b_{2} \\ ... \end{bmatrix} = \begin{bmatrix} 1 \\ 0 \\ 0 \\ ... \end{bmatrix}

Для більш детальної інформації зверніться до документації лінійної системи рівнянь (LSE).

Ми можемо знайти рішення для xx аналітично як x=A1bx = A^{-1}b; дивіться, наприклад, Посилання [1] або [2]. Однак це точне рішення може бути "погано обумовленим", що призводить до дуже великих L1-норм наших коефіцієнтів, xx, що може призвести до поганої продуктивності MPF. Замість цього можна також отримати апроксимоване рішення, яке мінімізує L1-норму xx, щоб спробувати оптимізувати поведінку MPF.

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

Тепер, коли ми вибрали наші значення kjk_j, ми повинні спочатку побудувати LSE, Ax=bAx=b, як пояснено вище. Матриця AA залежить не тільки від kjk_j, але й від нашого вибору PF, зокрема її порядку. Додатково, Ви можете врахувати, чи є PF симетричною чи ні (дивіться [1]), встановивши symmetric=True/False. Однак це не є обов'язковим, як показано у Посиланні [2].

from qiskit_addon_mpf.static import setup_static_lse

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

Давайте пропрацюємо значення, вибрані вище, щоб побудувати матрицю AA і вектор bb. З j=0,1,2j=0,1, 2 кроками Троттера kj=[1,2,4]k_j = [1, 2, 4], порядком χ=2\chi = 2 і вибором несиметричних кроків Троттера (s=1s=1), ми маємо, що матричні елементи AA нижче першого рядка визначаються виразом Ai>0,j=kj(2+1(i1))A_{i>0,j} = k_{j}^{-(2 + 1(i-1))}, конкретно:

A0,0=A0,1=A0,2=1A_{0,0} = A_{0,1} = A_{0,2} = 1 A1,j=kj1A1,0=112,  ,A1,1=122,  ,A1,2=142 A_{1,j} = k_{j}^{-1} \rightarrow A_{1,0} = \frac{1}{1^2}, \;, A_{1,1} = \frac{1}{2^2}, \;, A_{1,2} = \frac{1}{4^2} A2,j=kj2A2,0=113,  ,A2,1=123,  ,A2,2=143 A_{2,j} = k_{j}^{-2} \rightarrow A_{2,0} = \frac{1}{1^3}, \;, A_{2,1} = \frac{1}{2^3}, \;, A_{2,2} = \frac{1}{4^3}

або в матричній формі:

A=[11111221421123143]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}

Це можна побачити, перевіривши об'єкт lse:

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

Тоді як вектор обмежень bb має наступні елементи: b0=1b_{0} = 1 b1=b2=0b_1 = b_2 = 0

Таким чином,

b=[100]b = \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix}

І аналогічно в lse:

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

Об'єкт lse має методи для знаходження статичних коефіцієнтів xjx_j, які задовольняють систему рівнянь.

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]
Оптимізація для xx з використанням точної моделі

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

У наступному розділі стане зрозуміло, чому існує цей інтерфейс.

from qiskit_addon_mpf.costs import setup_exact_problem

model_exact, coeffs_exact = setup_exact_problem(lse)
model_exact.solve()
print(coeffs_exact.value)
[ 0.04761905 -0.57142857  1.52380952]

Як індикатор того, чи MPF, побудований з цими коефіцієнтами, дасть хороші результати, ми можемо використовувати L1-норму (див. також Посилання [1]).

print(
"L1 norm of the exact coefficients:",
np.linalg.norm(coeffs_exact.value, ord=1),
) # ord specifies the norm. ord=1 is for L1
L1 norm of the exact coefficients: 2.1428571428556378
Оптимізація для xx з використанням наближеної моделі

Може статися, що L1-норма для обраного набору значень kjk_j вважається занадто високою. Якщо це так, і Ви не можете вибрати інший набір значень kjk_j, Ви можете використовувати наближений розв'язок LSE замість точного.

Для цього просто використовуйте setup_approximate_model для побудови іншого екземпляра cvxpy.Problem, який обмежує L1-норму до обраного порогу, мінімізуючи різницю між AxAx і bb.

from qiskit_addon_mpf.costs import setup_sum_of_squares_problem

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, позначений як μS(t)\mu^S(t), не є оптимальною проекцією ρ(t)\rho(t) на підпростір, натягнутий станами формули добутку {ρki(t)}i=1r\{\rho_{k_i}(t)\}_{i=1}^r.

Щоб вирішити це, ми розглядаємо динамічний MPF (введений у Посиланні [2] і експериментально продемонстрований у Посиланні [3]), який дійсно мінімізує похибку апроксимації в нормі Фробеніуса. Формально, ми зосереджуємося на мінімізації

ρ(t)μD(t)F2  =  Tr[(ρ(t)μD(t))2],\|\rho(t) - \mu^D(t)\|_F^2 \;=\; \mathrm{Tr}\bigl[ \left( \rho(t) - \mu^D(t)\right)^2 \bigr],

відносно деяких коефіцієнтів xi(t)x_i(t) у кожний момент часу tt. Оптимальний проектор у нормі Фробеніуса тоді є μD(t)  =  i=1rxi(t)ρki(t)\mu^D(t) \;=\; \sum_{i=1}^r x_i(t)\,\rho_{k_i}(t), і ми називаємо μD(t)\mu^D(t) динамічним MPF. Підставляючи визначення вище:

ρ(t)μD(t)F2  =  =Tr[(ρ(t)μD(t))2]  =  =Tr[(ρ(t)i=1rxi(t)ρki(t))(ρ(t)j=1rxj(t)ρkj(t))]  =  =1  +  i,j=1rMi,j(t)xi(t)xj(t)    2i=1rLiexact(t)xi(t),\|\rho(t) - \mu^D(t)\|_F^2 \;=\; \\ = \mathrm{Tr}\bigl[ \left( \rho(t) - \mu^D(t)\right)^2 \bigr] \;=\; \\ = \mathrm{Tr}\bigl[ \left( \rho(t) - \sum_{i=1}^r x_i(t)\,\rho_{k_i}(t) \right) \left( \rho(t) - \sum_{j=1}^r x_j(t)\,\rho_{k_j}(t) \right) \bigr] \;=\; \\ = 1 \;+\; \sum_{i,j=1}^r M_{i,j}(t)\,x_i(t)\,x_j(t) \;-\; 2 \sum_{i=1}^r L_i^{\mathrm{exact}}(t)\,x_i(t),

де Mi,j(t)M_{i,j}(t) є матрицею Грама, визначеною як

Mi,j(t)  =  Tr[ρki(t)ρkj(t)]  =  ψin ⁣S(t/ki)kiS(t/kj)kj ⁣ψin2.M_{i,j}(t) \;=\; \mathrm{Tr}\bigl[\rho_{k_i}(t)\,\rho_{k_j}(t)\bigr] \;=\; \bigl|\langle \psi_{\mathrm{in}} \!\mid S\bigl(t/k_i\bigr)^{-k_i}\,S\bigl(t/k_j\bigr)^{k_j} \!\mid \psi_{\mathrm{in}} \rangle \bigr|^2.

і

Liexact(t)=Tr[ρ(t)ρki(t)]L_i^{\mathrm{exact}}(t) = \mathrm{Tr}[\rho(t)\,\rho_{k_i}(t)]

представляє перекриття між точним станом ρ(t)\rho(t) і кожною апроксимацією формули добутку ρki(t)\rho_{k_i}(t). У практичних сценаріях ці перекриття можуть вимірюватися лише наближено через шум або частковий доступ до ρ(t)\rho(t).

Тут ψin\lvert\psi_{\mathrm{in}}\rangle є початковим станом, а S()S(\cdot) є операцією, застосованою у формулі добутку. Вибираючи коефіцієнти xi(t)x_i(t), які мінімізують цей вираз (і обробляючи наближені дані перекриття, коли ρ(t)\rho(t) не повністю відомо), ми отримуємо "найкращу" (у сенсі норми Фробеніуса) динамічну апроксимацію ρ(t)\rho(t) у підпросторі MPF. Величини Li(t)L_i(t) і Mi,j(t)M_{i,j}(t) можуть бути обчислені ефективно з використанням методів тензорних мереж [3]. Доповнення MPF Qiskit надає кілька "бекендів" для виконання обчислень. Наведений нижче приклад показує найбільш гнучкий спосіб зробити це, а документація бекенду на основі шарів TeNPy також пояснює це дуже детально. Щоб використовувати цей метод, почніть з ланцюга, що реалізує бажану еволюцію в часі, і створіть моделі, які представляють ці операції з шарів відповідного ланцюга. Нарешті, створюється об'єкт Evolver, який можна використовувати для генерації величин, що еволюціонують у часі, Mi,j(t)M_{i,j}(t) і Li(t)L_i(t). Починаємо зі створення об'єкта Evolver, що відповідає наближеній еволюції в часі (ApproxEvolverFactory), реалізованій ланцюгами. Зокрема, приділіть особливу увагу змінній order, щоб вони збігалися. Зверніть увагу, що при генерації ланцюгів, що відповідають наближеній еволюції в часі, ми використовуємо значення-заповнювачі для time = 1.0 і кількості кроків Троттера (reps=1). Правильні апроксимуючі ланцюги потім створюються розв'язувачем динамічної задачі в setup_dynamic_lse.

from qiskit_addon_utils.slicing import slice_by_depth
from qiskit_addon_mpf.backends.tenpy_layers import LayerModel
from qiskit_addon_mpf.backends.tenpy_layers import LayerwiseEvolver
from functools import partial

# 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, які визначають деталі симуляції тензорної мережі, повинні бути обрані ретельно, щоб уникнути налаштування некоректно визначеної оптимізаційної задачі.

Потім ми налаштовуємо точний еволвер (наприклад, ExactEvolverFactory), який повертає об'єкт Evolver, що обчислює істинну або "референсну" еволюцію в часі. Реалістично, ми б апроксимували точну еволюцію, використовуючи формулу Сузукі-Троттера вищого порядку або інший надійний метод з малим кроком часу. Нижче ми апроксимуємо точний стан, що еволюціонує в часі, формулою Сузукі-Троттера четвертого порядку, використовуючи малий крок часу dt=0.1, що означає, що кількість кроків Троттера, використаних у момент часу tt, становить k=t/dtk=t/dt. Ми також вказуємо деякі специфічні для TeNPy опції усічення для обмеження максимальної розмірності зв'язку базової тензорної мережі, а також мінімальних сингулярних значень розділених зв'язків тензорної мережі. Ці параметри можуть впливати на точність очікуваного значення, обчисленого з динамічними коефіцієнтами MPF, тому важливо дослідити діапазон значень для знаходження оптимального балансу між обчислювальним часом і точністю. Зверніть увагу, що обчислення коефіцієнтів MPF не залежить від очікуваного значення PF, отриманого з апаратного виконання, і тому його можна налаштовувати в пост-обробці.

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,
},
)

Далі створіть початковий стан Вашої системи у форматі, сумісному з TeNPy (наприклад, MPS_neel_state=0101...01\vert 0101...01 \rangle). Це налаштовує багаточастинкову хвильову функцію, яку Ви будете еволюціонувати в часі ψin\lvert\psi_{\mathrm{in}}\rangle як тензор.

from qiskit_addon_mpf.backends.tenpy_tebd import MPOState
from qiskit_addon_mpf.backends.tenpy_tebd import MPS_neel_state

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

mps_initial_state = MPS_neel_state(models[0].lat)

Для кожного кроку часу tt ми налаштовуємо динамічну лінійну систему рівнянь за допомогою методу setup_dynamic_lse. Відповідний об'єкт містить інформацію про динамічну задачу MPF: lse.A дає матрицю Грама MM, тоді як lse.b дає перекриття LL. Потім ми можемо розв'язати LSE (якщо вона не є некоректно визначеною) для знаходження динамічних коефіцієнтів, використовуючи setup_frobenius_problem. Важливо відзначити різницю зі статичними коефіцієнтами, які залежать лише від деталей використаної формули добутку і не залежать від деталей еволюції в часі (гамільтоніана і початкового стану).

from qiskit_addon_mpf.dynamic import setup_dynamic_lse
from qiskit_addon_mpf.costs import setup_frobenius_problem

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

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

import matplotlib.pyplot as plt

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
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",
color="orange",
label="Static MPF - Approximate",
)

# # 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",
color="pink",
label="Dynamic MPF",
)

plt.plot(
exact_evolution_times,
exact_expvals,
lw=3,
color="red",
label="Exact time-evolution",
)

plt.title(
f"Expectation values for (ZZ,{(L//2-1, L//2)}) as a function of time"
)
plt.xlabel("Time")
plt.ylabel("Expectation Value")
plt.legend()
plt.grid()

Output of the previous code cell

У випадках, подібних до наведеного вище прикладу, де k=1k=1 PF поводиться погано у всі моменти часу, якість результатів динамічного MPF також суттєво постраждає. У таких ситуаціях корисно дослідити можливість використання окремих PF з вищою кількістю кроків Троттера для покращення загальної якості результатів. У цих симуляціях ми бачимо взаємодію різних типів похибок: похибки від кінцевого семплювання і похибки Троттера від формул добутку. MPF допомагає зменшити похибку Троттера через формули добутку, але має вищу похибку семплювання порівняно з формулами добутку. Це може бути вигідним, оскільки формули добутку можуть зменшити похибку семплювання зі збільшенням семплювання, але систематична похибка через апроксимацію Троттера залишається недоторканою.

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

Тепер розглянемо один момент часу t=1.0t=1.0 і обчислимо очікуване значення намагніченості різними методами, використовуючи один QPU. Конкретний вибір tt було зроблено так, щоб максимізувати різницю між різними методами та спостерігати їхню відносну ефективність. Щоб визначити вікно часу, для якого динамічний MPF гарантовано виробляє спостережувані величини з меншою похибкою, ніж будь-які з окремих формул Троттера в межах мультипродукту, ми можемо реалізувати "MPF тест" - див. рівняння (17) та навколишній текст у [3].

Налаштування схем Троттера

На цьому етапі ми знайшли наші коефіцієнти розкладання, xx, і все, що залишилось зробити - це згенерувати троттеризовані квантові схеми. Знову, модуль qiskit_addon_utils.problem_generators приходить на допомогу з корисною функцією для цього:

from qiskit.synthesis import SuzukiTrotter
from qiskit_addon_utils.problem_generators import (
generate_time_evolution_circuit,
)
from qiskit import QuantumCircuit

total_time = 1.0
mpf_circuits = []
for k in mpf_trotter_steps:
# Initial Neel state preparation
circuit = QuantumCircuit(L)
circuit.x([i for i in range(L) if i % 2 != 0])

trotter_circ = generate_time_evolution_circuit(
hamiltonian,
synthesis=SuzukiTrotter(order=order, reps=k),
time=total_time,
)

circuit.compose(trotter_circ, inplace=True)

mpf_circuits.append(pm.run(circuit))
mpf_circuits[-1].draw("mpl", fold=-1, scale=0.4)

Output of the previous code cell

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

Повернемось до обчислення очікуваного значення для одного моменту часу. Ми виберемо бекенд для виконання експерименту на обладнанні.

from qiskit_ibm_runtime import QiskitRuntimeService

service = QiskitRuntimeService()
backend = service.least_busy(min_num_qubits=127)
print(backend)

qubits = list(range(backend.num_qubits))

Потім ми видаляємо кубіти-викиди з мапи зв'язків, щоб забезпечити, що етап розміщення транспайлера не включає їх. Нижче ми використовуємо повідомлені властивості бекенду, збережені в об'єкті target, і видаляємо кубіти, які мають або помилку вимірювання, або двокубітний гейт вище певного порогу (max_meas_err, max_twoq_err), або час T2T_2 (який визначає втрату когерентності) нижче певного порогу (min_t2).

import copy
from qiskit.transpiler import Target, CouplingMap

target = backend.target
instruction_2q = "cz"

cmap = target.build_coupling_map(filter_idle_qubits=True)
cmap_list = list(cmap.get_edges())

max_meas_err = 0.012
min_t2 = 40
max_twoq_err = 0.005

# Remove qubits with bad measurement or t2
cust_cmap_list = copy.deepcopy(cmap_list)
for q in range(target.num_qubits):
meas_err = target["measure"][(q,)].error
if target.qubit_properties[q].t2 is not None:
t2 = target.qubit_properties[q].t2 * 1e6
else:
t2 = 0
if meas_err > max_meas_err or t2 < min_t2:
# print(q)
for q_pair in cmap_list:
if q in q_pair:
try:
cust_cmap_list.remove(q_pair)
except ValueError:
continue

# Remove qubits with bad 2q gate or t2
for q in cmap_list:
twoq_gate_err = target[instruction_2q][q].error
if twoq_gate_err > max_twoq_err:
# print(q)
for q_pair in cmap_list:
if q == q_pair:
try:
cust_cmap_list.remove(q_pair)
except ValueError:
continue

cust_cmap = CouplingMap(cust_cmap_list)

cust_target = Target.from_configuration(
basis_gates=backend.configuration().basis_gates
+ ["measure"], # or whatever new set of gates
coupling_map=cust_cmap,
)

sorted_components = sorted(
[list(comp.physical_qubits) for comp in cust_cmap.connected_components()],
reverse=True,
)
print("size of largest component", len(sorted_components[0]))
size of largest component 10

Ми хочемо встановити max_meas_err, min_t2 та max_twoq_err таким чином, щоб знайти достатньо велику підмножину кубітів, яка підтримує виконання схеми. У нашому випадку достатньо знайти 10-кубітний 1D ланцюг.

cust_cmap.draw()

Output of the previous code cell

Потім ми можемо відобразити схему та спостережувану величину на фізичні кубіти пристрою.

from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

transpiler = generate_preset_pass_manager(
optimization_level=3, target=cust_target
)

transpiled_circuits = [transpiler.run(circ) for circ in mpf_circuits]

qubits_layouts = [
[
idx
for idx, qb in circuit.layout.initial_layout.get_physical_bits().items()
if qb._register.name != "ancilla"
]
for circuit in transpiled_circuits
]

transpiled_circuits = []
for circuit, layout in zip(mpf_circuits, qubits_layouts):
transpiler = generate_preset_pass_manager(
optimization_level=3, backend=backend, initial_layout=layout
)
transpiled_circuit = transpiler.run(circuit)
transpiled_circuits.append(transpiled_circuit)

# transform the observable defined on virtual qubits to
# an observable defined on all physical qubits
isa_observables = [
observable.apply_layout(circ.layout) for circ in transpiled_circuits
]
print(transpiled_circuits[-1].depth(lambda x: x.operation.num_qubits == 2))
print(transpiled_circuits[-1].count_ops())
transpiled_circuits[-1].draw("mpl", idle_wires=False, fold=False)
51
OrderedDict([('sx', 310), ('rz', 232), ('cz', 132), ('x', 19)])

Output of the previous code cell

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

За допомогою примітиву Estimator ми можемо отримати оцінку очікуваного значення від QPU. Ми виконуємо оптимізовані AQC схеми з додатковими методами пом'якшення та придушення помилок.

from qiskit_ibm_runtime import EstimatorV2 as Estimator

estimator = Estimator(mode=backend)
estimator.options.default_shots = 30000

# Set simple error suppression/mitigation options
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, 3, 5)
estimator.options.resilience.zne.extrapolator = ("exponential", "linear")

estimator.options.environment.job_tags = ["mpf small"]

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

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

Єдиним етапом постобробки є комбінування очікуваного значення, отриманого від примітивів Qiskit Runtime на різних кроках Троттера, використовуючи відповідні коефіцієнти MPF. Для спостережуваної величини AA ми маємо:

Ampf=Tr[Aμ(t)]=jxjTr[Aρkj]=jxjAj \langle A \rangle_{\text{mpf}} = \text{Tr} [A \mu(t)] = \sum_{j} x_j \text{Tr} [A \rho_{k_j}] = \sum_{j} x_j \langle A \rangle_j

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

result_exp = job.result()
evs_exp = [res.data.evs for res in result_exp]
evs_std = [res.data.stds for res in result_exp]

print(evs_exp)
[array(-0.06361607), array(-0.23820448), array(-0.50271805)]

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

exact_mpf_std = np.sqrt(
sum(
[
(coeff**2) * (std**2)
for coeff, std in zip(coeffs_exact.value, evs_std)
]
)
)
print(
"Exact static MPF expectation value: ",
evs_exp @ coeffs_exact.value,
"+-",
exact_mpf_std,
)
approx_mpf_std = np.sqrt(
sum(
[
(coeff**2) * (std**2)
for coeff, std in zip(coeffs_approx.value, evs_std)
]
)
)
print(
"Approximate static MPF expectation value: ",
evs_exp @ coeffs_approx.value,
"+-",
approx_mpf_std,
)
dynamic_mpf_std = np.sqrt(
sum(
[
(coeff**2) * (std**2)
for coeff, std in zip(mpf_dynamic_coeffs_list[7], evs_std)
]
)
)
print(
"Dynamic MPF expectation value: ",
evs_exp @ mpf_dynamic_coeffs_list[7],
"+-",
dynamic_mpf_std,
)
Exact static MPF expectation value:  -0.6329590442738475 +- 0.012798249760406036
Approximate static MPF expectation value: -0.5690390035339492 +- 0.010459559917168473
Dynamic MPF expectation value: -0.4655579758795695 +- 0.007639139186720507

Нарешті, для цієї малої задачі ми можемо обчислити точне еталонне значення, використовуючи scipy.linalg.expm наступним чином:

from scipy.linalg import expm
from qiskit.quantum_info import Statevector

exp_H = expm(-1j * total_time * hamiltonian.to_matrix())

initial_state_circuit = QuantumCircuit(L)
initial_state_circuit.x([i for i in range(L) if i % 2 != 0])
initial_state = Statevector(initial_state_circuit).data

time_evolved_state = exp_H @ initial_state

exact_obs = (
time_evolved_state.conj() @ observable.to_matrix() @ time_evolved_state
)
print("Exact expectation value ", exact_obs.real)
Exact expectation value  -0.39909900734489434
sym = {1: "^", 2: "s", 4: "p"}
# Get expectation values at all times for each Trotter step
for k, step in enumerate(mpf_trotter_steps):
plt.errorbar(
k,
evs_exp[k],
yerr=evs_std[k],
alpha=0.5,
markersize=4,
marker=sym[step],
color="grey",
label=f"{mpf_trotter_steps[k]} Trotter steps",
) # , , )

plt.errorbar(
3,
evs_exp @ coeffs_exact.value,
yerr=exact_mpf_std,
markersize=4,
marker="o",
color="purple",
label="Static MPF",
)

plt.errorbar(
4,
evs_exp @ coeffs_approx.value,
yerr=approx_mpf_std,
markersize=4,
marker="o",
color="orange",
label="Approximate static MPF",
)

plt.errorbar(
5,
evs_exp @ mpf_dynamic_coeffs_list[7],
yerr=dynamic_mpf_std,
markersize=4,
marker="o",
color="pink",
label="Dynamic MPF",
)

plt.axhline(
y=exact_obs.real,
linestyle="--",
color="red",
label="Exact time-evolution",
)

plt.title(
f"Expectation values for (ZZ,{(L//2-1, L//2)}) 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()

Output of the previous code cell

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

def relative_error(ev, exact_ev):
return abs(ev - exact_ev)

relative_error_k = [relative_error(ev, exact_obs.real) for ev in evs_exp]
relative_error_mpf = relative_error(evs_exp @ mpf_coeffs, exact_obs.real)
relative_error_approx_mpf = relative_error(
evs_exp @ coeffs_approx.value, exact_obs.real
)
relative_error_dynamic_mpf = relative_error(
evs_exp @ mpf_dynamic_coeffs_list[7], exact_obs.real
)

print("relative error for each trotter steps", relative_error_k)
print("relative error with MPF exact coeffs", relative_error_mpf)
print("relative error with MPF approx coeffs", relative_error_approx_mpf)
print("relative error with MPF dynamic coeffs", relative_error_dynamic_mpf)
relative error for each trotter steps [0.33548293650112293, 0.16089452939226306, 0.10361904247828346]
relative error with MPF exact coeffs 0.2338600369291003
relative error with MPF approx coeffs 0.16993999618905486
relative error with MPF dynamic coeffs 0.06645896853467514

Частина II: масштабування

Давайте масштабуємо задачу за межі того, що можливо точно змодулювати. У цьому розділі ми зосередимося на відтворенні деяких результатів, показаних у посиланні [3].

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

Гамільтоніан

Для великомасштабного прикладу ми використовуємо модель XXZ на лінії з 50 вузлів:

H^XXZ=i=1L1Ji,(i+1)(XiX(i+1)+YiY(i+1)+2ZiZ(i+1)),\hat{\mathcal{H}}_{XXZ} = \sum_{i=1}^{L-1} J_{i,(i+1)}\left(X_i X_{(i+1)}+Y_i Y_{(i+1)}+ 2\cdot Z_i Z_{(i+1)} \right) \, ,

де Ji,(i+1)J_{i,(i+1)} є випадковим коефіцієнтом, що відповідає ребру (i,i+1)(i, i+1). Це гамільтоніан, розглянутий у демонстрації, представленій у посиланні [3].

L = 50
# Generate some coupling map to use for this example
coupling_map = CouplingMap.from_line(L, bidirectional=False)
graphviz_draw(coupling_map.graph, method="circo")

Output of the previous code cell

import numpy as np
from qiskit.quantum_info import SparsePauliOp, Pauli

# Generate random coefficients for our XXZ Hamiltonian
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,
)

print(hamiltonian)
SparsePauliOp(['IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXX', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYY', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZ', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'XXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'YYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'ZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXI', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYI', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZI', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII'],
coeffs=[1. +0.j, 2.09762701+0.j, 2.09762701+0.j, 4.19525402+0.j,
2.43037873+0.j, 2.43037873+0.j, 4.86075747+0.j, 2.20552675+0.j,
2.20552675+0.j, 4.4110535 +0.j, 2.08976637+0.j, 2.08976637+0.j,
4.17953273+0.j, 1.8473096 +0.j, 1.8473096 +0.j, 3.6946192 +0.j,
2.29178823+0.j, 2.29178823+0.j, 4.58357645+0.j, 1.87517442+0.j,
1.87517442+0.j, 3.75034885+0.j, 2.783546 +0.j, 2.783546 +0.j,
5.567092 +0.j, 2.92732552+0.j, 2.92732552+0.j, 5.85465104+0.j,
1.76688304+0.j, 1.76688304+0.j, 3.53376608+0.j, 2.58345008+0.j,
2.58345008+0.j, 5.16690015+0.j, 2.05778984+0.j, 2.05778984+0.j,
4.11557968+0.j, 2.13608912+0.j, 2.13608912+0.j, 4.27217824+0.j,
2.85119328+0.j, 2.85119328+0.j, 5.70238655+0.j, 1.14207212+0.j,
1.14207212+0.j, 2.28414423+0.j, 1.1742586 +0.j, 1.1742586 +0.j,
2.3485172 +0.j, 1.04043679+0.j, 1.04043679+0.j, 2.08087359+0.j,
2.66523969+0.j, 2.66523969+0.j, 5.33047938+0.j, 2.5563135 +0.j,
2.5563135 +0.j, 5.112627 +0.j, 2.7400243 +0.j, 2.7400243 +0.j,
5.48004859+0.j, 2.95723668+0.j, 2.95723668+0.j, 5.91447337+0.j,
2.59831713+0.j, 2.59831713+0.j, 5.19663426+0.j, 1.92295872+0.j,
1.92295872+0.j, 3.84591745+0.j, 2.56105835+0.j, 2.56105835+0.j,
5.12211671+0.j, 1.23654885+0.j, 1.23654885+0.j, 2.4730977 +0.j,
2.27984204+0.j, 2.27984204+0.j, 4.55968409+0.j, 1.28670657+0.j,
1.28670657+0.j, 2.57341315+0.j, 2.88933783+0.j, 2.88933783+0.j,
5.77867567+0.j, 2.04369664+0.j, 2.04369664+0.j, 4.08739329+0.j,
1.82932388+0.j, 1.82932388+0.j, 3.65864776+0.j, 1.52911122+0.j,
1.52911122+0.j, 3.05822245+0.j, 2.54846738+0.j, 2.54846738+0.j,
5.09693476+0.j, 1.91230066+0.j, 1.91230066+0.j, 3.82460133+0.j,
2.1368679 +0.j, 2.1368679 +0.j, 4.2737358 +0.j, 1.0375796 +0.j,
1.0375796 +0.j, 2.0751592 +0.j, 2.23527099+0.j, 2.23527099+0.j,
4.47054199+0.j, 2.22419145+0.j, 2.22419145+0.j, 4.44838289+0.j,
2.23386799+0.j, 2.23386799+0.j, 4.46773599+0.j, 2.88749616+0.j,
2.88749616+0.j, 5.77499231+0.j, 2.3636406 +0.j, 2.3636406 +0.j,
4.7272812 +0.j, 1.7190158 +0.j, 1.7190158 +0.j, 3.4380316 +0.j,
1.87406391+0.j, 1.87406391+0.j, 3.74812782+0.j, 2.39526239+0.j,
2.39526239+0.j, 4.79052478+0.j, 1.12045094+0.j, 1.12045094+0.j,
2.24090189+0.j, 2.33353343+0.j, 2.33353343+0.j, 4.66706686+0.j,
2.34127574+0.j, 2.34127574+0.j, 4.68255148+0.j, 1.42076512+0.j,
1.42076512+0.j, 2.84153024+0.j, 1.2578526 +0.j, 1.2578526 +0.j,
2.51570519+0.j, 1.6308567 +0.j, 1.6308567 +0.j, 3.2617134 +0.j])

Для спостережуваної величини ми вибираємо Z24Z25Z_{24}Z_{25}, як показано на нижній панелі Рис. 5 посилання [3].

observable = SparsePauliOp.from_sparse_list(
[("ZZ", (L // 2 - 1, L // 2), 1.0)], num_qubits=L
)
print(observable)
SparsePauliOp(['IIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIII'],
coeffs=[1.+0.j])

Вибір кроків Троттера

Експеримент, продемонстрований на Рис. 4 посилання [3], використовує kj=[2,3,4]k_j = [2, 3, 4] симетричні кроки Троттера порядку 22. Ми зосереджуємося на результатах для часу t=3t=3, де MPF і PF з більшою кількістю кроків Троттера (6 у цьому випадку) мають однакову похибку Троттера. Однак, очікуване значення MPF обчислюється з ланцюгів, що відповідають меншій кількості кроків Троттера і тому є більш мілкими. На практиці, навіть якщо MPF і ланцюг з більш глибокими кроками Троттера мають однакову похибку Троттера, ми очікуємо, що експериментальне очікуване значення, обчислене з ланцюгів MPF, буде ближчим до теоретичного, оскільки воно передбачає виконання більш мілких ланцюгів, менш схильних до апаратного шуму порівняно з ланцюгом, що відповідає PF з більшою кількістю кроків Троттера.

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

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

Тут ми розглядаємо статичні коефіцієнти MPF для цієї задачі.

lse = setup_static_lse(mpf_trotter_steps, order=order, symmetric=symmetric)
mpf_coeffs = lse.solve()
print(
f"The static coefficients associated with the ansatze are: {mpf_coeffs}"
)
print("L1 norm:", np.linalg.norm(mpf_coeffs, ord=1))
The static coefficients associated with the ansatze are: [ 0.26666667 -2.31428571  3.04761905]
L1 norm: 5.628571428571431
model_approx, coeffs_approx = setup_sum_of_squares_problem(
lse, max_l1_norm=2.0
)
model_approx.solve()
print(coeffs_approx.value)
print(
"L1 norm of the approximate coefficients:",
np.linalg.norm(coeffs_approx.value, ord=1),
)
[-0.24255546 -0.25744454  1.5       ]
L1 norm of the approximate coefficients: 2.0

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

# 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": 4,
},
)

# Create exact time-evolution circuits
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)
]

# Create the time-evolution object
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)

lse = setup_dynamic_lse(
mpf_trotter_steps,
total_time,
identity_factory,
exact_factory,
approx_factory,
mps_initial_state,
)
problem, coeffs = setup_frobenius_problem(lse)
try:
problem.solve()
mpf_dynamic_coeffs = coeffs.value
except Exception as error:
print(error, "Calculation Failed for time", total_time)
print("")

Побудова кожного з ланцюгів Троттера в нашому розкладі MPF

from qiskit.synthesis import SuzukiTrotter
from qiskit_addon_utils.problem_generators import (
generate_time_evolution_circuit,
)
from qiskit import QuantumCircuit

mpf_circuits = []
for k in mpf_trotter_steps:
# Initial state preparation |1010..>
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)

Побудова ланцюга Троттера з порівнянною похибкою Троттера до MPF

k = 6

# Initial state preparation |1010..>
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=k, order=order),
time=total_time,
)

comp_circuit.compose(trotter_circ, qubits=range(L), inplace=True)

mpf_circuits.append(comp_circuit)

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

import copy
from qiskit.transpiler import Target, CouplingMap

target = backend.target
instruction_2q = "cz"

cmap = target.build_coupling_map(filter_idle_qubits=True)
cmap_list = list(cmap.get_edges())

max_meas_err = 0.055
min_t2 = 30
max_twoq_err = 0.01

# Remove qubits with bad measurement or t2
cust_cmap_list = copy.deepcopy(cmap_list)
for q in range(target.num_qubits):
meas_err = target["measure"][(q,)].error
if target.qubit_properties[q].t2 is not None:
t2 = target.qubit_properties[q].t2 * 1e6
else:
t2 = 0
if meas_err > max_meas_err or t2 < min_t2:
# print(q)
for q_pair in cmap_list:
if q in q_pair:
try:
cust_cmap_list.remove(q_pair)
except ValueError:
continue

# Remove qubits with bad 2q gate or t2
for q in cmap_list:
twoq_gate_err = target[instruction_2q][q].error
if twoq_gate_err > max_twoq_err:
# print(q)
for q_pair in cmap_list:
if q == q_pair:
try:
cust_cmap_list.remove(q_pair)
except ValueError:
continue

cust_cmap = CouplingMap(cust_cmap_list)

cust_target = Target.from_configuration(
basis_gates=backend.configuration().basis_gates
+ ["measure"], # or whatever new set of gates
coupling_map=cust_cmap,
)

sorted_components = sorted(
[list(comp.physical_qubits) for comp in cust_cmap.connected_components()],
reverse=True,
)
print("size of largest component", len(sorted_components[0]))
size of largest component 73
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

transpiler = generate_preset_pass_manager(
optimization_level=3, target=cust_target
)

transpiled_circuits = [transpiler.run(circ) for circ in mpf_circuits]

qubits_layouts = [
[
idx
for idx, qb in circuit.layout.initial_layout.get_physical_bits().items()
if qb._register.name != "ancilla"
]
for circuit in transpiled_circuits
]

transpiled_circuits = []
for circuit, layout in zip(mpf_circuits, qubits_layouts):
transpiler = generate_preset_pass_manager(
optimization_level=3, backend=backend, initial_layout=layout
)
transpiled_circuit = transpiler.run(circuit)
transpiled_circuits.append(transpiled_circuit)

# transform the observable defined on virtual qubits to
# an observable defined on all physical qubits
isa_observables = [
observable.apply_layout(circ.layout) for circ in transpiled_circuits
]

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

from qiskit_ibm_runtime import EstimatorV2 as Estimator

estimator = Estimator(mode=backend)
estimator.options.default_shots = 30000

# Set simple error suppression/mitigation options
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 = ["mpf large"]

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

Крок 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.08034071), array(-0.00605026), array(-0.15345759), array(-0.18127293)]
[array(0.04482517), array(0.03438413), array(0.21540776), array(0.21520829)]
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.47510243192011536 +- 0.6613940032465087
Approximate static MPF expectation value: -0.20914170384216998 +- 0.32341567460419135
Dynamic MPF expectation value: -0.07994951978722761 +- 0.07423091963310202
sym = {2: "^", 3: "s", 4: "p"}
# Get expectation values at all times for each Trotter step
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="6 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"Expectation values for (ZZ,{(L//2-1, L//2)}) 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()

Output of the previous code cell

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

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

Також важливо зазначити, що в прикладі, який відтворює експеримент у посиланні [3], часова точка t=3t=3 виходить за межі, при якій очікується, що PF з k=2k=2 поводиться добре, що становить t/k>1t/k>1, як обговорюється в цьому посібнику.

Посилання

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

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

[3] Kim, Y., Eddins, A., Anand, S., Wei, K. X., Van Den Berg, E., Rosenblatt, S., ... & Kandala, A. (2023). Evidence for the utility of quantum computing before fault tolerance. Nature, 618(7965), 500-505.

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