【轉載】讓你的MATLAB代碼飛起來

原文地址:http://developer.51cto.com/art/201104/255128_all.htm數組

MATLAB語言是一種被稱爲是「演算紙」式的語言,所以追求的是方便性、靈活性以及交互性。在快速性上要比C語言這種性能強勁著稱的稍遜一籌。然而,經過一些手段,咱們也能讓MATLAB語言快起來,甚至和C差很少了!瀏覽器

 

MATLAB語言是一種被稱爲是「演算紙」式的語言,所以追求的是方便性、靈活性以及交互性。在快速性上要比C語言這種性能強勁著稱的稍遜一籌。然而,經過一些手段,咱們也能讓MATLAB語言快起來,甚至和C差很少了!函數

首先聲明:本文是一個初級教程,所以不少知識是假定你已經很熟悉了的;雖然我在討論讓代碼飛起來,但歷來不會說最快有多快,究竟有多快你要本身感受;做者水平不是很高,不免誤導你,當心甄別。工具


在正式討論以前,先看看這些好習慣你有沒有?性能

1. 使用 M-Lintspa

M-Lint是一個代碼分析檢查工具,它在你寫代碼的過程當中實時交互,發現你代碼的問題,按照最佳性能和最可維護性給出修改建議。命令行

注意:我可沒說是最正確!orm

若是沒有激活這個功能,依次使用File > Preferences > M-Lint,勾選Enable integrated M-Lint warning and error messages 。同時,還能夠設定你的偏好。視頻

激活後,在你寫代碼時就會實時交互了,錯誤的或者不推薦的部分會以紅色下劃線標出,鼠標通過紅色下劃線的語句或單詞,M-Lint給出提示信息。想一會兒看遍所有提示信息。使用Tools >M-Lint > (Save and) Show M-Lint Report2。htm

注:首次「觀看」先提示先保存一下。

2. 組織

給每個項目(project)創建一個單獨的文件夾。同屬於一個項目的文件保存在哪兒的都有,你找的時候就不費勁嗎!

寫頭部註釋,尤爲是H1。第一行就是H1。MATLAB中的內置函數的 help的內容其實就是讀取的這個函數的頭部註釋。怎麼寫,參照MATLAB內置函數。

將常常用到的控制檯命令存儲爲腳本(script)。若是有些命令反覆使用,仍是存爲腳本吧,沒別的意思,你要少敲多少次鍵盤啊!

3. 避免數據丟失

不要在腳本中使用 clear all。不幸的是這是一個你們經常使用的命令,有些書上還做爲一條規則確立起來,建議必須使用!要知道這個命令一執行,工做空間的數據可就不可逆轉的全沒了啊!

警告:注意呦! 

當心同名覆蓋。若是你在一個文件中,原本你的意思是兩個變量,你卻給他們起了相同的名字,那麼第一次的數據可就沒了。好比:

  1. result=max(a,b); %想求 a和 b之較大者  
  2. result=max(c,d); %想求 c 和d之較大者 

result結果是什麼?恐怕不是你想要的。不妨將其改成result1和 result2。相似的,也要當心文件重名的覆蓋,這個後果貌似更嚴重些。

下條內容請重視!

如何讓 MATLAB崩潰。

儘管 MATLAB是很穩定的,可是咱們仍然可讓它崩潰!使用第三方的MEX函數或者耗內存的操做好比視頻處理或者超大規模矩陣均可能會形成MATLAB崩潰。


 

若是你已經有這些好習慣,那麼恭喜,你要是還有其餘好習慣麻煩也告訴我一聲!若是沒有,相信你看完以後總該有了吧?好了,咱們開始!

1.使用profile

profile,Longman 給出的解釋是:a short description that gives important details about a person, a group of people, or a place。

MATLAB中內置了一個叫作profile的工具,來協助評估程序,也就是對程序運行過程的一個short description吧。主要命令有:

profile on 開啓

profile off 關閉

profile clear 清空數據

profile viewer 在profiler中看結果

下面咱們評估一下下面這個函數:

  1. function result =example1(count)   
  2. for k = 1:count   
  3.         result(k) = sin(k/50);   
  4.        
  5.         if result(k)<-0.9   
  6.            result(k) = gammaln(k);   
  7.         end   
  8. end  

