NumPyを使用してユークリッド距離を計算するにはどうすればよいですか?
-
05-07-2019 - |
質問
3Dには2つのポイントがあります:
(xa, ya, za)
(xb, yb, zb)
そして、距離を計算したい:
dist = sqrt((xa-xb)^2 + (ya-yb)^2 + (za-zb)^2)
これをNumPyで、またはPython全般で行う最良の方法は何ですか私が持っている:
a = numpy.array((xa ,ya, za))
b = numpy.array((xb, yb, zb))
解決
dist = numpy.linalg.norm(a-b)
他のヒント
SciPyにはそのための関数があります。 ユークリッドと呼ばれています。
例:
from scipy.spatial import distance
a = (1, 2, 3)
b = (4, 5, 6)
dst = distance.euclidean(a, b)
複数の距離を一度に計算することに興味がある人のために、 perfplot (私の小さなプロジェクト)。判明したこと
a_min_b = a - b
numpy.sqrt(numpy.einsum('ij,ij->i', a_min_b, a_min_b))
a
および b
の行の距離を最速で計算します。これは実際には1行だけにも当てはまります!
プロットを再現するコード:
import matplotlib
import numpy
import perfplot
from scipy.spatial import distance
def linalg_norm(data):
a, b = data
return numpy.linalg.norm(a-b, axis=1)
def sqrt_sum(data):
a, b = data
return numpy.sqrt(numpy.sum((a-b)**2, axis=1))
def scipy_distance(data):
a, b = data
return list(map(distance.euclidean, a, b))
def mpl_dist(data):
a, b = data
return list(map(matplotlib.mlab.dist, a, b))
def sqrt_einsum(data):
a, b = data
a_min_b = a - b
return numpy.sqrt(numpy.einsum('ij,ij->i', a_min_b, a_min_b))
perfplot.show(
setup=lambda n: numpy.random.rand(2, n, 3),
n_range=[2**k for k in range(20)],
kernels=[linalg_norm, scipy_distance, mpl_dist, sqrt_sum, sqrt_einsum],
logx=True,
logy=True,
xlabel='len(x), len(y)'
)
def dist(x,y):
return numpy.sqrt(numpy.sum((x-y)**2))
a = numpy.array((xa,ya,za))
b = numpy.array((xb,yb,zb))
dist_a_b = dist(a,b)
さまざまなパフォーマンスノートを使用して簡単な答えを説明したいと思います。 np.linalg.normはおそらくあなたが必要とする以上のことをします:
dist = numpy.linalg.norm(a-b)
まず-この関数はリストを処理し、すべての値を返すように設計されています。 pA
から一連のポイント sP
までの距離を比較するには:
sP = set(points)
pA = point
distances = np.linalg.norm(sP - pA, ord=2, axis=1.) # 'distances' is a list
いくつかのことを覚えておいてください:
- Python関数呼び出しは高価です。
- [通常] Pythonは名前検索をキャッシュしません。
そう
def distance(pointA, pointB):
dist = np.linalg.norm(pointA - pointB)
return dist
見た目ほど無邪気ではありません。
>>> dis.dis(distance)
2 0 LOAD_GLOBAL 0 (np)
2 LOAD_ATTR 1 (linalg)
4 LOAD_ATTR 2 (norm)
6 LOAD_FAST 0 (pointA)
8 LOAD_FAST 1 (pointB)
10 BINARY_SUBTRACT
12 CALL_FUNCTION 1
14 STORE_FAST 2 (dist)
3 16 LOAD_FAST 2 (dist)
18 RETURN_VALUE
まず-呼び出すたびに、「np」のグローバルルックアップ、「linalg」のスコープルックアップを行う必要があります。 " norm"のスコープルックアップ、および関数を単に呼び出すオーバーヘッドは、数十のPython命令に相当します。
最後に、結果を保存し、リロードのためにリロードするために2つの操作を無駄にしました...
改善の最初のパス:検索を高速化し、ストアをスキップします
def distance(pointA, pointB, _norm=np.linalg.norm):
return _norm(pointA - pointB)
はるかに合理化されました:
>>> dis.dis(distance)
2 0 LOAD_FAST 2 (_norm)
2 LOAD_FAST 0 (pointA)
4 LOAD_FAST 1 (pointB)
6 BINARY_SUBTRACT
8 CALL_FUNCTION 1
10 RETURN_VALUE
ただし、関数呼び出しのオーバーヘッドは依然としてある程度の作業になります。そして、自分で数学を行う方が良いかどうかを判断するためのベンチマークを行う必要があります。
def distance(pointA, pointB):
return (
((pointA.x - pointB.x) ** 2) +
((pointA.y - pointB.y) ** 2) +
((pointA.z - pointB.z) ** 2)
) ** 0.5 # fast sqrt
一部のプラットフォームでは、 ** 0.5
は math.sqrt
よりも高速です。走行距離は異なる場合があります。
****高度なパフォーマンスノート。
なぜ距離を計算するのですか?表示のみが目的の場合は、
print("The target is %.2fm away" % (distance(a, b)))
一緒に移動します。ただし、距離を比較したり、範囲チェックを行ったりする場合は、有用なパフォーマンス観測を追加したいと思います。
2つのケースを考えてみましょう:距離によるソート、または範囲の制約を満たすアイテムへのリストのカリング。
# Ultra naive implementations. Hold onto your hat.
def sort_things_by_distance(origin, things):
return things.sort(key=lambda thing: distance(origin, thing))
def in_range(origin, range, things):
things_in_range = []
for thing in things:
if distance(origin, thing) <= range:
things_in_range.append(thing)
最初に覚えておく必要があるのは、ピタゴラスを使用して計算することです距離( dist = sqrt(x ^ 2 + y ^ 2 + z ^ 2)
)ので、多くの sqrt
呼び出しを行っています。数学101:
dist = root ( x^2 + y^2 + z^2 )
:.
dist^2 = x^2 + y^2 + z^2
and
sq(N) < sq(M) iff M > N
and
sq(N) > sq(M) iff N > M
and
sq(N) = sq(M) iff N == M
要するに:実際にX ^ 2ではなくX単位の距離が必要になるまで、計算の最も難しい部分を排除できます。
# Still naive, but much faster.
def distance_sq(left, right):
""" Returns the square of the distance between left and right. """
return (
((left.x - right.x) ** 2) +
((left.y - right.y) ** 2) +
((left.z - right.z) ** 2)
)
def sort_things_by_distance(origin, things):
return things.sort(key=lambda thing: distance_sq(origin, thing))
def in_range(origin, range, things):
things_in_range = []
# Remember that sqrt(N)**2 == N, so if we square
# range, we don't need to root the distances.
range_sq = range**2
for thing in things:
if distance_sq(origin, thing) <= range_sq:
things_in_range.append(thing)
素晴らしい、両方の関数はもはや高価な平方根をしません。それははるかに高速です。また、ジェネレーターに変換することでin_rangeを改善できます。
def in_range(origin, range, things):
range_sq = range**2
yield from (thing for thing in things
if distance_sq(origin, thing) <= range_sq)
これは、次のようなことをしている場合に特に利点があります。
if any(in_range(origin, max_dist, things)):
...
しかし、あなたがやろうとしている非常に次のものが距離を必要とする場合、
for nearby in in_range(origin, walking_distance, hotdog_stands):
print("%s %.2fm" % (nearby.name, distance(origin, nearby)))
タプルの生成を検討する:
def in_range_with_dist_sq(origin, range, things):
range_sq = range**2
for thing in things:
dist_sq = distance_sq(origin, thing)
if dist_sq <= range_sq: yield (thing, dist_sq)
これは、範囲チェックを連鎖する可能性がある場合に特に役立ちます(距離を再度計算する必要がないため、「Xの近くでYのNm以内にあるものを見つける」)。
しかし、 things
の非常に大きなリストを検索していて、それらの多くが考慮に値しないと予想される場合はどうでしょうか。
実際には非常に単純な最適化があります:
def in_range_all_the_things(origin, range, things):
range_sq = range**2
for thing in things:
dist_sq = (origin.x - thing.x) ** 2
if dist_sq <= range_sq:
dist_sq += (origin.y - thing.y) ** 2
if dist_sq <= range_sq:
dist_sq += (origin.z - thing.z) ** 2
if dist_sq <= range_sq:
yield thing
これが有用かどうかは、「もの」のサイズに依存します。
def in_range_all_the_things(origin, range, things):
range_sq = range**2
if len(things) >= 4096:
for thing in things:
dist_sq = (origin.x - thing.x) ** 2
if dist_sq <= range_sq:
dist_sq += (origin.y - thing.y) ** 2
if dist_sq <= range_sq:
dist_sq += (origin.z - thing.z) ** 2
if dist_sq <= range_sq:
yield thing
elif len(things) > 32:
for things in things:
dist_sq = (origin.x - thing.x) ** 2
if dist_sq <= range_sq:
dist_sq += (origin.y - thing.y) ** 2 + (origin.z - thing.z) ** 2
if dist_sq <= range_sq:
yield thing
else:
... just calculate distance and range-check it ...
また、dist_sqの生成を検討してください。ホットドッグの例は次のようになります。
# Chaining generators
info = in_range_with_dist_sq(origin, walking_distance, hotdog_stands)
info = (stand, dist_sq**0.5 for stand, dist_sq in info)
for stand, dist in info:
print("%s %.2fm" % (stand, dist))
matplotlib.mlabで「dist」関数を見つけましたが、十分に便利だとは思いません。
参照用にここに投稿しています。
import numpy as np
import matplotlib as plt
a = np.array([1, 2, 3])
b = np.array([2, 3, 4])
# Distance between a and b
dis = plt.mlab.dist(a, b)
次のように実行できます。どれだけ高速かはわかりませんが、NumPyを使用していません。
from math import sqrt
a = (1, 2, 3) # Data point 1
b = (4, 5, 6) # Data point 2
print sqrt(sum( (a - b)**2 for a, b in zip(a, b)))
ベクトルを減算してから、innerproductを減算できます。
例に従って、
a = numpy.array((xa, ya, za))
b = numpy.array((xb, yb, zb))
tmp = a - b
sum_squared = numpy.dot(tmp.T, tmp)
result sqrt(sum_squared)
これは単純なコードであり、理解しやすいです。
np.dot
(ドット積)が好きです:
a = numpy.array((xa,ya,za))
b = numpy.array((xb,yb,zb))
distance = (np.dot(a-b,a-b))**.5
定義したとおりに a
と b
を使用すると、次も使用できます。
distance = np.sqrt(np.sum((a-b)**2))
素敵なワンライナー:
dist = numpy.linalg.norm(a-b)
ただし、速度が懸念される場合は、お使いのマシンで実験することをお勧めします。私のマシンでは、スクエアの **
演算子で math
ライブラリの sqrt
を使用すると、1ライナーのNumPyソリューションよりもはるかに高速であることがわかりました。 。
この単純なプログラムを使用してテストを実行しました:
#!/usr/bin/python
import math
import numpy
from random import uniform
def fastest_calc_dist(p1,p2):
return math.sqrt((p2[0] - p1[0]) ** 2 +
(p2[1] - p1[1]) ** 2 +
(p2[2] - p1[2]) ** 2)
def math_calc_dist(p1,p2):
return math.sqrt(math.pow((p2[0] - p1[0]), 2) +
math.pow((p2[1] - p1[1]), 2) +
math.pow((p2[2] - p1[2]), 2))
def numpy_calc_dist(p1,p2):
return numpy.linalg.norm(numpy.array(p1)-numpy.array(p2))
TOTAL_LOCATIONS = 1000
p1 = dict()
p2 = dict()
for i in range(0, TOTAL_LOCATIONS):
p1[i] = (uniform(0,1000),uniform(0,1000),uniform(0,1000))
p2[i] = (uniform(0,1000),uniform(0,1000),uniform(0,1000))
total_dist = 0
for i in range(0, TOTAL_LOCATIONS):
for j in range(0, TOTAL_LOCATIONS):
dist = fastest_calc_dist(p1[i], p2[j]) #change this line for testing
total_dist += dist
print total_dist
私のマシンでは、 math_calc_dist
は numpy_calc_dist
よりもはるかに高速に実行されます: 1.5秒対 23.5秒。
fastest_calc_dist
と math_calc_dist
の測定可能な差を得るには、 TOTAL_LOCATIONS
を6000に上げる必要がありました。その後、 fastest_calc_dist
〜50秒かかり、 math_calc_dist
は〜60秒かかります。
numpy.sqrt
と numpy.square
を試すこともできますが、どちらも私のマシンの math
の選択肢よりも低速です。
Python 2.6.6でテストを実行しました。
Pythonでリストとして表される2つのポイントが与えられた場合のPythonのユークリッド距離の簡潔なコードを次に示します。
def distance(v1,v2):
return sum([(x-y)**2 for (x,y) in zip(v1,v2)])**(0.5)
import numpy as np
from scipy.spatial import distance
input_arr = np.array([[0,3,0],[2,0,0],[0,1,3],[0,1,2],[-1,0,1],[1,1,1]])
test_case = np.array([0,0,0])
dst=[]
for i in range(0,6):
temp = distance.euclidean(test_case,input_arr[i])
dst.append(temp)
print(dst)
import math
dist = math.hypot(math.hypot(xa-xb, ya-yb), za-zb)
式は簡単に使用できます
distance = np.sqrt(np.sum(np.square(a-b)))
実際には、ピタゴラスの定理を使用して距離を計算するだけで、&#916; x、&#916; y、&#916; zの2乗を加算し、結果をルート化するだけです。
多次元空間のユークリッド距離を計算します:
import math
x = [1, 2, 6]
y = [-2, 3, 2]
dist = math.sqrt(sum([(xi-yi)**2 for xi,yi in zip(x, y)]))
5.0990195135927845
最初に2つの行列の差を見つけます。次に、numpyの乗算コマンドで要素ごとの乗算を適用します。その後、要素ごとに乗算された新しい行列の合計を見つけます。最後に、合計の平方根を見つけます。
def findEuclideanDistance(a, b):
euclidean_distance = a - b
euclidean_distance = np.sum(np.multiply(euclidean_distance, euclidean_distance))
euclidean_distance = np.sqrt(euclidean_distance)
return euclidean_distance