ハミルトニアンニューラルネットワーク（Hamiltonian Neural Network, HNN）は、物理系のエネルギー保存則をニューラルネットワークの構造に組み込む手法です。Greydanus et al.（2019）によって提案されました。力学系をデータから学習する際、通常のニューラルネットワークではエネルギーが保存されず長時間の予測で大きな誤差が蓄積しますが、HNNはハミルトン力学の構造を利用することで、エネルギー保存を自動的に保証します。
本記事では、ハミルトン力学の数学的基盤を丁寧に復習したうえで、HNNの理論、エネルギー保存の証明、学習の定式化、シンプレクティック積分法との組み合わせ、そしてPythonによる実装と通常のニューラルネットワークとの比較を行います。
本記事の内容
- ハミルトン力学の復習（正準方程式の導出）
- 通常のNNがエネルギー保存則を破る問題
- HNNの構造とエネルギー保存の数学的証明
- 学習の損失関数と訓練方法
- シンプレクティック積分法（リープフロッグ法）
- Pythonで単振り子のHNN実装（通常NNとのエネルギー保存精度比較）
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
ハミルトン力学の復習
ラグランジュ力学からの導出
ニュートン力学を一般化する枠組みとして、まずラグランジュ力学を考えます。一般化座標 $\bm{q} = (q_1, q_2, \dots, q_n)$ を使って系の状態を記述します。ラグランジアン $L$ は運動エネルギー $T$ とポテンシャルエネルギー $V$ の差として定義されます。
$$ L(\bm{q}, \dot{\bm{q}}) = T(\bm{q}, \dot{\bm{q}}) – V(\bm{q}) $$
オイラー＝ラグランジュ方程式は以下です。
$$ \frac{d}{dt}\frac{\partial L}{\partial \dot{q}_i} – \frac{\partial L}{\partial q_i} = 0, \quad i = 1, 2, \dots, n $$
一般化運動量とハミルトニアン
一般化運動量（正準運動量）を以下のように定義します。
$$ p_i = \frac{\partial L}{\partial \dot{q}_i} $$
ハミルトニアン $H$ は、ラグランジアンのルジャンドル変換として定義されます。
$$ H(\bm{q}, \bm{p}) = \sum_{i=1}^{n} p_i \dot{q}_i – L(\bm{q}, \dot{\bm{q}}) $$
ここで右辺の $\dot{\bm{q}}$ は $\bm{q}$ と $\bm{p}$ の関数として表されます（$p_i = \frac{\partial L}{\partial \dot{q}_i}$ を $\dot{q}_i$ について解く）。
多くの力学系では $T = \frac{1}{2}\dot{\bm{q}}^\top \bm{M} \dot{\bm{q}}$（$\bm{M}$ は質量行列）であり、このとき $\bm{p} = \bm{M}\dot{\bm{q}}$、$\dot{\bm{q}} = \bm{M}^{-1}\bm{p}$ となります。ハミルトニアンは以下になります。
$$ H(\bm{q}, \bm{p}) = \frac{1}{2}\bm{p}^\top \bm{M}^{-1} \bm{p} + V(\bm{q}) = T + V $$
つまり、ハミルトニアンは系の全エネルギー（運動エネルギー + ポテンシャルエネルギー）に等しくなります。
ハミルトンの正準方程式
ハミルトニアンの全微分を計算します。
$$ dH = \sum_{i} \frac{\partial H}{\partial q_i} dq_i + \sum_{i} \frac{\partial H}{\partial p_i} dp_i $$
一方、定義 $H = \sum_i p_i \dot{q}_i – L$ から、
$$ \begin{align} dH &= \sum_i (\dot{q}_i dp_i + p_i d\dot{q}_i) – \sum_i \frac{\partial L}{\partial q_i} dq_i – \sum_i \frac{\partial L}{\partial \dot{q}_i} d\dot{q}_i \\ &= \sum_i \dot{q}_i dp_i + \sum_i p_i d\dot{q}_i – \sum_i \frac{\partial L}{\partial q_i} dq_i – \sum_i p_i d\dot{q}_i \\ &= \sum_i \dot{q}_i dp_i – \sum_i \frac{\partial L}{\partial q_i} dq_i \end{align} $$
2行目で $p_i = \frac{\partial L}{\partial \dot{q}_i}$ を使い、$d\dot{q}_i$ の項が相殺されました。
$dq_i$ と $dp_i$ の係数を比較すると、
$$ \frac{\partial H}{\partial p_i} = \dot{q}_i, \quad \frac{\partial H}{\partial q_i} = -\frac{\partial L}{\partial q_i} $$
オイラー＝ラグランジュ方程式 $\frac{d}{dt}\frac{\partial L}{\partial \dot{q}_i} = \frac{\partial L}{\partial q_i}$ より $\dot{p}_i = \frac{\partial L}{\partial q_i}$ なので、
$$ \boxed{\frac{dq_i}{dt} = \frac{\partial H}{\partial p_i}, \quad \frac{dp_i}{dt} = -\frac{\partial H}{\partial q_i}} $$
これがハミルトンの正準方程式です。ベクトル表記では以下のように書けます。
$$ \frac{d}{dt}\begin{pmatrix} \bm{q} \\ \bm{p} \end{pmatrix} = \begin{pmatrix} \phantom{-}\frac{\partial H}{\partial \bm{p}} \\ -\frac{\partial H}{\partial \bm{q}} \end{pmatrix} = \bm{J} \nabla H $$
ここで $\bm{J} = \begin{pmatrix} \bm{0} & \bm{I} \\ -\bm{I} & \bm{0} \end{pmatrix}$ はシンプレクティック行列です。
エネルギー保存の証明
ハミルトニアンが時間に陽に依存しない（$\frac{\partial H}{\partial t} = 0$）とき、$H$ は時間とともに保存されます。証明は以下の通りです。
$$ \begin{align} \frac{dH}{dt} &= \sum_{i} \frac{\partial H}{\partial q_i}\frac{dq_i}{dt} + \sum_{i} \frac{\partial H}{\partial p_i}\frac{dp_i}{dt} \\ &= \sum_{i} \frac{\partial H}{\partial q_i}\frac{\partial H}{\partial p_i} + \sum_{i} \frac{\partial H}{\partial p_i}\left(-\frac{\partial H}{\partial q_i}\right) \\ &= \sum_{i} \frac{\partial H}{\partial q_i}\frac{\partial H}{\partial p_i} – \sum_{i} \frac{\partial H}{\partial p_i}\frac{\partial H}{\partial q_i} \\ &= 0 \end{align} $$
2行目で正準方程式 $\dot{q}_i = \frac{\partial H}{\partial p_i}$, $\dot{p}_i = -\frac{\partial H}{\partial q_i}$ を代入しました。各項が完全に相殺されるため、$\frac{dH}{dt} = 0$ が成立します。
この結果は正準方程式の構造から導かれる性質であり、$H$ の具体的な形によりません。これがHNNの核心的なアイデアにつながります。
通常のNNがエネルギー保存を破る問題
力学系のダイナミクスをデータから学習する場合、素朴なアプローチは以下のようにニューラルネットワークで時間発展を直接モデル化することです。
$$ \frac{d}{dt}\begin{pmatrix} \bm{q} \\ \bm{p} \end{pmatrix} = \text{NN}_{\bm{\theta}}(\bm{q}, \bm{p}) $$
このアプローチでは、ネットワークの出力 $\text{NN}_{\bm{\theta}}(\bm{q}, \bm{p}) = (\bm{f}_q, \bm{f}_p)$ に何の構造的制約もありません。
訓練データ付近では良い近似が得られても、長時間のシミュレーションではエネルギーがドリフト（系統的に増加または減少）します。これは、出力ベクトル場 $(\bm{f}_q, \bm{f}_p)$ がハミルトン系の構造 $\bm{J} \nabla H$ に従わないためです。
具体的に考えてみましょう。エネルギーの時間変化は、
$$ \frac{dH}{dt} = \frac{\partial H}{\partial \bm{q}} \cdot \bm{f}_q + \frac{\partial H}{\partial \bm{p}} \cdot \bm{f}_p $$
$\bm{f}_q$ と $\bm{f}_p$ が任意の関数の場合、この値は一般にゼロになりません。正準方程式のように $\bm{f}_q = \frac{\partial H}{\partial \bm{p}}$, $\bm{f}_p = -\frac{\partial H}{\partial \bm{q}}$ という特定の関係がある場合にのみ相殺が起こります。
HNNの構造
HNNの核心的なアイデアは、ベクトル場を直接学習する代わりに、ハミルトニアン関数 $H_{\bm{\theta}}(\bm{q}, \bm{p})$ をニューラルネットワークで学習し、正準方程式を通じて時間発展を導出することです。
$$ H_{\bm{\theta}}(\bm{q}, \bm{p}) = \text{NN}_{\bm{\theta}}(\bm{q}, \bm{p}) $$
時間発展は自動微分で計算した勾配から以下のように決まります。
$$ \frac{d\bm{q}}{dt} = \frac{\partial H_{\bm{\theta}}}{\partial \bm{p}}, \quad \frac{d\bm{p}}{dt} = -\frac{\partial H_{\bm{\theta}}}{\partial \bm{q}} $$
エネルギー保存の自動的保証
前節で証明したように、正準方程式に従う系では $\frac{dH}{dt} = 0$ が構造的に成立します。HNNは正準方程式の構造を保持しているため、学習されたハミルトニアン $H_{\bm{\theta}}$ がどのような関数であっても、$H_{\bm{\theta}}$ は時間とともに保存されます。
改めて確認しましょう。
$$ \begin{align} \frac{dH_{\bm{\theta}}}{dt} &= \frac{\partial H_{\bm{\theta}}}{\partial \bm{q}} \cdot \frac{d\bm{q}}{dt} + \frac{\partial H_{\bm{\theta}}}{\partial \bm{p}} \cdot \frac{d\bm{p}}{dt} \\ &= \frac{\partial H_{\bm{\theta}}}{\partial \bm{q}} \cdot \frac{\partial H_{\bm{\theta}}}{\partial \bm{p}} + \frac{\partial H_{\bm{\theta}}}{\partial \bm{p}} \cdot \left(-\frac{\partial H_{\bm{\theta}}}{\partial \bm{q}}\right) \\ &= \frac{\partial H_{\bm{\theta}}}{\partial \bm{q}} \cdot \frac{\partial H_{\bm{\theta}}}{\partial \bm{p}} – \frac{\partial H_{\bm{\theta}}}{\partial \bm{p}} \cdot \frac{\partial H_{\bm{\theta}}}{\partial \bm{q}} \\ &= 0 \end{align} $$
この証明は $H_{\bm{\theta}}$ の具体的な形に全く依存しません。ニューラルネットワークが $H_{\bm{\theta}}$ として何を学習しても、正準方程式の構造から自動的にエネルギーが保存されます。
重要な注意点として、ここでの「保存される量」は $H_{\bm{\theta}}$（学習されたハミルトニアン）であり、必ずしも真のハミルトニアン $H$ とは一致しません。しかし、学習が十分に進めば $H_{\bm{\theta}} \approx H$ となり、真のエネルギー保存も近似的に成立します。
学習の損失関数
HNNの学習には、状態 $(\bm{q}, \bm{p})$ と時間微分 $(\dot{\bm{q}}, \dot{\bm{p}})$ のペアからなる訓練データが必要です。
$$ \mathcal{D} = \left\{(\bm{q}^{(i)}, \bm{p}^{(i)}, \dot{\bm{q}}^{(i)}, \dot{\bm{p}}^{(i)})\right\}_{i=1}^{N} $$
損失関数は、正準方程式から導かれる時間微分と、訓練データの時間微分の二乗誤差です。
$$ \mathcal{L}(\bm{\theta}) = \frac{1}{N}\sum_{i=1}^{N}\left[\left\|\dot{\bm{q}}^{(i)} – \frac{\partial H_{\bm{\theta}}}{\partial \bm{p}}\bigg|_{(\bm{q}^{(i)}, \bm{p}^{(i)})}\right\|^2 + \left\|\dot{\bm{p}}^{(i)} + \frac{\partial H_{\bm{\theta}}}{\partial \bm{q}}\bigg|_{(\bm{q}^{(i)}, \bm{p}^{(i)})}\right\|^2\right] $$
第1項は $\dot{\bm{q}} = \frac{\partial H_{\bm{\theta}}}{\partial \bm{p}}$ の誤差、第2項は $\dot{\bm{p}} = -\frac{\partial H_{\bm{\theta}}}{\partial \bm{q}}$ の誤差です。
訓練データの生成
実際の実験では、以下の方法で訓練データを生成します。
- シミュレーション: 既知のハミルトニアンからODEソルバーで軌道を生成し、状態と時間微分を記録
- 有限差分: 軌道データ $(\bm{q}(t), \bm{p}(t))$ から $\dot{\bm{q}} \approx \frac{\bm{q}(t + \Delta t) – \bm{q}(t – \Delta t)}{2\Delta t}$ で近似
- 実験データ: センサーで測定された位置と速度のデータ
勾配の計算
損失関数の勾配 $\frac{\partial \mathcal{L}}{\partial \bm{\theta}}$ を計算するには、$\frac{\partial H_{\bm{\theta}}}{\partial \bm{q}}$ や $\frac{\partial H_{\bm{\theta}}}{\partial \bm{p}}$ のパラメータ $\bm{\theta}$ に関する勾配が必要です。これは2階の微分（$H_{\bm{\theta}}$ の入力に関する微分の、パラメータに関する微分）ですが、PyTorchの自動微分では create_graph=True を指定することで計算グラフを保持し、高階微分を自動的に処理できます。
シンプレクティック積分法
HNNで学習したハミルトニアンを使って長時間のシミュレーションを行う際、ODEソルバーの選択が重要です。通常のRunge-Kutta法は高精度ですが、ハミルトン系の幾何学的構造（シンプレクティック構造）を保存しません。そのため、長時間では微小なエネルギーのドリフトが蓄積します。
シンプレクティック構造
ハミルトン系の時間発展はシンプレクティック写像です。すなわち、位相空間の体積要素 $d\bm{q} \wedge d\bm{p}$（シンプレクティック2形式）を保存します。数学的には、状態 $\bm{z} = (\bm{q}, \bm{p})$ の時間発展 $\bm{z}(t) = \phi_t(\bm{z}(0))$ に対して、ヤコビアン $\bm{M} = \frac{\partial \phi_t}{\partial \bm{z}(0)}$ が以下を満たします。
$$ \bm{M}^\top \bm{J} \bm{M} = \bm{J} $$
ここで $\bm{J}$ はシンプレクティック行列です。
リープフロッグ法（Stormer-Verlet法）
最も広く使われるシンプレクティック積分法はリープフロッグ法（leapfrog method）です。ステップ幅 $\Delta t$ の更新式は以下です。
$$ \begin{align} \bm{p}_{n+1/2} &= \bm{p}_n – \frac{\Delta t}{2}\frac{\partial H}{\partial \bm{q}}\bigg|_{(\bm{q}_n, \bm{p}_{n+1/2})} \\ \bm{q}_{n+1} &= \bm{q}_n + \Delta t \frac{\partial H}{\partial \bm{p}}\bigg|_{(\bm{q}_n, \bm{p}_{n+1/2})} \\ \bm{p}_{n+1} &= \bm{p}_{n+1/2} – \frac{\Delta t}{2}\frac{\partial H}{\partial \bm{q}}\bigg|_{(\bm{q}_{n+1}, \bm{p}_{n+1/2})} \end{align} $$
ハミルトニアンが $H = T(\bm{p}) + V(\bm{q})$ と分離可能な場合（多くの力学系が該当）、$\frac{\partial H}{\partial \bm{q}} = \frac{\partial V}{\partial \bm{q}}$, $\frac{\partial H}{\partial \bm{p}} = \frac{\partial T}{\partial \bm{p}}$ となり、各ステップが陽的に計算できます。
$$ \begin{align} \bm{p}_{n+1/2} &= \bm{p}_n – \frac{\Delta t}{2}\nabla_{\bm{q}} V(\bm{q}_n) \\ \bm{q}_{n+1} &= \bm{q}_n + \Delta t \, \nabla_{\bm{p}} T(\bm{p}_{n+1/2}) \\ \bm{p}_{n+1} &= \bm{p}_{n+1/2} – \frac{\Delta t}{2}\nabla_{\bm{q}} V(\bm{q}_{n+1}) \end{align} $$
リープフロッグ法の重要な性質は以下です。
- シンプレクティック性: 位相空間のシンプレクティック構造を保存する
- 2次精度: 局所打ち切り誤差は $O(\Delta t^3)$、大域誤差は $O(\Delta t^2)$
- エネルギーの長期安定性: 真のハミルトニアン $H$ の周りで有界に振動し、系統的なドリフトが生じない
- 時間反転対称性: $\Delta t \to -\Delta t$ としても同じ方法が得られる
HNNの場合、$H_{\bm{\theta}}$ が分離可能な形でない場合もありますが、一般のシンプレクティック積分法（暗黙的中点法など）を適用するか、分離可能な構造を仮定してネットワークを設計する方法があります。
Pythonでの実装: 単振り子のHNN
ここからは、単振り子を例にHNNを実装し、通常のニューラルネットワークとエネルギー保存の精度を比較します。
問題設定: 単振り子
単振り子のハミルトニアンは以下です（質量 $m=1$, 長さ $l=1$, 重力加速度 $g=9.81$ とします）。
$$ H(q, p) = \frac{p^2}{2ml^2} + mgl(1 – \cos q) $$
ここで $q$ は角度、$p = ml^2 \dot{q}$ は角運動量です。正準方程式は以下になります。
$$ \dot{q} = \frac{\partial H}{\partial p} = \frac{p}{ml^2}, \quad \dot{p} = -\frac{\partial H}{\partial q} = -mgl\sin q $$
データの生成
import torch
import torch.nn as nn
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
# 再現性
torch.manual_seed(42)
np.random.seed(42)
device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')
# 物理パラメータ
m = 1.0 # 質量
l = 1.0 # 長さ
g = 9.81 # 重力加速度
def true_hamiltonian(q, p):
"""真のハミルトニアン"""
return p**2 / (2 * m * l**2) + m * g * l * (1 - np.cos(q))
def pendulum_dynamics(t, state):
"""単振り子のハミルトン方程式"""
q, p = state
dqdt = p / (m * l**2)
dpdt = -m * g * l * np.sin(q)
return [dqdt, dpdt]
def generate_pendulum_data(n_trajectories=30, t_span=(0, 5), dt=0.05):
"""訓練データの生成: 複数の初期条件から軌道を生成"""
t_eval = np.arange(t_span[0], t_span[1], dt)
all_states = []
all_derivatives = []
for _ in range(n_trajectories):
# ランダムな初期条件
q0 = np.random.uniform(-np.pi/2, np.pi/2)
p0 = np.random.uniform(-3, 3)
# ODEを解く
sol = solve_ivp(pendulum_dynamics, t_span, [q0, p0],
t_eval=t_eval, method='RK45', rtol=1e-10, atol=1e-10)
states = sol.y.T # (n_times, 2)
# 時間微分を解析的に計算
derivatives = np.array([pendulum_dynamics(0, s) for s in states])
all_states.append(states)
all_derivatives.append(derivatives)
states = np.vstack(all_states)
derivatives = np.vstack(all_derivatives)
return states, derivatives
# 訓練データの生成
states, derivatives = generate_pendulum_data(n_trajectories=30)
print(f'Training data: {states.shape[0]} samples')
# テンソルに変換
X_train = torch.tensor(states, dtype=torch.float32).to(device)
Y_train = torch.tensor(derivatives, dtype=torch.float32).to(device)
HNNモデルの定義
class HNN(nn.Module):
"""ハミルトニアンニューラルネットワーク"""
def __init__(self, hidden_dim=128):
super().__init__()
# ハミルトニアン H(q, p) をニューラルネットワークで表現
self.hamiltonian_net = nn.Sequential(
nn.Linear(2, hidden_dim),
nn.Tanh(),
nn.Linear(hidden_dim, hidden_dim),
nn.Tanh(),
nn.Linear(hidden_dim, hidden_dim),
nn.Tanh(),
nn.Linear(hidden_dim, 1),
)
# 重みの初期化
for mod in self.hamiltonian_net:
if isinstance(mod, nn.Linear):
nn.init.orthogonal_(mod.weight)
nn.init.zeros_(mod.bias)
def forward(self, qp):
"""ハミルトニアンの値を計算"""
return self.hamiltonian_net(qp)
def time_derivative(self, qp):
"""正準方程式から時間微分を計算"""
qp = qp.requires_grad_(True)
H = self.forward(qp)
# H の q, p に関する勾配
dH = torch.autograd.grad(
H, qp,
grad_outputs=torch.ones_like(H),
create_graph=True
)[0]
dH_dq = dH[:, 0:1] # ∂H/∂q
dH_dp = dH[:, 1:2] # ∂H/∂p
# 正準方程式: dq/dt = ∂H/∂p, dp/dt = -∂H/∂q
dqdt = dH_dp
dpdt = -dH_dq
return torch.cat([dqdt, dpdt], dim=1)
通常のNNモデル（比較用）
class BaselineNN(nn.Module):
"""通常のニューラルネットワーク（ベクトル場を直接学習）"""
def __init__(self, hidden_dim=128):
super().__init__()
self.net = nn.Sequential(
nn.Linear(2, hidden_dim),
nn.Tanh(),
nn.Linear(hidden_dim, hidden_dim),
nn.Tanh(),
nn.Linear(hidden_dim, hidden_dim),
nn.Tanh(),
nn.Linear(hidden_dim, 2), # 直接 (dq/dt, dp/dt) を出力
)
for mod in self.net:
if isinstance(mod, nn.Linear):
nn.init.orthogonal_(mod.weight)
nn.init.zeros_(mod.bias)
def time_derivative(self, qp):
"""時間微分を直接出力"""
return self.net(qp)
学習
from torch.utils.data import DataLoader, TensorDataset
# データローダー
dataset = TensorDataset(X_train, Y_train)
loader = DataLoader(dataset, batch_size=256, shuffle=True)
def train_model(model, loader, epochs=2000, lr=1e-3, model_name='Model'):
"""モデルの学習"""
optimizer = torch.optim.Adam(model.parameters(), lr=lr)
scheduler = torch.optim.lr_scheduler.CosineAnnealingLR(optimizer, T_max=epochs)
loss_fn = nn.MSELoss()
losses = []
for epoch in range(1, epochs + 1):
epoch_loss = 0
for batch_x, batch_y in loader:
optimizer.zero_grad()
pred = model.time_derivative(batch_x)
loss = loss_fn(pred, batch_y)
loss.backward()
torch.nn.utils.clip_grad_norm_(model.parameters(), max_norm=1.0)
optimizer.step()
epoch_loss += loss.item() * batch_x.size(0)
scheduler.step()
avg_loss = epoch_loss / len(loader.dataset)
losses.append(avg_loss)
if epoch % 500 == 0:
print(f'[{model_name}] Epoch {epoch}/{epochs}, Loss: {avg_loss:.8f}')
return losses
# HNNの学習
hnn = HNN(hidden_dim=128).to(device)
hnn_losses = train_model(hnn, loader, epochs=2000, lr=1e-3, model_name='HNN')
# 通常NNの学習
baseline = BaselineNN(hidden_dim=128).to(device)
baseline_losses = train_model(baseline, loader, epochs=2000, lr=1e-3, model_name='Baseline')
学習曲線の比較
plt.figure(figsize=(10, 5))
plt.semilogy(hnn_losses, label='HNN', linewidth=2)
plt.semilogy(baseline_losses, label='Baseline NN', linewidth=2)
plt.xlabel('Epoch', fontsize=12)
plt.ylabel('MSE Loss', fontsize=12)
plt.title('Training Loss Comparison', fontsize=14)
plt.legend(fontsize=12)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
シミュレーションとエネルギー保存の比較
def simulate(model, q0, p0, dt=0.01, n_steps=3000):
"""学習済みモデルでシミュレーション（オイラー法）"""
trajectory = [(q0, p0)]
q, p = q0, p0
for _ in range(n_steps):
qp = torch.tensor([[q, p]], dtype=torch.float32).to(device)
with torch.no_grad():
deriv = model.time_derivative(qp).cpu().numpy()[0]
dq, dp = deriv
q = q + dt * dq
p = p + dt * dp
trajectory.append((q, p))
return np.array(trajectory)
def simulate_leapfrog(model, q0, p0, dt=0.01, n_steps=3000):
"""HNN + リープフロッグ法でシミュレーション"""
trajectory = [(q0, p0)]
q, p = q0, p0
for _ in range(n_steps):
# 半ステップの運動量更新
qp = torch.tensor([[q, p]], dtype=torch.float32, requires_grad=True).to(device)
H = model(qp)
dH = torch.autograd.grad(H, qp, torch.ones_like(H))[0]
dH_dq = dH[0, 0].item()
p_half = p - 0.5 * dt * dH_dq
# 全ステップの座標更新
qp = torch.tensor([[q, p_half]], dtype=torch.float32, requires_grad=True).to(device)
H = model(qp)
dH = torch.autograd.grad(H, qp, torch.ones_like(H))[0]
dH_dp = dH[0, 1].item()
q_new = q + dt * dH_dp
# 半ステップの運動量更新
qp = torch.tensor([[q_new, p_half]], dtype=torch.float32, requires_grad=True).to(device)
H = model(qp)
dH = torch.autograd.grad(H, qp, torch.ones_like(H))[0]
dH_dq = dH[0, 0].item()
p_new = p_half - 0.5 * dt * dH_dq
q, p = q_new, p_new
trajectory.append((q, p))
return np.array(trajectory)
# 初期条件
q0, p0 = 1.0, 0.0
dt = 0.01
n_steps = 3000
t_sim = np.arange(n_steps + 1) * dt
# 真の軌道
sol_true = solve_ivp(pendulum_dynamics, (0, n_steps * dt), [q0, p0],
t_eval=t_sim, method='RK45', rtol=1e-12, atol=1e-12)
true_traj = sol_true.y.T
# 各モデルのシミュレーション
hnn_traj_euler = simulate(hnn, q0, p0, dt=dt, n_steps=n_steps)
baseline_traj = simulate(baseline, q0, p0, dt=dt, n_steps=n_steps)
hnn_traj_leapfrog = simulate_leapfrog(hnn, q0, p0, dt=dt, n_steps=n_steps)
軌道の比較
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# q(t) の比較
ax = axes[0, 0]
ax.plot(t_sim, true_traj[:, 0], 'k-', linewidth=2, label='True')
ax.plot(t_sim, hnn_traj_leapfrog[:, 0], 'b--', linewidth=1.5, label='HNN + Leapfrog')
ax.plot(t_sim, hnn_traj_euler[:, 0], 'g--', linewidth=1.5, label='HNN + Euler')
ax.plot(t_sim, baseline_traj[:, 0], 'r--', linewidth=1.5, label='Baseline NN')
ax.set_xlabel('Time', fontsize=12)
ax.set_ylabel('q (angle)', fontsize=12)
ax.set_title('Angle q(t)', fontsize=14)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
# p(t) の比較
ax = axes[0, 1]
ax.plot(t_sim, true_traj[:, 1], 'k-', linewidth=2, label='True')
ax.plot(t_sim, hnn_traj_leapfrog[:, 1], 'b--', linewidth=1.5, label='HNN + Leapfrog')
ax.plot(t_sim, hnn_traj_euler[:, 1], 'g--', linewidth=1.5, label='HNN + Euler')
ax.plot(t_sim, baseline_traj[:, 1], 'r--', linewidth=1.5, label='Baseline NN')
ax.set_xlabel('Time', fontsize=12)
ax.set_ylabel('p (momentum)', fontsize=12)
ax.set_title('Momentum p(t)', fontsize=14)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
# 位相空間の軌道
ax = axes[1, 0]
ax.plot(true_traj[:, 0], true_traj[:, 1], 'k-', linewidth=2, label='True')
ax.plot(hnn_traj_leapfrog[:, 0], hnn_traj_leapfrog[:, 1], 'b--',
linewidth=1.5, label='HNN + Leapfrog')
ax.plot(hnn_traj_euler[:, 0], hnn_traj_euler[:, 1], 'g--',
linewidth=1.5, label='HNN + Euler')
ax.plot(baseline_traj[:, 0], baseline_traj[:, 1], 'r--',
linewidth=1.5, label='Baseline NN')
ax.set_xlabel('q', fontsize=12)
ax.set_ylabel('p', fontsize=12)
ax.set_title('Phase Space', fontsize=14)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
ax.set_aspect('equal')
# エネルギーの時間発展
ax = axes[1, 1]
E_true = np.array([true_hamiltonian(s[0], s[1]) for s in true_traj])
E_hnn_leapfrog = np.array([true_hamiltonian(s[0], s[1]) for s in hnn_traj_leapfrog])
E_hnn_euler = np.array([true_hamiltonian(s[0], s[1]) for s in hnn_traj_euler])
E_baseline = np.array([true_hamiltonian(s[0], s[1]) for s in baseline_traj])
ax.plot(t_sim, E_true, 'k-', linewidth=2, label='True')
ax.plot(t_sim, E_hnn_leapfrog, 'b-', linewidth=1.5,
label='HNN + Leapfrog', alpha=0.8)
ax.plot(t_sim, E_hnn_euler, 'g-', linewidth=1.5,
label='HNN + Euler', alpha=0.8)
ax.plot(t_sim, E_baseline, 'r-', linewidth=1.5,
label='Baseline NN', alpha=0.8)
ax.set_xlabel('Time', fontsize=12)
ax.set_ylabel('Energy H(q, p)', fontsize=12)
ax.set_title('Energy Conservation', fontsize=14)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
plt.suptitle('Pendulum: HNN vs Baseline NN Comparison', fontsize=15)
plt.tight_layout()
plt.show()
エネルギー保存精度の定量比較
# エネルギー誤差の統計
E0 = E_true[0]
print("=" * 60)
print(f"{'Method':<25} {'Mean |ΔE|':>12} {'Max |ΔE|':>12} {'Std |ΔE|':>12}")
print("=" * 60)
for name, E in [('True (RK45)', E_true),
('HNN + Leapfrog', E_hnn_leapfrog),
('HNN + Euler', E_hnn_euler),
('Baseline NN', E_baseline)]:
dE = np.abs(E - E0)
print(f"{name:<25} {np.mean(dE):>12.6f} {np.max(dE):>12.6f} {np.std(dE):>12.6f}")
print("=" * 60)
# エネルギー誤差の時間発展
plt.figure(figsize=(10, 5))
plt.semilogy(t_sim, np.abs(E_hnn_leapfrog - E0) + 1e-16, 'b-',
linewidth=1.5, label='HNN + Leapfrog', alpha=0.8)
plt.semilogy(t_sim, np.abs(E_hnn_euler - E0) + 1e-16, 'g-',
linewidth=1.5, label='HNN + Euler', alpha=0.8)
plt.semilogy(t_sim, np.abs(E_baseline - E0) + 1e-16, 'r-',
linewidth=1.5, label='Baseline NN', alpha=0.8)
plt.xlabel('Time', fontsize=12)
plt.ylabel('|E(t) - E(0)|', fontsize=12)
plt.title('Energy Error Over Time', fontsize=14)
plt.legend(fontsize=12)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
学習されたハミルトニアンの可視化
# 学習されたHNNのハミルトニアンと真のハミルトニアンを比較
q_grid = np.linspace(-np.pi, np.pi, 100)
p_grid = np.linspace(-5, 5, 100)
Q_mesh, P_mesh = np.meshgrid(q_grid, p_grid)
# 真のハミルトニアン
H_true = true_hamiltonian(Q_mesh, P_mesh)
# HNNの予測
qp_flat = np.column_stack([Q_mesh.ravel(), P_mesh.ravel()])
qp_tensor = torch.tensor(qp_flat, dtype=torch.float32).to(device)
hnn.eval()
with torch.no_grad():
H_pred = hnn(qp_tensor).cpu().numpy().reshape(Q_mesh.shape)
# 可視化
fig, axes = plt.subplots(1, 3, figsize=(18, 5))
# 真のハミルトニアン
im0 = axes[0].contourf(Q_mesh, P_mesh, H_true, levels=30, cmap='viridis')
axes[0].set_xlabel('q', fontsize=12)
axes[0].set_ylabel('p', fontsize=12)
axes[0].set_title('True Hamiltonian $H(q,p)$', fontsize=14)
plt.colorbar(im0, ax=axes[0])
# 学習されたハミルトニアン
im1 = axes[1].contourf(Q_mesh, P_mesh, H_pred, levels=30, cmap='viridis')
axes[1].set_xlabel('q', fontsize=12)
axes[1].set_ylabel('p', fontsize=12)
axes[1].set_title('Learned Hamiltonian $H_\\theta(q,p)$', fontsize=14)
plt.colorbar(im1, ax=axes[1])
# 誤差（定数オフセットを補正）
offset = np.mean(H_pred) - np.mean(H_true)
error = np.abs(H_pred - offset - H_true)
im2 = axes[2].contourf(Q_mesh, P_mesh, error, levels=30, cmap='hot')
axes[2].set_xlabel('q', fontsize=12)
axes[2].set_ylabel('p', fontsize=12)
axes[2].set_title('|$H_\\theta - H_{true}$| (offset corrected)', fontsize=14)
plt.colorbar(im2, ax=axes[2])
plt.tight_layout()
plt.show()
ベクトル場の可視化
# HNNと通常NNのベクトル場を比較
q_vec = np.linspace(-2, 2, 20)
p_vec = np.linspace(-4, 4, 20)
Q_v, P_v = np.meshgrid(q_vec, p_vec)
# 真のベクトル場
DQ_true = P_v / (m * l**2)
DP_true = -m * g * l * np.sin(Q_v)
# HNNのベクトル場
qp_v = torch.tensor(
np.column_stack([Q_v.ravel(), P_v.ravel()]), dtype=torch.float32
).to(device)
hnn.eval()
with torch.no_grad():
hnn_deriv = hnn.time_derivative(qp_v).cpu().numpy()
DQ_hnn = hnn_deriv[:, 0].reshape(Q_v.shape)
DP_hnn = hnn_deriv[:, 1].reshape(Q_v.shape)
# 通常NNのベクトル場
baseline.eval()
with torch.no_grad():
base_deriv = baseline.time_derivative(qp_v).cpu().numpy()
DQ_base = base_deriv[:, 0].reshape(Q_v.shape)
DP_base = base_deriv[:, 1].reshape(Q_v.shape)
fig, axes = plt.subplots(1, 3, figsize=(18, 5))
for ax, DQ, DP, title in [
(axes[0], DQ_true, DP_true, 'True Vector Field'),
(axes[1], DQ_hnn, DP_hnn, 'HNN Vector Field'),
(axes[2], DQ_base, DP_base, 'Baseline NN Vector Field'),
]:
speed = np.sqrt(DQ**2 + DP**2)
ax.quiver(Q_v, P_v, DQ, DP, speed, cmap='coolwarm', alpha=0.8)
ax.set_xlabel('q', fontsize=12)
ax.set_ylabel('p', fontsize=12)
ax.set_title(title, fontsize=14)
ax.grid(True, alpha=0.3)
ax.set_aspect('equal')
plt.tight_layout()
plt.show()
まとめ
本記事では、ハミルトニアンニューラルネットワーク（HNN）の理論を数式の導出とともに解説し、Pythonで単振り子の実装と通常NNとの比較を行いました。
- ハミルトンの正準方程式 $\dot{q}_i = \frac{\partial H}{\partial p_i}$, $\dot{p}_i = -\frac{\partial H}{\partial q_i}$ は、$\frac{dH}{dt} = 0$（エネルギー保存）を構造的に保証する
- 通常のNNはベクトル場を直接学習するため、正準方程式の構造が保存されず、長時間シミュレーションでエネルギーがドリフトする
- HNNはハミルトニアン $H_{\bm{\theta}}(\bm{q}, \bm{p})$ を学習し、正準方程式を通じて時間発展を導出することで、$H_{\bm{\theta}}$ の保存を自動的に保証する
- 損失関数は $\|\dot{\bm{q}} – \frac{\partial H_{\bm{\theta}}}{\partial \bm{p}}\|^2 + \|\dot{\bm{p}} + \frac{\partial H_{\bm{\theta}}}{\partial \bm{q}}\|^2$ であり、自動微分で計算される
- シンプレクティック積分法（リープフロッグ法）と組み合わせることで、数値積分レベルでもエネルギーの長期安定性が得られる
- 単振り子の実装では、HNN + リープフロッグ法が通常NNと比較してエネルギー保存精度が大幅に優れていることを確認した
次のステップとして、以下の記事も参考にしてください。