Rotacja modów bozonowych

ModeRotation miesza mody bozonowe za pomocą macierzy unitarnej $U$, tak jak robi to sieć dzielników wiązki i przesuwników fazy. Nazywa się ją pasywną, ponieważ zachowuje całkowitą liczbę bozonów: przenosi bozony między modami i nigdy żadnego nie tworzy ani nie niszczy. (Ściskanie i przesunięcie to operacje aktywne, a rotacja modów nie potrafi ich wyrazić.) Odwzorowuje operatory kreacji jako

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

co na pełnej przestrzeni Focka jest unitarnością wielociałową

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

To bramka stojąca za interferometrem liniowo-optycznym: wysyłasz stany Focka, rotujesz mody i odczytujesz, jak bozony się redystrybuują. Jak każda bramka BosonicCircuit (zobacz Praca z obwodami bozonowymi) nie niesie ani obcięcia, ani kodowania. Oba wchodzą w czasie transpilacji, a ten przewodnik pokazuje, dlaczego obcięcie wymaga szczególnej uwagi w przypadku rotacji.

Zbuduj rotację

Konstruktor przyjmuje macierz unitarną $M \times M$. Bramka działa na $M$ modach. Klasa nie sprawdza, czy macierz jest unitarna, więc to twoja odpowiedzialność.

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

Typowe rotacje mają dedykowane konstruktory. ModeRotation.beamsplitter() przyjmuje kąt mieszania $\theta$ i opcjonalną fazę $\phi$, a $\theta = \pi/4$ daje zrównoważony dzielnik wiązki 50:50. ModeRotation.phase_shift() przyjmuje jedną fazę na mod.

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

Możesz też zacząć od antyhermitowskiego generatora $K$ i pozwolić bramce go wyeksponencjować, co daje rotację $\exp(K)$. Generator odzyskujesz za pomocą generator().

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

Mostem do algebry operatorów jest to_operator(), który zwraca hermitowski generator $H = \sum_{ij} \mathrm{i} \log(U)_{ij} \, b^\dagger_i b_j$ jako BosonOperator, tak że bramka to $e^{-\mathrm{i} H}$. Dla dzielnika wiązki 50:50 jest to para członów przeskoku z przeciwnymi współczynnikami urojonymi. Operator jest hermitowski, czego potrzebuje bramka Evolution, i zachowuje liczbę cząstek.

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

Każdy człon odczytuje (is_creation, mode) od lewej do prawej, więc drugi to $\mathrm{i} \tfrac{\pi}{4} \, b^\dagger_1 b_0$.

Działanie na stanach Focka

Rotacja przenosi bozony między modami, więc jej skutek najłatwiej zobaczyć na stanach Focka. Standardowym przykładem jest efekt Honga-Oua-Mandela. Wyślij po jednym bozonie do każdego wejściowego modu zrównoważonego dzielnika wiązki. Dwie ścieżki, które zostawiają bozony w różnych modach wyjściowych, interferują destruktywnie, więc stan $\lvert 1, 1 \rangle$ znika, a bozony zawsze wychodzą razem, w $\lvert 2, 0 \rangle$ lub $\lvert 0, 2 \rangle$, każde z prawdopodobieństwem jedna druga.

Zbuduj obwód na modach bozonowych i sprowadź go za pomocą predefiniowanego menedżera przebiegów. LevelSpec ustala w tym miejscu obcięcie i kodowanie. Obcięcie max_occupation=2 zachowuje poziomy 0, 1 i 2 każdego modu.

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

Potok zwraca zwykły obwód kubitowy. Odczytaj prawdopodobieństwa wyników ze stanu końcowego, tu za pomocą Statevector z Qiskit, który jest dokładny przy tym rozmiarze. Działa każdy symulator zwracający amplitudy lub zliczenia. Stan Focka $\lvert n_0, n_1 \rangle$ to pojedynczy stan bazy obliczeniowej, którego indeks wynika ze wzorców bitowych LevelEncoding.level_to_bits(); kubity drugiego modu następują po kubitach pierwszego.

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

Wynik $\lvert 1, 1 \rangle$ ma prawdopodobieństwo zero, a dwa wyniki zgrupowane dzielą prawdopodobieństwo po równo.

Wybór obcięcia

