Synteza bramki Evolution w przestrzeni bozonowej

Ewolucji czasowej $e^{-itH}$ hamiltonianu będącego sumą nieprzemiennych składników nie można zaimplementować dokładnie bramka po bramce. Wzór produktowy przybliża ją ciągiem krótkich ewolucji pod działaniem poszczególnych składników, a błąd tego przybliżenia nazywa się błędem Trottera.

Wzór produktowy można zastosować w dwóch miejscach, które nie są wymienne:

1. W przestrzeni bozonowej, zanim wybrano jakiekolwiek obcięcie lub kodowanie. Wzór widzi grupy operatora (zobacz Grupowanie członów operatora: wykorzystanie struktury operatora) i może je przeplatać.

2. W przestrzeni kubitów, wewnątrz PauliEvolutionGate z Qiskit, po zmapowaniu operatora. Wzór widzi jeden operator Pauliego na grupę i może jedynie dzielić człony wewnątrz niej.

Tylko pierwszy sposób może zmniejszyć błąd między grupami, które nie komutują, a który jest zwykle dominujący. Ten przewodnik pokazuje, jak go użyć.

Ważne

Oba kroki są opcjonalne. Bramka Evolution, która nigdy nie jest dekomponowana w przestrzeni bozonowej, trafia w całości do etapu bozon–kubit i tam jest mapowana i syntetyzowana. Nie powoduje to żadnego błędu na poziomie bozonowym, a jedynie ten, który wprowadza wzór na poziomie kubitów.

Konfiguracja

Dwumiejscowy łańcuch Bosego-Hubbarda, najmniejszy przykład, w którym grupy nie komutują (pełny model znajdziesz w Symulacja dynamiki Bosego-Hubbarda z wyborem kodowania). Jest pogrupowany w jedną grupę hermitowską dla wiązania i jedną grupę diagonalną dla oddziaływania na węźle:

>>> import numpy as np
>>> from qiskit.quantum_info import SparsePauliOp, Statevector
>>> from qiskit.transpiler import PassManager
>>> from scipy.sparse import csc_matrix
>>> from scipy.sparse.linalg import expm_multiply
>>> from qiskit_bosons.circuit import BosonicCircuit
>>> from qiskit_bosons.circuit.library import Evolution, InitializeModes
>>> from qiskit_bosons.mappers import LevelSpec
>>> from qiskit_bosons.mappers.library import occupation_encoding
>>> from qiskit_bosons.operators import BosonOperator, ann, cre
>>> from qiskit_bosons.transpiler.presets import generate_preset_occupation_pass_manager
>>>
>>> num_sites, max_occupation, total_time = 2, 2, 0.8
>>> bonds = [(0, 1)]
>>> terms = {}
>>> for left, right in bonds:
...     terms[(cre(left), ann(right))] = -1.0
...     terms[(cre(right), ann(left))] = -1.0
>>> for site in range(num_sites):
...     terms[(cre(site), cre(site), ann(site), ann(site))] = 1.0
>>> hamiltonian = BosonOperator.from_dict(terms)
>>> hamiltonian.groups = [
...     len(bonds)
...     if len({mode for _, mode in term}) == 1
...     else bonds.index(tuple(sorted({mode for _, mode in term})))
...     for term, _ in hamiltonian.iter_terms()
... ]
>>> hamiltonian.num_groups()
2

Aby zmierzyć błąd, potrzebujesz wartości odniesienia. Właściwa jest dokładna ewolucja hamiltonianu zmapowanego przy tym samym obcięciu, ponieważ obcięcie jest osobnym przybliżeniem:

