QunaSys Tech Blog

Quantum-Selected Configuration Interaction（QSCI）では、量子デバイスを用いてAnsatz状態をサンプリングし、対象とする状態にとって重要な電子配置、すなわち行列式を提案します。続いて古典ソルバーが、選択された部分空間内でハミルトニアンを構築し、対角化します。得られるエネルギーは変分的であり、選択した部分空間内では厳密で、部分空間を拡大するにつれて単調に減少します [1, 2]。したがって、より大きな部分空間を選択することで結果を改善していくことができますが、それは古典ソルバーがその規模の部分空間を扱える場合に限られます。

QunaSysでは、QURI SDK Enterprise向けに高性能なQSCIを開発しています。今回の開発以前は、その古典計算段階で、部分空間のハミルトニアンを1つのMessage Passing Interface（MPI）ランク上で構築・対角化する必要があり、扱える問題の規模は1プロセスで利用可能なメモリに制約されていました。そこで、各MPIランクが行列の一部だけを保持するように両段階を分散化しました。新しいソルバーでは、富岳上で400万行列式を扱うことができます。100万行列式では、ハミルトニアン構築が4ランクから96ランクへの拡張で11.8倍高速化し、また、実行結果から新たなボトルネックも明らかになりました。最大の部分空間では対角化が支配的となり、800万行列式へのすべてのケースで構築中にメモリ不足が発生しました。計算コストは量子ビット数よりもハミルトニアンの密度によって、より確実に予測することができます。

図1. 32ランクにおける、選択した行列式数に対するエネルギー誤差。H14では独立に計算したFCIエネルギーを基準とし、bare P450モデルでは18,150行列式からなる完全活性空間で計算したエネルギーを基準としています。縦軸が対数であるため、誤差がゼロとなるP450の終点は省略しています。

このソルバーはQURI SDK Enterprise向けに開発された分散QSCI実装であり、近日中にリリースされる予定です。ベンチマークは富岳上で実行しました。部分空間ハミルトニアンをランク間で行方向に分割して構築するため、単一のランクが行列全体を保持することはありません。その後、分散行列の保存にはPETScを、分散Krylov–Schur法による求解にはSLEPcを使用します [3, 4]。

実行にはFujitsu A64FXノード [5] を使用しました。各ノードは48コアと32 GiBの高帯域幅メモリ（HBM2）を備え、すべての実行で OMP_NUM_THREADS=1 としました。以降で具体的に示すメモリ上の理由から、ノードあたりのランク数は少なく設定しています。400万行列式のすべての測定点を含む32ランク固定の実行では、ノードあたり1ランク、すなわち32ノードを使用し、各ランクがノード全体の32 GiBを利用できるようにしました。これらの固定サイズ実行は、強スケーリングの測定とは別です。4〜16ランクのジョブではノードあたり1ランク、32〜96ランクのジョブではノードあたり4ランクを使用しました。したがって、96ランクの強スケーリング測定では24ノードを使用しています。100万行列式のH18には、約4ノード分の全メモリが必要です。

今回のベンチマークでは、9つの系を測定しました。結合を引き伸ばした直線状の水素鎖H14、H18、H20、H22の4系は28〜44量子ビットにわたり、強相関かつ比較的疎なストレスケースです。シトクロムP450活性部位の4つのモデルは22〜30量子ビットで、より密な現実的化学系です。これらは軸配位子が異なります。6-31G基底の水分子は、特性が完全に分かっている対照系です。量子ビット数よりハミルトニアン密度の影響が大きいため、疎な系と密な系を分けて報告します。

行列式集合には、一様ランダムサンプリングではなく、full configuration interaction（FCI）またはheat-bath configuration interaction（HCI）[2] から得たものを使用しました。HCIは推定重要度の基準に基づいて行列式を順位付けし、強く結合した配置を優先するため、これらの部分空間は行列構築と対角化にとって負荷の高い入力となります。以下の実行では量子ハードウェアを使用していません。これにより、サンプラーの統計的性質から切り離して、出自が明確な行列式集合に対するソルバー性能を評価できます。

以降では、古典計算ステップを、部分空間ハミルトニアンを組み立てる H-build と、それを対角化する diag に分け、部分空間サイズを独立変数として評価します。QSCIの各呼び出しでは2回の実行時間を記録し、合計時間が短い方を報告します。2回の実行はジョブ全体を繰り返すのではなく、1回のリソース割り当てを共有しているため、得られた比率は点推定として扱ってください。記事後半のSBD比較は、各測定点につき1回の計測です。

