「《高速微分法》」の版間の差分

提供: ORWiki
ナビゲーションに移動 検索に移動
(新しいページ: ''''【こうそくびぶんほう (fast differentiation)】'''  非線形関数の勾配, ヤコビ行列, ヘッセ行列等の値を数値的に計算する方...')
 
5行目: 5行目:
 
 以下,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 と 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>{\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>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> となる.  
+
 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">
 
<div align="center">
 
<table width="600" border="0" cellspacing="0" cellpadding="0">
 
<table width="600" border="0" cellspacing="0" cellpadding="0">
<tr>
+
<tr>
<td colspan="2" align="center">表1:基本演算と要素的偏導関数<br>
+
<td colspan="2" align="center">表1:基本演算と要素的偏導関数<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>
<td align="center"><table width="200" border="1" cellspacing="2" cellpadding="2">
+
</tr>
<tr>
+
<tr>
<td align="center" valign="middle"><math>\varphi_k\, </math></td>
+
<td align="center"><table width="200" border="1" cellspacing="2" cellpadding="2">
<td align="center" valign="middle"><math>\partial \varphi_k/ v_\alpha\, </math></td>
+
<tr>
</tr>
+
<td align="center" valign="middle"><math>\varphi_k</math></td>
<tr>
+
<td align="center" valign="middle"><math>\partial \varphi_k/ v_\alpha</math></td>
<td align="center" valign="middle"><math>v_k=\exp(v_\alpha)\, </math></td>
+
<td align="center" valign="middle"><math>\partial \varphi_k/ v_\beta</math></td>
<td align="center" valign="middle"><math>\exp(v_\alpha)\,\,(=v_k)</math></td>
+
</tr>
</tr>
+
<tr>
<tr>
+
<td align="center" valign="middle"><math>v_k=v_\alpha \pm v_\beta\, </math></td>
<td align="center" valign="middle"><math>v_k=\log(v_\alpha)\, </math></td>
+
<td align="center" valign="middle"><math>1\, </math></td>
<td align="center" valign="middle"><math>1/v_\alpha\, </math></td>
+
<td align="center" valign="middle"><math>\pm1</math></td>
</tr>
+
</tr>
<tr>
+
<tr>
<td align="center" valign="middle"><math>v_k=\sqrt{v_\alpha}\, </math></td>
+
<td align="center" valign="middle"><math>v_k=v_\alpha * v_\beta\, </math></td>
<td align="center" valign="middle"><math>1/(2\sqrt{v_\alpha})\, </math><br>
+
<td align="center" valign="middle"><math>v_\beta\, </math></td>
<p><math>(=0.5/v_k)\, </math><br>
+
<td align="center" valign="middle"><math>v_\alpha\, </math></td>
</p>
+
</tr>
</td>
+
<tr>
</tr>
+
<td align="center" valign="middle"><math>v_k=v_\alpha / v_\beta\, </math></td>
</table>
+
<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>
</tr>
+
<td align="center"><table width="200" border="1" cellspacing="2" cellpadding="2">
</table></div>
+
<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>回繰り返す必要がある.  
 
 先の例では, <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>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>(\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> についてだけ計算すればよい.  
+
 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>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回の計算で勾配の値は全て計算できることに注意.  

2007年7月15日 (日) 00:23時点における版

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

 以下,BUAD と TDAD による計算法を説明する.例として,2変数関数 について, の値を計算する代入文の列(プログラム), を考えよう. ただし, 各代入文の右辺には, 演算(基本演算とよぶ)が高々1回だけ現れるとする. , , に対応し, の値が計算される. 一般には, 変数関数 について, 回目の代入文には, 回目までに計算される変数が現れうるから, 延べ 回の演算を行なう代入文の列はと表される. これを計算過程といい, を中間変数という. のとき という入力定数の代入演算に相当する.

 BUADは, 補助変数 を導入し, 任意に を固定して, 合成関数の に関する偏微分則 に基づき, を計算する式を導出する. 基本演算 を四則演算や初等関数などの2項・単項の演算に限れば, 表1により, (これを要素的偏導関数という)を導出できる. , と初期設定すれば, について を計算済みとみなすことができ, の値を計算できる. 最終的に となる.

表1:基本演算と要素的偏導関数





 先の例では, に注意して, , , , , , という代入文の列を生成する. これを実行すると には の値が計算される( の計算の直後に を計算してもよい). 高々2項までの基本演算だけ使用するという条件の下では, BUADの手間は である. , と一部変更し, もう一度計算すれば, には, の値が計算される. 変数関数の勾配を計算するには, 同様の計算を 回繰り返す必要がある.

 TDADはこれとは異なり, 先の計算過程を と書き直し, これらを に関する制約式とみなす. この制約の下で, (の値) の停留点を考える. ラグランジュ関数 の停留点( かつ が成立する点)では, ラグランジュ乗数 は, 番目の制約式の摂動に対する関数値 の感度を与えるが, については に等しい. 入力 を定めると は一意に定まるが, は連立一次方程式を満たす. これを解くには, が実質的に単項・2項演算であることを考慮すると, と初期化しておき, の順に を計算する. 各 について高々2個の についてだけ計算すればよい.

 先の例では, を計算し, と初期化した後, , となる. 最終的に の値が計算される. 同じ条件の下で, TDADの手間は である. 1回の計算で勾配の値は全て計算できることに注意.

  変数 値関数 について, 全成分の値を計算するのに延べ 回の基本演算を実行したとする. ヤコビ行列 の列の線形結合はBUADで, 行についてはTDADで の手間で計算できる. 全成分については BUADでは , TDAD では である.

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



参考文献

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

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