準ニュートン法のソースを表示

'''準ニュートン法'''（じゅんニュートンほう、{{lang-en-short|quasi-Newton method}}）とは、非線形連立方程式の解、あるいは連続[[最適化問題]]の[[関数 (数学)|関数]]の[[極値|極大・極小]]解を見つけるための[[アルゴリズム]]である。準ニュートン法は[[ニュートン法]]を元にしており、非線形連立方程式の解を求めることが基本になるが、最適化問題においては、関数の停留点を見つけるために、関数の[[勾配 (ベクトル解析)|勾配]]=0の連立方程式を解くという立場から考えることができる。以下は、主に最適化問題の立場からの説明である。

== 概要 ==
通常のニュートン法は最適解の近傍が[[二次関数]]で近似できると仮定し、一階および二階の[[導関数]]を解の探索に用いる。高次元空間上で定義される関数に対しては、最小化（最大化）したい関数の[[勾配 (ベクトル解析)|勾配]]ベクトルおよび[[ヘッセ行列]]を用いる。一方で、準ニュートン法では[[ヘッセ行列]]を陽に計算する必要がない。その代わりにヘッセ行列が最適化の繰り返し計算の過程で得られる勾配ベクトルにより近似される。準ニュートン法は多次元関数の零点 (関数の値が0となる場所) を探すアルゴリズムの一種である[[セカント法]]（割線法）の一般化であると見ることも出来る。多次元の問題においてはセカント方程式は1次元の場合と違い一意に定まらず、劣決定問題となるが、準ニュートン法は近似の制約が異なっており、具体的には現在の推定ヘッセ行列を低ランク行列成分を用いて更新する。

