数学」カテゴリーアーカイブ

似ているgraphを探す方法

ある与えられたグラフに対して,全く同じ形をしているグラフを探す問題をグラフ同型判定問題という.数学の分野で同型という言葉が出てくるときは,一見同じようには見えない場合でも,見方(例えば記号の付け方)などを変えることで同じにできるような場合にそう呼ぶことが多い.この場合の同型も同じような意味合いである.グラフには大きく分けて有向グラフなのか無効グラフなのかという違いがあるが,ここでは話を簡単にするため無向グラフに限定する.実は,多くの実問題で計算量や問題設定の緩和を図るために有向グラフを無向グラフにすることは比較的ありふれた話である.したがって,以下の話は必ずしも有向グラフでは無駄というわけではない.

グラフの同型判定問題や,もしくはより難しい問題である,あるグラフの中に探したいグラフを一部として含むかどうかを探す部分グラフ同型判定問題はNPであることはわかっているものの,NP困難であるかどうかはわかっていない(Web上の記事によってNP完全だったりNP困難であったり,もしくはどちらかわかっていないと書かれていたり,実際のところどうなのかわからない).どちらにせよ,この問題を解くことが今の所容易ではないことは変わらない.多くの文脈(例えば化学式の検索など)で,このグラフ同型性判定問題を解決したいと考えられており,多くのアルゴリズムが今までに提案されてきた.

さて,私が今回,扱いたいのはこのグラフ同型性判定問題ではなく,より曖昧とした問題,すなわち2つのグラフの類似度を知りたいという問題である.

不完全マッチング(inexact matching)

今回扱いたいのは不完全マッチング(inexact matching)である.実はこの単語,論文中に出てきた表現であり,私の知る限り日本語訳を見たことがない.inexact matchingは直訳すれば不正確マッチングだが,不正確は日本語の意味合いを考えると,誤解を招きかねない.完全マッチングという別の用語はあるものの,不完全という用語がおそらく今回の意味ではピッタリのように感じられるので,ここでは便宜上,不完全マッチングと呼ぶことにする.

さて,不完全マッチングとは要するに,与えられた2つのグラフは同型ではないのだが,微妙に似ており,その上で2つのグラフ間の最も妥当なノードもしくはエッジの対応付けはどのようなものかという問題である.この対応付けがわかると,2つのグラフの間の距離のような概念まで導入することができ,一昔に流行ったSVNではカーネル法に帰着できると分類問題に帰着可能になるため,このグラフ間の距離を求めることにはビッグデータなどの応用も含めて,一定の需要があるらしいのだ.さて,不完全マッチングでは,この妥当な対応付けというところが味噌で,不完全マッチングを提供するアルゴリズムによって,妥当の意味が違ってくる.一番多く使われているのは,グラフ編集距離(graph edit distance)と呼ばれるものを用いるやり方である.これは文字列に対する編集距離に概念としては近い.グラフのノードもしくはエッジの追加や削除などを行って,片方のグラフをもう片方のグラフに近づけたときに,最も少ない手数で対応付けられたときのコストをもって,2つのグラフの距離とする.また,このときに対応付けられたノードもしくはエッジが今回の場合の答えになるというものだ.しかし,この最も最短のグラフ編集の手順を探す問題はNP困難であることが知られている.最も高速に厳密解を得る方法としては,A*アルゴリズムを用いる方法があるが,この方法ではたかだか15ノード程度が求められる関の山で,それ以上となると手も足も出ない.そこで,これまでに提案されたアルゴリズムでは,近似的にこの問題を解く方法が模索されてきた.

もう一つのこの問題の考え方は,2つのグラフのノード,エッジ同士の対応付に関する適応関数を定義し,その適応関数が最大となるような割当を探す問題に帰着させる方法である.この方法では,予め2つのグラフ間のノードもしくはエッジを対応付けした際の利益のようなものを設定しておく.例えば,グラフA, Bについて,Aの1番めのノードとBの2番めのノードを対応付けたら,10点だとか,Aの2番めとBの2番め同士なら5点みたいな感じである.その上で,すべてのノードとエッジを対応付けたあとのすべての適応値の合計が最大となるように対応付けを行ったときの対応をマッチングの結果とするのである.

