Bosonische Modi rotieren

Eine ModeRotation mischt bosonische Modi mit einer unitären Matrix $U$, wie es ein Netzwerk aus Strahlteilern und Phasenschiebern tut. Sie heißt passiv, weil sie die Gesamtzahl der Bosonen erhält: Sie verschiebt Bosonen zwischen Modi und erzeugt oder vernichtet nie eines. (Squeezing und Verschiebung sind die aktiven Operationen, und eine Modenrotation kann sie nicht ausdrücken.) Sie bildet die Erzeugungsoperatoren ab als

$$b^\dagger_i \mapsto \sum_j U_{ji} \, b^\dagger_j \, ,$$

was auf dem vollen Fock-Raum die Vielteilchen-Unitäre

$$\exp\left(\sum_{ij} \log(U)_{ij} \, b^\dagger_i b_j\right) \, .$$

ist. Das ist der Gate hinter einem linear-optischen Interferometer: Fock-Zustände hineinschicken, die Modi rotieren und auslesen, wie sich die Bosonen neu verteilen. Wie jeder Gate eines BosonicCircuit (siehe Mit bosonischen Circuits arbeiten) trägt er weder eine Abschneidung noch eine Kodierung. Beide kommen zur Transpile-Zeit hinzu, und diese Anleitung zeigt, warum die Abschneidung bei einer Rotation besondere Sorgfalt erfordert.

Eine Rotation aufbauen

Der Konstruktor nimmt die unitäre $M \times M$-Matrix. Der Gate wirkt auf $M$ Modi. Die Klasse prüft nicht, ob die Matrix unitär ist; das liegt in deiner Verantwortung.

>>> import numpy as np
>>> from qiskit_bosons.circuit.library import ModeRotation
>>>
>>> unitary = np.array([[1, -1], [1, 1]]) / np.sqrt(2)
>>> rotation = ModeRotation(unitary)
>>> rotation.num_modes  # the gate acts on two modes
2

Gängige Rotationen haben eigene Konstruktoren. ModeRotation.beamsplitter() nimmt einen Mischwinkel $\theta$ und eine optionale Phase $\phi$, und $\theta = \pi/4$ ergibt den symmetrischen 50:50-Strahlteiler. ModeRotation.phase_shift() nimmt eine Phase pro Modus.

>>> beamsplitter = ModeRotation.beamsplitter(np.pi / 4)
>>> np.round(beamsplitter.rotation_unitary.real, 4)
array([[ 0.7071, -0.7071],
       [ 0.7071,  0.7071]])
>>> bool(np.allclose(beamsplitter.rotation_unitary, unitary))
True
>>> shift = ModeRotation.phase_shift([0.0, np.pi / 2])
>>> np.round(np.diag(shift.rotation_unitary), 4)
array([1.+0.j, 0.+1.j])

Du kannst auch vom antihermiteschen Generator $K$ ausgehen und den Gate ihn exponenzieren lassen, was die Rotation $\exp(K)$ ergibt. Der Generator wird mit generator() zurückgewonnen.

>>> generator = np.array([[0, -0.3], [0.3, 0]], dtype=complex)
>>> small = ModeRotation.from_generator(generator)
>>> np.round(small.rotation_unitary.real, 4)
array([[ 0.9553, -0.2955],
       [ 0.2955,  0.9553]])
>>> bool(np.allclose(small.generator(), generator))
True

Die Brücke zur Operatoralgebra ist to_operator(), das den hermiteschen Generator $H = \sum_{ij} \mathrm{i} \log(U)_{ij} \, b^\dagger_i b_j$ als BosonOperator zurückgibt, sodass der Gate $e^{-\mathrm{i} H}$ ist. Für den 50:50-Strahlteiler ist das ein Paar von Hopping-Termen mit entgegengesetzten imaginären Koeffizienten. Der Operator ist hermitesch, was ein Evolution-Gate braucht, und er erhält die Teilchenzahl.

