演算子表現の設計原則

このガイドでは、BosonOperator がどのように構成されているか、またその扱い方を説明します。演算子の作成方法、演算子の調べ方と結合の仕方、そして2つの演算子を比較できる形にする方法です。ハミルトニアンを書き下すときは、ここから始めてください。

BosonOperator は、複素係数を持つボソンの生成演算子と消滅演算子の積の厳密な和です。占有数のカットオフは持ちません。ボソンモードは任意の数のボソンを保持できるため、そのフォック空間は無限次元であり、切り捨ては演算子が量子ビットにマッピングされるときにのみ導入されます（フォック空間を切り捨てるを参照）。したがって、このガイドのすべては、正準交換関係の厳密な代数の中で行われます。

$$[b_i, b_j^\dagger] = \delta_{ij}, \qquad [b_i, b_j] = [b_i^\dagger, b_j^\dagger] = 0 \, .$$

2つの帰結がインターフェースを形作っています。$b_i^\dagger b_i^\dagger$ は2つのボソンを生成する有効な演算子なので、モードのインデックスは1つの項の中で繰り返すことができます。また、2つの因子を並べ替えても符号は生じませんが、消滅演算子を同じモードの生成演算子の向こう側へ移動すると、$b_i b_i^\dagger = b_i^\dagger b_i + 1$ のように項が1つ加わります。

ストレージ

演算子は、疎行列の圧縮スパース行形式に似た、4つのフラットな配列に項を保持します。

coeffs	項ごとに1つの複素係数。
actions	因子ごとに1つのブール値。生成演算子なら True、消滅演算子なら False。
modes	因子ごとに1つのモードインデックス。
boundaries	actions と modes の中で、各項が始まる位置と終わる位置。

コンストラクターはこれら4つの配列を直接受け取ります。以下の演算子には、$2 b_0^\dagger b_1$、$-i b_1^\dagger b_1^\dagger$、定数 $0.5$ の3つの項があり、get_* メソッドで配列を読み戻すことができます。

>>> from qiskit_bosons.operators import BosonOperator, ann, cre
>>> op = BosonOperator(
...     coeffs=[2.0, -1.0j, 0.5],
...     actions=[True, False, True, True],
...     modes=[0, 1, 1, 1],
...     boundaries=[0, 2, 4, 4],
... )
>>> op.get_coeffs()
[(2+0j), (-0-1j), (0.5+0j)]
>>> op.get_actions()
[True, False, True, True]
>>> op.get_modes()
[0, 1, 1, 1]
>>> op.get_boundaries()
[0, 2, 4, 4]
>>> len(op)
3

最初の項は actions と modes の位置0から2にわたり、2番目は2から4、3番目は4から4にわたります。これが、因子を持たない項である定数項の保存方法です。2番目の項で繰り返されているモード 1 は、書かれたとおりに保持されます。モードインデックスが2回現れるからといって、項が削除されたり結合されたりすることはありません。

注

from_dict() で作成した演算子は、辞書の順序で項を保存しないため、実行ごとに配列内の項の順序が変わることがあります。演算子自体は同じであり、以下の演算子を比較したり出力したりする処理は、保存された順序に依存しません。

このレイアウトにより、構築と算術演算が低コストになりますが、その結果として、演算子が同じ項を同時に複数持つことがあります。和と積は項を追加するだけで結合しないため、一意な項が必要な場合は simplify() を呼び出してください（算術演算を行うを参照）。

演算子を作成する

演算子を作成する最も直接的な方法は from_dict() で、因子のタプルから係数への辞書を受け取ります。関数 cre() と ann() が因子を作成し、空のタプルは定数項です。

>>> b = BosonOperator.from_dict({(ann(0),): 1.0})
>>> b_dag = BosonOperator.from_dict({(cre(0),): 1.0})
>>> number = BosonOperator.from_dict({(cre(0), ann(0)): 1.0})
>>> print(format(number))
  1.000000e0 +0.000000e0j * (+0 -0)

出力形式では、+0 は $b_0^\dagger$、-0 は $b_0$ を表します。コンストラクター zero() と one() は加法単位元と乗法単位元を返し、from_terms() は iter_terms() のイテレーターから演算子を再構築します。

演算子ライブラリは、無次元の位置と運動量の直交位相成分を提供します。

$$\hat{q} = \frac{b^\dagger + b}{\sqrt{2}}, \qquad \hat{p} = \frac{i \left(b^\dagger - b\right)}{\sqrt{2}} \, ,$$

