SVMを理解するための準備として、最適化理論の制約なし最小値問題と、その解き方(ニュートン法)について。以下ノートとペンが必要です。

ゴール

最適化問題の概要を知る
1次元制約なし最適化問題をニュートン法で解く
高次元の制約なし最適化問題をニュートン法で解く

モチベーション

最適化理論

これから数回に分けて**SVM(サポートベクターマシン)**と呼ばれる、機械学習のアルゴリズムについて¹説明します。SVMで何ができるのかは後々説明するとして、まずはSVMを理解するための準備として、最適化理論のいくつかの話題について説明します。

最適化理論とは、与えられた条件の下で最善の選択をしなければならない問題に対し、それを数理的に解決する手助けとなる手段について深く考察する理論です。最適化理論で有名なものは「遺伝的アルゴリズム」でしょうか。遺伝的アルゴリズムでできることは、例えばこーじ氏のこの動画をご覧ください。非常に面白いです。

「最善の選択をしなければならない問題」とは、どのようなものがあるでしょうか。実は、これまでにも何度か遭遇していました。

Cycle2で扱った多項式回帰では、以下のような問題を解きました。

「$E(w_0, w_1) = \frac12 \sum_{i=0}^{N} (y(w_0, w_1, x_i) - t_i)^2$を最小にする$\vec{w}$は何か？」

またCycle4で扱ったk-means法では、以下のような問題を解きました。

「$E(q_{nk}, \vec{\mu_k}) = \sum_{k=1}^K \sum_{n=1}^N q_{nk} || \vec{x_n} - \vec{\mu_k} ||^2$を最小にする$q_{nk}, \vec{\mu_k}$は何か？」

これらの問題は以下のように定式化できるでしょう：

「関数$f(\vec{x})$を最小にする最適な選択(つまり$\vec{x}$)は何か？」

このような問題は**最適化問題**と呼ばれます。最適化問題は経済学や物理学など様々な分野で遭遇します。最適化理論は、これら異なる分野であっても適用できる、汎用性の高い解決手段を提供してくれます。また双対定理をはじめとする、純粋に数学的にも興味深い現象をも見ることができます。

そして機械学習もまた、最適化問題によく遭遇する分野の一つであり、さらに言えば、これから説明しようとしているSVMはこの最適化理論をふんだんに用いて体系化されています。そのため、SVMを理解するには最適化理論について理解しておくことが必要です。

ちなみに、次のように考える人もいるかもしれません：「最小値を最大値と置き換えた問題は最適化問題とは呼ばないのか？」---それも最適化問題と呼びます。ただ問題「関数$f(x)$が最大となるような$x$を求めよ」は、$f(x)$にマイナスをかけると「関数$-f(x)$が最小となるような$x$」を求める問題へと変わり、さらに二つの問題の答えが一致することがわかります。このような理由で、"最小値"に関する最適化問題だけ考えていても一般性は失われません。

最適化問題の種類

最適化問題には制約なし問題と制約あり問題があります。

制約あり問題とは以下のようなものです：

「関数$f(\vec{x})$を最小にする$\vec{x}$を求めよ。ただし、$\vec{x}$は$g(\vec{x}) \leq 0$を満たさねばならない。」

「ただし、$\vec{x}$は$g(\vec{x}) \leq 0$を満たさねばならない。」の部分がこの問題の制約条件になっています。制約条件が一つもない問題を制約なし問題と呼びます。

具体例をみます。例えば以下のようなシンプルな制約あり最適化問題を考えます。

「$f(x) = x^2$を最小にする$x$を求めよ。ただし$x$は$-(x-2) \leq 0$を満たさねばならない」

制約条件を無視すれば、$f(x)$を最小にする$x$は$0$になります。これは$\frac{df}{dx} = 0$を解けばでてきます。一方制約条件を考慮すると、$f(x)$を最小にする$x$は$2$になります。(★)

Cycle2とCycle4で遭遇した最適化問題は制約なし問題となる一方、SVMは制約あり問題と解くことになります。また上の例が示唆するように、制約あり問題のほうが考慮すべきことが多いため、一般的に難しいです。

最適化問題の解き方

最適化問題の解き方(計算アルゴリズム)は数多く提案されています。制約あり/なし問題それぞれで豊富にあり、また考えている関数によっても解き方が異なります。Wikipediaの記事をみると沢山あることがわかります。

今回と次回のCycleでは制約なし、制約ありそれぞれの計算アルゴリズムを一つずつ紹介します。具体的には、制約なし問題ではニュートン法、制約あり問題では内点法を扱います。SVMは制約あり問題なので内点法を使いますが、内点法の内部ではニュートン法を用いるため、いずれも説明します。