爲了分析這個函數的效率,首先開啓並清空 profiler,而後運行這個函數,接下來看結果報告。即依次輸入:

  1. >> profile on, profile clear   
  2. >> example1(5000);   
  3. >> profile viewer 

這就是 profile 的基本語法。也有使用鼠標操做的方法,這裏就不介紹了,那樣雖然直觀單遠不及使用,命令方便。

因爲系統的不一樣,報告的結果通常是不同的。如下是個人系統得出的結果。

1.先看profile summary:

profile summary

2.點擊example1連接,進入具體各小項的評估。

(1)調用函數(children)、被調用函數(parents)。本例中都沒有。若是被 profile 的對象有調用函數或者被調用函數的話,會給出相應的數據。

(2)時間在哪些行被消耗(Lines where the most time was spent):

Lines where the most time was spent

從數據中咱們能夠看出哪些行消耗了多少時間(總時間和相對時間),被調用了多少次,以及直觀的柱形圖。

(3)另外一個有用的項目是 M-Lint 結果,給出了錯誤(警告、提示)所在的行,以及對應的建議修改信息,這些建議對代碼的改進是頗有價值的信息:

 M-Lint 結果

(4)最下面還有一個函數列表,是(2)的另外一種形式。看圖:

函數列表

最右側是函數代碼,前有行號、每一行調用的次數和小號的時間。消耗時間最多的行被標示了出來。最紅的消耗時間最多。

profiler工具的時間分辨率不是很高,所以,若是你的代碼運行的時間很短,有時候恐怕不能感知到。這時候不妨人爲的加入幾個循環,讓程序所運行幾回,而後進行分析。

必須指出,profile工具的做用主要是分析程序,得到程序運行的信息。若是想要知道程序運行的精確時間,使用計時器 tic/toc。以上面程序爲例,在命令行輸入:

  1. >> tic;example1(5000);toc 

輸出是:

  1. Elapsed time is 0.058522 seconds. 

爲了得到更爲精準的結果,你最好把瀏覽器、殺毒軟件、防火牆等等佔用CPU時間片的程序先關了,只剩下不能關掉的系統進程。

注意:profile在新版本中不斷被增強,可以使用的參數也愈來愈多,不過大多數根本用不着,若是你以爲那些參數頗有用,我相信你根本用不找看我這個小冊子了,要真是這樣,麻煩您不吝賜教,分享一些經驗。更詳細的內容,您仍是去看文檔去吧!

 

2. 預分配矩陣

MATLAB中的矩陣變量能夠動態增加行和列。好比:

  1. >>x=2   
  2. x=   
  3. 2   
  4. >>x(2,3)=1   
  5. x=   
  6. 2   0   0   
  7. 0   0   1  

看到沒?MATLAB自動調整了矩陣的大小!從內部實現上看,矩陣數據存儲單元被從新分配了更大的單元。若是矩陣的大小被反覆的調整(好比在循環中),從新分配存儲空間帶來的額外開銷會是很顯著的。爲了不反覆的矩陣存儲從新分配,預分配矩陣的存儲單元是一個不錯的選擇。一個推薦的方法是使用 zeros 函數命令。看下面的代碼:

  1. a(1) = 1;   
  2. b(1) = 0;   
  3. for k = 2:8000   
  4. a(k) = 0.99803 * a(k-1)-0.06279 * b(k-1);   
  5. b(k) = 0.06279 * a(k-1) + 0.99803 * b(k-1);   
  6. end   
  7. tic/toc計時運行獲得:   
  8. Elapsed time is 0.013306 seconds.  

簡單分析上面的代碼,知道,每一次 for,矩陣 a 和 b 的大小都要被從新分配,最終的大小事 8000 的列向量。若是咱們提早就給它們分配好大小爲 8000的存儲空間,看看結果怎麼樣:

  1. a=zeros(1,8000);    %預分配矩陣存儲單元   
  2. b=zeros(1,8000);   
  3. a(1) = 1;   
  4. b(1) = 0;   
  5. for k = 2:8000   
  6. a(k) = 0.99803 * a(k-1)-0.06279 * b(k-1);   
  7. b(k) = 0.06279 * a(k-1) + 0.99803 * b(k-1);   
  8. end   
  9. 及時運行獲得:   
  10. Elapsed time is 0.000753 seconds.  

