Formula berbilang produk untuk mengurangkan ralat Trotter

Anggaran penggunaan QPU: Empat minit pada pemproses Heron r2 (NOTA: Ini hanyalah anggaran sahaja. Masa jalan sebenar anda mungkin berbeza.)

Latar belakang

Tutorial ini menunjukkan cara menggunakan Formula Berbilang Produk (MPF) untuk mencapai ralat Trotter yang lebih rendah pada observable kita berbanding ralat yang ditanggung oleh Circuit Trotter terdalam yang akan kita jalankan sebenarnya. MPF mengurangkan ralat Trotter bagi dinamik Hamiltonian melalui gabungan berwajaran bagi beberapa pelaksanaan Circuit. Pertimbangkan tugas mencari nilai jangkaan observable untuk keadaan kuantum $\rho(t)=e^{-i H t} \rho(0) e^{i H t}$ dengan Hamiltonian $H$. Kita boleh menggunakan Formula Produk (PF) untuk menghampiri evolusi masa $e^{-i H t}$ dengan cara berikut:

• Tulis Hamiltonian $H$ sebagai $H=\sum_{a=1}^d F_a,$ di mana $F_a$ ialah pengoperasi Hermitian supaya setiap unitary yang sepadan boleh dilaksanakan dengan cekap pada peranti kuantum.
• Hampirkan sebutan $F_a$ yang tidak bertukar ganti antara satu sama lain.

Kemudian, PF peringkat pertama (formula Lie-Trotter) ialah:

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

yang mempunyai sebutan ralat kuadratik $S_1(t)=e^{-i H t}+\mathcal{O}\left(t^{2}\right)$. Kita juga boleh menggunakan PF peringkat lebih tinggi (formula Lie-Trotter-Suzuki), yang menumpu lebih pantas, dan ditakrifkan secara rekursif sebagai:

$S_2(t):=\prod_{a=1}^d e^{-i F_a t/2}\prod_{a=1}^d e^{-i F_a 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,$

di mana $\chi$ ialah peringkat PF simetri dan $s_p = \left( 4 - 4^{1/(2p-1)} \right)^{-1}$. Untuk evolusi masa yang panjang, kita boleh membahagikan selang masa $t$ kepada $k$ selang, dipanggil langkah Trotter, berdurasi $t/k$ dan menghampiri evolusi masa dalam setiap selang dengan formula produk peringkat $\chi$, iaitu $S_{\chi}$. Oleh itu, PF peringkat $\chi$ untuk pengoperasi evolusi masa merentasi $k$ langkah Trotter ialah:

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

di mana sebutan ralat berkurang dengan bilangan langkah Trotter $k$ dan peringkat $\chi$ bagi PF.

Diberi integer $k \geq 1$ dan formula produk $S_{\chi}(t)$, keadaan yang dievolusi masa secara hampiran $\rho_k(t)$ boleh diperoleh daripada $\rho_0$ dengan menerapkan $k$ lelaran formula produk $S_{\chi}\left(\frac{t}{k}\right)$.

$$\rho_k(t)=S_{\chi}\left(\frac{t}{k}\right)^k \rho_0 S_{\chi}\left(\frac{t}{k}\right)^{-k}$$

$\rho_k(t)$ ialah hampiran bagi $\rho(t)$ dengan ralat hampiran Trotter ||$\rho_k(t)-\rho(t) ||$. Jika kita mempertimbangkan gabungan linear bagi hampiran Trotter untuk $\rho(t)$:

$$\mu(t) = \sum_{j}^{l} x_j \rho^{k_j}_{j}\left(\frac{t}{k_j}\right) + \text{beberapa ralat Trotter yang tinggal} \, ,$$

di mana $x_j$ ialah pekali pemberat kita, $\rho^{k_j}_j$ ialah matriks ketumpatan yang sepadan dengan keadaan tulen yang diperoleh dengan mengevolusi keadaan awal menggunakan formula produk, $S^{k_j}_{\chi}$, yang melibatkan $k_j$ langkah Trotter, dan $j \in {1, ..., l}$ mengindeks bilangan PF yang membentuk MPF. Semua sebutan dalam $\mu(t)$ menggunakan formula produk yang sama $S_{\chi}(t)$ sebagai asasnya. Matlamatnya adalah untuk mengatasi ||$\rho_k(t)-\rho(t) \|$ dengan mencari $\mu(t)$ yang mempunyai $\|\mu(t)-\rho(t)\|$ yang lebih rendah lagi.

• $\mu(t)$ tidak semestinya keadaan fizikal kerana $x_i$ tidak semestinya positif. Matlamat di sini adalah untuk meminimumkan ralat dalam nilai jangkaan observable dan bukan untuk mencari pengganti fizikal bagi $\rho(t)$.
• $k_j$ menentukan kedua-dua kedalaman Circuit dan tahap hampiran Trotter. Nilai $k_j$ yang lebih kecil menghasilkan Circuit yang lebih pendek, yang menanggung lebih sedikit ralat Circuit tetapi akan menjadi hampiran yang kurang tepat bagi keadaan yang dikehendaki.

Perkara utama di sini ialah ralat Trotter yang tinggal yang diberikan oleh $\mu(t)$ adalah lebih kecil daripada ralat Trotter yang akan diperoleh dengan hanya menggunakan nilai $k_j$ terbesar.

Kamu boleh melihat kegunaan ini dari dua perspektif:

1. Untuk bajet langkah Trotter yang tetap yang boleh kamu jalankan, kamu boleh mendapatkan keputusan dengan ralat Trotter yang lebih kecil secara keseluruhan.
2. Diberi beberapa sasaran bilangan langkah Trotter yang terlalu besar untuk dijalankan, kamu boleh menggunakan MPF untuk mencari koleksi Circuit berkedalaman rendah untuk dijalankan yang menghasilkan ralat Trotter yang serupa.

Keperluan

Sebelum memulakan tutorial ini, pastikan kamu telah memasang perkara berikut:

• Qiskit SDK v1.0 atau lebih baru, dengan sokongan visualisasi
• Qiskit Runtime v0.22 atau lebih baru (pip install qiskit-ibm-runtime)
• MPF Qiskit addons (pip install qiskit_addon_mpf)
• Qiskit addons utils (pip install qiskit_addon_utils)
• Pustaka Quimb (pip install quimb)
• Pustaka Qiskit Quimb (pip install qiskit-quimb)
• Numpy v0.21 untuk keserasian merentasi pakej (pip install numpy==0.21)

Bahagian I. Contoh skala kecil

Teroka kestabilan MPF

Tidak ada sekatan yang jelas pada pilihan bilangan langkah Trotter $k_j$ yang membentuk keadaan MPF $\mu(t)$. Walau bagaimanapun, ini mesti dipilih dengan teliti untuk mengelakkan ketidakstabilan dalam nilai jangkaan yang dikira daripada $\mu(t)$. Peraturan umum yang baik adalah untuk menetapkan langkah Trotter terkecil $k_{\text{min}}$ supaya $t/k_{\text{min}} \lt 1$. Jika kamu ingin mengetahui lebih lanjut tentang ini dan cara memilih nilai $k_j$ lain kamu, rujuk panduan Cara memilih langkah Trotter untuk MPF.

Dalam contoh di bawah, kita meneroka kestabilan penyelesaian MPF dengan mengira nilai jangkaan magnetisasi untuk pelbagai masa menggunakan keadaan yang dievolusi masa yang berbeza. Secara khusus, kita membandingkan nilai jangkaan yang dikira daripada setiap evolusi masa hampiran yang dilaksanakan dengan langkah Trotter yang sepadan dan pelbagai model MPF (pekali statik dan dinamik) dengan nilai tepat bagi observable yang dievolusi masa. Pertama, mari kita takrifkan parameter untuk formula Trotter dan masa evolusi

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

Untuk contoh ini kita akan menggunakan keadaan Neel sebagai keadaan awal $\vert \text{Neel} \rangle = \vert 0101...01 \rangle$ dan model Heisenberg pada garis 10 tapak untuk Hamiltonian yang mengawal evolusi masa

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

di mana $J$ ialah kekuatan gandingan untuk tepi jiran terdekat.

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

Observable yang akan kita ukur ialah magnetisasi pada sepasang Qubit di tengah-tengah rantai.

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

Kita takrifkan laluan Transpiler untuk mengumpulkan putaran XX dan YY dalam Circuit sebagai satu Gate XX+YY tunggal. Ini akan membolehkan kita memanfaatkan sifat pemuliharaan spin TeNPy semasa pengiraan MPO, yang mempercepatkan pengiraan dengan ketara.

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

Kemudian kita buat Circuit yang melaksanakan evolusi masa Trotter secara hampiran.

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)

