フェルミオン・量子ビット変換の最適化

Posted at 2025-12-24

$$
\def\bra#1{\mathinner{\left\langle{#1}\right|}}
\def\ket#1{\mathinner{\left|{#1}\right\rangle}}
\def\fbra#1{\mathinner{\left({#1}\right|}}
\def\fket#1{\mathinner{\left|{#1}\right)}}
\def\braket#1#2{\mathinner{\left\langle{#1}\middle|#2\right\rangle}}
$$

はじめに

連載ブログ記事ZX-calculusを用いたフェルミオン量子ビット変換（１）（２）（３）（４）において、フェルミオン・量子ビット変換として従来からよく知られた「Jordan-Wigner変換」「Parity変換」「Bravyi-Kitaev変換」以外にも、より汎用的にその変換を記述できる手法としてTernary Treeを用いた変換法があるということを学びました。そして、その変換を実行・最適化するためのPythonライブラリ¹を自作して以下githubに公開しました、ということも追記していました。

yaf2q - Yet another ferimion to qubit mapping library

yaf2qは以下の機能を提供しています。

フェルミオン・量子ビット変換の実行
フェルミオン・量子ビット変換の方法として従来からよく知られている「Jordan-Wigner変換」「Parity変換」「Bravyi-Kitaev変換」に対応するとともに、より汎用的にフェルミオン・量子ビット変換を記述できる「Ternary Tree（三分木）を用いた変換」にも対応しています。
フェルミオン・量子ビット変換の最適化
Ternary Treeを用いてフェルミオン・量子ビット変換を実行したときに、量子ビット演算子に関する何らかの指標（例えば、ハミルトニアンの各項を構成するパウリ重みや、量子位相推定などの量子アルゴリズムを実行する量子回路の深さなど）が最小の値をとるような最適なTernary Treeの構造を探索することができます。

せっかく作ったので、このライブラリを使ってどんなことができるのか、本ブログ記事で簡単に紹介したいと思います²。ちなみに、「フェルミオン・量子ビット変換の実行」については、基本的には参考にさせていただいた論文"From fermions to Qubits: A ZX-Calculus Perspective"に記載されているアルゴリズムを実装しており³、「フェルミオン・量子ビット変換の最適化」については、とりあえず簡単に実装できそうなシミュレーテッド・アニーリングを用いた方法で最適化機能を実装しています⁴。

それでは、はじめます！まずは、フェルミオン・量子ビット変換の実行です⁵

フェルミオン・量子ビット変換の実行

従来の変換法：Jordan-Wigner／Parity／Bravyi-Kitaev変換

フェルミオン・量子ビット変換を管理するF2QMapperクラスを以下のようにimportします。

from yaf2q.f2q_mapper import F2QMapper

従来の「Jordan-Wigner変換」「Parity変換」「Bravyi-Kitaev変換」を実行するためのF2QMapperのインスタンスは、以下のように作成することができます。コンストラクタの引数として変換手法の文字列と量子ビット数を指定します。

f2q_mapper = F2QMapper(kind = "jordan-wigner", num_qubits = 4)
f2q_mapper = F2QMapper(kind = "parity", num_qubits = 4)
f2q_mapper = F2QMapper(kind = "bravyi-kitaev", num_qubits = 4)

f2q_mapperをprintすると、以下のようになります（bravyi-kitaevを指定した場合）。yaf2qでは従来の3手法を"named mapper"と呼んでいます（一般的な用語ではありません）。

print(f2q_mapper)

出力:

named mapper - bravyi-kitaev (num_qubits:4)

フェルミオン・量子ビット変換の変換行列は、encoding_matrixプロパティとして格納されています。

print(f2q_mapper.encoding_matrix)

出力:

[[1 0 0 0]
 [1 1 0 0]
 [0 0 1 0]
 [1 1 1 1]]