上記の2つのやり方,すなわち
1. グラフの編集距離を最小とするやり方
2. グラフの各要素の対応付で最も適応関数の値が最大となるものを探すやり方
の2つは,本質的に違う問題のように感じられるだろうが,実はどちらについても二次割当問題(QAP)というNP困難な問題に帰着できることが知られている.したがって,結局の所,この問題を解くことはQAPを如何に高速に解くかという話になる.

手っ取り早く最新の結果を使いたい人へ

本当は,このあとにこの問題がどうして二次割当問題に帰着可能で,過去の研究ではいかにしてこれを解こうとしているのかを説明しようと思ったのだが,今回も例に漏れず面倒になってきたので,すでにあるライブラリだけ紹介して終わりにする.余裕があれば,そのあたりの話も追記したい.調べた感じでは,次の3つのgithubの実装が役立ちそうだ.
– https://github.com/bgauzere/graph-lib
– http://www.f-zhou.com/gm_code.html
– https://github.com/egbertbouman/fgm-cpp
上記の内,2つはC++で実装されており,残りの一つはMatlabだ.特にこれらは割と最近の結果(2017年)のものまで含まれており,その論文によれば500ノードくらいまではぼちぼち動くらしいので,それなりに使い物になるのではないかと思う.

内容についてどう思いますか?
  • いいね (2)
  • だめ (0)

べき乗則の話について

最近,「歴史は「べき乗則」で動く」という本を読んだ.amazonのURI

この本では,世の中の多くの減少における確率分布がべき乗則に従っており,その背景にはフラクタル構造のような再帰的な枠組みが隠されているのだということを終始語っている.これはべき乗則がスケールに対する普遍性を持っており,その背後にフラクタルのようないくら拡大縮小しても構造が変わらないものが存在していると仮定するとうまく説明可能であることに由来している.

本書では,地震の規模と発生頻度,株価の変動,その他生物の絶滅などに至る多くの事柄がべき乗則に従っていることを述べている.本書で取り扱われる議論をそのまま適用すれば,これらの現象を支配する要因にはフラクタル構造が隠れており,それゆえ大地震が発生した理由や,株価が大暴落することには特別な意味があるわけではないことになる.したがって,なぜ人がこれらの現象が起きるのか説明(予測)できないのかといえば,そもそも予測すること自体が不可能な事柄であるからという結論になるのである.

この話は一見,非常に的を得た説明のように感じられる.しかし,本書の最後にある訳者のコメントには,このべき乗則による説明は不十分であるという見解が述べられている.理由としては簡単で,たとえ現象の発生確率がべき乗則ではなく,正規分布に従っている場合でも,平均から大きくハズレた現象が発生することを予想することは困難であるからである.これも確かに正しいように感じられる.

さて,ではどちらの意見のほうが正しいのだろうか.この話を考える際には,おそらく次の2つのことに注意しなければならない.それはある分布からサンプリングすると必ずその分布にそったランダムな値が得られてるということ.またそのことと,予想可能であるかどうかは特に関連がないことである.すなわち,ある現象の発生確率がある分布に従っているとするとき,その情報だけでは,ある瞬間にその現象が起きる可能性は確率的にしか評価できない.したがって,その事実だけではどれほどの大きさで対象となる現象が発生するのか予測することは困難となる.しかし,この予測可能性は,分布のサンプリングだけに着目するのではなく,分布の形に着目すれば,その可否がわかるというのが本書の主題だろう.例えば,正規分布に従っている現象を対象とすれば,正規分布はスケールに対する普遍性を持たないのだから,平均に近い部分に対応する現象と,平均から外れた部分に対応する現象とでは,その発生原因は大きく異なるだろうという予想がたつ.しかし,べき乗則のような分布に従っていた場合には,スケール普遍性により,正規分布に当たるような特定の発生原因(すなわち予兆)を見つけることが困難なので,極端な事象(大地震など)を予想することができないというのが,本書の著者が言いたいことだったのではないだろうか.