今回やること

今回はニュートン法を用いて1次元、また高次元の制約なし最適化問題を解きます。

確認

(★)を計算して確かめてください。

準備

これから説明に用いる表現や、数学的知識の準備をします。

最適化問題の表現

最適化問題を今後以下のように表現します。

\begin{align}
\min_{\vec{x}} \ & f(\vec{x}) \\
\text{s.t.} \ & g_1(\vec{x}) \leq 0, \\
& g_2(\vec{x}) \leq 0, \\
& \vdots \\
& g_M(\vec{x})\leq 0 . 
\end{align}

これは「$f(\vec{x})$を最小にする$\vec{x}$を求めよ、ただし$\vec{x}$は$g_i(\vec{x}) \leq 0, , i = 1, 2, \cdots, M$を満たさねばならない」という意味です。$\text{s.t.}$以下が制約条件です($\text{s.t.}$は"such that"の略です)。制約条件は必ず0以下となるように問題を設定します。例えば制約条件が$2 \leq x$だったら、式変形して無理やり$-(x - 2) \leq 0$と表現を変えます。

十分条件と必要条件

十分条件、必要条件という言葉の意味を確認してください。以下の動画で簡単に学べます。

【高校　数学Ⅰ】　数と式７４　十分・必要条件　（１５分）

扱う関数

以後扱う関数は十分に滑らかな(つまり何回でも微分可能な)関数とします。

勾配ベクトル

勾配ベクトルの意味を理解してください。以下の動画が大変丁寧に勾配ベクトルについて説明してくれています。尚、記号の表記ですが、$\nabla f(\vec{x})$と$\frac{\partial f(\vec{x})}{\partial \vec{x}}$は同じ勾配ベクトルを意味しています。表記は統一せず、これらの記号を混合して使います。

とある八雲の科学解説『勾配（gradient）を工作で知ろう』

確認

以下の関数の勾配ベクトルを求めてください。

\begin{align}
\text{(1)} & f(x_1, x_2) = x_1^2 + x_2^2 \\
\text{(2)} & f(\vec{x}) = \vec{a}^T \vec{x}, \\
& \text{ただし} \  \vec{a} = (a_1, a_2, \cdots, a_N)^T, \ \vec{x} = (x_1, x_2, \cdots, x_N)^T.
\end{align}

1次元制約なし最適化問題をニュートン法で解く

説明

以下のような制約なし最適化問題を考えます：

\min_x f(x)

関数$f(x)$が1変数しかとらないため1次元制約なし問題と呼びました。

まず、そもそもですが、この最適化問題には解がないかもしれません。例えば$f(x) = x^3$を考えてみましょう。$x$を負の方向に動かしていった場合、$f(x)$は$-\infty$に近づいていき、いくらでも小さな値をとることができます。すなわちこの問題には解が存在しません。(★)

また解が存在したしても、複数存在するかもしれません。例えば$f(x) = (x+1)^2 (x - 1)^2$を考えてみます。これの最小値は0ですが、解は-1と1の二つ存在します。(★)

以後、考えるすべての最適化問題は解を持ち、かつ解の個数が1つであるものとします。具体的には考える関数$f$はすべて狭義凸関数であることを仮定します。凸関数のグラフイメージはグーグル画像検索をしてみるとおよそわかります(画像検索)。狭義凸関数であれば、$f'(x)$は狭義単調増加となるため、この関数の最適化問題はただ一つの解を持ちます(証明は略²)。さらに、実は$f'(x) = 0$が、$x$がこの最適化問題の解であることの必要十分条件となります。幸いSVMで考える最適化問題は狭義凸関数です。

関数$f$が狭義凸関数であることを仮定したので、この最適化問題の解を求めるには以下の方程式を解けばよいことがわかります：

f'(x) = 0

以下$g(x) = f'(x)$とします。

代数的に(つまり中学高校で学んだ左辺から右辺に移項して～という解き方)で解いてももちろん構いませんが、ここではあえて別の方法で解きます。その別の方法とはニュートン法です。

ニュートン法とは次のようなアルゴリズムを言います。

