2-1+%A6V%B6q%A4%C6%B9B%BA%E2%BBP+JIT+%A5[%B3t

在 6.5 版之前,MATLAB 程式環境是一個傳統的解譯器(Interpreter),在執行 MATLAB的程式碼時,會進行下列動作:

  1. 逐列對程式碼轉換為 p-code,這是 MATLAB 可以快速解讀的格式。
  2. 對產生的 p-code 進行逐列執行。

由於在執行每一列 p-code 時,都還包含一些「經常開銷」(Overhead),如果此列執行需要大量運算,這些經常開銷就會顯得很微小而不會拖累程式執行速度。反之,如果此列只是簡單的運算,那經常開銷的比例就會相對提高,更糟的是,如果此列是放在迴圈內,那麼程式碼的執行速度就會被這些經常開銷大幅拖慢。

因此簡單的說,在 MATLAB 6.5 版之前,若要加快執行速度,就要盡量不使用迴圈(例如 for 迴圈、while 迴圈等),而盡量改用向量化運算(Vectorized operations),以降低經常開銷的比例,尤其是要盡量使用 MATLAB 內建的指令來完成運算,盡量避免使用迴圈。

在 6.5 版後,MATLAB 引進了 JIT (Just-In-Time) 的編譯技術,簡稱 JIT 加速器(JIT-Accelerator),大幅地提高了迴圈的執行效率,也使得一般程式碼的撰寫,可以使用程式設計師最熟悉的方式,或是使用最能代表演算法精神的方式來呈現,而不必為了提高效率、勉強使用向量化運算,造成程式碼維護的困難。在這種情況下,判斷是否需要使用向量化程式碼,就變成不是那麼明顯。

JIT 加速器的主要功能,可以列出如下:

  1. JIT 加速器可將每一列 MATLAB 程式碼直接轉成「原機指令」(Native Machine Instructions),而不再直接使用 p-code,因此省去大部分 p-code 會遇到的經常開銷。
  2. JIT 加速器也會對 Intel X86 為主的 Windows 及 Linux 進行程式碼最佳化。

事實上,JIT 加速器可以分成兩部分:

  1. Just-In-Time Code Generation:將 p-code 轉成原機指令。
  2. Runtime Type Analysis:對資料型態的分析,以便用於 Just-In-Time Code Generation.

由於 MATLAB 的變數資料型態可以隨時改變,因此 6.5 版之前的 MATLAB 在執行 p-code 時,要花一些時間來決定變數的型態,這就是屬於 p-code 的經常開銷的一部份,會拖慢執行速度。而在 JIT 加速器在產生原機指令前,會先對每一列 MATLAB 程式碼進行資料型態的分析,其所把握的原則是:如果某一列程式碼已經被分析過,那麼若遇到同樣內容的一列程式碼,其變數的資料型態和維度都應該和已分析過的前一列程式碼一樣。根據這個原則,只要我們在迴圈中的變數都保持一樣的資料型態和維度,那麼 JIT 加速器就可以大幅地提高執行效率。

讀者或許會問:既然已經有了 JIT 加速器,我們是否還要使用向量化運算呢?基本的判斷原則如下:

  1. 如果所用的向量化運算並非根基於「以空間換取時間」,那還是要盡量使用(尤其是要盡量使用內建指令),因為這種向量化運算不需要額外的空間來暫存資訊,因此效率會比 JIT 加速器更好。(請見以下兩個範例。)
  2. 如果所用的向量化運算是根基於「以空間換取時間」,就必須經過仔細測試,才能決定是否值得使用。
  3. 碰到較複雜的資料結構(例如稀疏矩陣及維度超過 3 維的陣列),JIT 加速器並不會轉成原機指令,因此在這種情況下,也要盡量使用向量化運算。

以下將舉出幾個向量化運算的範例來進行探討,以便讓讀者瞭解「向量化運算」和「JIT 加速」這兩者對 MATLAB 執行效率的交互影響。(由於這些範例和 MATLAB 的版本有密切的關係,因此我們在每個範例都先印出 MATLAB 的版本,以供參考。)

首先我們使用「元素對元素」(Element by Element)的矩陣乘法,來說明如何使用「不需要額外空間」的向量化運算來提昇計算速度:

Example 1: 02-程式碼與記憶體之最佳化/matMultiply01.mfprintf('MATLAB version = %s\n', version); ns = 10*(1:30); for j=1:length(ns) n = ns(j); a = rand(n); b = rand(n); % 第一種方法:for-loop operation c1 = zeros(n); tic for p = 1:n for q = 1:n c1(p, q) = a(p, q)*b(p, q); end end time1(j)=toc; % 第二種方法:vectorized operation tic c2 = a.*b; time2(j)=toc; end subplot(2,1,1); plot(ns, time1, 'v-', ns, time2, '^-'); grid on legend('time1 for for-loop code', 'time2 for vectorized code', 'location', 'NorthWest'); xlabel('n'); ylabel('second'); subplot(2,1,2); plot(ns, time1./time2, '.-'); grid on xlabel('n'); ylabel('time1/time2'); fprintf('isequal(c1, c2) = %g\n', isequal(c1, c2));MATLAB version = 8.1.0.604 (R2013a) isequal(c1, c2) = 1

在上述範例中,我們計算兩個矩陣的「元素對元素」的乘積,在上圖中,第一個圖畫出迴圈運算和向量化運算所花的時間,第二個圖則是時間的比值。由第二個圖可以看出,加速的倍數大約都在 5 倍以上。

如果 for 迴圈的層數增加,計算速度會更慢,因此改用向量化運算的加速效能會更加顯著,下面這個範例,我們使用一般的矩陣乘法來進行驗證:

Example 2: 02-程式碼與記憶體之最佳化/matMultiply02.mfprintf('MATLAB version = %s\n', version); ns = 10*(1:20); for j=1:length(ns) n = ns(j); a = rand(n); b = rand(n); % 第一種方法:for-loop operation c1 = zeros(n); tic for p = 1:n for q = 1:n for r = 1:n c1(p, q) = c1(p, q)+a(p, r)*b(r, q); end end end time1(j)=toc; % 第二種方法:vectorized operation tic c2 = a*b; time2(j)=toc; end subplot(2,1,1); plot(ns, time1, 'v-', ns, time2, '^-'); grid on legend('time1 for for-loop code', 'time2 for vectorized code', 'location', 'NorthWest'); xlabel('n'); ylabel('second'); subplot(2,1,2); plot(ns, time1./time2, '.-'); grid on xlabel('n'); ylabel('time1/time2'); fprintf('isequal(c1, c2) = %g\n', isequal(c1, c2));MATLAB version = 8.1.0.604 (R2013a) isequal(c1, c2) = 1

由上述範例可看出,在不須要額外空間的前提下,使用向量化運算能夠將此範例(雙迴圈程式碼)有效地加速數十倍以上。

接著我們來看另外一類的向量化運算,這一類的範例通常須要額外的記憶體空間,已達到「以空間換取時間」的效果。但由於 JIT 加速氣的作用,「以空間換取時間」的效果可能不會很顯著,因此我們通常須要經過實際測試,才能得知在向量化後是否有實際加速功能。

在下面這個範例,我們要計算 n 項調和數列的總和,其中 n 是一個很大的整數,我們分別使用迴圈及向量化的方式來進行計算,並顯示加速倍數,如下:

Example 3: 02-程式碼與記憶體之最佳化/hsum01.mfprintf('MATLAB version = %s\n', version); n = 100000000; % === 第一種方法:for-loop operation tic total1 = 0; for i = 1:n total1 = total1+1/i; end time1 = toc; % === 第二種方法:vectorized operation tic total2 = sum(1./(1:n)); time2 = toc; fprintf('time1 = %g, time2 = %g, speedup factor = %g\n', time1, time2, time1/time2);MATLAB version = 8.1.0.604 (R2013a) time1 = 0.842769, time2 = 0.612797, speedup factor = 1.37528

由上述範例可知,使用向量化的運算,在執行速度的確有提升(這部分可能隨機器而有所不同),但並是非常不顯著。但是此種加速倍數是隨著 n 而變,我們可以畫出其關係曲線如下:

Example 4: 02-程式碼與記憶體之最佳化/hsum02.mfprintf('MATLAB version = %s\n', version); ns = 1000*(1:1000); for i=1:length(ns) n=ns(i); % 第一種方法:for-loop operation tic total1 = 0; for j = 1:n total1 = total1+1/j; end time1 = toc; % 第二種方法:vectorized operation tic total2 = sum(1./(1:n)); time2 = toc; % 計算倍數 speedupFactor(i)=time1/time2; end plot(ns, speedupFactor, '.-'); grid on xlabel('n'); ylabel('time1/time2');MATLAB version = 8.1.0.604 (R2013a)

由上述曲線可以看出,加速倍數並不是一條簡單遞增或遞減的曲線,不過只要 n 夠大,加速倍數的值都會大於一,代表向量化運算在這個範例的確會有加速的效果。事實上,如果您的程式是純粹的數值運算,而且沒有使用任何 for 迴圈或 while 迴圈,那麼其執行速度將會相當接近於純粹用 C 語言寫的程式碼。

但必須注意的是,上述向量化的運算式根據向量 sequence 來求和,但是 sequence 變數本身就要佔用記憶體空間,因此這種向量化運算的策略可以說是「以空間換取時間」,但是如果 n 太大,造成所需的空間超過機器的記憶體容量,此時機器就會將資料暫存到硬碟,造成整體運算速度的降低,而反而看不出來向量化運算的優點。

Hint
  • 若根據我的舊電腦(Pentium-450, 256 MB RAM)在 MATLAB 6.1(不支援 JIT)來進行測試,上述兩個範例的差異性會更加顯著,向量化程式碼的執行速度約是非向量化程式碼的 40 倍左右。
  • 在 MATLAB 6.5 版之後,採用了 JIT 編譯技術,因此「向量化」和「非向量化」的程式執行速度差異已經變小。

此外,為了使 JIT 加速器發揮最大效能,我們的程式碼也要盡量配合,以下是幾點注意事項:

  1. 迴圈中的索引變數盡量使用純量。
  2. 迴圈中的變數盡量是簡單的資料型態,並以不超過二維為主。
  3. 迴圈內所呼叫的函式僅限於 MATLAB 內建函式。

雖然向量化程式碼已經不再是提高 MATLAB 執行效率的唯一方法,在本節我們還是將針對向量化的運算來進行說明,主要著眼點在於:

  1. 不需要耗費額外記憶體的向量化程式碼,還是比非向量化的程式碼來得有效率。
  2. 若讀者還在使用 MATLAB 6.5 之前的版本,向量化的運算還是提高執行效率的主要關鍵。
  3. 碰到較複雜的資料結構(例如稀疏矩陣及維度超過 3 維的陣列),JIT 的效能並無法發揮,還是需要配合向量化的運算才可以讓程式碼更有效率。

一個 MATLAB 程式設計者的程度高低,可由其對「向量化之運算」的使用嫻熟程度來看出。要能夠熟練地運用向量化的運算,有下列三要素:

  1. 要對矩陣的索引(Indexing)非常熟悉,如此才能加以避免使用迴圈。
  2. 要對MATLAB可用的內建(Built-in)指令非常瞭解,如此才能盡量使用這些快速的指令。
  3. 要對問題本身認識清楚,才能將之轉化成使用內建指令即可解決的演算過程。

以下我們舉一些向量化運算的簡單(但卻常見)應用範例。例如,給定一向量 a 及一矩陣 x,若要將 a 的每一個元素乘上 x 的每一個直行,我們可使用三種方法:

  1. 以 for 迴圈來完成。
  2. 以內建的 diag 指令來達成此功能,此為向量化運算,但須要額外空間。
  3. 以內見的 bsxfun 指令來達成此功能,此為向量化運算,但不須要額外空間。

範例如下:

Example 6: 02-程式碼與記憶體之最佳化/colMultiply02.mfprintf('MATLAB version = %s\n', version); n = 10000; x = rand(3, n); y1 = zeros(size(x)); a = 1:n; % === 第 1 種方法:for-loop operation tic for i = 1:n y1(:,i)=x(:,i)*a(i); end time1 = toc; % === 第 2 種方法:vectorized operation which requires extra space tic y2 = x*diag(a); time2 = toc; % === 第 3 種方法:builtin command for vectorized operation tic y3 = bsxfun(@times, x, a); time3 = toc; fprintf('time1 = %g sec, time2 = %g sec, time3 = %g sec\n', time1, time2, time3); fprintf('time1/time2 = %g\n', time1/time2); fprintf('time1/time3 = %g\n', time1/time3);MATLAB version = 8.1.0.604 (R2013a) time1 = 0.00156713 sec, time2 = 0.310063 sec, time3 = 0.00122577 sec time1/time2 = 0.00505423 time1/time3 = 1.27849

由上述範例可以看出: 由此我們可以得到一個簡單的結論:若要提高運算速度,我們要盡量使用「不須要額外空間的內建指令」!

我們來看看幾個「不須要額外空間的內建指令」的加速效果。第一個指令是 cellfun,可用來對異質陣列的每一個元素進行運算。下面這個範例,我們讀入一個 CMU 發音辭典(包含約12萬個英文詞彙及其發音),然後計算每一個英語詞彙的長度:

Example 7: 02-程式碼與記憶體之最佳化/cellFun01.mfile='english.wpa'; [word, pa]=textread(file, '%s %s'); tic; len1=zeros(length(word),1); for i=1:length(word) len1(i)=length(word{i}); end time1=toc; tic; len2=cellfun(@length, word); time2=toc; fprintf('time1=%g sec, time2=%g sec, time1/time2=%g\n', time1, time2, time1/time2);time1=0.214358 sec, time2=0.0536109 sec, time1/time2=3.9984

由此可以看出,內建指令 cellfun 的確計算速度較快。(另一個類似的指令是 arrayfun,可以對一般陣列的每一個元素進行運算,請讀者自行參考線上說明,在此不再贅述。)

有關其他向量化的相關範例和問題,讀者們可以參考本章的習題和解答。

如果您所用的演算法無法向量化,花費的計算時間又很久,而且 MATLAB 的 JIT 也幫不上忙,那麼您就可以考慮下列方案:

相關細節請見本書「應用程式介面」之相關章節。
MATLAB程式設計:進階篇