Syntéza hradla Evolution v bosonovém prostoru

Časový vývoj $e^{-itH}$ hamiltoniánu, který je součtem nekomutujících částí, nelze implementovat přesně hradlo po hradle. Součinový vzorec (product formula) ho aproximuje posloupností krátkých vývojů pod jednotlivými částmi a chyba aproximace se nazývá Trotterova chyba.

Součinový vzorec můžeš použít na dvou místech a ta nejsou zaměnitelná:

1. V bosonovém prostoru, ještě před zvolením jakéhokoli omezení (truncation) nebo kódování. Vzorec vidí skupiny operátoru (viz Seskupování členů operátoru: využij strukturu operátoru) a může je prokládat.

2. V qubitovém prostoru, uvnitř Qiskit PauliEvolutionGate, poté co byl operátor zmapován. Vzorec vidí jeden Pauliho operátor na skupinu a může dělit jen členy uvnitř ní.

Pouze první z nich umí snížit chybu mezi skupinami, které nekomutují, a ta bývá obvykle dominantní. Tento návod ukazuje, jak ji použít.

Důležité

Oba kroky jsou volitelné. Hradlo Evolution, které se v bosonovém prostoru nikdy nerozloží, dorazí do fáze boson-na-qubit celé a je tam zmapováno a syntetizováno. To nezpůsobí žádnou chybu na bosonové úrovni, jen tu, kterou zavede vzorec na úrovni qubitů.

Příprava

Dvoumístný Bose-Hubbardův řetězec, nejmenší příklad, ve kterém skupiny nekomutují (celý model viz Simulace Bose-Hubbardovy dynamiky s volbou kódování). Je seskupen do jedné hermitovské skupiny pro vazbu a jedné diagonální skupiny pro interakci na místě:

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

K měření chyby potřebuješ referenci. Správnou referencí je přesný vývoj hamiltoniánu zmapovaného při stejném omezení (truncation), protože omezení je samostatná aproximace:

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

Poznámka

error záměrně simuluje převedený obvod bez rozkladu jeho instancí PauliEvolutionGate. Statevector aplikuje každou z nich jako přesnou exponenciálu její skupiny, takže zbývá jen součin exponenciál skupin, který zvolila bosonová syntéza.

Výchozí chování

Bez argumentu synthesis se hradlo Evolution rozkládá pomocí BosonicLieTrotter, vzorce prvního řádu: jeden faktor na skupinu, každý po celou dobu.

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

Bez rozkladu je celé hradlo předáno qubitové fázi a neexistuje žádná bosonová Trotterova chyba:

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

Po rozkladu se objeví skutečná chyba, protože obě skupiny nekomutují:

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

Snížení chyby

Předej hradlu jinou metodu syntézy. Vyšší řád skládá faktory symetricky, takže se vedoucí členy chyby vyruší; více opakování dělí čas na kratší kroky.

>>> 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ždé zvýšení řádu o dvě skládá pět průchodů řádu nižšího (růst $5^{(k-2)/2}$ ze Suzukiho rekurze, aplikovaný na $2n - 1$ faktorů řádu 2), což je počet v posledním řádku označený hvězdičkou. Průchody se setkávají na stejném faktoru, takže sousední hradla se vyvíjejí pod jedním operátorem a slučují se do jediného hradla se sečteným časem. To ušetří jedno hradlo na každém spoji a je to přesné. Předej merge_adjacent=False, chceš-li zachovat učebnicovou posloupnost. Tak či onak vyšší řád přináší menší chybu:

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

Opakování vyměňují hloubku za přesnost jinou rychlostí. Chyba prvního řádu klesá lineárně s velikostí kroku, chyba druhého řádu kvadraticky:

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

Které řešení je výhodnější, závisí na operátoru. Všimni si, že operátor s jedinou skupinou je přesný v libovolném řádu, protože osamělý faktor komutuje sám se sebou, a vzorec se pak zredukuje na jedno hradlo místo opětovného dělení.

Změř cenu

Menší chyba stojí více faktorů a každý faktor stojí provázací hradla, jakmile ho qubitová fáze syntetizuje. Kolik jich je, závisí na kódování, takže je spočítej přes předvolený pipeline s pevnou bází a bez optimalizace na úrovni qubitů (optimization_level=0), abys viděl cenu každé 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

Přesné počty závisejí na tom, jak Qiskit syntetizuje hradla vývoje, což se mezi verzemi mění, proto je ukazuje graf místo tabulky:

Každá skupina sloupců je jedna metoda: Lie-Trotter, Suzuki řádu 2 a Suzuki řádu 4. Řád kupuje přesnost za cenu, která roste faktorem rekurze, a unární kódování platí víc než Grayův kód u každé metody, protože každý jeho faktor se mapuje na více provázacích hradel. Metodu vol podle potřebné přesnosti a kódování podle hardwaru a obě volby drž odděleně: bosonový obvod výše se nikdy nemění.

Použití jedné metody na celý obvod

