冬籠りの三ヶ月間どっぷり自動微分(Automatic Differentiation)に浸かりました。俺の一世代前の貯留層シミュレーターiSMARTは導入の容易な差分商で近似する数値微分法を採用しているが、演算に大幅な手間が掛かる上に最適な微小値\(\epsilon\)を見出すのが難しく誤った解に収束することも起こりうる状況にある。近年研究対象のシミュレーターでは自動微分を用いてヤコビアン行列の組込みの手間の短縮化を図っていることを知り自動微分の組込みを夢見ていた。野外活動が鈍る北陸の冬の時期、iSMARTに自動微分の実装に向けてチャレンジを開始。Three months this winter are dedicated to the automatic differentiation (AD). My primitive reservoir simulator iSMART adopts the numerical differentiation method which is easy to implement but may end up converging to a wrong solution due to the loss of the search for an optimal \(\epsilon\). I wish to incorporate an AD function into iSMART, once I learn that the latest research simulators tend to use the AD function instead of spending huge time on constructing analytical Jacobian matrices. While my outdoor activities is reduced during the winter season in Hokuriku, I shall have challenge to vigorously implement the AD function in iSMART.

自動微分法は単純な原理に基づくが実用的効果の高い技法である。合成関数の導関数を微分の連鎖律(chain rule)に基づいて計算する手法である為、原理は単純であるが、\(n\)変数関数の勾配の値(\(n\)) 個の偏微分の値)を計算するのに、関数自身を計算するために実行する演算回数の高々定数倍の演算回数しか要しない。

inhouse貯留層シミュレーターSMARTは、市販の貯留層シミュレーター同様、煩雑で誤謬の恐れのあるものの演算速度を追求する為、解析微分法に基づくヤコビアン行列を組込んでいる。今回、自動微分法のイロハから学びながらiSMARTに対して自動微分法の実装を目指す。iSMARTはFortranでコーディングしているが、自動微分の実装が容易なPythonを用いる。実装できる自信はないがやれるだけやってみる腹積もり。

Partial Derivatives of Elementary Operators

先ず、基礎的な演算子の一次偏導関数に付き整理しておく。

