סינתזה של שער Evolution במרחב בוזוני

את ההתפתחות בזמן $e^{-itH}$ של המילטוניאן שהוא סכום של חלקים שאינם מתחלפים אי אפשר לממש במדויק שער אחרי שער. נוסחת מכפלה מקרבת אותה באמצעות סדרה של התפתחויות קצרות תחת החלקים הבודדים, ושגיאת הקירוב נקראת שגיאת טרוטר.

אפשר להחיל נוסחת מכפלה בשני מקומות, והם אינם ניתנים להחלפה:

1. במרחב הבוזוני, לפני שנבחרו קטימה או קידוד. הנוסחה רואה את הקבוצות של האופרטור (ראו קיבוץ איברי אופרטור: שימוש במבנה האופרטור) ויכולה לשלב ביניהן.

2. במרחב הקיוביטים, בתוך PauliEvolutionGate של Qiskit, אחרי שהאופרטור מופה. הנוסחה רואה אופרטור פאולי אחד לכל קבוצה ויכולה לחלק רק את האיברים בתוכה.

רק הראשונה יכולה להקטין את השגיאה בין קבוצות שאינן מתחלפות, שהיא בדרך כלל הדומיננטית. המדריך הזה מראה איך להשתמש בה.

חשוב

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

הכנה

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

>>> import numpy as np
>>> from qiskit.quantum_info import SparsePauliOp, Statevector
>>> from qiskit.transpiler import PassManager
>>> from scipy.sparse import csc_matrix
>>> from scipy.sparse.linalg import expm_multiply
>>> from qiskit_bosons.circuit import BosonicCircuit
>>> from qiskit_bosons.circuit.library import Evolution, InitializeModes
>>> from qiskit_bosons.mappers import LevelSpec
>>> from qiskit_bosons.mappers.library import occupation_encoding
>>> from qiskit_bosons.operators import BosonOperator, ann, cre
>>> from qiskit_bosons.transpiler.presets import generate_preset_occupation_pass_manager
>>>
>>> num_sites, max_occupation, total_time = 2, 2, 0.8
>>> bonds = [(0, 1)]
>>> terms = {}
>>> for left, right in bonds:
...     terms[(cre(left), ann(right))] = -1.0
...     terms[(cre(right), ann(left))] = -1.0
>>> for site in range(num_sites):
...     terms[(cre(site), cre(site), ann(site), ann(site))] = 1.0
>>> hamiltonian = BosonOperator.from_dict(terms)
>>> hamiltonian.groups = [
...     len(bonds)
...     if len({mode for _, mode in term}) == 1
...     else bonds.index(tuple(sorted({mode for _, mode in term})))
...     for term, _ in hamiltonian.iter_terms()
... ]
>>> hamiltonian.num_groups()
2

כדי למדוד שגיאה צריך ייחוס. ההתפתחות המדויקת של ההמילטוניאן כפי שמופה באותה קטימה היא הנכונה, כי הקטימה היא קירוב נפרד:

>>> spec = LevelSpec(max_occupation, "unary_projected")
>>> mapped = occupation_encoding(hamiltonian, spec).simplify()
>>> matrix = csc_matrix(SparsePauliOp.from_sparse_observable(mapped).to_matrix(sparse=True))
>>> pm = generate_preset_occupation_pass_manager(spec)
>>> pm.qubit = PassManager()  # keep each group as one exact PauliEvolutionGate
>>>
>>> def build(**evolution_kwargs):
...     circuit = BosonicCircuit(num_sites)
...     circuit.append(InitializeModes([1, 1]), circuit.modes)
...     circuit.append(
...         Evolution(num_sites, hamiltonian, time=total_time, **evolution_kwargs), circuit.modes
...     )
...     return circuit
>>>
>>> # Under the unary encodings, occupation ``n`` of mode ``m`` sets qubit ``m * d + n``, with
>>> # ``d = max_occupation + 1`` levels per mode. That places |1, 1> on the basis state below.
>>> levels = max_occupation + 1
>>> index = sum(1 << (mode * levels + n) for mode, n in enumerate([1, 1]))
>>> initial = Statevector.from_int(index, 2**(num_sites * levels)).data
>>> exact = expm_multiply(-1j * total_time * matrix, initial)

