《高速微分法》のソースを表示

'''【こうそくびぶんほう (fast differentiation)】'''

　非線形関数の[[勾配]], [[ヤコビ行列]], [[ヘッセ行列]]等の値を数値的に計算する方法のひとつ. 高速自動微分法(fast automatic differentiation), 計算微分法(computational differentiation), 単純に自動微分(automatic differentiation; 以下 AD)ともいう. 主なアルゴリズムは2種あり, ボトムアップ(前進)自動微分(bottom-up AD, forward AD; 以下 BUAD) と, トップダウン(逆行)自動微分(top-down AD, reverse AD, backward AD; 以下 TDAD) という [1, 2]. 高速微分法は狭義には, TDADを指す. AD は「関数の値を計算するプログラム」から「偏導関数の値を計算するプログラム」を生成する手順を与え, 生成物を(コンパイルし)実行すれば, 差分商近似のような打ち切り誤差無しで, 正確な偏導関数の値を計算できる. 大規模システムの数学モデル等の大規模プログラム(数千行以上)により表現される関数の偏導関数を計算できるのが特長. <math>n\,</math>変数関数の勾配の<math>n\,</math>個の値を関数計算の手間の定数倍で計算できる点が「高速」微分の由来である. 

　以下，BUAD と TDAD による計算法を説明する．例として，2変数関数 <math>f(x,y)=x/\sqrt{x^2+y}</math> について, <math>f(3,4)\,</math> の値を計算する代入文の列(プログラム), <math>x=3, y=4, v_1=x, v_2=y, v_3=v_1*v_1, v_4=v_3+v_2, v_5=\sqrt{v_4}, v_6=v_1/v_5</math> を考えよう. ただし, 各代入文の右辺には, 演算(基本演算とよぶ)が高々1回だけ現れるとする. <math>v_1\,</math>, <math>v_2\,</math> が <math>x\,</math>, <math>y\,</math> に対応し, <math>v_6\,</math> に <math>f(x,y)\,</math> の値が計算される. 一般には, <math>n\,</math>変数関数 <math>f(x_1,\cdots,x_n)</math> について, <math>k\,</math>回目の代入文には, <math>k-1\,</math>回目までに計算される変数が現れうるから, 延べ <math>r\,</math> 回の演算を行なう代入文の列は<math>\{v_k=\varphi_k(v_1,\cdots,v_{k-1})\}_{k=1}^r</math>と表される. これを計算過程といい, <math>v_k\,</math> を中間変数という. <math>k\leq n</math> のとき <math>\varphi_k</math> は<math>v_k=x_k</math> という入力定数の代入演算に相当する.

　BUADは, 補助変数 <math>\{s_k\}_{k=1}^r</math> を導入し, 任意に <math>j\,</math> <math>(1\leq j\leq n)</math>を固定して, 合成関数の<math>x_j\,</math> に関する偏微分則 <math>\textstyle {\partial v_k}/{\partial x_j} = \sum_{i=1}^{k-1}({\partial \varphi_k}/{\partial v_i})\cdot({\partial v_i}/{\partial x_j})</math> に基づき, <math>s_k\,</math> を計算する式を導出する. 基本演算 <math>\varphi_k</math> を四則演算や初等関数などの2項・単項の演算に限れば, 表1により, <math>{\partial \varphi_k}/{\partial v_i}</math> (これを要素的偏導関数という)を導出できる. <math>s_j=1\,</math>, <math>s_\ell=0</math> <math>(1\leq\ell\leq n, \ell\not=j)</math> と初期設定すれば, <math>k=n+1\,, n+2\,, \cdots</math> について<math>s_i=\partial v_i/\partial x_j</math> <math>(i=1,\cdots,k-1)</math>を計算済みとみなすことができ, <math>\textstyle s_k=\sum_{i=1}^{k-1}({\partial \varphi_k}/{\partial v_i})\cdot s_i</math> の値を計算できる. 最終的に <math>s_r=\partial f/\partial x_j</math> となる. 

