《大規模問題の分解法》のソースを表示

'''【だいきぼもんだいのぶんかいほう (decomposition method for large-scale problems)】'''

　現実世界で発生する複雑な問題を最適化問題としてモデル化すると，非常に多くの変数や制約条件をもつ大規模問題になることが多い．計算機を用いてその解を求める場合，目的関数や制約条件式がすべて線形であっても，変数や制約条件の数が増えるとともに計算時間は急激に増加する．また，それらの一部に非線形の項が含まれると，問題の解きにくさは飛躍的に増大する．そこで，大規模問題を直接解くのではなく，一部の変数や制約条件だけから成る小規模な，または，解きやすい部分問題を逐次解くことにより，もとの問題の解を得ようとするアルゴリズムが提案されている．それらを一般に[[大規模問題の分解法]] (decomposition method for large-scale problems) と呼ぶ [3]．

　現実の大規模問題は，しばしば特徴的なブロック構造をもつ．たとえば，<math>n\, </math>個の比較的独立な部分から成るシステムが共通の変数 <math>x_{0}\, </math>を含む場合は，


:<math>\begin{array}{llrlll}
    \mbox{min.} & f_{0}(x_{0}) + & \sum_{j=1}^{n} f_{j}(x_{j})\\
    \mbox{s. t.} & g_{0}(x_{0})  & & \leq 0, & & \quad \mbox{(1)}\\
                 & g_{j}(x_{0}) + & h_{j}(x_{j}) & \leq 0  & 
                                         (j = 1, \ldots, n), 
\end{array}\, </math>


という最適化問題が得られる．この問題は，変数 <math>x_{0}\, </math> を一時的に固定すると，


:<math>\begin{array}{llll}
\mbox{min.} \quad f_{j}(x_{j}) \quad 
  \mbox{s. t.} \quad g_{j}(x_{0}) + h_{j}(x_{j}) \leq 0,  & & \quad \mbox{(2)}
\end{array}\, </math>


という <math>n\, </math> 個の部分問題に分解される．[[ベンダース分解法]] (Benders decomposition method) はこのような性質を利用しており，関数 <math>f_{j}\, </math>, <math>h_{j}\, </math> が線形ならば，有限回の反復でもとの問題の解に到達できることが知られている．さらに，分解された部分問題はたがいに独立であり，それらは比較的大きい（粒度が粗い）問題となるため， MIMD (multiple instruction stream multiple data stream) 型の並列計算機で効率よく実行できる．

　また，システム全体にまたがる付加的な制約条件が存在する場合は，


:<math>\begin{array}{lrllll}
    \mbox{min.} & f_{0}(x_{0}) + \sum_{j=1}^{n} f_{j}(x_{j}) &         & \\
    \mbox{s. t.} & g_{0}(x_{0}) + \sum_{j=1}^{n} g_{j}(x_{j}) & \leq 0, & & \quad \mbox{(3)}\\
                 &                             h_{j}(x_{j})   & \leq 0  & (j = 1, \ldots, n),  
\end{array}\, </math>


という最適化問題が得られるが，<math>g_{j}\, </math> を含む制約条件をラグランジュ緩和により目的関数に組み込み, さらに変数 <math>x_0\, </math> を一時的に固定すると, 問題 (2) に類似した<math>n\, </math> 個の独立な部分問題に分解される．とくにすべての関数が線形ならば，問題 (3) は部分問題の解を用いて効率的に解けることが知られており，[[ダイツィク・ウルフ分解法]] (Dantzig-Wolfe decomposition method) と呼ばれている．

　一方，大規模で複雑な問題から取り扱いやすい構造をもつ部分のみを抽出して部分問題を構成し，これを逐次解くことによりもとの問題の解を得ようとする方法は，分割法 (splitting method) と呼ばれている．分割法は，問題 (1) や (3) のようなブロック構造をもたない問題にも適用可能である．また，並列処理可能な部分問題を構成すれば，それらはしばしば小規模な（粒度の細かい）問題となり，SIMD (single instruction stream multiple data stream) 型の大規模並列計算機を用いて効率的に実行できる．

　分割法は，線形方程式に対する反復法として，線形代数の分野において古くから研究されている．行列 <math>M\, </math> とベクトル <math>q\, </math> により定義される線形方程式


:<math>M x + q = 0\, </math>


に対して，条件 <math>M = B + C\, </math> を満たす行列 <math>B\, </math>, <math>C\, </math> を選び，方程式


:<math>\begin{array}{ll}
  B x + C x^{(k)} + q = 0 & \mbox{(4)}
\end{array}\, </math>


の解を <math>x^{(k+1)}\, </math> とおくことにより点列 <math>\{ x^{(k)} \}\, </math> を生成する方法は，[[行列分割法]] (matrix splitting method) と呼ばれる．行列 <math>B\, </math> を (ブロック) 対角行列に選べば方程式 (4) は並列的に解けるので，大規模問題に対する効率的な解法を得る．行列分割法は，線形相補性問題や線形変分不等式問題にも拡張できる．

　一般の凸計画問題に対しても，分割法に基づくアルゴリズムを構成できる．凸計画問題は，最適性条件を考慮すると，ある写像 <math>F: \mathbf{R}^{n} \rightarrow \mathbf{R}^{n}\, </math> を用いて


:<math>\begin{array}{ll}
  \mbox{find} \quad x \in \mathbf{R}^{n} \quad 
  \mbox{s. t.} \quad 0 \in F(x), & \quad \mbox{(5)}
\end{array}\, </math>


と記述される．条件 <math>F = G + H\, </math> を満たす写像 <math>G\, </math>, <math>H\, </math> を選び，部分問題


:<math>\begin{array}{ll}
  \mbox{find} \quad x \in \mathbf{R}^{n} \quad 
  \mbox{s. t.} \quad 0 \in G(x) + H(x^{(k)}), & \quad \mbox{(6)}
\end{array}\, </math>


の解を <math>x^{(k+1)}\, </math> とする反復法は，[[作用素分割法]] (operator splitting method) と呼ばれる．写像 <math>G\, </math> が分離可能ならば部分問題 (6) は並列的に解けるので，大規模な凸計画問題に対する効率的な解法が得られる．

　問題 (5) に対する有力な反復法に，[[近接点法]] (proximal point method) がある．近接点法では，単調非減少な正定数の列 <math>\{ \lambda^{(k)} \}\, </math> を定め，問題


:<math>\mbox{find} \quad x \in \mathbf{R}^{n} \quad 
\mbox{s. t.} \quad ( x^{(k)} - x ) \, / \, \lambda^{(k)} \in F(x),\, </math> 


の解を <math>x^{(k+1)}\, </math> とおく．作用素分割法や近接点法は，一般的な最適化問題のクラスである変分不等式問題にも拡張できる．また，これらの方法を組合せることによりさまざまな[[並列アルゴリズム(数理計画問題の)]](parallel algorithm (nonlinear programming)) を構成できることが知られている．このような考え方に基づく並列アルゴリズムについては，参考文献 [1, 2] に詳しい解説がある．



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

[1]  D. P. Bertsekas and J. N. Tsitsiklis, ''Parallel and Distributed Computation: Numerical Methods'', Prentice-Hall, 1989.  

[2] Y. Censor and S. A. Zenios, ''Parallel Optimization: Theory, Algorithms, and Applications'', Oxford University Press, 1997.  

[3] J. F. Shapiro, ''Mathematical Programming: Structures and Algorithms'', John Wiley & Sons, 1979.


●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●

●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●

'''【こうそくびぶんほう (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>{\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> となる. 

<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>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> についてだけ計算すればよい. 

　先の例では, <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.