H14では、diagの割合は1万行列式での25%から、400万行列式での79%まで上昇します。水素鎖では、H14は100K、H18とH20は400K、H22は100万行列式でdiagが過半のコストを占めるようになります。測定した最大部分空間では、すべての系でdiagが支配的です。この測定結果を受け、最適化対象を行列の組み立てから対角化へ移しました。

図2. 2つの水素鎖とP450 oxoモデルについて、32ランクにおける部分空間サイズに対するH-buildとdiagの割合。

上表の100Kの2行を比較してください。H18は36量子ビットの疎な水素鎖で、対角化に0.9秒かかります。一方、30量子ビットのP450 oxoモデルは、ランク数と部分空間サイズが同じでも3.4秒を要します。量子ビット数が少ないにもかかわらず、ハミルトニアンがより密であるため、コストは約4倍になります。

参照ケースである100万行列式のH18では、H-buildは4ランクから96ランクで11.8倍にスケールし、ランク数を24倍に増やすことで90.9秒から7.7秒まで短縮しました。テストした最大ランク数である96ランクでも、H-buildはスケールし続けていました。この曲線には2種類の構成が含まれています。4〜16ランクのジョブでは各ランクに1ノード全体を割り当て、32〜96ランクのジョブでは各ノードに4ランクを配置しました。したがって、これは均一な配置によるスケーリング測定ではなく、実際に使用した構成での経過時間を示しています。

図3. 100万行列式のH18について、ランク数に対するH-buildとdiagの実時間。配置は、16ランクまではノードあたり1ランク、32ランク以降はノードあたり4ランクに変わります。

対角化もほぼ同じようにスケールします。同じ4〜96ランクの範囲で11.8倍となり、163.8秒から13.9秒まで短縮しました。4〜96ランクの全系列に T = T_s + T_p/P をフィッティングすると、逐次実行部分の割合はH-buildで2.0%、diagで3.4%となります。したがって、2つのフェーズは並列化効率が異なるのではなく、絶対コストが異なります。測定したすべてのランク数で、diagはH-buildのおよそ2倍です。古典計算ステップがdiag律速なのは、diagのスケーリングが悪いからではなく、diag自体のコストが大きいからです。

今回の測定ではハミルトニアン構築と対角化のみを対象としており、サンプリング性能は対象外です。

32 MPIプロセスでは、各プロセスが保持する座標形式配列は、100万行列式で約0.5 GiB、200万行列式で1.2 GiB、400万行列式で3.0 GiBでした。直接記録した非ゼロ要素数は、系に応じて dim^1.32 から dim^1.39 で増加しました。これは1万〜400万行列式の範囲で測定したもので、決定係数はいずれも0.995以上です。部分空間サイズが大きくなるにつれて負荷分散が改善するため、ランクあたりのメモリ増加はこれより緩やかです。4つの水素鎖では、最も負荷の高いランクが保持する量は、25K行列式では均等配分時の2.4〜2.9倍、400万行列式では1.3〜2.2倍でした。32ビット整数のインデックスと64ビット実数を使用しているため、行列自体の非ゼロ要素あたりの容量は一律16バイトです。

これらの値は座標形式の行列に関するもので、必要量そのものではなく下限です。構築時には、座標形式配列、PETScがそこから構築する圧縮構造、PETScの座標インデックスマップが同時にメモリ上に存在する場合があります。バイト数からの見積もりでは、この一時的な使用量は記録された座標形式配列の約2.7倍です。富岳のプロセスレベルのメモリカウンターではこれらの割り当てを取得できなかったため、この倍率はピークメモリの直接測定ではなく推定値です。固有値ソルバーのKrylov基底はランクあたり数十MBを追加しますが、支配的な要因ではありません。測定した400万行列式の容量にこの推定を適用すると、ランクあたり約8 GiBとなります。ここでは、ノードのメモリをコア数で分割して考えることはできません。400万行列式の実行では、1ランクに48コアノード全体を割り当て、行列を配置する場所を確保するためだけに、ノードメモリの4分の1強と48コア中1コアのみを使用しています。

図4. 4つの水素鎖について、部分空間サイズに対するランクあたりの座標形式配列の測定容量。

800万行列式まで拡張することで、構築処理が次の改善対象であることを特定できました。試したすべての水素鎖で、定常状態の使用量なら収まるにもかかわらず、構築中にメモリ不足となりました。

ノードを追加すればランクあたりのメモリが減るという理由から、これをソフトリミットと呼ぶのは簡単です。しかし、私たちの弱スケーリングの測定結果は、その解釈を支持しません。ランクあたりの行列式数を25,000に固定し、ジョブを1ランクから96ランクまで拡大すると、H-build時間は4.8倍、ランクあたりのメモリは13.4倍、diag時間は34倍になりました。したがって、より大きな部分空間を保持するためにノードを追加しても、一定の実行時間でその部分空間を扱えるわけではありません。メモリの壁を、そのすぐ先にある時間の壁へ置き換えるだけです。

