仅使用NumPy计算马氏距离

Question

我正在寻找计算两个核阵列（x和y）之间的马氏距离的核方法。 下面的代码可以使用Scipy的cdist函数正确计算相同的值。因为这个函数在我的例子中计算不必要的matix，所以我想要更直接的方法，只使用NumPy计算它。

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

Answer 1

另一个简单的解决方案，它和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))

Answer 2

我认为你的问题在于协方差矩阵的构造。尝试：

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)))

输出：

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

要在不隐式创建中间数组的情况下执行此操作，您可能需要为Python牺牲一个C循环：

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]))

编辑：实际上，使用np.einsum巫术，您可以删除Python循环并大大加快它的速度（在我的系统中，从84.3微秒到2.9微秒）：

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

编辑：正如@Warren Weckesser指出的，einsum也可以用来消除中间的A和B数组：

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