>>> spec = LevelSpec(max_occupation, "unary_projected")
>>> mapped = occupation_encoding(hamiltonian, spec).simplify()
>>> matrix = csc_matrix(SparsePauliOp.from_sparse_observable(mapped).to_matrix(sparse=True))
>>> pm = generate_preset_occupation_pass_manager(spec)
>>> pm.qubit = PassManager()  # keep each group as one exact PauliEvolutionGate
>>>
>>> def build(**evolution_kwargs):
...     circuit = BosonicCircuit(num_sites)
...     circuit.append(InitializeModes([1, 1]), circuit.modes)
...     circuit.append(
...         Evolution(num_sites, hamiltonian, time=total_time, **evolution_kwargs), circuit.modes
...     )
...     return circuit
>>>
>>> # Under the unary encodings, occupation ``n`` of mode ``m`` sets qubit ``m * d + n``, with
>>> # ``d = max_occupation + 1`` levels per mode. That places |1, 1> on the basis state below.
>>> levels = max_occupation + 1
>>> index = sum(1 << (mode * levels + n) for mode, n in enumerate([1, 1]))
>>> initial = Statevector.from_int(index, 2**(num_sites * levels)).data
>>> exact = expm_multiply(-1j * total_time * matrix, initial)

>>> def error(circuit):
...     """Distance to the exact evolution, of the circuit's product of group exponentials."""
...     lowered = pm.run(circuit)
...     state = Statevector.from_int(0, 2**lowered.num_qubits).evolve(lowered).data
...     return float(np.max(np.abs(state - exact)))

Uwaga

error celowo symuluje zredukowany obwód bez dekomponowania jego instancji PauliEvolutionGate. Statevector stosuje każdą z nich jako dokładną eksponentę jej grupy, więc pozostaje tylko iloczyn eksponent grup wybrany przez syntezę bozonową.

Ustawienie domyślne

Bez argumentu synthesis bramka Evolution jest dekomponowana za pomocą BosonicLieTrotter, czyli wzoru pierwszego rzędu: jeden czynnik na grupę, każdy na cały czas.

>>> dict(build().decompose().count_ops())
{'Evolution': 2, 'InitializeModes': 1}

Bez dekompozycji cała bramka trafia do etapu kubitowego i nie ma żadnego bozonowego błędu Trottera:

>>> error(build()) < 1e-6
True

Po dekompozycji pojawia się prawdziwy błąd, ponieważ obie grupy nie komutują:

>>> error(build().decompose()) > 0.1
True

Zmniejszanie błędu

Przekaż bramce inną metodę syntezy. Wyższy rząd składa czynniki symetrycznie, tak że wiodące człony błędu się znoszą; więcej powtórzeń dzieli czas na krótsze kroki.

>>> from qiskit_bosons.circuit.library.synthesis import BosonicLieTrotter, BosonicSuzukiTrotter
>>>
>>> for label, synthesis in [
...     ("Lie-Trotter", BosonicLieTrotter()),
...     ("Suzuki, order 2", BosonicSuzukiTrotter(order=2)),
...     ("Suzuki, order 4", BosonicSuzukiTrotter(order=4)),
...     ("Suzuki, order 4*", BosonicSuzukiTrotter(order=4, merge_adjacent=False)),
... ]:
...     circuit = build(synthesis=synthesis).decompose()
...     print(f"{label:16s} {circuit.count_ops()['Evolution']:3d} factors")
Lie-Trotter        2 factors
Suzuki, order 2    3 factors
Suzuki, order 4   11 factors
Suzuki, order 4*  15 factors

Każde zwiększenie rzędu o dwa składa pięć przejść rzędu niższego (wzrost $5^{(k-2)/2}$ z rekursji Suzukiego, zastosowany do $2n - 1$ czynników rzędu 2), co odpowiada liczbie w ostatnim wierszu, oznaczonej gwiazdką. Przejścia spotykają się na tym samym czynniku, więc sąsiednie bramki ewoluują pod działaniem jednego operatora i są scalane w jedną bramkę z zsumowanym czasem. Oszczędza to jedną bramkę na każde złączenie i jest dokładne. Przekaż merge_adjacent=False, aby zachować podręcznikową sekwencję. Tak czy inaczej wyższy rząd daje mniejszy błąd:

>>> lie = error(build(synthesis=BosonicLieTrotter()).decompose())
>>> second = error(build(synthesis=BosonicSuzukiTrotter(order=2)).decompose())
>>> fourth = error(build(synthesis=BosonicSuzukiTrotter(order=4)).decompose())
>>> lie > second > fourth
True

