準ニュートン法

準ニュートン法（じゅんニュートンほう、英: 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年に発表され[1]、BFGS法を記憶制限準ニュートン法にした物としてL-BFGS法があり[2]、良く用いられるアルゴリズムの1つである。

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

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

手法の説明

ニュートン法と同様、関数 f ( x ) {\displaystyle f({\boldsymbol {x}})} の解を見つけるために二次の近似を用いる。 f ( x ) {\displaystyle f({\boldsymbol {x}})} の二階のテイラー展開は

f ( x k + Δ x ) ≈ f ( x k ) + ∇ f ( x k ) ⊺ Δ x + 1 2 Δ x ⊺ B Δ x {\displaystyle 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}}}

となる。この式で ∇ f {\displaystyle \nabla f} は勾配を表し、 B {\displaystyle {\boldsymbol {B}}} はヘッセ行列の近似を表す。勾配 ∇ f {\displaystyle \nabla f} はさらに次のように近似される。

∇ f ( x k + Δ x ) ≈ ∇ f ( x k ) + B Δ x {\displaystyle \nabla f({\boldsymbol {x}}_{k}+\Delta {\boldsymbol {x}})\approx \nabla f({\boldsymbol {x}}_{k})+{\boldsymbol {B}}\Delta {\boldsymbol {x}}}

この勾配が0になると仮定するとニュートン・ステップが次の式で計算される。

Δ x = − B − 1 ∇ f ( x k ) {\displaystyle \Delta {\boldsymbol {x}}=-{\boldsymbol {B}}^{-1}\nabla f({\boldsymbol {x}}_{k})}

そこでヘッセ行列の近似 B {\displaystyle {\boldsymbol {B}}} は次の式を満たすように行われる。

∇ f ( x k + Δ x ) = ∇ f ( x k ) + B Δ x {\displaystyle \nabla f({\boldsymbol {x}}_{k}+\Delta {\boldsymbol {x}})=\nabla f({\boldsymbol {x}}_{k})+{\boldsymbol {B}}\Delta {\boldsymbol {x}}}

この方程式はセカント方程式と呼ばれる。 f {\displaystyle f} が多次元空間上で定義される関数のとき、この式から B {\displaystyle {\boldsymbol {B}}} を求める問題は劣決定問題になる(式の数よりも未知数の数が多い問題のこと)。このとき B {\displaystyle {\boldsymbol {B}}} を求めて、ニュートン・ステップにより解を更新するという処理はセカント方程式を解くことに帰着される。多くの準ニュートン法はセカント方程式の解の選び方が異なっている。ほとんどの方法では B = B ⊺ {\displaystyle {\boldsymbol {B}}={\boldsymbol {B}}^{\intercal }} という対称性を仮定して解を求める。加えて、以下のリストに示す方法では新たな近似 B k + 1 {\displaystyle {\boldsymbol {B}}_{k+1}} を得るために、その前の近似 B k {\displaystyle {\boldsymbol {B}}_{k}} と、あるノルムの意味において近い解を選ぼうとする。すなわち何らかの正定値行列 V {\displaystyle {\boldsymbol {V}}} に対して B k + 1 = arg ⁡ min B ‖ B − B k ‖ V {\displaystyle {\boldsymbol {B}}_{k+1}=\arg \min _{\boldsymbol {B}}\|{\boldsymbol {B}}-{\boldsymbol {B}}_{k}\|_{\boldsymbol {V}}} と成るように B {\displaystyle {\boldsymbol {B}}} を更新する。近似ヘッセ行列の初期値としては単位行列などが用いられる[3]。最適化の解 x k {\displaystyle {\boldsymbol {x}}_{k}} は近似によって得られた B k {\displaystyle {\boldsymbol {B}}_{k}} を用いたニュートン・ステップにより更新される。

アルゴリズムをまとめると以下のようになる。

  • Δ x k = − α B k − 1 ∇ f ( x k ) {\displaystyle \Delta {\boldsymbol {x}}_{k}=-\alpha {\boldsymbol {B}}_{k}^{-1}\nabla f({\boldsymbol {x}}_{k})}
  • x k + 1 = x k + Δ x k {\displaystyle {\boldsymbol {x}}_{k+1}={\boldsymbol {x}}_{k}+\Delta {\boldsymbol {x}}_{k}}
  • 新しい点での勾配 ∇ f ( x k + 1 ) {\displaystyle \nabla f({\boldsymbol {x}}_{k+1})} を計算し y k = ∇ f ( x k + 1 ) − ∇ f ( x k ) {\displaystyle {\boldsymbol {y}}_{k}=\nabla f({\boldsymbol {x}}_{k+1})-\nabla f({\boldsymbol {x}}_{k})} とする
  • y k {\displaystyle {\boldsymbol {y}}_{k}} を用いてヘッセ行列の逆行列 B k + 1 − 1 {\displaystyle {\boldsymbol {B}}_{k+1}^{-1}} を直接近似する。近似の式は各手法ごとに以下のリストの通り。
