DBScan
DBScanという技をラーニング
論文
http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.71.1980
簡単な説明
http://ibisforest.org/index.php?DBSCAN
x-means
http://d.hatena.ne.jp/m-a-o/20100610#p1
などと同じく、クラスタ数を指定する必要がないクラスタリング手法
具体的には、点集合が与えられた時、点を二種類に分類する。データ間の距離が与えられているとして、距離eps以内にN個以上の点があるような点をcore pointと呼び、N個未満しかないような点をborder pointと呼ぶ。ただし、epsとNはクラスタリング時に指定するパラメータ。つまり大雑把には、ある程度高密度にデータが集積してる点と、そうでない例外的な点を分ける。互いに距離eps以内にあるcore point同士を同じクラスターに入れるようにすると、core pointのクラスタリングができる。
一方、border pointについては、次の3つの状況
がありえる
(1)距離eps以内にあるcore pointがない
(2)距離eps以内にあるcore pointは全て同一のクラスターに属する
(3)距離eps以内にあるcore pointが複数のクラスターに属する
(1)の場合、そのborder pointはNoiseとして扱われる。(2)の場合は、近くのcore pointが属するクラスターに入れてやればよい。(3)の場合は、どうすればよいか。論文では、このケースを考慮し忘れている気がする。(3)のようなことが起きると、一意にクラスタリングできない。
(3)のようなことが起きる例として、例えば
S0 = [(0.0 , 0.0)] S1 = [(0.99+0.1*n , 0) for n in xrange(10)] S2 = [(-0.99-0.1*n , 0) for n in xrange(10)] S=S0+S1+S2
のような点集合を考えて、(eps,N)=(1.0 , 3)とすればよい。距離は、普通にユークリッド距離で測る。もし、S1とS2の点が(論文Definition5の意味で)同じクラスターに入るとすると、S1の点とS2の点は互いにdensity-connectiveにならない。一方、これらがそれぞれ別のクラスターに入るとすると、S0の点は、どちらからもdensity-reachableなので、クラスターの極大性に反する。
まあ、こういうことは些細なケースであると言えなくもない。論文にあるアルゴリズムをそのまま実装すると、上記のようなケースで、S0がS1のクラスターに含まれるか、S2のクラスターに含まれるかは、点集合についてのループで最初に出現するcore pointがS1に属するかS2に属するかで決まる。些細なこととは言いつつ、(3)のような、どっちつかずが一番タチが悪くて、これがなければ別にN=1でよいともいえる。こうすれば、単に距離epsを閾値として、データを分類するだけ。実際にクラスタリングすると、大体は(3)のような奴のせいで、分かれてほしいクラスターがくっついたり、それを避けようとすると、くっついてほしいとこが分かれるジレンマに陥る。
x-meansなんかだと、各クラスターがある点の周りにランダムに分布するような時はうまくいって、クラスターの形状は超球状になるのだと思う。現実問題としては、各クラスターがある点の周りに正規分布しているものの、分散が方向に依存して、クラスター形状は超楕円的になる場合が多そうだけれど、分散が非等方的な場合は、どうすればよいのか知らない
一方で、距離関数の作り方の都合とかから、各クラスターがある点の周りに正規分布してくれるとは言えないようなケースも多々あると思う。DBScanは、そういう場合に有効であることを期待されている。
Nとepsをどう決定するのかというのは、当然の疑問で、論文にも最後の方に怪しいことが書いてあるが、とりあえず、自動的に決定できないらしい。
とりあえず、適当に実装
def dbscan(Pts , distance , eps , minPts):
visited = set([])
Noise = set([])
clusters = []
for P in Pts:
if P in visited:continue
visited.add(P)
seeds = [q for q in Pts if distance(P,q)<eps and P!=q]
if len(seeds)<minPts:
Noise.add(P)
else:
cluster = [P]
while 1:
if len(seeds)==0:break
Q = seeds.pop()
if not (Q in visited):
N = [q for q in Pts if distance(Q,q)<eps and Q!=q]
if len(N)>=minPts:
seeds += N
if not (Q in visited) or (Q in Noise):
cluster.append(Q)
visited.add(Q)
if Q in Noise:Noise.remove(Q)
clusters.append( cluster )
return (clusters,Noise)
if__name__=="__main__":
sqrt = __import__("math").sqrt
S0 = [(0.0 , 0.0)]
S1 = [(0.99+0.1*n , 0) for n in xrange(10)]
S2 = [(-0.99-0.1*n , 0) for n in xrange(10)]
S=S0+S1+S2
dbscan(S , lambda x,y:sqrt((x[0]-y[0])**2+(x[1]-y[1])**2) , 1.0 , 3)
謎のベンチマーク
個々の言語で最適化するんじゃなくて、なるべく"同じ"コードで性能を測定。速いコードがほしいなら、Cで書くよって話なので、教科書のアルゴリズムを、そのまま素朴に実装しました的な感じ。
Haskellとpythonのようにパラダイムの違う言語間では、同じ書き方をすると、一方が不利になったりするので、C/C++とSchemeとHaskellはついでで、主にLua/LuaJIT/V8/pythonあたりの差が見たかった。
| 処理系 | 空ループ | 空末尾再帰 | fib(35) | mandelbrot集合 | 挿入ソート | tarai |
|---|---|---|---|---|---|---|
| gcc-4.1.2(optionなし) | 44ms | 205ms | 46ms | 1ms以下 | ||
| gcc-4.1.2(-O3) | 72ms | 30ms | 1ms以下 | |||
| Lua-5.1.4 | 190ms | 1920ms | 4300ms | 559ms | 24390ms | 2939ms |
| LuaJIT 2.0.0 | 10ms | 30ms | 379ms | 31ms | 223ms | 280ms |
| python-2.4 | 795ms | 9768ms | 3703ms | 48185ms | 6756ms | |
| python-2.4+psyco | 19ms | 444ms | 398ms | 2472ms | 244ms | |
| V8 | 85ms | 498ms | 311ms | 49ms | 340ms | |
| Gauche 0.8.14 | 1209ms | 2551ms | 2234ms | 1805ms | ||
| GHC 6.10.1 | 941ms | <1ms | ||||
| LuaJIT 2.0.0(jit.off) | 66ms | 439ms | 1242ms | 160ms | 6000ms | 730ms |
注釈。
・空ループ(10^7回)
地味に面白い結果になった。
gcc(-O3)に結果がないのは、最適化オプションを付けると、空ループが消えてしまうから。で、LuaJITがgccを上回ってるのは、gccがループ変数をレジスターに割り当ててくれないせい。register修飾子とかつけても、あほなコードを吐くので、asmで書いたら、9msとかになった。多分LuaJITも同じようなコード生成しているのかね。速度的にはそんな感じだけど
pythonは若干謎で、関数の中にループを書かないと倍くらい遅いし、psycoも高速化してくれない。人の書いたpythonコード見てると、if __name__=="__main__"の中にコードを直書きするのではなく、main関数を別に作ってる人がいて、趣味の問題かと思ってたけど、こんな差があるのか。けどまあ、psycoだとV8より速い。
あと、Lua/python/JSのいずれでも、whileで回すより、forで回す方が速い。フーン
・空末尾再帰(10^7回)
空ループを末尾再帰版。gaucheの空ループを、どうしようかと思ったのだけど、なるべく類似性の高いコードで比較するという原則に従って別項目として扱うべきかと思った。しかし、殆どの処理系が再帰限界の前に屈したのだった。だらしない。LuaJITは、なんかイカサマしてるんじゃないかと思うくらいの速さ。
・fib(35)
35番目のFibonacci数。再帰とかどんな感じか的な。やっぱC/C++速い。V8/psyco/LuaJITあたりのネイティブコード野郎たちが、まあ大体同じ速度。GHCはだらしない。ネイティブコードの越えられない壁の後に、バイトコード野郎が続く。Gaucheが意外と速くて、まあでもpythonはLuaより遅いねっていう
・mandelbrot集合
よくある80x80のmandelbrot集合の描画。LuaJIT異常に速いな〜っていう。まあ、単純なループの中で、簡単な浮動小数点演算をするだけなので、ループさえ最適化できるなら、C/C++に匹敵できる感じなのかもしれない
・挿入ソート
JSには、ちゃんとした配列があるけど、Luaにはテーブルしかないという差なのかね
・tarai(12,6,0)
GHCに面目躍如させたかった。わけではなく、fib(35)でLuaがGaucheに負けてるのが気になったので、再帰をもう一個。やっぱり、gaucheの方が速い
psycoが結構速い。挿入ソートで、V8に大差で負けてるのは、Pythonのリストがindexアクセスが苦手とかいう点が大きい気がする。なので、pyscoはLanguage shootoutで、よい成績を出せそうではある。
それはともかく、そもそも何をしたかったかというと。LuaとかLuaJITが速いとかよく聞くんだけど、なんでなんかな？的な部分がよく分からんというのが動機だった。Lua/JS/pythonはよく似ていて、今回のコードも、{}がbegin~endになるとかの差はあっても、ほぼ同じコードで書けた(まあ、単にC/C++風ともいえる)。pythonなんかは、coroutineなし、再帰限界あり、制限されたlambdaと、機能的にはLuaに大分劣るとさえ言える。なのに、pythonやJSがLuaより遅いとしたら、その原因は、どこらへんにあるのか?(メモリ消費量もなんでそんな少なくて済むのか?というのも)という点が気になった。ていうか、本当に速いの?速いと言っても、何が得意で何が苦手なの？みたいなのを、もうちょっと知りたかった。
で、とりあえず、LuaJITは異常に速い。JITを使用しないですら、通常のLuaの3~4倍の速度。JITなしでもmandelbrot集合の描画とか、V8より速い。JITすれば最適化したC/C++並の速度が出る場合もある。公式のベンチマーク見ると、ここまで速くなったのは、2.0.0かららしい。再帰はループのようにはいかないっぽいけど、Luaはスタックオーバーフローしました！ってわけにはいかないので、そのへんで限界があるような気がしないでもない。
一体なにが起きているのかと思って、ソース見てみたところ、LuaJITは、単に普通のLuaにJIT機能を付け加えたものなんだと思ってたけど、(少なくとも、LuaJIT2.0.0では)VMからして全く別物っぽく、LuaJIT VMはLuaVMと異なる命令セットを持って、VMのレジスターをCPUのレジスターに割り当てるよう、VMは全部アセンブラで書くという気合の入り方。標準ライブラリもx86で速くなるようにチューニングされているっぽい。VM自体が違うので、実際のところ、何がどのくらい性能向上に貢献しているか分からないし、バイトコードを吐いてくれない/ポータビリティも失われたという点は大きなデメリットではある。
けどまあ、本気出せばここまで速くなるんだな〜というよい目標ではある。尤も、上のベンチマークだと、psycoで10~30倍くらい加速してるけど、実際のアプリケーションでは2~3倍にしか加速しないのが普通だし、Cで書いた部分で時間を食われてれば全く加速しないことも稀に良くあることなので、上のベンチマークに於ける2~3倍くらいの差は、実際のアプリケーションでは誤差程度の差にしかならないかもしれない。Luaはそもそも単体で動かすことはあまり多くないだろうし、LuaJITですらLuaの3~4倍程度だし、果たして、そこまで手間かける価値あるのか？という点は難しいとこであるけど
一方、Luaが速いという点については、どうなんだろう。確かにpythonよりは速そうではある。でも、末尾再帰や再帰の速度は、gaucheに負けている。schemeでは再帰が多用される一方、Luaで再帰を使うことは比較的少ないだろうので、、不利な勝負ではあるけれど。Language shootout見ても、多くの処理系がJITなりAOTなりでネイティブコード生成するので、それらより遅くても仕方ないという感がある。純粋にバイトコード生成して実行する処理系は、RubyとPythonくらい。なので、Luaが速いのではなく、RubyやPythonが遅いのでは?という可能性もある。
とはいえ。ループは速いので、何が起きてるか見てみる。
> cat loop.lua
for i=0,10000000 do end
> luac -l loop.lua
main <loop.lua:0,0> (6 instructions, 24 bytes at 0x9302d10)
0+ params, 4 slots, 0 upvalues, 4 locals, 3 constants, 0 functions
1 [1] LOADK 0 -1 ; 0
2 [1] LOADK 1 -2 ; 10000000
3 [1] LOADK 2 -3 ; 1
4 [1] FORPREP 0 0 ; to 5
5 [1] FORLOOP 0 -1 ; to 5
6 [1] RETURN 0 1
なんかfor文専用の命令がある。手続き型言語のVMだと、これが普通だったりするのだろうか。for文の中で、何かする時は、FORPREPとFORLOOPの間に色々挟まるらしい。一方、whileで書くと
> cat loop2.lua
local c=0
while 1 do c=c+1 if c==10000000 then break end end
> luac -l loop2.lua
main <loop2.lua:0,0> (7 instructions, 28 bytes at 0x8ab4d10)
0+ params, 2 slots, 0 upvalues, 1 local, 3 constants, 0 functions
1 [1] LOADK 0 -1 ; 0
2 [2] ADD 0 0 -2 ; - 1
3 [2] EQ 0 0 -3 ; - 10000000
4 [2] JMP -3 ; to 2
5 [2] JMP 1 ; to 7
6 [2] JMP -5 ; to 2
7 [2] RETURN 0 1
これの実行時間は、上のforで書いた時の3倍くらい遅い。ADDして、EQしてJMPするので、大雑把に一命令で済むのが3命令になる感じなのか。Pythonと比較しても、1.5倍くらい速い程度。結局Luaのループが速いのは、単にc++;if(c==10000000)break;という処理を1命令でやってくれるからという部分がでかい。まあ、Luaでループを多用するであろうことを考えれば、合理的な設計ではある。ちなみに、Squirrelを、ちらっと見た感じでは、FORPREP/FORLOOPに相当する命令は見当たらなかった。
ちなみに、cをlocal変数にしないと、無茶苦茶遅くなって、for文の10倍以上遅くなる。バイトコード見ると、SETGLOBAL/GETGLOBALが沢山挟まって、ステップ数が大幅に増える。pythonには、local修飾子とかないので、Luaと比べて不利なのもこの点かな、という気はする。
ついでに末尾再帰版
> cat loop3.lua
function loop(n)
if n<10000000 then
return loop(n+1)
else
return
end
end
loop(0)
> luac -l loop3.lua
main <loop3.lua:0,0> (6 instructions, 24 bytes at 0x97c2d10)
0+ params, 2 slots, 0 upvalues, 0 locals, 2 constants, 1 function
1 [7] CLOSURE 0 0 ; 0x97c2e40
2 [1] SETGLOBAL 0 -1 ; loop
3 [9] GETGLOBAL 0 -1 ; loop
4 [9] LOADK 1 -2 ; 0
5 [9] CALL 0 2 1
6 [9] RETURN 0 1
function <loop3.lua:1,7> (9 instructions, 36 bytes at 0x97c2e40)
1 param, 3 slots, 0 upvalues, 1 local, 3 constants, 0 functions
1 [2] LT 0 0 -1 ; - 10000000
2 [2] JMP 5 ; to 8
3 [3] GETGLOBAL 1 -2 ; loop
4 [3] ADD 2 0 -3 ; - 1
5 [3] TAILCALL 1 2 0
6 [3] RETURN 1 0
7 [3] JMP 1 ; to 9
8 [5] RETURN 0 1
9 [7] RETURN 0 1Gaucheの場合
gosh> (define (loop n) (if (< n 10000000) (loop (+ n 1)) '()))
loop
gosh> (disasm loop)
main_code (name=loop, code=0x9fd6ed0, size=11, const=1, stack=4):
args: #f
0 CONST 10000000
2 LREF-VAL0-BNLT(0,0) 9 ; (< n 10000000)
4 LREF0 ; n
5 NUMADDI(1) ; (+ n 1)
6 PUSH-GREF-TAIL-CALL(1) #<identifier user#loop>; (loop (+ n 1))
8 RET
9 CONSTN
10 RET
#<undef>Luaと同じような命令が生成される。実行速度も、Luaの1.5倍程度。まあでも、Schemeが、末尾再帰最適化があればループはいらんと言っても、速度的には厳しいものがあるのかな〜という感じではある。とはいえ、LuaJITではJITすると、whileより末尾再帰の方が速くなったりするけど
最後に、Pythonのfor文
>>> def test():
... for _ in xrange(10000000):pass
...
>>> __import__('dis').dis(test)
2 0 SETUP_LOOP 20 (to 23)
3 LOAD_GLOBAL 0 (xrange)
6 LOAD_CONST 1 (10000000)
9 CALL_FUNCTION 1
12 GET_ITER
>> 13 FOR_ITER 6 (to 22)
16 STORE_FAST 0 (_)
19 JUMP_ABSOLUTE 13
>> 22 POP_BLOCK
>> 23 LOAD_CONST 0 (None)
26 RETURN_VALUEもっと速くてもいい気はするけど、謎
色々見てると、Squirrelはどうなんだろうとか、LLVM-luaはどうなんだとか、気になってきたけど、キリがないので、このへんで。
全部貼ると多いので、pythonとLuaのソースコードだけ。pythonに慣れると、endしかない行とか、}しかない行とかが、とてもうざく感じる
続きを読む