trim SQD を活用した配置選択の改善

diagonalize_fermionic_hamiltonian のデフォルトの動作は、サブサンプルした各バッチを対角化し、最良の結果を保持することです。Trim SQD では代わりに、バッチを選別対象の候補として扱います。各バッチを対角化し、最大の重みを持つ配置までトリムし、すべてのバッチの生き残りをマージして一緒に対角化します。このマージされた対角化のエネルギーがそのイテレーションで報告され、そのトリムされた出力が次のイテレーションの種になります。

trim SQD の実装は、クリーブランド・クリニック財団、理化学研究所 (RIKEN)、IBM® による共同研究に基づいています。この機能を研究で使用する場合は、論文も引用してください。引用情報は CITATION.bib ファイルにあり、引用キーは merz2026crossing です。

この例で使用する SQD コードの詳細については、化学ハミルトニアンのチュートリアルをご覧ください。

分子を設定する

# Added by doQumentation — installs the packages this notebook needs if they are missing
import importlib.util

_needed = {"numpy": "numpy", "pyscf": "pyscf", "qiskit_addon_sqd": "qiskit-addon-sqd"}
_missing = [pip for module, pip in _needed.items()
            if importlib.util.find_spec(module) is None]
# One at a time, so a package that fails to install does not block the others
for _pip in _missing:
    %pip install -q {_pip}
if not _missing:
    print("\u2713 All required packages are installed")

import numpy as np
import pyscf
import pyscf.mcscf
from qiskit_addon_sqd.counts import generate_bit_array_uniform
from qiskit_addon_sqd.fermion import diagonalize_fermionic_hamiltonian
from qiskit_addon_sqd.trim import TrimPolicy

# Build N2 molecule
mol = pyscf.gto.Mole()
mol.build(
    atom=[["N", (0, 0, 0)], ["N", (1.0, 0, 0)]],
    basis="6-31g",
    symmetry="Dooh",
)

# Define active space
n_frozen = 2
active_space = range(n_frozen, mol.nao_nr())

# Get molecular integrals
scf = pyscf.scf.RHF(mol).run()
num_orbitals = len(active_space)
n_electrons = int(sum(scf.mo_occ[active_space]))
num_elec_a = (n_electrons + mol.spin) // 2
num_elec_b = (n_electrons - mol.spin) // 2
cas = pyscf.mcscf.CASCI(scf, num_orbitals, (num_elec_a, num_elec_b))
mo = cas.sort_mo(active_space, base=0)
hcore, nuclear_repulsion_energy = cas.get_h1cas(mo)
eri = pyscf.ao2mo.restore(1, cas.get_h2cas(mo), num_orbitals)

# Compute exact energy
exact_energy = cas.run().e_tot

# Create a seed to control randomness throughout this workflow
rng = np.random.default_rng(24)

# Generate random samples
bit_array = generate_bit_array_uniform(
    10_000, num_orbitals * 2, rand_seed=rng
)

WARN: Unable to to identify input symmetry using original axes.
Different symmetry axes will be used.

converged SCF energy = -108.835236570774
CASCI E = -109.046671778080  E(CI) = -32.8155692383188  S^2 = 0.0000000

ループを実行する

trim SQD の1回のイテレーションでは、次の処理を行います。

選別ラウンド。 ビット列をバッチに分割します。各バッチを対角化し、スピンセクターごとに、重みでランク付けした CI 文字列のうち trim_ratio の割合を保持します。
マージラウンド。 生き残った文字列をマージして1回対角化します。これがそのイテレーションで報告される結果であり、そのトリムされた文字列が次のイテレーションの種になります。

ループの sci_solver が両方のラウンドを実行し、デフォルトではこのパッケージの PySCF ソルバーを使用します。

次のように TrimPolicy を渡すことで、trim SQD を有効にします。

result = diagonalize_fermionic_hamiltonian(
    hcore,
    eri,
    bit_array,
    samples_per_batch=300,
    norb=num_orbitals,
    nelec=(num_elec_a, num_elec_b),
    num_batches=5,
    max_iterations=5,
    policy=TrimPolicy(trim_ratio=0.5, max_strings_per_trim=250),
    symmetrize_spin=True,
    max_dim=250,
    seed=rng,
)

sci_state = result.sci_state
print(
    f"Merged subspace: {len(sci_state.ci_strs_a)} x {len(sci_state.ci_strs_b)} CI strings"
)
print(f"Exact energy:     {exact_energy}")
print(f"Estimated energy: {result.energy + nuclear_repulsion_energy}")

Merged subspace: 148 x 148 CI strings
Exact energy:     -109.04667177808024
Estimated energy: -109.034369147063

選別が役立っているかを確認する

