ラグランジアンのこころ

この記事は物理学におけるハミルトンの原理、すなわちラグランジアンやハミルトンを用いた変分原理（停留条件）についてのイメージと、数値計算による簡易なシミュレーション方法の覚書である。とくに$L=K-U$の直感的な理解を示す。また物理情報機械学習との関連についてコメントする。

ハミルトンの原理のイメージ

ラグランジュ形式

ラグランジュ形式では軌道 $q(t)$ とポテンシャルエネルギー $U(q)$ について、作用汎関数が以下で定義される。また（拘束条件等を考えなければ）停留条件のもとでそれぞれの項の大小が調整される。

$$S[q] = \underbrace{\int_{t_0}^{t_1} \mathrm dt \left(\frac{1}{2}m\dot q^2\right)}_{小さくする} - \underbrace{\int_{t_0}^{t_1} \mathrm dt \, U(q)}_{大きくする}$$

ここで運動エネルギーの時間積分は線素の長さと理解できる。まず、$x(s), y(s)$ で与えられる曲線について、$\mathrm ds$あたりの微小線素の長さ $\mathrm \ell(s)$ は

$$\ell(s)^2 = \left(\frac{\mathrm dx}{\mathrm ds}\right)^2 + \left(\frac{\mathrm dy}{\mathrm ds}\right)^2$$

である。ここで $x=s=t$ また $y=q$ とすると、

$$\ell(s)^2 = 1 + \left(\frac{\mathrm dq}{\mathrm dt}\right)^2 = 1 + \dot q^2$$

となる。したがって運動エネルギーの時間積分

$$\int_{t_0}^{t_1} \mathrm dt \left(\frac{1}{2}m\dot q^2\right)$$

を最小化するということは、幾何学的には経路がたどった時空間上の（2乗）距離を最小化することに対応している。自由粒子、すなわちポテンシャルエネルギーが０の場合は2点 $q(t_0), q(t_1)$ の最短経路、つまり等速直線運動をする。

一方でポテンシャルエネルギーが存在する場合、その積分

$$\int_{t_0}^{t_1} \mathrm dt \, U(q)$$

を最大化するように軌道が補正される。重力ポテンシャル $U(q) = mgq$ の場合は以下のようになる。

いわゆる $L=K-U$ はこのようにイメージできる。

ハミルトン形式

ハミルトン形式の場合はラグランジュ形式とは異なり、停留点が極小値ではないことに注意が必要である。とくに運動量については極大を選択する必要があり、全体としては鞍点が選ばれる。したがってラグランジュ形式のように項毎の最大化・最小化でイメージするのは難しい。しいて書けば次のようになる。

$$S[q,p] = \underbrace{\int_{q_0}^{q_1} p\mathrm d q}_{\substack{pで大きくする\\ qで小さくする}} - \underbrace{\int_{t_0}^{t_1} \mathrm dt\left(\frac{p^2}{2m}\right)}_{pで小さくする} - \underbrace{\int_{t_0}^{t_1} \mathrm dt \, U(q)}_{qで大きくする}$$

$p$ については停留条件が極大値を与えることについて、簡単な証左を見てみよう。自由粒子において $p$ が一定であるとすると

$$\begin{split} S[q,p] &= \int_{q_0}^{q_1} p\mathrm d q - \int_{t_0}^{t_1} \mathrm dt\left(\frac{p^2}{2m}\right)\\ & = pL - \frac{p^2}{2m} T \\ &= - \frac{T}{2m} \left (p - \frac{mL}{T}\right)^2 + \cdots \end{split}$$

したがって $p = mL/T$ が作用の極大値を与える。ここで $L = q_1 - q_0$, $T = t_1 - t_0$とおいた。

ハミルトンの原理のシミュレーション

ハミルトンの原理は求める軌道に仮設(ansatz)を与えることで容易にシミュレーションできる。ここではポテンシャルエネルギーとして重力

$$U(q) = mgq$$

