読者です 読者をやめる 読者になる 読者になる

ゆるふわめも

in Kyoto or Tokyo

最適化:非線形計画について(最急降下、ニュートン、KKT条件)とNP困難問題に対する動的計画法

この記事の内容

 試験が近いからぱっとノートとかまとめた記事。細かい定義とかは省略してます。そして特に断りが無い限り{ x}はベクトルで、{ x_i }といった添字のあるものはベクトルではない単なる実数か定義域の内部にある数の一つです。あと微分可能性とかはほぼ触れてないので察してください。

参考文献

  • [1] 数理計画入門(1996) 福島雅夫著 朝倉書店


過去のメモ


非線形計画の最適化問題:最急降下法、ニュートン法、KKT条件まで - 雑なメモ
途中までは、上の記事とほぼ同じ内容。

このつぎのメモ(こっちの方が最新)


最適化:非線形計画+組み合わせ最適化のまとめのメモ - 雑なメモ

非線形計画問題の分類

問題の定式化

{
 目的関数: f(x) -> min \\
 制約条件: x \in S
}

一般的な非線形計画問題をとくための道具

以下でおおまかに分類する前に、どの問題でもたいがい使う道具とその性質。

勾配

{ \displaystyle {x} = (x_1 , x_2 , \dots , x_n)^T }について、このベクトル{ \displaystyle {x} }での(二次元での)傾きと呼ばれるものに対応する値は、

{ \displaystyle
\nabla f = \frac{d f({x})}{d {x}} = \left(
    \begin{array}{c}
      \frac{d f({x})}{d {x_1}} \\
      \frac{d f({x})}{d {x_2}} \\
      \vdots \\
      \frac{d f({x})}{d {x_n}} \\
    \end{array}
  \right)
}
であり勾配と呼ばれる。

ヘッセ行列

{ \displaystyle
H = \Delta f = \nabla^2 f = \frac{d^2 f({x})}{d {x}^2} = \left(
    \begin{array}{c}
      \frac{d^2 f({x})}{d {x_1}{x_1}} \dots \frac{d^2 f({x})}{d {x_1}{x_n}}\\
      \frac{d^2 f({x})}{d {x_2}{x_1}} \dots \frac{d^2 f({x})}{d {x_2}{x_n}}\\
      \vdots \\
      \frac{d^2 f({x})}{d {x_n}{x_1}} \dots \frac{d^2 f({x})}{d {x_n}{x_n}}\\
    \end{array}
  \right)
}

と記述されて、これをヘッセ行列と呼ぶ。つまり{ H_{i,j} = \frac{d^2 f({x})}{d {x_i}{x_j}} }である。ここで偏微分{ \frac{d^2 f({x})}{d {x_i}{x_j}} 、\frac{d^2 f({x})}{d {x_j}{x_i}} }がともに連続だとするとこの二式は等しい。

勾配とヘッセ行列の性質
  • 勾配は上で記述されているようにベクトルであり、その向きは関数がもっとも急に増加する方向を向いているから、その逆の向きは関数がもっとも急に減少する方向。そしてこの性質は最急勾配法で使うはず。
  • 目的関数のヘッセ行列について、目的関数が凸関数ならば任意の点でヘッセ行列は半正定値行列。(半正定値行列、正定値行列に付いては以下を参照)

行列の定値性 - Wikipedia

問題の分類

すべての問題は目的関数に-1をかけることで最小化問題に書き換えることができる。なので、問題の種類を分類分けするとしたら制約条件のS、実行可能領域によって分類。

制約無問題

 つまり S が { R^n }全体の場合で目的関数にとくに制限はない。以下のような性質がある。

  • 最適解{x^*}では{\nabla f(x^*) = 0}が成り立つ。

 これは{x^*}が局所最適解であるための必要条件。

  • 最適解{x^*}では{\nabla^2 f(x^*)}は半正定値行列である。

 この証明は参考文献の第四章が詳しい。
証明:
 背理法で示す。「最適解{x^*}では{\nabla^2 f(x^*)}は半正定値行列である。」が成り立たないと仮定する。