この変換行列は、フォック状態を量子ビット状態に変換するためのものです。例えば、$\fket{1100}$というフォック状態は[1,1,0,0]というリストとして表現され、この変換行列を適用することで[1,0,0,0]というリストを得ることができ、これは量子ビット状態$\ket{1000}$を表しています。この行列をフォック状態を表すリストに行列として乗算することで量子ビット状態を得ることはできるのですが、fock_to_qubit_state()というメソッドを使う方が簡単です。以下のようにします。

fock_state = [1, 1, 0, 0]
qubit_state = f2q_mapper.fock_to_qubit_state(fock_state)
print(qubit_state)

出力:

[1, 0, 0, 0]

逆変換のメソッドもあります。

qubit_state = [1, 0, 0, 0]
fock_state = f2q_mapper.qubit_to_fock_state(qubit_state)
print(fock_state)

出力:

[1, 1, 0, 0]

フェルミオン演算子を量子ビット演算子に変換する場合は、まずOpenFermionを使って解析したいフェルミオン演算子を作成することから始めます。例えば、$H_2$分子のフェルミオン・ハミルトニアンは、OpenFermionを使って以下のように作成することができます。仕様の詳細はOpenFermionのドキュメントをご参照ください。

from openfermion import MolecularData
from openfermionpyscf import run_pyscf

molecule = MolecularData(
    geometry = [('H', (0, 0, 0)), ('H', (0, 0, 0.65))]
    basis = "sto-3g",
    multiplicity = 1
    charge = 0,
)
molecule = run_pyscf(molecule, run_scf=1, run_fci=1)
fermion_hamiltonian = molecule.get_molecular_hamiltonian()

このfermion_hamiltonian（フェルミオン演算子）をfermion_to_qubit_operator()メソッドの引数に指定することで、以下のように、量子ビット演算子を得ることができます。

qubit_hamiltonian = f2q_mapper.fermion_to_qubit_operator(fermion_hamiltonian)

このqubit_hamiltonianはyaf2qのQubitOperatorSetクラスのインスタンスになっていて、この中に、OpenFermionで定義されている量子ビット演算子クラス（QubitOperator）やqiskitで定義されている量子ビット演算子クラス（SparsePauliOp）がメンバとして格納されています。各々openfermion_formおよびqiskit_formプロパティとして、以下のように取得することができます⁶。

print(qubit_hamiltonian.openfermion_form)

出力:

0.03775110394645509 [] +
0.04407961290255181 [X0 Z1 X2] +
0.04407961290255181 [X0 Z1 X2 Z3] +
0.04407961290255181 [Y0 Z1 Y2] +
0.04407961290255181 [Y0 Z1 Y2 Z3] +
0.18601648886230604 [Z0] +
0.18601648886230604 [Z0 Z1] +
0.1699209784826151 [Z0 Z1 Z2] +
0.1699209784826151 [Z0 Z1 Z2 Z3] +
0.12584136558006329 [Z0 Z2] +
0.12584136558006329 [Z0 Z2 Z3] +
0.17297610130745106 [Z1] +
-0.26941693141631995 [Z1 Z2 Z3] +
0.17866777775953394 [Z1 Z3] +
-0.26941693141631995 [Z2]

print(qubit_hamiltonian.qiskit_form)

出力:

SparsePauliOp(['IIII', 'IIZI', 'IZII', 'IZIZ', 'IZZZ', 'XZXI', 'XZXZ', 'YZYI', 'YZYZ', 'ZIII', 'ZIZI', 'ZIZZ', 'ZZII', 'ZZZI', 'ZZZZ'], 
	coeffs=[ 0.0377511 +0.j, -0.26941693+0.j,  0.1729761 +0.j,  0.17866778+0.j, -0.26941693+0.j,  0.04407961+0.j,  0.04407961+0.j,  0.04407961+0.j, 0.04407961+0.j,  0.18601649+0.j,  0.12584137+0.j,  0.12584137+0.j, 0.18601649+0.j,  0.16992098+0.j,  0.16992098+0.j])

QubitOperatorSetクラスには、固有値や固有ベクトルを求めるメソッドが定義されています。eigenvalues()メソッドの引数numに指定された数の固有値を小さいものからnum個並んだリストとして得ることができます⁷。