Seterusnya, kita kira nilai jangkaan yang dievolusi masa daripada Circuit Trotter.

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)

Kita juga kira nilai jangkaan tepat untuk perbandingan.

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)

Pekali MPF statik

MPF statik ialah yang nilai $x_j$-nya tidak bergantung pada masa evolusi, $t$. Mari kita pertimbangkan PF peringkat $\chi = 1$ dengan $k_j$ langkah Trotter, ini boleh ditulis sebagai:

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

di mana $A_n$ ialah matriks yang bergantung pada perolak $F_a$ dalam penguraian Hamiltonian. Adalah penting untuk diperhatikan bahawa $A_n$ sendiri adalah bebas daripada masa dan bilangan langkah Trotter $k_j$. Oleh itu, adalah mungkin untuk membatalkan sebutan ralat peringkat rendah yang menyumbang kepada $\mu(t)$ dengan pilihan berat $x_j$ bagi gabungan linear yang teliti. Untuk membatalkan ralat Trotter bagi $l-1$ sebutan pertama (ini akan memberikan sumbangan terbesar kerana ia sepadan dengan bilangan langkah Trotter yang lebih rendah) dalam ungkapan untuk $\mu(t)$, pekali $x_j$ mesti memenuhi persamaan berikut:

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

dengan $n=0, ... l-2$. Persamaan pertama menjamin tiada berat sebelah dalam keadaan yang dibina $\mu(t)$, manakala persamaan kedua memastikan pembatalan ralat Trotter. Untuk PF peringkat lebih tinggi, persamaan kedua menjadi $\sum_{j=1}^{l-1} \frac{x_j}{k_j^{\eta}} = 0$ di mana $\eta = \chi + 2n$ untuk PF simetri dan $\eta = \chi + n$ sebaliknya, dengan $n=0, ..., l-2$. Ralat yang terhasil (Ruj. [1],[2]) ialah kemudian

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

Menentukan pekali MPF statik untuk set nilai $k_j$ yang diberikan adalah sama dengan menyelesaikan sistem linear persamaan yang ditakrifkan oleh dua persamaan di atas untuk pemboleh ubah $x_j$: $Ax=b$. Di mana $x$ ialah pekali yang diminati, $A$ ialah matriks yang bergantung pada $k_j$ dan jenis PF yang kita gunakan ($S$), dan $b$ ialah vektor kekangan. Secara khusus:

$A_{0,j} = 1$ $A_{i>0,j} = k_{j}^{-(\chi + s(i-1))}$ $b_0 = 1$ $b_{i>0} = 0$

di mana $\chi$ ialah order, $s$ ialah $2$ jika symmetric adalah True dan $1$ sebaliknya, $k_{j}$ ialah trotter_steps, dan $x$ ialah pemboleh ubah untuk diselesaikan. Indeks $i$ dan $j$ bermula pada $0$. Kita juga boleh memvisualisasikan ini dalam bentuk matriks:

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

dan

$$b = \begin{bmatrix} b_{0} \\ b_{1} \\ b_{2} \\ ... \end{bmatrix} = \begin{bmatrix} 1 \\ 0 \\ 0 \\ ... \end{bmatrix}$$

Untuk maklumat lanjut, rujuk dokumentasi Sistem Linear Persamaan (LSE).

