Pagbawas ng Trotter error ng Hamiltonian dynamics gamit ang multi-product formulas

Sa notebook na ito, matututuhan mo kung paano gamitin ang isang Multi-Product Formula (MPF) upang makamit ang mas mababang Trotter error sa ating observable kumpara sa kung ano ang nakukuha ng pinakamalalim na Trotter circuit na talagang isasagawa natin. Gagawin mo ito sa pamamagitan ng pagdaan sa mga hakbang ng isang Qiskit pattern:

• Hakbang 1: I-map sa quantum problem
  • I-initialize ang Hamiltonian ng ating problema
  • Gumamit ng MPF upang bumuo ng Trotterized time-evolution circuits
• Hakbang 2: I-optimize ang problema
  • Dito ay tina-transpile natin ang ating mga circuit para sa isang GenericBackendV2
• Hakbang 3: Isagawa ang mga eksperimento
  • Gumamit ng StatevectorEstimator para sa pagiging simple sa notebook na ito
• Hakbang 4: I-reconstruct ang mga resulta
  • Kalkulahin ang MPF expectation value

Hakbang 1: Map sa quantum problem

1a: Pagse-set up ng ating Hamiltonian

Ginagamit natin ang Ising model sa isang linya ng 10 site:

$$\hat{\mathcal{H}}_{\text{Ising}} = \sum_{i=1}^{9} J_{i,(i+1)} Z_i Z_{(i+1)} + \sum_{i=1}^{10} h_i X_i \, ,$$

kung saan ang $J$ ay ang coupling strength sa pagitan ng dalawang site at ang $h$ ay ang panlabas na magnetic field. Ang qiskit_addon_utils package ay nagbibigay ng ilang reusable na functionalities para sa iba't ibang layunin.

Ang qiskit_addon_utils.problem_generators module nito ay nagbibigay ng mga function para bumuo ng mga Heisenberg-like Hamiltonian sa ibinigay na connectivity graph. Ang graph na ito ay maaaring rustworkx.PyGraph o CouplingMap na ginagawang madaling gamitin sa mga Qiskit-centric workflow.

Sa mga sumusunod, gumagawa tayo ng simpleng linya ng 10 qubits gamit ang CouplingMap.from_line method.

# Added by doQumentation — required packages for this notebook
!pip install -q numpy qiskit qiskit-addon-mpf qiskit-addon-utils rustworkx scipy

from qiskit.transpiler import CouplingMap

# Generate some coupling map to use for this example
coupling_map = CouplingMap.from_line(10, bidirectional=False)

from rustworkx.visualization import graphviz_draw

graphviz_draw(coupling_map.graph, method="circo")

Susunod, bumubuo tayo ng SparsePauliOp sa ibinigay na connectivity gamit ang nais na mga constant.

from qiskit_addon_utils.problem_generators import generate_xyz_hamiltonian

# Get a qubit operator describing the Ising field model
hamiltonian = generate_xyz_hamiltonian(
    coupling_map,
    coupling_constants=(0.0, 0.0, 1.0),
    ext_magnetic_field=(0.4, 0.0, 0.0),
)
print(hamiltonian)

SparsePauliOp(['IIIIIIIZZI', 'IIIIIZZIII', 'IIIZZIIIII', 'IZZIIIIIII', 'IIIIIIIIZZ', 'IIIIIIZZII', 'IIIIZZIIII', 'IIZZIIIIII', 'ZZIIIIIIII', 'IIIIIIIIIX', 'IIIIIIIIXI', 'IIIIIIIXII', 'IIIIIIXIII', 'IIIIIXIIII', 'IIIIXIIIII', 'IIIXIIIIII', 'IIXIIIIIII', 'IXIIIIIIII', 'XIIIIIIIII'],
              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, 0.4+0.j, 0.4+0.j, 0.4+0.j, 0.4+0.j, 0.4+0.j, 0.4+0.j, 0.4+0.j,
 0.4+0.j, 0.4+0.j, 0.4+0.j])

Ang observable na susukatin natin ay ang kabuuang magnetization na maaari nating buuin nang simple tulad ng ipinapakita sa ibaba:

