MATLAB 程式設計入門篇 程式碼與記憶體之最佳化 張智星 jang@mirlab.org http://mirlab.org/jang 清大資工系 多媒體檢索實驗室

2-1 程式碼的向量化 重點說明 本章範例之執行環境 6.5 版之前，MATLAB 執行效率非常倚賴程式碼的向量化 6.5 （含）版後，MATLAB 引進了JIT 加速器，程式碼向量化的重要性相對降低。 同學們仍須瞭解程式碼向量化的實作方法、優缺點。 本章範例之執行環境 OS: Windows 7 CPU: Intel Core i7-2670QM, 2.20GHz RAM: 8GB MATLAB: 7.12.0.635 (R2011a)

MATLAB 程式運作流程 (I) 6.5 版之前，MATLAB 程式環境是一個傳統的解譯器（Interpreter），在執行 MATLAB的程式碼時，會進行下列動作: 逐列對程式碼轉換為 p-code，這是 MATLAB 可以解讀的格式。 對產生的 p-code 進行逐列執行。

MATLAB 程式運作流程 (II) 在執行每一列 p-code 時，都還包含一些「經常開銷」（Overhead） 若該列執行需要大量運算，經常開銷就會顯得很渺小而不會拖累程式執行速度。 若該列只是簡單的運算，那經常開銷的比例就會相對提高。 若該列是放在迴圈內，那麼程式碼的執行速度就會被這些經常開銷大幅拖慢。

MATLAB 程式碼的向量化 若要加快執行速度 對向量化運算來進行說明，主要著眼點在於： MATLAB 6.5 版（含）之後，已經加入了 JIT-Accelerator，有效地降低了經常開銷的比例，同時也使得向量化程式碼的重要性越來越低。 JIT: Just-In-Time 對向量化運算來進行說明，主要著眼點在於： 很多舊的程式碼還是包含向量化的運算。 MATLAB 6.5 之前的版本，向量化的運算還是提高執行效率的主要關鍵。 某些情況下，JIT-Accelerator 的效能無法完全發揮，還是需要配合向量化的運算才可以讓程式碼更有效率。

傳統範例 (I) 計算 n 項調和數列的總和，若使用 for 迴圈，當 n = 100,000 時，可計算其執行時間如下： 範例2-1：hsum01.m tic n = 100000; total = 0; for i = 1:n total = total+1/i; end toc elapsed_time = 2.63000

傳統範例 (II) 由上可知，當 n = 100,000 時，使用 for 迴圈的程式碼約需 2.63 秒才能執行完畢。若改用向量化的運算，可計時如下： 範例2-2：hsum02.m tic n = 100000; sequence = 1:n; total = sum(1./sequence); toc elapsed_time = 0.0600

提示 上述的測試時間是根據Pentium-450, 256 MB RAM在 MATLAB 6.1 所得到的結果。如果你的電腦更快，計算時間的差異性就不會那麼顯著。 若使用我目前環境進行測試： hsum01  0.009773 sec hsum02  0.008938 sec 在 MATLAB 6.5 版之後，採用了 JIT (Just-In-Time) 的編譯技術，因此「向量化」和「非向量化」的程式執行速度差異較小。

程式碼向量化的要訣 要能夠熟練地運用向量化的運算，有下列三個訣竅： 對矩陣的索引（Indexing）非常熟悉。 對 MATLAB 可用的內建（Built-in）指令非常熟悉。 對問題本身的解法非常瞭解。

矩陣直行乘法 若要對一個矩陣 x 的每一個直行乘上向量 a 的每一個元素，我們可使用內建的 diag 指令來達成此功能，例如： 範例2-3： colMultiply01.m 此為「以空間換取時間」的範例。 x = [1 2 3; 1 1 1]; a = [3 2 1]; y = x*diag(a) y = 3 4 3 3 2 1

矩陣橫列乘法 若要對一個矩陣 x 的每一個橫列乘上向量 a 的每一個元素，也可以使用 diag 指令，例如： 範例2-4： rowMultiply01.m 此為「以空間換取時間」的範例。 x = [1 2 3;4 5 6]; a = [3 1]; y = diag(a)*x y = 3 6 9 4 5 6

行向量拷貝 將一個行向量拷貝 n 份，可執行如下： 範例2-5：copyCol.m a*ones(1, n) 也可以產生相同的效果。 rep = a(:, ones(1,n)) rep = 1 1 1 1 1 1 1 1 1 1 3 3 3 3 3 3 3 3 3 3 7 7 7 7 7 7 7 7 7 7

列向量拷貝 將一個列向量拷貝 n 份，可執行如下： 範例2-6：copyRow.m ones(n, 1)*a 也可以產生相同的效果。 rep = a(ones(1,n), :) rep = 1 4 8

矩陣拷貝 將一個 m×n 矩陣 A 拷貝 M×N 次,並安排成維度大小為(M*m)×(N*n)的大矩陣，可輸入如下： 範例2-7：copyMatrix.m 亦可使用內建指令 repmat 來達成同樣的功能。 A = [1 2 3; 4 5 6]; m = size(A, 1); n = size(A, 2); M = 3; N = 2; mIndex = (1:m)'*ones(1,M); nIndex = (1:n)'*ones(1,N); B = A(mIndex, nIndex) B = 1 2 3 1 2 3 4 5 6 4 5 6

2-2 JIT加速器 6.5 版後，MATLAB 引進了 JIT (Just-In-Time) 的編譯技術，簡稱 JIT 加速器（JIT-Accelerator） 大幅地提高了迴圈的執行效率。 使得一般程式碼的撰寫，可以使用程式設計師最能代表演算法精神的方式來呈現。 不必為了提高效率、勉強使用向量化運算，造成程式碼維護的困難。

