matlab程式設(shè)計(jì)入門(mén)篇程式碼與記憶體之最佳化

1、MATLAB 程式設(shè)計(jì)入門(mén)篇程式碼與記憶體之最佳化,張智星jang@mirlab.orghttp://mirlab.org/jang清大資工系 多媒體檢索實(shí)驗(yàn)室,重點(diǎn)說(shuō)明6.5 版之前，MATLAB 執(zhí)行效率非常倚賴程式碼的向量化6.5 （含）版後，MATLAB 引進(jìn)了JIT 加速器，程式碼向量化的重要性相對(duì)降低。同學(xué)們?nèi)皂毑t解程式碼向量化的實(shí)作方法、優(yōu)缺點(diǎn)。本章範(fàn)例之執(zhí)行環(huán)境OS: Windows 7CPU: In

2、tel Core i7-2670QM, 2.20GHzRAM: 8GBMATLAB: 7.12.0.635 (R2011a),2-1 程式碼的向量化,6.5 版之前，MATLAB 程式環(huán)境是一個(gè)傳統(tǒng)的解譯器（Interpreter），在執(zhí)行 MATLAB的程式碼時(shí)，會(huì)進(jìn)行下列動(dòng)作:逐列對(duì)程式碼轉(zhuǎn)換為 p-code，這是 MATLAB 可以解讀的格式。對(duì)產(chǎn)生的 p-code 進(jìn)行逐列執(zhí)行。,MATLAB 程式運(yùn)作流程 (I),在執(zhí)

3、行每一列 p-code 時(shí)，都還包含一些「經(jīng)常開(kāi)銷(xiāo)」（Overhead）若該列執(zhí)行需要大量運(yùn)算，經(jīng)常開(kāi)銷(xiāo)就會(huì)顯得很渺小而不會(huì)拖累程式執(zhí)行速度。若該列只是簡(jiǎn)單的運(yùn)算，那經(jīng)常開(kāi)銷(xiāo)的比例就會(huì)相對(duì)提高。若該列是放在迴圈內(nèi)，那麼程式碼的執(zhí)行速度就會(huì)被這些經(jīng)常開(kāi)銷(xiāo)大幅拖慢。,MATLAB 程式運(yùn)作流程 (II),若要加快執(zhí)行速度MATLAB 6.5 版之前，就要盡量不使用迴圈 ，而盡量改用向量化的運(yùn)算，以降低經(jīng)常開(kāi)銷(xiāo)的比例。 MATLAB

4、 6.5 版（含）之後，已經(jīng)加入了 JIT-Accelerator，有效地降低了經(jīng)常開(kāi)銷(xiāo)的比例，同時(shí)也使得向量化程式碼的重要性越來(lái)越低。JIT: Just-In-Time對(duì)向量化運(yùn)算來(lái)進(jìn)行說(shuō)明，主要著眼點(diǎn)在於：很多舊的程式碼還是包含向量化的運(yùn)算。MATLAB 6.5 之前的版本，向量化的運(yùn)算還是提高執(zhí)行效率的主要關(guān)鍵。某些情況下，JIT-Accelerator 的效能無(wú)法完全發(fā)揮，還是需要配合向量化的運(yùn)算才可以讓程式碼更有效率

5、。,MATLAB 程式碼的向量化,傳統(tǒng)範(fàn)例 (I),計(jì)算 n 項(xiàng)調(diào)和數(shù)列的總和，若使用 for 迴圈，當(dāng) n = 100,000 時(shí)，可計(jì)算其執(zhí)行時(shí)間如下：範(fàn)例2-1：hsum01.m,ticn = 100000;total = 0;for i = 1:ntotal = total+1/i;endtoc,elapsed_time = 2.63000,由上可知，當(dāng) n = 100,000 時(shí)，使用 for 迴圈的程

6、式碼約需 2.63 秒才能執(zhí)行完畢。若改用向量化的運(yùn)算，可計(jì)時(shí)如下： 範(fàn)例2-2：hsum02.m,ticn = 100000;sequence = 1:n;total = sum(1./sequence);toc,elapsed_time = 0.0600,傳統(tǒng)範(fàn)例 (II),提示,上述的測(cè)試時(shí)間是根據(jù)Pentium-450, 256 MB RAM在 MATLAB 6.1 所得到的結(jié)果。如果你的電腦更快，計(jì)算時(shí)間的差