しかし,このような説明もあまりに感覚に訴えかけすぎていて,その妥当性がそこまで納得できるものになっていない.結局の所どういうふうに考えればよいのだろうか(誰か教えてください).

内容についてどう思いますか?
  • いいね (1)
  • だめ (0)

配列の中で2番めに大きい要素

最近ブログを書いてないって言われたので,何かを書こうとおもったけれど,何を書くのが全く思いつかないので,今日,聞かれた質問について書こうと思う.

聞かれたのは,与えられた配列の中から2番めに大きい要素を知る方法である.1番大きい要素の取得の仕方はよくある話だが2番めはあまり聞いたことがない気がする.どうやって解くのが最も効率が良いだろうか.ちなみにプログラム演習の課題のようなので,あまり高度なライブラリを使用するのは避けなければいけないという制約がついている.

最も簡単な解法は,一番最大の要素をまず見つけ出し,次にその要素を取り除いた配列から最大の要素を取り出す方法だ.これは愚直だがわかりやすく書きやすい.ただし,最大値を求めるループを2回使う必要があり,いくらか冗長な気がする.次に,思いつくのはソートして2番めの要素を取り出す方法である.スクリプト言語を使っているのであればこれが一番ありがちなのではないだろうか.ただし,一から全て実装する場合は,これは少し面倒だ(ライブラリがまったくないそんな状態は普通ないのだが).

さてここではよりフォーマルな形で定式化することで,この問題を効率よく解く方法を考えたい.まず与えられた配列を\(a_n, \dots, a_1\)という数列であるとしよう.\(a_k, \dots, a_1\)の中で一番大きい要素を返す関数を\(f(k)\)とする.このとき明らかに
$$f(k) = \begin{cases}
-\infty & (k=0) \\
f(k-1) & (f(k-1) \geq a_{k}) \\
a_{k} & (otherwise)
\end{cases}$$
が成り立つ.次に,\(a_k, \dots, a_1\)の中で2番めに大きい要素を返す関数を\(g(k)\)としよう.定義から明らかに\(g(k) \leq f(k)\)である.\(g(k)\)の計算は\(f(k)\)と非常に近いはずである.しかし,2番めであるという制約が,問題を難しくしている.まず,\(g(k)\)と\(g(k-1)\),そして\(a_k\)の関係を考えよう.まず,\(a_k\leq g(k-1)\)のときは明らかに\(g(k)=g(k-1)\)である.次に,\(a_k> g(k-1)\)の場合は,\(g(k)=a_k\)としたいところだが,\(f(k)\)との兼ね合いがあるため,それだけでは条件として不十分だ.実際には,\(f(k-1) < a_k\)のときは,\(g(k)\ne a_k\)としなければ,2番めに大きいという条件を満たせない.逆に言えば,これだけで必要条件を満たすことができる.したがって,まとめれば
$$g(k) = \begin{cases}
-\infty & (k=0) \\
g(k-1) & (g(k-1) \geq a_{k}) \\
a_{k} & (f(k-1) \geq a_{k}) \\
g(k-1) & (otherwise)
\end{cases}$$
である.思ったより簡単だった.

内容についてどう思いますか?
  • いいね (2)
  • だめ (1)

NP困難とNP完全

最近,証明書が切れたみたいでこのページが危険だって表示されるらしい.そろそろ更新しなきゃと思いつつ,何もしてない(めんどくさい).

さて,今回はNP困難とNP完全の話.よく間違えるので書こうと思います(間違えるのはもしかして俺だけ?).

計算量を議論する際に,知ってる人も多いと思うけれどPとNPという問題のクラスが存在する.Pは多項式時間で解ける問題を表していて,NPは答えを知っているときに問題を多項式時間で検証できる問題を意味してる.多項式時間が英語でPolynomial Timeなので,頭文字を取ってきてPだと思うのはいいけれど,NPをNon Polynomial Timeと考えて,多項式時間じゃない問題だと考えるのは間違い.

NPはNon Polynomial Timeじゃなくて,Non-deterministic Polynomial Timeのこと.意味が違うから勘違いしてはいけない.ただこの紛らわしい命名規則がそもそもの混乱の元のような気がする.