Kita boleh mencari penyelesaian untuk $x$ secara analitik sebagai $x = A^{-1}b$; lihat contohnya Ruj. [1] atau [2]. Walau bagaimanapun, penyelesaian tepat ini boleh "terkondisi buruk", menghasilkan norma-L1 yang sangat besar bagi pekali $x$ kita, yang boleh menyebabkan prestasi MPF yang buruk. Sebaliknya, seseorang boleh juga mendapatkan penyelesaian hampiran yang meminimumkan norma-L1 bagi $x$ dalam usaha mengoptimumkan kelakuan MPF.

Sediakan LSE

Kini setelah kita memilih nilai $k_j$ kita, kita mesti terlebih dahulu membina LSE, $Ax=b$ seperti yang diterangkan di atas. Matriks $A$ bergantung bukan sahaja pada $k_j$ tetapi juga pilihan PF kita, khususnya peringkat-nya. Selain itu, kamu mungkin mengambil kira sama ada PF simetri atau tidak (lihat [1]) dengan menetapkan symmetric=True/False. Walau bagaimanapun, ini tidak diperlukan, seperti yang ditunjukkan oleh Ruj. [2].

from qiskit_addon_mpf.static import setup_static_lse

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

Mari kita kerjakan nilai yang dipilih di atas untuk membina matriks $A$ dan vektor $b$. Dengan $j=0,1, 2$ langkah Trotter $k_j = [1, 2, 4]$, peringkat $\chi = 2$ dan pilihan langkah Trotter tidak simetri ($s=1$), kita dapati elemen matriks $A$ di bawah baris pertama ditentukan oleh ungkapan $A_{i>0,j} = k_{j}^{-(2 + 1(i-1))}$, khususnya:

$A_{0,0} = A_{0,1} = A_{0,2} = 1$ $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}$ $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}$

atau dalam bentuk matriks:

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

Ini boleh dilihat dengan memeriksa objek lse:

lse.A

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

Manakala vektor kekangan $b$ mempunyai elemen berikut: $b_{0} = 1$ $b_1 = b_2 = 0$

Oleh itu,

$$b = \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix}$$

Dan begitu juga dalam lse:

lse.b

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

Objek lse mempunyai kaedah untuk mencari pekali statik $x_j$ yang memenuhi sistem persamaan.

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]

Optimumkan $x$ menggunakan model tepat

Sebagai alternatif kepada pengiraan $x=A^{-1}b$, kamu juga boleh menggunakan setup_exact_model untuk membina contoh cvxpy.Problem yang menggunakan LSE sebagai kekangan dan penyelesaian optimalnya akan menghasilkan $x$.

Dalam bahagian seterusnya, akan jelas mengapa antara muka ini wujud.

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]

