Efektywne Obliczanie euklidesowej macierzy odległości za pomocą Numpy

Mam zbiór punktów w przestrzeni 2-wymiarowej i muszę obliczyć odległość od każdego punktu do drugiego punktu.

Mam stosunkowo małą liczbę punktów, może najwyżej 100. Ale ponieważ muszę to robić często i szybko, aby określić relacje między tymi ruchomymi punktami, a ponieważ jestem świadomy, że iteracja przez punkty może być tak zła, jak złożoność O (N^2), szukam sposobów, aby skorzystać z magii matrycy numpy (lub scipy).

Jak to jest w moim kodzie, współrzędne każdego obiektu są przechowywane w jego klasie. Jednakże, mogę również zaktualizować je w tablicy numpy, gdy aktualizuję współrzędną klasy.

class Cell(object):
    """Represents one object in the field."""
    def __init__(self,id,x=0,y=0):
        self.m_id = id
        self.m_x = x
        self.m_y = y

Przyszło mi do głowy stworzenie euklidesowej macierzy odległości, aby zapobiec powielaniu, ale być może masz mądrzejszą strukturę danych.

Jestem otwarty na wskazówki do sprytnych algorytmów.

Zauważam również, że istnieją podobne pytania dotyczące odległości euklidesowej i numpy, ale nie znalazłem żadnego, który bezpośrednio odnosi się do tej kwestii skutecznego wypełniania macierzy pełnej odległości.

Author: Wes Modes, 2014-03-28
Source

4 answers

Możesz skorzystać z complex typu: 

# build a complex array of your cells
z = np.array([complex(c.m_x, c.m_y) for c in cells])

Pierwsze rozwiązanie 

# mesh this array so that you will have all combinations
m, n = np.meshgrid(z, z)
# get the distance via the norm
out = abs(m-n)

Drugie rozwiązanie

Siatkowanie jest główną ideą. Ale numpy jest sprytny, więc nie musisz generować m & n. Wystarczy obliczyć różnicę używając transponowanej wersji z. Siatka jest wykonywana automatycznie: 

out = abs(z[..., np.newaxis] - z)

Trzecie rozwiązanie

I jeśli {[11] } jest bezpośrednio ustawiona jako tablica 2-wymiarowa, możesz użyć z.T zamiast dziwnej z[..., np.newaxis]. Więc w końcu, Twój kod będzie wyglądał tak: 

z = np.array([[complex(c.m_x, c.m_y) for c in cells]]) # notice the [[ ... ]]
out = abs(z.T-z)

Przykład 

>>> z = np.array([[0.+0.j, 2.+1.j, -1.+4.j]])
>>> abs(z.T-z)
array([[ 0.        ,  2.23606798,  4.12310563],
       [ 2.23606798,  0.        ,  4.24264069],
       [ 4.12310563,  4.24264069,  0.        ]])

Jako dopełnienie możesz później usunąć duplikaty, przyjmując górny trójkąt: 

>>> np.triu(out)
array([[ 0.        ,  2.23606798,  4.12310563],
       [ 0.        ,  0.        ,  4.24264069],
       [ 0.        ,  0.        ,  0.        ]])

Niektóre benchmarki 

>>> timeit.timeit('abs(z.T-z)', setup='import numpy as np;z = np.array([[0.+0.j, 2.+1.j, -1.+4.j]])')
4.645645342274779
>>> timeit.timeit('abs(z[..., np.newaxis] - z)', setup='import numpy as np;z = np.array([0.+0.j, 2.+1.j, -1.+4.j])')
5.049334864854522
>>> timeit.timeit('m, n = np.meshgrid(z, z); abs(m-n)', setup='import numpy as np;z = np.array([0.+0.j, 2.+1.j, -1.+4.j])')
22.489568296184686
 35
Author: Kiwi,
Warning: date(): Invalid date.timezone value 'Europe/Kyiv', we selected the timezone 'UTC' for now. in /var/www/agent_stack/data/www/doraprojects.net/template/agent.layouts/content.php on line 54
2014-03-29 08:26:37

Jeśli nie potrzebujesz pełnej macierzy odległości, lepiej będzie użyć KD-tree. Rozważmy scipy.spatial.cKDTree LUB sklearn.neighbors.KDTree. Dzieje się tak dlatego, że KD-drzewo kan znajduje k-najbliżej sąsiadów w czasie O(n log n) i dlatego unikasz O(n**2) złożoności obliczeń wszystkich n przez N odległości.

 9
Author: Sturla Molden,
Warning: date(): Invalid date.timezone value 'Europe/Kyiv', we selected the timezone 'UTC' for now. in /var/www/agent_stack/data/www/doraprojects.net/template/agent.layouts/content.php on line 54
2014-03-29 20:37:09

Oto jak możesz to zrobić używając numpy: 

import numpy as np

x = np.array([0,1,2])
y = np.array([2,4,6])

# take advantage of broadcasting, to make a 2dim array of diffs
dx = x[..., np.newaxis] - x[np.newaxis, ...]
dy = y[..., np.newaxis] - y[np.newaxis, ...]
dx
=> array([[ 0, -1, -2],
          [ 1,  0, -1],
          [ 2,  1,  0]])

# stack in one array, to speed up calculations
d = np.array([dx,dy])
d.shape
=> (2, 3, 3)

Teraz pozostaje obliczenie normy L2 wzdłuż osi 0 (jak omówiono tutaj):

(d**2).sum(axis=0)**0.5
=> array([[ 0.        ,  2.23606798,  4.47213595],
          [ 2.23606798,  0.        ,  2.23606798],
          [ 4.47213595,  2.23606798,  0.        ]])
 8
Author: shx2,
Warning: date(): Invalid date.timezone value 'Europe/Kyiv', we selected the timezone 'UTC' for now. in /var/www/agent_stack/data/www/doraprojects.net/template/agent.layouts/content.php on line 54
2014-03-28 19:21:49

Jake Vanderplas podaje ten przykład używając broadcastingu w Python Data Science Handbook , który jest bardzo podobny do tego, co zaproponował @shx2. 

import numpy as np
rand = random.RandomState(42)
X = rand.rand(3, 2)  
dist_sq = np.sum((X[:, np.newaxis, :] - X[np.newaxis, :, :]) ** 2, axis = -1)

dist_sq
array([[0.        , 0.18543317, 0.81602495],
       [0.18543317, 0.        , 0.22819282],
       [0.81602495, 0.22819282, 0.        ]])
 8
Author: Rich Pauloo,
Warning: date(): Invalid date.timezone value 'Europe/Kyiv', we selected the timezone 'UTC' for now. in /var/www/agent_stack/data/www/doraprojects.net/template/agent.layouts/content.php on line 54
2019-01-06 00:28:19