Powtórzenia wymieniają głębokość na dokładność w innym tempie. Błąd pierwszego rzędu maleje liniowo wraz z rozmiarem kroku, a drugiego rzędu kwadratowo:

>>> lie_errors = [error(build(synthesis=BosonicLieTrotter(reps=r)).decompose()) for r in (1, 2, 4, 8)]
>>> bool(all(a > b for a, b in zip(lie_errors, lie_errors[1:])))
True
>>> bool(lie_errors[0] / lie_errors[-1] > 6)  # about 8, i.e. linear in 1 / reps
True
>>> suzuki_errors = [
...     error(build(synthesis=BosonicSuzukiTrotter(order=2, reps=r)).decompose()) for r in (1, 2, 4)
... ]
>>> bool(suzuki_errors[0] / suzuki_errors[-1] > 10)  # about 16, i.e. quadratic in 1 / reps
True

To, który kompromis jest lepszy, zależy od operatora. Zauważ, że operator z jedną grupą jest dokładny w dowolnym rzędzie, ponieważ pojedynczy czynnik komutuje sam ze sobą, a wzór redukuje się wtedy do jednej bramki zamiast dzielić ją ponownie.

Pomiar kosztu

Mniejszy błąd kosztuje więcej czynników, a każdy czynnik kosztuje bramki splątujące, gdy etap kubitowy go zsyntetyzuje. Ile ich jest, zależy od kodowania, więc policz je w preset pipeline, przy ustalonej bazie i bez optymalizacji na poziomie kubitów (optimization_level=0), aby zobaczyć cenę każdej metody:

>>> basis = ["rz", "sx", "x", "cx"]
>>>
>>> def cx_count(circuit, encoding):
...     spec = LevelSpec(max_occupation, encoding)
...     lowered = generate_preset_occupation_pass_manager(
...         spec, optimization_level=0, basis_gates=basis
...     ).run(circuit)
...     return lowered.count_ops().get("cx", 0)
>>>
>>> methods = [
...     ("Lie-Trotter", BosonicLieTrotter()),
...     ("Suzuki, order 2", BosonicSuzukiTrotter(order=2)),
...     ("Suzuki, order 4", BosonicSuzukiTrotter(order=4)),
... ]
>>> encodings = ("gray", "unary_unprojected")
>>> counts = {
...     encoding: [
...         cx_count(build(synthesis=method).decompose(), encoding) for _, method in methods
...     ]
...     for encoding in encodings
... }
>>> # the order buys accuracy at a price that grows with it
>>> bool(counts["gray"][0] < counts["gray"][1] < counts["gray"][2])
True

Dokładne liczby zależą od tego, jak Qiskit syntetyzuje bramki ewolucji, co zmienia się między wydaniami, dlatego zamiast tabeli pokazuje je wykres:

Każda grupa słupków to metoda: Lie-Trotter, Suzuki rzędu 2 i Suzuki rzędu 4. Rząd kupuje dokładność za cenę rosnącą o czynnik rekursji, a kodowanie unarne płaci więcej niż kod Graya przy każdej metodzie, ponieważ każdy jego czynnik mapuje się na więcej bramek splątujących. Wybierz metodę ze względu na potrzebną dokładność, a kodowanie ze względu na sprzęt, i rozdziel oba wybory: powyższy obwód bozonowy nigdy się nie zmienia.

Zastosowanie jednej metody do całego obwodu

Aby ustawić synthesis na każdej bramce, musisz przeprowadzić ten wybór przez wszystko, co ją buduje. Przebieg B2BSynthesis dokonuje tego wyboru raz dla całego pipeline'u. Jego config wskazuje metodę dla każdego typu bramki, a przebieg zarówno wybiera metodę, jak i rozwija bramki:

>>> from qiskit_bosons.transpiler import BosonicPassManager
>>> from qiskit_bosons.transpiler.passes import B2BSynthesis
>>>
>>> trotterize = BosonicPassManager([B2BSynthesis({"Evolution": BosonicSuzukiTrotter(order=2)})])
>>> dict(trotterize.run(build()).count_ops())
{'Evolution': 3, 'InitializeModes': 1}