position() と momentum() です。同じモジュールの交換子関数と組み合わせて、代数を直接確認できます。

>>> from qiskit_bosons.operators.library import commutator, momentum, position
>>> q, p = position(0), momentum(0)
>>> q.is_hermitian() and p.is_hermitian()
True
>>> commutator(q, p).equiv(1j * BosonOperator.one())
True
>>> commutator(b, b_dag).equiv(BosonOperator.one())
True

このモジュールには anti_commutator() と double_commutator() もあります。3つとも同じ型の演算子を受け取り、新しい演算子を返します。

積の順序に注意する

2つの演算子は2通りの方法で掛け合わせることができ、混同しやすいので注意が必要です。@ 演算子は数学的な積なので、a @ b は $AB$ です。& 演算子は合成で、Qiskit の左から右へ読む規約に従います。a & b は最初に a を適用し、次に b を適用するので、$BA$ になります。交換子関数は @ を使うので、commutator(a, b) は $[A, B]$ です。& を使って自分でブラケットを作ると、気づかないうちに符号が反転します。

>>> (b @ b_dag).equiv(b_dag @ b + BosonOperator.one())
True
>>> (b & b_dag).equiv(b_dag @ b)
True
>>> ((b @ b_dag) - (b_dag @ b)).equiv(BosonOperator.one())
True
>>> ((b & b_dag) - (b_dag & b)).equiv(-BosonOperator.one())
True

演算子を調べる

演算子の項は上記のフラットな配列に保持されているため、演算子は直接イテレートできません。代わりに iter_terms() を使ってください。各項について、(is_creation, mode) のペアのリストと係数を返します。イテレーターはコピーを返すため、それらを変更しても演算子は変わりません。

>>> hopping = BosonOperator.from_dict(
...     {
...         (cre(0), ann(1)): -1.0,
...         (cre(1), ann(0)): -1.0,
...         (cre(2), cre(2), ann(2), ann(2)): 0.5,
...     }
... )
>>> len(hopping)
3
>>> sorted(hopping.iter_terms())
[([(True, 0), (False, 1)], (-1+0j)), ([(True, 1), (False, 0)], (-1+0j)), ([(True, 2), (True, 2), (False, 2), (False, 2)], (0.5+0j))]

いくつかのメソッドが演算子の構造を要約します。

>>> hopping.get_support()
{0, 1, 2}
>>> hopping.min_num_modes()
3
>>> hopping.max_rank()
4

get_support() は演算子が作用するモードの集合、min_num_modes() はそれを保持できる最小のレジスタ、max_rank() は最も長い項の長さです。

注

モードは、ボソンの1つの自由度を表す抽象的なインデックスです。それ自体に意味はありません。モードが格子のどのサイト、どの振動座標、どの光学ポートを表すかは、あなたが決めます。BosonicCircuit は同じインデックスを使うため（ボソン回路を扱うを参照）、演算子とそれを運ぶゲートは、どちらもモード 2 が何を表すかを知らなくても、それが何であるかについて一致します。演算子はモードのレベル数を何も固定しないため、レジスタのサイズは下限にすぎません。min_num_modes() は演算子を保持する最小のレジスタであり、演算子の性質ではありません。

占有数に関するメソッドは、1つの項が1つのモードに置けるボソンの数を報告します。これはマッピングに必要なカットオフの下限です。時間発展のように演算子を繰り返し適用すると、より高いレベルが占有される可能性があるため、十分なカットオフではありません（フォック空間を切り捨てるを参照）。

>>> hopping.mode_occupations()
[1, 1, 2]
>>> hopping.mode_occupations([2, 0])
[2, 1]
>>> hopping.max_occupation()
2

さらに2つのチェックで保存量を調べます。すべての項で生成演算子と消滅演算子の数が同じであれば、演算子は粒子数を保存します。ブロックサイズのリストを渡すと、連続するモードの各ブロック内で同じ条件を検査します。

>>> hopping.conserves_particle_number()
True
>>> move = BosonOperator.from_dict({(cre(0), ann(2)): 1.0})
>>> move.conserves_sector([3])
True
>>> move.conserves_sector([2, 1])
False

最後に、relabel_modes() はモードの名前を付け替えます。モード i は permutation[i] になります。リストは演算子のすべてのモードをインデックスできる長さで、エントリーを繰り返してはいけません。

>>> print(format(move.relabel_modes([5, 6, 7])))
  1.000000e0 +0.000000e0j * (+5 -7)