容量を追加するだけでは問題規模は拡大できません。そもそも行列を構築しないことが必要です。行列フリーの演算子適用であれば、構築時のピークメモリ使用量と定常状態の行列容量の両方を回避できます。また、これは固有値ソルバーの改善でも通る必要があるコードパスです。

同じ問題を、Sample-based Quantum Diagonalization（SQD）の qiskit-community/qiskit-c-api-demo 実装で使われている固有値ソルバー、Selected Basis Diagonalization（SBD）でも実行しました [6, 8]。SQDは、QSCIを基盤とするワークフローに、configuration recoveryとそれに関連する後処理を加えたものです。したがってこのベンチマークは、部分空間の次元をおおむね揃えた上で、選択された部分空間を古典的に対角化するステップについて、2つの実装を比較するものです。両者とも同じランク数、OMP_NUM_THREADS=1 で実行し、時間にはハミルトニアンのファイルI/Oとセットアップを含めていません。

テストしたすべてのP450部分空間で、私たちのソルバーの方が高速でした。2,000行列式では1.8〜3.5倍、各モデルの最大部分空間では、water配位子で8.5倍、oxoで11.5倍、pyridineで16.9倍、18,150行列式が完全活性I空間となるbareモデルで15.5倍高速でした。

両者とも完全なFCI順位付き行列式プールから開始しますが、部分空間の構築方法が異なります。QSCIは、順位の高いα-β行列式ペアを1つのフラットなリストとして保持します。一方、SBDへのブリッジではαストリングとβストリングを別々に順位付けし、それらの直積を対角化します。100,000行列式では、2つの空間の重なりはwater配位子で58%、pyridineとoxoで47%です。完全な18,150行列式のCI空間に到達するbareモデルでのみ、2つの行列が一致します。P450の全22測定点における2つのエネルギーの差は最大0.36 mHaで、25,000行列式以上のすべての部分空間では0.01 mHa未満です。したがって、実行時間の比は、bareモデルの完全空間の測定点を除けば、同一行列ではなく、ほぼ同じ次元で近いエネルギーを与える異なる部分空間同士を比較したものです。

図5. 4つのP450モデルについて、32ランクにおける目標部分空間サイズに対するQSCIとSBDの実時間。次元はおおむね一致させていますが、完全空間の終点を除き、行列式集合は異なります。

• 実装。 比較対象は qiskit-c-api-demo のCPU実装（commit 47dbfb7、SBD submodule commit f0b85ad）です [8]。事前計算した部分空間を実行し、αとβの行列式ファイルを別々に受け取り、処理別の実行時間を記録するため、ローカルのベンチマーク用パッチを適用しました。パッチはご要望に応じて提供できます。SBDにはGPUアクセラレーション版もあります [7]。今回の実行はA64FX上のCPUのみを使用しているため、GPU版との比較は対象外です。
• 最適化。 今回の測定では、私たちの実装を継続的にプロファイリングしましたが、SBDは同程度には調整していません。SBDにも私たちのコードで見つかったものと同程度の改善余地があれば、この差は縮まる可能性があります。
• SBDの設定。 制限時間のデフォルトは600秒で、未完了の結果を通知せずに返す場合があるため、P450の一連の測定では7,200秒に設定しました。また、反復回数の上限を200、ブロックサイズを40、許容誤差を 1e-9 に設定しました。大きな部分空間では、ブロックサイズ40はデフォルトの10に比べて約2倍高速でした。
• 許容誤差。 固有値ソルバーの許容誤差は完全には一致していません。SLEPcのデフォルトは 1e-8、SBDはより厳しい 1e-9 です。別途実施した4ランクでの対照実験では、QSCIを 1e-9 で再実行すると、H18と2つのP450ケースで対角化時間が4〜14%増加しましたが、エネルギーに意味のある変化はありませんでした。観測された最大の変化は1.14倍で、報告した8.5〜16.9倍という差よりはるかに小さい値です。

6-31G基底の水分子では、私たちのSLEPcによる求解は、200,000行列式のフラットな部分空間から-76.1208 Haを返します。一方、SBDは199,809行列式を含む447×447の直積を対角化し、-74.108 Haを返します。粒子数を固定した完全CI空間には1,656,369行列式が含まれるため、この447×447の直積は選択された部分空間にすぎず、完全CI空間ではありません。2つの部分空間で共通する行列式は55,693個、約28%であるため、QSCIのエネルギーはSBDの行列に対する変分上界ではありません。