を考えよう。そして軌道の候補として3次関数を考えてみる。すなわち

$$q(t) = c_3t^3 + c_2t^2 + c_1t + c_0$$

とおく。変分原理では端点を固定するため、固定された端点を通るような係数を次で与える。

$$\left( \begin{matrix} c_3 \\ c_2 \\ c_1 \\ c_0 \end{matrix} \right) = \left( \begin{matrix} t_0{}^3 & t_0{}^2 & t_0 & 1 \\ t_1{}^3 & t_1{}^2 & t_1 & 1 \\ t_a{}^3 & t_a{}^2 & t_a & 1 \\ t_b{}^3 & t_b{}^2 & t_b & 1 \\ \end{matrix} \right)^{-1} \left( \begin{matrix} q_0 \\ q_1 \\ q_a \\ q_b \end{matrix} \right)$$

ここで $t_0, t_1, q_0, q_1$ は変分の際に固定される点であり、$t_a, t_b, q_a, q_b$ が経路を決めるための自由パラメータになる。

さて、このように経路に仮設を与えてしまえばあとは Sympy などを用いることで作用汎関数

$$S[q] = \int_{t_0}^{t_1} \mathrm dt \left(\frac{1}{2}m\dot q\right) - \int_{t_0}^{t_1} \mathrm dt \left(mgq\right)$$

を陽に書き下せる。停留点の導出には

$$\mathrm{loss}(t_a,t_b,q_a,q_b) = \left(\frac{\partial S}{\partial t_a}\right)^2 + \left(\frac{\partial S}{\partial t_b}\right)^2 + \left(\frac{\partial S}{\partial q_a}\right)^2 + \left(\frac{\partial S}{\partial q_b}\right)^2 = 0$$

を勾配法等で解けばよい。下はこの考え方で鉛直投げ上げの軌道を求めた際の最適化プロセスである。点線は厳密解を表す。

ハミルトン形式

$$S[q,p] = \int_{q_0}^{q_1} p\mathrm d q - \int_{t_0}^{t_1} \mathrm dt\left(\frac{p^2}{2m}\right) - \int_{t_0}^{t_1} \mathrm dt \left(mgq\right)$$

の場合も同様だが、運動量 $p$ についても $q$ と同様に ansatz を置く必要があるため変数が増える。結果例は以下のとおりである。ラベルのaction は $\int p\mathrm dq$ 項である。

コメント: 物理情報機械学習との関連

ここでの発想は以下のサイクルを回すことであった。

軌道の生成にはとうぜんニューラルネットワークやガウス過程を用いてもよい。積分は数値的に評価せざるを得ないが、運動エネルギーに含まれる微分の計算等には自動微分が使える。loss 関数の計算にはデータとの残差を含めることもでき、その場合、このアルゴリズムは物理情報機械学習と呼べるものになる。ふつう物理情報機械学習では運動方程式を評価に用い、わざわざ変分原理を使うことはない。一方で運動方程式が煩雑なときなど、変分原理によるアプローチが有効な場合もあるかもしれない。また終点を固定できることはpath planning など制御理論における応用がありそうである。

コメント：自由エネルギーとの関連

ラグランジュ形式の式

$$S[q] = \underbrace{\int_{t_0}^{t_1} \mathrm dt \left(\frac{1}{2}m\dot q^2\right)}_{小さくする} - \underbrace{\int_{t_0}^{t_1} \mathrm dt \, U(q)}_{大きくする}$$

から、変分ベイズ法の式や、それに基づく自由エネルギーの式

$$F_{\mathrm{free\text{-}energy}} = -L_{\mathrm{ELBO}} = \underbrace{-\log p}_{\mathrm{suprize}} + \underbrace{\mathbb{KL}[q\parallel p]}_{\mathrm{divergence}}$$

を彷彿とした人もあるかもしれない。これらの背景を「ベイズ力学」という概念で理解・統一しようとしている人たちもいる。