>>> def error(circuit):
...     """Distance to the exact evolution, of the circuit's product of group exponentials."""
...     lowered = pm.run(circuit)
...     state = Statevector.from_int(0, 2**lowered.num_qubits).evolve(lowered).data
...     return float(np.max(np.abs(state - exact)))

הערה

error מדמה בכוונה את המעגל שהורד בלי לפרק את מופעי PauliEvolutionGate שלו. Statevector מחיל כל אחד מהם כאקספוננט המדויק של הקבוצה שלו, כך שמה שנשאר הוא רק מכפלת אקספוננטים של קבוצות שהסינתזה הבוזונית בחרה.

ברירת המחדל

בלי הארגומנט synthesis, שער Evolution מפורק עם BosonicLieTrotter, הנוסחה מהסדר הראשון: גורם אחד לכל קבוצה, כל אחד לכל משך הזמן.

>>> dict(build().decompose().count_ops())
{'Evolution': 2, 'InitializeModes': 1}

בלי פירוק, השער כולו מועבר לשלב הקיוביטים ואין שגיאת טרוטר בוזונית בכלל:

>>> error(build()) < 1e-6
True

אחרי פירוק מופיעה שגיאה אמיתית, כי שתי הקבוצות אינן מתחלפות:

>>> error(build().decompose()) > 0.1
True

הקטנת השגיאה

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

>>> from qiskit_bosons.circuit.library.synthesis import BosonicLieTrotter, BosonicSuzukiTrotter
>>>
>>> for label, synthesis in [
...     ("Lie-Trotter", BosonicLieTrotter()),
...     ("Suzuki, order 2", BosonicSuzukiTrotter(order=2)),
...     ("Suzuki, order 4", BosonicSuzukiTrotter(order=4)),
...     ("Suzuki, order 4*", BosonicSuzukiTrotter(order=4, merge_adjacent=False)),
... ]:
...     circuit = build(synthesis=synthesis).decompose()
...     print(f"{label:16s} {circuit.count_ops()['Evolution']:3d} factors")
Lie-Trotter        2 factors
Suzuki, order 2    3 factors
Suzuki, order 4   11 factors
Suzuki, order 4*  15 factors

כל הגדלה של הסדר בשניים מרכיבה חמישה מעברים של הסדר שמתחתיו (הגידול של $5^{(k-2)/2}$ ברקורסיה של סוזוקי, המוחלת על $2n - 1$ הגורמים של סדר 2), וזה הספירה בשורה האחרונה, המסומנת בכוכבית. המעברים נפגשים באותו גורם, ולכן השערים השכנים מתפתחים תחת אופרטור אחד ומתמזגים לשער יחיד עם הזמן המסוכם. זה חוסך שער אחד בכל תפר והוא מדויק. העבירו merge_adjacent=False כדי לשמור על הרצף מספרי הלימוד. בכל מקרה, סדר גבוה יותר קונה שגיאה קטנה יותר:

>>> lie = error(build(synthesis=BosonicLieTrotter()).decompose())
>>> second = error(build(synthesis=BosonicSuzukiTrotter(order=2)).decompose())
>>> fourth = error(build(synthesis=BosonicSuzukiTrotter(order=4)).decompose())
>>> lie > second > fourth
True

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

>>> lie_errors = [error(build(synthesis=BosonicLieTrotter(reps=r)).decompose()) for r in (1, 2, 4, 8)]
>>> bool(all(a > b for a, b in zip(lie_errors, lie_errors[1:])))
True
>>> bool(lie_errors[0] / lie_errors[-1] > 6)  # about 8, i.e. linear in 1 / reps
True
>>> suzuki_errors = [
...     error(build(synthesis=BosonicSuzukiTrotter(order=2, reps=r)).decompose()) for r in (1, 2, 4)
... ]
>>> bool(suzuki_errors[0] / suzuki_errors[-1] > 10)  # about 16, i.e. quadratic in 1 / reps
True

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

מדידת העלות

שגיאה קטנה יותר עולה ביותר גורמים, וכל גורם עולה בשערי הסתבכות לאחר ששלב הקיוביטים סינתז אותו. כמה תלוי בקידוד, ולכן ספרו אותם דרך צינור ההגדרה המוכנה, עם בסיס קבוע וללא אופטימיזציה ברמת הקיוביט (optimization_level=0), כדי לראות את המחיר של כל שיטה:

