Numpy Broadcast to performing Euclidean distance vectorized

Oto oryginalne zmienne wejściowe:

A = np.array([[1,1,1,1],[2,2,2,2]])
B = np.array([[1,2,3,4],[1,1,1,1],[1,2,1,9]])
A
# array([[1, 1, 1, 1],
#        [2, 2, 2, 2]])
B
# array([[1, 2, 3, 4],
#        [1, 1, 1, 1],
#        [1, 2, 1, 9]])

A jest tablicą 2x4. B to tablica 3x4.

Chcemy obliczyć euklidesową operację macierzy odległości w jednej całkowicie wektoryzowanej operacji, gdzie dist[i,j] zawiera odległość między I-tą instancją w A I J-tą instancją w B. zatem dist jest 2x3 W tym przykładzie.

Odległość

Może być zapisany z numpy jako

dist = np.sqrt(np.sum(np.square(A-B))) # DOES NOT WORK
# Traceback (most recent call last):
#   File "<stdin>", line 1, in <module>
# ValueError: operands could not be broadcast together with shapes (2,4) (3,4)

Jednakże, jak pokazano powyżej, problem polega na tym, że elementowa operacja odejmowania A-B obejmuje niezgodne rozmiary tablic, w szczególności 2 i 3 w pierwszym wymiarze.

A has dimensions 2 x 4
B has dimensions 3 x 4

Aby wykonać elementowe odejmowanie, musimy pad A lub B, aby spełnić reguły numpy ' ego. Wybieram pad a z dodatkowym wymiarem, aby stał się 2 x 1 x 4, co pozwala na wyrównanie wymiarów tablic do nadawania. Aby dowiedzieć się więcej o numpy broadcasting, zobacz samouczek w instrukcji scipy i ostatni przykład w ten tutorial .

Możesz wykonać wypełnienie za pomocą wartości np.newaxis lub polecenia np.reshape. Pokazuję oba poniżej:

# First approach is to add the extra dimension to A with np.newaxis
A[:,np.newaxis,:] has dimensions 2 x 1 x 4
B has dimensions                     3 x 4

# Second approach is to reshape A with np.reshape
np.reshape(A, (2,1,4)) has dimensions 2 x 1 x 4
B has dimensions                          3 x 4

Jak widzisz, użycie obu metod pozwoli na wyrównanie wymiarów. Użyję pierwszego podejścia z np.newaxis. Tak więc teraz będzie to działać, aby utworzyć tablicę A-B, która jest tablicą 2x3x4:

diff = A[:,np.newaxis,:] - B
# Alternative approach:
# diff = np.reshape(A, (2,1,4)) - B
diff.shape
# (2, 3, 4)

Teraz możemy umieścić to wyrażenie różnicy w dist równaniu, aby uzyskać ostateczne wynik:

dist = np.sqrt(np.sum(np.square(A[:,np.newaxis,:] - B), axis=2))
dist
# array([[ 3.74165739,  0.        ,  8.06225775],
#        [ 2.44948974,  2.        ,  7.14142843]])

Zauważ, że sum jest nad axis=2, co oznacza, że pobieramy sumę nad trzecią osią tablicy 2x3x4 (gdzie ID osi zaczyna się od 0).

Jeśli macierze są małe, powyższe polecenie będzie działać poprawnie. Jeśli jednak masz duże tablice, możesz napotkać problemy z pamięcią. Zauważ, że w powyższym przykładzie numpy wewnętrznie stworzył tablicę 2x3x4 do wykonywania transmisji. Jeśli uogólnimy A na wymiary a x z i B na Wymiary b x z, to numpy wewnętrznie utworzy tablicę a x b x z do nadawania.

Możemy uniknąć tworzenia tej pośredniej tablicy, wykonując pewną matematyczną manipulację. Ponieważ obliczasz odległość euklidesową jako sumę kwadratów różnic, możemy skorzystać z matematycznego faktu, że suma kwadratów różnic może zostać przepisana.

