Monthly Archives: 2月 2016

scikit-learnのk最近傍の結果の加工について

sciki-learnの最近傍法について

scikit-learnにはk最近傍(K-Nearest Neighbor)を計算するNearestNeighborsというクラスがあるが、これはこんな感じで使う:

>>> from sklearn.neighbors import NearestNeighbors
>>> nn = NearestNeighbors(n_neighbors=3)
>> x = np.array([[0,0],[0,1],[2,2],[3,3]])
>> nn.fit(x)
earestNeighbors(algorithm='auto', leaf_size=30, metric='minkowski',
        metric_params=None, n_neighbors=3, p=2, radius=1.0)
>> distances, indices  = nn.kneighbors(x)
>> distances
rray([[ 0.        ,  1.        ,  2.82842712],
      [ 0.        ,  1.        ,  2.23606798],
      [ 0.        ,  1.41421356,  2.23606798],
      [ 0.        ,  1.41421356,  3.60555128]])
>> indices
rray([[0, 1, 2],
      [1, 0, 2],
      [2, 3, 1],
      [3, 2, 1]], dtype=int64)

ここで入力は、(0,0), (0,1), (2,2), (3,3)(順に点0,…,点3とする)として与えられている。クラスのインスタンス化のときに与えられているn_neighbors=3は、k最近傍のkを表していて、つまりここでは各点から一番近い3点ずつを計算する。

結果は2つの変数distancesとindicesに入れられていて、indicesは各点から近い順に点のインデックスが入っている。行i(以下、行列の行を上から順に行0, 行1と表現する)には点iから近い順に点のインデックスが入っている。ここで行1が[1,0,2]なのは、点1から近い順に点1、点0、点2であるということである。

問題提起

レコメンデーションシステムや情報検索の分野では、自分自身を除いて最近傍をとりたいことがよくある。つまり点Xから近いもののリストに点Xが入っていてほしくない。そういうものを取得するにはどうすればよいだろうか。distanceもindicesも両方とも自分自身を除いたものにしたい。

一瞬、自分自身が必ず一番近いから一番左の列を削ればいいのではと思うかもしれないが、話はそう単純ではない。座標が同じ点が含まれる場合、このようになる:

>>> nn = NearestNeighbors(n_neighbors=3)
>>> x=[[0,0],[0,0],[1,1]]
>>> nn.fit(x)
NearestNeighbors(algorithm='auto', leaf_size=30, metric='minkowski',
         metric_params=None, n_neighbors=3, p=2, radius=1.0)
>>> distances, indices = nn.kneighbors(x)
>>> indices
array([[0, 1, 2],
       [0, 1, 2],
       [2, 1, 0]], dtype=int64)

この場合は、行1の一番左に0がくるが、この行からも1を取り除きたい。

さらにやっかいなのは、行iにiが含まれないこともある。

>>> x=np.zeros((4,2))
>>> nn.fit(x)
NearestNeighbors(algorithm='auto', leaf_size=30, metric='minkowski',
         metric_params=None, n_neighbors=3, p=2, radius=1.0)
>>> distances, indices = nn.kneighbors(x)
>>> indices
array([[0, 1, 2],
       [0, 1, 2],
       [0, 1, 2],
       [0, 1, 2]], dtype=int64)

これは極端な例だが、行3には3が含まれていない。つまりここで2近傍を計算するときは行3は一番右を消して[0,1]にならないといけない。

やりたいことの仕様をまとめると次のようになる。

  • 行列indicesの各行について、行iにiという値があればその値を削除する
  • 行iにiという値がなければ一番右の値を削除する
  • 以上の削除により元のindicesより1列少ない行列を取得する
  • 行列distancesについても、indicesで削除された要素と同じ位置の要素を削除して1列少ない行列を取得する

このことをNumpyの機能を駆使してfor文を使わずにやりたい。何しろ、for文使ったら負けだと思っているので。

解説

ここではなぜそのようなコードになるのかという理由を解説するので、コードだけを知りたい場合は読み飛ばして「結論」を参照されたい。

次のような人工的なindicesとdistancesを用意して計算してみる。

>>> indices=np.array([[0,1,3,4] for _ in range(3)])
>>> distances=np.arange(10,22).reshape(3,4)
>>> indices
array([[0, 1, 3, 4],
       [0, 1, 3, 4],
       [0, 1, 3, 4]])
>>> distances
array([[10, 11, 12, 13],
       [14, 15, 16, 17],
       [18, 19, 20, 21]])

方針としては、削除するべき要素がFalseでそれ以外をTrueにしたマスク行列を用意して、インデキシングを利用して抽出する。

行iで要素がiであるものをFalseにするには次のような計算をする。

>>> mask = indices != np.arange(indices.shape[0])[:,np.newaxis]
>>> mask
array([[False,  True,  True,  True],
       [ True, False,  True,  True],
       [ True,  True,  True,  True]], dtype=bool)

ここでは演算子!=がブロードキャストされて各要素の比較が行われている。indicesの行数をnとすると、演算子!=の右辺はn行1列で0~n-1が順に入った行列を意味している。indicesの各要素と、右辺の行列の対応する行の値が比較されている。

次に、一番下の行がすべてTrueなのでこの一番右をFalseにしたい。つまり一般にはすべてがTrueな行の一番右をFalseに書き換えたい。まず各行について「すべてがTrueか?」はこれで判定できる。

>>> mask.all(axis=1)
array([False, False,  True], dtype=bool)

ここでallは要素の論理積をとるメソッドで、引数axis=1が指定されているのは横方向に計算することを意味している。