さて,多項式時間で解けるとは一体何を指しているのだろうか.これは簡単に言えば,問題のサイズに対して多項式に比例する時間で解けるという意味.例えばソートの問題であれば,性能の悪いソートアルゴリズムは要素数\(n\)に対して,大体\(n^2\)に比例する時間でソートができる.多項式ってのが例えば\(x^3+2x^2+x+-1\)みたいな形をしてるものを指してることを考えれば,できるの悪いソートアルゴリズムの計算時間が多項式時間であることがわかる.対して,各家庭に宅配をする問題を考えてみる.宅配をしなくちゃいけない家のリストとそれぞれの家の間の距離が与えられたときに,最短で回れるルートを考える問題はめちゃくちゃ難しくて,愚直にやると家の数\(x\)に対して,\(x!\)の時間がかかる.これは多項式じゃないので多項式時間で解けるとは言わない.

基本的には,問題のサイズ(例えば要素数が\(n\)とか)に対して,\(n^a+n^b+…+n^z\)(ただしa,b,…,zは定数)で解ければ,多項式時間で解けるってことになる.多項式時間の場合は別に,このaとかbの値がどのくらいの大きさであってもかまわないので,これがめちゃくちゃ大きければ解くのが大変なわけだけど,それでも宅配の問題のような\(n!\)みたいなのが出てこない限り,ときやすいと考える(たぶん).

Pの問題は多項式時間で解けるものを指しているけれど,NPは答えを与えられたときに多項式時間で検証できるものを表している.多くの場合,答えが与えられたらそれがあってるかどうかを検証することは,答えを見つけることより簡単なことが多い.例えばソートの問題だったら,ソートするよりソートが終わってるかどうかを調べるほうが断然簡単なのは明らかだ.だから,ソートの問題はPにもNPにも属していることになる.

他の問題に部分和問題がある.これは,複数個の数が与えられたときに,そこから幾つか選んで足したものが,ある数\(N\)に等しくなるように選ぶことができるかっていう問題.簡単な例として,1,5,7が与えられてたとしよう.ここからいくつか選んで足し合わせたら6にできるかって聞かれたら,即座に1と5を選べばいいってわかると思う.似たように今度は9にできるかって聞かれたら,できないってこともわかる.ただし,これは個数が少ないからわかるのであって,これが膨大な数から選択する問題だったら,簡単じゃない.ただし,膨大な数が渡されていたとしても,答えがわかっていてそれがあってるか確認するのならきっと直ぐにできる.実際,この問題はNPに属していて,さらにNPの問題の中では特に難しい問題だとわかってる.こういうのをNP完全という.

問題が難しいか簡単かという話は,ある問題を多項式時間で別の問題に帰着できるかどうか判断できる.例えば,幾つかの要素の中から最小のものを探したいとしよう.これは要素すべてを調べればわかるから,要素数\(n\)に計算時間が比例する.だけど,要素全体をソートしてから,一番端にある要素を取ってきても同じことができる.ソートするには一番速くて\(n \log(n)\)の時間がかかることがわかってるので,ソートするほうが遠回りになる.でもソートすることでも解けるということで,最小値を求める問題は,ソート問題に帰着できると言える.この例の通り,帰着元の問題は,帰着先の問題より簡単なので,帰着関係で問題の難しさを言うことができる.似たように難しさがわからない問題があったら,計算量が分かってる問題に帰着させる方法を考えればいい.

NP困難の問題は,クラスNPの中で特に難しい問題だと言ったけど,そういう問題よりも難しい問題がある.要するに多項式時間で検証できないような問題が存在するということ.そういったNP完全な問題と同等か,それ以上に難しい問題をNP困難という.NP困難はNPってついてるけど,必ずしもクラスNPに属してるわけじゃない.NP完全より(以上)は難しいという意味.

部分和問題の話の中で,この問題はNP完全だという話をした.この問題はNPには属していそうだけれど,Pには属してなさそうな感じがするかもしれない.しかし残念ながら,部分和問題は多項式時間で解けるかどうかはわかってない.もっというとPとNPが同じかどうかはミレニアム問題という多額の賞金が出る問題となっていて,今現在も全力を上げて研究者が調べてる.なので,P=NPかP≠NPのどっちであるのかわからないけれど,計算機科学を取り扱ってる多くの人の見解は,同じではないだろうということで一致してるのだと思う(おそらく).