算術演算を行う

演算子は、加算、減算、符号反転、数による乗算と除算をサポートしており、通常の演算子形式でも、インプレースでも使えます。

>>> total = 2 * b_dag + b_dag - 3j * b
>>> len(total)
3
>>> len(total.simplify())
2
>>> total.simplify().equiv(3 * b_dag - 3j * b)
True

total の3つの項は結合されたことがなかったため、simplify() がそれらを結合するまで、演算子は b_dag を2回保持していました。このメソッドは同一の項の係数を合計し、合計の大きさが許容誤差を下回るものを削除します。これは構造的な操作です。書かれた因子を比較するだけで、交換関係は適用しません。関連メソッドの ichop() は、先に結合することなく小さな係数をインプレースで削除します。そのため、演算子が重複した項を持つ可能性がある場合は、ichop() の前に simplify() を呼び出してください。

積には、上で説明したとおり @ または & を使い、** は演算子を整数乗します。

>>> x = b_dag + b
>>> print(format((x ** 2).simplify()))
  1.000000e0 +0.000000e0j * (-0 -0)
  1.000000e0 +0.000000e0j * (-0 +0)
  1.000000e0 +0.000000e0j * (+0 -0)
  1.000000e0 +0.000000e0j * (+0 +0)

結果は、書かれたとおりの4つの順序を列挙します。それらを単一の形にまとめるには交換関係が必要で、次のセクションで適用します。

随伴は、各項の因子の順序を逆にし、生成演算子と消滅演算子を入れ替え、係数を複素共役にします。メソッド is_hermitian() は、演算子がその随伴と等しいかどうかを検査します。

>>> hop = BosonOperator.from_dict({(cre(0), ann(1)): 1j})
>>> print(format(hop.adjoint()))
  0.000000e0 -1.000000e0j * (+1 -0)
>>> hop.is_hermitian()
False
>>> (hop + hop.adjoint()).is_hermitian()
True

単独のホッピング項はエルミートではありません。時間発展のために項をグループ化する場合、各グループはそれ自体でエルミートでなければならないため、この点が重要になります（演算子の項をグループ化する：演算子の構造を使うを参照）。

正規順序と演算子の比較

代数では $b\,b^\dagger$ と $b^\dagger b + 1$ は等しいため、梯子演算子の2つの積は、同じ演算子でありながら見た目が異なることがあります。この意味で等しい演算子が同じように保存されているとは限らないため、比較には交換関係が必要です。メソッド normal_ordered() がそれを適用します。すべての生成演算子をすべての消滅演算子の左に移動し、同じモード上のペアが入れ替わるたびに項を追加します。

>>> product = BosonOperator.from_dict({(ann(0), cre(0), ann(0), cre(0)): 1.0})
>>> ordered = product.normal_ordered()
>>> len(ordered)
5
>>> print(format(ordered.simplify()))
  1.000000e0 +0.000000e0j * ()
  3.000000e0 +0.000000e0j * (+0 -0)
  1.000000e0 +0.000000e0j * (+0 +0 -0 -0)

順序付けされた演算子は結合されていない和として返され、5つの保存された項は、simplify() で結合すると3つに減ります。正規順序によって項が消えることはありません。$b^\dagger b^\dagger$ やその他の繰り返された因子はそのまま残ります。

sandwich 引数は、2つの因子グループ内でモードをどのように並べるかを設定します。デフォルトでは、両方のグループが降順に並べ替えられます。sandwich=True の場合、生成演算子は昇順、消滅演算子は降順になり、大きなインデックスが中央寄りになります。

>>> pair = BosonOperator.from_dict({(cre(0), cre(1), ann(0), ann(1)): 1.0})
>>> list(pair.normal_ordered().iter_terms())
[([(True, 1), (True, 0), (False, 1), (False, 0)], (1+0j))]
>>> list(pair.normal_ordered(sandwich=True).iter_terms())
[([(True, 0), (True, 1), (False, 1), (False, 0)], (1+0j))]

== 演算子は保存された項の完全な等価性を、許容誤差も簡約化もなしで検査するため、数学的に等しい2つの演算子の違いを報告することがあります。メソッド equiv() は、まず2つの演算子の差を正規順序にしてから、すべての係数が許容誤差を下回ることを検査します。ほとんどの場合、求めている比較はこれです。

