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

機械学習について思うこと

ブログだけに書くのもいいかなと思っていたのだけれど,やはり文章を書く以上は人の目に触れたほうが良いと思ったので,ブログとQiitaで同時に投稿することにした.どちらもmarkdownでかけるので共有するのも楽だしね.

以下の文章は極めて長いが,お付き合いいただければ幸いである.

機械学習について思うこと

これは私が常日頃思っていることで,もしかしたら私がここで言及しなくてもすでに多くの場所で語られていることもしれない.また,私自身は別にこの分野の専門家ではないので間違ったことを言っているかもしれないが,その場合は意見をいただければと思う.

さて,機械学習,特に深層学習は目覚ましい進歩を遂げたが,この技術はいつも絶対的な保証ができないという問題を抱えている.そしてそれを理由に工学的には利用が難しいと考えられている.この文章ではなぜそのように考えられているのか,そしてそれを解決することが原理的に可能なのかを議論したい.

医療における信頼度

DNNは医用画像の分野でもよく使われている.これは例えば,従来,内視鏡画像やレントゲンなどをもとにガンなどの様々な疾患を発見することが主な目的である.このような仕事は,疾患を見逃せば重篤な問題を引き起こす可能性がある上に,高い精度を求めるためには熟練した医師が必要なため,今までは人間の医師が行ってきた.

しかし,近年注目を浴びているDNNを用いることで,こういった問題を解決する兆しが見えてきている.DNNはタスクによっては人間を超える性能を叩き出しており,特に今まで自動化が難しいと考えられていた,上記のような診断を医師無しでできることは,医師に対する負担が減る上に,経済的にも良いことであると考えられる.

一方で,このような機械学習に基づく診断は,その診断結果が高い識別率を持っていたとしても,その信頼性に疑問が呈されることがある.これは,診断結果が間違っていた場合の危険性を考えると,例え識別率が人間を上回っていたとしても,どうしても不安が残るためだ.そして,可能であれば,理論的に識別結果に対する信頼性を定義することで,診断結果に説得力をもたせられないかという考えに至るようだ.実際に,この話は他の研究室にいる知り合いが,教授から提示された課題の一つである.

さて,ではこのような理論的な裏付けはつけられるだろうか.最初に,これに対する私の考えを述べさせてもらうと,おそらく不可能だと思う.以下ではその理由について説明したい.

そもそも機械学習とはなにか

すでに機械学習は広域な分野を有しており,その目的を一括りに一言で述べることはできない.しかし,医用画像に対する診断のような,教師あり学習に限れば,それは簡単に言うことができる.すなわち,未知の入力が与えられた際に,その予測を得ることである.

多項式回帰の例

これの最も簡単な例は最小二乗法による回帰問題に見えることができる.次の図を見てほしい,
回帰直線.png
今図中で,バツじるしがついているところが与えられたデータ点である.ここで,知りたいことは,データ点が与えられていない\(x\)について,対応する\(y\)を知りたいということである.そしてこのような問題に対して,最小二乗法では,適当な近似曲線(ここでは直線)を引き,これをもとに\(y\)の値を予測する.これは非常に簡単なタスクだが,例えば先ほど説明した医用画像の例も,横軸を画像,縦軸をガンを発症している確率とすれば,やはり近似曲線を引く問題と同じものとみなせる.

さて,上の図はあくまで直線で近似しているため,これでは教師データを真似ることすらできていない.多くの文脈では,ここで直線のような一次式ではなく,より高い次数の次元で近似する.すると例えば,
回帰直線2.png
のような,赤や緑の曲線を得ることもできるだろう.ここで,赤い曲線は教師データをいい感じに通っている感じがするが,緑の曲線はやりすぎ感がある.特にこのようなやりすぎな場合を普通オーバーフィッティングと言ったりする.対して,先程の直線回帰のようなあまり近似できていない場合をアンダーフィッティングという.

機械学習の分野では,このような問題をバイアス・バリアンス分解によって理解することが多い.多項式回帰では,このような高次の多項式を用いることによる振動問題をルンゲ現象と呼び,次数が大きい多項式の微分値が大きくなりやすいことをその理由とする.

さて,そもそもの目的は,この教師データを生み出した”真なる曲線”を予測することだった.今までの例で挙げたこれらの曲線は,この問題を本質的に解決しているのだろうか.答えはわからないである.オーバーフィッティングにしろ,アンダーフィッティングにしろ,それがだめかどうかは,実際に教師データ以外の実データで検証を行う以外にほかはない.

確率的解釈

多項式回帰は簡単だが,そこから多くのものを得られる.例えば多項式は入力データから作られる(例えば\(x\)や\(x^2\)など)が,別に多項式である必要はない.任意の場合に目を向けると実はカーネル法につながっている.加えて,この回帰問題は比較的簡単に確率的な解釈を導入することができる.例えば正規分布を使ってみよう.この場合,\(x\)と\(y\)の関係は,\(\mathcal{N}(y|\mu, \sigma^2)\)とかける.ここで\(\mu\)と\(\sigma\)は\(x\)で与えられることにすれば,確率的な解釈が導入されたことになる.そして,これを応用することでベイズ的な取扱を行うことすら可能となる.しかし,そのような高度な解釈を導入したところで,上記の真の曲線(分布)を予測できているかどうかはわからないのである.

なぜうまくフィッティングしていることが言えないのか

どのようなモデルが与えられたところで,それが真に求めるべき曲線を近似しているかどうかは理論的に知る方法はない.これはなぜなら,データ点がないところについては,我々はなんの知識も用いることができないからだ.

