有效更新点之间的距离答案

【问题标题】：Efficiently updating distance between points有效更新点之间的距离
【发布时间】：2017-08-05 06:14:10
【问题描述】：

我有一个包含 n 行（观察）和 p 列（特征）的数据集：

import numpy as np
from scipy.spatial.distance import pdist, squareform
p = 3
n = 5
xOld = np.random.rand(n * p).reshape([n, p])

我有兴趣在真正具有 n x (n-1)/2 唯一值的 nxn 矩阵中获取这些点之间的距离

sq_dists = pdist(xOld, 'sqeuclidean')
D_n = squareform(sq_dists)

现在假设我收到了N 的其他观察结果并想更新D_n。一种非常低效的方法是：

N = 3
xNew = np.random.rand(N * p).reshape([N, p])
sq_dists = pdist(np.row_stack([xOld, xNew]), 'sqeuclidean')
D_n_N = squareform(sq_dists)

但是，考虑到 n ~ 10000 和 N ~ 100，这将是多余的。我的目标是使用D_n 更有效地获得D_n_N。为了做到这一点，我将 D_n_N 划分如下。我已经有D_n 并且可以计算B [N x N]。但是，我想知道是否有一种好方法可以在没有一堆 for 循环的情况下计算 A（或 A 转置）并最终构造 D_n_N

D_n (n x n)    A [n x N]
A.T [N x n]    B [N x N]

提前致谢。

【问题讨论】：

标签： python arrays numpy scipy euclidean-distance

【解决方案1】：

相当有趣的问题！好吧，在获得解决方案的过程中，我必须在这里学习一些新东西。

涉及的步骤：

首先，我们在这里介绍新的点。因此，我们需要使用cdist 来获得新旧点之间的平方欧几里得距离。这些将被容纳在新输出的两个块中，一个在旧距离的正下方，一个在旧距离的右侧。
我们还需要计算新点中的pdist，并将其square-formed块放在新对角区域的尾部。

新D_n_N 的示意图如下所示：

[   D_n      cdist.T
  cdist      New pdist squarefomed]

总而言之，实现看起来是这样的 -

cdists = cdist( xNew, xOld, 'sqeuclidean')

n1 = D_n.shape[0]
out = np.empty((n1+N,n1+N))
out[:n1,:n1] = D_n
out[n1:,:n1] = cdists
out[:n1,n1:] = cdists.T
out[n1:,n1:] = squareform(pdist(xNew, 'sqeuclidean'))

运行时测试

方法-

# Original approach
def org_app(D_n, xNew):
    sq_dists = pdist(np.row_stack([xOld, xNew]), 'sqeuclidean')
    D_n_N = squareform(sq_dists)
    return D_n_N    

# Proposed approach
def proposed_app(D_n, xNew, N):
    cdists = cdist( xNew, xOld, 'sqeuclidean')    
    n1 = D_n.shape[0]
    out = np.empty((n1+N,n1+N))
    out[:n1,:n1] = D_n
    out[n1:,:n1] = cdists
    out[:n1,n1:] = cdists.T
    out[n1:,n1:] = squareform(pdist(xNew, 'sqeuclidean'))
    return out

时间安排 -

In [102]: # Setup inputs
     ...: p = 3
     ...: n = 5000
     ...: xOld = np.random.rand(n * p).reshape([n, p])
     ...: 
     ...: sq_dists = pdist(xOld, 'sqeuclidean')
     ...: D_n = squareform(sq_dists)
     ...: 
     ...: N = 3000
     ...: xNew = np.random.rand(N * p).reshape([N, p])
     ...: 

In [103]: np.allclose( proposed_app(D_n, xNew, N), org_app(D_n, xNew))
Out[103]: True

In [104]: %timeit org_app(D_n, xNew)
1 loops, best of 3: 541 ms per loop

In [105]: %timeit proposed_app(D_n, xNew, N)
1 loops, best of 3: 201 ms per loop

【讨论】：

太棒了。这正是我想要的。一个单独的问题。定义一个诸如out = np.empty((n1+N,n1+N)) 之类的空矩阵并对其进行填充，或者使用hstack 和vstack 函数将这四个部分放在一起是否更有效？在 MATLAB 中，预先分配一个空矩阵是有意义的。我想知道哪种方法在 Python 中是最佳的。感谢您的详细回答。
的原因是空的更好！最好的办法是在开始的时候预留一个很大的空矩阵，然后在每一步填充。
@ahoosh 如果您有很多分配步骤，那么初始化一次然后执行这些分配是有意义的。我想即使与这三个 hstack/vstack 组合进行比较，初始化和分配可能会更好。

【解决方案2】：

只需使用 cdist ：

D_OO=cdist(xOld,xOld)

D_NN=cdist(xNew,xNew)
D_NO=cdist(xNew,xOld)
D_ON=cdist(xOld,xNew) # or D_NO.T

最后：

D_=vstack((hstack((D_OO,D_ON)),(hstack((D_NO,D_NN)))))

【讨论】：