《高速微分法》

提供:ORWiki
移動: 案内, 検索

【こうそくびぶんほう (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 は「関数の値を計算するプログラム」から「偏導関数の値を計算するプログラム」を生成する手順を与え, 生成物を(コンパイルし)実行すれば, 差分商近似のような打ち切り誤差無しで, 正確な偏導関数の値を計算できる. 大規模システムの数学モデル等の大規模プログラム(数千行以上)により表現される関数の偏導関数を計算できるのが特長. n\,変数関数の勾配のn\,個の値を関数計算の手間の定数倍で計算できる点が「高速」微分の由来である.

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

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

表1:基本演算と要素的偏導関数
\varphi_k \partial \varphi_k/ v_\alpha \partial \varphi_k/ v_\beta
v_k=v_\alpha \pm v_\beta\, 1\, \pm1
v_k=v_\alpha * v_\beta\, v_\beta\, v_\alpha\,
v_k=v_\alpha / v_\beta\, 1/v_\beta\, -v_\alpha/({v_\beta}^2)\,

(=-v_k/v_\beta)\,

\varphi_k\, \partial \varphi_k/ v_\alpha\,
v_k=\exp(v_\alpha)\, \exp(v_\alpha)\,\,(=v_k)
v_k=\log(v_\alpha)\, 1/v_\alpha\,
v_k=\sqrt{v_\alpha}\, 1/(2\sqrt{v_\alpha})\,

(=0.5/v_k)\,


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

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

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

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

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



参考文献

[1] M.Bucker, G.Corliss, P.Hovland, U.Naumann, and B.Norris (eds.), Automatic Differentiation: Applications, Theory, and Implementations, Lecture Notes in Computational Science and Engineering, Springer, Vol.50, 2006.

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

個人用ツール