クリロフ時間発展回路(SKQD)の生成

重要

このガイドの概念は、現在Python APIでのみ利用できます。同等の機能は、今後のリリースでC APIでも提供される予定です。

サンプルベース・クリロフ量子対角化(SKQD)は、部分空間対角化の量子中心型の変種です。単一の変分アンザッツをサンプリングする(通常のサンプルベース量子対角化(SQD)のように)、またはランダム化された時間発展回路のアンサンブルをサンプリングする(SqDRIFTのように)のではなく、クリロフ基底、つまり次の状態の族をサンプリングします。

$$\lvert \psi_k \rangle = e^{-i k \, \Delta t \, H} \lvert \psi_0 \rangle, \qquad k = 0, 1, \dots, D - 1,$$

参照状態$\lvert \psi_0 \rangle$を、固定の時間ステップ$\Delta t$の倍数を増やしながら時間発展させて得られます。これら$D$個の回路からサンプリングされたビット列は、その後ハミルトニアンを古典的に対角化するクリロフ部分空間を構成します。

このガイドでは、SKQDの論文に従い、非相互作用の浴に結合した磁気不純物である単一不純物アンダーソン模型(SIAM)について、SKQDアルゴリズムの回路構築ステップを再現します。

注

SKQDのサンプリングステップが機能するのは、基底状態が計算基底で疎である場合に限られます。つまり、古典的な対角化によって、扱える数のサンプリングされたビット列から基底状態を再構成できる必要があります。以下で示すように、この疎性は自動的に得られるものではなく、問題を表現する1粒子基底の性質です。基底を適切に選ぶことが、この構築の核心です。

1. SIAMハミルトニアンをフェルミオン演算子として構築する

単一不純物アンダーソン模型は、1つの相互作用する不純物軌道(インデックス$0$)を、鎖状に配置された$L$個のサイトからなる非相互作用の浴に接触させたものです。そのハミルトニアンは、1体部分$H_1$と2体部分$H_2$に分けられます。

$$H = H_1 + H_2,$$

1体部分は、不純物のオンサイトエネルギー(化学ポテンシャル$\mu$)、不純物と浴の混成$V$、および浴の最近接ホッピング$t$をまとめたものです。

$$H_1 = \mu \sum_\sigma a^\dagger_{0\sigma} a_{0\sigma} - V \sum_\sigma \left( a^\dagger_{0\sigma} a_{1\sigma} + \text{h.c.} \right) - t \sum_{\sigma} \sum_{i=1}^{L-1} \left( a^\dagger_{i\sigma} a_{i+1,\sigma} + \text{h.c.} \right),$$

2体部分は、不純物のみにオンサイトのクーロン反発$U$を持ちます。

$$H_2 = U \, a^\dagger_{0\uparrow} a_{0\uparrow} a^\dagger_{0\downarrow} a_{0\downarrow}.$$

このガイドでは、粒子・正孔対称点$\mu = -U/2$で扱います。これら2つの部分は、1体積分h1e(浴のホッピング、不純物の混成、不純物のオンサイトエネルギー)と、単一の2体積分h2e(不純物のクーロン項)に直接対応します。まず、模型が自然に定義される位置基底でこれらを組み立てます。

