Обертання бозонних мод

ModeRotation змішує бозонні моди унітарною матрицею $U$, як це робить мережа світлоподільників і фазозсувників. Його названо пасивним, бо воно зберігає повну кількість бозонів: переміщує бозони між модами і ніколи не створює та не знищує жодного. (Стиснення та зсув — це активні операції, і обертання мод не може їх виразити.) Воно відображає оператори народження як

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

що на повному просторі Фока є багаточастинковим унітарним оператором

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

Це вентиль, що стоїть за лінійно-оптичним інтерферометром: подай стани Фока, оберни моди й прочитай, як бозони перерозподілилися. Як і кожен вентиль BosonicCircuit (див. Робота з бозонними схемами), він не несе ні обрізання, ні кодування. Обидва з'являються під час транспіляції, і цей посібник показує, чому обрізання потребує особливої уваги для обертання.

Побудова обертання

Конструктор приймає унітарну матрицю $M \times M$. Вентиль діє на $M$ мод. Клас не перевіряє, чи матриця унітарна, тож це твоя відповідальність.

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

Поширені обертання мають окремі конструктори. ModeRotation.beamsplitter() приймає кут змішування $\theta$ і необов'язкову фазу $\phi$, а $\theta = \pi/4$ дає збалансований світлоподільник 50:50. ModeRotation.phase_shift() приймає одну фазу на моду.

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

Можна також почати з антиермітового генератора $K$ і дозволити вентилю його експоненціювати, що дає обертання $\exp(K)$. Генератор відновлюється за допомогою 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

Місток до алгебри операторів — to_operator(), який повертає ермітів генератор $H = \sum_{ij} \mathrm{i} \log(U)_{ij} \, b^\dagger_i b_j$ як BosonOperator, так що вентиль дорівнює $e^{-\mathrm{i} H}$. Для світлоподільника 50:50 це пара доданків стрибка з протилежними уявними коефіцієнтами. Оператор ермітів, чого й потребує вентиль Evolution, і зберігає кількість частинок.

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

Кожен доданок читається як (is_creation, mode) зліва направо, тож другий — це $\mathrm{i} \tfrac{\pi}{4} \, b^\dagger_1 b_0$.

Дія на стани Фока

Обертання переміщує бозони між модами, тож його дію найлегше побачити на станах Фока. Ефект Хонга-Оу-Менделя — стандартний приклад. Подай по одному бозону на кожну вхідну моду збалансованого світлоподільника. Два шляхи, що залишають бозони в різних вихідних модах, інтерферують деструктивно, тому стан $\lvert 1, 1 \rangle$ зникає, і бозони завжди залишають світлоподільник разом, у $\lvert 2, 0 \rangle$ або $\lvert 0, 2 \rangle$, кожен з імовірністю одна друга.

Побудуй схему на бозонних модах і знизь її рівень за допомогою передналаштованого менеджера проходів. LevelSpec на цьому етапі фіксує обрізання та кодування. Обрізання max_occupation=2 зберігає рівні 0, 1 і 2 кожної моди.

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

Конвеєр повертає звичайну кубітну схему. Прочитай ймовірності результатів із кінцевого стану, тут за допомогою Statevector із Qiskit, який точний за такого розміру. Підійде будь-який симулятор, що повертає амплітуди або лічильники. Стан Фока $\lvert n_0, n_1 \rangle$ — це один обчислювальний базисний стан, чий індекс задається бітовими шаблонами LevelEncoding.level_to_bits(); кубіти другої моди йдуть після кубітів першої.

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

Результат $\lvert 1, 1 \rangle$ має нульову ймовірність, а два згруповані результати ділять ймовірність порівну.

Вибір обрізання

Попередній результат потребує max_occupation=2. Запусти ту саму схему з обрізанням на 1. Цього досить, щоб вмістити вхідний стан $\lvert 1, 1 \rangle$, але не стан, який з нього створює обертання, що має вагу на $\lvert 2, 0 \rangle$ і $\lvert 0, 2 \rangle$:

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

Схема виконується без жодної помилки, а відповідь хибна: стан залишається в $\lvert 1, 1 \rangle$ з імовірністю одиниця. Обертання змішує сектори заповнення всередині моди. Два бозони можуть опинитися в одній вихідній моді, тож точний стан має вагу на заповненні 2 у моді, що починала з одного бозона. З обрізанням на 1 цього рівня не існує, і обрізаний генератор не може перемістити туди стан.

Тому обмеження на заповнення початкового стану недостатнє. Обрізання має покривати повну кількість бозонів, що може потрапити в будь-яку одну моду обертання. Для $N$ бозонів загалом безпечний вибір — обрізання $N$.

Те саме міркування застосовується до кожної моди окремо, тож суміш різних обрізань на моду (див. Опис обрізання за допомогою LevelSpec) робить обертання легким для хибного відображення. Мода, що може отримувати бозони, потребує обрізання, яке покриває те, що вона може отримати, а не лише те, що вона містить спочатку. Перевірка LevelSpec.truncation() порівнює специфікацію з нижньою межею заповнення оператора. Вона не може позначити випадок вище. Кожен доданок генератора світлоподільника містить один оператор народження, тож межа дорівнює одиниці, і перевірка не повідомляє про нестачу за обрізання на 1, хоча результат там хибний:

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