Sebagai petunjuk sama ada MPF yang dibina dengan pekali ini akan menghasilkan keputusan yang baik, kita boleh menggunakan norma-L1 (lihat juga Ruj. [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

Optimumkan $x$ menggunakan model hampiran

Mungkin berlaku bahawa norma-L1 untuk set nilai $k_j$ yang dipilih dianggap terlalu tinggi. Jika demikian dan kamu tidak boleh memilih set nilai $k_j$ yang berbeza, kamu boleh menggunakan penyelesaian hampiran kepada LSE dan bukannya penyelesaian tepat.

Untuk berbuat demikian, cukup gunakan setup_approximate_model untuk membina contoh cvxpy.Problem yang berbeza, yang mengehadkan norma-L1 kepada ambang yang dipilih sambil meminimumkan perbezaan $Ax$ dan $b$.

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

Perhatikan bahawa kamu mempunyai kebebasan penuh dalam cara menyelesaikan masalah pengoptimuman ini, yang bermakna kamu boleh menukar penyelesai pengoptimuman, ambang penumpuannya, dan sebagainya. Semak panduan masing-masing tentang Cara menggunakan model hampiran.

Pekali MPF dinamik

Dalam bahagian sebelumnya, kita memperkenalkan MPF statik yang mengatasi hampiran Trotter standard. Walau bagaimanapun, versi statik ini tidak semestinya meminimumkan ralat hampiran. Secara konkrit, MPF statik, dilambangkan $\mu^S(t)$, bukan unjuran optimum bagi $\rho(t)$ ke atas subruang yang direntang oleh keadaan formula produk $\{\rho_{k_i}(t)\}_{i=1}^r$.

Untuk menangani perkara ini, kita mempertimbangkan MPF dinamik (diperkenalkan dalam Ruj. [2] dan dibuktikan secara eksperimen dalam Ruj. [3]) yang meminimumkan ralat hampiran dalam norma Frobenius. Secara formal, kita fokus pada meminimumkan

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

berkenaan dengan beberapa pekali $x_i(t)$ pada setiap masa $t$. Unjuran optimum dalam norma Frobenius ialah kemudian $\mu^D(t) \;=\; \sum_{i=1}^r x_i(t)\,\rho_{k_i}(t)$, dan kita panggil $\mu^D(t)$ sebagai MPF dinamik. Dengan memasukkan definisi di atas:

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

di mana $M_{i,j}(t)$ ialah matriks Gram, ditakrifkan oleh

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

dan

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

mewakili pertindihan antara keadaan tepat $\rho(t)$ dan setiap hampiran formula produk $\rho_{k_i}(t)$. Dalam senario praktikal, pertindihan ini mungkin hanya boleh diukur secara hampiran disebabkan oleh hingar atau akses separa kepada $\rho(t)$.

Di sini, $\lvert\psi_{\mathrm{in}}\rangle$ ialah keadaan awal, dan $S(\cdot)$ ialah operasi yang diterapkan dalam formula produk. Dengan memilih pekali $x_i(t)$ yang meminimumkan ungkapan ini (dan mengendalikan data pertindihan hampiran apabila $\rho(t)$ tidak diketahui sepenuhnya), kita mendapatkan hampiran dinamik "terbaik" (dalam erti norma Frobenius) bagi $\rho(t)$ dalam subruang MPF. Kuantiti $L_i(t)$ dan $M_{i,j}(t)$ boleh dikira dengan cekap menggunakan kaedah rangkaian tensor [3]. MPF Qiskit addon menyediakan beberapa "Backend" untuk menjalankan pengiraan. Contoh di bawah menunjukkan cara yang paling fleksibel untuk berbuat demikian, dan dokumentasi Backend berasaskan lapisan TeNPy juga menerangkan secara terperinci. Untuk menggunakan kaedah ini, mulakan dari Circuit yang melaksanakan evolusi masa yang dikehendaki dan buat model yang mewakili operasi ini daripada lapisan Circuit yang sepadan. Akhir sekali, objek Evolver dibuat yang boleh digunakan untuk menghasilkan kuantiti yang dievolusi masa $M_{i,j}(t)$ dan $L_i(t)$. Kita mulakan dengan mencipta objek Evolver yang sepadan dengan evolusi masa hampiran (ApproxEvolverFactory) yang dilaksanakan oleh Circuit. Khususnya, perhatikan dengan teliti pemboleh ubah order supaya ia sepadan. Perhatikan bahawa dalam menjana Circuit yang sepadan dengan evolusi masa hampiran, kita menggunakan nilai pemegang tempat untuk time = 1.0 dan bilangan langkah Trotter (reps=1). Circuit hampiran yang betul kemudiannya dihasilkan oleh penyelesai masalah dinamik dalam 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,
    },
)

amaran

Pilihan LayerwiseEvolver yang menentukan butiran simulasi rangkaian tensor mesti dipilih dengan teliti untuk mengelakkan persediaan masalah pengoptimuman yang tidak terdefinisi.