<div align="center">
<table width="600" border="0" cellspacing="0" cellpadding="0">
<tr>
<td colspan="2" align="center">表１：基本演算と要素的偏導関数<br>
</td>
</tr>
<tr>
<td align="center"><table width="200" border="1" cellspacing="2" cellpadding="2">
<tr>
<td align="center" valign="middle"><math>\varphi_k</math></td>
<td align="center" valign="middle"><math>\partial \varphi_k/ v_\alpha</math></td>
<td align="center" valign="middle"><math>\partial \varphi_k/ v_\beta</math></td>
</tr>
<tr>
<td align="center" valign="middle"><math>v_k=v_\alpha \pm v_\beta\, </math></td>
<td align="center" valign="middle"><math>1\, </math></td>
<td align="center" valign="middle"><math>\pm1</math></td>
</tr>
<tr>
<td align="center" valign="middle"><math>v_k=v_\alpha * v_\beta\, </math></td>
<td align="center" valign="middle"><math>v_\beta\, </math></td>
<td align="center" valign="middle"><math>v_\alpha\, </math></td>
</tr>
<tr>
<td align="center" valign="middle"><math>v_k=v_\alpha / v_\beta\, </math></td>
<td align="center" valign="middle"><math>1/v_\beta\, </math></td>
<td align="center" valign="middle"><math>-v_\alpha/({v_\beta}^2)\, </math><br>
<p><math>(=-v_k/v_\beta)\, </math><br>
</p>
</td></tr></table>
</td>
<td align="center"><table width="200" border="1" cellspacing="2" cellpadding="2">
<tr>
<td align="center" valign="middle"><math>\varphi_k\, </math></td>
<td align="center" valign="middle"><math>\partial \varphi_k/ v_\alpha\, </math></td>
</tr>
<tr>
<td align="center" valign="middle"><math>v_k=\exp(v_\alpha)\, </math></td>
<td align="center" valign="middle"><math>\exp(v_\alpha)\,\,(=v_k)</math></td>
</tr>
<tr>
<td align="center" valign="middle"><math>v_k=\log(v_\alpha)\, </math></td>
<td align="center" valign="middle"><math>1/v_\alpha\, </math></td>
</tr>
<tr>
<td align="center" valign="middle"><math>v_k=\sqrt{v_\alpha}\, </math></td>
<td align="center" valign="middle"><math>1/(2\sqrt{v_\alpha})\, </math><br>
<p><math>(=0.5/v_k)\, </math><br>
</p>
</td>
</tr>
</table>
</td>
</tr>
</table></div>


　先の例では, <math>\partial v_1/\partial x=1, \partial v_2/\partial x=0</math> に注意して, <math>s_1=1\,</math>, <math>s_2=0\,</math>, <math>s_3=2*v_1*s_1\, </math>, <math>s_4=s_3+s_2\, </math>, <math>s_5=0.5/v_5*s_4\, </math>, <math>s_6=(1/v_5)*s_1+(-v_6/v_5)*s_5\, </math> という代入文の列を生成する. これを実行すると <math>s_6\,</math> には <math>(\partial f/\partial x)(3,4)\, </math> の値が計算される(<math>v_k\,</math> の計算の直後に <math>s_k\,</math> を計算してもよい). 高々2項までの基本演算だけ使用するという条件の下では, BUADの手間は <math>\mbox{O}(r)\, </math> である. <math>s_1=0\, </math>, <math>s_2=1\, </math> と一部変更し, もう一度計算すれば, <math>s_6\,</math> には, <math>(\partial f/\partial y)(3,4)</math> の値が計算される. <math>n\,</math>変数関数の勾配を計算するには, 同様の計算を <math>n\, </math>回繰り返す必要がある. 

