Terasaki's blog

勉強したことをまとめたブログ

Newton法

今回は最適化手法の1つであるNewton法の原理と、Newton法の収束速度について説明したいと思います。

1変数関数 $f(x)$ の最小化問題を考えます。まず、 $\bar{x}+\Delta x$ まわりで関数 $f(x)$ を次のようにTaylor展開します。
$$ f(\bar{x}+\Delta x) = f(\bar{x}) + f'( \bar{x} ) \Delta x + \frac{1}{2}f''(\bar{x}) \Delta x^2 + \cdots $$

ここで、 $\Delta x$ に関する3次以上の項を十分小さいとみなして無視すると、
$$ f(\bar{x}+\Delta x) \approx f(\bar{x}) + f'(\bar{x})\Delta x + \frac{1}{2}f''(\bar{x})\Delta x^2 $$

この式は $\Delta x$ に関する2次関数になっています。そこで、極値を求めるために上式を $\Delta x$ に関して微分し、 $0$ とおきます。
$$ \begin{aligned} f'(\bar{x}) + f''(\bar{x})\Delta x & = 0 \\ \\ \Delta x & = -\frac{f'(\bar{x})}{f''(\bar{x})} \end{aligned} $$

したがって、初期値として $x_0$ を与え、次式を使って $ x_1 , x_2 , \cdots $ を逐次求めることで、最適解を近似的に求めることが出来ます。