qubit_hamiltonian.eigenvalues(num = 3)

同様にeigenvectors()メソッドの引数numに指定された数の固有ベクトルを固有値の小さいものからnum個並んだリストとして得ることができます。

qubit_hamiltonian.eigenvectors(num = 3)

固有値と固有ベクトルを同時に得たい場合は、eigsh()メソッドを使います。

eigenvalues, eigenvectors = qubit_hamiltonian.eigsh(num = 3)

また、パウリ重み(pauli weight)⁸のリストを得るpauli_weights()メソッドもあります。実行すると、

print(qubit_hamiltonian.pauli_weights())

出力:

[3, 1, 4, 1, 2, 2, 3, 4, 4, 0, 3, 1, 3, 2, 3]

のようなリストが得られます⁹。パウリ重みの平均値や最大値・最小値などを求める元データとして利用することができます。

Ternary Treeを用いた変換法

Ternary Treeは、一般には各ノードが高々3個の子ノードを持つようになっている木構造のことです。が、フェルミオン・量子ビット変換で対象としているのは、以下に示すような特別なTernary Treeです。

ここで、ノードは楕円で表現されておりノード番号がその中に記載されています。ノードから伸びているエッジは3本であり、各々X,Y,Zというラベルがついています。これは各々パウリX,Y,Z演算子に対応しています。一般的なTernary Treeと違うのは子ノードがないエッジも許されているということです。詳細は論文"From fermions to Qubits: A ZX-Calculus Perspective"やそれを解説した連載ブログ記事ZX-calculusを用いたフェルミオン量子ビット変換（１）（２）（３）（４）をご参照いただくとして、ここでは、とりあえず、このようなTernary Treeによって任意のフェルミオン・量子ビット変換を規定することができて、その変換のアルゴリズムもわかっているということをおさえておけば十分です。

さて、それでは、yaf2qでどのようにこのTernary Treeを定義するかを見ていきます。それには、TernaryTreeSpecクラスを使います。そのコンストラクタに、以下のように、indicesとedgesという引数を指定することで、そのインスタンスを生成します。

from yaf2q.ternary_tree_spec import TernaryTreeSpec

ttspec = TernaryTreeSpec(
    indices = [0, 1, 2, 3],
    edges = {1: (0, 'X'), 2: (1, 'Y'), 3: (1, 'Z')},
)

ここで、indicesはノード番号のリストを表しています。4個の要素からなっているので4個のノードからなるTernary Treeということになります。edgesはエッジ集合を表しているのですが、少々説明が必要です。上の例では1,2,3をキーとして、各々に対応した値が(0, 'X'),(1, 'Y'),(1, 'Z')というタプルである、辞書データとして表現されています。キーは子ノード番号を表しています。対応したタプルの1番目の要素はそれが接続している親ノード番号で、2番目の要素は、親ノードに接続するエッジのラベル文字列です。ルートノード番号は0と決められている前提です。なので、子ノードを表すキーに0は含まれずそれ以外のすべてのノード番号が含まれます¹⁰。そして、親ノード番号は必ず子ノード番号より小さく、値のタプルに重複があってはならないという前提もあります¹¹。これでTernary Treeの構造は一意に決定できます（とりあえずyaf2qで独自定義したフォーマットです）。

ttspec.show()

とすると、Ternary Tree図を表示してくれます。実は、上に図示したTernary Treeはこれで作成したものです。

Ternary Treeの木構造はこのままにして、ノード番号のみを入れ替えたものを定義したい場合もあるかもしれません。その場合は、indicesの要素の並び順のみを変えます。例えば、

ttspec = TernaryTreeSpec(
    indices = [1, 0, 3, 2],
    edges = {1: (0, 'X'), 2: (1, 'Y'), 3: (1, 'Z')},
)
ttspec.show()

のようにすると、ノード番号だけ異なるTernary Treeを以下のように得ることができます。このとき、edgesは変更しないでください。変更すると木構造の形も変わってしまいます。indicesはedegesに記載されているノード番号がそれぞれを何番目の量子ビットとするかを表しています。