Kemudian kita sediakan pengevolusi tepat (contohnya, ExactEvolverFactory), yang mengembalikan objek Evolver yang mengira evolusi masa sebenar atau "rujukan". Secara realistiknya, kita akan menghampiri evolusi tepat dengan menggunakan formula Suzuki-Trotter peringkat lebih tinggi atau kaedah boleh dipercayai lain dengan langkah masa yang kecil. Di bawah, kita menghampiri keadaan yang dievolusi masa secara tepat dengan formula Suzuki-Trotter peringkat keempat menggunakan langkah masa kecil dt=0.1, yang bermaksud bilangan langkah Trotter yang digunakan pada masa $t$ ialah $k=t/dt$. Kita juga menentukan beberapa pilihan pemotongan khusus TeNPy untuk mengehadkan dimensi bon maksimum rangkaian tensor yang mendasari, serta nilai singular minimum bagi bon rangkaian tensor yang dipecah. Parameter ini boleh mempengaruhi ketepatan nilai jangkaan yang dikira dengan pekali MPF dinamik, jadi adalah penting untuk meneroka pelbagai nilai bagi mencari keseimbangan optimum antara masa pengiraan dan ketepatan. Perhatikan bahawa pengiraan pekali MPF tidak bergantung pada nilai jangkaan PF yang diperoleh daripada pelaksanaan perkakasan, dan oleh itu ia boleh ditala dalam pasca-pemprosesan.

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

Seterusnya, buat keadaan awal sistem kamu dalam format yang serasi dengan TeNPy (contohnya, MPS_neel_state=$\vert 0101...01 \rangle$). Ini menyediakan fungsi gelombang berbilang zarah yang akan kamu evolusi dalam masa $\lvert\psi_{\mathrm{in}}\rangle$ sebagai tensor.

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)

Untuk setiap langkah masa $t$ kita sediakan sistem persamaan linear dinamik dengan kaedah setup_dynamic_lse. Objek yang sepadan mengandungi maklumat tentang masalah MPF dinamik: lse.A memberikan matriks Gram $M$ manakala lse.b memberikan pertindihan $L$. Kita kemudiannya boleh menyelesaikan LSE (apabila tidak terdefinisi buruk) untuk mencari pekali dinamik menggunakan setup_frobenius_problem. Adalah penting untuk diperhatikan perbezaan dengan pekali statik, yang hanya bergantung pada butiran formula produk yang digunakan dan bebas daripada butiran evolusi masa (Hamiltonian dan keadaan awal).

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

Akhir sekali, plot nilai jangkaan ini sepanjang masa evolusi.

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

Dalam kes seperti contoh di atas, di mana PF $k=1$ berkelakuan buruk pada semua masa, kualiti keputusan MPF dinamik juga sangat terjejas. Dalam situasi sedemikian, adalah berguna untuk menyiasat kemungkinan menggunakan PF individu dengan bilangan langkah Trotter yang lebih tinggi untuk meningkatkan kualiti keseluruhan keputusan. Dalam simulasi ini, kita melihat interaksi pelbagai jenis ralat: ralat daripada pensampelan terhingga, dan ralat Trotter daripada formula produk. MPF membantu mengurangkan ralat Trotter akibat formula produk tetapi menanggung ralat pensampelan yang lebih tinggi berbanding formula produk. Ini boleh menjadi kelebihan, kerana formula produk boleh mengurangkan ralat pensampelan dengan peningkatan pensampelan, tetapi ralat sistematik akibat hampiran Trotter tetap tidak berubah.

Langkah 1: Petakan input klasik ke masalah kuantum

Sekarang kita akan pertimbangkan satu masa $t=1.0$ dan kira nilai jangkaan magnetisasi menggunakan pelbagai kaedah dengan satu QPU. Masa $t$ ini dipilih khas supaya perbezaan antara pelbagai kaedah dapat dimaksimumkan dan kita boleh lihat keberkesanan relatif mereka. Untuk tentukan julat masa di mana MPF dinamik dijamin menghasilkan pemboleh ubah dengan ralat lebih rendah berbanding sebarang formula Trotter secara individu dalam multi-product, kita boleh laksanakan "ujian MPF" — rujuk persamaan (17) dan teks sekitarnya dalam [3].

Sediakan Circuit Trotter