from qiskit.quantum_info import SparsePauliOp

L = coupling_map.size()
observable = SparsePauliOp.from_sparse_list([("Z", [i], 1 / L / 2) for i in range(L)], num_qubits=L)
print(observable)

SparsePauliOp(['IIIIIIIIIZ', 'IIIIIIIIZI', 'IIIIIIIZII', 'IIIIIIZIII', 'IIIIIZIIII', 'IIIIZIIIII', 'IIIZIIIIII', 'IIZIIIIIII', 'IZIIIIIIII', 'ZIIIIIIIII'],
              coeffs=[0.05+0.j, 0.05+0.j, 0.05+0.j, 0.05+0.j, 0.05+0.j, 0.05+0.j, 0.05+0.j,
 0.05+0.j, 0.05+0.j, 0.05+0.j])

1b: Multi-Product Formulas

Binabawasan ng mga MPF ang Trotter error ng Hamiltonian dynamics sa pamamagitan ng weighted combination ng ilang circuit execution.

Para gawing mas konkreto ito, tinutukoy natin ang isang MPF bilang:

$$\mu(t) = \sum_{j} x_j \rho^{k_j}_{j}\left(\frac{t}{k_j}\right) + \text{some remaining Trotter error} \, ,$$

kung saan ang $x_j$ ay ang ating weighting coefficients, ang $\rho^{k_j}_j$ ay ang density matrix na tumutugma sa pure state na nakuha sa pamamagitan ng pag-evolve ng inisyal na state gamit ang product formula, $S^{k_j}$, na may kasamang $k_j$ Trotter steps, at ang $j$ ay nag-iindex sa bilang ng mga product formula na bumubuo sa MPF.

Ang susi dito ay ang natitirang Trotter error ay mas maliit kaysa sa Trotter error na makukuha sa pamamagitan ng simpleng paggamit ng pinakamalaking $k_j$ value!

Maaari mong tingnan ang pagiging kapaki-pakinabang ng MPF mula sa dalawang perspektibo:

1. Para sa fixed na badyet ng mga Trotter step na maaari mong isagawa, makakakuha ka ng mga resulta na may mas maliit na Trotter error sa kabuuan.
2. Para sa isang bilang ng mga Trotter step na nagreresulta sa malalim na mga circuit, maaari mong gamitin ang MPF upang makahanap ng ilang mas maikling-depth na mga circuit na patatakbuhin na nagreresulta sa katulad na Trotter error.

Isang panimula sa mga static na MPF

Ang mga static na MPF ay yaong ang mga $x_j$ value ay HINDI nakadepende sa evolution time, $t$.

Ang pagtukoy ng static MPF coefficients para sa ibinigay na set ng mga $k_j$ value ay nangangahulugan ng paglutas ng linear system of equations: $Ax=b$, kung saan ang $x$ ay ang ating mga coefficient na interesado, ang $A$ ay isang matrix na nakadepende sa $k_j$ at sa uri ng PF na ginagamit natin ($S$), at ang $b$ ay isang vector ng mga constraint. Para sa kaiklian, hindi tayo magpapakadetalyado pa dito at sa halip ay isasangguni ka natin sa documentation ng LSE.

Maaari nating mahanap ang solusyon para sa $x$ nang analytical bilang $x = A^{-1}b$, tingnan hal. Carrera Vazquez et al., 2023 o Zhuk et al., 2023. Gayunpaman, ang exact na solusyon na ito ay maaaring "ill-conditioned" na nagreresulta sa napakalaking L1-norms ng ating mga coefficient, $x$, na maaaring humantong sa masamang performance ng MPF. Sa halip, maaari ring makakuha ng approximate na solusyon na nag-mi-minimize ng L1-norm ng $x$ upang subukang i-optimize ang gawi ng MPF.

Sa mga sumusunod, matututuhan mo kung paano gawin ang lahat ng ito.

Pagpili ng $k_j$

