概要
機械学習を調べていると自動微分という技術があることを知りました。内容をよく理解するために、自動微分を実装してみたいと思います。最後に、自動微分を使った減衰振動の微分を求めてみようと思います。
自動微分とは
任意の関数$f$に引数$x$を代入して戻り値を得るとします。
y=f(x)
この時、同時に以下のような関数$f$の$x$に対する微分を求めるのが自動微分です。
g=\frac{\partial f(x)}{\partial x}
なお、$y$も勾配$g$も具体的な数値$x$に対する数値であり、抽象的な数式を得られる訳ではありません。
実装上のルール
自動微分を実装する1つの方法として、あらゆる演算において、オリジナル値の計算と微分値を同時に行うという方法があります。
実装上では、あらゆる値は$\mathbb R^2$を、あらゆる関数は$\mathbb R^2\to\mathbb R^2$の形を取ることとし、要素0番目に0階微分のオリジナル値を、要素1番目に1階微分した勾配値を入れます。例えば
Z = f(Y)
とすると、
-
Y[0]$=y$ ←オリジナル値 -
Y[1]$=\partial y/\partial x$ ←微分値 -
Z[0]$=f(y)$ ←オリジナル値 -
Z[1]$=\partial f(y)/\partial x$ ←微分値
となるように頑張って実装します。また、変数名は大文字とします。
実装Level1
微分対象にある変数
変数$x$のオリジナル値X[0]は任意の数でいいですが、微分値X[1]$=\partial x/\partial x$は常に$1$です。
X = (3, 1)
微分対象にない定数
定数$a$のオリジナル値A[0]は任意の数でいいですが、微分値A[1]$=\partial a/\partial x$は常に$0$です。
A = (5, 0)
B = (7, 0)
足し算
$i+j$を行う関数add(I,J)を考えます。オリジナル値add(I,J)[0]は$i+j$ですが、微分値add(I,J)[1]は
\frac{\partial(i+j)}{\partial x}=\frac{\partial i}{\partial x}+\frac{\partial j}{\partial x}
となります。
def add(I, J): return I[0]+J[0], I[1]+J[1]
引き算
$i-j$を行う関数sub(I,J)を考えます。オリジナル値sub(I,J)[0]は$i-j$ですが、微分値sub(I,J)[1]は
\frac{\partial(i-j)}{\partial x}=\frac{\partial i}{\partial x}-\frac{\partial j}{\partial x}
となります。
def sub(I, J): return I[0]-J[0], I[1]-J[1]
掛け算
$ij$を行う関数mul(I,J)を考えます。オリジナル値mul(I,J)[0]は$ij$ですが、微分値mul(I,J)[0]は
\frac{\partial(ij)}{\partial x}=i\frac{\partial j}{\partial x}+j\frac{\partial i}{\partial x}
となります。
def mul(I, J):return I[0]*J[0], I[0]*J[1] + J[0]*I[1]
割り算
$i/j$を行う関数div(I,J)を考えます。微分値は
\frac{\partial(ij^{-1})}{\partial x} = j^{-2}(j\frac{\partial i}{\partial x}-i\frac{\partial j}{\partial x})
となります。
def div(I, J): return I[0]/J[0], J[0]**(-2) * (J[0]*I[1] - I[0]*J[1])
テストLevel1
1
$y=ax+b\ (a=5,x=3,b=7)$を考えます。解析的には
y=5\times 3+7=22\\
\frac{\partial y}{\partial x}=a=5
となりますが、ADでも
Y = mul(A,X)
Y # (22, 5)
と正しく求まります。
2
$y=a/x\ (a=5,x=3)$を考えます。解析的には
y=5/3\approx 1.67\\
\frac{\partial y}{\partial x}=a(-x^{-2})\approx -0.56
となりますが、Pythonでも
Y = div(A,X)
Y # (1.6666666666666667, -0.5555555555555556)
と正しく求まります。
実装Level2
少し高難易度な関数等を実装します。
sin
$\sin(i)$の微分は合成関数の微分より
\frac{\partial \sin(i)}{\partial x}=\frac{\partial i}{\partial x}\cos(i)
となるので
def sin(I): return math.sin(I[0]), I[1]*math.cos(I[0])
cos
$\cos(i)$の微分は合成関数の微分より
\frac{\partial \cos(i)}{\partial x}=-\frac{\partial i}{\partial x}\sin(i)
となるので
def cos(I): return math.cos(I[0]), - I[1]*math.sin(I[0])
exp
$\exp(i)$の微分は合成関数の微分より
\frac{\partial \exp(i)}{\partial x}=\frac{\partial i}{\partial x}\exp(i)
となるので
def exp(I): return math.exp(I[0]), I[1]*math.exp(I[0])
テスト - 減衰振動
最後に、減衰振動の式を組み、複雑な式でもちゃんと自動微分できるかを検証してみます。
y(t)=ce^{at}\sin(\omega t+\gamma)\\
\dot y(t)=c[ae^{at}\sin(\omega t+\gamma)+\omega e^{at}\cos(\omega t+\gamma)]\\
c=5, a=-0.4, \omega=3, \gamma=1
今回の微分対象にある変数は$t$となります。また、微分によって求まる$\dot y(t)$は減衰振動の速度にあたります。
まずは、様々な$t$における関数の値$y(t)$及び微分$\dot y(t)$を求めるメソッドを組みます。
def get_Y(t):
# 変数を作る
T = (t, 1)
# 定数を作る
C = (5, 0)
A = (-0.4, 0)
OMEGA = (3, 0)
GAMMA = (1, 0)
Y = \
mul(
C,
mul(
exp(mul(A, T)),
sin(add(mul(OMEGA, T), GAMMA))
)
)
return Y
次に、$y(t)$及び$\dot y(t)$の解析解(Calc)及び自動微分結果(AD)をそれぞれ以下のように作成します。
import numpy as np
import matplotlib.pyplot as plt
# 解析解を作成する
c, a, omega, gamma = 5, -0.4, 3, 1
t_s = np.linspace(0, 5)
calc_y_s = c*np.exp(a*t_s)*np.sin(omega*t_s+gamma)
calc_doty_s = c*(
a*np.exp(a*t_s)*np.sin(omega*t_s+gamma) +
omega*np.exp(a*t_s)*np.cos(omega*t_s+gamma)
)
# 自動微分結果を作成する
ad_y_s = []
ad_doty_s = []
for t in t_s:
Y = get_Y(t)
ad_y_s.append(Y[0])
ad_doty_s.append(Y[1])
最後に両者を比較します。
まずは$y(t)$です。これは微分をする前の数値なので当たり前ですが、両者は一致しています。
次は$y(t)$を微分した$\dot y(t)$です。
自動微分の結果が、解析解と一致しました🎉🎉🎉
私の実装した自動微分のコード達が正しく動いていることを確認できました。
さいごに
これで勾配法とか実装してみようかな