Ein Evolution-Gate im bosonischen Raum synthetisieren

Die Zeitentwicklung $e^{-itH}$ eines Hamiltonians, der eine Summe nicht kommutierender Teile ist, lässt sich nicht exakt Gate für Gate umsetzen. Eine Produktformel nähert sie durch eine Folge kurzer Entwicklungen unter den einzelnen Teilen an, und der Näherungsfehler heißt Trotter-Fehler.

Du kannst eine Produktformel an zwei Stellen anwenden, und die beiden sind nicht austauschbar:

1. Im bosonischen Raum, bevor Trunkierung oder Kodierung gewählt wurden. Die Formel sieht die Gruppen des Operators (siehe Operatorterme gruppieren: die Operatorstruktur nutzen) und kann sie ineinander verschachteln.

2. Im Qubit-Raum, innerhalb von Qiskits PauliEvolutionGate, nachdem der Operator abgebildet wurde. Die Formel sieht pro Gruppe einen Pauli-Operator und kann nur die Terme innerhalb davon weiter unterteilen.

Nur der erste Weg kann den Fehler zwischen nicht kommutierenden Gruppen verringern, der meist der dominante ist. Diese Anleitung zeigt, wie du ihn nutzt.

Wichtig

Beide Schritte sind optional. Ein Evolution-Gate, das im bosonischen Raum nie zerlegt wird, erreicht die Stufe Boson-zu-Qubit als Ganzes und wird dort abgebildet und synthetisiert. Das verursacht keinen Fehler auf bosonischer Ebene, sondern nur den, den die Formel auf Qubit-Ebene einführt.

Einrichtung

Eine Bose-Hubbard-Kette mit zwei Plätzen, das kleinste Beispiel, in dem die Gruppen nicht kommutieren (das vollständige Modell findest du unter Bose-Hubbard-Dynamik mit wählbarer Kodierung simulieren). Sie wird in eine hermitesche Gruppe für die Bindung und eine diagonale Gruppe für die Wechselwirkung am Platz gruppiert:

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

Um einen Fehler zu messen, brauchst du eine Referenz. Die exakte Entwicklung des Hamiltonians so, wie er bei derselben Trunkierung abgebildet wurde, ist die richtige, denn die Trunkierung ist eine eigene Näherung:

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

Hinweis

error simuliert den abgesenkten Circuit bewusst ohne seine PauliEvolutionGate-Instanzen zu zerlegen. Statevector wendet jede als exakte Exponentialfunktion ihrer Gruppe an, sodass nur noch das Produkt der Gruppenexponentiale übrig bleibt, das die bosonische Synthese gewählt hat.

Der Standard

Ohne das Argument synthesis wird ein Evolution-Gate mit BosonicLieTrotter zerlegt, der Formel erster Ordnung: ein Faktor pro Gruppe, jeweils für die gesamte Zeit.

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

Ohne Zerlegung wird das gesamte Gate an die Qubit-Stufe übergeben, und es gibt überhaupt keinen bosonischen Trotter-Fehler:

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

Zerlegt, tritt ein echter Fehler auf, weil die beiden Gruppen nicht kommutieren:

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

Den Fehler verringern

Übergib dem Gate eine andere Synthesemethode. Eine höhere Ordnung setzt die Faktoren symmetrisch zusammen, sodass sich die führenden Fehlerterme aufheben; mehr Wiederholungen teilen die Zeit in kürzere Schritte.

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

Jede Erhöhung der Ordnung um zwei setzt fünf Durchläufe der darunterliegenden Ordnung zusammen (das Wachstum $5^{(k-2)/2}$ der Suzuki-Rekursion, angewandt auf die $2n - 1$ Faktoren der Ordnung 2), was die Anzahl in der letzten, mit einem Stern markierten Zeile ist. Die Durchläufe treffen sich am selben Faktor, sodass die benachbarten Gates unter einem Operator entwickeln und zu einem einzigen Gate mit der summierten Zeit verschmolzen werden. Das spart ein Gate pro Naht und ist exakt. Übergib merge_adjacent=False, um die Lehrbuchfolge beizubehalten. In beiden Fällen erkauft eine höhere Ordnung einen kleineren Fehler:

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