看出來沒?速度提升了近 18 倍!像這種只需添加幾行代碼就能作到的狀況是不少的。這個例子也有特殊性,就是最後的結果大小已知,若是結果的大小可變、未知呢?不要緊,咱們能夠估計一下,最終結果最大能是多少?比估計到的最大再留出一些餘量就成了!若是你估計的仍是不夠大,那超出的部分還要反覆從新分配,不過這樣節省下來的時間也是很可觀的,畢竟能夠少分配不少次了! 最後呢,還要處理一下後事,好比你分配給變量 a 有 1000 個單元,但最終它只佔了300個,那你還要將那700個給收回來。看下面的代碼:

  1. a = zeros(1,10000);     %  預分配   
  2. count = 0;   
  3. for k = 1:10000   
  4. v = exp(rand*rand);   
  5. if v > 0.5        %  增加結果不肯定的來源   
  6. count = count + 1; a(count) = v;   
  7. end   
  8. end   
  9. a = a(1:count);    %調整矩陣大小   
  10. 未預分配時:Elapsed time is 0.052395 seconds.   
  11. 預分配後:Elapsed time is 0.008935 seconds. 

感慨:些微時間的意義在哪呢?背後是你對 MATLAB 的理解深度。哥玩的不是時間,是技術。

 

3. 向量化

不少狀況下,程序中的某些代碼能夠被向量化,向量化先後的速度每每在10 倍以上!向量化是最基本和最有效的讓代碼快起來的技巧,我都不肯意在後面叫「之一」了。

(1)向量化的計算

不少常規函數都是向量化的,它們做用於數組時,就好像是做用於數組中的每個元素。例如:

  1. >> sqrt([1,4,9,16])   
  2. ans =   
  3.       1     2     3     4   
  4. 考慮下面的函數:   
  5. function d = minDistance(x,y,z)    %尋找點集中距離遠點最近點   
  6. nPoints = length(x);   
  7. d = zeros(nPoints,1);    %  預分配   
  8. for k = 1:nPoints    %  計算每個點的距離   
  9. d(k) = sqrt(x(k)^2 + y(k)^2 + z(k)^2);   
  10. end   
  11. d = min(d);    %  獲得最小距離   
  12. 取  x=[1 2 3 4 5 6]; y=[2 3 5 2 1 4];z=[9 2 3 2 1 5];   
  13. 計時運行:Elapsed time is 0.008006 seconds. 

若是你寫出上面相似的代碼,說明你認真看了前面的內容。爲d預分配空間確實爲本例節省了很多時間。若是採用向量化計算,咱們能夠去掉for循環,直接計算向量。這裏要隆重推出「.」運算符,它表示的是對應元素進行運算。有.*和./和.\和.'和.^等。分別表示不帶.運算的對應元素運算。假設A是方陣,A^2是矩陣的 2 次乘冪,而 A.^2 表示矩陣 A 中的元素各自求平方組成新的矩陣。考慮下面的代碼:

  1. function d = minDistance(x,y,z)   
  2. d = sqrt(x.^2 + y.^2 + z.^2);    %  計算每一點的距離   
  3. d = min(d);     
  4. 計時運行:   
  5. Elapsed time is 0.005326 seconds. 

貌似差異不大?這就對了,別忘了,咱可就計算了6個值啊!這麼幾個值就有了這樣的差距,那x、y、z向量要是大一點,結果的差別就可想而知了!

更進一步的,咱們可使用d = sqrt(min(x.^2 + y.^2 + z.^2))取代後兩行語句,讓程序更加簡潔。

一下函數使用向量化的計算會更爲節省時間:min, max, repmat, meshgrid,sum, diff, prod等等。

(2)向量化邏輯

上面討論了計算的向量化,其實MATLAB的邏輯運算也是向量化的。好比:

  1. >> [1 4 2]>[2 3 1]   
  2. ans =   
  3.       0     1     1 

兩個數組「按元素」進行比較。向量的邏輯操做返回二進制的邏輯結果向量,即用0表明假,用1表明真。這爲何有用呢?由於MATLAB中有一些強勁的針對邏輯向量的函數。例如:

  1. >> find([1,5,3] < [2,2,4])     
  2. ans =   
  3. 1      3   
  4. >> any([1,5,3] < [2,2,4])     
  5. ans =   
  6. 1   
  7. >> all([1,5,3] < [2,2,4])     
  8. ans =   