関数値及び導関数の算定に当たっては、ベクトル表示\(\vec{u}=(u, u’)\)を用いるのが一般的である。ここで、第一項は関数値\(f(\vec{x})\)、第二項は導関数\(f'(\vec{x})\)を示す。

Arithmetic Operators

$$\begin{align}
\vec{u} \pm \vec{v} &= (u±v, u’±v’) \\
\vec{u}\times\vec{v} &= (uv, u’v+uv’) \\
\vec{u}\div\vec{v} &= (u/v, (u’-(u/v)v’)/v)　where　v\neq{0} \\
\end{align} \tag{A-1}$$

Multivariate Functions

一変数の合成関数に於いては連鎖律を用いて拡張定義できる事は周知の事実。

$$\vec{f}(\vec{u}) = \vec{f}(u, u’ ) = (f(u), u’f'(u)) \tag{A-2}$$

上式に基づき下記演算子の一次偏導関数は次式で与えられる。

$$\begin{align}
\sin \vec{u} &=\sin (u,u’) = (\sin u, u’\cos u) \\
\cos \vec{u} &=\cos (u,u’) = (\cos u, −u’\sin u) \\
\tan \vec{u} &=\tan (u,u’) = (\tan u, u’/cos^2u) \\
e^{\vec{u}} &=e^{(u,u’)} = (e^{u} , u’e^{u}) \\
\log \vec{u} &=\log (u,u’) = (\log u, u’/u)
\end{align} \tag{A-3}$$

ここで、上式\((A-2)\)に於いて二変数以上の合成関数とその偏導関数の自動微分の定義について振り返る。

二変数の偏導関数の自動微分についても、一変数\(u\)を\(u=(x,y)\)として連鎖律を用いて、元の集合\(f:\mathbb {R^1} \to \mathbb{R^2}\)に拡張したものと考えられるので、一般化すると、新たな集合\(\mathbb{R}^n\)に属する\(n\)変数についても同様に自動微分が成立することを理解できよう。

Toward Implementation

出来合いのオープンソースライブラリーも多数あるが、何せ時間はたっぷりあるので自前でコーディングすることにする。

Fortranによる実装が理想であるが実例は少なく、又Fortranの最新機能にも疎い俺には無理。C++、Pythonが選択肢となるが、C++はEasy-Readableな言語でない故触手がわかない。結局Pythonにならざるを得ない。

演算子のオーバーロード機能を利用して実装する。演算子のオーバーロード機能とは、「ユーザー定義型で演算や比較を行えるようにするための機能」で、ユーザー定義型でも同じような振る舞いを持たせることが可能となる。基本演算子には次の特殊メソッドが用意されている。これらを用いてクラスに自動微分を組込むことができる。

Arithmetic Operator: Left Operand

算術演算子:左辺の特殊メソッドは計算時にself(自身)が左辺の場合呼び出される。A+Bの時は、Aが左辺でBが右辺となる。

Operator	Function	Expression
+	足し算	__add__(self, other)
–	引き算	__sub__(self, other)
*	掛け算	__mul__(self, other)
/	割り算	__truediv__(self, other)
**	べき乗	__pow__(self, other)

Arithmetic Operator: Right Operand

計算時にselfが右辺の場合に呼び出される特殊メソッドで、特殊メソッドが左辺と右辺のどちらも定義されている場合は、左辺の特殊メソッドが優先される。

Operation	Function	Expression
+	足し算	__radd__(self, other)
–	引き算	__rsub__(self, other)
*	掛け算	__rmul__(self, other)
/	割り算	__rtruediv__(self, other)
**	べき乗	_rpow__(self, other)

Relational Operator

比較演算子用の特殊メソッドによって比較演算子を使えるようになる。

Operator	Function	Expression
<	小なり	__lt__(self, other)
<=	小なりイコール	__le__(self, other)
>	大なり	__gt__(self, other)
>=	大なりイコール	__ge__(self, other)
==	等しい	__eq__(self, other)
!=	等しくない	__ne__(self, other)

Partial Differentiation by Hand Calculation

手計算で自動微分を実行してみる。手計算で出来ればコンピュータに実行させる過程を直ぐに理解できよう。
変数値と変数の偏導分値のペアを(\(x,1)\)で示す。\(“1″\)は\(x\)の偏微分値が\(“1″\)であることによる。尚、定数\(a\)に対しては(\(a,0)\)となることはお分かり頂けよう。このペアを使って忍耐のいる作業をコンピュータで淡々と演算させるだけだ。これが先に触れた「単純な原理に基づく」に当たるが、思いついた御仁は偉い！

Four Elementary Arithmetic Operators

以下の4つの課題で練習してみよう。

$$\begin{align}
f(x) &= (x+2)(x-1) \\
\vec{f}(\vec{x}) &=(\vec{x}+\vec{2}) \times (\vec{x}-\vec{1}) \\
&= ((x,1)+(2,0)) \times ((x,1)-(1,0)) \\
f(3,1) &= ((3,1)+(2,0)) \times ((3,1)-(1,0)) = (5,1) \times (2,1) \\
&=(5×2,(5,1)'(2,1)+(2,1)'(5,1))=(10,1\times(2,1)+1\times(5,1)) \\
&=(10,2+5)=(10,7) //
\end{align}$$

$$\begin{align}
f(x) &= \frac{x+2}{x+1} \\
\vec{f}(\vec{x}) &=\frac{\vec{x}+\vec{2}}{\vec{x}+\vec{1}} \\
&= \frac{(x,1)+(2,0)}{(x,1)+(1,0)} \\
f(3,1) &=\frac{(3,1)+(2,0)}{(3,1)+(1,0)} = \frac{(5,1)}{(4,1)} \\
&=( \frac{5}{4}, \frac{(5,1)’-((5,1)/(4,1))(4,1)’}{(4,1)})=( \frac{5}{4}, \frac{(5,1)'(4,1)-(5,1)(4,1)’}{(4,1)^{2}}) \\
&=(\frac{5}{4}, \frac{1\times(4,1)-(5,1)\times1}{(4,1)^{2}})=(\frac{5}{4}, \frac{4-5}{16}) =(\frac{5}{4}, \frac{-1}{16}) = (1.25,-0.0625) //
\end{align}$$

$$\begin{align}
f(x) &= \frac{(x+2)(x-1)}{x+1} \\
\vec{f}(\vec{x}) &= \frac{(\vec{x}+\vec{2})(\vec{x}-\vec{1})}{\vec{x}+\vec{1}} \\
&= \frac{ ((x,1)+(2,0)) \times ((x,1)-(1,0)) }{(x,1)+(1,0)} \\
f(3,1) &= \frac{ ((3,1)+(2,0)) \times ((3,1)-(1,0)) }{(3,1)+(1,0)} \\
&= \frac{(5,1) \times (2,1)} {(4,1)} = \frac{(10,7)}{(4,1)} \\
&=( \frac{5}{2}, \frac{(10,7)’-((10,7)/(4,1))(4,1)’}{(4,1)})=( \frac{5}{2}, \frac{(10,7)'(4,1)-(10,7)(4,1)’}{(4,1)^{2}}) \\
&= (\frac{5}{2},\frac{(7\times(4,1)-(10,7)\times1}{(4,1)^{2}})=(\frac{5}{2},\frac{28-10}{16})=(\frac{5}{2},\frac{9}{8}) = (2.5,1.125) //
\end{align}$$

Elementary Composite Function

$$\begin{align}
f(x) &= \frac{\sin x}{1+e^{-x}} \\
\vec{f}(\vec{x}) &=\frac{\sin \vec{x}}{1+e^{-\vec{x}}} \\
&=( \frac{\sin (3,0)}{(1,0)+e^{-(3,0)}}, \frac{ (3,1)’\sin (3,1)’ – \sin (3,1)/((1,0)+e^{-(3,0)})((1,0)’+(e^{-(3,1)})’ ) }{(1,0)+e^{-(3,0)}}) \\
&=( \frac{\sin (3)}{1+e^{-3}}, \frac{ 1\times \cos (3) – \sin (3)/(1+e^{-3})(-1\times e^{-3} ) }{1+e^{-3}}) \\
&=( \frac{\sin (3)}{1+e^{-3}}, \frac{ \cos (3)}{1+e^{-3}}+\frac{e^{-3} \sin (3)}{(1+e^{-3})^{2}}) \\
&= (0.134, -0.937) //
\end{align}$$

お分かり頂けたかな。俺はこれを理解するのに多くの時間が掛かりました。

Closing

今回は手計算によって自動微分を実行する醍醐味？を味わった処で終わりとします。
次回は実装に向けたクラスの構築並びに実証テストにつき紹介する予定。単なるプログラミング知見を駆使するだけでこれ以上地頭は要りません。

さて、北陸にも春がやって来ました。

「開花待つ寒梅の鉢　紅椿花一輪で傍らともす」

References

最初の頃は自動微分の理論は取っ付きにくく難儀したが、手計算で算出できるようになってから一気に理解が進むようになりました。参考資料は山ほどありますが最も参考になった資料を下記に列挙します。特に、資料[2]はMATLAB向けなるも実装の事例もあって大変参考になる。

1. コラム：自動微分 — DeZero Book
2. Neidinger, R. 2010. Introduction to automatic differentiation and MATLAB object-oriented programming. SIAM Review, 52(3), 545–563. doi:10.1137/080743627
3. An Introduction to Reservoir Simulation Using MATLAB/GNU Octave: User Guide for the MATLAB Reservoir Simulation Toolbox (MRST), Cambridge University Press, c2019

〆