>>> basis = ["rz", "sx", "x", "cx"]
>>>
>>> def cx_count(circuit, encoding):
...     spec = LevelSpec(max_occupation, encoding)
...     lowered = generate_preset_occupation_pass_manager(
...         spec, optimization_level=0, basis_gates=basis
...     ).run(circuit)
...     return lowered.count_ops().get("cx", 0)
>>>
>>> methods = [
...     ("Lie-Trotter", BosonicLieTrotter()),
...     ("Suzuki, order 2", BosonicSuzukiTrotter(order=2)),
...     ("Suzuki, order 4", BosonicSuzukiTrotter(order=4)),
... ]
>>> encodings = ("gray", "unary_unprojected")
>>> counts = {
...     encoding: [
...         cx_count(build(synthesis=method).decompose(), encoding) for _, method in methods
...     ]
...     for encoding in encodings
... }
>>> # the order buys accuracy at a price that grows with it
>>> bool(counts["gray"][0] < counts["gray"][1] < counts["gray"][2])
True

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

כל קבוצת עמודות היא שיטה: Lie-Trotter, סוזוקי מסדר 2 וסוזוקי מסדר 4. הסדר קונה דיוק במחיר שגדל בגורם הרקורסיה, והקידוד האונארי משלם יותר מקוד גריי בכל שיטה כי כל אחד מהגורמים שלו ממופה ליותר שערי הסתבכות. בחרו את השיטה לפי הדיוק שאתם צריכים ואת הקידוד לפי החומרה, והפרידו בין שתי הבחירות: המעגל הבוזוני שלמעלה אף פעם לא משתנה.

החלת שיטה אחת על מעגל שלם

כדי להגדיר synthesis על כל שער, צריך להעביר את הבחירה דרך כל מה שבונה שער. מעבר B2BSynthesis עושה את הבחירה פעם אחת לצינור שלם. ה-config שלו נותן שם לשיטה לכל סוג שער, והוא גם בוחר את השיטה וגם מרחיב את השערים:

>>> from qiskit_bosons.transpiler import BosonicPassManager
>>> from qiskit_bosons.transpiler.passes import B2BSynthesis
>>>
>>> trotterize = BosonicPassManager([B2BSynthesis({"Evolution": BosonicSuzukiTrotter(order=2)})])
>>> dict(trotterize.run(build()).count_ops())
{'Evolution': 3, 'InitializeModes': 1}

אפשר גם לתת שם לשיטה, עם הארגומנטים ל-__init__ שלה, במקום לבנות אותה: B2BSynthesis({"Evolution": ("SuzukiTrotter", (), {"order": 2})}). השמות הזמינים הם אלה של B2BSynthesisPluginManager, ושיטה שאינה רשומה בשום מקום אפשר למסור כמופע, כמו למעלה.

עם apply=False המעבר רק בוחר את השיטה ומשאיר את ההרחבה לנקודה שבה רוצים אותה, למשל כדי לבדוק קודם את המעגל. שער אטומי אף פעם לא נוגעים בו, אפילו לא מתייגים אותו:

>>> selected = BosonicPassManager(
...     [B2BSynthesis({"Evolution": ("SuzukiTrotter", (), {"order": 2})}, apply=False)]
... ).run(build())
>>> selected.count_ops()["Evolution"]  # selected, not yet expanded
1
>>> selected.decompose().count_ops()["Evolution"]  # expanded on demand
3

שער שלא הורחב מגיע שלם לשלב מבוזון לקיוביט, שם הוא ממופה בלי שה-synthesis שלו נבדק, ולכן כל השיטות מפיקות אז פלט זהה.

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

הגבלת המעבר לשערים שנהנים מכך

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

>>> # one group per term, each Hermitian; the order is fixed because order two evolves the first
>>> # group twice, so which group comes first changes the cost
>>> diagonal = BosonOperator.from_terms_with_groups(
...     [([cre(0), ann(0), cre(1), ann(1)], 1.0, 0), ([cre(0), ann(0)], 0.5, 1)]
... )
>>>
>>> def diagonal_cost(synthesis):
...     circuit = BosonicCircuit(num_sites)
...     circuit.append(Evolution(num_sites, diagonal, time=total_time, synthesis=synthesis), circuit.modes)
...     return cx_count(circuit.decompose(), "gray")
>>>
>>> bool(diagonal_cost(BosonicSuzukiTrotter(order=2)) > diagonal_cost(BosonicLieTrotter()))
True

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