其實,對通常向量(非邏輯向量)也是適用的!

  1. >> find(eye(4)==1)   
  2. ans =   
  3.       1   
  4.       6   
  5.      11   
  6.      16 

以上函數的用法請本身查閱函數說明。

 

4. 示例

(1)向量歸一標準化

將一個向量v歸一標準化,咱們但是使用v = v/norm(v),norm函數的做用是求模(範數)。

若是對一組向量 v(:,1), v(:,2),…進行歸一標準化,可使用一個循環計算v(:,k)/norm(v(:,k))。更好的策略是向量化計算:

  1. vMag = sqrt(sum(v.ˆ2));   
  2. v = v./vMag(ones(1,size(v,1)),:); 

(2)剔除元素

有時候,咱們須要將矩陣中的符合某些條件的元素剔除,固然可使用條件判斷加循環。咱們使用向量化剔除矩陣中的NaN和無窮兩類數:

  1. i = find(isnan(x) | isinf(x));   %在x中找到符合條件的數的位置   
  2. x(i) = [];    %剔除它   
  3. 或者,一樣的功能:   
  4. i = find(˜isnan(x) & ˜isinf(x));     %找到不符合的數   
  5. x = x(i);     %保留它   
  6. 進一步的,咱們能夠更加簡化,省略中間變量:   
  7. x(isnan(x) | isinf(x)) = [];     
  8. 以及   
  9. x = x(˜isnan(x) & ˜isinf(x)); 

(3)分段函數

信號分析中十分重要的 sinc(x)函數是分段的:x=0 時的值是 1,x!=0 時,sinc(x)=sin(x)/x。下面的代碼使用向量化方法處理分段:

  1. function y = sinc(x)   
  2. y = ones(size(x));          %  先設全部的y都是1   
  3. i = find(x ˜= 0);          %  找到非零x值   
  4. y(i) = sin(x(i)) ./ x(i);      %  計算非零值處的函數值   
  5. 更簡潔的,能夠寫成:   
  6. y = (sin(x) + (x == 0))./(x + (x == 0)) 

能看出來嗎?裏面用到了邏輯運算,實在是巧妙的很!

(4)其餘

還有些不經常使用的,算了,知道也八輩子用不着,珍惜腦細胞吧!

感慨:向量、矢量、相量、複數、數組、矩陣,這些名詞能分清楚麼?能分清楚知道內涵也就是爲何要這樣規定麼?不會也別問我!

 

5. 內嵌簡單函數

內嵌函數的意思就是將函數調用的函數的代碼直接寫到這個函數裏面來。因爲函數調用要作保護現場以及恢復現場等工做,也會額外增長一些時間消耗。若是調用的次數不是不少,這些時間是能夠忽略的,可是當調用次數不少的時候(好比500次),這個時間就很可觀了!

什麼樣的被調用函數適合內嵌呢?正如標題所說,是簡單的函數,特徵呢就是這個函數只有幾行代碼。若是這個函數很複雜,代碼很長,仍是死了這個心吧,內嵌是內嵌了,但是你看不懂代碼了,得不償失。程序的可讀性是很是重要的!

注意:必須是 M-File 實現的函數才能內嵌!

下面的代碼演示一個反覆調用median函數的內嵌方法。原代碼:

  1. y = zeros(size(x));      %  預分配   
  2. for k = 3:length(x)-2   
  3. y(k) = median(x(k-2:k+2));   
  4. end   
  5. 取  x=rand(1,2500);   
  6. 計時運行:Elapsed time is 0.030949 seconds. 

下面咱們試試內嵌。首先,要研究一下你要內嵌的函數,本例中就是median。在命令行中輸入:edit median,發現它是使用sort進行工做的。將核心代碼內嵌:

  1. y = zeros(size(x));     
  2. for k = 3:length(x)-2   
  3. tmp = sort(x(k-2:k+2));   
  4. y(k) = tmp(3); ;   
  5. end   
  6. 仍取x=rand(1,2500);   
  7. 計時運行:Elapsed time is 0.011379 seconds. 

以上就是一個演示,可見時間確實省去了很多。爲了確認你想內嵌的函數是不是用M-File實現的,你可使用「edit 函數名」命令試試看。

相關文章
相關標籤/搜索