Metodę można też wskazać nazwą, wraz z argumentami jej __init__, zamiast ją budować: B2BSynthesis({"Evolution": ("SuzukiTrotter", (), {"order": 2})}). Dostępne nazwy to te z B2BSynthesisPluginManager, a metodę, która nigdzie nie jest zarejestrowana, można przekazać jako instancję, jak wyżej.

Z apply=False przebieg tylko wybiera metodę i zostawia rozwinięcie do momentu, w którym jest potrzebne, na przykład aby najpierw przejrzeć obwód. Bramka atomowa nigdy nie jest ruszana, nawet oznaczana:

>>> selected = BosonicPassManager(
...     [B2BSynthesis({"Evolution": ("SuzukiTrotter", (), {"order": 2})}, apply=False)]
... ).run(build())
>>> selected.count_ops()["Evolution"]  # selected, not yet expanded
1
>>> selected.decompose().count_ops()["Evolution"]  # expanded on demand
3

Nierozwinięta bramka trafia w całości do etapu bozon–kubit, gdzie jest mapowana bez sprawdzania jej synthesis, więc każda metoda daje wtedy identyczny wynik.

Opcjonalny filter ogranicza przebieg do bramek, które na tym zyskują. Ewolucje, których grupy wszystkie komutują, jak hamiltonian złożony wyłącznie z operatorów liczby cząstek, są syntetyzowane dokładnie w dowolnym rzędzie, więc wyższy rząd tylko zwiększyłby tam głębokość.

Ograniczenie przebiegu do bramek, które na tym zyskują

Czysto diagonalny operator dobrze to ilustruje, ponieważ operatory liczby cząstek komutują ze sobą. Jego grupy są dokładne w pierwszym rzędzie, a drugi rząd zużywa bramki splątujące na darmo:

>>> # one group per term, each Hermitian; the order is fixed because order two evolves the first
>>> # group twice, so which group comes first changes the cost
>>> diagonal = BosonOperator.from_terms_with_groups(
...     [([cre(0), ann(0), cre(1), ann(1)], 1.0, 0), ([cre(0), ann(0)], 0.5, 1)]
... )
>>>
>>> def diagonal_cost(synthesis):
...     circuit = BosonicCircuit(num_sites)
...     circuit.append(Evolution(num_sites, diagonal, time=total_time, synthesis=synthesis), circuit.modes)
...     return cx_count(circuit.decompose(), "gray")
>>>
>>> bool(diagonal_cost(BosonicSuzukiTrotter(order=2)) > diagonal_cost(BosonicLieTrotter()))
True

Oba obwody są dokładne, więc drugi zużywa więcej bramek splątujących bez korzyści. (Symetryczny wzór drugiego rzędu stosuje pierwszą grupę dwa razy, każdą w połowie czasu, a zatem koszt zależy od tego, która grupa jest pierwsza.) filter przyjmuje predykat na DAGOpNode i pozostawia odrzucone bramki nietknięte, więc jeden przebieg może podnieść rząd dla bramek, które na tym zyskują, a resztę zostawić w pierwszym rzędzie. Poniższy predykat akceptuje bramkę tylko wtedy, gdy jej operator przenosi bozony między trybami, czyli gdy grupy nie komutują:

>>> def is_diagonal(operator):
...     """Whether every term only counts bosons, rather than moving them between modes."""
...     return all(
...         sorted(mode for created, mode in actions if created)
...         == sorted(mode for created, mode in actions if not created)
...         for actions, _ in operator.iter_terms()
...     )
>>>
>>> def mixed_circuit():
...     circuit = BosonicCircuit(num_sites)
...     circuit.append(Evolution(num_sites, hamiltonian, time=total_time), circuit.modes)
...     circuit.append(Evolution(num_sites, diagonal, time=total_time), circuit.modes)
...     return circuit
>>>
>>> order_2 = BosonicSuzukiTrotter(order=2)
>>> everywhere = cx_count(
...     BosonicPassManager([B2BSynthesis({"Evolution": order_2})]).run(mixed_circuit()), "gray"
... )
>>> restricted = cx_count(
...     BosonicPassManager(
...         [
...             B2BSynthesis(
...                 {"Evolution": order_2},
...                 filter=lambda node: not is_diagonal(node.op.operator),
...             )
...         ]
...     ).run(mixed_circuit()),
...     "gray",
... )
>>> bool(everywhere > restricted)  # the filter saves entangling gates
True

