ボソン空間でのEvolutionゲートの合成

非可換な部分の和であるハミルトニアンの時間発展 $e^{-itH}$ は、ゲートごとに厳密に実装することはできません。積公式は、個々の部分による短い発展の列でこれを近似し、その近似誤差をトロッター誤差と呼びます。

積公式は2か所に適用でき、両者は交換可能ではありません。

1. ボソン空間で、トランケーションやエンコーディングが選ばれる前に適用します。公式は演算子のグループ(演算子の項のグループ化:演算子の構造を使うを参照)を認識し、それらを交互に配置できます。

2. 量子ビット空間で、演算子がマッピングされた後に、QiskitのPauliEvolutionGate内部で適用します。公式はグループごとに1つのパウリ演算子を認識し、その内部の項しか細分できません。

非可換なグループ間の誤差を減らせるのは最初の方法だけであり、これは通常、支配的な誤差です。このガイドでは、その使い方を示します。

重要

どちらのステップも任意です。ボソン空間で分解されないEvolutionゲートは、そのままボソンから量子ビットへの段階に到達し、そこでマッピングおよび合成されます。この場合、ボソンのレベルでの誤差はまったく生じず、量子ビットレベルの公式が導入する誤差だけが残ります。

セットアップ

2サイトのボース・ハバード・チェーンは、グループが可換でない最小の例です(完全なモデルについてはエンコーディングを選択してボース・ハバード動力学をシミュレートするを参照)。これは、ボンド用のエルミートなグループ1つと、サイト内相互作用用の対角なグループ1つにグループ化されています。

>>> 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ゲートは1次の公式であるBosonicLieTrotterで分解されます。グループごとに1つの因子があり、それぞれが全時間分です。

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

分解しない場合、ゲート全体が量子ビット段階に渡されるため、ボソンのトロッター誤差はまったく生じません。

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

分解すると、2つのグループが可換でないため、実際の誤差が現れます。

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

次数が2増えるたびに、1つ下の次数のスイープが5回合成されます(次数2の $2n - 1$ 個の因子に適用される、鈴木の再帰の $5^{(k-2)/2}$ の増加)。これが、星印のついた最終行の数です。スイープは同じ因子で接するため、隣接するゲートは1つの演算子のもとで発展し、時間を合計した単一のゲートに融合されます。これにより継ぎ目ごとに1ゲート節約でき、しかも厳密です。教科書どおりの列を保持するには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

繰り返し回数は、深さと精度を異なる比率で交換します。1次の誤差はステップ幅に対して線形に減少し、2次の誤差は2乗で減少します。

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

どちらがより良い選択かは演算子によります。グループが1つだけの演算子は、単独の因子が自分自身と可換であるため、どの次数でも厳密であり、その場合、公式は再分割されずに1つのゲートにまとめられることに注意してください。

コストを測定する

誤差を小さくすると因子が増え、各因子は量子ビット段階で合成されると、もつれゲートを消費します。その数はエンコーディングに依存します。そこで、基底を固定し量子ビットレベルの最適化なし(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がEvolutionゲートを合成する方法に依存し、それはリリース間で変わるため、表ではなくプロットで示します。

各棒のグループは1つの方法を表します。リー・トロッター、次数2の鈴木、次数4の鈴木です。次数を上げると精度が得られますが、そのコストは再帰の係数ずつ増加します。また、ユナリー・エンコーディングは、どの方法でも、各因子がより多くのもつれゲートにマッピングされるため、グレイコードより多く支払います。必要な精度に応じて方法を選び、ハードウェアに応じてエンコーディングを選んでください。この2つの選択は切り離しておいてください。上のボソン回路は変わりません。

1つの方法を回路全体に適用する

各ゲートに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を使うと、恩恵を受けるゲートだけに限定できます。数演算子だけで構成されたハミルトニアンのように、すべてのグループが可換なEvolutionは、どの次数でも厳密に合成されるため、次数を上げてもそこでは深さが増えるだけです。

恩恵を受けるゲートにパスを限定する

純粋に対角な演算子は、数演算子どうしが可換なので、この点を具体的に示します。そのグループは1次で厳密であり、2次ではもつれゲートを無駄に消費します。

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

どちらの回路も厳密なので、2つ目はもつれゲートを余分に消費するだけで利得がありません。(対称な2次の公式は最初のグループをそれぞれ半分の時間で2回適用するため、どのグループが先に来るかによってコストが変わります。)filterはDAGOpNodeに対する述語を受け取り、拒否されたゲートには手を加えません。そのため、1つのパスで、恩恵を受けるゲートでは次数を上げ、残りは1次のままにできます。以下の述語は、演算子がモード間でボソンを移動させる場合にのみゲートを受け入れます。これはグループが可換でなくなる場合です。

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

パスをもう一度実行しても何も変わりません。出力される因子はアトミックだからです。

アトミックな因子

atomicとマークされたEvolutionには定義がないため、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で拒否します。偶数次の対称化は非エルミートな因子の誤差を部分的に打ち消すため、誤ってグループ化された演算子は、依然として誤りであるにもかかわらず、次数1より次数2のほうが良く見えることがあります。

短いステップを繰り返す

回路レベルの代替手段として、短いトロッター・ステップを1つ構築して繰り返します。repeat()は繰り返しを1つの回路に平坦化し、境界を隣接したままにするため、パスはそれらをまたいで作用できます。

>>> 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を渡します。エンコーディングは1モードにつき複数の量子ビットを必要とし、モードには量子ビットごとの像がないため、バリアは、それが触れるレジスターのすべての量子ビットにマッピングされます。

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

注記

repeatは合成方法のrepsとは無関係です。合成は単一のEvolutionを、より多くのより短い因子に細分します。1つのEvolutionを含む回路を繰り返しても、その総時間が増えるだけです。

量子ビットレベルの公式と組み合わせる

2つのレベルは合成でき、両方が近似である場合、結果の精度は2つのうち弱いほうで決まります。量子ビット側の公式が1次のままでボソンの次数を上げても、得られるものは少なくなります。量子ビット側の公式についてはMapperFnEvolutionSynthesisを参照してください。これはマッパーが選んだパウリ項の順序を保持するため、group_wise()は各グループの項を互いに隣接させて出力でき、simplify()は代わりに正準な順序に固定できます。

次のステップ

• ボソン回路のトランスパイルでは、ボソン合成に続く段階と、B2BSynthesisがその中のどこに位置するかを説明します。

• 演算子の項のグループ化:演算子の構造を使うでは、積公式が交互に配置するグループの割り当て方と確認方法を説明します。

• エンコーディングを選択してボース・ハバード動力学をシミュレートするでは、繰り返しの方式を完全な例に適用します。