>>> operator = beamsplitter.to_operator()
>>> # The stored order of the terms is not specified, so sort them before printing.
>>> sorted((term[0], round(term[1].imag, 4)) for term in operator.iter_terms())
[([(True, 0), (False, 1)], -0.7854), ([(True, 1), (False, 0)], 0.7854)]
>>> operator.is_hermitian()
True
>>> operator.conserves_particle_number()
True

Jeder Term liest sich von links nach rechts als (is_creation, mode), sodass der zweite $\mathrm{i} \tfrac{\pi}{4} \, b^\dagger_1 b_0$ ist.

Auf Fock-Zustände wirken

Eine Rotation verschiebt Bosonen zwischen Modi, daher ist ihre Wirkung an Fock-Zuständen am leichtesten zu sehen. Der Hong-Ou-Mandel-Effekt ist das Standardbeispiel. Schicke je ein Boson in jeden Eingangsmodus eines symmetrischen Strahlteilers. Die beiden Pfade, die die Bosonen in verschiedenen Ausgangsmodi lassen, interferieren destruktiv, sodass der Zustand $\lvert 1, 1 \rangle$ verschwindet und die Bosonen immer zusammen den Strahlteiler verlassen, in $\lvert 2, 0 \rangle$ oder $\lvert 0, 2 \rangle$, jeweils mit Wahrscheinlichkeit einhalb.

Baue den Circuit auf bosonischen Modi auf und senke ihn mit dem Preset-Pass-Manager ab. Die LevelSpec legt an dieser Stelle die Abschneidung und die Kodierung fest. Eine Abschneidung von max_occupation=2 behält die Niveaus 0, 1 und 2 jedes Modus.

>>> from qiskit_bosons.circuit import BosonicCircuit
>>> from qiskit_bosons.circuit.library import InitializeModes
>>> from qiskit_bosons.mappers import LevelEncoding, LevelSpec
>>> from qiskit_bosons.transpiler.presets import generate_preset_occupation_pass_manager
>>>
>>> circuit = BosonicCircuit(2)
>>> circuit.append(InitializeModes([1, 1]), circuit.modes)
>>> circuit.append(ModeRotation.beamsplitter(np.pi / 4), circuit.modes)

Die Pipeline gibt einen gewöhnlichen Qubit-Circuit zurück. Lies die Ergebniswahrscheinlichkeiten aus dem Endzustand ab, hier mit Qiskits Statevector, der bei dieser Größe exakt ist. Jeder Simulator, der Amplituden oder Counts liefert, funktioniert. Ein Fock-Zustand $\lvert n_0, n_1 \rangle$ ist ein einzelner Basiszustand der Rechenbasis, dessen Index sich aus den Bitmustern von LevelEncoding.level_to_bits() ergibt; die Qubits des zweiten Modus folgen auf die des ersten.

>>> from qiskit.quantum_info import Statevector
>>> from qiskit.transpiler import PassManager
>>>
>>> def fock_probabilities(circuit, max_occupation):
...     """Maps each pair of occupations to the probability of that outcome."""
...     encoding = LevelEncoding.UNARY_UNPROJECTED
...     spec = LevelSpec(max_occupation, encoding)
...     pass_manager = generate_preset_occupation_pass_manager(spec)
...     pass_manager.qubit = PassManager()  # keep the evolutions as the pipeline emitted them
...     qubit_circuit = pass_manager.run(circuit)
...     state = Statevector(qubit_circuit)
...     width = encoding.num_qubits(max_occupation + 1)
...     probabilities = {}
...     for n0 in range(max_occupation + 1):
...         for n1 in range(max_occupation + 1):
...             index = encoding.level_to_bits(n0) | (encoding.level_to_bits(n1) << width)
...             probability = round(float(abs(state.data[index]) ** 2), 4)
...             if probability > 0:
...                 probabilities[(n0, n1)] = probability
...     return probabilities
>>>
>>> fock_probabilities(circuit, max_occupation=2)
{(0, 2): 0.5, (2, 0): 0.5}