7、異性就不會(huì)那麼顯著。若使用我目前環(huán)境進(jìn)行測(cè)試：hsum01 ? 0.009773 sechsum02 ? 0.008938 sec在 MATLAB 6.5 版之後，採(cǎi)用了 JIT (Just-In-Time) 的編譯技術(shù)，因此「向量化」和「非向量化」的程式執(zhí)行速度差異較小。,要能夠熟練地運(yùn)用向量化的運(yùn)算，有下列三個(gè)訣竅：對(duì)矩陣的索引（Indexing）非常熟悉。對(duì) MATLAB 可用的內(nèi)建（Built-in）指令非常熟悉。

8、對(duì)問(wèn)題本身的解法非常瞭解。,程式碼向量化的要訣,若要對(duì)一個(gè)矩陣 x 的每一個(gè)直行乘上向量 a 的每一個(gè)元素，我們可使用內(nèi)建的 diag 指令來(lái)達(dá)成此功能，例如：範(fàn)例2-3： colMultiply01.m此為「以空間換取時(shí)間」的範(fàn)例。,x = [1 2 3; 1 1 1]; a = [3 2 1]; y = x*diag(a),y = 3 4 3 3 2 1

9、,矩陣直行乘法,若要對(duì)一個(gè)矩陣 x 的每一個(gè)橫列乘上向量 a 的每一個(gè)元素，也可以使用 diag 指令，例如： 範(fàn)例2-4： rowMultiply01.m此為「以空間換取時(shí)間」的範(fàn)例。,x = [1 2 3;4 5 6];a = [3 1];y = diag(a)*x,y = 3 6 9 4 5 6,矩陣橫列乘法,將一個(gè)行向量拷貝 n 份，可執(zhí)行如下：範(fàn)例2-5：

10、copyCol.ma*ones(1, n) 也可以產(chǎn)生相同的效果。,n = 10;a = [1 3 7]';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

11、 7 7 7 7 7 7 7 7 7 7,行向量拷貝,將一個(gè)列向量拷貝 n 份，可執(zhí)行如下：範(fàn)例2-6：copyRow.mones(n, 1)*a 也可以產(chǎn)生相同的效果。,n = 5;a = [1 4 8];rep = a(ones(1,n), :),rep = 1 4 8 1 4 8

12、 1 4 8 1 4 8 1 4 8,列向量拷貝,將一個(gè) m×n 矩陣 A 拷貝 M×N 次,並安排成維度大小為(M*m)×(N*n)的大矩陣，可輸入如下：範(fàn)例2-7：copyMatrix.m亦可使用內(nèi)建指令 repmat 來(lái)達(dá)成同樣的功能。,A = [1 2 3; 4 5 6];m = size(A, 1);n =

13、 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 1 2 3 1 2

14、 3 4 5 6 4 5 6 1 2 3 1 2 3 4 5 6 4 5 6,矩陣拷貝,2-2 JIT加速器,6.5 版後，MATLAB 引進(jìn)了 JIT (Just-In-Time) 的編譯技術(shù)，簡(jiǎn)稱(chēng) JIT 加速器（JIT-Accelerator）大幅地提高了迴圈的執(zhí)行效率。使

15、得一般程式碼的撰寫(xiě)，可以使用程式設(shè)計(jì)師最能代表演算法精神的方式來(lái)呈現(xiàn)。不必為了提高效率、勉強(qiáng)使用向量化運(yùn)算，造成程式碼維護(hù)的困難。,在JIT加速器的運(yùn)作下，向量化運(yùn)算的光環(huán)不如往日耀眼，測(cè)試範(fàn)例如下：範(fàn)例2-8：hsum03.m,JIT加速器可以加速迴圈運(yùn)算,ns = 1000*(1:100);for i=1:length(ns)% 第一種方法：for-loop operationtictotal = 0;for

16、j = 1:ns(i)total = total+sin(1/j);endtime1 = toc;,由此圖可以看出，time1/time2 的值都大約在 4 左右，顯示 JIT 加速器的效能已經(jīng)大幅提升。（若無(wú)JIT，此值約在40左右。）,% 第二種方法：vectorized operationtictotal = sum(sin(1./(1:n)));time2 = toc;% 計(jì)算倍數(shù)ratio(i

17、)=time1/time2;endplot(ns, ratio, '.-'); grid onxlabel('n'); ylabel('time1/time2');,2-2 JIT加速器,JIT 加速器的主要功能,JIT 加速器的主要功能，列出如下：JIT 加速器可將每一列 MATLAB 程式碼直接轉(zhuǎn)成「原機(jī)指令」，而不再直接使用 p-code，因此省去大部分 p-code 會(huì)遇到

