ボソンモードの回転

ModeRotation は、ビームスプリッターと位相シフターのネットワークと同様に、ユニタリ行列 $U$ でボソンモードを混合します。これは全ボソン数を保存するためパッシブと呼ばれます。ボソンをモード間で移動させますが、ボソンを生成したり消滅させたりすることはありません。(スクイージングと変位はアクティブな操作であり、モード回転ではそれらを表現できません。)これは生成演算子を次のように写します。

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

これは、完全なFock空間では多体ユニタリになります。

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

これは線形光学干渉計の背後にあるゲートです。Fock状態を入力し、モードを回転させ、ボソンがどのように再分配されるかを読み出します。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() はモードごとに1つの位相を受け取ります。

>>> 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) と読むため、2つ目は $\mathrm{i} \tfrac{\pi}{4} \, b^\dagger_1 b_0$ です。

Fock状態に作用させる

回転はボソンをモード間で移動させるため、その効果はFock状態で最も見やすくなります。標準的な例はHong-Ou-Mandel効果です。平衡なビームスプリッターの各入力モードに1つずつボソンを送ります。ボソンを異なる出力モードに残す2つの経路は破壊的に干渉するため、状態 $\lvert 1, 1 \rangle$ は消え、ボソンは常に一緒に、それぞれ確率2分の1で $\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)

パイプラインは通常の量子ビット回路を返します。最終状態から結果の確率を読み取ります。ここではQiskitの Statevector を使い、このサイズでは厳密です。振幅またはカウントを返すシミュレーターであれば何でも使えます。Fock状態 $\lvert n_0, n_1 \rangle$ は単一の計算基底状態であり、そのインデックスは LevelEncoding.level_to_bits() のビットパターンによって与えられます。2番目のモードの量子ビットは、1番目のモードの量子ビットの後に続きます。

>>> 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$ の結果の確率はゼロで、2つのバンチした結果が確率を等しく分け合います。

カットオフを選ぶ

前の結果には 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}

回路はエラーなしで実行されますが、答えは誤っています。状態は確率1で $\lvert 1, 1 \rangle$ にとどまります。回転は、モード内の占有セクターを混合します。2つのボソンが同じ出力モードに入りうるため、厳密な状態は、1つのボソンから始まったモードの占有数2に重みを持ちます。カットオフが1だとそのレベルは存在せず、切り捨てられた生成子は状態をそこへ移動させることができません。

したがって、初期状態の占有数に上限があるだけでは不十分です。カットオフは、回転の任意の単一のモードに到達しうる全ボソン数をカバーしなければなりません。合計 $N$ 個のボソンの場合、安全な選択はカットオフ $N$ です。

同じ推論がすべてのモードに個別に当てはまるため、モードごとに異なるカットオフを混在させると(LevelSpecで切り捨てを記述するを参照)、回転を誤ってマッピングしやすくなります。ボソンを受け取りうるモードには、最初に保持している数だけでなく、受け取りうる数をカバーするカットオフが必要です。チェック LevelSpec.truncation() は、スペックを演算子の占有数の下限と比較します。上のケースは検出できません。ビームスプリッターの生成子の各項は単一の生成演算子を持つため、下限は1であり、カットオフ1での結果が誤っているにもかかわらず、チェックは不足を報告しません。

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

さらなる注意点は、トランスパイラーが回転を構築する方法(次のセクション)から生じます。これは回転を2モード因子の列として適用し、プラグインのドキュメントは、単一の因子が完全な回転では決して到達しないレベルを占有させうると警告しています。その場合、切り捨てられた列と切り捨てられた単一の時間発展は、同じユニタリの異なる近似になります。2つの入力モードにそれぞれ1つのボソンを置き、出力に二重占有がない1,080個のランダムな3モード回転を探索しましたが、そのようなケースは見つからなかったため、これは一般的な失敗ではなく予防策として扱ってください。他の回路と同様に、結果が変化しなくなるまでカットオフを増やして収束させてください。

トランスパイラーが回転を合成する方法

プリセットパイプラインの optimization ステージは、まず B2BSynthesis パスと BosonicClementsDecomposition メソッド(ModeRotation.Clements として登録)を使って、ボソン空間で ModeRotation を分解します。これはプリセットの選択であり、フレームワークの要件ではありません。後のセクションでは、分解を一切行わずに同じゲートを下位に変換します。この分解は、完全な生成子を一度に指数化するのではありません。まず、Clementsらの矩形メッシュに従って、ユニタリを最大 $M(M-1)/2$ 個の隣接するモード間の回転と $M$ 個の単一モード位相に因数分解します [1]。各因子はアトミックな Evolution ゲートになり、ボソンから量子ビットへのステージが他のゲートと同様にそれをマッピングします。

因数分解そのものは givens_decomposition() で、分解処理がこれを呼び出します。ユニタリから多体回転への写像は群準同型であるため、各因子の回転の積は、切り捨てを行わない空間での全体の回転と厳密に等しくなります。したがって、メッシュ自体が新たな誤差を加えることはありません。利点は局所性です。出力される各時間発展は、混合する2つのモード（位相の場合は1つのモード）のパウリ項だけを含み、すべてのモードにまたがる項は含みません。

ランダムな3モード回転をトランスパイルし、パイプラインが出力した内容を調べます。一般的な $3 \times 3$ ユニタリには、最大3つの回転と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]

各モードに3つのレベルがある場合、各モードは3つの量子ビットを占めます。2つのモード間のビームスプリッターは、それらのモードの6つの量子ビットに作用し、位相は1つのモードの量子ビットに作用します。

各時間発展が適切に合成されて初めて、因子は厳密になります。ビームスプリッター生成子の各項は可換ではないため、デフォルトの1ステップの積公式では、すべての因子に誤差が残ります。誤差を減らすには、因子を対応付ける 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に入った1つのボソンは、確率 $|U_{j0}|^2$ でモード $j$ から出ていきます。繰り返し回数を増やしながら、合成された回路をこの厳密な結果と比較します。積公式が適用されるように、最初に回路を分解してください。分解しないと、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

繰り返し回数を4倍にするごとに誤差は約4分の1になります。これは1次の公式として想定どおりです。すべての因子が繰り返されるため、コストは比例して増加します。

分解せずに回転を低レベル化する

ボソン空間での分解は、回転を低レベル化する方法の1つです。（ライブラリには他の方法は含まれていませんが、ボソンから量子ビットへの段階はそれに依存しません。）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])

回転は、メッシュでは6つなのに対し、ここでは9つすべての量子ビットにわたる1つの時間発展になります。メッシュの各因子は6量子ビットを超えて広がることはありません。この2つは同じ近似ではありません。切り捨ては、ここでは回転全体の生成子に1回だけ適用されますが、メッシュでは各因子に個別に適用されます（上記のカットオフに関する警告はそのためです）。直接的な経路は、どのカットオフでも演算子の 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.