最小作用の原理をTensorflowで確認

# 解いた式

ラグランジアン

$$
L = K - U
$$

を時間積分した

$$
S = \int L \mathrm{d}t
$$

を馬鹿正直に数値的に最小化したら、それっぽい軌跡が得られました。
２次元平面での重力ポテンシャルの公転運動です。

# お断り

今回も機械学習の話ではありません。

# もうちょいくやしく

数値的に解くので、位置 $\vec{r}(t)$ は時間について離散的な関数となります。$N$ ステップの軌跡だと思うと、軌跡はベクトルみたいに書けます。

$$
x = \Bigl(\vec{r}(0), \vec{r}(\Delta t), \cdots, \vec{r}(N\Delta t) \Bigr)^\mathrm{T}
$$

速度は単純な差分で書けます。
ある時間ステップ $i$ では、その前後の成分の差分で速度が計算できます。テイラー展開の1次までとる発想です。

```math
\vec{p}(i\Delta t) = \frac{1}{2\Delta t}\Bigl(\vec{r}\left((i + 1)\Delta t\right) - \vec{r}\left((i - 1)\Delta t\right)\Bigr)
```

$p$は普通は運動量に使う記号ですが、速度に使っています。これを行列表記にしますと

```math
v = \frac{1}{2\Delta t}\left(
\begin{matrix}
-2 & 2 & 0 & 0 & \cdots \\
-1 & 0 & 1 & 0 & \cdots \\
0 & -1 & 0 & 1 & \cdots \\
  &    &  \vdots & & \\
0 &  \cdots & -1 & 0 & 1  \\
0 &  \cdots & 0 & -2 & 2  \\

\end{matrix}
\right)
\left(
\begin{matrix}
\vec{r}(0) \\
\vec{r}(\Delta t) \\
\vec{r}(2\Delta t) \\
\vdots \\
\vec{r}((N-1)\Delta t) \\
\vec{r}(N\Delta t) \\
\end{matrix}
\right)
```

1番目と$N$ 番目の成分（つまり端、最初と最後）はとなりの成分との差とします。
差分を取る行列を $D = \frac{1}{2\Delta t}(\cdots) $としますと、

```math
v = Dx
```

になります。

以降、$x$, $v$の$i$ステップ目の成分は添え字で表記することにします。

```math
x_i = \vec{r}(i\Delta t)\\
v_i = \vec{p}(i\Delta t)
```

また運動の途中で質量 $m$ は一定だとします。
位置と速度が定義されたので、各時間ステップでの運動エネルギー

```math
K_i = \frac{1}{2}mv_i^2
```

とポテンシャル

```math
U_i = G\frac{mM}{\left(x_i - \vec{R}\right)^2}
```

が書けます。
ここで$G$は重力定数、$M$は公転中心の質量、$\vec{R}$は公転中心の座標とします。今回は$G=1$、$M=1$、$\vec{R}=0$としてしまいます。

最後に、各時間ステップごとのラグランジアン

```math
L_i = K_i - U_i
```

を求めた後、それを全ステップにわたって合計した作用

```math
S = \sum_i \Delta t L_i
```

が計算できます。

さて、最初の位置 $x_1$ と最後の位置 $x_N$ を固定し、その間の $x_i$ をいろいろ動かして、$S$を最小化させてみることを考えます。

# 最小作用の原理

変分的に書くと、

```math
\delta S = 0
```

となり、自然界はこれを満たすような位相空間中の軌跡を取るというのが、解析力学の教えるところです。
ここから運動方程式などのハミルトン力学が派生していきます。つまりこの原理は普通は理論の礎になるもので、そのまま解く対象にはならないなという印象です。

$S$ は普通は位置と共役運動量 $r$, $p$の関数なのでこのように変分的に書きますが、
今回の話の依存関係を追ってみますと、