18、的經(jīng)常開(kāi)銷(xiāo)。 JIT 加速器也會(huì)對(duì) Intel X86 為主的 Windows 及 Linux 進(jìn)行程式碼最佳化。 JIT 加速器可以分成兩部分： Just-In-Time Code Generation：將 p-code 轉(zhuǎn)成原機(jī)指令。Runtime Type Analysis：對(duì)資料型態(tài)的分析，以便用於 Just-In-Time Code Generation,有了 JIT 加速器，是否還要使用向量化運(yùn)算呢？基本的判斷原則如

19、下：所用的向量化運(yùn)算是根基於「以空間換取時(shí)間」，就必須經(jīng)過(guò)測(cè)試，才能決定是否值得使用。 所用的向量化運(yùn)算並非根基於「以空間換取時(shí)間」，那還是要盡量使用（尤其是要盡量使用內(nèi)建指令），因?yàn)檫@種向量化運(yùn)算不需要額外的空間來(lái)暫存資訊，因此效率會(huì)比 JIT 加速器更好。碰到較複雜的資料結(jié)構(gòu)（例如稀疏矩陣及維度超過(guò) 3 維的陣列），JIT 加速器並不會(huì)轉(zhuǎn)成原機(jī)指令，因此在這種情況下，也要盡量使用向量化運(yùn)算。,JIT 加速器與向量化運(yùn)算,使用「

20、不需要額外空間」的向量化運(yùn)算來(lái)提昇效率： 範(fàn)例2-9： matMultiply01.m,「不需要額外空間」的向量化運(yùn)算 (I),ns = 50*(1:20);for j=1:length(ns)n = ns(j);a = rand(n);b = rand(n);% 第一種方法：for-loop operationc = zeros(n);ticfor p = 1:nfor q = 1:nc(

21、p, q) = a(p, q)*b(p, q);end,endtime1(j)=toc;% 第二種方法：vectorized operationticc = a.*b;time2(j)=toc;endsubplot(2,1,1);plot(ns, time1, 'v-', ns, time2, '^-'); grid onlegend('time1 for for

22、-loop code', 'time2 for vectorized code');xlabel('n'); ylabel('second');subplot(2,1,2);plot(ns, time1./time2, '.-'); grid onxlabel('n'); ylabel('time1/time2');,「不需要

23、額外空間」的向量化運(yùn)算 (II),「不需要額外空間」的向量化運(yùn)算 (III),,第一個(gè)圖畫(huà)出迴圈運(yùn)算和向量化運(yùn)算所花的時(shí)間。第二個(gè)圖則是時(shí)間的比值，有此可以看出加速的倍數(shù)大約都在10～20之間。,為了使 JIT 加速器發(fā)揮最大效能，我們的程式碼也要盡量配合，以下是幾點(diǎn)注意事項(xiàng) 迴圈中的索引變數(shù)盡量使用純量。迴圈中的變數(shù)盡量是簡(jiǎn)單的資料型態(tài)，並以不超過(guò)二維為主。迴圈內(nèi)所呼叫的函式僅限於 MATLAB 內(nèi)建函式。,如何發(fā)揮JIT加速

24、器的最大功能,MATLAB 可隨時(shí)改變矩陣的維度，以容納新的資料，但是每一次增大矩陣的容量時(shí)，可能造成下列現(xiàn)象: MATLAB 必須要向電腦的作業(yè)系統(tǒng)索取記憶體，因此會(huì)造成程式執(zhí)行效率的降低 造成記憶體使用的分散現(xiàn)象，可能導(dǎo)致「實(shí)際雖仍有充足的記憶體，但卻未有連續(xù)空間以處理較大矩陣」的現(xiàn)象?；兑陨戏N種原因，若預(yù)先知道使用矩陣的維度大小，則最好實(shí)施矩陣的預(yù)先配置,2-3 矩陣的預(yù)先配置,計(jì)算亂數(shù)矩陣的各個(gè)元素平方範(fàn)例2-10：

25、preAllocate01.m,n = 10000000;x = rand(1, n);clear y;% 清除變數(shù) yticfor i = 1:ny(i) = x(i)^2;endtime1=toc;,矩陣預(yù)先配置的範(fàn)例 (I),y = zeros(1, n);% 預(yù)先配置所須矩陣 yticfor i = 1:ny(i) = x(i)^2;endtime2=toc;fprintf('tim

26、e1 = %g, time2 = %g, time1/time2 = %g\n', time1, time2, time1/time2);,矩陣預(yù)先配置的範(fàn)例 (II),time1 = 2.86088, time2 = 0.152491, time1/time2 = 18.761,上述範(fàn)例中，使用 y = zeros(1, n) 來(lái)事先宣告 y 所佔(zhàn)用的空間，因此程式碼的執(zhí)行效率快了將近 20 倍。,2-4 記憶體管理,MATL