Das Ergebnis $\lvert 1, 1 \rangle$ hat die Wahrscheinlichkeit null, und die beiden gebündelten Ergebnisse teilen sich die Wahrscheinlichkeit gleichmäßig.

Die Abschneidung wählen

Das vorherige Ergebnis braucht max_occupation=2. Führe denselben Circuit mit einer Abschneidung von 1 aus. Das genügt, um den Eingangszustand $\lvert 1, 1 \rangle$ zu fassen, aber nicht den Zustand, den die Rotation daraus erzeugt und der Gewicht auf $\lvert 2, 0 \rangle$ und $\lvert 0, 2 \rangle$ hat:

>>> fock_probabilities(circuit, max_occupation=1)
{(1, 1): 1.0}

Der Circuit läuft ohne Fehler durch, und die Antwort ist falsch: Der Zustand bleibt mit Wahrscheinlichkeit eins bei $\lvert 1, 1 \rangle$. Eine Rotation mischt Besetzungssektoren innerhalb eines Modus. Zwei Bosonen können im selben Ausgangsmodus landen, sodass der exakte Zustand Gewicht auf der Besetzung 2 in einem Modus hat, der mit einem Boson gestartet ist. Bei einer Abschneidung von 1 existiert dieses Niveau nicht, und der abgeschnittene Generator kann den Zustand nicht dorthin bewegen.

Eine Schranke für die Besetzung des Anfangszustands genügt daher nicht. Die Abschneidung muss die Gesamtzahl der Bosonen abdecken, die jeden einzelnen Modus der Rotation erreichen können. Bei insgesamt $N$ Bosonen ist eine Abschneidung von $N$ die sichere Wahl.

Dieselbe Überlegung gilt für jeden Modus einzeln, sodass eine Mischung unterschiedlicher Abschneidungen pro Modus (siehe Die Abschneidung mit einer LevelSpec beschreiben) eine Rotation leicht falsch abbilden lässt. Ein Modus, der Bosonen empfangen kann, braucht eine Abschneidung, die abdeckt, was er empfangen kann, nicht nur, was er anfangs hält. Die Prüfung LevelSpec.truncation() vergleicht eine Spezifikation mit der unteren Besetzungsschranke eines Operators. Den obigen Fall kann sie nicht erkennen. Jeder Term des Strahlteiler-Generators enthält einen einzelnen Erzeugungsoperator, sodass die Schranke eins ist und die Prüfung bei einer Abschneidung von 1 kein Defizit meldet, obwohl das Ergebnis dort falsch ist:

>>> LevelSpec(1, "unary_unprojected").truncation(beamsplitter.to_operator())
[0, 0]

Eine weitere Warnung ergibt sich aus der Art, wie der Transpiler eine Rotation aufbaut (nächster Abschnitt). Er wendet die Rotation als Folge von Zwei-Modus-Faktoren an, und die Plugin-Dokumentation warnt, dass ein einzelner Faktor ein Niveau bevölkern kann, das die vollständige Rotation nie erreicht. Die abgeschnittene Folge und die abgeschnittene einzelne Evolution sind dann unterschiedliche Näherungen derselben Unitären. Eine Suche über 1.080 zufällige Drei-Modus-Rotationen, mit einem Boson in jedem von zwei Eingangsmodi und ohne Doppelbesetzung im Ausgang, brachte keinen solchen Fall zutage; betrachte dies also als Vorsichtsmaßnahme und nicht als häufigen Fehler. Lass die Abschneidung konvergieren wie bei jedem anderen Circuit, indem du sie erhöhst, bis sich das Ergebnis nicht mehr ändert.

Wie der Transpiler eine Rotation synthetisiert