Method B k + 1 = {\displaystyle {\boldsymbol {B}}_{k+1}=} H k + 1 = B k + 1 − 1 = {\displaystyle {\boldsymbol {H}}_{k+1}={\boldsymbol {B}}_{k+1}^{-1}=}
DFP法 ( I − y k Δ x k ⊺ y k ⊺ Δ x k ) B k ( I − Δ x k y k ⊺ y k ⊺ Δ x k ) + y k y k ⊺ y k ⊺ Δ x k {\displaystyle \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}}}} H k + Δ x k Δ x k ⊺ y k T Δ x k − H k y k y k ⊺ H k ⊺ y k ⊺ H k y k {\displaystyle {\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}}}}
BFGS法 B k + y k y k ⊺ y k ⊺ Δ x k − B k Δ x k ( B k Δ x k ) ⊺ Δ x k T B k Δ x k {\displaystyle {\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}}}} ( I − y k Δ x k ⊺ y k ⊺ Δ x k ) ⊺ H k ( I − y k Δ x k ⊺ y k ⊺ Δ x k ) + Δ x k Δ x k ⊺ y k ⊺ Δ x k {\displaystyle \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}}}}
ブロイデン法 B k + y k − B k Δ x k Δ x k ⊺ Δ x k Δ x k ⊺ {\displaystyle {\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 }} H k + ( Δ x k − H k y k ) Δ x k ⊺ H k Δ x k ⊺ H k y k {\displaystyle {\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}}}}
Broyden family ( 1 − φ k ) B k + 1 BFGS + φ k B k + 1 DFP , φ ∈ [ 0 , 1 ] {\displaystyle (1-\varphi _{k}){\boldsymbol {B}}_{k+1}^{\text{BFGS}}+\varphi _{k}{\boldsymbol {B}}_{k+1}^{\text{DFP}},\qquad \varphi \in [0,1]}
SR1法（英語版） B k + ( y k − B k Δ x k ) ( y k − B k Δ x k ) ⊺ ( y k − B k Δ x k ) ⊺ Δ x k {\displaystyle {\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}}}} H k + ( Δ x k − H k y k ) ( Δ x k − H k y k ) ⊺ ( Δ x k − H k y k ) ⊺ y k {\displaystyle {\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}}}}

実装

準ニュートン法は現在、汎用的に用いられている最適化アルゴリズムの1つであり、多くのプログラミング言語による実装が存在する。

  • NAG数値計算ライブラリには準ニュートン法を用いた関数の最大化・最小化アルゴリズムが存在する。
  • PythonのライブラリであるSciPyでは scipy.optimize.minimize のデフォルトは BFGS である[4]。
  • GNU Octaveは関数'funcmin'においてBFGS法を用いている。
  • MATLABのOptimization Toolboxにある関数fminuncでもBFGS法が実装されている。
  • MathematicaやRにも一般的な関数の最適化法としてBFGS法が実装されている。

脚注

[脚注の使い方]
  1. ^ Nocedal, J. (1980). “Updating Quasi-Newton Matrices with Limited Storage”. Mathematics of Computation 35 (151): 773–782. doi:10.1090/S0025-5718-1980-0572855-7. 
  2. ^ Liu, D. C.; Nocedal, J. (1989). “On the Limited Memory Method for Large Scale Optimization”. Mathematical Programming B 45 (3): 503–528. doi:10.1007/BF01589116. http://www.ece.northwestern.edu/~nocedal/PSfiles/limited-memory.ps.gz. 
  3. ^ William H. Pressほか (2007). Numerical Recepes. Cambridge Press. p. 521-526. ISBN 978-0-521-88068-8 
  4. ^ scipy.optimize.minimize — SciPy Reference Guide

関連項目

  • 表
  • 話
  • 編
  • 歴
非線形(無制約)
… 関数　
勾配法
収束性
準ニュートン法
その他の求解法
… ヘッセ行列
  • 最適化におけるニュートン法（英語版）
The graph of a strictly concave quadratic function is shown in blue, with its unique maximum shown as a red dot. Below the graph appears the contours of the function: The level sets are nested ellipses.
Optimization computes maxima and minima.
非線形(制約付き)
一般
微分可能
凸最適化
凸縮小化
  • 切除平面法（英語版、デンマーク語版、ドイツ語版、スペイン語版）
  • 簡約勾配法
  • 劣勾配法（英語版）
線型 および
二次
内点法
ベイズ-交換
  • 単体法
  • 改訂シンプレックス法（英語版）
  • 十字法（英語版）
  • レムケの主ピボット操作法（英語版）
組合せ最適化
系列範例
(Paradigms)
グラフ理論
最小全域木
最大フロー問題
メタヒューリスティクス
  • カテゴリ（最適化 • アルゴリズム） • ソフトウェア（英語版）