《高速微分法》

【こうそくびぶんほう (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 _{k}$ は$v_{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}$ となる.

表１：基本演算と要素的偏導関数
$\varphi _{k}$ $\partial \varphi _{k}/v_{\alpha }$ $\partial \varphi _{k}/v_{\beta }$ $v_{k}=v_{\alpha }\pm v_{\beta }\,$ $1\,$ $\pm 1$ $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}\cdot 1,\lambda _{2}\gets \lambda _{2}+\lambda _{4}\cdot 1$, $\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.