Jeśli uruchomisz przebieg ponownie, nic się nie zmieni: czynniki, które emituje, są atomowe.

Czynniki atomowe

Bramka Evolution oznaczona jako atomic nie ma definicji, więc decompose() zostawia ją na miejscu, bez względu na to, ile razy zostanie zastosowane. Synteza oznacza w ten sposób emitowane przez siebie czynniki, podobnie jak BosonicQDrift dla próbkowanych członów.

To coś więcej niż porządek. Dalszy podział grupy hermitowskiej, człon po członie, daje czynniki, które z osobna nie są hermitowskie, a eksponenta operatora niehermitowskiego nie jest unitarna:

>>> circuit = build()
>>> once = circuit.decompose()
>>> dict(circuit.decompose(reps=5).count_ops()) == dict(once.count_ops())
True

Ostrzeżenie

Każdy czynnik musi być hermitowski sam w sobie. Operator hermitowski podzielony na grupy może dać grupy niehermitowskie, takie jak człon przeskoku oddzielony od swojego sprzężenia. Sprawdź za pomocą groups_are_hermitian(). Operator bez grup jest dzielony człon po członie, a etap kubitowy odrzuca wtedy czynniki z komunikatem Operator contains complex coefficients. Symetryzacja rzędu parzystego częściowo znosi błąd czynnika niehermitowskiego, więc błędnie pogrupowany operator może wyglądać lepiej w rzędzie 2 niż w rzędzie 1, a mimo to być błędny.

Powtarzanie krótkiego kroku

Jako alternatywę na poziomie obwodu zbuduj jeden krótki krok Trottera i powtórz go. repeat() spłaszcza powtórzenia do jednego obwodu i zostawia granice sąsiadujące, więc przebiegi mogą działać ponad nimi:

>>> step = BosonicCircuit(num_sites)
>>> step.append(Evolution(num_sites, hamiltonian, time=total_time / 4), step.modes)
>>> dict(step.repeat(4).count_ops())
{'Evolution': 4}

Przekaż insert_barriers=True, aby zachować powtórzenia oddzielone. Bariera jest mapowana na każdy kubit rejestrów, których dotyka, ponieważ kodowanie potrzebuje kilku kubitów na tryb, a tryb nie ma własnego obrazu na poziomie pojedynczego kubitu:

>>> [inst.operation.name for inst in step.repeat(3, insert_barriers=True)._inner.data]
['Evolution', 'barrier', 'Evolution', 'barrier', 'Evolution']

Uwaga

repeat nie ma związku z reps metody syntezy. Synteza dzieli pojedynczą ewolucję na więcej krótszych czynników. Powtórzenie obwodu zawierającego jedną ewolucję jedynie mnoży jej całkowity czas.

Połączenie ze wzorem na poziomie kubitów

Oba poziomy składają się ze sobą, a tam, gdzie oba przybliżają, o dokładności wyniku decyduje słabszy z nich. Podniesienie rzędu bozonowego przy pierwszorzędowym wzorze po stronie kubitów daje niewiele. Wzór po stronie kubitów opisuje MapperFnEvolutionSynthesis. Zachowuje on kolejność członów Pauliego wybraną przez mapper, więc group_wise() może emitować człony każdej grupy obok siebie, a simplify() może zamiast tego narzucić kanoniczną kolejność.

Dalsze kroki

• Transpilacja obwodów bozonowych wyjaśnia etapy następujące po syntezie bozonowej i miejsce, w którym znajduje się wśród nich B2BSynthesis.

• Grupowanie członów operatora: wykorzystanie struktury operatora wyjaśnia, jak przypisywać i sprawdzać grupy, które wzór produktowy przeplata.

• Symulacja dynamiki Bosego-Hubbarda z wyborem kodowania stosuje schemat powtórzeń w kompletnym przykładzie.