非線形最適化(1) - 最急降下法, Newton法, GN法
初めに
工学の分野で遭遇するモデルの多くは非線形なモデルが殆どで、それらのモデルによって構築される制約条件や最小化したい目的関数も非線形となる。そして、構築したモデルを用いて、条件を満たす若しくは目的関数を最小化するパラメータの推定が行われる。非線形最適化の手法は、それらの場面で活躍する。
非線形のモデルで表される問題を解くためには、それを解くためのソルバーを実装する。しかし、実装を一から構築する必要はなく、オープンソースの非線形ソルバーを使えばよい。問題を解くための道具は既に用意はされている。ただ、全ての問題に対して万能な解法はなく、また十分に性能を高めるには問題の構造に応じて適切な方法・パラメータを選択する必要がある。ソルバーの関数群はコンパクトにまとまっていて、使い易く整備はされているが、それらを正しく使いこなすには非線形最適化についての基本的な理解が必要となる。
本エントリーでは、まず非線形最適化の概要について簡単に述べる。続いて、基本的な手法である最急降下法、Newton法、Gauss-Newton (GN)法についてまとめる。これらの方法について詳細について解説した記事は、優秀な解説が数多く散見される。正確性や具体的な理論の解説はそれらに譲るとして、なるべく感覚的に分かるような説明を目指してまとたいと思う。また、より実用的に用いられているGauss-Newton (GN)法を改良したLevenberg Marquardt(LM)法については少々内容も複雑なので、別のエントリーにてまとめようと思います。
非線形最適化とは
非線形最適化とは、文字通り線形でない関数$f({\bf x})$で表された条件を満たす変数を求める問題である。具体的には、変数${\bf x} \in R ^n$について
$$ f({\bf x})=0 \tag{1}$$
を満たす${\bf x ^ *}$を求めたいとする。この${\bf x ^ *}$を求めるための方法が非線形最適化という分野によって扱われる。非線形最適化にも、さらに制約ありの場合となしの場合に大別される。本エントリーは、基本となる制約なしの場合についてとりあげる。
非線形最適化が必要になる背景として、関数$f({\bf x})$が線形関数である場合の例を最初に確認する。最も単純な具体例として、行列$\bf{A} \in R ^{n \times n}$とベクトル$\bf{b} \in R ^n$を用いて式(2)の様に条件が表される場合を考える。
$${\bf A} {\bf x} - {\bf b} = 0 \tag{2}$$
この時、行列$\bf{A}$が逆行列をもつ、つまり正則であれば条件(2)を満たす$\bf{x ^ *}$は逆行列を使って式(3)のように表すことができる。 $${\bf x ^ *}={\bf A} ^{-1}{\bf b} \tag{3}$$ したがって、条件式に登場する行列${\bf A}$とベクトル${\bf b}$をそのまま使って解析的に${\bf x ^ *}$を求めることが出来る。
世の中の様々な分野で扱う数理モデルは式(2)のような線形なモデルである事は珍しく、非線形である事の方が多い。その場合、式(3)のように解析的に式(1)を厳密に満たす$\bf{x ^ *}$を求めることは困難となる。そこで、解析的に厳密な解を求める代わりに左辺の$f({\bf x})$を最小にする${\bf x}$を求めるということを行う。
$${\bf x ^ *}=\underset{{\bf x} \in R}{argmin} \space f({\bf x}) \tag{4}$$
この最小となる${\bf x}$を求めるための理論を与えてくれるのが非線形最適化という分野である。続いて、非線形最適化について具体的に見ていく。
全体像
制約なしの非線形最適化についてその全体像を確認する。非線形最適化の基本的なアイデアは、微分を使って関数$f({\bf x})$の極小値を求めることで最小の$\bf{x ^ *}$ を求める。関数$f({\bf x})$が極小である時$f({\bf x})$の導関数$f'({\bf x})$は0となる。つまり、$f'({\bf x})=0$となる$\bf{x ^ *}$を求める。
しかし、高校数学でも習うとおり、$f'({\bf x})=0$となる${\bf x}$は図1に示すとおり3通りあり得る。なので、${\bf x ^ *}$を求めるためには、さらに$f'({\bf x})$がどう変化するのかも同時に調べる必要がある。同じ$f'({\bf x})=0$となる${\bf x}$であっても、傾き$f'({\bf x})$が減少し正から負へ転じる場合は極大、符号が変化しない場合は鞍点であり、増加し負から正へと転じる場合が極小である。二次導関数$f''({\bf x})$を用いればそれぞれ、$f''({\bf x}) <0$、$f''({\bf x})=0$、$f''({\bf x})>0$の場合である。
ここで更に、極小$\neq$最小である事に注意しなければならない。図2に描いたとおり、$f'({\bf x})=0 \land f''(x) > 0$で合っても大域的に最小であるとは限らない。
以上で述べたとおりで、非線形最適化では$f'({\bf x})=0$となる${\bf x}$を求めるが、その際に傾き$f'({\bf x})$の符号と大域的にそれが最小であるかをケアしながら行うことが重要となる。非線形最適化では、後述するとおりいくつかの求解法があるが共通して探索的に${\bf x}$を微小変化させながら$f'({\bf x})=0$となる${\bf x}$を求める。どの方法も共通して、
$${\bf x}_{k+1} = {\bf x}_k + \delta {\bf x} \tag{5}$$
で表される漸化式でステップ毎に$\delta {\bf x}$を求めながら${\bf x}_k$を繰り返し更新し、収束判定を行った段階での${\bf x}_k$を解とする。この$\delta {\bf x}$をどう求めるかが方法によって異なる。
$\delta {\bf x}$を求める最も基本的な方法として最急降下(Gradient Discent)法があり、その他Newton法、Gauss-Newton (GN)法、Levenberg Marquardt(LM)法といった方法がある。どれも先述したとおり”関数$f{{\bf (x)}}$の傾き”の情報を使って$\delta {\bf x}$を求め、式(5)の漸化式によって”繰り返し”${\bf x}$を変化させながら最小値を求める。図３にそれぞれの方法の特徴についてまとめた。
勾配だけをみて更新方向・量を決めるのが勾配降下法で、最もシンプルで演算量も少なく膨大なパラメータを持つモデルにも適用出来る。昨今注目を集めているニューラルネットワークもこの勾配降下法を用いている。
ニュートン法は、勾配だけでなくグラフの形状をさらにみて更新方向・量を決める方法で勾配降下法より収束性がよい。一方で2次導関数を求める必要があって演算量が多く、また局所最適解に陥りやすいという問題がある。最小2乗問題に限定してこれらの問題を解決したのが、ガウス・ニュートン法およびレーベンバーグ・マーカート法である。
次の項からは、これらの方法について最急降下法から順に解説をしていく。
最急降下法
最急降下法のアイデアはシンプルで、傾きが減少する方向に${\bf x}$を更新し続ければ最小値を与える${\bf x}$に辿り着くはず、と考える。具体的には、各ステップ${k}$において最も$f({\bf x})$が減少する方向$\delta {\bf x}$で${\bf x}$を更新し続ける。図４(a)は最急降下法のイメージで、各ステップの${\bf x}_k$において$f({\bf x})$が減少する方向つまり勾配$f({\bf x}_k)$の符号を反転させた方向に係数$\alpha$を乗じた分だけ${\bf x}_k$を更新している。
つまり、最急降下法では式(5)の漸化式が式(6)のようになる。 $$ {\bf x}_{k+1} = {\bf x}_k - \alpha \left.\frac{\partial f}{\partial {\bf x}}\right|_{{\bf x}_{k}} \tag{6} $$ 係数$\alpha$はハイパーパラメータで、毎回の更新の度合いを決める値であり、問題毎に適当な値を設定する必要がある。最急降下法はシンプルに勾配$\left.\frac{\partial f}{\partial {\bf x}}\right|_{{\bf x}_{k}}$さえ求まれば良いという扱いやすさはあるが、大雑把に${\bf x}_k$の周りで1次の線形関数に近似して減少方向を求めているため関数の具体的な形状をあまり考慮しておらず、図4(b)のように極小点周りで振動してしまい中々収束しないという問題がある。これに対し、次に述べるNewton法では2次の項まで考慮して$\delta {\bf x}$を求めることで収束速度をより速くすることが出来ます。
Newton法
Newton法では、各ステップの${\bf x}_k$周りに$f({\bf x})$を2次の項まで含めてより正確に近似し、それを最小にし続けるように${\bf x}_k$を更新し続け、全体の最小にたどり着くという方法をとる。最急降下法と同様に傾きの方向に応じて更新方向と大きさが変わるのと同時に、各${\bf x}_k$周りの微小領域においてより関数の形状を精度良く近似し、勾配$f'({\bf x}_k+\delta {\bf x})=0$となるように更新し続けるので図5(a)のように極小値周りで振動することなく高速に極小値を求めることが出来る。
簡単のためにまず1次元の場合から具体的な式を追ってみる。${\bf x}_k$周りに$f({\bf x})$をTaylor展開すると、
$$
f(x+dx) = f(x)+f'(x)dx+\frac{1}{2}f''(x)dx^2+O(dx^3) \tag{7}
$$
となる。$O(dx^3)$を無視した上で、この左辺を最小にする$dx$は、単純に2次関数の最小値を求めるために$dx$について平方完成し、
$$
f(x+dx) = \frac{1}{2}f''(x) \left( dx+\frac{f'(x)}{f''(x)} \right) +f(x)-\frac{f'(x)^2}{f''(x)} \tag{8}
$$
となるので、漸化式(5)は式(9)のようになる。 $$ {\bf x}_{k+1} = {\bf x}_{k} - \frac{f'(x)}{f''(x)} \tag{9} $$ イメージとしては図5(a)のように各${\bf x}_k$で2次関数を近似してフィッティングしてその最小値を求めることを繰り返すのがNewton法である。平方完成してその中心を求めているだけなので、最小であるためには当然下に凸であること、即ち$f''({\bf x})>0$であることが必要となる。
一般の多次元の場合も同様で、多変数関数のTaylor展開を用いると${\bf x} \in R ^n,{\bf \delta x} \in R ^n$について
$$
\begin{align}
f({\bf x} + \delta{\bf x}) &= \sum\_{i=0} ^n \frac{1}{n!} \left ( \delta {\bf x} \cdot \nabla \right ) ^i f({\bf x}) \\\
&=f({\bf x}) + \delta {\bf x} \cdot \nabla f+\frac{1}{2} \left ( \delta {\bf x} \cdot \nabla \right ) ^2f + O(\delta {\bf x} ^3) \\\
&=f({\bf x}) + \delta {\bf x} \cdot \nabla f+\frac{1}{2}\delta {\bf x} ^T (\nabla ^2 f)\delta {\bf x} + O(\delta {\bf x} ^3)
\end{align} \tag{10}
$$
となる。 ここで、${\bf H}({\bf x})$：ヘシアン行列(Hessian Matrix)を以下で定義する。
$$
{\bf H}({\bf x}) = \begin{bmatrix}
\frac{\partial ^2 f}{\partial x_1 ^2} & \frac{\partial ^2 f}{\partial x_1 \partial x_2} & \cdots & \frac{\partial ^2 f}{\partial x_1 \partial x_n} \\
\frac{\partial ^2 f}{\partial x_2 \partial x_1} & \frac{\partial ^2 f}{\partial x_2 ^2} & \cdots & \frac{\partial ^2 f}{\partial x_2 \partial x_n} \\
\vdots & \vdots & \ddots & \vdots \\
\frac{\partial ^2 f}{\partial x_n \partial x_1} & \frac{\partial ^2 f}{\partial x_n \partial x_2} & \cdots & \frac{\partial ^2 f}{\partial x_n ^2}
\end{bmatrix} \tag{11}
$$
${\bf \delta x}$について3次以降の項を無視しヘシアン行列を用いると、$f({\bf x} + \delta{\bf x})$を${\bf x}$で展開した結果は式(12)のようになる。
$$
f({\bf x} + \delta{\bf x})=f({\bf x}) + \delta {\bf x} \cdot \nabla f+\frac{1}{2}\delta {\bf x}^T {\bf H}({\bf x}) \delta {\bf x} \tag{12}
$$
右辺を${\bf \delta x}$について微分し、極小値を与える${\bf \delta x}$を求める。*1
$$ \nabla f+\frac{1}{2}({\bf H({\bf x})}+{\bf H({\bf x}) ^T}){\bf \delta x} = {\bf 0} \tag{13} $$ ${\bf H}({\bf x})$が対称行列であり、${\bf H}({{\bf x}}) = {\bf H ^T}({{\bf x}})$であることを用いると $$ \nabla f+{\bf H({\bf x})}{\bf \delta x} = {\bf 0} \tag{14} $$ となる。よって${\bf \delta x}=-{\bf H ^{-1}({\bf x})} \nabla f$となるので、漸化式(5)は以下のようになる。 $$ {\bf x}_{k+1} = {\bf x}_{k} - {\bf H ^{-1}({\bf x}_{k})} \nabla f \tag{15} $$ 一次関数の時と同様にこの更新を繰り返すことで最小値にたどり着くためには、$f({\bf x})$が下に凸の関数であることが必要となるが、これは${\bf H}({\bf x})$が半正定値行列(positive semi-definite)であればよく、即ち $$ \delta {\bf x} ^T {\bf H}({\bf x}) \delta {\bf x} \geq 0 \tag{16} $$ であればよい。
一方でNewton法には問題点があり、ひとつは図5(b)のように局所最適解に陥りやすいという問題点がある。もう一つは、ヘシアン行列：$H({\bf x})$を求めるために計算量が大きくなってしまう点である。 非線形最小二乗法を解く場合において、後者の問題を解決したのがGauss-Newton(GN)法であり、前者の問題を解決したのがLevenberg Marquardt(LM)法である。
Gauss-Newton(GN)法
非線形最適化で扱う問題の多くは、2乗誤差の最小化つまり$f({\bf x})$が以下のようにあらわされる場合を扱うことが多い。 $$ f({\bf x})=\frac{1}{2}\sum_{i=i} ^{m} {r_i({\bf x})} ^2 \tag{17} $$ このとき勾配$\nabla f$は、 $$ \nabla f = \sum_{i=i} ^{m} {r_i({\bf x})}\frac{\partial r_i({\bf x}) }{\partial {\bf x}} \tag{18} $$ となる。ヘシアン行列：$H({\bf x})$は式(18)を更に微分することによって求めることができ、$H({\bf x})$の$(j,k)$成分を${\bf H_{jk}}$とすると、 $$ {\bf H}_{jk}=\sum_{i=i} ^{m} \left( \frac{\partial r_i }{\partial x_k} \frac{\partial r_i }{\partial x_j} + r_i \frac{\partial ^2 r_i }{\partial x_k {\partial x_j}} \tag{19} \right) $$ となる。更にここで、第2項の2階微分項を無視すると、 $$ {\bf H}_{jk}=\sum_{i=i} ^{m} \frac{\partial r_i }{\partial x_k} \frac{\partial r_i }{\partial x_j} \tag{20} $$ を得る。 ヤコビ行列${\bf J} \in {\bf R ^{m \times n}}$および残渣ベクトル${\bf r} \in {\bf R ^m}$を用いると、勾配とヘシアン行列：$H({\bf x})$はそれぞれ、 $$ \nabla f = {\bf J} ^T {\bf r} \tag{21} $$ $$ H({\bf x})={\bf J} ^T {\bf J} \tag{22} $$ となる。よって漸化式(5)は、 $$ {\bf x}_{k+1} = {\bf x}_{k} - {\bf ({\bf J} ^T {\bf J}) ^{-1}} {\bf J} ^T {\bf r} \tag{23} $$ となる。 導出過程であったとおり、GN法では2階微分項を無視するという近似を用いているため収束が保証されていない。最小値布付近で初期値を与えて、十分に小さい残差ベクトル${\bf r}$と非線形性が比較的穏やかである必要がある。
まとめ
非線形最適化についてその概要と各種法の基本と特徴についてまとめた。最急降下法とNewton法の違い、Newton法とGauss Newton法の違い、辺りがクリアに分かっていることが重要だと思う。
機械学習やコンピュータービジョンにおいて頻出の2乗誤差においては、Gauss Newton法よりもそれを改良したLevenberg Marquardt(LM)法の方が実用的には用いられる。LM法は、軽く触れた通り最急降下法とGauss Newton法を組み合わせ、双方の欠点を解決し、性能と収束性の両立を目指したものと言える。詳細はまた別のエントリーで紹介したいと思う。