JIT加速器可以加速迴圈運算 在JIT加速器的運作下，向量化運算的光環不如往日耀眼，測試範例如下： 範例2-8：hsum03.m ns = 1000*(1:100); for i=1:length(ns) % 第一種方法：for-loop operation tic total = 0; for j = 1:ns(i) total = total+sin(1/j); end time1 = toc;

2-2 JIT加速器 % 第二種方法：vectorized operation tic total = sum(sin(1./(1:n))); time2 = toc; % 計算倍數 ratio(i)=time1/time2; end plot(ns, ratio, '.-'); grid on xlabel('n'); ylabel('time1/time2'); 由此圖可以看出，time1/time2 的值都大約在 4 左右，顯示 JIT 加速器的效能已經大幅提升。（若無JIT，此值約在40左右。）

JIT 加速器的主要功能 JIT 加速器的主要功能，列出如下： JIT 加速器可以分成兩部分： JIT 加速器可將每一列 MATLAB 程式碼直接轉成「原機指令」，而不再直接使用 p-code，因此省去大部分 p-code 會遇到的經常開銷。 JIT 加速器也會對 Intel X86 為主的 Windows 及 Linux 進行程式碼最佳化。 JIT 加速器可以分成兩部分： Just-In-Time Code Generation：將 p-code 轉成原機指令。 Runtime Type Analysis：對資料型態的分析，以便用於 Just-In-Time Code Generation

JIT 加速器與向量化運算 有了 JIT 加速器，是否還要使用向量化運算呢？基本的判斷原則如下： 所用的向量化運算是根基於「以空間換取時間」，就必須經過測試，才能決定是否值得使用。 所用的向量化運算並非根基於「以空間換取時間」，那還是要盡量使用（尤其是要盡量使用內建指令），因為這種向量化運算不需要額外的空間來暫存資訊，因此效率會比 JIT 加速器更好。 碰到較複雜的資料結構（例如稀疏矩陣及維度超過 3 維的陣列），JIT 加速器並不會轉成原機指令，因此在這種情況下，也要盡量使用向量化運算。

「不需要額外空間」的向量化運算 (I) 使用「不需要額外空間」的向量化運算來提昇效率： 範例2-9： matMultiply01.m ns = 50*(1:20); for j=1:length(ns) n = ns(j); a = rand(n); b = rand(n); % 第一種方法：for-loop operation c = zeros(n); tic for p = 1:n for q = 1:n c(p, q) = a(p, q)*b(p, q); end

「不需要額外空間」的向量化運算 (II) end time1(j)=toc; % 第二種方法：vectorized operation tic c = a.*b; time2(j)=toc; subplot(2,1,1); plot(ns, time1, 'v-', ns, time2, '^-'); grid on legend('time1 for for-loop code', 'time2 for vectorized code'); xlabel('n'); ylabel('second'); subplot(2,1,2); plot(ns, time1./time2, '.-'); grid on xlabel('n'); ylabel('time1/time2');

「不需要額外空間」的向量化運算 (III) 第一個圖畫出迴圈運算和向量化運算所花的時間。 第二個圖則是時間的比值，有此可以看出加速的倍數大約都在10～20之間。

如何發揮JIT加速器的最大功能 為了使 JIT 加速器發揮最大效能，我們的程式碼也要盡量配合，以下是幾點注意事項 迴圈中的索引變數盡量使用純量。 迴圈中的變數盡量是簡單的資料型態，並以不超過二維為主。 迴圈內所呼叫的函式僅限於 MATLAB 內建函式。

2-3 矩陣的預先配置 MATLAB 可隨時改變矩陣的維度，以容納新的資料，但是每一次增大矩陣的容量時，可能造成下列現象: 造成記憶體使用的分散現象，可能導致「實際雖仍有充足的記憶體，但卻未有連續空間以處理較大矩陣」的現象。 基於以上種種原因，若預先知道使用矩陣的維度大小，則最好實施矩陣的預先配置

矩陣預先配置的範例 (I) 計算亂數矩陣的各個元素平方 範例2-10： preAllocate01.m n = 10000000; x = rand(1, n); clear y; % 清除變數 y tic for i = 1:n y(i) = x(i)^2; end time1=toc;

矩陣預先配置的範例 (II) y = zeros(1, n); % 預先配置所須矩陣 y tic for i = 1:n y(i) = x(i)^2; end time2=toc; fprintf('time1 = %g, time2 = %g, time1/time2 = %g\n', time1, time2, time1/time2); time1 = 2.86088, time2 = 0.152491, time1/time2 = 18.761 上述範例中，使用 y = zeros(1, n) 來事先宣告 y 所佔用的空間，因此程式碼的執行效率快了將近 20 倍。

2-4 記憶體管理 MATLAB 提供 4 個和記憶體管理直接相關的指令： 指令 說明 clear 從記憶體中清除變數或函數 pack 重新安排記憶體的使用情況，使記憶體分散使用的情況降至最低 doc memory 顯示如何處理「Out of memory」的相關措施 quit 離開 MATLAB，將所有變數佔用的記憶體歸還給作業系統

「記憶體不足」之對策 執行 MATLAB 程式碼時，發生「Out of memory」的警告訊息，即表示已無足夠空間可容納新變數，可採取的解決方案為： 隨時使用 clear 指令清除不再須要的變數。 執行 pack 指令，減小記憶體分散使用的現象。 增加虛擬記憶體的大小。 選用適合的資料型態，例如 doubleuint8。 如果上述方法仍無法解決，那麼只有減少所用資料量的大小。