誰か頭のいい人証明して~

内容についてどう思いますか?
  • いいね (6)
  • だめ (2)

モンゴメリ乗算

#このページがGoogle検索の上の方に引っかかってるみたいなので,慌てて内容を充実させた.何かおかしな点やわかりづらい点があればコメントを頂ければと思う

#2017/9/30 わかりやすいように編集

\(N > 0\)を法とした合同式を計算する場合,\(N\)による割り算は避けては通れない問題となる.四則演算の中で,割り算が最も遅いことを考えれば,この剰余を求める操作を高速に行うのが簡単ではないことは明らかだ.
しかし,世の中には不思議なことにこれを比較的高速に演算する方法がある.その一つとして今回はモンゴメリ乗算を紹介しよう.

モンゴメリ乗算の形式的説明

まずモンゴメリ乗算を説明するための前準備として,次のモンゴメリ写像\(f\)とモンゴメリリダクション\({\rm MR}\)を導入する.

$$f(x)=xR \mod N \\
{\rm MR}(x) = xR^{-1} \mod N$$

このとき,\(N\)は\(R\)と互いに素であり,かつ\(R > N\)が成り立つとする.
\(R^{-1}\)は,\(R\times R^{-1} = 1 \mod N\)となるような整数である.これは,例えば5が法の場合は,\(R=2\)であれば,\(R^{-1}=3\)となる.また,上式の関係から,モンゴメリ写像とモンゴメリリダクションは逆写像の関係にあるといえる.

次に\(N\)より小さな整数,\(a\),\(b\)を用意し,その積の剰余\(ab \mod N\)を考えたい.ここではモンゴメリリダクションを使って,うまいこと表現できないかを考える.そのためにまず

$$
\begin{align}
f(ab) &= abR \mod N \\
&= aR \times bR \times R^{-1} \mod N \\
&= f(a)f(b)R^{-1} \mod N
\end{align}
$$

であることに注意する.今,明らかに\(f(a) = {\rm MR}(aR^2)\)が成り立つので,

$$
f(ab) = {\rm MR}({\rm MR}(aR^2){\rm MR}(bR^2))
$$

両辺に対して,モンゴメリリダクションを取れば

$$
ab \mod R = {\rm MR}({\rm MR}({\rm MR}(aR^2){\rm MR}(bR^2)))
$$

さて,ここまでは単純に式の形を変えたにすぎない.当初の目的は如何に\(\mod N\)を高速に計算するかだった.それをモンゴメリリダクションを用いた形に直したところで,その中に剰余演算があるのでは意味がないように思われる.ところがこのモンゴメリリダクションを高速に計算するアルゴリズムがある.すなわち

$${\rm MR} (T) = \begin{cases} g(T) & (g(T) < N) \\ g(T)-N & (otherwise) \end{cases}$$

$$g(T)=\frac{T+(TN’ \mod R)N}{R}$$

ただし$$NN’\equiv -1 \mod R$$が成り立つ.

モンゴメリリダクションの計算方法を見ると,今までネックだった\(\mod N\)がいつの間にかなくなり,代わりに\(\mod R\)だけになっている.他にも\(R\)による除算が増えているが,結局どちらも\(R\)を恣意的に選べることから,\(R\)を2のべき乗にすれば,アンド演算とシフト演算のみでそれぞれ計算できる.残りは除算以外の演算なので,それほど重くない.よって,モンゴメリリダクションで\(ab\)を高速に計算する方法がわかったことになる.

さて,ここまではモンゴメリ乗算についての形式的な説明である.モンゴメリ写像から始まり,モンゴメリリダクションについて紹介した.また,その高速な計算方法を天下り的に与えたので,上の話は多くの人にとって納得できるものではないだろう.そこでここからは,上記の内容について1つずつ整理して説明していく.

モンゴメリ乗算の意味