Ang pagpili ng $k_j$ ay nasa end-user. Sa prinsipyo, anumang halaga ang maaaring piliin ngunit ang ilang $k_j$ ay hahantong sa mas malaking noise amplification sa mga totoong device kaysa sa ibang mga pagpili. Kaya, mahalagang subukan ng isa na maghanap ng "magagandang" value ng $k_j$.

Dito, kukunin lang natin ang ilang fixed na value para sa $k_j$. Ang pinakamaliit na value ay binigyang-katwiran ng target evolution time na $t=8.0$ na karaniwang nagsasabi sa atin na masiyahan sa $t/k_{\text{min}} \lt 1$ ngunit empirikal na alam natin na ang pagtatakda nito sa katumbas ng $1$ ay karaniwang gumagana rin. Kung gusto mong matuto pa tungkol dito at kung paano pumili ng iba pang $k_j$ value, sumangguni sa kaukulang guide: How to choose the Trotter steps for an MPF.

time = 8.0
trotter_steps = (8, 12, 19)

Pagse-set up ng LSE

Ngayon na napili na natin ang ating mga $k_j$, dapat muna nating buuin ang LSE, $Ax=b$ tulad ng ipinaliwanag sa itaas. Ang matrix na $A$ ay nakadepende hindi lamang sa $k_j$ kundi pati na rin sa ating pagpili ng product formula (PF) -- partikular ang order nito. Bukod pa rito, maaaring isaalang-alang ng isa kung ang PF ay symmetric o hindi (tingnan ang Carrera Vazquez et al., 2023), sa pamamagitan ng pagtatakda ng symmetric=True. Gayunpaman, hindi ito kinakailangan tulad ng ipinakita ng Zhuk et al., 2023.

Dito, gagamit tayo ng second-order Suzuki-Trotter formula na nag-yi-yield ng order=2 at itatakda natin ang symmetric=True.

from qiskit_addon_mpf.static import setup_static_lse

lse = setup_static_lse(trotter_steps, order=2, symmetric=True)
print(lse)

LSE(A=array([[1.00000000e+00, 1.00000000e+00, 1.00000000e+00],
       [1.56250000e-02, 6.94444444e-03, 2.77008310e-03],
       [2.44140625e-04, 4.82253086e-05, 7.67336039e-06]]), b=array([1., 0., 0.]))

Paglutas para sa $x$ nang analytical

Tulad ng nabanggit dati, maaari nating mahanap ang $x$ nang analytical:

import numpy as np

coeffs_analytical = lse.solve()
print(coeffs_analytical)

[ 0.17239057 -1.19447005  2.02207947]

Pag-optimize para sa $x$ gamit ang exact model

Bilang alternatibo sa pagkalkula ng $x=A^{-1}b$, maaari mo ring gamitin ang setup_exact_problem upang bumuo ng cvxpy.Problem instance na gumagamit ng LSE bilang mga constraint at ang optimal na solusyon nito ay magbibigay ng $x$.

Sa susunod na seksyon, magiging malinaw kung bakit umiiral ang interface na ito.

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.17239057 -1.19447005  2.02207947]

Bilang indikasyon kung ang isang MPF na binuo sa mga coefficient na ito ay magbibigay ng magagandang resulta, maaari nating gamitin ang L1-norm (tingnan din ang Carrera Vazquez et al., 2023).

print(np.linalg.norm(coeffs_exact.value, ord=1))

3.3889400921655914

Pag-optimize para sa $x$ gamit ang approximate model

Maaaring mangyari na ang L1 norm para sa napiling set ng mga $k_j$ value ay tinuturing na masyadong mataas. Kung ganoon ang kaso at hindi ka makakapili ng ibang set ng mga $k_j$ value, maaari kang gumamit ng approximate na solusyon sa LSE sa halip na exact.

Para gawin ito, gamitin lang ang setup_sum_of_squares_problem upang bumuo ng ibang cvxpy.Problem instance na nagpapaliit sa L1-norm sa napiling threshold habang minimi-minimize ang pagkakaiba ng $Ax$ at $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=3.0)
model_approx.solve()
print(coeffs_approx.value)
print(np.linalg.norm(coeffs_approx.value, ord=1))

