C++」タグアーカイブ

似ている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)

関数型言語に関して思うこと

ここ最近というわけではないが,関数型言語に関する話題をよく目にするようになったような気がする.何を持ってして関数型言語であるかというのは非常に悩ましい問題ではあるが,おそらく最も古くに生み出された関数型言語がLispであろうことを考えれば,関数型言語の歴史は非常に古い.それに対して,世の中のコンピュータの基盤であるOSなどの低レイヤを支えるc言語がLispより20年近くあとに生み出されたことを考えれば,c言語やその影響を受けた後続の言語(JavaやLLなど)がいかに新しいかがわかるはずだ.

私はもともとc言語からプログラミングを入門し,次にC++を学び,Java,Pythonと進めてきた.特にC++に触れてきた時間は非常に長く,プログラミングの思考回路がC++をもとに組み上がっていると言っても過言ではない.したがって,初めて関数型言語に触れたときは強い隔たりを感じた.特に違和感を感じたのは,関数型言語の記述がアセンブリから完全に乖離している点だった.CやC++などの手続き型言語が,何をしたいのかを書くのに対して,関数型言語が何が成り立つのかを書くというのはよく言われる話だが,コンピュータの根幹をなすアセンブリ言語が前者の立場にあることを考えれば,関数型言語がそのままの形では全く相容れないことは明らかである.

ではなぜそのような関数型言語が今になって注目されているのだろうか.その一つの答えは,関数型言語を支える非常に強力な仕組みである型システムが,プログラムの記述を簡単にし,更にはバグさえも防ぎ得るということにあるではないだろうか.そもそもc言語やそこから派生した言語は,人間がミスをすると簡単にバグが挿入される可能性がある.これは例えれば鋭いナイフのようなもので,その切れ味をうまく活かせれば大抵のことは効率よく行うことができるが,使い方を誤れば自身を傷つけてしまう諸刃の剣となる.これはc言語の前身であるB言語の更に前身であるBCPLの哲学である

BCPLの哲学は、最良を知ると皆が考える暴君がいるのではなく、何が許されて何が許されないのかという決まりが定められているのでもない。どちらかといえばBCPLは、たとえ明らかに馬鹿げた事態に直面したときでさえ、不平を言わずに自らの能力を最大限に活かしてサービスを提供せんとする召使いとして振舞う。彼が何をしていて、細かい制限に制約されないということを、プログラマーは常に理解しているものと考える。
「BCPL, the language and its compiler」より

にあらわれている.これはプログラマが過ちを侵さないという絶対の信頼のもとで成り立っている.しかし,実際にはそのようなことは絶対にありえない.むしろ世の中にこれほど人為的ミスとしてのバグがはびこっているのは,人間がそこまで賢くないことを示唆しているのではなかろうか.

さて,人の書くプログラムには絶対にバグが入り込むと考えたとき,これを防ぐためにはどのようにしたら良いだろうか.関数型言語における,この問題に対する答えはバグが入り込む可能性があるような書き方は,そもそも許さないというものである.関数型言語では,関数の副作用を極力許さないことと,型に基づいた統一的な記述により制約をかけている.このような制約は,動的型付けであり副作用がてんこ盛りのLLから見れば,非常に強いものだが,バグを防ぐ上では非常によく機能することがわかっている.

ある言語が,手続き型言語と関数型言語のどちらであるかは程度の問題であって,完全な線引があるわけではない.pythonやperlなどであっても記述の仕方によっては関数型言語のように書くことができるだろうし,LispやOCamlでも手続き型言語のような記述ができるはずである.しかし,言語のサポートによって関数型言語として代表されるような記述がどこまで可能であるかが変わってくる.特に大きな線引は高階関数が扱えるかどうかではないかと思う.関数型言語が特に流行りだした頃から,多くの言語でラムダ式(無名関数)をサポートするようになったが,部分適用が言語レベルでサポートされているのは殆ど見たことがない.しかし,関数型言語にとって高階関数が無くてはならない存在であることを考えれば,これをサポートしていない多くの言語は一般には関数型言語とは呼べないのではないか.

とにかく,高階関数はもちろんとして,完全に純粋すなわち副作用を持たない言語としてHaskellが有名である.Haskellのすごいところは,その徹底ぶりにある.多くの関数型言語は,副作用を全く持たないことを保証できていない.これはおそらく,利便性の観点から完全に副作用を除外すると,非常に書きづらくなる恐れがあるからではないだろうか.画面に文字を表示するということですら,簡約のやり方によっては表示の順番が変わってしまう恐れがあるため,そのままでは副作用は避けられない.しかし,Haskellではモナドを用いることで,理論的にも妥当なやり方で副作用を追いやることに成功している.さらに,Haskellは遅延評価を用いることで便宜的とはいえ,加算無限の概念を記述することも可能である.Haskellのすごいところは,そういった,一見達成困難なことを理論的に美しい方法で,コンピュータ上の実装のレベルまで落とし込んでいる点にある.これにより,Haskellは数学と非常に相性が良い.特に関数型言語の利点である,性質を記述することでプログラムが可能であるという点がうまく引き出されており,これは数学の式の記述に非常に近いものになっている.

ここまで思ったことをつらつら書いてみたが,あまりに何も考えてなさすぎてとりとめのないものになってしまった.とにかく,最後に私が言えることは,関数型言語を全く知らない人はどのようなものかを一度勉強してみるべきだということだ.そこには,手続き型言語だけでは感じることができない,理論によって整理されたきれいな世界が広がっているのだから.

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