最初の準ニュートン法のアルゴリズムは物理学者のW.C. Davidonが[[アルゴンヌ国立研究所]]で研究中に提案したものである。彼が1959年に提案した最初の準ニュートン法は、1963年にFletcherとPowellが世に広め、後に[[DFP法]]（Davidon-Fletcher-Powell法）とよばれるようになったが、1970年後半、より効率的な[[BFGS法]]の登場により、今日ではあまり用いられていない。現在最も用いられている準ニュートン法のアルゴリズムは[[SR1法]]（対称ランクワン法、symmetric rank-one）、[[BHHH法]]、そして[[BFGS法]]（提案者であるBroyden, Fletcher, Goldfarb, Shannoの頭文字から）である。大規模問題に対応させる方法の一つとして記憶制限準ニュートン法が1980年に発表され<ref>{{cite journal|doi=10.1090/S0025-5718-1980-0572855-7|first=J. |last=Nocedal|title= Updating Quasi-Newton Matrices with Limited Storage |year=1980|journal=Mathematics of Computation |volume=35|pages=773–782|issue=151}}</ref>、BFGS法を記憶制限準ニュートン法にした物として[[L-BFGS法]]があり<ref>{{cite journal|doi= 10.1007/BF01589116|first=D. C.|last1= Liu |first2= J.|last2= Nocedal|url=http://www.ece.northwestern.edu/~nocedal/PSfiles/limited-memory.ps.gz|title= On the Limited Memory Method for Large Scale Optimization|year=1989|journal= Mathematical Programming B|volume=45|issue= 3|pages= 503–528}}</ref>、良く用いられるアルゴリズムの1つである。

SR1法は行列の更新が行列の[[行列の定値性|正定値性]]を保存しないため、不定値の行列のヘッセ行列に対しても用いることが出来る。また[[ブロイデン法]]は行列が対称行列でなくとも良く、通常の連立方程式の解を求めるのにも使うことが出来る。

通常のニュートン法と比べたときの準ニュートン法の最大の利点はヘッセ行列の逆行列を計算する必要がない更新法があるという点である。ニュートン法や、それを部分的に用いる[[内点法]]のアルゴリズムはこのヘッセ行列の逆行列を計算する部分が計算のボトルネックとなることが多い。これに対して準ニュートン法はヘッセ行列の逆行列自体を直接近似できる。

== 手法の説明 ==
ニュートン法と同様、関数 <math>f(\boldsymbol{x})</math> の解を見つけるために二次の近似を用いる。<math>f(\boldsymbol{x})</math> の二階の[[テイラー展開]]は
:<math>f(\boldsymbol{x}_k + \Delta \boldsymbol{x}) \approx f(\boldsymbol{x}_k) + \nabla f(\boldsymbol{x}_k)^\intercal \Delta \boldsymbol{x} + \frac{1}{2} \Delta \boldsymbol{x}^\intercal \boldsymbol{B} \Delta \boldsymbol{x} </math>
となる。この式で <math>\nabla f</math> は勾配を表し、<math>\boldsymbol{B}</math> はヘッセ行列の近似を表す。勾配 <math>\nabla f</math> はさらに次のように近似される。
:<math>\nabla f(\boldsymbol{x}_k + \Delta \boldsymbol{x}) \approx \nabla f(\boldsymbol{x}_k) + \boldsymbol{B} \Delta \boldsymbol{x}</math>
この勾配が0になると仮定するとニュートン・ステップが次の式で計算される。
:<math> \Delta \boldsymbol{x} = - \boldsymbol{B}^{-1} \nabla f(\boldsymbol{x}_k)</math>
そこでヘッセ行列の近似<math>\boldsymbol{B}</math>は次の式を満たすように行われる。
:<math>\nabla f(\boldsymbol{x}_k + \Delta \boldsymbol{x}) = \nabla f(\boldsymbol{x}_k) + \boldsymbol{B} \Delta \boldsymbol{x}</math>
この方程式はセカント方程式と呼ばれる。<math>f</math> が多次元空間上で定義される関数のとき、この式から <math>\boldsymbol{B}</math> を求める問題は劣決定問題になる(式の数よりも未知数の数が多い問題のこと)。このとき <math>\boldsymbol{B}</math> を求めて、ニュートン・ステップにより解を更新するという処理はセカント方程式を解くことに帰着される。多くの準ニュートン法はセカント方程式の解の選び方が異なっている。ほとんどの方法では <math>\boldsymbol{B} = \boldsymbol{B}^\intercal</math> という対称性を仮定して解を求める。加えて、以下のリストに示す方法では新たな近似 <math>\boldsymbol{B}_{k+1}</math> を得るために、その前の近似 <math>\boldsymbol{B}_k</math> と、ある[[ノルム]]の意味において近い解を選ぼうとする。すなわち何らかの正定値行列 <math>\boldsymbol{V}</math> に対して <math>\boldsymbol{B}_{k+1} = \arg \min_\boldsymbol{B} \| \boldsymbol{B} - \boldsymbol{B}_k \|_{\boldsymbol{V}}</math> と成るように <math>\boldsymbol{B}</math> を更新する。近似ヘッセ行列の初期値としては単位行列などが用いられる<ref>{{Cite book|author=William H. Pressほか|year=2007|title=Numerical Recepes|publisher=Cambridge Press|page=521-526|isbn=978-0-521-88068-8}}</ref>。最適化の解 <math>\boldsymbol{x}_k</math> は近似によって得られた <math>\boldsymbol{B}_k</math> を用いたニュートン・ステップにより更新される。

アルゴリズムをまとめると以下のようになる。
*<math>\Delta \boldsymbol{x}_k = -\alpha \boldsymbol{B}_k^{-1} \nabla f(\boldsymbol{x}_k) </math>
*<math>\boldsymbol{x}_{k+1} = \boldsymbol{x}_k + \Delta \boldsymbol{x}_k</math>
* 新しい点での勾配 <math>\nabla f(\boldsymbol{x}_{k+1})</math> を計算し <math>\boldsymbol{y}_k = \nabla f(\boldsymbol{x}_{k+1}) - \nabla f(\boldsymbol{x}_k)</math> とする
*<math>\boldsymbol{y}_k</math> を用いてヘッセ行列の逆行列 <math>\boldsymbol{B}_{k+1}^{-1}</math> を直接近似する。近似の式は各手法ごとに以下のリストの通り。

{| class="wikitable"

|-
! Method
! <math>\boldsymbol{B}_{k+1} = </math>
!<math>\boldsymbol{H}_{k+1} = \boldsymbol{B}_{k+1}^{-1} = </math>
|-

| [[DFP法]]

|<math>\left( \boldsymbol{I} - \frac{\boldsymbol{y}_k \, \Delta \boldsymbol{x}_k^\intercal} {\boldsymbol{y}_k^\intercal \, \Delta \boldsymbol{x}_k} \right) \boldsymbol{B}_k \left( \boldsymbol{I} - \frac{\Delta \boldsymbol{x}_k \boldsymbol{y}_k^\intercal}{\boldsymbol{y}_k^\intercal \, \Delta \boldsymbol{x}_k} \right) + \frac{\boldsymbol{y}_k \boldsymbol{y}_k^\intercal} {\boldsymbol{y}_k^\intercal \, \Delta \boldsymbol{x}_k}</math>

|<math> \boldsymbol{H}_k + \frac{\Delta \boldsymbol{x}_k \Delta \boldsymbol{x}_k^\intercal}{\boldsymbol{y}_k^{T} \, \Delta \boldsymbol{x}_k} - \frac{\boldsymbol{H}_k \boldsymbol{y}_k \boldsymbol{y}_k^\intercal \boldsymbol{H}_k^\intercal} {\boldsymbol{y}_k^\intercal \boldsymbol{H}_k \boldsymbol{y}_k}</math>

|-

| [[BFGS法]]

| <math> \boldsymbol{B}_k + \frac{\boldsymbol{y}_k \boldsymbol{y}_k^\intercal}{\boldsymbol{y}_k^\intercal \Delta \boldsymbol{x}_k} - \frac{\boldsymbol{B}_k \Delta \boldsymbol{x}_k (\boldsymbol{B}_k \Delta \boldsymbol{x}_k)^\intercal} {\Delta \boldsymbol{x}_k^{T} \boldsymbol{B}_k \, \Delta \boldsymbol{x}_k}</math>

| <math> \left( \boldsymbol{I} - \frac{\boldsymbol{y}_k \Delta \boldsymbol{x}_k^\intercal}{\boldsymbol{y}_k^\intercal \Delta \boldsymbol{x}_k} \right)^\intercal \boldsymbol{H}_k \left( \boldsymbol{I} - \frac{\boldsymbol{y}_k \Delta \boldsymbol{x}_k^\intercal}{\boldsymbol{y}_k^\intercal \Delta \boldsymbol{x}_k} \right) + \frac{\Delta \boldsymbol{x}_k \Delta \boldsymbol{x}_k^\intercal}{\boldsymbol{y}_k^\intercal \, \Delta \boldsymbol{x}_k}</math>

|-
| [[ブロイデン法]]

| <math>\boldsymbol{B}_k + \frac{\boldsymbol{y}_k - \boldsymbol{B}_k \Delta \boldsymbol{x}_k}{\Delta \boldsymbol{x}_k^\intercal \, \Delta \boldsymbol{x}_k} \, \Delta \boldsymbol{x}_k^\intercal</math>

|<math>\boldsymbol{H}_{k} + \frac{(\Delta \boldsymbol{x}_k - \boldsymbol{H}_k \boldsymbol{y}_k) \Delta \boldsymbol{x}_k^\intercal \boldsymbol{H}_k}{\Delta \boldsymbol{x}_k^\intercal \boldsymbol{H}_k \boldsymbol{y}_k}</math>

|-

| ブロイデン公式族

| <math>(1 - \varphi_k) \boldsymbol{B}_{k+1}^{\text{BFGS}}+ \varphi_k \boldsymbol{B}_{k+1}^{\text{DFP}}, \qquad \varphi \in [0,1]</math>

|
|-

| [[SR1法]]

| <math>\boldsymbol{B}_{k} + \frac{(\boldsymbol{y}_k - \boldsymbol{B}_k \, \Delta \boldsymbol{x}_k) (\boldsymbol{y}_k - \boldsymbol{B}_k \, \Delta \boldsymbol{x}_k)^\intercal}{(\boldsymbol{y}_k - \boldsymbol{B}_k \, \Delta \boldsymbol{x}_k)^\intercal \, \Delta \boldsymbol{x}_k}</math>

| <math>\boldsymbol{H}_{k} + \frac{(\Delta \boldsymbol{x}_k - \boldsymbol{H}_k \boldsymbol{y}_k) (\Delta \boldsymbol{x}_k - \boldsymbol{H}_k \boldsymbol{y}_k)^\intercal}{(\Delta \boldsymbol{x}_k - \boldsymbol{H}_k \boldsymbol{y}_k)^\intercal \boldsymbol{y}_k}</math>

|}

== 実装 ==
準ニュートン法は現在、汎用的に用いられている最適化アルゴリズムの1つであり、多くの[[プログラミング言語]]による実装が存在する。
* [[NAG数値計算ライブラリ]]には準ニュートン法を用いた関数の最大化・最小化アルゴリズムが存在する。
* [[Python]]のライブラリである[[SciPy]]では scipy.optimize.minimize のデフォルトは BFGS である<ref>[http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html scipy.optimize.minimize — SciPy Reference Guide]</ref>。
* [[GNU Octave]]は関数'funcmin'においてBFGS法を用いている。
* [[MATLAB]]のOptimization Toolboxにある関数fminuncでもBFGS法が実装されている。
* [[Mathematica]]や[[R]]にも一般的な関数の最適化法としてBFGS法が実装されている。

== 脚注 ==
{{脚注ヘルプ}}
{{Reflist}}

== 関連項目 ==
* [[ニュートン法]]

{{最適化アルゴリズム}}
{{デフォルトソート:しゆんにゆうとんほう}}
[[Category:最適化アルゴリズム]]
[[Category:数学に関する記事]]
[[Category:数学のエポニム]]