ここで、目的関数fを任意の点{x'}テイラー展開すると、

{
f(x) \\
= f(x') + \sum_{k=1}^n \frac{d f(x')}{d x_i} (x_i - x'_i) + \frac{1}{2} \sum_{k=1}^n \sum_{l=1}^n \frac{d f(x')}{d x_k} \frac{d f(x')}{d x_l} (x_k - x'_k)(x_l - x'_l) + (三次以上の項)\\
= f(x') + \nabla f(x')^T (x - x') + \frac{1}{2} (x - x')^T \nabla^2 f(x') (x - x') + (三次以上の項)\\
}

と記述できる。今、任意の点{x}は最適解{x^*}と適当な実数{\alpha}とベクトル{d}を用いて{x = x^* + \alpha d}とかける。今点{x}を最適解{x^*}にてテイラー展開すると、上の式に代入して

{
f(x^* + \alpha d) \\
= f(x^*) + \nabla f(x^*)^T (x^* + \alpha d - x^*) + \frac{1}{2} (x^* + \alpha d - x^*)^T \nabla^2 f(x^*) (x^* + \alpha d - x^*) + (三次以上の項)\\
= f(x^*) + \frac{1}{2} \alpha^2 d^T \nabla^2 f(x^*) d + O(\alpha^3)\\
}

とかける。「最適解{x^*}では{\nabla^2 f(x^*)}は半正定値行列である。」が成り立たないと仮定しているのだから、
{\alpha^2 d^T \nabla^2 f(x^*) d < 0}となるようなベクトル{d}が存在するはずである。ここで上のテイラー展開の結果でαの値を十分小さくすると{O(\alpha^3)}の項を無視できる。
すると{\alpha^2 d^T \nabla^2 f(x^*) d < 0}なのだから、

{
f(x^* + \alpha d) = f(x^*) + \frac{1}{2} \alpha^2 d^T \nabla^2 f(x^*) d < f(x^*)
}

であるがこれは最適解が{x^*}であることに矛盾している。■

凸計画問題

 Sが凸領域であり、さらに目的関数fが凸関数であるならばその問題は凸計画問題と呼ぶ。つまり、

  1. {x,y \in S , 0 \leq \alpha \leq 1}ならば{ \alpha x + (1 - \alpha)y \in S}が成り立ち (実行可能領域が凸領域であり)
  2. {x,y \in S , 0 \leq \alpha \leq 1}ならば{f(\alpha x + (1 - \alpha)y) \leq \alpha f(x) + (1-\alpha)f(y) }が成り立つ (目的関数が凸関数)

の場合がこれにあたる。
凸最適化 - Wikipedia

この問題は以下の性質をもつ。

  • 局所最適解が見つかったならそれは大域最適解である

証明:背理法で示す。局所最適解と大域最適解が一致しないとして局所最適解を{x^*}、大域最適解を{y^*}とする。すると目的関数は凸関数なので
{f(\alpha x^* + (1 - \alpha)y^*) \leq \alpha f(x^*) + (1-\alpha)f(y^*) }
が成り立ち、実行可能領域は凸領域だから
{x^*,y^* \in S , 0 \leq \alpha \leq 1}ならば{ \alpha x^* + (1 - \alpha)y^* \in S}
が成り立つ。よって最適解{ \alpha x^* + (1 - \alpha)y^* }が見つかったがこれは大域最適解が{y^*}であることに矛盾する。

  • {m+n}数の制約条件を以下のように定義する。

{ c_1(x) = 0 , \dots , c_m(x) = 0}
{ d_1(x) \leq 0, \dots , d_n(x) \leq 0}
このとき、 { d_1(x) \leq 0, \dots , d_n(x) \leq 0}がすべて凸関数ならば実行可能領域は凸領域になる。

非線形計画問題の解の発見法

ここでも、実行可能領域の種類(制約なし、制約あり)の場合で分けて考える。

制約無問題

  1. 最急勾配法
  2. ニュートン法
  3. 準ニュートン法

が授業で触れた範囲。以下、{x_i}はベクトル。

最急降下法 - Wikipedia
ニュートン法 - Wikipedia
準ニュートン法 - Wikipedia

  1. 適当な点{x_0}を決める。
  2. {x_k + \alpha (\nabla f(x_k))}が最小となるような{\alpha}を求める。
  3. {x_{k+1} = x_k + \alpha (\nabla f(x_k))}とする。
  4. {|\nabla f(x_{k+1})| < \epsilon }(事前に定めた閾値を下回る)なら終了、そうでない場合はステップ2へ


今、点{x_0}から最急降下法で到達しうる範囲では常にヘッセ行列が正定値行列であるとする(つまり、解に近い適当な点をとったらそこが凸領域で、そこでは目的関数も凸関数であるように?)。このときこの領域では{x^T \nabla^2 f(x) x}は常に正の値だからこれを用いてノルム{ ||x||_G = \sqrt{x^T G x} }を定義すると条件数{\tau}を用いて、

{

 ||x_{k+1} - x^*||_G \leq (\frac{\tau - 1}{\tau + 1})||x_{k} - x^*||_G

}

が成り立ち、ある比率{(\frac{\tau - 1}{\tau + 1})}に従って収束するから一次収束と言う。ここで条件数は

条件数(じょうけんすう、英: condition number)は、問題のコンピュータでの数値解析しやすさの尺度であり、その問題がどれだけ数値解析に適しているかを表す。条件数が小さい問題は「良条件 (well-conditioned)」であり、条件数が大きい問題は「悪条件 (ill-conditioned)」である。

条件数 - Wikipedia

問題を振り返ると、悪条件 (ill-conditioned)」ならば比率{(\frac{\tau - 1}{\tau + 1})}は1に近づくから収束が遅いと考えられる。


  1. 適当な点{x_0}を決める。
  2. {|\nabla f(x_{k+1})| < \epsilon }(事前に定めた閾値を下回る)なら終了。
  3. {x_{k+1} = x_k - \frac{\nabla f(x_k)}{\nabla^2 f(x_k)} }とする。
  4. ステップ2へ。


目的関数をテイラー展開すれば
{
f(x) = f(x_k) + \nabla f(x_k)^T (x - x_k) + \frac{1}{2} (x - x_k)^T \nabla^2 f(x_k) (x - x_k) + (三次以上の項)\\
}
となることが上の方に書かれていた(気がする)。
今、この式の最後の項を無視して
{
f(x) \simeq f(x_k) + \nabla f(x_k)^T (x - x_k) + \frac{1}{2} (x - x_k)^T \nabla^2 f(x_k) (x - x_k) 
}
とする。今、点{x_0}からニュートン法で到達しうる範囲では常にヘッセ行列が正定値行列であるとする(さっきと同じ)。
そして、このテイラー展開した式がなるべく小さくなうように{(x - x_k)}を定めて一番小さくなる点を{x_{k+1}}にすればいい。だから上の式を両辺{(x - x_k)}について微分して、
{

\frac{d f(x)}{d(x - x_k)} = \nabla f(x_k) + \nabla^2 f(x_k) (x - x_k) = 0

}
となる点を探す。なぜなら、「常にヘッセ行列が正定値行列であるとする」のだから「凸関数」になりそのときは「局所最適解は大域最適解に一致」するから。この式から{ x_{k+1} = - (\nabla^2 f(x_k))^{-1} \nabla f(x_k) }が得られる。
 ニュートン法の収束が二次収束であることを示すためには(参考文献[1]p111)やっぱり目的関数を{x_k}まわりでテイラー展開する。さらに、ニュートン法では{\frac{\nabla f(x_k)}{\nabla^2 f(x_k)}}があるのでこの部分に関係する項が出てくるように{\nabla f(x)}も最適解でテイラー展開する。
{
f(x^*) = f(x_k) + \nabla f(x_k)^T (x^* - x_k) + \frac{1}{2} (x^* - x_k)^T \nabla^2 f(x_k) (x^* - x_k) + O<a href= この式を整理すると が二つめのテイラー展開から導けて二次収束とわかる。 準ニュートン法(参考文献[1]p111より) 適当な点を決める。 適当な正定値対称行列を示す行列とする。 (事前に定めた閾値を下回る)なら終了。 {x_k + \alpha (B_k)^{-1} (\nabla f(x_k*1}"/>が最小となるような{\alpha}を求める。

  • {x_{k+1} = x_k + \alpha (B_k)^{-1} (\nabla f(x_k))}
  • {b_{k+1} (x_{k+1} - x_k) = \nabla f(x_{k+1}) - \nabla f(x_k)} とする。
  • ステップ2へ。

  • ニュートン法は一ステップごとにヘッセ行列を計算するからすごい時間がかかるので、準ニュートン法ではヘッセ行列を近似した代替をもちいていてその働きをしているのが行列B。なのでこのBをうまく近似すれば性能も良くなる。

    a.普通の方法
    普通の微分の定義式を思い出して
    {
 \nabla^2 f(x) \simeq = \frac{\nabla x_{k+1} - \nabla x_{k}}{x_{k+1} - x_{k}}
}

    b,BFGS法
    NTTデータシステムさんのページにすばらしい記事がありました。以下参照。

    BFGS 法 - 数理計画用語集

    ルーシュ・キューン・タッカー条件(KKT条件)

    カルーシュ・キューン・タッカー条件 - Wikipedia


    KKT 条件 - 数理計画用語集

    上記のページに詳しく記述されている(ただし出典はこの記事の参考文献におなじらしい)。

    問題の定式化

    { 目的関数: f(x) -> min }
    { c_1(x) = 0 , \dots , c_m(x) = 0}
    { d_1(x) \leq 0, \dots , d_n(x) \leq 0}

    制約なし問題では、{\nabla f(x) = 0}が最適であるための必要条件のひとつだが、制約ありの場合は{\nabla f(x) = 0}ではなくても実行可能領域の境界上などで最適解がみつかる場合が多いから{\nabla f(x)}という情報だけでは最適解は見つからない。

    そこで、最適解である点がいくつもの制約を満たすために制約を満たすような方向へ”引っ張られて釣り合った結果”の点であると考えたら、

    {
 \nabla f(x^*) +( u_1 \nabla c_1(x^*) + \dots + u_{m} \nabla c_{m} (x^*) ) + ( u_{m+1} \nabla d_1(x^*) + \dots + u_{m} \nabla d_{m+n} (x^*) ) = \vec{0} 
}

    となるような非負の実数{u_1,\dots,u_{(m+n)}}が存在すると言える(証明は参考文献に書いてないし実際学部の範囲ではなさそうなので考え方だけ理解)。{\nabla}をつけるのは勾配は関数が最も増加する方向の方向ベクトルなのだから、最適解において各関数(制約)がもっとも大きく引っ張る方向と個人的に解釈、釣り合うのだからその和は零ベクトルで表される。

    そして、有効制約でない制約(最適解である{x^*}を選んだ時に{c_k (x^*) = 0}とならない制約)はこの文脈で言えば最適解の点を制約がなんとか満たすように"引っ張っていない(まだいくらか動ける余裕がある)"のだから

    {
 \nabla f(x^*) +( u_1 \nabla c_1(x^*) + \dots + u_{m} \nabla c_{m} (x^*) ) + ( u_{m+1} \nabla d_1(x^*) + \dots + u_{m} \nabla d_{m+n} (x^*) ) = \vec{0} 
}

    にて最適解を移動させるための役割を持たないはずだと考えると「有効制約でないものに対応する係数uは0である」と考えられる。参考文献にはもっとしっかりベクトルの一次従属などを用いて説明があった。

    そしてもしも制約{d_1(x),\dots,d_n(x)}が存在しないとすると、この"引っぱりあい"を示す式は

    {
 \nabla f(x^*) +( u_1 \nabla c_1(x^*) + \dots + u_{m} \nabla c_{m} (x^*) )  = \vec{0} 
}

    と書くことができる。これはまさにラグランジュ未定乗数法の式そのもので実際この最適化の文脈でもuをラグランジュ係数と呼ぶらしい。

    ラグランジェの未定乗数法 [物理のかぎしっぽ]
    ラグランジュの未定乗数法 - Wikipedia

    KKT条件の定式化

     以上の考えをもって、KKT条件を以下に記述(参考文献[1]p121)。

    KKT条件
    { \displaystyle
  \left\{
    \begin{array}{l}
     \nabla f(x^*) +( u_1 \nabla c_1(x^*) + \dots + u_{m} \nabla c_{m} (x^*) ) + ( u_{m+1} \nabla d_1(x^*) + \dots + u_{m} \nabla d_{m+n} (x^*) ) = \vec{0} \\
\\
     u_k = \left\{
    \begin{array}{lll}
      u_k \geq 0 & if & c_k({x^*}) = 0 \\
      0 & if & d_k({x^*}) < 0 \\
      u_k \geq 0 & if & d_k({x^*}) \leq 0 \\
    \end{array}
  \right.
    \end{array}
  \right.
}

    最適解であるための必要・十分条件(参考文献[1]p123)

     KKT条件を満たすすべての点が最適解であるとは限らないのは適当に制約と目的関数を作ってためしてみると分かる。なのでさらに目的関数と制約が凸関数であるという制約を加えることで局所最適解が大域最適解と一致するようにするとKKT条件を満たす点が最適解であると言えるようになる。さらに以下の式Lを考える。

    {

L = f(x) + \sum_{k=1}^m u_k c_k (x) + \sum_{k=m+1}^n u_k d_k (x)

}

    このLの勾配は

    {

\nabla L = \nabla f(x) + \sum_{k=1}^m u_k \nabla c_k (x) + \sum_{k=m+1}^n u_k \nabla d_k (x)

}

    であり、まさにKKT条件の”引っ張り合い”の式であって

    {

\nabla^2 L = \nabla^2 f(x) + \sum_{k=1}^m u_k \nabla^2 c_k (x) + \sum_{k=m+1}^n u_k \nabla^2 d_k (x)

}

    がLのヘッセ行列。そして有効制約の勾配ベクトルと直行するベクトル{  y }{ y^T \nabla^2 L y \geq 0}が成り立つことが知られている。この式は以前の制約なしの凸計画問題の「ヘッセ行列が半正定値行列である」ことと対応しているとみると、見つけた点が最適解であるための十分条件は凸計画問題の「ヘッセ行列が正定値行列である」ことに対応して「{ y^T \nabla^2 L y > 0}」と見当がつく。

    ラグランジュの双対理論

    以下に詳しい記述がありました。

    • 社団法人 日本オペレーションズ・リサーチ学会 OR事典さまより参照

    双対性理論 - ORWiki

    講義では最後に触れましたが、試験と関係しているかは不明なのでいったん飛ばして時間の許す時に追記予定。

    NP困難な問題に対する解の求め方

    ここからは最適化の授業の後半のメモ、だけど図が多かったのでこの記事を(外部の人が)読んだだけではわからないかも。すいません。授業後半の主な話題は...

    だったはず(まだ講義は終わってないのでわからない)。

    たぶんこの章の範囲はアルゴリズム理論などの本が詳しい。あと

    • エージェントアプローチ 人工知能 スチュワート=ラッセル著

    この本の途中に局所探索やヒューリステリックについていろいろ記述してあるが本が網羅的でしかも高いので図書館で借りたほうがいいかも。プログラムコードはなく、基本疑似プログラミングと文章での記述。

    組み合わせ計画(参考文献[1]p141)

     これまでの問題は制約関数(連続値をとる)に基づいて最適解を見つけていたが、ここからは有限この要素のなかから最適となる値や組み合わせを探す問題を議論。有限個とあるように、制約や目的関数はもはや連続関数ではなく勾配が0になる点を解の候補にするといったアプローチは無理。しかし制約が有限個ということは全通りを計算すれば必ず解がみつかる(かないこと)が分かる。だけどこれは計算時間の観点から不可能に近い。そのため、「最適となる値や組み合わせを探す」ために「要素の中から解(の一部を構成するもの)を選ぶ」基準に基づいて早い段階で最適解(だと思われる)ものを見つけることを考える。

    貪欲法

    貪欲法 - Wikipedia


    最強最速アルゴリズマー養成講座:そのアルゴリズム、貪欲につき――貪欲法のススメ (1/3) - ITmedia エンタープライズ

    その名の通り、各時点で最も最適と思われる点を(先のことを考えずに)選ぶ方法。たとえば目的地に向かう途中の道で分かれ道に出会ったとして道には「次の分かれ道まで○×km」と書かれているとするとその道の次の分岐や方向を考えずに「最もつぎの分岐までの距離が小さい道」を選ぶ。貪欲法に分類されるアルゴリズムには

    1. クラスカルアルゴリズム
    2. ダイクストラ

    がある。以下、クラスカルアルゴリズムの疑似コードをwikipediaより引用。


    クラスカル法 - Wikipedia

    グラフの各頂点がそれぞれの木に属するように、森(木の集合) F を生成する(つまり頂点1個だけからなる木が頂点の個数だけ存在する)
    S ← グラフの全ての辺を含む集合
    while (S が空集合ではない)
    S から重みが最小の辺 e を削除する
    if (e が2つの木を連結するもの)
    e を森に加え、2つの木を1つにまとめる

    分枝限定法

     以下の説明を参照。要は全通り試さなくても明らかにだめだと思ったらその選択肢とそこから派生する組み合わせは試さずに他を試していく。


    分枝限定法とは - OR事典 Weblio辞書

    • 授業でやったこと:ナップザック問題に対する分枝限定法

    ナップサック問題 - Wikipedia
    以下、(試験での)解を求めるまでのステップ。

    1.問題を定式化
    {
 目的関数 : max  \sum_{1 \leq k \leq n} c_i x_i \\
 制約条件 : \sum_{1 \leq k \leq n} a_i x_i \leq b  かつ x_i \in \{ 0 , 1 \} (\forall k) 
}

    つまり、a=各アイテムの重さ、b=ナップザックに入る重さの限界、c=各アイテムの価値、x=ナップザックに入っているときのみ1、目的関数はナップザックに入るすべてのアイテムの価値の総和を最大にすることを示している。

    2.アイテムを(何らかの基準に基づき)並べ替え
    このステップはなくてもいい。だけどはじめのほうにナップザックに入れるアイテムは「最適解に含まれていそうなもの」がいいので、そういうものがはじめに入りそうな基準を何でもいいから作って並び替え。ここでは参考文献に基づいて「重さあたりの価値が最大のもの」から順にソートして、
    {
 \frac{c_1}{a_1} \leq \frac{c_2}{a_2} \leq \frac{c_3}{a_3} \leq \frac{c_4}{a_4} 
}
    であるとする。(簡単のためにアイテムは4個に限定)

    3.{x_i}に基づいた木を作成する
     はじめ、すべての{x_i = 0}とする、つまり何も入っていないことを示している。以下、”重みあたりの価値が高いアイテム”から順にナップザックに入れる、入れないことが試されていくけど文書だけでは説明できないので他参照。とりあえず分枝限定法で作られる木についての性質は

    • 木は{x_i }の値によって分かれていく二分木
    • 各ノードは以下の値を記憶する。
      • 自分より上のノードの{x_i}に対応して、現在の残りのナップザックの空きスペースb_{empty}
      • b_{empty}に基づいて
        • (b-1)b_{empty}を超えない値で最も”重みあたりの価値が高いアイテム”を入れた場合の目的関数{max  \sum_{1 \leq k \leq n} c_i x_i }の値を記憶
        • (b-2)b_{empty}が0でない場合はb_{empty}を超えて一個だけ入れることを許して最も”重みあたりの価値が高いアイテム”を入れた場合の目的関数{max  \sum_{1 \leq k \leq n} c_i x_i }の値を記憶


    こうすることで、(b-1)はその時点での「最適解は最低でもこの値以上」である境界を示し(b-2)はその時点での「最適解がこの値よりも上である可能性がない」境界が木に示される。なので、分岐点の(b-1)の値が分岐先の点の(b-2)の値よりも大きな値ならばその分岐をすすんでも最適解が見つかることは絶対ないから分岐をすすまない。

    動的計画法(ダイナミックプログラミング)

    動的計画法 - Wikipedia

    同じナップザック問題を動的計画法で解くことも可能。

    (途中です)

    {}
    {
}

    役立つページ

    *1:x^* - x_k)^3)\\ \nabla f(x^*) = \nabla f(x^*_k) + \nabla^2 f(x_k) ^T (x^* - x_k) + O((x^* - x_k)^2)\\ }"/>
    この式を整理すると
    {
 f(x^*) = f(x_k) + \nabla f(x_k)^T (x^* - x_k) + \frac{1}{2} (x^* - x_k)^T \nabla^2 f(x_k) (x^* - x_k) + O((x^* - x_k)^3) \\ 
 0 = (x_{k+1} - x^*) + O((x^* - x_k)^2)  
}
    が二つめのテイラー展開から導けて二次収束とわかる。


    1. 適当な点{x_0}を決める。
    2. 適当な正定値対称行列を示す行列{B_0 = (n,n)単位行列 }とする。
    3. {|\nabla f(x_k)|< \epsilon }(事前に定めた閾値を下回る)なら終了。