Poprzedni wynik wymaga max_occupation=2. Uruchom ten sam obwód z obcięciem równym 1. To wystarcza, by pomieścić stan wejściowy $\lvert 1, 1 \rangle$, ale nie stan, który rotacja z niego wytwarza i który ma wagę na $\lvert 2, 0 \rangle$ i $\lvert 0, 2 \rangle$:

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

Obwód działa bez żadnego błędu, a odpowiedź jest błędna: stan pozostaje w $\lvert 1, 1 \rangle$ z prawdopodobieństwem jeden. Rotacja miesza sektory obsadzeń wewnątrz modu. Dwa bozony mogą skończyć w tym samym modzie wyjściowym, więc dokładny stan ma wagę na obsadzeniu 2 w modzie, który zaczynał z jednym bozonem. Przy obcięciu równym 1 ten poziom nie istnieje, a obcięty generator nie może tam przenieść stanu.

Ograniczenie obsadzenia stanu początkowego nie wystarcza więc. Obcięcie musi obejmować całkowitą liczbę bozonów, które mogą trafić do dowolnego pojedynczego modu rotacji. Dla łącznie $N$ bozonów bezpiecznym wyborem jest obcięcie równe $N$.

To samo rozumowanie dotyczy każdego modu osobno, więc mieszanka różnych obcięć na mod (zobacz Opisz obcięcie za pomocą LevelSpec) sprawia, że rotację łatwo zmapować niepoprawnie. Mod, który może otrzymać bozony, potrzebuje obcięcia obejmującego to, co może otrzymać, a nie tylko to, co zawiera na początku. Sprawdzenie LevelSpec.truncation() porównuje specyfikację z dolnym ograniczeniem obsadzenia operatora. Nie potrafi wychwycić powyższego przypadku. Każdy człon generatora dzielnika wiązki zawiera pojedynczy operator kreacji, więc ograniczenie wynosi jeden, a sprawdzenie nie zgłasza niedoboru przy obcięciu równym 1, choć wynik jest tam błędny:

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

Dalsze ostrzeżenie wynika ze sposobu, w jaki transpilator buduje rotację (następna sekcja). Stosuje rotację jako sekwencję czynników dwumodowych, a dokumentacja wtyczki ostrzega, że pojedynczy czynnik może obsadzić poziom, którego pełna rotacja nigdy nie osiąga. Obcięta sekwencja i obcięta pojedyncza ewolucja są wtedy różnymi przybliżeniami tej samej unitarności. Przeszukanie 1080 losowych rotacji trzymodowych, z jednym bozonem w każdym z dwóch modów wejściowych i bez podwójnego obsadzenia na wyjściu, nie wykazało takiego przypadku, więc traktuj to jako środek ostrożności, a nie częstą awarię. Doprowadź obcięcie do zbieżności tak jak dla każdego innego obwodu, zwiększając je, aż wynik przestanie się zmieniać.

Jak transpilator syntezuje rotację

Etap optimization predefiniowanego potoku najpierw rozkłada ModeRotation w przestrzeni bozonowej, za pomocą przebiegu B2BSynthesis i metody BosonicClementsDecomposition (zarejestrowanej jako ModeRotation.Clements). To wybór presetu, a nie wymóg frameworka: późniejsza sekcja sprowadza tę samą bramkę bez żadnego rozkładu. Rozkład nie eksponencjuje pełnego generatora w jednym kroku. Najpierw rozkłada unitarność na co najwyżej $M(M-1)/2$ rotacji między sąsiednimi modami oraz $M$ faz jednomodowych, zgodnie z prostokątną siatką Clementsa i in. [1]. Każdy czynnik staje się atomową bramką Evolution, którą etap bozon-kubit mapuje jak każdą inną.

Samą faktoryzacją jest givens_decomposition(), którą wywołuje dekompozycja. Odwzorowanie macierzy unitarnej na rotację wielociałową jest homomorfizmem grup, więc iloczyn rotacji czynnikowych dokładnie równa się pełnej rotacji w nieobciętej przestrzeni. Siatka nie wnosi zatem żadnego własnego błędu. Jej zaletą jest lokalność. Każda wygenerowana ewolucja obejmuje tylko człony Pauliego dwóch mieszanych modów (lub jednego modu w przypadku fazy), a nie człony ze wszystkich modów.