[-0.40454257  0.57553173  0.8290123 ]
1.8090865903790838

Tandaan, na mayroon kang ganap na kalayaan kung paano malulutas ang optimization problem na ito, ibig sabihin maaari mong baguhin ang optimization solver, ang convergence thresholds nito, at iba pa. Tingnan ang kaukulang guide sa How to use the approximate model.

1c: Pagse-set up ng Trotter circuits

Sa puntong ito, nahanap na natin ang ating expansion coefficients, $x$, at ang natitira na lang gawin ay buuin ang Trotterized quantum circuits. Muli, ang qiskit_addon_utils.problem_generators module ay tumulong para gawin ito:

from qiskit.synthesis import SuzukiTrotter
from qiskit_addon_utils.problem_generators import generate_time_evolution_circuit

circuits = []
for k in trotter_steps:
    circ = generate_time_evolution_circuit(
        hamiltonian,
        synthesis=SuzukiTrotter(order=2, reps=k),
        time=time,
    )
    circuits.append(circ)

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

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

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

Hakbang 2: I-optimize ang problema

Karaniwan, ito ang hakbang sa pattern kung saan ino-optimize mo ang iyong mga circuit para sa execution sa hardware. Dito, dahil gumagamit lang tayo ng noiseless simulator, simpleng tina-transpile lang natin ang ating circuit para sa isang GenericBackendV2.

from qiskit.providers.fake_provider import GenericBackendV2
from qiskit.transpiler import generate_preset_pass_manager

backend = GenericBackendV2(num_qubits=10)
transpiler = generate_preset_pass_manager(optimization_level=2, backend=backend)

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

Hakbang 3: Isagawa ang quantum experiments

Tulad ng ipinaliwanag sa simula pa lang, ila-laktawan natin ang optimization step 2 dahil simpleng kakalkulahin lang natin ang mga expectation value ng ating target observable gamit ang noise-free simulator, partikular ang StatevectorEstimator.

from qiskit.primitives import StatevectorEstimator

estimator = StatevectorEstimator()
job = estimator.run([(circ, observable) for circ in transpiled_circuits])
result = job.result()

Hakbang 4: I-reconstruct ang mga resulta

Una, kinukuha natin ang mga indibidwal na expectation value na nakuha para sa bawat isa sa mga Trotter circuit:

evs = [res.data.evs for res in result]
print(evs)

[array(0.23799162), array(0.35754312), array(0.38649906)]

Susunod, simpleng pinagsasama-sama lang natin ang mga ito sa ating mga MPF coefficient upang makuha ang kabuuang mga expectation value ng MPF. Sa ibaba, ginagawa natin ito para sa bawat isa sa iba't ibang paraan kung paano natin nakuha ang $x$.

print("Analytical    solution:", evs @ coeffs_analytical)
print("Exact model   solution:", evs @ coeffs_exact.value)
print("Approx. model solution:", evs @ coeffs_approx.value)

Analytical    solution: 0.3954847855980006
Exact model   solution: 0.39548478559800204
Approx. model solution: 0.42991214253489807

Sa wakas, para sa maliit na problemang ito ay maaari nating kalkulahin ang exact reference value gamit ang scipy.linalg.expm tulad ng sumusunod:

from scipy.linalg import expm

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

initial_state = np.zeros(exp_H.shape[0])
initial_state[0] = 1.0

time_evolved_state = exp_H @ initial_state

exact_obs = time_evolved_state.conj() @ observable.to_matrix() @ time_evolved_state
print(exact_obs.real)

0.40060242487899755

Malinaw nating makikita na ang MPF ay nakabawas ng Trotter error kumpara sa nakuha sa pinakamalalim na indibidwal na PF na may $k_j=19$. Gayunpaman, nakikita rin natin na ang approximate model ay hindi perpekto dahil aktwal itong nagresulta sa expectation value na mas masama kaysa sa exact na solusyon. Ipinapakita nito ang kahalagahan ng paggamit ng masisigasig na convergence criteria sa approximate model tulad ng matututuhan mo sa guide na How to use the approximate model.