Przyspieszanie obliczeń w MATLABie

Przyspieszanie wykonywania algorytmów w środowisku MATLAB to w zasadzie temat rzeka. Drogi do obranego celu są różne. Parallel Computing Toolbox umożliwia zastąpienie pętli for przez parfor, zrównoleglającej obliczenia na wszystkich dostępnych w komputerze rdzeniach obliczeniowych, a także pozwala na wykorzystanie w obliczeniach kart graficznych firmy nVidia. Przy naprawdę zasobożernych zadaniach MATLAB Distributed Computing Server daje możliwość zastosowania do obliczeń zewnętrznego klastra obliczeniowego. Można też wykorzystać MATLAB Coder do wygenerowania kodu C z części naszych algorytmów, co również powinno przełożyć się na prędkość ich działania. Zanim jednak sięgniemy po te zaawansowane i nierzadko kosztowne narzędzia, warto zacząć od przeglądnięcia m-kodu, którym dysponujemy. Nierzadko skrócenie czasu wykonania algorytmu można osiągnąć modyfikując nieznacznie m-kod.

Alokacja pamięci

To prawda, że jedną z zalet programu MATLAB jest to, że coś takiego jak alokacja pamięci może użytkownika nie obchodzić. Rzeczywiście, w „codziennym” użyciu można o tym, w jaki sposób MATLAB operuje na pamięci, w ogóle nie myśleć. Jednak, jeśli czas wykonania naszego algorytmu niebezpiecznie się wydłuża, warto się nad tą kwestią zastanowić. Bardzo wygodna jest w MATLABie swobodna zmiana wielkości macierzy i dopisywanie do niej nowych wartości. Przy różnych okazjach atrakcyjny może okazać się kod typu:

wyniki = [];
For i = 1:n
	wynik = superFajnyAlgorytm();
	wyniki = [wyniki wynik];
end

Kod jest wygodny w zapisie i łatwy do interpretacji przez człowieka, jednak może być problematyczny dla komputera. Jeśli nasze n wynosi powiedzmy 100, to w zasadzie wszystko w porządku. Jeśli jednak n wynosi 10e6 to może pojawić się problem. Rzecz w tym, że jeśli MATLAB ma uaktualnić macierz wyniki o nowy wpis, to za każdym razem musi utworzyć w pamięci miejsce na nową macierz wyniki, skopiować tam zawartość starej macierzy powiększonej o nowy wynik cząstkowy i w końcu zwolnić miejsce zajmowane przez starą macierz wyniki. To naprawdę może bardzo spowolnić działanie algorytmu. Dla przykładu, załóżmy, że przeprowadzamy symulację rzutu 5 kośćmi do gry w taki sposób, że powtarzamy rzuty do czasu, aż po raz 10 wypadnie pięć szóstek. Kod można zapisać tak:

tic
wyniki = [];

for i = 1:10
    wynik = 0;
    while sum(wynik)~=30
        wynik = randi(6,5,1);
        wyniki = [wyniki wynik];
    end
end
toc   

Elapsed time is 43.526072 seconds.

Czas wykonania algorytmu może się oczywiście zmieniać w zależności od naszego szczęścia do gry w kości, niemniej jednak jest on irytująco długi. Zobaczmy co się stanie, jeśli zaalokujemy pamięć na macierz w wynikami.

tic
komplet_szostek = 0;
wyniki = zeros(5,100000);
for i = 1:length(wyniki)
    wyniki(:,i) = randi(6,5,1);
    if sum(wyniki(:,i)) == 30
        komplet_szostek = komplet_szostek + 1;
        if komplet_szostek == 10
            break;
        end
    end
end
wyniki = wyniki(:, 1:i);
toc

Elapsed time is 0.209002 seconds.
Różnica jest kolosalna!

Unikanie pętli for

Unikanie pętli for może w pierwszej chwili wydawać się pomysłem niedorzecznym. Że niby co, wszystko przerabiać na pętle while? Kompletnie nie o to chodzi. MATLAB jest językiem stworzonym do pracy na macierzach i okazuje się, że w niektórych sytuacjach można zrezygnować ze stosowania pętli for, jeśli tylko zaczniemy operować na macierzach w sensowny sposób. Załóżmy, że jest potrzeba wygenerowania sygnału będącego sumą kolejnych, dajmy na to 50 harmonicznych. Można to zrobić oczywiście tak:

fp = 1000000; %1MHz
f0 = 2; %częstotliwość podstawowa
n = 50; % liczba harmonicznych

dt = 1/fp;
t = 0:dt:1-dt;
tic;
y = zeros(1, fp); %wektor wynikowy
for i=1:2:2*n
    x = (1/i) * sin(2*pi*f0*i*t);
    y = y + x;
end
toc

plot(t, y); xlabel('Czas [s]'); ylabel('Amplituda');

Elapsed time is 0.928059 seconds.
A można też inaczej:

fp = 1000000; %1MHz
f0 = 2; %częstotliwość podstawowa
n = 50; % liczba harmonicznych

dt = 1/fp;
t = 0:dt:1-dt;
tic;
i = 1:2:2*n;
A = 1./i;     %wektor kolejnych amplitud
f = (f0.*i)'; % wektor kolejnych częstotliwości
y = A*sin(2*pi*f*t);
plot(t,y);
toc

Elapsed time is 0.781626 seconds.

Różnica na korzyść rozwiązania „nieforowanego” ujawni się tylko wtedy, kiedy korzystać będziemy z komputera wieloprocesorowego (czy też procesora wielordzeniowego), uzbrojonego przynajmniej w 4 niezależne rdzenie obliczeniowe. Im więcej rdzeni, tym lepiej oczywiście. W obliczeniach, w których przekształcane są macierze, angażowane są wszystkie dostępne w komputerze rdzenie obliczeniowe niezależnie od tego, czy mamy zainstalowany Parallel Computing Toolbox czy też nie. Wymaga to nieco wprawy, ale czasami warto z tego skorzystać.