Zauważ, że średnioterminowy termin obejmuje sumę ponadpierwiastkiem mnożenie. Ta suma ponad multipleksacje są lepiej znane jako produkt kropkowy. Ponieważ A i B są macierzami, to operacja ta jest w rzeczywistości mnożeniem macierzy. Możemy więc przepisać powyższe jako:

Możemy wtedy napisać następujący kod numpy:

threeSums = np.sum(np.square(A)[:,np.newaxis,:], axis=2) - 2 * A.dot(B.T) + np.sum(np.square(B), axis=1)
dist = np.sqrt(threeSums)
dist
# array([[ 3.74165739,  0.        ,  8.06225775],
#        [ 2.44948974,  2.        ,  7.14142843]])

Zauważ, że powyższa odpowiedź jest dokładnie taka sama jak poprzednia implementacja. Ponownie, zaletą jest to, że nie musimy tworzyć pośredniej tablicy 2x3x4 do nadawania.

Dla kompletności, sprawdźmy dwukrotnie, czy wymiary każdego summanda w threeSums pozwoliły na nadawanie.

np.sum(np.square(A)[:,np.newaxis,:], axis=2) has dimensions 2 x 1
2 * A.dot(B.T) has dimensions                               2 x 3
np.sum(np.square(B), axis=1) has dimensions                 1 x 3

Tak więc, zgodnie z oczekiwaniami, ostateczna dist tablica ma wymiary 2x3.

To użycie iloczynu kropkowego zamiast sumy pierwiastkowego mnożenia jest również omówione w ten tutorial .

Miałem ten sam problem ostatnio pracując z Deep learning(stanford cs231n, Assignment1), ale kiedy używałem

 np.sqrt((np.square(a[:,np.newaxis]-b).sum(axis=2)))

Wystąpił błąd

MemoryError

To znaczy, że zabrakło mi pamięci(w rzeczywistości,to wyprodukowało szereg 500*5000*1024 w middle.It jest taki ogromny!)

Aby zapobiec temu błędowi, możemy użyć formuły upraszczającej:

Kod:

import numpy as np
aSumSquare = np.sum(np.square(a),axis=1);
bSumSquare = np.sum(np.square(b),axis=1);
mul = np.dot(a,b.T);
dists = np.sqrt(aSumSquare[:,np.newaxis]+bSumSquare-2*mul)

Ta funkcjonalność jest już zawarta w modułu przestrzennego scipy i polecam jej użycie, ponieważ będzie ona wektoryzowana i wysoce zoptymalizowana pod maską. Ale, jak wynika z drugiej odpowiedzi, są sposoby, w jakie możesz to zrobić sam.

import numpy as np
a = np.array([[1,1,1,1],[2,2,2,2]])
b = np.array([[1,2,3,4],[1,1,1,1],[1,2,1,9]])
np.sqrt((np.square(a[:,np.newaxis]-b).sum(axis=2)))
# array([[ 3.74165739,  0.        ,  8.06225775],
#       [ 2.44948974,  2.        ,  7.14142843]])
from scipy.spatial.distance import cdist
cdist(a,b)
# array([[ 3.74165739,  0.        ,  8.06225775],
#       [ 2.44948974,  2.        ,  7.14142843]])

Używając numpy.linalg.Norma działa również dobrze z nadawaniem. Podanie wartości całkowitej dla axis użyje normy wektorowej, która domyślnie jest normą euklidesową.

import numpy as np

a = np.array([[1,1,1,1],[2,2,2,2]])
b = np.array([[1,2,3,4],[1,1,1,1],[1,2,1,9]])
np.linalg.norm(a[:, np.newaxis] - b, axis = 2)

# array([[ 3.74165739,  0.        ,  8.06225775],
#       [ 2.44948974,  2.        ,  7.14142843]])