明示的にindicesとedgesを指定せずにランダムにTernary Treeを作成することもできます。以下のように、random()メソッドにノード数を引数に与えることで得ることができます。

ttspec = TernaryTreeSpec.random(4)
print(ttspec)

出力:

indices:[3, 0, 2, 1], edges:{1: (0, 'Z'), 2: (0, 'Y'), 3: (2, 'Y')

このようにTernary Treeが作成できたところで、これに基づいたフェルミオン・量子ビット変換を実行するため、以下のようにF2QMapperインスタンスを作成します。

f2q_mapper = F2QMapper(ttspec=ttspec)

これ以降は、「従来の変換法：Jordan-Wigner／Parity／Bravyi-Kitaev変換」で述べたのと同様にして、量子ビット状態を取得したり、量子ビット演算子やそれに対応した固有値・固有ベクトルやパウリ重みのリストを得ることができます。

フェルミオン・量子ビット変換の最適化

次に、フェルミオン・量子ビット変換の最適化について説明します。

いま、4個の水素原子からなる1次元水素鎖（$H_4$）を考えて、その量子ビット・ハミルトニアンを構成するパウリ重みの平均を最小化するためのフェルミオン・量子ビット変換（＝Ternary Tree）を求めたいとします。そのために、まず、以下のように、OpenFerimionを使ってフェルミオン・ハミルトニアンfermion_hamiltonianを作成することから始めます。

from openfermion import MolecularData
from openfermionpyscf import run_pyscf

from yaf2q.ternary_tree_spec import TernaryTreeSpec
from yaf2q.f2q_mapper import F2QMapper
from yaf2q.optimizer.sa import SAParams, SA

distance = 0.65
molecule = MolecularData(
    geometry = [('H', (0, 0, 0)), ('H', (0, 0, distance)), ('H', (0, 0, 2*distance)), ('H', (0, 0, 3*distance))],
    basis = 'sto-3g',
    multiplicity = 1,
    charge = 0,
)
molecule = run_pyscf(molecule, run_scf=1, run_fci=1)
fermion_hamiltonian = molecule.get_molecular_hamiltonian()
num_qubits = fermion_hamiltonian.one_body_tensor.shape[0]

そして、TernaryTreeSpecのみを引数としてfloatを返す目的関数を以下のように定義します（関数内関数として定義してください）。

def objective_func(ttspec: TernaryTreeSpec):
    f2q_mapper = F2QMapper(ttspec=ttspec)
    qubit_hamiltonian = f2q_mapper.fermion_to_qubit_operator(fermion_hamiltonian)
    weights = qubit_hamiltonian.pauli_weights()
    weight_ave = sum(weights) / len(weights)
    return weight_ave

フェルミオン・量子ビット変換を実行してできた量子ビット演算子qubit_hamiltonianからパウリ重みリストを取得して、その平均値を返すようになっています。この目的関数は、TernaryTreeSpecを引数としてfloatを返すものであれば、どんな関数を定義しても良いです。例えば、qubit_hamiltonianに基づき量子位相推定の量子回路を作って、その回路深さを(floatとして)返すようにしても良いです。その場合、量子位相推定の深さを最小化するフェルミオン・量子ビット変換（=ternary tree）が得られます。

目的関数が定義できたら、これを最小化する解(ternary tree)を求めます。それにはyaf2qのSAクラスを使います。SAはシミュレーテッド・アニーリング(Simulated Annealing)の略です。yaf2qでは、このアルゴリズムを使って目的関数を最小化する最適解(正確には近似解)を求めるようにしています¹²。SAのコンストラクタには、量子ビット数(num_qubit)、目的関数(objective_func)、パラメータ(params)(後述)、詳細表示するか否かを表すブーリアン(verbose)を指定します。

ttspec_opt = SA(
        num_qubits = num_qubits,
        objective_func = objective_func,
        params = SAParams(init_sampling=10, num_steps=30, cooling_factor=1.0),
        verbose = False,
).optimize()

ここで、paramsはシミュレーテッド・アニーリングを制御するパラメータSAParamsのインスタンスなのですが、これは初期サンプリング(init_sampling)、アニーリングステップ数(num_steps)、冷却因子(cooling_factor)という属性値をもっており、それらをコンストラクタで指定できるようになっています。init_samplingは初期温度を決定するためにランダムサンプリングを最初に実行するのですが、そのサンプル数です。num_stepsはアニーリングの温度を徐々に段階的に下げていく、その段階数です。cooling_factorは温度を下げていくカーブをどれだけ急峻にするかを表す因子です。0.0より大きい値を指定します。大きくなるほど急激に温度を下げる形になります。各々明示的に指定しない場合、init_sampling=10,num_steps=10,cooling_factor=1.0がデフォルト値として設定されます。どんな目的関数を定義するかによって探索の挙動は変わるので、デフォルト値で収束が悪いと感じた場合、各パラメータ値を調整する必要があります。

これで、SAによる最適化器が作成できたので、実際の最適化を実行します。すでに上に示されていますが、optimize()メソッドを使います。これにより(近似的に)最適なternary tree(ttspec_opt)が得られます。

そして、このttspec_optを使って、改めてフェルミオン・量子ビット変換を実行して、量子ビット演算子とパウリ重みを、以下のように計算して表示します。

print(f"* ternary tree:\n{ttspec_opt}")
f2q_mapper = F2QMapper(ttspec=ttspec_opt)
qubit_hamiltonian = f2q_mapper.fermion_to_qubit_operator(fermion_hamiltonian)
weights = qubit_hamiltonian.pauli_weights()
weight_ave = sum(weights) / len(weights)

print(f"* pauli weight (ave): {weight_ave}")

これを実行してみると、例えば、以下の結果を得ることができます。

出力:

[jordan-wigner]
* pauli wieght (ave): 4.583783783783784
[parity]
* pauli weight (ave): 4.691891891891892
[bravyi-kitaev]
* pauli weight (ave): 4.562162162162162
[ternary tree optimization]
* ternary tree:
indices:[3, 7, 2, 0, 5, 4, 1, 6], edges:{1: (0, 'Y'), 2: (1, 'Y'), 3: (2, 'Z'), 4: (3, 'Z'), 5: (4, 'Y'), 6: (3, 'X'), 7: (4, 'X')
* pauli product length (ave): 4.4324324324324325

従来の3手法と比べ、最適Ternary Treeのパウリ重みが一番小さいという結果になっています¹³。

シミュレーテッド・アニーリングは確率的な手法なので、実行のたびに[ternary tree optimization]の結果は変わります。が、乱数のシードを固定することで、結果を固定することもできます。以下のようにseedを指定します。

ttspec_opt = SA(
        num_qubits = num_qubits,
        objective_func = objective_func,
        params = SAParams(init_sampling=10, num_steps=30, cooling_factor=1.0, seed=123),
        verbose = False,
).optimize()

最適化やってみた！

$H_4$分子で最適化の結果が得られたので、別の分子でもパウリ重みがどの程度最適化されるのかやってみたくなります。というわけで、やってみました。$H_4$,$BeH_2$,$H_2 O$,$CH_4$での結果を以下に示します。前述の通り、シミュレーテッド・アニーリングは確率的な手法なので、実行のたびに結果は変わります。本当は何回も試行して平均をとるべきなのですが、今回は、とりあえず1回実行した結果を示しています。すべての分子において、Ternary Treeを最適化した結果がもっとも良いということがわかりました。

ちなみに、ここで使用したコードは、pauli_weight_optimization.pyに置いてあります。ご参考まで。

おわりに

今回は、量子ビット・ハミルトニアンのパウリ重みをなるべく小さくしたいというモチベーションで最適化計算を実行してみましたが、本文中にも記載した通り、目的関数としては、TernaryTreeSpecを引数にして実数値を返す関数であれば何でも良いです。例えば、量子ビット・ハミルトニアンに基づき構成された量子位相推定の量子回路の深さや2量子ビット演算子の数を目的関数にしても良いですし、あるいは、それを特定の量子チップのトポロジー向けにトランスパイルした回路の深さや2量子ビット演算子の数でも良いです。いろいろ最適化してみて結果がどうなるか見てみたい気もしていますが、まだそこまで実行できていません。ご興味ある人がもしいれば、いろいろお試しいただけるとよろしいかと思います。

ただし、まだまだ改善の余地は沢山あります。例えば、シミュレーテッド・アニーリングの近傍探索の方法はまだ工夫できるかもしれませんし、そもそもシミュレーテッド・アニーリング以外の手法を採用した方が良いのかもしれません。それから、本文中では述べませんでしたが、近傍探索するたびに目的関数の評価を愚直に実行しているため、目的関数の中身が複雑な場合、最適化計算の時間が非常にかかってしまうという問題がありそうです。これに対して、並列処理の可能性や目的関数をいちいち評価しないで済ませるための何らかのサロゲートモデルを構築して目的関数の代わりにするということを検討しても良いかもしれません。

とにかく、まだまだ当たり実験的な初期バージョンといったレベルのものなので、今後少しずつ改良していきたいと思っています。

以上

なのですが、中身の処理の大半はRustで書いています。実はRust初心者です。Rust勉強がてらこのライブラリを作成しました。 ↩
本ブログ記事は、yaf2qの日本語チュートリアルをベースにしてなるべくわかりやすく解説したものです。 ↩
ただし、Ternary TreeをPython APIとして扱うためにどんなデータの持ち方をすれば良いかについては、独自に悩んだ（汗）結果がライブラリの仕様に反映されています。 ↩
「簡単に実装できそう」という理由だけでシミュレーテッド・アニーリングを採用しているので、もっと良い最適化の方法はあるのだろうなーと思っています。 ↩
yaf2qの動作はUbuntu24.04(Wsl2)で確認しています。Ubuntu環境が手元にある人は、README.mdに記載されているようにインストールしてみていただけると、以降で説明する内容を確認することができます。…のはずです。が、不具合などあればお知らせください（汗）。 ↩
いまのところOpenfermion形式とqiskit形式のみしかサポートしてないのですが、その他の形式にもそのうち対応するかもしれません。 ↩
固有値の大きいものから取り出すようなメソッドは用意していません。が、そのうち対応するかもしれません。 ↩
pauli weightを「パウリ重み」ととりあえず訳しましたが、日本語として定着した用語はないように思います。本ブログ記事では「パウリ重み」で通します。で、パウリ重みというのは、パウリ演算子の積の中に含まれるパウリ演算子の数のことです。$X_0 Z_2 Y_3$の場合、パウリ重みは$3$です。 ↩
もし、量子ビット演算子が$0.1 X_0 Z_2 Y_3 + 0.2 Z_1 Z_2$だとすると、pauli_weights()メソッドの返り値は例えば[3,2]のようになります（ただし、リストに含まれるパウリ重みの順番は不定です）。 ↩
従って、edgesの辞書サイズは常に(ノード数-1)で、そのキーには1からノード数-1までの整数値が重複のない形で入ることになります。 ↩
edgesに記載されている整数値は、いま考えているTernary Treeの最上位の親ノードから深さ方向に昇順に並んだ整数値であるというルールです。深さが同じノード同士の大小関係はどうでも良いです（自由に決めて良いです）。 ↩
シミュレーテッド・アニーリングを実行する場合、近傍探索をどのように定義するかが重要です。yaf2qの場合、(1)indicesのランダムに選んだ任意の2つの要素をスワップする操作と(2)edgesのランダムに選んだ子ノードの親ノードを（重複がないように）ランダムに変更する操作の2つを交互に実行して近傍探索を行うようにしています。詳細はソースコードをご参照ください。 ↩
従来の3手法に関するコードは省略しました。pauli_weight_optimization.pyをご覧ください。 ↩

You get articles that match your needs
You can efficiently read back useful information
You can use dark theme

What you can do with signing up