Dokonaj transpilacji losowej rotacji trzech modów i sprawdź, co wygenerował potok. Ogólna macierz unitarna $3 \times 3$ wymaga maksymalnie trzech rotacji oraz trzech faz, a każdy czynnik jest osobną bramką PauliEvolutionGate. Wszystkie są zadeklarowane na całym rejestrze, ale operatory w ich wnętrzu działają tylko na kilku kubitach:

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

Przy trzech poziomach na mod każdy mod zajmuje trzy kubity. Dzielnik wiązki między dwoma modami dotyka sześciu kubitów tych modów, a faza dotyka kubitów jednego modu.

Czynniki są dokładne dopiero wtedy, gdy każda ewolucja jest dobrze zsyntetyzowana. Człony generatora dzielnika wiązki nie komutują, więc domyślna jednokrokowa formuła produktowa pozostawia błąd w każdym czynniku. Przekaż formułę produktową z większą liczbą powtórzeń do wtyczki Evolution, która odwzorowuje czynniki, aby go zmniejszyć:

>>> 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(),
...     )

Pojedynczy bozon wchodzący do modu 0 wychodzi w modzie $j$ z prawdopodobieństwem $|U_{j0}|^2$. Porównaj zsyntetyzowany obwód z tym dokładnym wynikiem dla rosnącej liczby powtórzeń. Najpierw zdekomponuj obwód, aby zastosowana została formuła produktowa; bez tego dokładny symulator, taki jak Statevector, stosuje bramkę jako jej pełną eksponentę i ukrywa błąd.

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

Błąd maleje około czterokrotnie przy każdym czterokrotnym zwiększeniu liczby powtórzeń, zgodnie z oczekiwaniami dla formuły pierwszego rzędu. Koszt rośnie proporcjonalnie, ponieważ każdy czynnik jest powtarzany.

Obniż rotację bez jej dekomponowania

Dekompozycja w przestrzeni bozonowej to jeden ze sposobów obniżenia rotacji. (Biblioteka nie zawiera innego sposobu, ale etap bozon-kubit od niego nie zależy.) B2QSynthesis wyszukuje wtyczkę po nazwie bramki, a ModeRotation pozostawiona w całości trafia do niej jak każda inna bramka. Wtyczka dla niej odwzorowuje generator ModeRotation.to_operator() za jednym razem i emituje pojedynczą bramkę PauliEvolutionGate w czasie jednostkowym. To wszystko:

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

Zarejestruj ją pod nazwą bramki w menedżerze przebiegów, który jej nie dekomponuje:

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

Rotacja jest teraz jedną ewolucją na wszystkich dziewięciu kubitach, podczas gdy siatka ma sześć, z których żadna nie obejmuje więcej niż sześciu. Oba podejścia nie są tym samym przybliżeniem. Obcięcie jest tu stosowane jednokrotnie do generatora całej rotacji, a w siatce osobno do każdego czynnika (stąd ostrzeżenie o odcięciu powyżej). Ścieżka bezpośrednia równa się occupation_encoding() operatora przy dowolnym odcięciu, więc jest punktem odniesienia do porównania siatki. Jej kosztem jest lokalność: jej człony obejmują wszystkie mody, a formuła produktowa ma więcej członów do podziału.

Żadne z nich nie jest ogólnie lepsze, i właśnie dlatego warto zostawić wybór otwarty. Napisz wtyczkę podobną do tej, aby porównać obniżenia, użyć rotacji specyficznej dla kodowania albo przekazać bramkę czemuś, co rozumie ją natywnie.

Kolejne kroki

• Transpilacja obwodów bozonowych wyjaśnia etapy uruchamiane przez predefiniowany menedżer przebiegów oraz sposób zastępowania wtyczki obsługującej bramkę.

• Odwzorowywanie reprezentacji operatorów omawia obcięcie, kodowania oraz LevelSpec, który ustala układ każdego modu.

• Praca z obwodami bozonowymi pokazuje, jak składać i łączyć obwody bozonowe.

• Dokumentacja API opisuje qiskit_bosons.circuit.library dla bramek, qiskit_bosons.linalg dla dekompozycji oraz qiskit_bosons.transpiler.passes dla wtyczek syntezy.

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