まず,そもそもモンゴメリ乗算とは,ある数の積をモンゴメリ写像したものを指す.例えば\(N\)未満の数,\(a,b\)を用意して,
$$
\begin{align}
f(ab) &= abR \mod N \\
&= (aR)(bR)R^{-1} \mod N \\
&= f(a)f(b)R^{-1} \mod N
\end{align}
$$
はモンゴメリ乗算である.この定義から明らかだが,別にモンゴメリ写像された2つの変数の積は,モンゴメリ乗算ではない.その後, \(R^{-1}\)を掛けてから\(\mod N\)を取って,はじめてモンゴメリ乗算となる.そしてその「\(R^{-1}\)を掛けてから\(\mod N\)」を取る演算がちょうどモンゴメリリダクション\(\rm MR\)の定義になっている.

また別の言い方をすれば,モンゴメリ表現(モンゴメリ写像されたもののこと)された2つの変数の積をモンゴメリ表現(モンゴメリ乗算の結果)にするには,積に対してモンゴメリリダクションを取ればいい.よって,例えば\(abc \mod N\)を求めたいとすれば,まず\(a,b,c\)それぞれのモンゴメリ表現\(A,B,C\)を用意し,次に
$$
D = {\rm MR}(AB)
$$
を求める.\(D\)はモンゴメリ表現なので,\(D\)と\(C\)のモンゴメリ乗算を考えることが出来て,
$$
E = {\rm MR}(DC)
$$
となる.モンゴメリ写像の逆変換はモンゴメリリダクションなので,結局\(abc \mod N\)は
$$
{\rm MR}(E)={\rm MR}(f(abc))=abc
$$
で与えられる.

このように,与えられた変数をすべてモンゴメリ表現にしてから,掛けてはモンゴメリリダクションを行っていけば,最終的にすべての変数の積のモンゴメリ表現を得ることが出来る.あとはそれをモンゴメリリダクションすることで知りたい答えが得られるというのがモンゴメリ乗算の全容である.

そうすると,モンゴメリ乗算は,モンゴメリリダクションが肝であることがわかる.もしこのモンゴメリリダクションが高速に行えればモンゴメリ表現を素早く求められることになるので,結果的に剰余演算が高速できることになる.そして,そのモンゴメリリダクションを高速に計算する方法が先程示した計算法だったわけだ.ではそこにある仕組みは一体何なのだろうか.それを探るためにもう一度モンゴメリリダクションの計算式を以下に示す.

$${\rm MR} (T) = \begin{cases} g(T) & (g(T) < N) \\ g(T)-N & (otherwise) \end{cases}$$

$$g(T)=\frac{T+(TN’ \mod R)N}{R}$$

ただし$$NN’\equiv -1 \mod R$$が成り立つ.

まずは \(g\)の式に注目しよう.こいつの分子は与えられた引数\(T\)に\(N\)の定数倍を足し合わせている.したがって,この分子の式自体は\(\mod N\)において,値を保存することは確かだろう.次に,面白い事実として\((TN’ \mod R)N\)は\(T\)の\(R\)以下のビットの桁をすべて繰り上げて0にするような数となっている.これは例えば\(R\)が128であれば,\(T+(TN’ \mod R)N\)は下位7ビットまでがすべて0となるということだ.試しにこのことを確認するためのC++のコードを次に示す(一部wikipediaのパクリ).

上記のコードのT,R,Nの値を色々と変更してみてほしい.どのパターンでも述べたとおりになっているだろう.要するに,これは\(T\)の下位桁のビットがすべて0になるような数を足し合わせることで,下位桁の値を上位桁に寄せ集めている事になっている.この寄せ集める操作は,あくまで\(\mod N\)においては保存するので,結果的にビットの広がりを小さくする効果のみを持つ.最後に寄せ集めた数を\(R\)で割って,右にシフトさせることで,\(N\)までの範囲に収めようとしているのである.ただし,\(T\)の値によっては,\(g(T)\)が\(N\)を少しだけはみ出すことがあり得るので,それを最後の条件式で補正している.

以上がモンゴメリリダクションの仕組みである.

内容についてどう思いますか?
  • いいね (15)
  • だめ (0)

ニュートン・ラフソン法と極値