　TDADはこれとは異なり, 先の計算過程を <math>\{-v_k+\varphi_k(v_1,\cdots,v_{k-1})=0\}_{k=1}^r</math> と書き直し, これらを <math>v_1, \cdots, v_r</math> に関する制約式とみなす. この制約の下で, <math>v_r\,</math> (<math>f\,</math>の値) の停留点を考える. ラグランジュ関数<math>\textstyle L(v_1,\cdots,v_r; \lambda_1,\cdots,\lambda_r)=v_r+\sum_{k=1}^r\lambda_k(-v_k+\varphi_k(v_1,\cdots,v_{k-1}))</math> の停留点(<math>\partial L/\partial \lambda_k=0</math> かつ<math>\partial L/\partial v_k=0</math> が成立する点)では, ラグランジュ乗数 <math>\lambda_k\, </math> は, <math>k\,</math>番目の制約式の摂動に対する関数値 <math>v_r\,</math> の感度を与えるが, <math>j=1,\cdots, n</math> については<math>\lambda_j\, </math> は <math>\partial f/\partial x_i</math> に等しい. 入力 <math>x_1, \cdots, x_n</math> を定めると<math>v_{1}, \cdots, v_r</math> は一意に定まるが, <math>\lambda_k\, </math> は連立一次方程式<math>\textstyle (\partial L/\partial v_r=)1+\lambda_r\cdot (-1)=0,(\partial L/\partial v_k=)\sum_{j=k+1}^r\lambda_j\cdot(\partial\varphi_j/\partial v_k) + \lambda_k\cdot(-1)=0 (k=r-1,\cdots,1)</math>を満たす. これを解くには, <math>\varphi_k</math> が実質的に単項・2項演算であることを考慮すると, <math>\lambda_r\gets 1, \lambda_{r-1}\gets 0,\cdots, \lambda_1\gets 0</math> と初期化しておき, <math>k=r-1,r-2,\cdots,1</math> の順に <math>\lambda_i\gets\lambda_i+\lambda_k\cdot(\partial \varphi_k/\partial v_i)(i=1,\cdots,k-1)</math> を計算する. 各 <math>k\,</math> について高々2個の <math>i\,</math> についてだけ計算すればよい. 

　先の例では, <math>v_1, \cdots, v_6</math> を計算し, <math>\lambda_6=1, \lambda_5=0, \cdots, \lambda_1=0</math> と初期化した後, <math>\lambda_1\gets\lambda_1+\lambda_6\cdot(1/v_5),</math><math>\lambda_5\gets\lambda_5+\lambda_6\cdot(-v_6/v_5),</math><math>\lambda_4\gets\lambda_4+\lambda_5\cdot(0.5/v_5),</math><math>\lambda_3\gets\lambda_3+\lambda_4\cdot1,\lambda_2\gets\lambda_2+\lambda_4\cdot1</math>, <math>\lambda_1\gets\lambda_1+\lambda_3\cdot(2v_1)</math> となる. 最終的に <math>\lambda_1, \lambda_2\, </math> に <math>(\partial f/\partial x)(3,4), (\partial f/\partial y)(3,4)</math> の値が計算される. 同じ条件の下で, TDADの手間は <math>\mbox{O}(r)\, </math>である. 1回の計算で勾配の値は全て計算できることに注意. 

　<math>n\,</math> 変数 <math>m\,</math> 値関数 <math>[f_1(x_1,\cdots,x_n),\cdots,f_m(x_1,\cdots,x_n)]^{\top}</math> について, 全成分の値を計算するのに延べ <math>r\,</math> 回の基本演算を実行したとする. ヤコビ行列 <math>J=(\partial f_i/\partial x_j)\, </math> の列の線形結合はBUADで, 行についてはTDADで <math>\mbox{O}(r)\, </math>の手間で計算できる. 全成分については BUADでは <math>\mbox{O}(nr)\, </math>, TDAD では <math>\mbox{O}(mr)\, </math> である. 

　実際には, 基本演算は表1に限らず, 代入文(やその列)を一つの基本演算とみなしてよい. また, プログラム中に条件分岐があっても, 与えられた入力値に関する関数の合成は上記の形で書けるから, ADを適用できる. ただし, 分岐の境目では, ADの結果は, 真の偏導関数値と異なることがある. たとえば, <math>\mbox{if(x=1.0)}\{\mbox{y=x*x}\}\mbox{else}\{\mbox{y=1.0}\}\, </math>の様なプログラムを自動微分すると, <math>x\,</math> の値が1.0 のときには不具合が起こりうるので注意が必要である. 



----
'''参考文献'''

[1] M. Berz, C. Bischof, G. Corliss and A. Griewank, ''Computational Differentiation: Techniques, Applications, and Tools'', SIAM, 1996.

[2]久保田光一, 伊理正夫, 『アルゴリズムの自動微分と応用』, コロナ社, 1998.