Pythonでxy座標上の2点間の距離をforループを使わずに計算する方法
機械学習でカーネル法やらクラスタリングをやっていると、何かと「2サンプル(点)間の距離」を計算することが多いです。ここではより一般的に「Pythonで2点間の距離をforループを使わずに行列(テンソル)計算として求める方法」を見ていきます。
目次
動機
カーネル法やクラスタリングをしていると、ニューラルネットワークの出力層(潜在空間)に対して、サンプル間の距離を求めるということがよく目にします。ディープラーニングになると損失関数の中で計算しないといけないので、forループが使うのが厳しい、行列(テンソル)計算でなんとかできればめっちゃ嬉しいというシチュエーションがあります。
実はKaggleのとあるカーネルを見ていたら、TensorFlowの場合ですが、行列計算であっさり求めていたので目からうろこでした。今回はTensorFlowではなくNumpyの問題として一般化します。
問題:日本の都市間の距離
今カーネル法やら機械学習の話はおいておくとしましょう。もっと一般的な問題の設定にします。
問題:以下のデータは、日本の主要都市の「緯度経度」を表したものです(データはWikipediaから)。ここからforループを使わずに、行列・テンソル計算を使って、都市間の距離の行列を求めてください。距離は2点間のユークリッド距離として計算してください。
import numpy as np
cities = np.array([[43.062083, 141.354389], # 札幌
[38.268056, 140.869722], # 仙台
[35.689472, 139.69175], # 東京
[35.181389, 136.906389], # 名古屋
[34.69375, 135.502111], # 大阪
[34.38525, 132.455306], # 広島
[33.590139, 130.401722], # 福岡
[26.213083, 127.678056]]) # 那覇
ただし、今問題を単純にするため、緯度1度あたりの距離と経度1度あたりの距離が同じとします(これは本来赤道以外異なります)。緯度1度あたりの距離は約111km、東京での経度1度あたりの距離は約91kmなので、今回は平均を取って緯度・経度とも1度あたり101kmで計算するものとします。
腕に覚えがある方はぜひ自力で解いてみてください。ちなみに途中の計算はテンソルになります。
答え
自分もよくわからなかったのでいきなり答え書きます。自力で解きたい人はがんばってください。自力でできた人はとてもえらい。
all_diffs = np.expand_dims(cities, axis=1) - np.expand_dims(cities, axis=0)
degree_distance = np.sqrt(np.sum(all_diffs ** 2, axis=-1))
distance = (101 * degree_distance).astype(np.int32)
print(" 札幌 仙台 東京 名古屋 大阪 広島 福岡 那覇")
print(distance)
出力は次のようになります。こんな距離表を飛行機で見たことあるのではないでしょうか。
札幌 仙台 東京 名古屋 大阪 広島 福岡 那覇
[[ 0 486 763 913 1031 1255 1462 2191]
[ 486 0 286 507 651 935 1158 1804]
[ 763 286 0 285 434 742 961 1545]
[ 913 507 285 0 150 456 676 1299]
[1031 651 434 150 0 309 527 1165]
[1255 935 742 456 309 0 222 956]
[1462 1158 961 676 527 222 0 794]
[2191 1804 1545 1299 1165 956 794 0]]
うーん、沖縄遠い。福岡から行っても、東京-札幌ぐらいの距離あるんですね。
ちなみにこれはだいぶいい加減な計算です。なぜなら日本のような中緯度地域では、緯度と経度の1度は違うからです。正しい値は国土地理院のデータを見ましょう。
正しい値は、東京~札幌が831.0km、東京~仙台が304.9km、東京~名古屋が259.1km、東京~大阪が395.9km、東京~広島が675.1km、東京~福岡が880.6km、東京~那覇が1553.6kmだそうです。東京~福岡のようにほぼ経度方向にだけ移動する場合は、緯度1度=経度1度の仮定による誤差がかなり大きいですね。逆に、東京~那覇のような地図上でほぼ対角線での移動は、距離が長くても誤差が小さいのがわかります。
プログラムの解説
若干脱線してしまいましたが、答えの解説をしましょう。このプログラムのキーは差分の計算です。
all_diffs = np.expand_dims(cities, axis=1) - np.expand_dims(cities, axis=0)
これはPythonのブロードキャスティングを使ったとても美しい書き方なので、ここがミソです。ぜひ味わってください。元のcitiesという配列は(8,2)というshapeなので、左側のexpand_dimsは(8,1,2)、右側のexpand_dimsは(1,8,2)になります。
なんとなく見えてきたでしょうか? Pythonブロードキャスティングは、shapeが1の場合は多いほうの次元分コピーされるので、(8,1,2)→(8,8,2)、(1,8,2)→(8,8,2)という次元にブロードキャスティングされます。軸を増やす→ブロードキャスティングだけを使って、転置(transpose)や内積(matmul)を一切使わずに点同士の差分を求めることができました。
あとは特に大丈夫ですよね。要素を2乗して、3つ目の軸の和を取って(8,8)の行列にする。あとは平方根を取れば無事、(緯度経度単位での)距離が出てくるわけです。あとはこれに101を掛け、キロメートル単位に変換すれば終わりです。
まとめ
np.expand_dims()とブロードキャスティングを使えば、forループを使わずに点同士の距離を求めることができました。ベクトル同士の場合は、内積や転置でそれっぽくできますが、行列以上の点(サンプル)同士の計算はしんどくなってくるので、ぜひこのやり方をマスターしてみてください。
Shikoan's ML Blogの中の人が運営しているサークル「じゅ~しぃ~すくりぷと」の本のご案内
技術書コーナー
北海道の駅巡りコーナー