*赤げふさんからコメントを頂いたので修正 2017/3/26

ちょっと混乱したのでそのまとめ。

世の中には非線形方程式に解を与えなければならないことがある。この手の問題は解析的に解を求められるのが最善であることは確かだが、実際にはそうかんたんではない場合も多い。そのような場合、コンピュータで数値的に解を求められることがある。その方法の一つにニュートンラフソン法と呼ばれるものがある。

ニュートンラフソン法の場合、対象とする問題は$$f:\mathbb{R}^n \to \mathbb{R}^n$$のような関数、すなわち

$$ f(\bf{x}) = \bf{0} $$

のような非線形連立方程式の解を求めたい場合に使用する。注意してほしいのは、入力と出力の次元が等しい非線形方程式でないといけないことである。そうでないと使用できない(と思われる)。さて、この問題を解く場合、ニュートンラフソン法ではある入力\(\bf{x_0}\)における一次近似式を求め、これが\( 0\)となる場合の\(\bf x\)を求めるというのを繰り返すことで解を求める。これに対する直感的な説明は、他のページに譲るとして、この方程式を解いてみよう。$$f(\bf{x})$$の一次近似式は

$$ f(\bf{x}_0)+ Jf(\bf{x_0}) (\bf{x-x_0})$$

ここで\(Jf(\bf{x_0})\)はヤコビ行列である。この一次近似式が\(0 \)となる点を求めればいいので

$$\bf x_1 = x_0- Jf^{-1}(\bf{x_0})f(\bf{x}_0)$$

となる。ここで\(\bf{x_1}\)が次の点となる。これを、\(\bf x_1\)を次の点として同様の操作を逐次繰り返せば求めたい解となる。ここでヤコビ行列の逆行列を求める必要があるので、関数\(f\)は入力と出力の次元が等しい必要がある。

さて、本来の場合、出力が\(0 \)となる入力を求めるのが目的だが、場合によって関数の極値を求めるのにも使用される。その場合入力の次元が多次元で、出力が1次元の場合を扱う。というのも、極値となる点(候補の点)は、\(\nabla f(\bf x)\)が\( 0\)となる点を求めることに等しいため、ナブラを使用する都合上、出力が入力の次元と等しくなる。すると、そのままニュートンラフソン法が使用できるからである。

そして多くの文献ではこれをテイラー展開によって説明する。すなわち関数\(f( \bf x)\)のある点\(\bf x_0\)を中心とした二次までのテイラー展開

$$f(\bf x) \approx f(\bf x_o) + {}^{t}(\bf x – x_0) \nabla f(\bf x_0) + \frac{1}{2} {}^{t}(\bf{x – x_0}) \nabla \nabla f(\bf x_0) (x-x_0) $$

を考え、この両辺の勾配を取る。

$$\nabla f(\bf x) \approx \nabla f(\bf x_o) +  \nabla \nabla f(\bf x_0) (x – x_0)$$

ここで\( \nabla \nabla f(\bf x_0)\)はヘッセ行列*1である。さて、テイラー展開はn次まででやめると、n次のテイラー多項式と呼ばれ、それ自体がn次多項式近似となること、また多項式展開は基本的に一意的であること*2を考えれば、先程の勾配を取ったものは\(\nabla f(\bf x)\)の一次近似式であることがわかる。したがって、先程と同様の議論で

$$\bf x_1 = x_0- \nabla \nabla f^{-1}(\bf{x_0})f(\bf{x}_0)$$

が得られる。基本的に最初の例と原理は同じなので、これもニュートンラフソン法と呼ばれる。さて、これが一体何の役に立つのかということがわからない読者も多いかもしれない。私は役に立たないとやる気が起きないという物臭な性格をしているので、そういうあたりが気になるところだが、例えば複数のパラメータを調整することで、与えられた点を近似するような関数を求めたい場合に使用する。具体的にはそれはカメラの内部パラメータを求める場合であったり、またあるいは機械学習でパラメータ推定する場合等がある。注意しなければならないのは、この調整しようとする関数の出力は別に多次元でもかまわないという点である。なぜなら大抵の場合これらの近似における指標はスカラーの誤差関数で与えられるので、多変数を取るスカラー値を返す関数の極値を求める問題に帰着できるからである。