初期値$x_0$と終了条件$\delta > 0$を選ぶ。
$x_{n+1}$を次のように更新する：$x_{n+1} = x_n - \frac{g(x_n)}{g'(x_{n})}$
$|g(x_n)| \geq \delta$である限り2.を繰り返す。
$|g(x_n)| < \delta$ならば終了。$x_n$を返す。

ひとつ注釈として、ニュートン法によって得られた解$x_n$は、$g(x)=0$の厳密な解ではありません。$g(x)=0$の真の解を$\hat{x}$とするならば$\hat{x} \neq x_n$です。しかしながら、両者の値は限りなく近い状態($\hat{x} \approx x_n$)となっています。このように真の解に限りなく近い値を近似解と呼びます。厳密に正しい解ではなく近似解を得ることも、私たちは「最適化問題を解いた」ことと見なします。実際、最適化問題を解く実用的な場面では近似解で十分です。

ニュートン法は次に説明するアイディアから来ています。

まず$g(x_{n+1})$を$x_n$を中心とした1次のテイラー展開で近似します。

g(x_{n+1}) \approx g(x_n) + g'(x_n) (x - x_n).

ここで$g(x_{n+1})=0$となる$x_{n+1}$が知りたいので、1次近似した右辺もそれが成立しているとみなして、右辺$=0$の式をたてます。

g(x_n) + g'(x_n) (x_{n+1} - x_n) = 0.

これを$x_{n+1}$について解くと以下のようになります。

x_{n+1} = x_0 - \frac{g(x_0)}{g'(x_0)}.

ここで得た$x_{n+1}$はあくまで本来の関数$g(x)$を１次近似した式の解に過ぎず、ほとんどの場合本来の$g(x)=0$の真の解とは大分異なることでしょう。しかしながら、もとの$x_n$よりも真の解には近づいてはいます(●)。つまりこの操作を繰り返すことで、真の解にどんどん近づくというわけです³。

確認

(★)を記した２つの説明を計算して確かめてください。増減表を書いてグラフの概形をかくとわかります。
(●)に関して、$g(x)=x^2 - 1$というグラフで、$x_n$が更新されるたびに真の解に近づいている様子を描写してください。
ニュートン法を実装して、以下の最適化問題をニュートン法を使って解いてください。$\delta = 0.001$くらいにしておいてください。

\min_{x > 0} \ x^3 - 6x - \log x \

n次元制約なし最適化問題をニュートン法で解く

説明

ここでは$n$次元の制約なし最適化問題と、それを解くためのニュートン法を考えます。次元があがりますが、基本となるアイディアは１次元のときと同じです。

以下のような最適化問題を考えます。ただし、$\vec{x} = (x_1, x_2, \cdots, x_n)^T$とします。

\min_{\vec{x}} \ f(\vec{x}).

$f$を狭義凸関数としているので、やはり以下の(連立)方程式の解が求める答えとなります。

\nabla f(\vec{x}) = \vec{0}.

これは連立方程式となり、代数的に解くことが一般に難しくなります。ニュートン法はここで効果を発揮します。

以下、$G(\vec{x}) = \nabla f(\vec{x})$とします。$G(\vec{x})$はn次元のベクトルであることに注意してください。$G(\vec{x}) = (g_1 (\vec{x}), g_2 (\vec{x}), \cdots, g_n (\vec{x}))^T$とします。また$\partial G(\vec{x})$を$G(\vec{x})$のヤコビ行列とします：

\partial G(\vec{x}) = 
\left(
\begin{matrix}
\frac{\partial g_1 (\vec{x})}{\partial x_1} & \frac{\partial g_1 (\vec{x})}{\partial x_2} & \cdots & \frac{\partial g_1 (\vec{x})}{\partial x_n} \\
\frac{\partial g_2 (\vec{x})}{\partial x_1} & \frac{\partial g_2 (\vec{x})}{\partial x_2} & \cdots & \frac{\partial g_2 (\vec{x})}{\partial x_n} \\
\vdots & & \ddots & \vdots \\
\frac{\partial g_n (\vec{x})}{\partial x_1} & \frac{\partial g_n (\vec{x})}{\partial x_2} & \cdots & \frac{\partial g_n (\vec{x})}{\partial x_n}
\end{matrix}
\right).

n次元のニュートン法のアルゴリズムは以下のようになります。

初期値$\vec{x_0}$と終了条件$\delta > 0$を選ぶ。
$\vec{x_{n+1}}$を次のように更新する：$\vec{x_{n+1}} = \vec{x_n} - \partial G(\vec{x_n})^{-1} G(\vec{x_n})$
$|| G(\vec{x_n}) || \geq \delta$である限り2.を繰り返す。
$||G(\vec{x_n})|| < \delta$ならば終了。$\vec{x_n}$を返す。

この更新式の導出も、１次元のときと同じアイディアでやります。

まず$g_i(\vec{x_{n+1}})$を$\vec{x_n}$を中心とした１次のテイラー展開で近似します(★1)。

\begin{align}
g_i(\vec{x_{n+1}}) & \approx g_i(\vec{x_n}) + \nabla g_i(\vec{x_n})^T (\vec{x_{n+1}} - \vec{x_n}) \\
& = g_i(\vec{x_n}) + \sum_{j=1}^{n} \frac{\partial g_i(\vec{x_n})}{\partial x_j} (x_{n+1, \, j} - x_{n, \, j}),
\end{align}

ただし、$\vec{x_k} = (x_{k, , 1}, x_{k, , 2}, \cdots, x_{k, , n})$としました。

ここで１次元のときと同様、右辺$=0$を考え、さらに、それをすべての$i$に対して式を立て、連立方程式として記述すると、以下のようになります。(★2)

G(\vec{x_n}) + \partial G(\vec{x_n}) (\vec{x_{n+1}} - \vec{x_n}) = \vec{0}.

これを$\vec{x_{n+1}}$について解けば以下のようになります。(★3)

\vec{x_{n+1}} = \vec{x_n} - \partial G(\vec{x_n})^{-1} G(\vec{x_n}).

このようにして$\vec{x_{n}}$を更新してゆくと、$G(\vec{x}) = \vec{0}$の解に近づいていきます。(これもやはりすべての考えている関数が凸であることから、近似解を得られることが保証されます。)

確認

(★1)にある多変数関数におけるテイラー展開を手元の数学書で確認してください。(下の補足にも書きました。)
(★2)と(★3)を計算して確かめてください。

実践問題

★ n次元のニュートン法を実装し、Cycle2の最適化問題をニュートン法で解いてください。Cycle2のときに出した$\vec{w}$と近いベクトルがでてくれば成功です。

参考

書籍

補足

多変数のテイラー展開の思い出し方

1次しか書いていないけど、同じ要領で何次でもできる。

$g(t) = f(\vec{x} + t \vec{y})$とする、ただし$\vec{x} = (x_1, x_2, \cdots, x_n)^T, \ \vec{y} = (y_1, y_2, \cdots, y_n)^T$.

$g(t)$の1次のマクローリン展開を考える。

g(t) = g(0) + g'(0)t.

ここで$g(0)$と$g'(0)$をそれぞれ計算する。$g'(0)$は連鎖律を用いる。

$g(0)$は

g(0) = f(\vec{x})

$g'(0)$は$\vec{u} = \vec{x} + t \vec{y}$, $\vec{u} = (u_1, u_2, \cdots, u_n )^T$とおくと、下の補足の補足を考慮して計算していくと、

\begin{align}
g'(t) &= \sum_{i=1}^{n} \frac{\partial f(\vec{u})}{\partial u_i} \cdot \frac{d}{dt} (x_i + y_it) \\
&= \sum_{i=1}^{n} \frac{\partial f(\vec{x} + \vec{y}t)}{\partial x_i} y_i .\\
g'(0) &= \sum_{i=1}^{n} \frac{\partial f(\vec{x})}{\partial x_i} y_i \\
&= \nabla f(\vec{x})^T \vec{y}.
\end{align}

なので、

g(t) = f(\vec{x}) + \nabla f(\vec{x})^T \vec{y} t.

$t=1$を代入すると、

\begin{align}
f(\vec{x} + \vec{y}) &= g(1) \\
&= f(\vec{x}) + \nabla f(\vec{x})^T \vec{y} .
\end{align}

$\vec{x} = \vec{x_n}$, $\vec{y} = \vec{x_{n+1}} - \vec{x_n}$とおきなおせば、

f(\vec{x_{n+1}}) = f(\vec{x_n}) + \nabla f(\vec{x_n})^T (\vec{x_{n+1}} - \vec{x_n}).

補足の補足

$g(x) = f(x+a)$, $u=x+a$とすると、

\frac{df(u)}{du} = \frac{g(x)}{dx} = \frac{df(x+a)}{dx}.

なぜなら

\begin{align}
\frac{df(u)}{du} &= \lim_{h \rightarrow 0} \frac{f(u+h) - f(u)}{h} \\
&= \lim_{h \rightarrow 0} \frac{f(x + a + h)- f(x + a)}{h} \\
&= \lim_{h \rightarrow 0} \frac{g(x + h) - g(x)}{h} \\
&= \frac{dg(x)}{dx} \\
&= \frac{df(x+a)}{dx}.
\end{align}

具体的にはSV分類のアルゴリズムについて。 ↩
杉浦解析入門Iを参照 ↩
これも杉浦解析入門I参照。今$g(x)$も凸関数を仮定していることに注意 ↩