Pada ketika ini, kita sudah jumpa pekali pengembangan, $x$, dan yang tinggal hanyalah menjana Circuit kuantum Trotterisasi. Sekali lagi, modul qiskit_addon_utils.problem_generators hadir membantu dengan fungsi yang berguna untuk melakukan ini:

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)

Langkah 2: Optimumkan masalah untuk pelaksanaan perkakasan kuantum

Mari kita kembali kepada pengiraan nilai jangkaan untuk satu titik masa. Kita akan pilih satu Backend untuk menjalankan eksperimen pada perkakasan sebenar.

from qiskit_ibm_runtime import QiskitRuntimeService

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

qubits = list(range(backend.num_qubits))

Kemudian kita buang Qubit terpencil dari peta gandingan supaya peringkat tataletak dalam Transpiler tidak memasukkan mereka. Di bawah kita gunakan sifat Backend yang dilaporkan yang tersimpan dalam objek target dan buang Qubit yang mempunyai ralat pengukuran atau Gate dua-Qubit melebihi ambang tertentu (max_meas_err, max_twoq_err) atau masa $T_2$ (yang menentukan kehilangan keherensi) di bawah ambang tertentu (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

Kita perlu tetapkan max_meas_err, min_t2, dan max_twoq_err supaya kita jumpa subset Qubit yang cukup besar untuk menyokong Circuit yang hendak dijalankan. Dalam kes kita, cukup sekadar mencari rantaian 1D 10 Qubit.

cust_cmap.draw()

Kemudian kita boleh petakan Circuit dan pemboleh ubah ke Qubit fizikal pada peranti.

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

Langkah 3: Jalankan menggunakan primitif Qiskit

Dengan primitif Estimator kita boleh dapatkan anggaran nilai jangkaan dari QPU. Kita jalankan Circuit AQC yang telah dioptimumkan bersama teknik pengurangan dan penindasan ralat tambahan.

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

Langkah 4: Proses pasca dan kembalikan keputusan dalam format klasik yang dikehendaki

Satu-satunya langkah pemprosesan pasca ialah menggabungkan nilai jangkaan yang diperoleh dari primitif Qiskit Runtime pada langkah Trotter yang berbeza menggunakan pekali MPF masing-masing. Untuk pemboleh ubah $A$ kita ada:

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

Pertama, kita ekstrak nilai jangkaan individu yang diperoleh untuk setiap Circuit Trotter:

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

Seterusnya, kita gabungkan semula dengan pekali MPF kita untuk menghasilkan jumlah nilai jangkaan MPF. Di bawah, kita buat ini untuk setiap cara berbeza yang kita gunakan untuk mengira $x$.

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

Akhirnya, untuk masalah kecil ini kita boleh kira nilai rujukan tepat menggunakan scipy.linalg.expm seperti berikut:

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

Dalam contoh di atas, kaedah MPF dinamik memberikan prestasi terbaik dari segi nilai jangkaan, mengatasi apa yang kita akan dapat sekiranya kita hanya menggunakan bilangan langkah Trotter tertinggi sahaja. Walaupun pelbagai teknik MPF tidak sentiasa mencapai nilai jangkaan yang lebih baik berbanding bilangan langkah Trotter tertinggi (seperti model tepat dan anggaran dalam plot di atas), sisihan piawai nilai-nilai ini menangkap dengan baik peningkatan varians yang berlaku apabila menggunakan teknik MPF. Ini menonjolkan ketidaktentuan sekitar nilai jangkaan yang diperoleh, yang sentiasa merangkumi nilai jangkaan yang kita jangkakan dari evolusi masa tepat sistem tersebut. Sebaliknya, nilai jangkaan yang dikira dengan bilangan langkah Trotter yang lebih rendah gagal menangkap nilai jangkaan tepat dalam ketidaktentuan mereka, justeru secara yakin memberikan keputusan yang salah.

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

Bahagian II: skala lebih besar

Jom kita tingkatkan skala masalah ni sampai melepasi apa yang boleh disimulasi dengan tepat. Dalam bahagian ni kita akan fokus untuk meniru semula sebahagian hasil yang ditunjukkan dalam Rujukan [3].

Langkah 1: Petakan input klasik kepada masalah kuantum

Hamiltonian

Untuk contoh berskala besar, kita guna model XXZ pada barisan 50 tapak:

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

di mana $J_{i,(i+1)}$ ialah pekali rawak yang sepadan dengan tepi $(i, i+1)$. Ini adalah Hamiltonian yang dipertimbangkan dalam demonstrasi yang ditunjukkan dalam Rujukan [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")

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

Untuk pemboleh ukur, kita pilih $Z_{24}Z_{25}$, seperti yang dapat dilihat dalam panel bawah Rajah 5 dalam Rujukan [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])

Pilih langkah Trotter

Eksperimen yang ditunjukkan dalam Rajah 4 Rujukan [3] menggunakan $k_j = [2, 3, 4]$ langkah Trotter simetri peringkat $2$. Kita fokus pada keputusan untuk masa $t=3$, di mana MPF dan PF dengan lebih banyak langkah Trotter (6 dalam kes ini) mempunyai ralat Trotter yang sama. Namun, nilai jangkaan MPF dikira daripada Circuit yang sepadan dengan bilangan langkah Trotter yang lebih rendah dan oleh itu lebih cetek. Dalam praktik, walaupun MPF dan Circuit langkah Trotter yang lebih dalam mempunyai ralat Trotter yang sama, kita jangkakan nilai jangkaan eksperimen yang dikira daripada Circuit MPF akan lebih hampir dengan nilai teori, kerana ia melibatkan pelaksanaan Circuit yang lebih cetek dan kurang terdedah kepada hingar perkakasan berbanding Circuit PF langkah Trotter yang lebih tinggi.

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

Sediakan LSE

Di sini kita tengok pekali MPF statik untuk masalah ini.

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

Pekali dinamik

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

Bina setiap Circuit Trotter dalam penguraian MPF kita

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)