この結果のnotをとって、maskの一番右の列との論理積をとり、maskの一番右の列に格納しなおせばよい。つまりこうする。

>> mask[:,-1] &= np.logical_not(mask.all(axis=1))
>> mask
rray([[False,  True,  True,  True],
      [ True, False,  True,  True],
      [ True,  True,  True, False]], dtype=bool)

これでマスク行列ができたので、これをindicesとdistancesに適用すればいいのだが、2次元配列にインデキシングすると1次元配列になってしまうのでreshapeする必要がある。

>>> shape = (indices.shape[0], indices.shape[1] - 1)
>>> indices = indices[mask].reshape(shape)
>>> distances = distances[mask].reshape(shape)
>>> indices
array([[1, 3, 4],
       [0, 3, 4],
       [0, 1, 3]])
>>> distances
array([[11, 12, 13],
       [14, 16, 17],
       [18, 19, 20]])

これでほしい結果が得られた。

結論

次のコードでやりたいことができる。ただし、indicesとdistancesを上書きしている。

mask = indices != np.arange(indices.shape[0])[:,np.newaxis]
mask[:,-1] &= np.logical_not(mask.all(axis=1))
shape = (indices.shape[0], indices.shape[1] - 1)
indices = indices[mask].reshape(shape)
distances = distances[mask].reshape(shape)

スクリプト言語の流儀〜C言語、Python、Ruby速度比較

2年くらい前の記事だが、こんな記事が気になったので補足する。

俺の言語がこんなに遅いわけがない!? 〜C, Java, PHP, Python, Rubyによるプログラミング言語 速度比較〜

確かにスクリプト言語は遅い。ただし、C言語で一般的なアルゴリズムをそのまま使って計算したとしたら。

スクリプト言語にはスクリプト言語なりの流儀があり、慣れている人は違う方法で実装する。

ちょうどPythonについて、同じようなことを先日ソフトウェアジャパンというイベントでしゃべって来たのでこのスライドも参照されたい。

ここで改めてベンチマークをとってみたい。まずはPythonから。Pythonではこの手の数値計算をするのはNumpyというライブラリを使うのが普通である。Numpyで1からnが入った配列を用意し、そのメソッドsum()で合計を計算する。

結果の格納先になぜ辞書型を使うのかなど、気に入らないところもあるが、和の計算以外は上記リンク先のベンチマークコードをそのまま使う。

N = 10000
import numpy as np


def sumup(n):
    return np.arange(1, n + 1).sum()


def main():
    print("python with numpy start.")
    result = {}
    for count in range(1, N + 1):
        result[count - 1] = sumup(count)
    print("python with numpy end.")

main()

Pythonの数値計算でNumpyを使うのは充分に一般的だと思うが、Python本体に含まれてないライブラリを使って「Pythonのベンチマーク」と主張することに違和感を感じる人がいるかもしれないので、Numpyを使わない版も実験してみる。その場合、1からnが入ったリストを用意してsumという関数を適用するのが普通かと思う。

N = 10000


def sumup(n):
    return sum(range(1, n + 1))


def main():
    print("python without numpy start.")
    result = {}
    for count in range(1, N + 1):
        result[count - 1] = sumup(count)
    print("python without numpy end.")

main()

次にRubyで実験してみる。Rubyの配列にはinjectというメソッドがあるのでそれを利用してみる。

N = 10000

def sumup(n)
  return (1..n).inject{ |sum, x| sum + x }
end

def main()
  result = []
  puts "ruby start."
  (1..N).each do |count|
    result[count-1] = sumup(count)
  end
  puts "ruby end."
end

main()

C言語は、参照先ブログのコードをそのまま使い、最適化なし(オプション-O0)と最適化あり(-O3)を試してみる。

#include <stdio.h>
#define N 10000

long sumup(int n) {
  int i;
  long sum = 0;
  for (i=1; i<=n; i++) {
    sum += i;
  }
  return sum;
}

int main() {
  int count;
  long result[N];
  printf("c start.\n");
  for (count = 1; count <= N; count++) {
    result[count-1] = sumup(count);
  }
  printf("c end.\n");
  return 0;
}

参照先ブログにあるJavaとPHPはここでは無視する。

ベンチマーク結果

AWSのt2.microで、それぞれのバージョンは以下のとおり
OS: Ubuntu 14.04.3 LTS
Python: 3.5.1
Ruby: 1.9.3p484
gcc: 4.8.4

参照先にあるようにtimeコマンドで測定した。

 

言語 計算時間(msec)
Python(Numpyあり) 167
Python(Numpyなし) 739
Ruby 4159
C言語(最適化なし) 185
C言語(最適化あり) 1

考察

最適化付きのC言語が速すぎるのは、最適化しすぎて何もやっていないからである。アセンブリコードを吐かせて確認したが、main関数からsumup関数の呼び出し自体が消えている。なので、これと他を比べるのはあまり意味がないかと思う。

Python+Numpyは爆速であることがわかる。Numpyなしでもそれなりの速さは出ている。この結果だけをみるとRubyが遅いようだが、にしても参照先ブログからかなり速くなっており、受ける印象が随分違うと思う。

まとめ

  • C言語と同じような発想で同じアルゴリズムをそのまま実装するとスクリプト言語は非常に遅い
  • しかしスクリプト言語はスクリプト言語なりの流儀があり、慣れている人はどうすれば速くなるか知っている。そのコードはC言語とかけ離れることもある。

更新履歴:
2016/2/11 誤植の修正と若干の加筆