Векторизованное пространственное расстояние в питоне с использованием NumPy - PullRequest
Новая
       37

Векторизованное пространственное расстояние в питоне с использованием NumPy

0 голосов
/ 27 августа 2018

У меня есть пустой массив в python, который содержит множество (10k +) трехмерных точек вершин (векторы с координатами [x, y, z]).Мне нужно рассчитать расстояние между всеми возможными парами этих точек.

это легко сделать с помощью scipy: 

import scipy
D = spdist.cdist(verts, verts)

, но я не могу использовать это из-за политики проекта по введению новых зависимостей.

Итак, я придумал этот наивный код: 

def vert_dist(self, A, B):
    return ((B[0]-A[0])**2+(B[1]-A[1])**2+(B[2]-A[2])**2)**(1.0/2)

# Pairwise distance between verts
#Use SciPy, otherwise use fallback
try:
    import scipy.spatial.distance as spdist
    D = spdist.cdist(verts, verts)
except ImportError:
    #FIXME: This is VERY SLOW:
    D = np.empty((len(verts), len(verts)), dtype=np.float64)
    for i,v in enumerate(verts):
        #self.app.setStatus(_("Calculating distance %d of %d (SciPy not installed => using SLOW AF fallback method)"%(i,len(verts))), True)
        for j in range(i,len(verts)):
            D[j][i] = D[i][j] = self.vert_dist(v,verts[j])

vert_dist () вычисляет трехмерное расстояние между двумя вершинами, а остальная часть кода просто перебирает вершины в одномерном массиве и для каждого еговычисляет расстояние друг от друга в одном и том же массиве и создает двумерный массив расстояний.

Но это очень медленно (в 1000 раз) по сравнению с собственным C-кодом Сципи.Интересно, смогу ли я ускорить его, используя чистый NumPy.по крайней мере, в некоторой степени.

Дополнительная информация: https://github.com/scipy/scipy/issues/9172

Кстати, я пробовал компилятор PyPy JIT, и он был даже медленнее (в 10 раз), чем чистый Python.

ОБНОВЛЕНИЕ: я смог немного ускорить процесс: 

    def vert_dist_matrix(self, verts):
            #FIXME: This is VERY SLOW:
            D = np.empty((len(verts), len(verts)), dtype=np.float64)
            for i,v in enumerate(verts):
                    D[i] = D[:,i] = np.sqrt(np.sum(np.square(verts-verts[i]), axis=1))
            return D

Это устраняет внутренний цикл, вычисляя сразу весь ряд, что делает вещи довольно быстрыми, но все же заметно медленнее, чем scipy.Так что я все еще смотрю на решение @ Divakar

Ответы [ 2 ]

0 голосов
/ 27 августа 2018

Вы можете использовать numpy.linalg.norm: 

from numpy.linalg import norm

a = np.random.rand(10000, 3)
b = np.random.rand(10000, 3)

c = norm(a-b, axis=1)  # will return a np.array of distances

Я не тестировал его, но n=10K случай сработал для меня мгновенно.

0 голосов
/ 27 августа 2018

Существует eucl_dist пакет (отказ от ответственности: я его автор), который в основном содержит два метода для решения проблемы вычисления квадратов евклидовых расстояний, которые более эффективны, чем SciPy's cdist, особенно для больших массивов (с большим количеством столбцов).

Мы будем использовать некоторые коды из его source code, чтобы приспособиться к нашей проблеме, чтобы дать нам два подхода.

Подход № 1

Следуя wiki contents, мы могли бы использовать matrix-multiplication и некоторые NumPy specific implementations для нашего первого подхода,вот так - 

def pdist_squareformed_numpy(a):
    a_sumrows = np.einsum('ij,ij->i',a,a)
    dist = a_sumrows[:,None] + a_sumrows -2*np.dot(a,a.T)
    np.fill_diagonal(dist,0)
    return dist

Подход № 2

Еще один способ - создать «расширенные» версии входных массивов, которые снова подробно обсуждаются в этом исходном коде github.используйте кодовую ссылку для нашего второго подхода, который лучше для меньших столбцов, как здесь, например, так: 

def ext_arrs(A,B, precision="float64"):
    nA,dim = A.shape
    A_ext = np.ones((nA,dim*3),dtype=precision)
    A_ext[:,dim:2*dim] = A
    A_ext[:,2*dim:] = A**2

    nB = B.shape[0]
    B_ext = np.ones((dim*3,nB),dtype=precision)
    B_ext[:dim] = (B**2).T
    B_ext[dim:2*dim] = -2.0*B.T
    return A_ext, B_ext

def pdist_squareformed_numpy_v2(a):
    A_ext, B_ext = ext_arrs(a,a)
    dist = A_ext.dot(B_ext)
    np.fill_diagonal(dist,0)
    return dist

Обратите внимание, что они дают нам евклидовы расстояния в квадрате.Таким образом, для фактических расстояний мы хотим использовать np.sqrt(), если это конечный результат.

Примерные прогоны - 

In [380]: np.random.seed(0)
     ...: a = np.random.rand(5,3)

In [381]: from scipy.spatial.distance import cdist

In [382]: cdist(a,a)
Out[382]: 
array([[0.  , 0.29, 0.42, 0.2 , 0.57],
       [0.29, 0.  , 0.58, 0.42, 0.76],
       [0.42, 0.58, 0.  , 0.45, 0.9 ],
       [0.2 , 0.42, 0.45, 0.  , 0.51],
       [0.57, 0.76, 0.9 , 0.51, 0.  ]])

In [383]: np.sqrt(pdist_squareformed_numpy(a))
Out[383]: 
array([[0.  , 0.29, 0.42, 0.2 , 0.57],
       [0.29, 0.  , 0.58, 0.42, 0.76],
       [0.42, 0.58, 0.  , 0.45, 0.9 ],
       [0.2 , 0.42, 0.45, 0.  , 0.51],
       [0.57, 0.76, 0.9 , 0.51, 0.  ]])

In [384]: np.sqrt(pdist_squareformed_numpy_v2(a))
Out[384]: 
array([[0.  , 0.29, 0.42, 0.2 , 0.57],
       [0.29, 0.  , 0.58, 0.42, 0.76],
       [0.42, 0.58, 0.  , 0.45, 0.9 ],
       [0.2 , 0.42, 0.45, 0.  , 0.51],
       [0.57, 0.76, 0.9 , 0.51, 0.  ]])

Время на 10k точках - 

In [385]: a = np.random.rand(10000,3)

In [386]: %timeit cdist(a,a)
1 loop, best of 3: 309 ms per loop

# Approach #1
In [388]: %timeit pdist_squareformed_numpy(a) # squared eucl distances
1 loop, best of 3: 668 ms per loop

In [389]: %timeit np.sqrt(pdist_squareformed_numpy(a)) # actual eucl distances
1 loop, best of 3: 812 ms per loop

# Approach #2
In [390]: %timeit pdist_squareformed_numpy_v2(a) # squared eucl distances
1 loop, best of 3: 237 ms per loop

In [391]: %timeit np.sqrt(pdist_squareformed_numpy_v2(a)) # actual eucl distances
1 loop, best of 3: 395 ms per loop

Второй подход кажется близким к cdist по производительности!

...