Bina Circuit Trotter dengan ralat Trotter yang setanding dengan 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)

Langkah 2: Optimumkan masalah untuk pelaksanaan perkakasan kuantum

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
]

Langkah 3: Laksanakan menggunakan primitif 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)
    ]
)

Langkah 4: Proses pasca dan kembalikan hasil dalam format klasik yang dikehendaki

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

Bila melaksanakan Circuit pada perkakasan, kita mungkin menghadapi cabaran tambahan dalam mendapatkan nilai jangkaan yang tepat disebabkan kehadiran hingar perkakasan. Perkara ini tidak diambil kira dalam formalisme MPF dan boleh menjadi halangan kepada penyelesaian MPF. Sebagai contoh, inilah mungkin sebabnya pekali dinamik gagal memberikan anggaran nilai jangkaan yang lebih baik berbanding pekali statik anggaran dalam plot. Iaitu, pengevolusi anggaran yang mensimulasikan Circuit anggaran, tidak mencerminkan dengan tepat keputusan yang diperoleh dengan melaksanakan Circuit anggaran dalam kehadiran hingar perkakasan. Atas sebab-sebab ini, adalah disyorkan untuk menggabungkan pelbagai teknik pengurangan ralat bagi mendapatkan keputusan yang sedekat mungkin dengan nilai ideal bagi setiap formula hasil darab. Ini akan menunjukkan manfaat yang konsisten daripada pendekatan MPF.

Secara keseluruhan, pekali statik anggaran masih memberikan penyelesaian yang lebih tepat berbanding formula hasil darab dengan bilangan langkah Trotter yang lebih tinggi dengan jumlah ralat Trotter yang sama dalam tetapan tanpa hingar.

Penting juga untuk diambil perhatian bahawa dalam contoh yang meniru eksperimen dalam Rujukan [3], titik masa $t=3$ adalah melebihi had di mana PF dengan $k=2$ dijangka berfungsi dengan baik, iaitu $t/k>1$ seperti yang dibincangkan dalam panduan ini.

Rujukan

[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] Robertson, N. F., et al. (2024). Tensor network enhanced dynamic multiproduct formulas. arXiv preprint arXiv:2407.17405.