>>> import numpy as np
>>>
>>> norb = 8
>>> nelec = (norb // 2, norb // 2)  # half filling
>>> hopping, onsite, hybridization = 1.0, 10.0, 1.0
>>> chemical_potential = -0.5 * onsite
>>>
>>> # one-body integrals: bath hopping, impurity hybridization, and impurity on-site energy
>>> h1e = np.zeros((norb, norb))
>>> np.fill_diagonal(h1e[:, 1:], -hopping)
>>> np.fill_diagonal(h1e[1:, :], -hopping)
>>> h1e[0, 1] = h1e[1, 0] = -hybridization
>>> h1e[0, 0] = chemical_potential
>>>
>>> # two-body integrals: on-site Coulomb repulsion on the impurity orbital
>>> h2e = np.zeros((norb,) * 4)
>>> h2e[0, 0, 0, 0] = onsite

2. 運動量基底への変換

SIAMの浴の部分は自由フェルミオン鎖であるため、運動量基底への変換によって対角化されます。これはSKQDにとって決定的なステップです。浴は並進対称性を持つため、基底状態は位置基底では疎ではありませんが、浴のハミルトニアンが対角になる運動量基底では疎になります。この変換を行う1粒子軌道回転は、問題固有の線形代数です。これは、不純物軌道には手を加えずに浴のホッピング行列を対角化し、その後、不純物を中心のサイトに移動させます。

>>> def momentum_basis(norb):
...     """Orbital rotation from the position to the momentum basis."""
...     n_bath = norb - 1
...     hopping_matrix = np.zeros((n_bath, n_bath))
...     np.fill_diagonal(hopping_matrix[:, 1:], -1)
...     np.fill_diagonal(hopping_matrix[1:, :], -1)
...     _, vecs = np.linalg.eigh(hopping_matrix)
...     orbital_rotation = np.zeros((norb, norb))
...     orbital_rotation[0, 0] = 1.0
...     orbital_rotation[1:, 1:] = vecs
...     new_index = n_bath // 2
...     perm = np.r_[1:(new_index + 1), 0, (new_index + 1):norb]
...     return orbital_rotation[:, perm]
...
>>> orbital_rotation = momentum_basis(norb)

変分ワークフローとは異なり、ここではハミルトニアンの積分を運動量基底に回転させます。そのため、回路全体(ハミルトニアンも状態も)が基底状態の疎な基底で記述され、サンプリングされたビット列が有用な情報を持ちます。1体積分と2体積分の回転は、電子構造テンソルの標準的な基底変換です。

>>> rotation = orbital_rotation.T.conj()
>>> h1e_momentum = np.einsum("ab,Aa,Bb->AB", h1e, rotation, rotation.conj(), optimize="greedy")
>>> h2e_momentum = np.einsum(
...     "abcd,Aa,Bb,Cc,Dd->ABCD",
...     h2e,
...     rotation,
...     rotation.conj(),
...     rotation,
...     rotation.conj(),
...     optimize="greedy",
... )

この回転は不純物軌道を固定したまま(浴にのみ作用するため)、オンサイト相互作用は変化しません。h2e_momentumはh2eと同様に、単一の数・数項のままです。唯一の変化は、不純物が中心のインデックスに移動したことです。これは積分テンソルから直接確認でき、それぞれが非ゼロの要素を1つだけ持ちます。

>>> print(f"h2e nonzero at:          {np.argwhere(np.abs(h2e) > 1e-12).tolist()}")
h2e nonzero at:          [[0, 0, 0, 0]]
>>> print(f"h2e_momentum nonzero at: {np.argwhere(np.abs(h2e_momentum) > 1e-12).tolist()}")
h2e_momentum nonzero at: [[3, 3, 3, 3]]

この局所性が、後で相互作用を安価に合成できる理由です。

1体積分は、下三角順に詰めて格納されていることを前提とする電子積分コンストラクターfrom_1body_tril_spin_sym()を使って、FermionOperatorに組み立てます。2体部分は単一の数・数項なので、2体積分の仕組みを経由せず、cre()/ann()で直接構築します。これは不純物モード$p$に対するオンサイトのクーロン演算子$U \, n_{p\uparrow} n_{p\downarrow}$です(ブロックスピン配置では、$\uparrow$モードはインデックスp、$\downarrow$モードはインデックスp + norb)。どちらのコンストラクターも、norb個の空間軌道をスピン倍化して2 * norb個のフェルミオンモードにします。

1体部分$H_1$と2体部分$H_2$は、すぐに足し合わせるのではなく、別々の演算子として保持します。これは、セクション1の式の項$H_1$と$H_2$に対応します。各部分を(両方の基底で)一度だけ構築しておけば、以降で再利用できます。足し合わせると、どちらの基底でも厳密対角化のための完全なハミルトニアンになります。個別には、$H_2$(運動量基底)が、ステップ4でトロッターステップが時間発展させる演算子です。

>>> from qiskit_fermions.operators import FermionOperator, ann, cre
>>>
>>> # U n_imp_up n_imp_down: the impurity's up/down modes sit at `imp` and `imp + norb` in the
>>> # block-spin layout (the impurity is orbital 0 in the position basis, orbital 3 in the
>>> # momentum basis -- the single relocated index found above)
>>> imp_pos, imp_mom = 0, 3
>>>
>>> one_body_position = FermionOperator.from_1body_tril_spin_sym(
...     h1e[np.tril_indices(norb)], norb
... )
>>> two_body_position = FermionOperator.from_dict(
...     {(cre(imp_pos), ann(imp_pos), cre(imp_pos + norb), ann(imp_pos + norb)): onsite}
... )
>>> one_body_momentum = FermionOperator.from_1body_tril_spin_sym(
...     h1e_momentum[np.tril_indices(norb)], norb
... )
>>> two_body_momentum = FermionOperator.from_dict(
...     {(cre(imp_mom), ann(imp_mom), cre(imp_mom + norb), ann(imp_mom + norb)): onsite}
... )
>>>
>>> hamiltonian_position = (
...     (one_body_position + two_body_position).normal_ordered().simplify(atol=1e-14)
... )
>>> hamiltonian_momentum = (
...     (one_body_momentum + two_body_momentum).normal_ordered().simplify(atol=1e-14)
... )

基底変換はユニタリ相似変換なので、スペクトルは変わりません。$(\text{norb}, \text{nelec})$セクターにおける最低固有値を厳密対角化と比較して、演算子が正しいことを確認します。FermionOperatorはSciPyのLinearOperatorを公開しており(ネイティブな完全配置間相互作用(FCI)の行列・ベクトル積カーネルを基盤とします)、scipy.sparse.linalg.eigsh()にそのまま渡せます。

>>> import ffsim
>>> import scipy.sparse.linalg
>>>
>>> linop = ffsim.linear_operator(hamiltonian_momentum, norb, nelec)
>>> reference_energy, ground_state = scipy.sparse.linalg.eigsh(linop, k=1, which="SA")
>>> reference_energy = reference_energy[0]
>>> ground_state = ground_state[:, 0]
>>> print(f"exact ground-state energy: {reference_energy:.6f}")
exact ground-state energy: -13.422492

基底変換の利点は疎性です。これを具体的に示すために、基底状態の重みの99%を捉えるのに必要な計算基底の行列式の数を数え、運動量基底と位置基底を比較します。エネルギーは基底に依存しませんが、行列式の数は依存します。

>>> def determinants_for(vec, fraction=0.99):
...     weights = np.sort(np.abs(vec) ** 2)[::-1]
...     return int(np.searchsorted(np.cumsum(weights), fraction) + 1)
...
>>> # the same Hamiltonian in the position basis has an identical spectrum, but a dense ground state
>>> position_linop = ffsim.linear_operator(hamiltonian_position, norb, nelec)
>>> position_energy, position_ground = scipy.sparse.linalg.eigsh(position_linop, k=1, which="SA")
>>>
>>> print(f"position-basis energy: {position_energy[0]:.6f}")
position-basis energy: -13.422492
>>> print(f"determinants for 99% weight: {determinants_for(position_ground[:, 0])} (position)")
determinants for 99% weight: 2563 (position)
>>> print(f"determinants for 99% weight: {determinants_for(ground_state)} (momentum)")
determinants for 99% weight: 21 (momentum)

2つの演算子のエネルギーは同じですが、位置基底の基底状態は数千個の行列式に広がっているのに対し、運動量基底ではおよそ20個に集中しています。これが、SKQDのサンプリングが位置基底ではなく運動量基底で機能する理由です。

3. 参照状態の準備

SKQDの参照状態$\lvert \psi_0 \rangle$は、フェルミ準位に最も近い電子を、最も近い空の運動量モードへ励起したすべての状態の重ね合わせです。SKQDの論文では、これを手書きのXXPlusYYGateのネットワークで構築しています。強く相関しているように見えます(計算基底で非ゼロの振幅が数百個あります)が、粒子数を保存する2モード回転のネットワークは1粒子(フェルミオン・ガウス)操作であるため、生成される状態は単一のスレーター行列式です。つまり、このネットワークが実装する1粒子回転が分かれば、1つのPrepareSlaterDeterminantゲートで準備できます。

各XXPlusYYGateはモードのペアに対するギブンス回転として作用するため、2量子ビットゲートを手作業で出力する代わりに、その回転を1粒子の描像で直接構築できます(ゲートごとに1つのギブンス行列)。これらのゲートは、フェルミ準位から拡大する層状に広がる最近接回転のレンガ状配置を形成し、フェルミ準位付近の各電子を、そのすぐ上の空のモードに広げます。

>>> def fermi_level_rotation(norb, nocc, n_excited=3):
...     """Orbital rotation exciting the near-Fermi electrons into the nearest empty modes."""
...     cos, sin = np.cos(np.pi / 4), np.sin(np.pi / 4)
...     fermi = nocc - 1  # highest occupied mode; the Fermi level sits between it and mode `nocc`
...     # Nearest-neighbor Givens rotations fanning out from the Fermi level in expanding layers:
...     # layer `depth` couples the mode `depth` steps below the level to the one above it, mixing
...     # the `n_excited` highest-occupied modes with the empty modes just above them. A final
...     # rotation carries amplitude up into the topmost empty mode of the window.
...     network = [
...         (lower, lower + 1)
...         for depth in range(n_excited)
...         for lower in range(fermi - depth, fermi + depth + 1, 2)
...     ]
...     network.append((fermi + n_excited, fermi + n_excited + 1))
...     unitary = np.eye(norb, dtype=complex)
...     for i, j in network:
...         givens = np.eye(norb, dtype=complex)
...         givens[i, i] = givens[j, j] = cos
...         givens[i, j], givens[j, i] = sin, -sin
...         unitary = givens @ unitary
...     return unitary
...
>>> reference_rotation = fermi_level_rotation(norb, nelec[0])

注

上記の回転は、ハミルトニアンに合わせて運動量基底で表現されています。フェルミ準位付近の電子は、占有された運動量モードのうち最もエネルギーの高いものであり、ネットワークはそれらを最もエネルギーの低い空のモードへ励起します。このため、参照状態は前のステップの基底変換の前ではなく、後に構築する必要があります。

スレーター行列式の準備はスピンセクターごとに行います。各セクターに固有の(電子数、軌道数)の形状があるためです。単一のPrepareSlaterDeterminantゲート(占有状態として最も低いnocc個のモードを埋めたものと、1粒子回転)が1つのセクターを準備します。ここでは両方のスピンセクターが同じ参照占有状態と回転を共有するため、同じゲートをアルファモード0..norbに1回、ベータモードnorb..2*norbに1回適用します。

>>> from qiskit_fermions.circuit.library import PrepareSlaterDeterminant
>>>
>>> occupation = [i < nelec[0] for i in range(norb)]
>>> reference_preparation = PrepareSlaterDeterminant(occupation, reference_rotation)

この宣言的なゲートは、手書きのXXPlusYYGateネットワークを置き換えます。ギブンス回転の合成はトランスパイラーが生成します。

4. クリロフ回路の組み立て

各クリロフ回路は、参照状態の準備と、それに続く$e^{-i k \, \Delta t \, H}$による時間発展で構成されます。ハミルトニアンを1体部分と2体部分$H = H_1 + H_2$に分割し、その分割に沿って時間発展をトロッター分解します。これは回路の深さにとって決定的な選択です。1体の時間発展$e^{-i t H_1}$は自由フェルミオン(フェルミオン・ガウス)操作なので、項ごとにパウリ・トロッター分解する必要はなく、1粒子基底回転(OrbitalRotation、ユニタリは$e^{-i t \, \mathbf{h}_1}$)になります。2体の時間発展$e^{-i t H_2}$は、単一のオンサイト相互作用項にのみ作用します。このようにステップを構築すると、トランスパイラーは浅いギブンス回転のレンガ状配置と、単一の2量子ビット相互作用ゲートを出力できます。これは、SKQDの論文が手作業で構築している浅い構造です。

2体部分はすでにあります。これはステップ2で構築したtwo_body_momentum演算子(単一のオンサイト項)です。1体の時間発展$e^{-i t H_1}$はスピンに依存しないため、各スピンセクターに1回ずつ適用する、セクターごとの単一のOrbitalRotation(ユニタリは$e^{-i t \, \mathbf{h}_1}$)になります。これは参照状態の準備と同じ並列構造です。

各クリロフ回路は、合計時間$k \, \Delta t$だけ時間発展させます。これは、サイズ$\Delta t$の$k$ステップの2次トロッター積によって実現されるため、クリロフ次元が増えてもステップごとの誤差は一定に保たれます。各ステップでは、(安価な単一項の)2体時間発展の2つの半ステップの間に、完全な1体回転(スピンセクターごとに1つ)を挟みます。クリロフ次元$D$はそのような回路の数であり、dimは0(参照状態のみ)から$D - 1$までの値を取ります。

>>> import scipy.linalg
>>>
>>> from qiskit_fermions.circuit import FermionicCircuit
>>> from qiskit_fermions.circuit.library import Evolution, OrbitalRotation
>>>
>>> num_modes = 2 * norb
>>> time_step = 0.2
>>>
>>> half_exp_h2 = Evolution(num_modes, two_body_momentum, time_step / 2)
>>> full_exp_h1 = OrbitalRotation(scipy.linalg.expm(-1j * time_step * h1e_momentum))
>>>
>>> def krylov_circuit(dim):
...     # closes over the fermionic gates defined in the outer scope:
...     # `reference_preparation`, `half_exp_h2` and `full_exp_h1`
...     circuit = FermionicCircuit(num_modes)
...     circuit.append(reference_preparation, range(norb))
...     circuit.append(reference_preparation, range(norb, num_modes))
...     for _ in range(dim):  # second-order Trotter product of `dim` steps
...         circuit.append(half_exp_h2, circuit.modes)
...         circuit.append(full_exp_h1, range(norb))
...         circuit.append(full_exp_h1, range(norb, num_modes))
...         circuit.append(half_exp_h2, circuit.modes)
...     return circuit

注

ここで構築する回路には、まだ測定が含まれていません。ステップ6の厳密シミュレーションでも同じフェルミオン回路を使用し、測定はユニタリ操作ではないためです。ステップ5で、モードiを古典ビットiに読み出すFermionicCircuit.measure_all()で測定を追加します。

トランスパイル後の回路ではなく、フェルミオン回路を測定することをお勧めします。フェルミオンの測定は読み出すモードを指定するため、レイアウトとルーティングによって、そのモードが配置された量子ビットとともに引き継がれ、サンプリングされたビット列はモードでインデックス付けされます。代わりにトランスパイル後の回路を測定すると、final_index_layout()を自分で元に戻す必要があり、忘れると、置換された結果が気付かれないまま返されます。これは、プリセットが使用する占有基底のジョルダン・ウィグナー符号化(各モードが1つの量子ビットで表される)に依存しています。

5. 量子ビット回路へのトランスパイル

FermionicCircuitをトランスパイルすると、そのフェルミオンゲートが量子ビットゲートにマッピングされます。ここではジョルダン・ウィグナーjordan_wigner()マッピングを使用します。OrbitalRotationゲートはギブンス回転のレンガ状配置に合成され、単一項のEvolutionは単一の2量子ビット相互作用ゲートに分解されます。

ここでのすべてのフェルミオンゲートにはデフォルトの合成があるため、各ステージを手作業で組み立てるのではなく、既成のgenerate_preset_jw_pass_manager()を直接使用できます。キーワード引数は量子ビットステージに転送されます。合成された深さを最適化なしで忠実に確認するには、optimization_level=0を渡します。クリロフ回路の全体は、クリロフ次元に対するループで生成できます。トランスパイルしていないFermionicCircuitを、トランスパイル後の量子ビット回路と一緒に保持してください。フェルミオン回路はステップ6の厳密(ステートベクトル)シミュレーションで使用するため、測定を含めてはいけません。measure_all(inplace=False)は、トランスパイルする測定付きのコピーを返し、元の回路には手を加えません。このコピーが、バックエンドが実行するものになります。

>>> from qiskit_fermions.transpiler.presets import generate_preset_jw_pass_manager
>>>
>>> pass_manager = generate_preset_jw_pass_manager(optimization_level=0)
>>>
>>> def krylov_circuits(krylov_dim):
...     fermionic_circuits, circuits = [], []
...     for dim in range(krylov_dim):
...         fermionic_circuit = krylov_circuit(dim)
...         fermionic_circuits.append(fermionic_circuit)
...         measured = fermionic_circuit.measure_all(inplace=False)
...         circuits.append(pass_manager.run(measured))
...     return fermionic_circuits, circuits
...
>>> krylov_dim = 5
>>> fermionic_circuits, circuits = krylov_circuits(krylov_dim)

クリロフ次元は、回路の2量子ビット深さに直接現れます。参照状態の準備は固定コストであり、クリロフのべきが1つ増えるごとに、時間発展演算子の2次トロッターステップがもう1つ加わります。

>>> two_qubit = lambda instr: instr.operation.num_qubits == 2
>>> [circuit.depth(two_qubit) for circuit in circuits]
[4, 12, 21, 30, 39]

注

ハミルトニアンを分割することで、これらの深さは控えめに抑えられます。1体の時間発展は純粋な軌道回転であるため、隣接する量子ビット上のXXPlusYYGateのレンガ状配置に合成されます。オンサイト相互作用の各半ステップは、不純物の2つのスピンモードを結合する単一のRZZGateです。これは、2次の積が回転の両側に挟む安価な項です。代わりにハミルトニアン全体を1つのEvolutionゲートに渡すと、1体部分が項ごとにパウリ・トロッター分解され、(長距離の)不純物・浴の結合にわたって長いジョルダン・ウィグナーのZストリングが伝播し、深さが桁違いに増大します。

最初の回路は参照行列式のみを準備し、最後の回路は最も多くのトロッターステップ(この場合は4つ)を含みます。

>>> circuits[0].draw("mpl", fold=-1, measure_arrows=False)
<Figure size ... with 1 Axes>

>>> circuits[-1].draw("mpl", fold=-1, measure_arrows=False)
<Figure size ... with 1 Axes>

6. クリロフ回路からのビット列のサンプリング

ハードウェアでは、トランスパイルした回路を実行して測定します。ここでは、そのステップのノイズなし(ステートベクトル)版を示し、SKQDが収集するビット列を見ます。トランスパイルしていない各FermionicCircuitをffsim.apply_unitary()でシミュレーションし、得られたステートベクトルをffsim.sample_state_vector()でサンプリングします。

PrepareSlaterDeterminantゲートは、シミュレーションでは検証してから回転を適用します(入力状態が参照行列式を占有していることを確認し、その後で回転を適用します)。したがって、入力は単純な占有行列式、つまりスピンごとに最も低いnocc個のモードを埋めた、回転なしの状態です。ゲートがreference_rotation自体を適用します。

>>> from collections import Counter
>>>
>>> reference = ffsim.slater_determinant(
...     norb, (list(range(nelec[0])), list(range(nelec[1])))
... )
>>>
>>> shots = 1000
>>> counts = Counter()
>>> for dim, fermionic_circuit in enumerate(fermionic_circuits):
...     statevector = ffsim.apply_unitary(reference, fermionic_circuit, norb=norb, nelec=nelec)
...     samples = ffsim.sample_state_vector(
...         statevector, norb=norb, nelec=nelec, shots=shots, seed=dim
...     )
...     counts += Counter(samples)

5つのクリロフ回路すべてからのサンプルをまとめると、サポートは非常に小さくなります。(4, 4)セクターの$\binom{8}{4}^2 = 4900$個の行列式のうち、異なるビット列は数百個にすぎません。この集中は、ステップ2で確立した運動量基底の疎性の直接的な帰結であり、その後の古典的な対角化を扱いやすくします。

qiskit.visualization.plot_histogram()でカウントをプロットすると、スパース性が見て取れます。参照行列式を筆頭に、ごく少数の配置が支配的です。

>>> from qiskit.visualization import plot_histogram
>>>
>>> plot_histogram(dict(counts), number_to_keep=25)
<Figure size ... with 1 Axes>

7. サンプリングした部分空間での対角化

サンプリングしたビット列は、SKQDの古典側の処理への入力になります。ハミルトニアンを、サンプリングした配置が張る部分空間に射影し、その空間で対角化します。この処理はqiskit-addon-sqdパッケージに渡され、そのdiagonalize_fermionic_hamiltonian()がサンプルベース・クォンタム対角化のループ全体を実行します。サンプリングした配置のバッチから部分空間を構築し、その中でハミルトニアンを対角化し、配置回復を用いて部分空間を反復的に洗練します。配置回復とは、既知の粒子数対称性に違反する配置を、占有を反転させて修復するものです。配置回復は、ハードウェアのノイズによってサンプリングされたビット列が誤った粒子数セクターに壊されてしまう影響を取り除くために設計されています。ここでのサンプリングはノイズがないため、その点で修復すべきものはありません。

ソルバーは、ステップ6でプールしたcountsから直接構築したBitArrayとしてサンプリングされたカウントを受け取ります。num_bitsは明示的に渡してください。渡さない場合、幅は最も幅の広いサンプリング済みビット列から推定されます。最後のモードを占有するサンプルがないときは、その幅が2 * norbより狭くなり、ソルバーは配列を拒否します。ソルバーは、構成の他の部分と同じ運動量基底で対角化するため、ステップ2のh1e_momentumとh2e_momentumの積分を、再パッキングなしでそのまま受け取ります。回復されたエネルギーは、わずか数百のサンプリング配置から、厳密な参照値と1ミリ・ハートリー未満の精度で一致します。SQDのエネルギーは真の基底状態に対する変分上界なので、常にreference_energy(ステップ2から引き続きスコープ内にあります)以上になります。

>>> from qiskit.primitives import BitArray
>>> from qiskit_addon_sqd.fermion import diagonalize_fermionic_hamiltonian
>>>
>>> bit_array = BitArray.from_counts(dict(counts), num_bits=2 * norb)
>>>
>>> # record each iteration's batch results so we can watch the energy converge
>>> result_history = []
>>> result = diagonalize_fermionic_hamiltonian(
...     h1e_momentum,
...     h2e_momentum,
...     bit_array,
...     samples_per_batch=100,
...     norb=norb,
...     nelec=nelec,
...     num_batches=3,
...     max_iterations=3,
...     symmetrize_spin=True,
...     callback=result_history.append,
...     seed=np.random.default_rng(24),
... )
>>> print(f"SQD estimate: {result.energy:.4f} (exact: {reference_energy:.4f})")
SQD estimate: -13.4224 (exact: -13.4225)

実行を可視化します。反復ごとのエネルギーの収束と、回復された基底状態における各空間軌道の平均占有数です。反復ごとのエネルギーは、その反復のバッチの中で最も低い値で、占有数は2つのスピンセクターにわたるorbital_occupanciesを合計したものです。補正すべきノイズがなくても、エネルギーは反復ごとに低下します。各反復は、プールからサブサンプリングしたsamples_per_batch個だけの配置からなるnum_batches個のバッチ内で対角化を行います。ビット列の引き継ぎによる反復的な洗練が、どの配置がそれらのバッチに入るかを着実に改善するため、回復されるエネルギーも改善します。

>>> import matplotlib.pyplot as plt
>>> from matplotlib.ticker import MaxNLocator
>>>
>>> min_energies = [min(batch, key=lambda r: r.energy).energy for batch in result_history]
>>> occupancies = np.sum(result.orbital_occupancies, axis=0)
>>>
>>> fig, axs = plt.subplots(1, 2, figsize=(12, 5))
>>> _ = axs[0].plot(range(len(min_energies)), min_energies, marker="o", label="SQD energy")
>>> _ = axs[0].axhline(reference_energy, color="gray", linestyle="--", label="exact energy")
>>> _ = axs[0].set_xlabel("iteration")
>>> _ = axs[0].set_ylabel("energy")
>>> _ = axs[0].set_title("Energy vs SQD iteration")
>>> _ = axs[0].xaxis.set_major_locator(MaxNLocator(integer=True))
>>> _ = axs[0].legend()
>>> _ = axs[1].bar(range(norb), occupancies)
>>> _ = axs[1].set_xlabel("spatial orbital")
>>> _ = axs[1].set_ylabel("average occupancy")
>>> _ = axs[1].set_title("Occupancy per spatial orbital")
>>> fig
<Figure size ... with 2 Axes>

これでSKQDのループが完結します。ステップ1〜5の運動量基底の回路がスパースなサンプルを生成し(ステップ6)、上記のサンプルベース対角化がそのサンプルを基底状態エネルギーに変換します。

次のステップ

このガイドでは、ノイズのない状態ベクトル・シミュレーターでSKQDパイプライン全体を実行しました。ハードウェアでは、ステップ6の状態ベクトルによるサンプリングの代わりに、トランスパイルしたcircuitsを実行して測定し、得られたカウントをそのままステップ7の対角化に渡します。回路の実行についてはQiskitドキュメントを、部分空間対角化の後処理(配置回復が最も力を発揮するノイズを含むサンプルでの動作を含む)についてはSQDアドオンのチュートリアルを参照してください。

ffsim.apply_unitary()とffsim.sample_state_vector()がこのパッケージのフェルミオン回路でネイティブに動作する理由と、ステップ6で使う固定粒子数セクターによって、ステップ2で確立した運動量基底のスパース性がサンプル効率の高いものになる理由を理解するには、ffsimの関係ガイドをお読みください。