![Screenshot_20190330_214439.png](https://qiita-image-store.s3.amazonaws.com/0/258019/fad6cc2a-222a-8573-63c4-b1751c31fddf.png)

というグラフになっています。ここで$M, m, G, D$は系が決まれば定数となるものです。となれば独立に変更できそうなのは$x$だけですね。そして$x$は２次元中の$N$個の点列だと思うと、$2(N-2)$ 個のパラメータがあるように見えます。すると、

```math
x(軌跡) = \mathop{\rm arg~min}\limits_x \left[S\right]
```

として、なんかオプティマイザがつかえそうですね。

Tensorflowの`trainer`の出番です。

# 実装開始

## ライブラリ

```python
import tensorflow as tf
import numpy as np
import matplotlib.pyplot as plt
```

### 注意

ここでの`tensorflow`のバージョンは`1.x`です。Eagerモードがデフォルトになる`2.x`では実装を変えなければなりません。


## 定数

```python
G = 1.0
m = 1.0
Mc = 1.0
Mm = Mc*m

fp = tf.float64
Nstep = 500
Ndim = 2
Time_all = 6.4
dt = Time_all/Nstep
```

`Time_all`は軌跡の始点と終点の間にかかった時間で、`dt`はそれを総ステップ数`Nstep`で割ったものになります。

## 位置

```python
r_init = tf.get_variable("position_init", shape=[1, Ndim], dtype=fp, trainable=False)
r_fin = tf.get_variable("position_final", shape=[1, Ndim], dtype=fp, trainable=False)
r_path = tf.get_variable("position_path", shape=[Nstep - 2, Ndim], dtype=fp)

r = tf.concat([r_init, r_path, r_fin], 0)
```

## 微分用行列

```python
d_mat = np.eye(Nstep, k=1)
d_mat = d_mat - d_mat.T
d_mat[0, 0] = -1
d_mat[-1, -1] = 1
d_mat = tf.constant(d_mat/2)/dt
```

上で出てきた式を頑張ってnumpy行列にし、最後に`tf.constant`に渡しています。

## 速度

```python
dr = tf.matmul(d_mat, r)
```

ここで、`r`の形状は`(N, 2)`で、`d_mat`は`(N, N)`なので、`dr`は`(N, 2)`の形状の`tf.Tensor`です。

## 運動エネルギー

```python
vr_i = tf.reduce_sum(dr**2, axis=1)
K_i = 0.5*vr_i*m
```

各成分を2乗して足したらベクトルの長さの2乗ですね。
`reduce_sum`に`axis`を渡しているので、`dr**2`の形状`(N, 2)`のうちの`2`のほうが総和の対象となります。そのため`vr_i`の形状は`(N, )`です。

## 位置エネルギー

```python
rr_i = tf.rsqrt(tf.reduce_sum(r**2, axis=1))
U_i = - G*Mm*rr_i
```

## ラグランジアン

```python
L_i = K_i - U_i
```

## 作用

```python
S = tf.reduce_sum(L_i)*dt
```

`dt`はなくてもいいんじゃないかな。

## 最小化

```python
optimizer = tf.train.AdamOptimizer(0.05)
minimize = optimizer.minimize(S)
```

困った時のAdam。結果のところにありますが、一応`GradientDescent`も試しています。

## 初期座標

### 代入演算

```python
r_path_0 = tf.placeholder(fp, shape=[Nstep - 2, Ndim])
r_init_0 = tf.placeholder(fp, shape=[1, Ndim])
r_fin_0 = tf.placeholder(fp, shape=[1, Ndim])
init_assign = [tf.assign(r_path, r_path_0), tf.assign(r_init, r_init_0), tf.assign(r_fin, r_fin_0)]
initializer = tf.global_variables_initializer()
```

`r`は最小化の対象化で`tf.Variable`です。これに所定の初期値を与える方法として、

* `tf.get_variables`の`initializer`オプションにて、`tf.constant(value=hogehoge)`とする。
* 受け皿の`tf.placeholder`を用意し、`tf.assign`で上書き

が考えられます。今回は`r_path`の定義が初期値を決める処理より先に来てしまったので前者が使いにくく、後者としました。

### 実際の値

```python
r1 = np.array([1.0, 0.0])
r2 = np.array([0.0, 1.0])
r3 = np.array([-1.0, 0.0])
r4 = np.array([0.0, -1.0])
r5 = np.array([1.0, 0.0])
rs = np.c_[r1, r2, r3, r4, r5].T
ts = np.arange(5)/4
ts_interp = np.arange(Nstep)/(Nstep - 1)
x_interp = np.interp(ts_interp, ts, rs[:, 0])
y_interp = np.interp(ts_interp, ts, rs[:, 1])
r0_np = np.c_[x_interp, y_interp].reshape(Nstep, Ndim)
```

$(1,0), (0, 1), (-1, 0), (0, -1), (1, 0)$ の順に原点周りを１周し、その間を適当に補間した感じです。つまり原点周りに菱型に回る軌跡を初期状態としました。
`rs`がその通って欲しいチェックポイント５点を並べたもので、`ts`がその時刻です。ここでは開始を$t=0$, 終点を$t=1$にしました。注意：ここでの$t$はMD中での$t$とは異なり、ただの`np.interp`用の横軸です。

## 実行

### 可視化のmatplotlibとその初期化

```python
fig = plt.figure()
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
lines1,  = ax1.plot(x_interp, y_interp)
lines2,  = ax2.plot(x_interp, y_interp)
ax2.set_xlim((-1.5, 1.5))
ax2.set_ylim((-1.5, 1.5))
ax1.set_xlim((-1.5, 1.5))
ax1.set_ylim((-1.5, 1.5))
```

`ax1`は左半分で座標を、`ax2`は右半分で速度をプロットします。
が、ここでは速度の計算がまだなされていないので、座標をプロットしてお茶を濁します。

### セッション

```python
with tf.Session() as sess:
    sess.run(initializer)
```

### 初期値代入

```python
    sess.run(init_assign, feed_dict={r_path_0: r0_np[1:-1, :], r_init_0: r0_np[0:1, :], r_fin_0: r0_np[-1:, :]})
```

### 最小化実行

```python
    step = 0
    while True:
        sess.run(minimize)
        step += 1
```

### 可視化

```python
        if step%2 == 0:
            r_t = sess.run(r)
            lines1.set_data(r_t[:, 0], r_t[:, 1])

            v_t = sess.run(dr)
            lines2.set_data(v_t[:, 0], v_t[:, 1])

            plt.pause(0.001)
```

２イテレーションごとに、matplotlibの`ax2`にセットされた値を`set_data`で更新していきます。
たまに初期ウィンドウから大きく外れるので、`set_xlim` `set_ylim`で適宜調整します。
終わるときは、端末をアクティブにしてCtrl+Cです。

# 結果

左が位置、右が速度の図です。

* Adam

学習率は $0.05$ です。

![adam_opt.gif](https://qiita-image-store.s3.amazonaws.com/0/258019/f9090a87-6500-b630-f3fa-2d9d6edac666.gif)

GIF動画１周が1000イテレーションです。
収束が近づいた時に暴れるのが難点ですね。真剣にやるなら、収束が進んだ段階で学習率を変えたり、オプティマイザを変えるなどの工夫がいることでしょう。

* GradientDecent

まず、Adamの時と同じ学習率 $0.05$ にすると発散しました。

![GD005.gif](https://qiita-image-store.s3.amazonaws.com/0/258019/acd7d01d-a0ce-3c10-06c3-bb4a9eaba86e.gif)

動画１周が２５０イテレーションくらい。

学習率を $0.01$にすると発散はしませんでしたが、非常に収束が遅くなりました。

![GD001.gif](https://qiita-image-store.s3.amazonaws.com/0/258019/b81f4feb-476c-0844-6e67-5941cedf8aff.gif)

動画一周で１０００イテレーション。フレームレートはAdamのものと同じです。

* おまけ

Adamで`Time_all = 10.0`としてみました。

![movie.gif](https://qiita-image-store.s3.amazonaws.com/0/258019/3cfecc2b-1938-1334-7c1e-fff80e89d125.gif)

1000イテレーションでは若干収束には届きませんでしたが、１周にかかる時間を増やしたのに合わせて先ほどの軌跡よりも原点から遠くに行って、速度も小さくなっているのがわかると思います。

# 余談

今回は経路の最初と最後を固定しました。
もしここも動くようにすると、運動エネルギーを小さく、位置エネルギーを大きくすれば作用も小さくなりますが、これは今回の系では無限遠で静止というのになってしまいそうですね。

まあ動かないまでもある程度の束縛（一周する＝最初と最後の位置が一緒である）とかを課したもとでの最小化とかも考えられますが、今回はこの辺にしておきます。

# 野望

経路積分