そもそも数学的には,このような回帰曲線を引く問題は,\(x\)が与えられた際に,それらしい出力||(y\)を出す関数\(y=f(x)\)を求めることに等しい.そのような関数は,教師データがない部分に関しては,何を返しても良い.例えば,教師データがある部分については,デルタ関数的に振る舞い,それ以外の部分については恒等的に0であっても何も問題がない.

直感的には,そのような関数ではだめなことはわかるが,では何がどのくらいだめなのかを知るすべはない.いい感じで回帰曲線を引くという意味では,別に多項式回帰を使おうが,カーネル法を用いようが,SVMを利用しようが,DNNに投げようが,どれでもそれらしい曲線が得られるはずである.この手の議論でよく引き合いに出されるのが,あるモデルが任意の関数を近似可能かどうかという問題である.近似できない関数が存在すれば,もしかしたら与えられた教師データを通るような関数を作り出せない可能性がある.したがって,パラメータを十分に増やすことで,任意の関数を表現できることを示すことには一定の価値がある.しかし,任意の関数を近似できることは,今言ったデータがないところを予測することに関しては何も言及していない.

例えば,最も簡単な多項式回帰は数学的には非常に強力なモデルである.というのも,多項式回帰はテイラー展開,もっと言えば複素数空間ではローラン展開であるとみなせる.しかし,ローラン展開はフーリエ級数展開と同等であるので,ほぼ任意の関数を近似できるはずである.しかし,それによって教師データを近似した結果が,良いものであるかどうかについては何も教えてくれない.

DNNがよくわかっていない最大の理由は実はそこにあるのではないかと私は思っている.高い予測性能を叩き出せるということは,言い換えれば,データ点がないところの予測が優れていることにほかならない.しかし,その部分について数学的に妥当性を言うことはできない.機械学習において,あるモデルが優れているかどうかは,究極的には,人間にとって,予測結果がどれほど自然であるかというその一点に集約されるのだ.

もし解決策があるとしたら

私は上記の問題に対して全く解決策がないとは考えていない.ぱっと思いつくだけで2つほど解決策があると考えている.

1.真の正解に制限をつける

フーリエ級数展開を考えよう.フーリエ級数展開は周波数成分に分解するための手法だが,もう一つの捉え方は曲線を近似するための方法でもある.例えばローパスフィルタは低周波成分だけを残すことで,データ点を通る近似曲線を求めているという解釈ができる.実はこのような考えは更に進めると,ウィナーフィルタやカルマンフィルタ,パーティクルフィルタのような話につながっている.

さて,フーリエ級数展開にはサンプリング定理が柱として存在する.これは,現信号の周波数がある特定の範囲までにしか分布していないのであれば,現信号を復元可能であることを保証する.フーリエ級数展開に利用価値があるのは,サンプリング定理によりその理論的背景が存在し,加えて,私達の周りにある多くの信号は,最大周波数がわかっているためである.

このようなサンプリング定理のようなことが,機械学習で用いられるデータに対しても同様に言えるのであれば,理論的な裏付けができるだろう.というのも,予測したいデータがどの範囲まで分布しているのかがわかれば,予測が最悪でどこまで外れているかを定性的に評価できるからだ.しかし,これが解けてしまうと,そもそも過学習のような問題は起きなくなるのではないかという気もするため,もしかしたら解けないのかもしれない.

2.脳を理解する

機械学習で対応しなければならないタスクの多くは,我々人間がその正解を決めていることが多い.例えば,音声認識を考えよう.この場合,入力された音声に対する正解は我々人間が聞いたときの感覚でラベル付をしている.音声合成も同様である.

これは人間にとって自然とは何かがわかれば,正解も自ずと分かることを示唆している.したがって,未だよくわかっていない脳科学が進展すれば,真の正解を知ることができるかもしれない.

しかし,これは極めて茨な道である.この理由の一つは生きている脳細胞をリアルタイムで計測する方法には,どれも問題が多いということである.例えば,脳波を取る方法にしろ,fMRIなどを用いる方法にしろ,大抵は時間的,もしくは空間的分解能が低すぎる.脳の活動が大量の,微小なシナプスやニューロンによって構成されていることを考えれば,現状の計測方法では十分とは言えない.マウスの脳を計測する方法に限れば,侵襲的な計測方法も用いることができる.例えば,電極アレイを脳に挿して電位を図るなどである.しかし,これも拒絶反応が起きるなどの様々な問題を抱えている.現状他に期待されている方法には,マウスの頭蓋骨をガラスに置き換え,光学的な方法で観測する方法も検討されているが,これも完璧とは言えない.

ただ,最近は脳細胞を人工的に作り出すことが可能になりつつあるようなので,そういった方法も使えるかもしれない.

関連する話

機械学習の問題は実は統計で出てくる問題と本質的に変わらない.それは機械学習がそもそも統計学の影響を強く受けているからだろう.統計学で出てくる問題の多くは母集団がわからないことに端を発する.

例えば,偏差値を考えよう.学生にテストを受けさせて採点を行う.採点結果は標本となり,その標本から母集団を予測する.母集団が予測されたら,そこから偏差値を計算する.しかし,実際のところ試験を受けていない学生については,統計は何も教えてくれない.もし,模試によって偏差値を算出したとしても,その模試を受けなかった学生の中に極めて成績の悪い,もしくは良い学生が大量に存在したら,実際の母集団の平均や分散は予測されるものより大きく異なることになるだろう.

では,統計学ではそのような問題をどのように解決するのだろうか.それは要は数の暴力である.統計で現れるt検定やχ二乗検定にしろ,標本による母集団の推定にしろ,十分なサンプル数があれば,確率的にその結果が怪しい可能性が十分に低くなる.逆に言えば,そのような手段でしか統計学では保証するすべがない.それにたとえ十分なサンプル数が存在しても,確率的には万が一を避けられない.この問題が機械学習における限界と同等のものであるように感じられる.

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

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