選別が最も効果を発揮するのは、サンプルされた配置の数が一度に対角化できる数を大きく上回る場合です。その場合、より大きなプールをランク付けすることで、単一のバッチでは見逃していたであろう配置を本当に見つけられます。逆に、現実的なサイズの部分空間がすでに波動関数の大部分を捉えている場合は、効果が最も小さくなります。

選別は、バッチ同士が互いに異なっていることに依存します。TrimPolicy は、diagonalize_fermionic_hamiltonian の動作を変更して samples_per_batch * num_batches 個のビット列からなる1つのプールを抽出し、それをバッチ間で分割することで、これを管理します。そのため、複数のバッチに現れるビット列はありません。

部分空間を制限する

ここでの部分空間は、スピンアルファとスピンベータの CI 文字列配列のデカルト積によって張られます。

トリムは配置をランク付けしますが、引き継がれるのはスピンごとの2つの文字列リストです。部分空間の重みが最も高い150個の配置を保持すると、通常は約150個の異なるアルファ文字列と150個の異なるベータ文字列が含まれ、その積は150個ではなく約22,500個の配置を張ります。そのため、次のイテレーションに渡される部分空間は、ランク付けされた集合よりも一般にはるかに大きくなり、保持する数に対して2乗で増加します。

トリムと引き継ぎは異なる時点で行われるため、別々に制限されます。トリムは選別ラウンドの各バッチを縮小し、引き継ぎは次のイテレーションに持ち越されるものです。

4つの引数でこれを制限します。少なくとも1つは設定しておくと便利です。

trim_ratio は、スピンセクターごとに、各バッチの文字列のうち選別を通過するものの割合です。値を小さくするとマージされた部分空間は縮小しますが、次のイテレーションがサンプリングで再発見しなければならない配置を捨てることになります。
max_strings_per_trim は、比率にかかわらず、各バッチがセクターごとに提供する文字列の上限を設定します。両方のスピンセクターを1つのリストにマージし、各バッチが提供するものをほぼ2倍にする symmetrize_spin と併用すると便利です。この制限はバッチごとです。マージされた部分空間は、最大で num_batches 倍の文字列を保持でき、バッチ間で共通していた分だけ少なくなりますが、実際にはそれほど多くありません。
carryover_threshold、carryover_ratio、max_carryover は、次のイテレーションの種になる文字列を制限します。これらはバッチごとではなく、マージされた部分空間全体に適用されます。carryover_threshold は StandardPolicy が使用する振幅のカットオフで、デフォルトも同じです。他の2つは数を制限しますが、カットオフだけでは数を制限できません。
max_dim は diagonalize_fermionic_hamiltonian 自体の引数で、マージされたものを含め、ループが対角化するすべての部分空間の各スピンセクターを制限します。これは安全策として機能し、引き継ぎが新しいサンプルとマージされた後で、ポリシーが返すものに適用されます。

上記の実行では max_strings_per_trim=250 と max_dim=250 を併せて設定しているため、各バッチが提供するものも、それから構築される部分空間も、際限なく増えることはありません。

外部の固有値ソルバーを使用する

sci_solver 引数を渡すことで、外部の固有値ソルバーと一緒に trim SQD を使用できます。

from sbd.sbd_solver import create_sbd_solver

result = diagonalize_fermionic_hamiltonian(
    hcore,
    eri,
    bit_array,
    samples_per_batch=300,
    norb=num_orbitals,
    nelec=(num_elec_a, num_elec_b),
    num_batches=5,
    policy=TrimPolicy(trim_ratio=0.1),
    sci_solver=create_sbd_solver(sbd_config={"method": 0, "eps": 1e-10}),
    seed=rng,
)

前のスニペットを実行するには、PyPI から sbd-eigensolver パッケージをインストールする必要があります。

作業を分散する

TrimPolicy を使うと、ループは1回のイテレーションで sci_solver を2回呼び出します。1回はすべてのバッチを渡して、もう1回はマージされた部分空間を渡して呼び出します。どちらの呼び出しも通常の sci_solver インターフェースを通り、このパッケージのマルチプロセス対応が、すべての MPI ランクで集合的に呼び出します。

利用できる並列性は潜在的に2つのレベルがあり、互いに独立しています。

1回の対角化の内部。 SBD などのソルバーは、1回の対角化を複数のランクに分散できます。
バッチ間。 ソルバーは原理的には、複数のバッチを同時に対角化できます。

バッチの集合が与えられると、ソルバーがワーカープロセスへの作業の委任方法を決定します。

分子を設定する​

ループを実行する​

選別が役立っているかを確認する​

部分空間を制限する​

外部の固有値ソルバーを使用する​

作業を分散する​