Szybkie algorytmy do wyszukiwania parami odległości euklidesowej

Question

Wiem, że program matlab ma wbudowaną funkcję pdist, która obliczy odległości parami. Jednak moja macierz jest tak duża, że ​​jej 60000 na 300 i matlab zabraknie pamięci. To jest kontynuacja Matlab euclidean pairwise square distance function.

Czy istnieje sposób obejścia tej nieefektywności obliczeniowej. Próbowałem ręcznie kodować obliczenia odległości parami i zwykle trwa to cały dzień (czasami od 6 do 7 godzin).

Każda pomoc jest bardzo doceniana!

Answer 1

Komputery nie są nieskończenie duże lub nieskończenie szybkie. Ludzie myślą, że mają dużo pamięci, szybki procesor, więc po prostu tworzą coraz większe problemy, a potem zastanawiają się, dlaczego ich problem działa powoli. Faktem jest, że nie jest to nieefektywność obliczeniowa. TO JEST przeciążony procesor.

Jak zaznacza Oli w komentarzu, istnieje coś takiego jak 2e9 wartości do obliczenia, nawet przy założeniu, że obliczasz tylko górną lub dolną połowę macierzy odległości. (6e4^2/2 to około 2e9.) Będzie to wymagać przechowywania około 16 gigabajtów pamięci RAM, przy założeniu, że tylko JEDNA kopia tablicy zostanie utworzona w pamięci. Jeśli twój kod jest niechlujny, możesz go z łatwością podwoić lub potroić. Jak tylko wejdziesz do pamięci wirtualnej, rzeczy stają się znacznie wolniejsze.

Nie wystarczy mieć duży problem, aby szybko uruchomić. Aby naprawdę pomóc, musimy wiedzieć, ile pamięci RAM jest dostępne. Czy to jest problem z pamięcią wirtualną? Czy używasz 64-bitowego MATLAB na procesorze, który może obsłużyć wszystkie potrzebne pamięci RAM?

Answer 2

Nie mogłem się oprzeć pokusie grania. Stworzyłem Matlab mex C file o nazwie pdistc, który implementuje parową odległość euklidesową dla pojedynczej i podwójnej precyzji. Na moim komputerze używającym Matlab R2012b i R2015a jest on o 20% szybszy niż pdist (i podstawową funkcję pomocniczą pdistmex) dla dużych wejść (np. 60 000 na 300).

Jak już podkreślono, problem ten jest zasadniczo ograniczony przez pamięć i prosisz go o wiele. Mój kod mex C wykorzystuje minimalną pamięć poza tym, co jest potrzebne do wyjścia. Porównując jego użycie pamięci do tego z pdist, wygląda na to, że oba są praktycznie takie same. Innymi słowy, pdist nie używa dużo dodatkowej pamięci. Twój problem z pamięcią jest prawdopodobnie w pamięci zużytej przed wywołaniem pdist (możesz użyć clear, aby usunąć duże tablice?) Lub po prostu dlatego, że próbujesz rozwiązać duży problem obliczeniowy na małym sprzęcie.

Tak więc, moja funkcja prawdopodobnie nie będzie w stanie całkowicie zaoszczędzić pamięci, ale być może będziesz mógł użyć innej wbudowanej funkcji. Możesz obliczyć fragmenty ogólnego wektora odległości parami. Coś takiego:

m = 6e3; 
n = 3e2; 
X = rand(m,n); 
sz = m*(m-1)/2; 

for i = 1:m:sz-m 
    D = pdistc(X', i, i+m); % mex C function, X is transposed relative to pdist 
    ...      % Process chunk of pairwise distances 
end 

ten jest znacznie wolniejszy (10 razy lub więcej) i ta część mojego kodu C nie jest zoptymalizowany dobrze, ale pozwoli znacznie mniejsze zużycie pamięci – zakładając, że nie potrzebują cała tablica w tym samym czasie. Zauważ, że możesz zrobić to samo znacznie skuteczniej z pdist (lub pdistc), tworząc pętlę, w której przekazałeś bezpośrednio podzestawy X, a nie wszystkie.

Jeśli masz 64-bitowy Intel Mac, nie musisz kompilować, ponieważ umieściłem plik binarny .mexmaci64, ale w przeciwnym razie musisz dowiedzieć się, jak skompilować kod dla swojego komputera. Nie mogę ci z tym pomóc.Możliwe, że nie możesz go skompilować lub że wystąpią problemy ze zgodnością, które musisz rozwiązać, samodzielnie edytując kod. Możliwe też, że są błędy, a kod zawiesi Matlab. Zauważ, że możesz uzyskać nieco inne wartości wyjściowe w stosunku do pdist z różnicami między nimi w zakresie epsilon maszyny (eps). pdist może, ale nie musi, robić fantazyjne rzeczy, aby uniknąć przepełnień dla dużych danych wejściowych i innych problemów numerycznych, ale należy pamiętać, że mój kod nie.

Dodatkowo stworzyłem prosty pure Matlab implementation. Jest znacznie wolniejszy niż kod mex, ale wciąż jest szybszy niż naiwna implementacja lub kod znaleziony w pdist.

Wszystkie pliki can be found here. Archiwum ZIP zawiera wszystkie pliki. Jest licencjonowany na BSD. Zapraszam do optymalizacji (wypróbowałem wywołania BLAS i OpenMP w kodzie C bez żadnego pożytku, – Może jakiś wskaźnik magii lub GPU/OpenCL mógłby go jeszcze przyspieszyć). Mam nadzieję, że może być pomocna dla ciebie lub kogoś innego.

Answer 3

W moim systemie Poniżej znajduje się najszybciej (nawet szybciej niż kod C pdistc przez @horchler):

function [ mD ] = CalcDistMtx (mX)  
    vSsqX = sum(mX .^ 2); 
    mD = sqrt(bsxfun(@plus, vSsqX.', vSsqX) - (2 * (mX.' * mX)));  
end 

Musisz bardzo dobrze dostrojony kod C, aby pokonać tę, myślę.

P. S.
Korzystanie Matlaba pdist dla porównania: squareform(pdist(mX.')) jest równoważna CalcDistMtx(mX).
Mianowicie wejście powinno zostać transponowane.