Ще одне застереження випливає зі способу, яким транспайлер будує обертання (наступний розділ). Він застосовує обертання як послідовність двомодових множників, і документація плагіна попереджає, що один множник може заповнити рівень, якого повне обертання ніколи не досягає. Тоді обрізана послідовність і обрізана єдина еволюція є різними наближеннями того самого унітарного оператора. Пошук серед 1080 випадкових тримодових обертань з одним бозоном у кожній із двох вхідних мод і без подвійного заповнення на виході не виявив такого випадку, тож сприймай це як застережний захід, а не поширений збій. Досягай збіжності обрізання, як і для будь-якої іншої схеми, збільшуючи його, доки результат не перестане змінюватися.

Як транспайлер синтезує обертання

Етап optimization передналаштованого конвеєра спочатку розкладає ModeRotation у бозонному просторі за допомогою проходу B2BSynthesis і методу BosonicClementsDecomposition (зареєстрованого як ModeRotation.Clements). Це вибір передналаштування, а не вимога фреймворку: пізніший розділ знижує рівень того самого вентиля без жодного розкладу. Розклад не експоненціює повний генератор за один крок. Спочатку він розкладає унітарний оператор на щонайбільше $M(M-1)/2$ обертань між сусідніми модами плюс $M$ однодомових фаз за прямокутною сіткою Клементса та ін. [1]. Кожен множник стає атомарним вентилем Evolution, який етап відображення бозонів на кубіти відображає як і будь-який інший.

Сама факторизація — це givens_decomposition(), яку викликає декомпозиція. Відображення унітарної матриці в багаточастинкове обертання є гомоморфізмом груп, тому добуток обертань-множників точно дорівнює повному обертанню на неусіченому просторі. Отже, сітка не додає власної похибки. Її перевага — локальність. Кожна згенерована еволюція містить лише паулієві члени двох мод, які вона змішує (або однієї моди для фази), а не члени за всіма модами.

Транспілюй випадкове обертання трьох мод і подивися, що згенерував конвеєр. Типова унітарна матриця $3 \times 3$ потребує щонайбільше трьох обертань плюс трьох фаз, і кожен множник є окремим PauliEvolutionGate. Усі вони оголошені на всьому регістрі, але оператори всередині діють лише на кілька кубітів:

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

За трьох рівнів на моду кожна мода займає три кубіти. Світлоподільник між двома модами зачіпає шість кубітів цих мод, а фаза зачіпає кубіти однієї моди.

Множники точні лише тоді, коли кожна еволюція синтезована добре. Члени генератора світлоподільника не комутують, тому типова однокрокова формула добутку залишає похибку в кожному множнику. Щоб зменшити її, передай формулу добутку з більшою кількістю повторень до плагіна Evolution, який відображає множники:

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

Один бозон, що входить у моду 0, виходить у моді $j$ з імовірністю $|U_{j0}|^2$. Порівняй синтезовану схему з цим точним результатом для зростаючої кількості повторень. Спочатку декомпонуй схему, щоб застосувалася формула добутку; без цього точний симулятор, як-от Statevector, застосовує вентиль як повну експоненту й приховує похибку.

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

Похибка зменшується приблизно в чотири рази за кожне чотириразове збільшення кількості повторень, як і очікується для формули першого порядку. Вартість зростає пропорційно, бо кожен множник повторюється.

Знизь рівень обертання без декомпозиції

Декомпозиція в бозонному просторі — це один зі способів знизити рівень обертання. (Бібліотека не містить іншого способу, але етап «бозон у кубіт» від нього не залежить.) B2QSynthesis шукає плагін за назвою вентиля, і ModeRotation, залишений цілим, потрапляє до нього як і будь-який інший вентиль. Плагін для нього відображає генератор ModeRotation.to_operator() за один раз і видає єдиний PauliEvolutionGate за одиничного часу. Це все:

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

Зареєструй його за назвою вентиля в менеджері проходів, який його не декомпонує:

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

Тепер обертання — це одна еволюція на всіх дев'яти кубітах, тоді як сітка має шість, жодна з яких не охоплює більше шести. Це не однакові наближення. Усічення тут застосовується один раз до генератора всього обертання, а в сітці — до кожного множника окремо (звідси попередження про обрізання вище). Прямий шлях дорівнює occupation_encoding() оператора за будь-якого обрізання, тому він є еталоном для порівняння сітки. Ціна — локальність: його члени охоплюють усі моди, і формулі добутку є більше що підрозділяти.

Жоден із них не кращий загалом, і саме тому варто залишити вибір відкритим. Напиши плагін, подібний до цього, щоб порівнювати зниження рівня, використати обертання, специфічне для кодування, або передати вентиль чомусь, що розуміє його нативно.

Наступні кроки

• Транспіляція бозонних схем пояснює етапи, які виконує попередньо налаштований менеджер проходів, і як замінити плагін, що обробляє вентиль.

• Відображення представлень операторів охоплює усічення, кодування та LevelSpec, що фіксує розкладку кожної моди.

• Робота з бозонними схемами показує, як збирати та поєднувати бозонні схеми.

• Документація API містить qiskit_bosons.circuit.library для вентилів, qiskit_bosons.linalg для декомпозиції та qiskit_bosons.transpiler.passes для плагінів синтезу.

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