>>> reordered = BosonOperator.from_dict({(cre(0), ann(0)): 1.0, (): 1.0})
>>> contracted = BosonOperator.from_dict({(ann(0), cre(0)): 1.0})
>>> contracted == reordered
False
>>> contracted.simplify() == reordered.simplify()
False
>>> contracted.equiv(reordered)
True

メソッド simplify だけでは、書かれたとおりの同一の項しか結合しないため、2つの形は分かれたままです。同様に、is_hermitian() は演算子をその随伴と正規順序で比較するため、追加の作業なしに $b\,b^\dagger + b^\dagger b$ をエルミートと認識します。

注

正規順序は切り捨てのない代数では厳密な恒等式であり、equiv はそれに依存しています。切り捨ての後では恒等式ではありません。マッピングされた演算子はマッピングする形に依存するため、上記の contracted と reordered のマッピング像は、保持される最高レベルでのみ異なります。違いの大きさと形の選び方については、フォック空間を切り捨てるを参照してください。

項の順序付けとフィルタリング

演算子の項を並べ替えても（項の内部の因子ではありません）演算子は変わりませんが、ランダム化積公式のように項を順に読み取るアルゴリズムの挙動は変わることがあります。qiskit_bosons.operators.terms の関数は新しい演算子を返し、項に付いているグループインデックスを保持します。

canonical_order() は、書かれた因子のみに依存する固定の順序に項を並べ替えるため、演算子がどのように組み立てられても同じように読めます。order_terms() は、指定したキーで並べ替えます。このキーは、iter_terms() が返すタプルとして各項を受け取ります。

>>> from qiskit_bosons.operators.terms import canonical_order, order_terms
>>> mixed = BosonOperator.from_dict({(cre(1), ann(0)): 1.0, (cre(0), ann(1)): 3.0, (): 2.0})
>>> list(canonical_order(mixed).iter_terms())
[([], (2+0j)), ([(True, 0), (False, 1)], (3+0j)), ([(True, 1), (False, 0)], (1+0j))]
>>> by_weight = order_terms(mixed, key=lambda term: abs(term[1]), reverse=True)
>>> [coeff for _, coeff in by_weight.iter_terms()]
[(3+0j), (2+0j), (1+0j)]

演算子がグループインデックスを持つ場合（演算子の項をグループ化する：演算子の構造を使うを参照）、関数 group_order() は、各グループを項の連続した1つの区間にします。

>>> from qiskit_bosons.operators.terms import group_order
>>> grouped = BosonOperator.from_dict({(cre(1), ann(0)): 1.0, (cre(0), ann(1)): 3.0, (): 2.0})
>>> grouped.groups = [1, 0, 1]
>>> group_order(grouped).groups
[0, 1, 1]

最後のヘルパーは、項を並べ替えるのではなく削除します。数演算子の積である項 $b_i^\dagger b_i$ は、各フォック状態に位相を掛けるだけなので、計算する対象にとってその位相が無関係な場合、filter_diagonal_terms() は定数、数演算子、およびそれらの積を削除できます。演算子をインプレースで変更し、正規順序の演算子を必要としますが、これは検査しません。

>>> from qiskit_bosons.operators.terms import filter_diagonal_terms
>>> model = BosonOperator.from_dict(
...     {
...         (cre(0), ann(1)): -1.0,
...         (cre(1), ann(0)): -1.0,
...         (cre(0), cre(0), ann(0), ann(0)): 0.5,
...         (cre(0), ann(0)): 0.2,
...         (): 0.1,
...     }
... )
>>> filter_diagonal_terms(model)
>>> len(model)
2
>>> model.is_hermitian()
True

2つのホッピング項だけが残ります。サイト内相互作用と数演算子は占有数基底で対角であるため、定数とともに削除されます。相互作用で繰り返されるモードは正しく扱われます。$b_0^\dagger b_0^\dagger b_0 b_0$ は対角ですが、$b_0^\dagger b_0^\dagger b_0$ は対角ではありません。

次のステップ

• 演算子の項をグループ化する：演算子の構造を使うでは、項にグループインデックスを付ける方法を示します。

• 演算子表現をマッピングするでは、エンコーディングの選択や切り捨ての影響を含め、演算子を量子ビット演算子に変換します。

• ボソン回路を扱うでは、演算子を Evolution ゲートとして回路に配置します。

• エンコーディングを選択してボース・ハバード動力学をシミュレートするでは、ボース・ハバード鎖のハミルトニアンを構築し、エンドツーエンドでシミュレートします。

• qiskit_bosons.operators と qiskit_bosons.operators.terms モジュールに、完全なインターフェースが記載されています。