Oktawa: oblicz odległość między dwiema macierzami wektorów

Załóżmy, że mam dwie macierze Nx2, Mx2 reprezentujące odpowiednio wektory N, M 2d. Czy istnieje prosty i dobry sposób obliczenia odległości między każdą parą wektorów (n, m)?

Łatwy, ale nieefektywny sposób to oczywiście:

d = zeros(N, M);
for i = 1:N,
  for j = 1:M,
    d(i,j) = norm(n(i,:) - m(j,:));
  endfor;
endfor;

Najbliższa odpowiedź, jaką znalazłem, jest następująca bsxfun:

bsxfun(inline("x-y"),[1,2,3,4],[3;4;5;6])

ans =
  -2 -1  0  1
  -3 -2 -1  0
  -4 -3 -2 -1
  -5 -4 -3 -2

performance octave vectorization

— Kelley van Evert
źródło

Spojrzałem na to i nie mogłem nic lepszego niż wektoryzować obliczenia. Myślę, że to obliczenie jest całkiem dobrym kandydatem do napisania zewnętrznej funkcji C / Fortran.

— Aron Ahmadia

Założę się, że możesz zrobić matrycę 2xNxM, którą wypełniasz produktem zewnętrznym, a następnie kwadrat każdego z wpisów i zsumować wzdłuż osi zerowej i pierwiastka kwadratowego. W Pythonie wyglądałoby to następująco: distance_matrix = (n [:,:, nexaxis] * m [:, newaxis ,:]); Distance_matrix = Distance_matrix ** 2; distance_matrix = sqrt (distance_matrix.sum (oś = 1)); Jeśli potrzebujesz tylko poznać najbliższe n-wektory, są na to znacznie lepsze sposoby!

— meawoppl

@meawoppl (nowy w Octave) Dowiedziałem się, jak korzystać z pakietu algebry liniowej w Octave, który zapewnia cartprod, więc teraz mogę pisać: (1) x = cartprod(n(:,1), m(:,1)); (2) y = cartprod(n(:,2), m(:,2)); (3) .. który d = sqrt((x(:,1)-x(:,2)).^2+(y(:,1)-y(:,2)).^2) działa znacznie szybciej!

— Kelley van Evert

Co powiesz na octave.sourceforge.net/statistics/function/pdist.html

— Nemo

Odpowiedzi:

Wektoryzacja jest prosta w takich sytuacjach przy użyciu strategii takiej jak ta:

eN = ones(N,1);
eM = ones(M,1);
d  = sqrt(eM*n.^2' - 2*m*n' + m.^2*eN');

Oto przykład, który wektoryzuje pętlę for z 15-krotnym przyspieszeniem dla M = 1000 i N = 2000.

n = rand(N,2);
m = rand(M,2);
eN = ones(N,2);
eM = ones(2,M);

tic;
d_vect  = sqrt(eN*m.^2' - 2*n*m' + n.^2*eM);
vect_time = toc;

tic;
for i=1:N
  for j=1:M
     d_for(i,j) = norm(n(i,:)-m(j,:));
  end
end
for_time = toc; 

assert(norm(d_vect-d_for) < 1e-10*norm(d_for))

— David Bindel
źródło

David, miło cię widzieć na scicomp! Bezwstydnie zredagowałem twój fragment kodu i rozwinąłem go trochę. Cofnij, jeśli moje zmiany poszły w złym kierunku od wyjaśnienia tego, co zamierzałeś.

— Aron Ahmadia

Od wersji Octave 3.4.3 i późniejszych operator wykonuje automatyczną transmisję (używa bsxfun wewnętrznie). Możesz więc postępować w ten sposób.

Dx = N(:,1) - M(:,1)';
Dy = N(:,2) - M(:,2)';
D = sqrt (Dx.^2 + Dy.^2);

Możesz zrobić to samo za pomocą matrycy 3D, ale myślę, że to jest bardziej jasne. D jest macierzą odległości NxM, każdy wektor w N względem każdego wektora w M.

Mam nadzieję że to pomoże

— JuanPi
źródło