NVIDIA CUDA-Q チュートリアル＆H100ベンチマーク３： 鈴木-トロッター近似による磁化の計算
（計算速度が線形に変化してるのがあっているのかちょっと未確認ではありますが）
量子コンピュータの代表的な応用のひとつに、量子多体系の時間発展のシミュレーションがあります。特に、スピン系モデルにおける磁化などの物理量を高精度かつ高速に再現することは、材料科学や量子情報の分野で極めて重要です。
本記事では、CUDA-Qを用いてハイゼンベルク模型に従うスピン鎖の磁化を計算するチュートリアルを紹介します。このシミュレーションでは、時間依存ハミルトニアンに対して鈴木–トロッター（Suzuki–Trotter）近似を適用し、量子状態の時間発展を効率的に追跡します。
CUDA-Qが提供する状態管理機能を活用することで、各ステップでの量子状態をGPUメモリ上に保持し、再帰的かつ高速に時間発展を継続できます。これにより、従来の「毎ステップでの再シミュレーション」が不要になり、最大24倍の計算高速化が実現されることを、本記事の例では示します。
ハミルトニアンの構築から、初期状態の準備、磁化期待値の計算、Trotter展開による反復的シミュレーション、さらには期待値の時間変化の可視化までを、CUDA-Qの機能を通じて一貫して実装する方法を、ステップごとに詳しく解説します。
プログラムの準備
最初にライブラリの読み込みです。
import cudaq
import time
import numpy as np
from typing import List
次に今回のモデルのハミルトニアンを確認します。
ハイゼンベルグハミルトニアンは下記の通り与えられています。
最初にハミルトニアンをPythonの形に書き下します。
g = 1.0
Jx = 1.0
Jy = 1.0
Jz = g
dt = 0.05
n_steps = 10
n_spins = 11
omega = 2 * np.pi
def heisenbergModelHam(t: float) -> cudaq.SpinOperator:
tdOp = cudaq.SpinOperator(num_qubits=n_spins)
for i in range(0, n_spins - 1):
tdOp += (Jx * cudaq.spin.x(i) * cudaq.spin.x(i + 1))
tdOp += (Jy * cudaq.spin.y(i) * cudaq.spin.y(i + 1))
tdOp += (Jz * cudaq.spin.z(i) * cudaq.spin.z(i + 1))
for i in range(0, n_spins):
tdOp += (np.cos(omega * t) * cudaq.spin.x(i))
return tdOp
最初に初期状態を決めていきます。
@cudaq.kernel
def getInitState(numSpins: int):
q = cudaq.qvector(numSpins)
for qId in range(0, numSpins, 2):
x(q[qId])
次に時間発展を記述します。
時間発展ではそのままではハミルトニアンが展開できないので、鈴木トロッター展開を適用します。
こちらではトロッター展開された時間発展を1ステップずつ記録しながら計算をします。
@cudaq.kernel
def trotter(state: cudaq.State, coefficients: List[complex],
words: List[cudaq.pauli_word], dt: float):
q = cudaq.qvector(state)
for i in range(len(coefficients)):
exp_pauli(coefficients[i].real * dt, q, words[i])
次は、ハミルトニアンの係数と演算子を抜き出すプログラムです。
def termCoefficients(op: cudaq.SpinOperator) -> List[complex]:
result = []
ham.for_each_term(lambda term: result.append(term.get_coefficient()))
return result
def termWords(op: cudaq.SpinOperator) -> List[str]:
result = []
ham.for_each_term(lambda term: result.append(term.to_string(False)))
return result
そして最後は計算したい平均磁化で、Z方向のスピンの平均をとってます。
average_magnetization = cudaq.SpinOperator(num_qubits=n_spins)
for i in range(0, n_spins):
average_magnetization += ((1.0 / n_spins) * cudaq.spin.z(i))
average_magnetization -= 1.0
実際のシミュレーション
上記準備したものを利用してシミュレーションします。
まず初期状態です。
state = cudaq.get_state(getInitState, n_spins)
次にシミュレーションです。
results = []
times = []
for i in range(0, n_steps):
start_time = time.time()
ham = heisenbergModelHam(i * dt)
coefficients = termCoefficients(ham)
words = termWords(ham)
magnetization_exp_val = cudaq.observe(trotter, average_magnetization, state,
coefficients, words, dt)
result = magnetization_exp_val.expectation()
results.append(result)
state = cudaq.get_state(trotter, state, coefficients, words, dt)
stepTime = time.time() - start_time
times.append(stepTime)
print(f"Step {i}: time [s]: {stepTime}, result: {result}")
print(f"Step times: {times}")
print(f"Results: {results}")
これで、CPUは11スピン、GPUは25スピンくらいまで計算できるようです。ステップも多くなると限界がきますが、それも比較してみます。
以下CPUとGPUでの比較です。
- CPU: AMD EPYC 9654
- GPU: NVIDIA H100
計算がかなり重たいので、CPUは6-11スピン、GPUは6-25スピンを計算しました。
GPUはかなり余裕ですね。
次にCPUは重たいので、GPUでスピンを20に固定して、ステップ数を10から100まで10刻みで測定しました。
結構高速です。以上です。