Wiederholungen tauschen Tiefe gegen Genauigkeit zu einer anderen Rate. Der Fehler erster Ordnung fällt linear mit der Schrittweite, der zweiter Ordnung quadratisch:

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

Welcher Tausch besser ist, hängt vom Operator ab. Beachte, dass ein Operator mit nur einer Gruppe bei jeder Ordnung exakt ist, weil ein einzelner Faktor mit sich selbst kommutiert, und die Formel dann zu einem Gate zusammenfällt, statt es erneut aufzuteilen.

Die Kosten messen

Ein kleinerer Fehler kostet mehr Faktoren, und jeder Faktor kostet verschränkende Gates, sobald die Qubit-Stufe ihn synthetisiert hat. Wie viele, hängt von der Kodierung ab. Zähle sie daher über die Preset-Pipeline, mit fester Basis und ohne Optimierung auf Qubit-Ebene (optimization_level=0), um den Preis jeder Methode zu sehen:

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

Die genauen Zahlen hängen davon ab, wie Qiskit die Evolutions-Gates synthetisiert, was sich zwischen Releases ändert, deshalb zeigt die Grafik sie statt einer Tabelle:

Jede Balkengruppe ist eine Methode: Lie-Trotter, Suzuki der Ordnung 2 und Suzuki der Ordnung 4. Die Ordnung erkauft Genauigkeit zu einem Preis, der um den Faktor der Rekursion wächst, und die unäre Kodierung zahlt bei jeder Methode mehr als der Gray-Code, weil jeder ihrer Faktoren auf mehr verschränkende Gates abgebildet wird. Wähle die Methode nach der Genauigkeit, die du brauchst, und die Kodierung nach der Hardware, und halte die beiden Entscheidungen getrennt: Der bosonische Circuit oben ändert sich nie.

Eine Methode auf einen ganzen Circuit anwenden

Um synthesis an jedem Gate zu setzen, musst du die Wahl durch alles schleusen, was eines baut. Der B2BSynthesis-Pass trifft die Wahl stattdessen einmal für eine Pipeline. Sein config benennt die Methode für jeden Gate-Typ, und er wählt die Methode aus und expandiert die Gates:

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

Eine Methode kann auch benannt werden, mit den Argumenten für ihr __init__, statt gebaut zu werden: B2BSynthesis({"Evolution": ("SuzukiTrotter", (), {"order": 2})}). Die verfügbaren Namen sind die des B2BSynthesisPluginManager, und eine Methode, die nirgends registriert ist, kann wie oben als Instanz übergeben werden.

Mit apply=False wählt der Pass nur die Methode aus und überlässt die Expansion dem Zeitpunkt, an dem sie gewünscht ist, zum Beispiel um den Circuit zuerst zu inspizieren. Ein atomares Gate wird nie angefasst, nicht einmal markiert:

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

Ein nicht expandiertes Gate erreicht die Stufe Boson-zu-Qubit als Ganzes, wo es abgebildet wird, ohne dass sein synthesis herangezogen wird, sodass jede Methode dann identische Ausgabe erzeugt.

Ein optionaler filter beschränkt ihn auf die Gates, die profitieren. Evolutionen, deren Gruppen alle kommutieren, etwa ein Hamiltonian, der nur aus Anzahloperatoren besteht, werden bei jeder Ordnung exakt synthetisiert, sodass eine höhere Ordnung dort nur Tiefe hinzufügen würde.

Den Pass auf die Gates beschränken, die profitieren

Ein rein diagonaler Operator macht den Punkt konkret, denn Anzahloperatoren kommutieren miteinander. Seine Gruppen sind in erster Ordnung exakt, und die zweite Ordnung verbraucht verschränkende Gates für nichts:

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