SBDの直積にはHartree–Fock行列式が含まれ、このハミルトニアンに対するそのエネルギーは-75.984 Haです。その部分空間における最小固有値は、そこに含まれる行列式の期待値より高くなることはありません。SBDの結果はこの上界より1.88 Ha高いため、それ自身の行列の最小固有値ではあり得ません。

SBDは許可された200回の反復のうち15回で停止し、要求された 1e-9 に対して残差4.7e-10を報告して正常終了します。SBD自身の基準では収束していますが、基底状態には収束していません。反復上限を200回から2,000回に増やしても、報告されるすべての桁は変わりません。Davidson法のブロックサイズを増やしたり許容誤差を緩めたりすると異なるエネルギーが得られますが、いずれもHartree–Fock上界より高い値です。これらの設定を変更しても、行列やその基底状態エネルギーは変わりません。

SBDは初期推定をハードコードしており、前処理を外部から指定する手段もないため、どちらの改善策も試すことができませんでした。SBDが正しい基底状態エネルギーを返さないため、比較可能な有効なSBD実行時間はありません。

次の実装目標は、行列フリー固有値ソルバーです。これは今回明らかになった2つの制約に対処します。最大の部分空間では対角化が実行時間の大半を占める一方、ハミルトニアンを実体化するとピークメモリが増大し、800万行列式で計算が停止します。行列を取り除けば、主にメモリ確保のためにノード全体を割り当てるのではなく、各ノードの計算能力をより有効に活用できるようになります。

実用上の利点は、単一プロセスの容量上限を超える選択部分空間を扱えることです。今回のベンチマークでは、量子サンプリングではなくFCIまたはHCIから得た行列式を使用しました。次の測定では、行列フリーソルバーと量子ハードウェアからサンプリングした行列式を組み合わせ、強相関系を対象に、精度とコストをHCI、密度行列繰り込み群（DMRG）、coupled clusterと比較することで、このエンジニアリング上の成果を化学へつなげます。

ソフトウェアとバージョン：

• QURI SDK Enterprise revision 8314ff0
• Python 3.11.11
• PETSc 3.24.5
• SLEPc 3.24.2
• qiskit-c-api-demo commit 47dbfb7
• SBD commit f0b85ad

本研究は、理化学研究所計算科学研究センターが提供するスーパーコンピュータ「富岳」の計算資源（課題番号 hp260359）を利用して実施しました。

1. K. Kanno, M. Kohda, R. Imai, S. Koh, K. Mitarai, W. Mizukami, Y. O. Nakagawa. Quantum-selected configuration interaction: classical diagonalization of Hamiltonians in subspaces selected by quantum computers. Phys. Rev. Research 8, 023268 (2026). https://doi.org/10.1103/dmn4-snfx.
2. A. A. Holmes, N. M. Tubman, C. J. Umrigar. Heat-bath configuration interaction: an efficient selected configuration interaction algorithm inspired by heat-bath sampling. J. Chem. Theory Comput. 12(8), 3674-3680 (2016). https://doi.org/10.1021/acs.jctc.6b00407.
3. S. Balay et al. PETSc/TAO Users Manual. Argonne National Laboratory, ANL-21/39, Revision 3.24.
4. V. Hernandez, J. E. Roman, V. Vidal. SLEPc: a scalable and flexible toolkit for the solution of eigenvalue problems. ACM Trans. Math. Software 31(3), 351-362 (2005). https://doi.org/10.1145/1089014.1089019.
5. M. Sato et al. Co-design for A64FX manycore processor and Fugaku. SC20: International Conference for High Performance Computing, Networking, Storage and Analysis (2020). https://doi.org/10.1109/SC41405.2020.00051.
6. J. Robledo-Moreno et al. Chemistry beyond the scale of exact diagonalization on a quantum-centric supercomputer. Sci. Adv. 11, eadu9991 (2025). https://doi.org/10.1126/sciadv.adu9991.
7. J. Doi, T. Shirakawa, Y. Kawashima, S. Yunoki, H. Horii. GPU-accelerated selected basis diagonalization with Thrust for SQD-based algorithms. arXiv:2601.16637 (2026).
8. qiskit-community. qiskit-c-api-demo. GitHub, commit 47dbfb70c785fc5465ac60f737e4899318e96741, with r-ccs-cms/sbd commit f0b85ada8f8ad76e123d9e1ad4ed7d769dd08518. https://github.com/qiskit-community/qiskit-c-api-demo/tree/47dbfb70c785fc5465ac60f737e4899318e96741