Pokud chceš nastavit synthesis u každého hradla, musíš volbu protáhnout vším, co nějaké hradlo vytváří. Průchod B2BSynthesis místo toho provede volbu jednou pro celý pipeline. Jeho config pojmenovává metodu pro každý typ hradla a zároveň metodu vybírá i hradla rozkládá:

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

Metodu lze také pojmenovat spolu s argumenty pro její __init__, místo aby se vytvářela: B2BSynthesis({"Evolution": ("SuzukiTrotter", (), {"order": 2})}). Dostupná jména jsou ta z B2BSynthesisPluginManager a metodu, která není nikde zaregistrována, lze předat jako instanci, jako výše.

S apply=False průchod pouze vybere metodu a rozklad nechá na místo, kde je žádoucí, například abys nejdřív obvod prohlédl. Atomické hradlo se nikdy nedotkne, ani se neoznačí:

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

Nerozložené hradlo dorazí do fáze boson-na-qubit celé, kde se zmapuje bez ohledu na jeho synthesis, takže všechny metody pak dávají shodný výstup.

Volitelný filter ho omezí na hradla, kterým to prospěje. Vývoje, jejichž skupiny všechny komutují, například hamiltonián tvořený jen operátory počtu, se syntetizují přesně v libovolném řádu, takže vyšší řád by tam jen přidal hloubku.

Omezení průchodu na hradla, kterým to prospěje

Čistě diagonální operátor to ukáže názorně, protože operátory počtu spolu komutují. Jeho skupiny jsou přesné v prvním řádu a druhý řád zbytečně spotřebuje provázací hradla:

>>> # 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 obvody jsou přesné, takže druhý spotřebuje více provázacích hradel bez zisku. (Symetrický vzorec druhého řádu aplikuje první skupinu dvakrát v polovičním čase, a cena proto závisí na tom, která skupina je první.) filter přijímá predikát nad DAGOpNode a zamítnutá hradla nechává nedotčená, takže jeden průchod může zvýšit řád u hradel, kterým to prospěje, a zbytek nechat v prvním řádu. Níže uvedený predikát přijme hradlo jen tehdy, když jeho operátor přesouvá bosony mezi módy, tedy když skupiny nekomutují:

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

Spustíš-li průchod znovu, nic se nezmění: faktory, které vydává, jsou atomické.

Atomické faktory

Evolution označené jako atomic nemá žádnou definici, takže decompose() ho nechá na místě, kolikrát se aplikuje. Syntéza takto označuje faktory, které vydává, a stejně tak BosonicQDrift členy, které vzorkuje.

Není to jen úhlednost. Další dělení hermitovské skupiny člen po členu dává faktory, které jednotlivě hermitovské nejsou, a exponenciála nehermitovského operátoru není unitární:

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

Upozornění

Každý faktor musí být hermitovský sám o sobě. Hermitovský operátor rozdělený do skupin může dát nehermitovské skupiny, například hopping člen oddělený od svého adjungovaného. Zkontroluj pomocí groups_are_hermitian(). Operátor bez skupin se dělí člen po členu a qubitová fáze pak faktory odmítne s chybou Operator contains complex coefficients. Symetrizace sudého řádu částečně ruší chybu nehermitovského faktoru, takže špatně seskupený operátor může v řádu 2 vypadat lépe než v řádu 1, a přesto být chybný.

Opakování krátkého kroku

Jako alternativu na úrovni obvodu sestav jeden krátký Trotterův krok a opakuj ho. repeat() zploští opakování do jednoho obvodu a ponechá hranice sousedící, takže průchody mohou působit přes ně:

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

Předej insert_barriers=True, chceš-li opakování udržet oddělená. Bariéra se mapuje na každý qubit registrů, kterých se dotýká, protože kódování potřebuje několik qubitů na mód a mód nemá vlastní obraz po jednotlivých qubitech:

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

Poznámka

repeat nesouvisí s reps metody syntézy. Syntéza rozděluje jediný vývoj na více kratších faktorů. Opakování obvodu, který obsahuje jeden vývoj, jen násobí jeho celkový čas.

Kombinace s vzorcem na úrovni qubitů

Obě úrovně se skládají a tam, kde obě aproximují, určuje přesnost výsledku slabší z nich. Zvyšování bosonového řádu při vzorci na straně qubitů prvního řádu přináší málo. Vzorec na straně qubitů viz MapperFnEvolutionSynthesis. Zachovává pořadí Pauliho členů, které zvolil mapper, a tak group_wise() může vydávat členy každé skupiny vedle sebe a simplify() může místo toho pevně stanovit kanonické pořadí.

Další kroky

• Transpilace bosonových obvodů vysvětluje fáze následující po bosonové syntéze a kam mezi ně patří B2BSynthesis.

• Seskupování členů operátoru: využij strukturu operátoru vysvětluje, jak přiřadit a zkontrolovat skupiny, které součinový vzorec prokládá.

• Simulace Bose-Hubbardovy dynamiky s volbou kódování aplikuje schéma opakování na úplný příklad.