*このヘッセ行列は、当然だが\(\nabla f(\bf x)\)のヤコビ行列になっている。
*この事実は強力だ。というのも例えば\(\frac{1}{1-r}\)のテイラー展開を知りたいとしよう。さて、勉強熱心な読者はこれは等比級数で表現できることを知っているだろう。すなわち
$$\frac{1}{1-r} = 1 – r + r^2 – r^3 + \cdots $$
が成り立つ。そして多項式展開の一意性よりこれがテイラー展開の結果でもあることがわかる。更にこれを両辺積分すれば
$$\log(1-r) = r – \frac{1}{2}r^2 + \frac{1}{3}r^3 + \cdots $$
であることまでわかり、また同様にこれもテイラー展開した結果でもある。ちなみに積分定数はr=0のとき0であるので省略した。

内容についてどう思いますか?
  • いいね (4)
  • だめ (2)

カルマンフィルタとEMアルゴリズム

この辺りの話についての記事が恐ろしくWEB上では枯渇しているようだ。中にはMATLABによる実装なども見つけることができるが、私がMATLABが好きではないという至極個人的な理由で、そのようなものは受け入れられない(そもそもあんなに言語として欠陥のある物を、高い金を出して購入するという風習を理解できない)。

今回、諸事情でカルマンフィルタと、それのオフライン実装、加えてEMアルゴリズムによるパラメータの推定までを行う必要があった。そこで、そのコードに多少の説明を加えて公開しようと思ったわけである。ただし、私は別に専門家でもなんでもないので、場合によっては内容に不備や間違いが含まれている可能性がある。したがって、ここに記載された内容により何か不利益を被るようなことになったとしても、私は責任を取ることは出来ない。よって、読者は内容を鵜呑みにせず自分の目で正しいかどうかをよく確かめてほしい。また、なにか気づいたことがあれば、些細な事でも連絡をいただけるとありがたい。

まずは定義式を確認しよう。

$${\bf z}_n = {\bf A}{\bf z}_{n-1} + {\bf w}_n$$

$${\bf x}_n = {\bf C}{\bf z}_{n} + {\bf v}_n$$

$${\bf z}_1 = {\bf \mu}_{0} + {\bf u}$$

ただしここで

$$p({\bf w}) = {\cal N}({\bf w|0,\Gamma})$$

$$p({\bf v}) = {\cal N}({\bf v|0,\Sigma})$$

$$p({\bf u}) = {\cal N}({\bf u|0},{\bf P}_0)$$

さてカルマンフィルタの目的は、観測値$${\bf x}_n$$から内部状態$${\bf z}_n$$を導くことである。詳しい説明は他のページに譲るとして,この問題はガウス性ノイズが邪魔となり簡単には求めることは出来ない。しかし、うまいことやることでこれは計算可能である.

ここまでの話はどこにでも乗っていることだ。重要なのはここからで、今までの枠組みでは$$\bf \Gamma$$や$$\bf \Sigma$$などは事前に知っている必要がある。しかし、世の中そんな都合のいい話ばかりではない。私たちは問題に取り組む際、パラメータを都合よく知っているといつでも仮定できるわけではない。実際問題としては、観測値$${\bf x}_n$$から内部パラメータを推定できることが最善だろう。そのための解決策としてここではEMアルゴリズムを使用する。これについても詳細は省略するが(知りたければコメントを書いてほしい)、要約すれば初期値を与えたあとに、そこから反復計算を行い、収束した頃にはだいたいちょうどよいパラメータ値になるというものである(超雑)。詳細は置いといて、実際のコードを下に示そう。

これがクラスファイルである。と言ってもstructかつすべての変数がpulicなので到底まともなクラスとはいえないが。使い方は次の通り

ここではノイズの乗ったsin関数について、そのノイズ除去を行っている。実際の出力結果をgnuplotで表示したものがこれ

test

色が気持ち悪いことを除けば、特に問題はなさそうだ。

内容についてどう思いますか?
  • いいね (3)
  • だめ (0)