>>> def is_diagonal(operator):
...     """Whether every term only counts bosons, rather than moving them between modes."""
...     return all(
...         sorted(mode for created, mode in actions if created)
...         == sorted(mode for created, mode in actions if not created)
...         for actions, _ in operator.iter_terms()
...     )
>>>
>>> def mixed_circuit():
...     circuit = BosonicCircuit(num_sites)
...     circuit.append(Evolution(num_sites, hamiltonian, time=total_time), circuit.modes)
...     circuit.append(Evolution(num_sites, diagonal, time=total_time), circuit.modes)
...     return circuit
>>>
>>> order_2 = BosonicSuzukiTrotter(order=2)
>>> everywhere = cx_count(
...     BosonicPassManager([B2BSynthesis({"Evolution": order_2})]).run(mixed_circuit()), "gray"
... )
>>> restricted = cx_count(
...     BosonicPassManager(
...         [
...             B2BSynthesis(
...                 {"Evolution": order_2},
...                 filter=lambda node: not is_diagonal(node.op.operator),
...             )
...         ]
...     ).run(mixed_circuit()),
...     "gray",
... )
>>> bool(everywhere > restricted)  # the filter saves entangling gates
True

אם מריצים את המעבר שוב, שום דבר לא משתנה: הגורמים שהוא פולט הם אטומיים.

גורמים אטומיים

שער Evolution המסומן atomic אין לו הגדרה, ולכן decompose() משאיר אותו במקומו לא משנה כמה פעמים מחילים אותו. סינתזה מסמנת כך את הגורמים שהיא פולטת, וכך עושה גם BosonicQDrift עבור האיברים שהוא דוגם.

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

>>> circuit = build()
>>> once = circuit.decompose()
>>> dict(circuit.decompose(reps=5).count_ops()) == dict(once.count_ops())
True

זהירות

כל גורם חייב להיות הרמיטי בפני עצמו. אופרטור הרמיטי שמפוצל לקבוצות יכול להפיק קבוצות שאינן הרמיטיות, כמו איבר דילוג שהופרד מהצמוד שלו. בדקו עם groups_are_hermitian(). אופרטור בלי קבוצות מפוצל איבר אחר איבר, ואז שלב הקיוביטים דוחה את הגורמים עם Operator contains complex coefficients. הסימטריזציה של סדר זוגי מבטלת חלקית את השגיאה של גורם שאינו הרמיטי, ולכן אופרטור מקובץ לא נכון יכול להיראות טוב יותר בסדר 2 מאשר בסדר 1 ועדיין להיות שגוי.

חזרה על צעד קצר

כחלופה ברמת המעגל, בנו צעד טרוטר קצר אחד וחזרו עליו. repeat() משטח את החזרות למעגל אחד ומשאיר את הגבולות צמודים, כך שמעברים יכולים לפעול מעבר להם:

>>> step = BosonicCircuit(num_sites)
>>> step.append(Evolution(num_sites, hamiltonian, time=total_time / 4), step.modes)
>>> dict(step.repeat(4).count_ops())
{'Evolution': 4}

העבירו insert_barriers=True כדי להפריד בין החזרות. מחסום ממופה לכל קיוביט באוגרים שהוא נוגע בהם, כי קידוד צריך כמה קיוביטים לכל אופן, ולאופן אין תמונה משלו לכל קיוביט:

>>> [inst.operation.name for inst in step.repeat(3, insert_barriers=True)._inner.data]
['Evolution', 'barrier', 'Evolution', 'barrier', 'Evolution']

הערה

repeat אינו קשור ל-reps של שיטת סינתזה. סינתזה מחלקת התפתחות אחת ליותר גורמים קצרים יותר. חזרה על מעגל שמכיל התפתחות אחת רק מכפילה את הזמן הכולל שלה.

שילוב עם נוסחה ברמת הקיוביט

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

הצעדים הבאים

• Transpile של מעגלים בוזוניים מסביר את השלבים שבאים אחרי הסינתזה הבוזונית, ואיפה B2BSynthesis ממוקם ביניהם.

• קיבוץ איברי אופרטור: שימוש במבנה האופרטור מסביר איך להקצות ולבדוק את הקבוצות ששילוב נוסחת מכפלה משלב ביניהן.

• סימולציה של דינמיקת בוז-האברד עם בחירת קידוד מחיל את תכנית החזרות על דוגמה שלמה.