Die Stufe optimization der Preset-Pipeline zerlegt eine ModeRotation zuerst im bosonischen Raum, mit dem Pass B2BSynthesis und der Methode BosonicClementsDecomposition (registriert als ModeRotation.Clements). Das ist eine Wahl des Presets und keine Anforderung des Frameworks: Ein späterer Abschnitt senkt denselben Gate ohne jede Zerlegung ab. Die Zerlegung exponenziert den vollen Generator nicht in einem Schritt. Sie faktorisiert die Unitäre zuerst in höchstens $M(M-1)/2$ Rotationen zwischen benachbarten Modi plus $M$ Ein-Modus-Phasen, nach dem rechteckigen Netz von Clements et al. [1]. Jeder Faktor wird zu einem atomaren Evolution-Gate, den die Boson-zu-Qubit-Stufe wie jeden anderen abbildet.

Die Faktorisierung selbst ist givens_decomposition(), die von der Zerlegung aufgerufen wird. Die Abbildung von einer unitären Matrix auf die Vielteilchen-Rotation ist ein Gruppenhomomorphismus, sodass das Produkt der Faktor-Rotationen im nicht abgeschnittenen Raum exakt der vollständigen Rotation entspricht. Das Netz fügt also keinen eigenen Fehler hinzu. Sein Vorteil ist die Lokalität. Jede erzeugte Evolution betrifft nur die Pauli-Terme der beiden Moden, die sie mischt (bzw. einer Mode bei einer Phase), statt Terme über alle Moden hinweg.

Transpiliere eine zufällige Rotation über drei Moden und untersuche, was die Pipeline erzeugt hat. Eine allgemeine unitäre $3 \times 3$-Matrix benötigt höchstens drei Rotationen plus drei Phasen, und jeder Faktor ist sein eigenes PauliEvolutionGate. Alle sind auf dem gesamten Register deklariert, aber die darin enthaltenen Operatoren wirken nur auf wenige Qubits:

>>> from qiskit.quantum_info import random_unitary
>>>
>>> random = np.asarray(random_unitary(3, seed=7).data)
>>> spec = LevelSpec(2, "unary_unprojected")
>>> rotation_circuit = BosonicCircuit(3)
>>> rotation_circuit.append(InitializeModes([1, 0, 0]), rotation_circuit.modes)
>>> rotation_circuit.append(ModeRotation(random), rotation_circuit.modes)
>>> pass_manager = generate_preset_occupation_pass_manager(spec)
>>> pass_manager.qubit = PassManager()
>>> qubit_circuit = pass_manager.run(rotation_circuit)
>>> dict(qubit_circuit.count_ops())
{'PauliEvolution': 6, 'x': 3}
>>> for instruction in qubit_circuit.data:
...     if instruction.operation.name == "PauliEvolution":
...         terms = instruction.operation.operator.to_sparse_list()
...         print(len(terms), sorted({qubit for _, indices, _ in terms for qubit in indices}))
64 [0, 1, 2, 3, 4, 5]
64 [3, 4, 5, 6, 7, 8]
64 [0, 1, 2, 3, 4, 5]
2 [1, 2]
2 [4, 5]
2 [7, 8]

Bei drei Niveaus pro Mode belegt jede Mode drei Qubits. Ein Strahlteiler zwischen zwei Moden berührt die sechs Qubits dieser Moden, und eine Phase berührt die Qubits einer Mode.

Die Faktoren sind nur dann exakt, wenn jede Evolution gut synthetisiert wird. Die Terme eines Strahlteiler-Generators kommutieren nicht, daher hinterlässt die standardmäßige einstufige Produktformel in jedem Faktor einen Fehler. Übergib dem Evolution-Plugin, das die Faktoren abbildet, eine Produktformel mit mehr Wiederholungen, um ihn zu verringern:

>>> from qiskit.passmanager import MultiStagePassManager
>>> from qiskit.synthesis import LieTrotter
>>> from qiskit_bosons.circuit.library.synthesis import BosonicClementsDecomposition
>>> from qiskit_bosons.mappers.library import occupation_encoding
>>> from qiskit_bosons.transpiler import (
...     BosonicCircuitToDAG,
...     BosonicPassManager,
...     QuantumDAGToCircuit,
... )
>>> from qiskit_bosons.transpiler.passes import (
...     B2BSynthesis,
...     B2QSynthesis,
...     MapperFnEvolutionSynthesis,
...     OccupationB2QLayout,
...     OccupationInitializeModesSynthesis,
... )
>>>
>>> def with_product_formula(product_formula):
...     synthesis = B2QSynthesis()
...     synthesis.methods["Evolution"] = MapperFnEvolutionSynthesis(
...         lambda operator, num_qubits: occupation_encoding(operator, spec, num_qubits),
...         product_formula,
...     )
...     synthesis.methods["InitializeModes"] = OccupationInitializeModesSynthesis(spec)
...     return MultiStagePassManager(
...         input=BosonicCircuitToDAG(),
...         optimization=BosonicPassManager(
...             B2BSynthesis({"ModeRotation": BosonicClementsDecomposition()})
...         ),
...         layout=OccupationB2QLayout(spec),
...         synthesis=synthesis,
...         output=QuantumDAGToCircuit(),
...     )

Ein einzelnes Boson, das in Mode 0 eintritt, verlässt das System mit der Wahrscheinlichkeit $|U_{j0}|^2$ in Mode $j$. Vergleiche den synthetisierten Circuit für steigende Wiederholungszahlen mit diesem exakten Ergebnis. Zerlege den Circuit zuerst, damit die Produktformel angewendet wird; andernfalls wendet ein exakter Simulator wie Statevector das Gate als vollständige Exponentialfunktion an und verbirgt den Fehler.

>>> encoding = LevelEncoding.UNARY_UNPROJECTED
>>> width = encoding.num_qubits(3)
>>> exact = np.abs(random[:, 0]) ** 2
>>>
>>> def one_boson_probabilities(qubit_circuit):
...     """Probability that the single boson sits in mode j, for each j."""
...     state = Statevector(qubit_circuit)
...     probabilities = []
...     for j in range(3):
...         occupations = [1 if mode == j else 0 for mode in range(3)]
...         index = sum(
...             encoding.level_to_bits(n) << (width * mode) for mode, n in enumerate(occupations)
...         )
...         probabilities.append(abs(state.data[index]) ** 2)
...     return np.array(probabilities)
>>>
>>> for repetitions in (1, 4, 16):
...     lowered = with_product_formula(LieTrotter(reps=repetitions)).run(rotation_circuit)
...     error = np.max(np.abs(one_boson_probabilities(lowered.decompose(reps=2)) - exact))
...     print(f"{repetitions:2d} repetitions: error {error:.1e}")
 1 repetitions: error 1.2e-03
 4 repetitions: error 7.6e-05
16 repetitions: error 4.7e-06

Der Fehler sinkt um etwa den Faktor vier bei jeder Vervierfachung der Wiederholungen, wie es für eine Formel erster Ordnung zu erwarten ist. Die Kosten wachsen proportional, weil jeder Faktor wiederholt wird.

Eine Rotation ohne Zerlegung absenken

Eine Zerlegung im bosonischen Raum ist eine Möglichkeit, eine Rotation abzusenken. (Die Bibliothek enthält keine andere Möglichkeit, aber die Stufe von Bosonen zu Qubits hängt nicht davon ab.) B2QSynthesis sucht ein Plugin anhand des Namens des Gates, und eine ModeRotation, die ganz belassen wurde, erreicht es wie jedes andere Gate. Ein Plugin dafür bildet den Generator von ModeRotation.to_operator() in einem Zug ab und erzeugt ein einzelnes PauliEvolutionGate bei Einheitszeit. Das ist alles:

>>> from qiskit.circuit.library import PauliEvolutionGate
>>> from qiskit_bosons.transpiler.passes.utils import map_node_single_register
>>>
>>> class DirectModeRotationSynthesis:
...     """Maps the generator of the whole rotation, with no bosonic decomposition."""
...
...     def __init__(self, mapper_fn):
...         self.mapper_fn = mapper_fn
...
...     def run(self, in_node, out_dag, *, b2q_layout):
...         modes, qreg = map_node_single_register(in_node, b2q_layout)
...         generator = in_node.op.to_operator().relabel_modes(list(modes))
...         observable = self.mapper_fn(generator, len(qreg))
...         out_dag.apply_operation_back(PauliEvolutionGate(observable, time=1.0), tuple(qreg))

Registriere es für den Namen des Gates in einem Pass Manager, der es nicht zerlegt:

>>> def mapper_fn(operator, num_qubits):
...     return occupation_encoding(operator, spec, num_qubits)
>>>
>>> synthesis = B2QSynthesis()
>>> synthesis.methods["Evolution"] = MapperFnEvolutionSynthesis(mapper_fn)
>>> synthesis.methods["InitializeModes"] = OccupationInitializeModesSynthesis(spec)
>>> synthesis.methods["ModeRotation"] = DirectModeRotationSynthesis(mapper_fn)
>>> direct = MultiStagePassManager(
...     input=BosonicCircuitToDAG(),
...     layout=OccupationB2QLayout(spec),
...     synthesis=synthesis,
...     output=QuantumDAGToCircuit(),
... )
>>> direct_circuit = direct.run(rotation_circuit)
>>> dict(direct_circuit.count_ops())
{'x': 3, 'PauliEvolution': 1}
>>> terms = direct_circuit.data[-1].operation.operator.to_sparse_list()
>>> len(terms), sorted({qubit for _, indices, _ in terms for qubit in indices})
(198, [0, 1, 2, 3, 4, 5, 6, 7, 8])

Die Rotation ist nun eine Evolution über alle neun Qubits, während das Netz sechs hat, von denen keine mehr als sechs umfasst. Die beiden sind nicht dieselbe Näherung. Die Abschneidung wird hier einmal auf den Generator der gesamten Rotation angewendet, im Netz dagegen auf jeden Faktor einzeln (daher die Warnung zum Cutoff oben). Der direkte Weg entspricht bei jedem Cutoff occupation_encoding() des Operators, ist also die Referenz, mit der du ein Netz vergleichst. Der Preis ist die Lokalität: Seine Terme erstrecken sich über alle Moden, und eine Produktformel hat mehr davon zu unterteilen.

Keines von beiden ist allgemein besser, und genau deshalb lohnt es sich, die Wahl offen zu halten. Schreibe ein Plugin wie dieses, um Absenkungen zu vergleichen, eine kodierungsspezifische Rotation zu verwenden oder das Gate an etwas zu übergeben, das es nativ versteht.

Nächste Schritte

• Bosonische Circuits transpilieren erklärt die Stufen, die der voreingestellte Pass Manager ausführt, und wie du das Plugin ersetzt, das ein Gate behandelt.

• Operatordarstellungen abbilden behandelt die Abschneidung, die Kodierungen und die LevelSpec, die das Layout jeder Mode festlegt.

• Mit bosonischen Circuits arbeiten zeigt, wie du bosonische Circuits zusammenstellst und kombinierst.

• Die API-Dokumentation listet qiskit_bosons.circuit.library für die Gates, qiskit_bosons.linalg für die Zerlegung und qiskit_bosons.transpiler.passes für die Synthese-Plugins auf.

[1]

W. R. Clements, P. C. Humphreys, B. J. Metcalf, W. S. Kolthammer, and I. A. Walmsley, Optimal design for universal multiport interferometers, Optica 3, 1460 (2016), arXiv:1603.08788.