Beide Circuits sind exakt, also verbraucht der zweite mehr verschränkende Gates ohne Gewinn. (Die symmetrische Formel zweiter Ordnung wendet die erste Gruppe zweimal mit halber Zeit an, und daher hängen die Kosten davon ab, welche Gruppe zuerst kommt.) Der filter nimmt ein Prädikat über den DAGOpNode und lässt abgelehnte Gates unberührt, sodass ein Pass die Ordnung bei den Gates erhöhen kann, die profitieren, und den Rest in erster Ordnung belässt. Das folgende Prädikat akzeptiert ein Gate nur, wenn sein Operator Bosonen zwischen Moden bewegt, was genau dann der Fall ist, wenn die Gruppen nicht kommutieren:

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

Wenn du den Pass erneut ausführst, ändert sich nichts: Die Faktoren, die er ausgibt, sind atomar.

Atomare Faktoren

Ein als atomic markiertes Evolution hat keine Definition, sodass decompose() es unabhängig davon, wie oft es angewendet wird, an Ort und Stelle lässt. Eine Synthese markiert die Faktoren, die sie ausgibt, auf diese Weise, und ebenso BosonicQDrift für die Terme, die es samplet.

Das ist mehr als Ordentlichkeit. Eine hermitesche Gruppe weiter Term für Term aufzuteilen, ergibt Faktoren, die einzeln nicht hermitesch sind, und die Exponentialfunktion eines nicht hermiteschen Operators ist nicht unitär:

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

Achtung

Jeder Faktor muss für sich hermitesch sein. Ein hermitescher Operator, der in Gruppen aufgeteilt wird, kann nicht hermitesche erzeugen, etwa einen Hopping-Term, der von seinem Adjungierten getrennt ist. Prüfe mit groups_are_hermitian(). Ein Operator ohne Gruppen wird Term für Term aufgeteilt, und die Qubit-Stufe lehnt die Faktoren dann mit Operator contains complex coefficients ab. Die Symmetrisierung einer geraden Ordnung hebt den Fehler eines nicht hermiteschen Faktors teilweise auf, sodass ein falsch gruppierter Operator bei Ordnung 2 besser aussehen kann als bei Ordnung 1 und trotzdem falsch ist.

Einen kurzen Schritt wiederholen

Als Alternative auf Circuit-Ebene baust du einen kurzen Trotter-Schritt und wiederholst ihn. repeat() glättet die Wiederholungen zu einem Circuit und lässt die Grenzen aneinanderstoßen, sodass Passes über sie hinweg wirken können:

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

Übergib insert_barriers=True, um die Wiederholungen getrennt zu halten. Eine Barriere wird auf jedes Qubit der Register abgebildet, die sie berührt, weil eine Kodierung mehrere Qubits pro Mode braucht und eine Mode kein eigenes Abbild pro Qubit hat:

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

Hinweis

repeat hat nichts mit den reps einer Synthesemethode zu tun. Eine Synthese unterteilt eine einzelne Evolution in mehr, kürzere Faktoren. Einen Circuit, der eine Evolution enthält, zu wiederholen, vervielfacht lediglich ihre Gesamtzeit.

Mit einer Formel auf Qubit-Ebene kombinieren

Die beiden Ebenen setzen sich zusammen, und wo beide nähern, wird die Genauigkeit des Ergebnisses von der schwächeren der beiden bestimmt. Die bosonische Ordnung zu erhöhen, während die Formel auf Qubit-Seite erster Ordnung bleibt, bringt wenig. Siehe MapperFnEvolutionSynthesis für die Formel auf Qubit-Seite. Sie behält die vom Mapper gewählte Reihenfolge der Pauli-Terme bei, und so kann group_wise() die Terme jeder Gruppe nebeneinander ausgeben, und simplify() kann stattdessen eine kanonische Reihenfolge festlegen.

Nächste Schritte

• Bosonische Circuits transpilieren erklärt die Stufen, die auf die bosonische Synthese folgen, und wo B2BSynthesis darin sitzt.

• Operatorterme gruppieren: die Operatorstruktur nutzen erklärt, wie du die Gruppen zuweist und prüfst, die eine Produktformel verschachtelt.

• Bose-Hubbard-Dynamik mit wählbarer Kodierung simulieren wendet das Wiederholungsschema auf ein vollständiges Beispiel an.