$$ x_{i+1} = x_{i}-\frac{f'(x_i)}{f''(x_i)} $$

この手法をNewton法といいます。実際のプログラムではある程度反復したら終了する必要がありますが、数列が収束したかどうかを判定する条件として例えば $ \|x_{i+1}-x_{i}\| < \delta $ などを用います( $\delta$ は小さな正の数)。

Newton法は多変数関数に拡張できます。 二階微分可能な$n$ 変数関数 $f(x)$ の勾配を $ \nabla f(x) $ 、 Hesse行列を $H(x)$とすると、更新式は

$$ x_{i+1} = x_{i} - H(x_i)^{-1} \nabla f(x_i) $$

と表されます。ただし、Hesse行列の逆行列を計算するのは非効率的であるため、$d_{i} \in \mathbb{R}^n $ についての連立方程式

$$ H(x_i)\ d_{i} = -\nabla f(x) $$

を解くことで探索方向 $d_{i}$ を求め、

$$ x_{i+1} = x_{i} + d_{i} $$

で解を更新します。

Newton法は2次収束するという重要な性質を持っています。2次収束とは、 $i+1$ 回目の誤差が $i$ 回目の誤差の2乗に比例して減少していくような収束のことです。すなわち、 $x^*$ に収束するある点列 $\{x_i\}$ が次式を満たすとき、その点列は2次収束するといいます。

$$ \| x_{k+1} - x^* \| \leq M \| x_{k} - x^* \|^2 $$

ここで $ M $ は実数です。例えば $i$ ステップ目で誤差が $0.1$ のとき、 $i+1$ ステップ目ではほぼ $0.01$ 、その次のステップではほぼ $0.0001$ , $\cdots $ というように、非常に速く収束していきます。

証明

真の解を $ x^* $ として $ i $ 回目の解と真の解の差を $ \varepsilon_i = x_i - x^*$ とおきます。Newton法の更新の式の両辺から $x^*$ を減じると、
$$ \varepsilon_{i+1} = \varepsilon_{i}-\frac{f'(x_i)}{f''(x_i)} $$

ここで、 $f'(x_i) = f'(x^*+\varepsilon_i)$ 、 $f''(x_i) = f''( x^* + \varepsilon_i) $ のそれぞれのTaylor展開は、 $f'(x^*)=0$ であることを用いると、

$$ \begin{aligned} f'(x_i) & = f''(x^*)\varepsilon_i +\frac{1}{2}f'''(x^*)\varepsilon_i^2 + O\left(\varepsilon_i^3\right) \approx f''(x^*)\varepsilon_i +\frac{1}{2}f'''(x^*)\varepsilon_i^2 \\ \\ f''(x_i) & = f''(x^*) +f'''(x^*)\varepsilon_i + O\left(\varepsilon_i^2\right) \approx f''(x^*) +f'''(x^*)\varepsilon_i \end{aligned} $$

のようになります。この2式を誤差の漸化式に代入すると、
$$ \varepsilon_{i+1} \approx \varepsilon_{i}-\frac{f''(x^*)\varepsilon_i +\frac{1}{2}f'''(x^*)\varepsilon_i^2}{f''(x^*) +f'''(x^*)\varepsilon_i} $$

が得られます。右辺を通分し、分母を $f''(x^*) +f'''(x^*)\varepsilon_i \approx f''(x^*)$ と近似すると、最終的に以下の式が得られます。

$$ \varepsilon_{i+1} \approx \frac{f'''(x^*)}{2f''(x^*)}\varepsilon_i^2 $$

この式は、 $i+1$ 回目の誤差が $i$ 回目の誤差の2乗に比例することを示しています。したがって、Newton法が2次収束することが導かれました。

Kamada-Kawaiのアルゴリズム

グラフが与えられたときに、それを綺麗に描画するアルゴリズムグラフ描画アルゴリズムといいます。今回はグラフ描画アルゴリズムの一つである、Kamada-Kawaiのアルゴリズムについて説明したいと思います。

Kamada-Kawaiのアルゴリズムは、辺の重みを理想的な2点間のユークリッド距離としてエネルギー関数 $E$ を定義し、それを最小にするような節点の座標 $p_i = (x_i,y_i)^T\ (i=1,\cdots,n)$ を求めることで、綺麗なグラフを描画します。

まず、 $d_{i,j}$ を2点間の最短経路長とします。この $d_{i,j}$ を求めるアルゴリズムとして、例えばFloyd-Warshall法などがあります。次に$d_{i,j}$を用いて2点間の理想的なバネの長さ $l_{i,j}$ を定義します。
$$ l_{i,j} = Ld_{i,j} $$

$ L $は次式で計算されます。

$$ L = \max_{ i < j } \frac{ L_0 }{ d_{i,j} } $$

$L_0$は表示する画面の1辺の長さです。

また、バネの強さを次式で定義します。( $K$ は定数 )

$$ k_{i,j} = \frac{K}{d_{i,j}} $$

エネルギー関数の最小化

これで各値が求まったので、エネルギー関数 $E$ を定義し、各点の座標 $p_i=(x_i,y_i)^T $ を求めていきます。

エネルギー関数 $E$ を次式で定義します。

$$ \begin{aligned} E & \equiv \sum_{i=1}^{n-1}\sum_{j=i+1}^{n} \frac{1}{2}k_{i,j} \left( \|p_i-p_j\|-l_{i,j} \right)^2 \\ & = \sum_{i=1}^{n-1}\sum_{j=i+1}^{n} \frac{1}{2}k_{i,j} \left( (x_i-x_j)^2+(y_i-y_j)^2+l_{i,j}^2-2l_{i,j}\sqrt{(x_i-x_j)^2+(y_i-y_j)^2} \right) \end{aligned} $$

これを $p_i = (x_i,y_i)^T$ について最小にする手法はいくつかありますが、参考pdfと同じように、Newton-Raphson法で最小化したいと思います。
Newton-Raphson法では目的関数の1,2階導関数が必要なので、実際に計算をしていきます。まずは1階導関数です。

$$ \begin{aligned} \frac{\partial E}{\partial x_m} & = \sum_{i=1}^{n-1}\sum_{j=i+1}^{n} \frac{1}{2}k_{i,j} \left( \frac{\partial}{\partial x_m}(x_i-x_j)^2-2l_{i,j}\frac{\partial}{\partial x_m}\sqrt{(x_i-x_j)^2+(y_i-y_j)^2} \right) \\ & = \sum_{ i=1 }^{ m - 1 } \frac{1}{2} k_{i,j} \left( -2(x_i-x_m)-2 l_{i,m} \frac{ -( x_i - x_m ) }{ \sqrt{ (x_i - x_m )^2+( y_i - y_m )^2 } } \right) \\ & \quad +\sum_{j=m+1}^{n} \frac{1}{2}k_{i,j} \left( 2(x_m-x_j)-2l_{m,j}\frac{(x_m-x_j)}{\sqrt{(x_m-x_j)^2+(y_m-y_j)^2}} \right) \\ & = \sum_{i=1,i \neq m}^n k_{m,i}(x_m-x_i)\left( 1 - \frac{l_{m,i}}{\| p_m - p_i \|} \right) \end{aligned} $$

同様にして、$\frac{\partial E}{\partial y_m}$は、
$$ \frac{\partial E}{\partial y_m} = \sum_{i=1,i \neq m}^n k_{m,i}(y_m-y_i)\left( 1 - \frac{l_{m,i}}{\| p_m - p_i \|} \right) $$

となります。2階導関数は次式のようになります。

$$ \begin{aligned} \frac{\partial^2 E}{\partial x_m\partial y_m} & = \sum_{i=1,i \neq m}^n k_{m,i} \left(- l_{m,i} \frac{ 2( x_m - x_i )( y_m - y_i ) }{\| p_m - p_i \|^3}\times \left(-\frac{1}{2}\right) \right) \\ & = \sum_{i=1,i \neq m}^n k_{m,i}l_{m,i}\left( \frac{( x_m - x_i )( y_m - y_i )}{\| p_m - p_i \|^3 } \right) \\ \\ \frac{\partial^2 E}{\partial x_m^2} & = \sum_{i=1,i \neq m}^n k_{m,i} \left( 1 - l_{m,i} \frac{ \| p_m - p_i \| - \frac{ ( x_m - x_i )^2 }{ \| p_m - p_i \| } }{\| p_m - p_i \|^2} \right) \\ & = \sum_{i=1,i \neq m}^n k_{m,i}\left( 1 - \frac{ l_{m,i}( y_m - y_i )^2}{\| p_m - p_i \|^3 } \right) \\ \\ \frac{\partial^2 E}{\partial y_m^2} & = \sum_{i=1,i \neq m}^n k_{m,i} \left( 1- \frac{l_{m,i}(x_m-x_i)^2 }{\| p_m - p_i \|^3} \right) \end{aligned} $$

これで必要な式が求まりました。Newton-Raphson法では以下の式を繰り返し計算して近似的に極小値を求めます。 $$ \begin{aligned} H\Delta p_m & = - \nabla E \\ p_m & \leftarrow p_m + \Delta p_m \end{aligned} $$ ただし、$H$はヘッシアンといい、2変数関数の場合

$$ H = \left( \begin{array}{cc} \frac{\partial^2 E}{\partial x_m^2} & \frac{\partial^2 E}{\partial x_m\partial y_m} \\ \frac{\partial^2 E}{\partial y_m\partial x_m} & \frac{\partial^2 E}{\partial y_m^2} \end{array} \right) $$

で定義されます。 この反復計算によって求まった各節点の座標を用いることで、元々の各エッジの重みを考慮したグラフを描画することが出来ます。このとき、反復計算における終了条件として、次の値を使用します。 $$ \Delta_i = \sqrt{\left( \frac{\partial E}{\partial x_m} \right)^2 + \left( \frac{\partial E}{\partial y_m} \right)^2} $$ 擬似コードは次のようになります。
  1. Floyd-Warshall法で$d_{i,j}$を求める
  2. $l_{i,j}$を求める
  3. $k_{i,j}$を求める
  4. $p_i = (x_i,y_i)^T$をランダムに初期化する
  5. $\mathrm{while}(\max_{i} {\Delta_i} > \varepsilon)$
    1. $m \leftarrow $ $\Delta_i$を最大化する$i$
    2. $\mathrm{while}(\Delta_m > \varepsilon)$
      1. $ H\Delta p_m = - \nabla E $ から$ \Delta p_m $を求める
      2. $ p_m \leftarrow p_m + \Delta p_m $ で座標を更新する

プログラム例

今回はC++で実装しました。

実行結果

下のような無向グラフに対してプログラムを実行させました。結果はgnuplotで表示しています。
5 10
0 1 1
0 2 1
0 3 1
0 4 1
1 2 1
1 3 1
1 4 1
2 3 1
2 4 1
3 4 1
ただし、このデータは最初の1行が 節点数$n$、辺の本数$ m $を表し、以降の$ m $行が始点$ a_i $ 終点$ b_i $ 重み $ c_i $を表しています。このプログラムでは、コマンドライン引数で無向グラフか有向グラフかを指定できるようになっています。 f:id:tera_kun:20170519005007p:plain
今回は$L_0=1,K=10$としました。結果を見ると、各辺の重みを考慮したグラフが描画できていることがわかります。

まとめ

このように、Kamada-Kawaiアルゴリズムは比較的単純な方法で良い節点座標を求めることが出来ます。しかし、注意すべきこととして、局所最適解の問題があります。Newton-Raphson法で求めた解は極小解なので、それよりもエネルギー関数が小さくなるような解が存在することもあります。したがって、あまり質の良くない結果を得ることがあり、何らかの工夫が必要です。例えば、Kamada-Kawaiアルゴリズムによって求めた解を、Fruchterman-Reingoldのアルゴリズムで改善する、という方法が挙げられます。

参考

以下のpdfを参考にしました。
https://cs.brown.edu/~rt/gdhandbook/chapters/force-directed.pdf

KaTeXメモ

ブログに$\KaTeX$を導入したので、使い方などについてのメモ。

$\KaTeX$は、Khan Academyが開発している、Webブラウザ上で動作する数式描画ライブラリです。有名な数式描画ライブラリとしてMathJaxがありますが、描画速度が若干遅いという欠点があります。それに対して$\KaTeX$は、高速に描画できることを売りにしているようです。

  • Fast: KaTeX renders its math synchronously and doesn’t need to reflow the page.
  • Print quality: KaTeX’s layout is based on Donald Knuth’s TeX, the gold standard for math typesetting.
  • Self contained: KaTeX has no dependencies and can easily be bundled with your website resources.
  • Server side rendering: KaTeX produces the same output regardless of browser or environment, so you can pre-render expressions using Node.js and send them as plain HTML.
KaTeX – The fastest math typesetting library for the web

はてなブログで$\KaTeX$を使用する

自分のブログの「詳細設定」「headに要素を追加」の部分に以下のコードを記述するだけで、$\KaTeX$を使用することが出来ます。

<link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/KaTeX/0.7.1/katex.min.css" integrity="sha384-wITovz90syo1dJWVh32uuETPVEtGigN07tkttEqPv+uR2SE/mbQcG7ATL28aI9H0" crossorigin="anonymous">
<script src="https://cdnjs.cloudflare.com/ajax/libs/KaTeX/0.7.1/katex.min.js" integrity="sha384-/y1Nn9+QQAipbNQWU65krzJralCnuOasHncUFXGkdwntGeSvQicrYkiUBwsgUqc1" crossorigin="anonymous"></script>
<script src="https://cdnjs.cloudflare.com/ajax/libs/KaTeX/0.7.1/contrib/auto-render.min.js"></script>
<script>// <![CDATA[
document.addEventListener("DOMContentLoaded", function(){
  renderMathInElement(
    document.body,{
      delimiters: [
        {left: "$$", right: "$$", display: true},
        {left: "$", right: "$", display: false}]})});
// ]]></script>

使い方

$で数式を囲むとインライン数式モードとなり、文章中に数式を表示できます。ディスプレイ数式モードで表示したいときは$$で囲みます。

$$ f(x) = x^2 $$

$$ f(x) = x^2 $$

これらの区切り文字を変えたいときは、headに記述したdelimitersの中身を書き換えます。

$\KaTeX$では、aligned環境を使うことができ、複数行の数式を書けます。例えば、

$$ \begin{aligned} f(x) & = x^2+2x+1 \\ & = (x+1)^2 \end{aligned} $$

と書くと、

$$ \begin{aligned} f(x) & = x^2+2x+1 \\ & = (x+1)^2 \end{aligned} $$

と表示されます。

array環境などを使うことで、行列も表示できます。

$$
A = \left( \begin{array}{cccc} 
a_{11} & a_{12} & \ldots & a_{1n} \\ 
a_{21} & a_{22} & \ldots & a_{2n} \\
\vdots & \vdots & \ddots & \vdots \\
a_{m1} & a_{m2} & \ldots & a_{mn} 
\end{array} \right) $$

$$ A = \left( \begin{array}{cccc} a_{11} & a_{12} & \ldots & a_{1n} \\ a_{21} & a_{22} & \ldots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{m1} & a_{m2} & \ldots & a_{mn} \end{array} \right) $$


なお、はてなブログでは改行すると<br />が挿入されてしまうので、LaTeXのコードを複数行にわたって記述する場合はpタグ停止記法を使う必要があります。

$\KaTeX$は開発途中であり、サポートしていない機能もあるようです。Function Support in KaTeX · Khan/KaTeX Wiki · GitHubに、現在使用可能な機能が書かれています。

また、はてなブログでKaTeXを使う - kivantium活動日記にもあるように、はてなブログキーワードリンク機能によって一部表示できない文字があるので、注意が必要です(一応スペースを入れて$ m $と書くことで$ m $と表示できるけど、使い勝手は悪そうです)。