2014-12-29 5 views
4

두 numpy 배열 (x와 y) 사이의 Mahalanobis 거리를 계산하는 NumPy 방법을 찾고 있습니다. 다음 코드는 Scipy의 cdist 함수를 사용하여 정확하게 계산할 수 있습니다. 이 함수는 필자의 경우 불필요한 matix를 계산하므로 NumPy를 사용하여 더 직접 계산할 수 있습니다.NumPy를 사용하여 Mahalanobis 거리를 계산하십시오.

import numpy as np 
from scipy.spatial.distance import cdist 

x = np.array([[[1,2,3,4,5], 
       [5,6,7,8,5], 
       [5,6,7,8,5]], 
       [[11,22,23,24,5], 
       [25,26,27,28,5], 
       [5,6,7,8,5]]]) 
i,j,k = x.shape 

xx = x.reshape(i,j*k).T 


y = np.array([[[31,32,33,34,5], 
       [35,36,37,38,5], 
       [5,6,7,8,5]], 
       [[41,42,43,44,5], 
       [45,46,47,48,5], 
       [5,6,7,8,5]]]) 


yy = y.reshape(i,j*k).T 

results = cdist(xx,yy,'mahalanobis') 
results = np.diag(results) 
print results 



[ 2.28765854 2.75165028 2.75165028 2.75165028 0.   2.75165028 
    2.75165028 2.75165028 2.75165028 0.   0.   0.   0. 
    0.   0.  ] 

내 재판 :

VI = np.linalg.inv(np.cov(xx,yy)) 

print np.sqrt(np.dot(np.dot((xx-yy),VI),(xx-yy).T)) 

이 사람이이 방법을 수정할 수 있을까요? 여기

그것을 위해 공식은 :

http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.spatial.distance.mahalanobis.html#scipy.spatial.distance.mahalanobis

+0

[1,11]과 [31,41] 사이의 마하라 누스 거리를 계산하고 싶었습니다. [2,22]와 [32,42], ... 등등. – Borys

+0

'scipy'에서의 구현은 순수한 파이썬 코드입니다. 그들의 접근 방식을 비교할 수 있습니다. 두 벡터 사이의 마할 라 노비스 거리를 계산하려면 여기를 참조하십시오. https://github.com/scipy/scipy/blob/6a7327e8bb8248b2ea165180bc602edf1ab33dda/scipy/spatial/distance.py#L508-541 거리 계산은 각 관측 벡터를 통해 반복해야하는 관측 행렬. – cel

+0

예, 저 소스에서 계산을 시도했지만 파이썬에 대한 지식이 부족하기 때문에 아직 완성되지 않았습니다. 내 재판 좀 봐 주 시겠어요? – Borys

답변

9

내가 당신의 문제가 당신의 공분산 행렬의 건설에있다 생각합니다. 시도 :

X = np.vstack([xx,yy]) 
V = np.cov(X.T) 
VI = np.linalg.inv(V) 
print np.diag(np.sqrt(np.dot(np.dot((xx-yy),VI),(xx-yy).T))) 

출력 :

A = np.dot((xx-yy),VI) 
B = (xx-yy).T 
n = A.shape[0] 
D = np.empty(n) 
for i in range(n): 
    D[i] = np.sqrt(np.sum(A[i] * B[:,i])) 

편집 :

[ 2.28765854 2.75165028 2.75165028 2.75165028 0.   2.75165028 
    2.75165028 2.75165028 2.75165028 0.   0.   0.   0. 
    0.   0.  ] 

이 암시 적으로 여기에서 생성 된 중간 배열없이이 작업을 수행하려면, 당신은 파이썬 하나에 대한 C 루프를 희생해야 할 수도 있습니다 실제로, np.einsum 부두를 사용하면 파이썬 루프를 제거하고 (시스템에서 84.3 μs에서 2.9 μs까지) 많은 속도를 낼 수 있습니다 :

D = np.sqrt(np.einsum('ij,ji->i', A, B)) 

편집은 : @Warren Weckesser가 지적 하듯이, einsum도 중간 AB 배열을 멀리 할 수 ​​있습니다 :

delta = xx - yy 
D = np.sqrt(np.einsum('nj,jk,nk->n', delta, VI, delta)) 
+0

시도해 주셔서 감사합니다, upvoted. 사실, 속도 향상을 위해 불필요한 계산을 줄이기 위해 np.diag를 사용하지 않기로했습니다. – Borys

+0

아, 알겠습니다. 내 게시물을 수정했습니다. – xnx

+0

정말 고마워요. C 루프가 무슨 뜻인지 아십니까? – Borys

0

einsum만큼이나 빠른 또 다른 간단한 해결책

e = xx-yy 
X = np.vstack([xx,yy]) 
V = np.cov(X.T) 
p = np.linalg.inv(V) 
D = np.sqrt(np.sum(np.dot(e,p) * e, axis = 1)) 
관련 문제