보손 모드 회전하기

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()는 모드마다 하나의 위상을 받습니다.

>>> 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$입니다.

Fock 상태에 작용시키기

회전은 보손을 모드 사이로 옮기므로, 그 효과는 Fock 상태에서 가장 쉽게 볼 수 있습니다. Hong-Ou-Mandel 효과가 대표적인 예입니다. 균형 잡힌 빔 스플리터의 각 입력 모드에 보손을 하나씩 보냅니다. 보손이 서로 다른 출력 모드에 남게 되는 두 경로는 상쇄 간섭을 일으키므로, 상태 $\lvert 1, 1 \rangle$은 사라지고 보손은 항상 함께 $\lvert 2, 0 \rangle$ 또는 $\lvert 0, 2 \rangle$로 나가며, 각각의 확률은 2분의 1입니다.

보손 모드에서 회로를 만들고 프리셋 패스 매니저로 낮추세요. 이 시점에서 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()의 비트 패턴으로 주어집니다. 두 번째 모드의 큐비트는 첫 번째 모드의 큐비트 뒤에 옵니다.

>>> 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$ 결과의 확률은 0이며, 뭉친 두 결과가 확률을 똑같이 나눠 가집니다.

절단 한계 선택하기

앞의 결과에는 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이면 그 준위가 존재하지 않으며, 절단된 생성자는 상태를 그곳으로 옮길 수 없습니다.

따라서 초기 상태의 점유수에 대한 한계만으로는 충분하지 않습니다. 절단 한계는 회전의 어느 단일 모드에든 도달할 수 있는 총 보손 수를 포괄해야 합니다. 총 $N$개의 보손에 대해 안전한 선택은 절단 한계 $N$입니다.

같은 논리가 모든 모드에 각각 적용되므로, 모드마다 서로 다른 절단 한계를 섞어 쓰면(LevelSpec으로 절단 기술하기 참고) 회전을 잘못 매핑하기 쉽습니다. 보손을 받을 수 있는 모드는 처음에 가진 것만이 아니라 받을 수 있는 것을 포괄하는 절단 한계가 필요합니다. 검사 LevelSpec.truncation()은 스펙을 연산자의 점유수 하한과 비교합니다. 이는 위의 경우를 잡아내지 못합니다. 빔 스플리터 생성자의 각 항은 생성 연산자를 하나만 가지므로 하한은 1이고, 이 검사는 절단 한계 1에서 부족분이 없다고 보고하지만, 거기서의 결과는 틀립니다.

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

트랜스파일러가 회전을 구성하는 방식(다음 절)에서 한 가지 추가 주의 사항이 따라옵니다. 트랜스파일러는 회전을 2모드 인자의 시퀀스로 적용하며, 플러그인 문서는 하나의 인자가 전체 회전이 도달하지 않는 준위를 채울 수 있다고 경고합니다. 그러면 절단된 시퀀스와 절단된 단일 발전은 같은 유니터리에 대한 서로 다른 근사가 됩니다. 두 입력 모드에 보손이 하나씩 있고 출력에 이중 점유가 없는 1,080개의 무작위 3모드 회전을 탐색했지만 그런 경우는 나오지 않았으므로, 이를 흔한 실패가 아니라 예방 조치로 받아들이세요. 다른 회로와 마찬가지로 결과가 더 이상 변하지 않을 때까지 절단 한계를 늘려 가며 수렴시키세요.

트랜스파일러가 회전을 합성하는 방식

프리셋 파이프라인의 optimization 단계는 ModeRotation을 먼저 보손 공간에서, B2BSynthesis 패스와 BosonicClementsDecomposition 방법(ModeRotation.Clements로 등록됨)으로 분해합니다. 이는 프레임워크의 요구 사항이 아니라 프리셋의 선택입니다. 뒤의 절에서는 어떤 분해도 없이 같은 게이트를 낮춥니다. 이 분해는 전체 생성자를 한 번에 지수화하지 않습니다. 먼저 유니터리를 인접한 모드 사이의 최대 $M(M-1)/2$개의 회전과 $M$개의 단일 모드 위상으로 분해하며, 이는 Clements 등의 직사각형 메시를 따릅니다 [1]. 각 인자는 원자적인 Evolution 게이트가 되며, 보손-큐비트 단계가 이를 다른 게이트와 똑같이 매핑합니다.

인수분해 자체는 givens_decomposition()이며, 분해 과정에서 이 함수를 호출합니다. 유니타리에서 다체 회전으로 가는 사상은 군 준동형사상이므로, 인자 회전들의 곱은 절단되지 않은 공간에서의 전체 회전과 정확히 같습니다. 따라서 메시는 자체적인 오차를 더하지 않습니다. 이 방식의 장점은 국소성입니다. 생성되는 각 시간 발전은 모든 모드에 걸친 항이 아니라, 섞이는 두 모드(위상의 경우 한 모드)의 파울리 항만 포함합니다.

임의의 3모드 회전을 트랜스파일하고 파이프라인이 내보낸 결과를 살펴보세요. 일반적인 $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으로 들어온 단일 보손은 확률 $|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

1차 공식에서 예상되는 대로, 반복 횟수가 4배 늘어날 때마다 오차는 약 4분의 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])

이제 회전은 아홉 개 큐비트 전체에 걸친 하나의 시간 발전이 되며, 메시는 여섯 개이고 그중 어느 것도 여섯 개를 넘는 범위에 걸치지 않습니다. 두 방식은 같은 근사가 아닙니다. 절단은 여기서는 전체 회전의 생성자에 한 번 적용되고, 메시에서는 모든 인자에 각각 적용됩니다(위의 컷오프에 대한 경고가 이 때문입니다). 직접 경로는 어떤 컷오프에서도 연산자의 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.