סיבוב מצבים בוזוניים

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 (ראו עבודה עם מעגלים בוזוניים), הוא לא נושא קטימה ולא קידוד. שניהם נכנסים בזמן ה-transpile, והמדריך הזה מראה למה הקטימה דורשת זהירות מיוחדת בסיבוב.

בניית סיבוב

הבנאי מקבל את המטריצה האוניטרית $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$.

פעולה על מצבי פוק

סיבוב מזיז בוזונים בין מצבים, ולכן ההשפעה שלו קלה לראות על מצבי פוק. אפקט Hong-Ou-Mandel הוא הדוגמה הסטנדרטית. שלחו בוזון אחד לכל מצב קלט של מפצל קרן מאוזן. שני המסלולים שמשאירים את הבוזונים במצבי פלט שונים מתאבכים הרסנית, ולכן המצב $\lvert 1, 1 \rangle$ נעלם והבוזונים תמיד יוצאים יחד, ב-$\lvert 2, 0 \rangle$ או ב-$\lvert 0, 2 \rangle$, כל אחד בהסתברות חצי.

בנו את המעגל על מצבים בוזוניים והורידו אותו עם מנהל המעברים (pass manager) המוגדר מראש. 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() משווה spec לחסם התחתון של אכלוס של אופרטור. היא לא יכולה לסמן את המקרה שלמעלה. כל איבר של מחולל מפצל הקרן מחזיק אופרטור יצירה בודד, ולכן החסם הוא אחד והבדיקה לא מדווחת על חסר בקטימה של 1, אף שהתוצאה שם שגויה:

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

אזהרה נוספת נובעת מהאופן שבו ה-Transpiler בונה סיבוב (בסעיף הבא). הוא מיישם את הסיבוב כסדרה של גורמים דו-מצביים, ותיעוד התוסף מזהיר שגורם בודד יכול לאכלס רמה שהסיבוב השלם לעולם לא מגיע אליה. הסדרה הקטומה וההתפתחות הבודדת הקטומה הן אז קירובים שונים של אותו אוניטרי. חיפוש על פני 1,080 סיבובים אקראיים של שלושה מצבים, עם בוזון אחד בכל אחד משני מצבי קלט וללא אכלוס כפול בפלט, לא מצא מקרה כזה, ולכן התייחסו לזה כאמצעי זהירות ולא ככשל נפוץ. הביאו את הקטימה להתכנסות כמו בכל מעגל אחר, על ידי הגדלתה עד שהתוצאה מפסיקה להשתנות.

איך ה-Transpiler מסנתז סיבוב

שלב ה-optimization של הצינור המוגדר מראש מפרק ModeRotation קודם במרחב הבוזוני, עם המעבר B2BSynthesis והשיטה BosonicClementsDecomposition (רשומה כ-ModeRotation.Clements). זו בחירה של ההגדרה המוקדמת ולא דרישה של המסגרת: סעיף מאוחר יותר מוריד את אותו שער בלי שום פירוק. הפירוק לא מעלה בחזקה את המחולל המלא בצעד אחד. הוא מפרק קודם את האוניטרי לכל היותר $M(M-1)/2$ סיבובים בין מצבים סמוכים בתוספת $M$ פאזות חד-מצביות, לפי הרשת המלבנית של Clements ועמיתיו [1]. כל גורם הופך לשער Evolution אטומי, ששלב מצב-לקיוביט ממפה כמו כל אחד אחר.

הפירוק עצמו הוא givens_decomposition(), שאליו הפירוק קורא. ההעתקה מאוניטרי לסיבוב הרב-גופי היא הומומורפיזם של חבורות, ולכן מכפלת סיבובי הגורמים שווה בדיוק לסיבוב המלא במרחב הלא קטוע. כלומר, הרשת לא מוסיפה שגיאה משלה. היתרון שלה הוא לוקליות. כל אבולוציה שנפלטת כוללת רק את איברי Pauli של שני המודים שהיא מערבבת (או של מוד אחד במקרה של פאזה), במקום איברים על פני כל המודים.

בצעו Transpile לסיבוב אקראי של שלושה מודים ובדקו מה ה-pipeline פלט. אוניטרי כללי בגודל $3 \times 3$ דורש לכל היותר שלושה סיבובים, בתוספת שלוש פאזות, וכל גורם הוא PauliEvolutionGate משלו. כולם מוגדרים על כל הרגיסטר, אבל האופרטורים שבתוכם פועלים רק על כמה Qubit:

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

עם שלוש רמות לכל מוד, כל מוד תופס שלושה Qubit. beamsplitter בין שני מודים נוגע בשישה ה-Qubit של אותם מודים, ופאזה נוגעת ב-Qubit של מוד אחד.

הגורמים מדויקים רק כאשר כל אבולוציה מסונתזת היטב. האיברים של היוצר של beamsplitter אינם מתחלפים, ולכן נוסחת המכפלה ברירת המחדל בצעד יחיד משאירה שגיאה בכל גורם. העבירו נוסחת מכפלה עם יותר חזרות ל-plugin של 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

השגיאה יורדת בגורם של כארבעה על כל הגדלה פי ארבעה במספר החזרות, כמצופה מנוסחה מסדר ראשון. העלות גדלה באופן פרופורציונלי, כי כל גורם חוזר על עצמו.

הורדת סיבוב לרמה נמוכה בלי לפרק אותו

פירוק במרחב בוזוני הוא דרך אחת להוריד סיבוב לרמה נמוכה. (הספרייה לא כוללת דרך אחרת, אבל השלב של בוזון-ל-Qubit לא תלוי בזה.) B2QSynthesis מחפש plugin לפי שם השער, ו-ModeRotation שנשאר שלם מגיע אליו כמו כל שער אחר. plugin עבורו ממפה את היוצר של 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))

רשמו אותו תחת שם השער ב-pass manager שאינו מפרק אותו:

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

הסיבוב הוא עכשיו אבולוציה אחת על כל תשעת ה-Qubit, בעוד שברשת יש שש, ואף אחת מהן לא פורשת על יותר משישה. שתי הגישות אינן אותו קירוב. הקטיעה מיושמת פעם אחת על היוצר של כל הסיבוב כאן, ולכל גורם בנפרד ברשת (ומכאן האזהרה על הקטיעה למעלה). המסלול הישיר שווה ל-occupation_encoding() של האופרטור בכל קטיעה, ולכן הוא ההתייחסות להשוואת רשת. המחיר הוא לוקליות: האיברים שלו פורשים על כל המודים, ולנוסחת מכפלה יש יותר מהם לחלק.

אף אחת מהן אינה טובה יותר באופן כללי, וזו הסיבה להשאיר את הבחירה פתוחה. כתבו plugin כמו זה כדי להשוות בין הורדות לרמה נמוכה, כדי להשתמש בסיבוב ייעודי לקידוד מסוים, או כדי להעביר את השער לכלי שמבין אותו באופן טבעי.

השלבים הבאים

• Transpile למעגלים בוזוניים מסביר את השלבים ש-pass manager המוגדר מראש מריץ, ואיך להחליף את ה-plugin שמטפל בשער.

• מיפוי ייצוגי אופרטורים מכסה את הקטיעה, את הקידודים, ואת LevelSpec שקובע את הפריסה של כל מוד.

• עבודה עם מעגלים בוזוניים מראה איך להרכיב ולשלב מעגלים בוזוניים.

• תיעוד ה-API מפרט את qiskit_bosons.circuit.library עבור השערים, את qiskit_bosons.linalg עבור הפירוק, ואת qiskit_bosons.transpiler.passes עבור ה-plugins של הסינתזה.

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