Oktaf: hitung jarak antara dua matriks vektor

12

Misalkan saya memiliki dua matriks Nx2, Mx2 masing-masing mewakili N, M 2d vektor. Apakah ada cara sederhana dan bagus untuk menghitung jarak antara masing-masing pasangan vektor (n, m)?

Cara yang mudah namun tidak efisien tentu saja:

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

Jawaban terdekat yang saya temukan adalah bsxfun, digunakan seperti:

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
sumber
Saya melihat ini dan saya tidak bisa melakukan lebih baik daripada membuat vektor perhitungan. Saya pikir perhitungan ini adalah kandidat yang cukup baik untuk menulis fungsi C / Fortran eksternal.
Aron Ahmadia
1
Saya yakin Anda bisa membuat matriks 2xNxM yang Anda isi dengan produk luar, lalu kuadrat setiap entri dan jumlah sepanjang sumbu nol dan akar kuadrat. Dalam Python ini akan terlihat seperti: distance_matrix = (n [:,:, nexaxis] * m [:, newaxis ,:]); distance_matrix = distance_matrix ** 2; distance_matrix = sqrt (distance_matrix.sum (sumbu = 1)); Jika Anda hanya perlu tahu n-vektor terdekat ada banyak cara yang lebih baik untuk melakukan ini!
meawoppl
3
@meawoppl (Baru ke Oktaf) Saya menemukan cara menggunakan paket aljabar linier dalam Oktaf, yang menyediakan cartprod, jadi sekarang saya dapat menulis: (1) x = cartprod(n(:,1), m(:,1)); (2) y = cartprod(n(:,2), m(:,2)); (3) d = sqrt((x(:,1)-x(:,2)).^2+(y(:,1)-y(:,2)).^2) ..yang berjalan jauh lebih cepat!
Kelley van Evert
Bagaimana dengan octave.sourceforge.net/statistics/function/pdist.html
Nemo

Jawaban:

6

Vektorisasi mudah dalam situasi ini menggunakan strategi seperti ini:

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

Berikut adalah contoh yang memvariasikan loop untuk dengan speedup 15x untuk M = 1000 dan 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
sumber
David, senang melihatmu di scicomp! Saya tanpa malu-malu mengedit fragmen kode Anda dan mengembangkannya, tolong kembali jika suntingan saya salah arah dengan menjelaskan apa yang Anda maksud.
Aron Ahmadia
2

Dari Octave 3.4.3 dan kemudian operator - melakukan siaran otomatis (menggunakan bsxfun secara internal). Jadi Anda bisa melanjutkan dengan cara ini.

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

Anda dapat melakukan hal yang sama menggunakan matriks 3d tapi saya rasa ini lebih jelas. D adalah matriks NxM jarak, setiap vektor dalam N terhadap setiap vektor dalam M.

Semoga ini membantu

JuanPi
sumber