やりたいこと

やりたいと思っていることがいくつかある.箇条書きにすると

  • FPGAのCPU実装をもっと頑張る(具体的にはUnix系OSを動かす)
  • Rustをやる(できれば研究とかバイトのコードも徐々に移行したい)
  • 日本語力を上げる
  • もっといろんなことを勉強する(暗号,音声,機械学習,etc)

のような感じ.別に最後のやつは普段のタスクの中に自然と組み込まれている感じなのであまり意識する必要はないだろうけど,上から2つは個人的に時間を作ってやらないとできるとは思えない.3つ目は文章が下手くそなのでうまくなりたいなということ(願望).

FPGAのやつはできればCPUの実装のコード(verilog)をきれいにしたというのもある.研修Aという研究室の自由課題でやったものだから,締切に合わせて結構その場しのぎでやってしまった感があって,それをどうにかしないと次に進む上で障害になりそう.かと言ってあれってそんな簡単に直せるとは思えない… 抜本的に書き直す羽目になりそう.そもそもverilogってアセンブリに近いものがあって全然きれいに書けてる自信がわかないんだけどどうすりゃいいんだろう.高位合成とかと違って結構低レイヤだから記述量が多くなりがちなのが影響して,とてもじゃないけど見やすくならない.あれも結局経験値で決まってくるやつなんだろうな.

Rustについて言えば,とにかく実践で使っていくしかないだろうなって思ってる.入門書とか初心者向けの記事とかをただ追って読んでるだけだと全然頭に入ってこないたちなので,実際に使って鍛えていくしかない.自分もそんなに詳しいわけじゃないんであれだけど,c言語とかc++みたいなシステム周りとか結構低いレイヤのプログラミング言語って,D言語とか他にもあるんだろうけど結局のところc言語 or C++がデファクトスタンダードって感じがあるように思う(マイコンとかの開発環境でC++/C以外の言語を使う場面ってなくはないけど少ない気がする).LLだとか関数型言語とかだと選択肢がありすぎるくらいにはあるのに,低レイヤ向けの,アセンブリより抽象度が一つ上なメジャーな高級言語ってのが,もう何十年も前に発明されたCとかC++のまま変わってないってのは,あまりに時代遅れな感じがする.と言いながらもう10年近くCだとかC++を書いてきた人間がこんなことを言うのは説得力にかけるのかもしれないけど,流石に限界というかこのあたりにもイノベーションがあってもいいよねってのは薄々思っていたところ.確かにC++は最近になってC++14とかC++17とか新しくなってきてるけど,もともと難解だったC++に更に付け足しつけたしで色々くっつけたもんだから複雑怪奇になって来てて,「C++を理解してます」ってのがジョークになるという笑うにも笑えない状態になってる.

だから流石にそろそろ新しくそういうシステム周りもかける言語を習得したいなと思ってたわけだけど,ちょうどMozillaがRustを開発してるという話を知って,これは良さそうと感じた.最近はライブラリも結構充実してきてて特に問題もなさそう.他に最近話題のコンパイルするような静的型付け言語としてはGoとかもあるけど,あれはポインタを叩くような処理をするというよりは,コンパイルできるLLみたいなイメージが有る.Google的には速度が出るLLが作りたかったってことなんかなという印象(素人目線).個人的にはやっぱりLLとシステム記述向けの言語って住み分けができてると思ってて,前者が簡単にかける代わりに詳細な記述をしたり速度を出したりすることが難しいもので,後者がその反対という感じ.よくGoとRustを対比させる記事を見るけど,それはこの2つがちょうどLLとシステム記述向け言語の間くらいに位置する言語だからということなんだろうな.ただGoはどちらかと言えばLL寄りで,Rustがシステム記述向けよりなんだろうと思うし,だからこの2つは用途が衝突するところもあるけど,全くかぶらないところもあるんだと思ってる.

確かにLLって半端じゃないくらいに遅いから早くしたいと思う人が多いのは当然だよね.pythonくらいしか経験がないけど,それでもC++なら一瞬で終わる処理が半端じゃないくらいに時間かかったりするのを見て,何とも言えない感じなることはある.プログラミング界隈で有名な話?に80:20則ってのがあって,プログラム全体の20%が実行時間の80%を占めてるって経験則なんだけれど,LLの哲学的には,重い20%のコードをCとかにまかせて,速度がでなくても大丈夫な80%ところをLLでどうにかするって話だったんだと思う.確かにそれは合理的に見えるし,全部を無駄に書くのが大変な言語で修行のごとく記述していくのは意味がないってのは大いに賛同するところなんだけど,現実はやっぱりなかなか厳しいものがあって,そう簡単じゃない.実際に使ってみると,LLがおそくてその80%が馬鹿みたいに大きくなったりして,僕の書き方が悪いってのはあるんだろうけど,なんだかなーって結果になることも往々にしてある.そんなこんなあって,みんな早いLLを求めてたんだろうから,Goみたいな言語が出て来たんだろうなという予想.本当のところどういうことなのかはGoを勉強しないとわからないけど,今のところGoが重要そうなWeb界隈の方には全く関わりがないので,実用的な観点でやることはなさそう.

とまあ前置きが長くなったけどRustを勉強したいなと思ってる.関数型言語の流れを受けて色々と拡張があるみたいだし,加えてC++とかと比べて安全だと言われればやらない理由も特にない.速度も結構早いみたいだしね.

ということで今年度も頑張っていきたいと思います.

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

モンゴメリ乗算

#このページが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)

カルマンフィルタと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)