【转】提高MATLAB运行效率
用過Matlab的人都知道,Matlab是一種解釋性語言,存在計算速度慢的問題,為了提高程序的運行效率,matlab提供了多種實用工具及編碼技巧。
?
1. 循環矢量化
Matlab是為矢量和矩陣操作而設計的,因此,可以通過矢量化方法加速M文件的運行。矢量化是指將for循環和while循環轉換為等價的矢量或矩陣操作。下面給出一個循環的例子:
i=0;
for n = 0:0.1:1000
????i=i+1;
????y(i)=cos(n);
end
那么我們可以矢量化為:
n= 0:0.1:1000;
y=cos(n);
我們可以用tic和toc函數來查看上述各代碼運行的時間,采用for循環的程序0.39秒(具體時間和計算機配置有關),而矢量化后幾乎耗時為0。
2. 給數組或矩陣預分配內存
????特別是使用大型數組或矩陣時,Matlab進行動態內存分配和取消時,可能會產生內存碎片,這將導致大量閑置內存產生,預分配可通過提前給大型數據結構預約足夠空間來避免這個問題。
3. 用函數代替腳本文件
????因為每次調用MATLAB的腳本文件都需要將不必要的中間變量加載到內存中,每執行一次,就加載一次。函數在調用時被編譯成了偽代碼,只需要加載到內存一次。當多次調用同一個函數時會運行快一些。因此盡量多使用函數文件而少使用腳本文件,也是提高執行效率的一種方法。
4. 用Mex文件編寫循環代碼
????Matlab提供了與C和C++的接口,那么我們可以在用C或C++語言編寫耗時的循環代碼,然后通過接口程序在Matlab中轉換成dll文件,這就是我們所要的Mex文件,通過這種方法可以極大地提高計算速率。
?
1.?盡量避免使用循環結構
MATLAB變量的基本類型是矩陣,當對矩陣的每個元素循環處理時,運算速度很慢。因此編程時應盡量把數組和矩陣看作一個整體來進行編程,而不是像其他的程序設計語言那樣,使用循環結構對矩陣的元素循環進行處理。利用MATLAB提供的用于矢量化操作的函數,把循環矢量化,這樣既可以提高編程效率,也可以提高程序的執行效率。下面給出一個循環的例子:
i=0;
for n = 0:0.1:100
??????????i=i+1;
????????y(i)=cos(n)
end
上述程序段把數組中的每個元素都進行函數值計算,這樣會耗費大量的運算時間,我們可以把數組看作一個整體來處理,計算函數值,可以修改這個程序段如下。
n = 0:0.1:100;
y = cos(n)
通過使用MATLAB專門提供的測試程序運行時間的函數,可以發現,把數組看作一個整體,進行操作后,執行效率提高約300倍。
另外,在必須使用多重循環的情況下,建議在循環的外環執行循環次數少的,內環執行循環次數多的,這樣也可以顯著提高程序執行速度。
2.?在使用數組或矩陣之前先定義維數
MATLAB中的變量在使用之前不需要明確地定義和指定維數。但當未預定義數組或矩陣的維數時,當需賦值的元素下標超出現有的維數時,MATLAB?就為該數組或矩陣擴維一次,這樣就會大大降低程序的執行效率。因此,在使用數組或矩陣之前,預定義維數可以提高程序的執行效率。
3.?對矩陣元素使用下標或者索引操作
在MATLAB中,矩陣元素的引用可用兩個下標來表示。例如:A(i,j) 表示矩陣的第i行第j列的元素;A(1:k,j)表示矩陣A的第j列的前k個元素;A(:,j) 表示矩陣的第j列的所有元素。求矩陣A的第j列元素的平均值的表達式為mean(A(:,j))。
4.?盡量多使用函數文件少使用腳本文件
因為每次調用MATLAB的腳本文件都需要將不必要的中間變量加載到內存中,每執行一次,就加載一次。函數在調用時被編譯成了偽代碼,只需要加載到內存一次。當多次調用同一個函數時會運行快一些。因此盡量多使用函數文件而少使用腳本文件,也是提高執行效率的一種方法。
5.?在必須使用循環時,可以考慮轉換為C-MEX
當必須使用耗時的循環時,可以考慮將循環體中的語句轉換為C-MEX。C-MEX是將M文件通過MATLAB的編譯器轉換為可執行文件,是按照 MEX 技術要求的格式編寫相應的程序,通過編譯連接,生成擴展名為.dll的動態鏈接庫文件,可以在MATLAB環境下直接執行。這樣,循環體中的語句在執行時不必每次都解釋(interpret)。一般來說,C-MEX 文件的執行速度是相同功能的M文件執行速率的20~40倍。編寫C-MEX不同于M文件,需要了解MATLAB?C-MEX規范。幸運的是MATLAB提供了將M文件轉換為C-MEX的工具。
6.?內存優化
MATLAB在進行復雜的運算時需要占用大量的內存。合理使用內存和提高內存的使用效率,可以加快運行速度,減少系統資源的占用。
7.?內存管理函數和命令
(1)Clear variablename:從內存中刪除名稱為variablename的變量。
(2)Clear all:從內存中刪除所有的變量。
(3)Save:將指令的變量存入磁盤。
(4)Load:將save命令存入的變量載入內存。
(5)Quit:退出MATLAB,并釋放所有分配的內存。
(6)Pack:把內存中的變量存入磁盤,再用內存中的連續空間載回這些變量。考慮到執行效率問題,不能在循環中使用。
8.?節約內存的方法
(1)避免生成大的中間變量,并刪除不再需要的臨時變量。
(2)當使用大的矩陣變量時,預先指定維數并分配好內存,避免每次臨時擴充維數。
(3)當程序需要生成大量變量數據時,可以考慮定期將變量寫到磁盤,然后清除這些變量。
當需要這些變量時,再重新從磁盤加載。
(4)當矩陣中數據極少時,將全矩陣轉換為稀疏矩陣。
?
提高MATLAB程序效率的幾點原則,這些都是俺在這兩年中參加四次數模編寫大量m程序總結的經驗,加之網上很多英雄也是所見略同。
1.“計算向量、矩陣化,盡量減少for循環。”[/B]
因為MATLAB本來就是矩陣實驗室的意思,他提供了極其強大而靈活的矩陣運算能力,你就沒必要自己再用自己編寫的for循環去實現矩陣運算的功能了。另外由于matlab是一種解釋性語言,所以最忌諱直接使用循環語句。但在有些情況下,使用for循環可以提高程序的易讀性,在效率提高不是很明顯的情況下可以選擇使用for循環。
口說無憑,下面是利用tic與toc命令計算運算所用時間的方法,測試兩種編程的效率。需要說明的是沒有辦法精確計算程序執行時間,matlab幫助這樣寫到“Keep in mind that tic and toc measure overall elapsed time. Make sure that no other applications are running in the background on your system that could affect the timing of your MATLAB programs.”意思是說在程序執行的背后很可能有其他程序在執行,這里涉及到程序進程的的問題,m程序執行的過程中很可能有其他進程中斷m程序來利用cup,所以計算出來的時間就不僅是m程序的了,在這種情況下我把那些寄點去掉,進行多次計算求他的平均時間。
?
n = 100;
A(1:1000,1:1000) = 13;
C(1:1000,1) = 15;
D(1:1000,1) = 0;
for k = 1:n
????D(:) = 0;
????tic
????for i = 1:1000
????????for j = 1:1000
????????????D(i) = D(i) + A(i,j)*C(j);
????????end
????end
????t1(k) = toc;
????%------------------
????D(:) = 0;
????tic
????D = A*C;
????t2(k) = toc;
end
u = t1./t2;
u(u<0) = [];
plot(u)
p = mean(u)
t1、t2分別代表用for循環編程和矩陣化編程計算矩陣乘向量所用時間,u代表時間的比值。u(u<0) = [];是認為t1不可能小于t2,所以去掉不可能出現的情況。然后畫出圖形計算平均值。
經多次試驗結果大致相同,其中一次結果如下:
p =
????9.6196
------------t1時間是t2的十倍左右。
?
2.“循環內大數組預先定義--預先分配空間”[/U]
這一點原則是極其重要的,以至于在編寫m程序時編輯器會給出提示“'ver' might be growing inside a loop.Consider prealloacting for speed.”
clear
n = 50;
m = 1000;
for k = 1:n
????A = [];
????tic
????A(1:m,1:m) = 3;
????for i = 1:m
????????A(i,i) = i;
????end
????t1(k) = toc;
????%------------
????A = [];
????tic
????for j = 1:m
????????A(j,j) = j;
????end
????t2(k) = toc;
end
t2(t1>10^9) = [];
t1(t1>10^9) = [];
plot([t1;t2]')
t1、t2分別表示預先分配空間和循環中分配空間的時間,下圖上面一條t2、下面t1
3.“盡可能利用matlab內部提供的函數”[/U]
因為matlab內部提供的函數絕對是各種問題的最優算法,那寫程序都是他們大師級人物寫出來的,程序應該說相當高效,有現成的為什么不用那!???? 這個原則就不用實際的程序測試了。
?
?
關于MATLAB程序提速的問題,可以參考網上很多智者的文章,都比較經典。也可以看看我的上一篇文章,和網上大部分帖子有點不同,我是以實際的測試程序作為依據對如何提高MATLAB程序速度進行介紹的。??????這里我再補充幾點大家需要注意的。下面是我在國內一個比較出名的論壇看到的關于m程序提速的帖子,開始還真以為他們談論的都應該遵循。(盡信書不如無書)
帖子的一部分這樣說道:“當要預分配一個非double型變量時使用repmat函數以加速,如將以下代碼:
A = int8(zeros(100));
換成:
A = repmat(int8(0), 100, 100);”
凡事不能只憑自己的感覺,沒有一點實際的例子,對于權威我們還要有挑戰精神那,就不用說現在還不是經典的觀點了。下面是我寫的測試程序,我本來是想得到這位網友大哥的結果,但是實事不是我們想象的那么簡單。
n = 100;
m = 1000;
for k=1:n
????tic
????A = int8(ones(m));
????t1(k) = toc;
????tic
????B = repmat(int8(1),m,m);
????t2(k) = toc;
end
plot(1:n,t1,'r',1:n,t2)
isequal(A,B)
可以看出下面的紅線是t1,而且最后的一句返回1,說明兩種形式返回的值完全一樣。
由此我想說的是,不管是在我們做論文,還是寫博客的時候,別直接從網上或者別人文章那兒找點知識定理之類的補充自己那蒼白無力的文章。最好是自己動手編一下,“實踐是檢驗真理的唯一標準”。
?
經過這一測試,我感覺有必要,也有責任對這個論壇上的一大批經典談論加以測試。盡管這個結論是錯誤的但這還不足以證明論壇上的帖子都不是經典。
還有一點關于m程序提速的這樣說到:“在必須使用多重循環時下,如果兩個循環執行的次數不同,則在循環的外環執行循環次數少的,內環執行循環次數多的。這樣可以顯著提高速度。”
n=1000;
A = ones(1000)*13;
for k=1:n
????tic
????for i=1:10
????????for j=1:1000
????????????A(i,j)=A(i,j)*15;
????????end
????end
????t1(k)=toc;
????tic
????for i=1:1000
????????for j=1:10
????????????A(i,j)=A(i,j)*16;
????????end
????end
????t2(k)=toc;
end
t2(t1>10^9)=[];
t1(t1>10^9)=[];
t1(t2>10^9)=[];
t2(t2>10^9)=[];%去除外界因素影響所產生的寄點
plot(1:size(t1,2),t1,'r',1:size(t1,2),t2)
?
由這個時間圖可以看出for循環的嵌套順序對于速度是有影響的,雖然相對來說差別不是很大,但是畢竟論壇上的觀點是正確的。至于他所說的“顯著”二字就沒必要加上了。
此論壇還有一些提速的觀點列舉如下:
“遵守Performance Acceleration的規則
???? 關于什么是“Performance Acceleration”請參閱matlab的幫助文件。我只簡要的將
其規則總結如下7條:
1、只有使用以下數據類型,matlab才會對其加速:
logical,char,int8,uint8,int16,uint16,int32,uint32,double 而語句中如果使用了非以上的數據類型則不會加速,如:numeric,cell,structure,single,function handle,java classes,user classes,int64,uint64
2、matlab不會對超過三維的數組進行加速。
3、當使用for循環時,只有遵守以下規則才會被加速:a、for循環的范圍只用標量值來表示;
b、for循環內部的每一條語句都要滿足上面的兩條規則,即只使用支持加速的數據類型,只使用三維以下的數組;c、循環內只調用了內建函數(build-in function)。
4、當使用if、elseif、while和switch時,其條件測試語句中只使用了標量值時,將加速運行。
5、不要在一行中寫入多條操作,這樣會減慢運行速度。即不要有這樣的語句:
x =?a.name; for k=1:10000, sin(A(k)), end;
6、當某條操作改變了原來變量的數據類型或形狀(大小,維數)時將會減慢運行速度。
7、應該這樣使用復常量x = 7 + 2i,而不應該這樣使用:x = 7 + 2*i,后者會降低運行速度。”
“盡量用向量化的運算來代替循環操作。如將下面的程序:
i=0;
for t = 0:.01:10
????i = i+1;
????y(i) = sin(t);
end
替換為:
t = 0:.01:10;
y = sin(t);
速度將會大大加快。最常用的使用vectorizing技術的函數有:All、diff、ipermute、permute、reshape、ueeze、y、find、logical、prod、shiftdim、sub2ind、cumsum、ind2sub、ndgrid、repmat、sort、sum 等。”
“優先使用matlab內建函數,將耗時的循環編寫進MEX-File中以獲得加速。
b、使用Functions而不是Scripts 。”
“ 絕招:你也許覺得下面兩條是屁話,但有時候它真的是解決問題的最好方法。
1、改用更有效的算法
2、采用Mex技術,或者利用matlab提供的工具將程序轉化為C語言、Fortran語言。關于如何將M文件轉化為C語言程序運行,可以參閱本版帖子:“總結:m文件轉化為c/c++語言文件,VC編譯”。 ”
?
?
除了m程序提速的問題,這里還列出了《MATLAB代碼矢量化指南(譯)》
一、基本技術?
-----------------------------------------------------?
1)MATLAB索引或引用(MATLAB Indexing or Referencing)?
在MATLAB中有三種基本方法可以選取一個矩陣的子陣。它們分別是 下標法,線性法和邏輯法(subscripted,linear, and logical)。?
如果你已經熟悉這個內容,請跳過本節
1.1)下標法?
非常簡單,看幾個例子就好。?
A = 6:12;?
A([3,5])?
ans =?
8 10?
A([3:2:end])?
ans =?
8 10 12
A =?
[11 14 17;?
12 15 18;
13 16 19];?
A(2:3,2)?
ans =?
15?
16
1.2)線性法?
二維矩陣以列優先順序可以線性展開,可以通過現行展開后的元素序號?
來訪問元素。?
A =?
[11 14 17;
12 15 18;
13 16 19];?
A(6)?
ans =?
16
A([3,1,8])?
ans =?
13 11 18?
A([3;1;8])?
ans =?
13?
11?
18
1.3)邏輯法?
用一個和原矩陣具有相同尺寸的0-1矩陣,可以索引元素。在某個位置上為1表示選取元素,否則不選。得到的結果是一個向量。?
A = 6:10;?
A(logical([0 0 1 0 1]))?
ans =?
8 10?
A =?
[1 2?
3 4];?
B = [1 0 0 1];?
A(logical(B))?
ans =?
1 4?
-----------------------------------------------------?
2)數組操作和矩陣操作(Array Operations vs. Matrix Operations)?
對矩陣的元素一個一個孤立進行的操作稱作數組操作;而把矩陣視為一個整體進行的運算則成為矩陣操作。MATLAB運算符*,/,\,^都是矩陣運算,而相應的數組操作則是.*, ./, .\, .^?
A=[1 0 ;0 1];?
B=[0 1 ;1 0];?
A*B % 矩陣乘法?
ans =?
0 1?
1 0?
A.*B % A和B對應項相乘?
ans =?
0 0?
0 0?
------------------------------------------------------?
3)布朗數組操作(Boolean Array Operations)?
對矩陣的比較運算是數組操作,也就是說,是對每個元素孤立進行的。因此其結果就不是一個“真”或者“假”,而是一堆“真假”。這個結果就是布朗數組。?
D = [-0.2 1.0 1.5 3.0 -1.0 4.2 3.14];?
D >= 0?
ans =?
0 1 1 1 0 1 1?
如果想選出D中的正元素:?
D = D(D>0)?
D =?
1.0000 1.5000 3.0000 4.2000 3.1400?
除此之外,MATLAB運算中會出現NaN,Inf,-Inf。對它們的比較參見下例?
Inf==Inf返回真?
Inf<1返回假?
NaN==NaN返回假?
同時,可以用isinf,isnan判斷,用法可以顧名思義。在比較兩個矩陣大小時,矩陣必須具有相同的尺寸,否則會報錯。這是 你用的上size和isequal,isequalwithequalnans(R13及以后)。?
------------------------------------------------------?
4)從向量構建矩陣(Constructing Matrices from Vectors)?
在MATLAB中創建常數矩陣非常簡單,大家經常使用的是:?
A = ones(5,5)*10?
但你是否知道,這個乘法是不必要的??
A = 10;?
A = A(ones(5,5))?
A =?
10 10 10 10 10?
10 10 10 10 10?
10 10 10 10 10?
10 10 10 10 10?
10 10 10 10 10?
類似的例子還有:?
v = (1:5)';?
n = 3;?
M = v(:,ones(n,1))?
M =
1 1 1?
2 2 2?
3 3 3?
4 4 4?
5 5 5?
事實上,上述過程還有一種更加容易理解的實現方法:?
A = repmat(10,[5 5]);?
M = repmat([1:5]', [1,3]);?
其中repmat的含義是把一個矩陣重復平鋪,生成較大矩陣。更多詳細情況,參見函數repmat和meshgrid。?
-----------------------------------------------------?
5)相關函數列表(Utility Functions)?
ones 全1矩陣?
zeros 全0矩陣?
reshape 修改矩陣形狀?
repmat 矩陣平鋪?
meshgrid 3維plot需要用到的X-Y網格矩陣?
ndgrid n維plot需要用到的X-Y-Z...網格矩陣?
filter 一維數字濾波器,當數組元素前后相關時特別有用。?
cumsum 數組元素的逐步累計?
cumprod 數組元素的逐步累計?
eye 單位矩陣?
diag 生成對角矩陣或者求矩陣對角線?
spdiags 稀疏對角矩陣?
gallery 不同類型矩陣庫?
pascal Pascal 矩陣?
hankel Hankel 矩陣?
toeplitz Toeplitz 矩陣
==========================================================?
二、擴充的例子?
------------------------------------------------------?
6)作用于兩個向量的矩陣函數?
假設我們要計算兩個變量的函數F?
F(x,y) = x*exp(-x^2 - y^2)?
我們有一系列x值,保存在x向量中,同時我們還有一系列y值。我們要對向量x上的每個點和向量y上的每個點計算F值。換句話說,我們要計算對于給定向量x和y的所確定的網格上的F值。?
使用meshgrid,我們可以復制x和y來建立合適的輸入向量。然后?
可以使用第2節中的方法來計算這個函數。?
x = (-2:.2:2);?
y = (-1.5:.2:1.5)';?
[X,Y] = meshgrid(x, y);?
F = X .* exp(-X.^2 - Y.^2);?
如果函數F具有某些性質,你甚至可以不用meshgrid,比如?
F(x,y) = x*y ,則可以直接用向量外積?
x = (-2:2);?
y = (-1.5:.5:1.5);?
x'*y?
在用兩個向量建立矩陣時,在有些情況下,稀疏矩陣可以更加有?
效地利用存儲空間,并實現有效的算法。我們將在第8節中以一個?
實例來進行更詳細地討論.?
--------------------------------------------------------?
7)排序、設置和計數(Ordering, Setting, and Counting Operations)?
在迄今為止討論過的例子中,對向量中一個元素的計算都是獨立于同一向量的其他元素的。但是,在許多應用中,你要做的計算則可能與其它元素密切相關。例如,假設你用一個向量x來表示一 個集合。不觀察向量的其他元素,你并不知道某個元素是不是一個冗余元素,并應該被去掉。如何在不使用循環語句的情況下刪除冗余元素,至少在現在,并不是一個明顯可以解決的問題。?
解決這類問題需要相當的智巧。以下介紹一些可用的基本工具
max 最大元素?
min 最小元素?
sort 遞增排序?
unique 尋找集合中互異元素(去掉相同元素)?
diff 差分運算符[X(2) - X(1), X(3) - X(2), ... X(n) - X(n-1)]?
find 查找非零、非NaN元素的索引值?
union 集合并?
intersect 集合交?
setdiff 集合差?
setxor 集合異或
繼續我們的實例,消除向量中的多余元素。注意:一旦向量排序后,任何多余的元素就是相鄰的了。同時,在任何相等的相鄰元素在向量diff運算時變為零。這是我們能夠應用以下策略達到目的。我們現在在已排序向量中,選取那些差分非零的元素。
% 初次嘗試。不太正確!?
x = sort(x(:));?
difference = diff(x);?
y = x(difference~=0);
這離正確結果很近了,但是我們忘了diff函數返回向量的元素個數比輸入向量少1。在我們的初次嘗試中,沒有考慮到最后一個元素也可能是相異的。為了解決這個問題,我們可以在進行差分之前給向量x加入一個元素,并且使得它與以前的元素一定不同。一種實現的方法是增加一個NaN。
% 最終的版本。?
x = sort(x(:));?
difference = diff([x;NaN]);?
y = x(difference~=0);
我們使用(:)運算來保證x是一個向量。我們使用~=0運算,而不用find函數,因為find函數不返回NaN元素的索引值,而我們操作中差分的最后元素一定是NaN。這一實例還有另一種實現方式:
y=unique(x);
后者當然很簡單,但是前者作為一個練習并非無用,它是為了練習使用矢量化技術,并示范如何編寫你自己的高效代碼。此外,前者還有一個作用:Unique函數提供了一些超出我們要求的額外功能,這可能降低代碼的執行速度。?
??????假設我們不只是要返回集合x,而且要知道在原始的矩陣里每個相異元素出現了多少個“復本”。一旦我們對x排序并進行了差分,我們可以用find來確定差分變化的位置。再將這個變化位置進行差分,就可以得到復本的數目。這就是"diff of find of diff"的技巧。基于以上的討論,我們有:
% Find the redundancy in a vector x?
x = sort(x(:));?
difference = diff([x;max(x)+1]);?
count = diff(find([1;difference]));?
y = x(find(difference));?
plot(y,count)
這個圖畫出了x中每個相異元素出現的復本數。注意,在這里我們避開了NaN,因為find不返回NaN元素的索引值。但是,作為特例,NaN和Inf 的復本數可以容易地計算出來:
count_nans = sum(isnan(x(:)));?
count_infs = sum(isinf(x(:)));
另一個用于求和或者計數運算的矢量化技巧是用類似建立稀疏矩陣的方法實現的。這還將在第9節中作更加詳細的討論.?
-------------------------------------------------------?
8)稀疏矩陣結構(Sparse Matrix Structures)
在某些情況下,你可以使用稀疏矩陣來增加計算的效率。如果你構造一個大的中間矩陣,通常矢量化更加容易。在某些情況下,你可以充分利用稀疏矩陣結構來矢量化代碼,而對于這個中間矩陣不需要大的存儲空間。?
假設在上一個例子中,你事先知道集合y的域是整數的子集,?
{k+1,k+2,...k+n};即,?
y = (1:n) + k
例如,這樣的數據可能代表一個調色板的索引值。然后,你就可以對集合中每個元素的出現進行計數(構建色彩直方圖?譯者)。這是對上一節中"diff of find of diff"技巧的一種變形。現在讓我們來構造一個大的m x n矩陣A,這里m是原始x向量中的元素數, n是集合y中的元素數。?
A(i,j) = 1 if x(i) = y(j)?
0 otherwise
回想一下第3節和第4節,你可能認為我們需要從x和y來構造矩陣A。如果當然可以,但要消耗許多存儲空間。我們可以做得更好,因為我們知道,矩陣A中的多數元素為0,x中的每個元素對應的行上只有一個值為1。?
以下就是構造矩陣的方法(注意到y(j) = k+j,根據以上的公式):?
x = sort(x(:));?
A = sparse(1:length(x), x+k, 1, length(x), n);
現在我們對A的列進行求和,得到出現次數。?
count = sum(A);?
在這種情況下,我們不必明確地形成排序向量y,因為我們事先知道?
y = 1:n + k.
這里的關鍵是使用數據,(也就是說,用x控制矩陣A的結構)。由于x在一個已知范圍內取整數值,我們可以更加有效地構造矩陣。 假設你要給一個很大矩陣的每一列乘以相同的向量。使用稀疏矩陣,不僅可以節省空間,并且要比在第5節介紹的方法更加快速. 下面是它的工作方式:?
F = rand(1024,1024);?
x = rand(1024,1);?
% 對F的所有行進行點型乘法.?
Y = F * diag(sparse(x));?
% 對F的所有列進行點型乘法.?
Y = diag(sparse(x)) * F;
我們充分利用矩陣乘法算符來執行大規模運算,并使用稀疏矩陣以防止臨時變量變得太大。?
--------------------------------------------------------?
9)附加的例子(Additional Examples)?
下面的例子使用一些在本技術手冊中討論過的方法,以及其它一些相關方法。請嘗試使用tic 和toc(或t=cputime和cputime-t),看一下速度加快的效果。?
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>?
用于計算數組的雙重for循環。?
使用的工具:數組乘法?
優化前:?
A = magic(100);?
B = pascal(100);?
for j = 1:100?
?? for k = 1:100;?
?????? X(j,k) = sqrt(A(j,k)) * (B(j,k) - 1);?
?? end?
end?
優化后:?
A = magic(100);?
B = pascal(100);?
X = sqrt(A).*(B-1);?
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>?
用一個循環建立一個向量,其元素依賴于前一個元素使用的工具:FILTER, CUMSUM, CUMPROD?
優化前:?
A = 1;?
L = 1000;?
for i = 1:L?
?? A(i+1) = 2*A(i)+1;?
end?
優化后:?
L = 1000;?
A = filter([1],[1 -2],ones(1,L+1));?
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>?
如果你的向量構造只使用加法或乘法,你可使用cumsum或cumprod函數。?
優化前:?
n=10000;?
V_B=100*ones(1,n);?
V_B2=100*ones(1,n);?
ScaleFactor=rand(1,n-1);?
for i = 2:n?
????V_B(i) = V_B(i-1)*(1+ScaleFactor(i-1));?
end?
for i=2:n?
????V_B2(i) = V_B2(i-1)+3;?
end?
優化后:?
n=10000;?
V_A=100*ones(1,n);?
V_A2 = 100*ones(1,n);?
ScaleFactor=rand(1,n-1);?
V_A=cumprod([100 1+ScaleFactor]);?
V_A2=cumsum([100 3*ones(1,n-1)]);?
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>?
向量累加,每5個元素進行一次:?
工具:CUMSUM , 向量索引?
優化前:?
% Use an arbitrary vector, x?
x = 1:10000;?
y = [];?
for n = 5:5:length(x)?
????y = [y sum(x(1:n))];?
end?
優化后(使用預分配):?
x = 1:10000;?
ylength = (length(x) - mod(length(x),5))/5;?
% Avoid using ZEROS command during preallocation?
y(1:ylength) = 0;?
for n = 5:5:length(x)?
????y(n/5) = sum(x(1:n));?
end?
優化后(使用矢量化,不再需要預分配):?
x = 1:10000;?
cums = cumsum(x);?
y = cums(5:5:length(x));?
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>?
操作一個向量,當某個元素的后繼元素為0時,重復這個元素:?
工具:FIND, CUMSUM, DIFF?
任務:我們要操作一個由非零數值和零組成的向量,要求把零替換成為它前面的非零數值。例如,我們要轉換下面的向量:?
a=2; b=3; c=5; d=15; e=11;?
x = [a 0 0 0 b 0 0 c 0 0 0 0 d 0 e 0 0 0 0 0];?
為:?
x = [a a a a b b b c c c c c d d e e e e e e];?
解(diff和cumsum是反函數):?
valind = find(x);?
x(valind(2:end)) = diff(x(valind));?
x = cumsum(x);?
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>?
將向量的元素累加到特定位置上?
工具:SPARSE?
優化前:?
% The values we are summing at designated indices?
values = [20 15 45 50 75 10 15 15 35 40 10];?
% The indices associated with the values are summed cumulatively?
indices = [2 4 4 1 3 4 2 1 3 3 1];?
totals = zeros(max(indices),1);?
for n = 1:length(indices)?
????totals(indices(n)) = totals(indices(n)) + values(n);?
end?
優化后:?
indices = [2 4 4 1 3 4 2 1 3 3 1];?
totals = full(sparse(indices,1,values));
注意:這一方法開辟了稀疏矩陣的新用途。在使用sparse命令創建稀疏矩陣時,它是對分配到同一個索引的所有值求和,而不是替代已有的數值。這稱為"向量累加",是MATLAB處理稀疏矩陣的方式。
?
關于MATLAB的效率問題,很多文章,包括我之前寫的一些,主要集中在使用向量化以及相關的問題上。但是,最近我在實驗時對代碼進行profile的過程中,發現在新版本的MATLAB下,for-loop已經得到了極大優化,而效率的瓶頸更多是在函數調用和索引訪問的過程中。
由于MATLAB特有的解釋過程,不同方式的函數調用和元素索引,其效率差別巨大。不恰當的使用方式可能在本來不起眼的地方帶來嚴重的開銷,甚至可能使你的代碼的運行時間增加上千倍(這就不是多買幾臺服務器或者增加計算節點能解決的了,呵呵)。
下面通過一些簡單例子說明問題。(實驗選在裝有Windows Vista的一臺普通的臺式機(Core2 Duo 2.33GHz + 4GB Ram)進行,相比于計算集群, 這可能和大部分朋友的環境更相似一些。實驗過程是對某一個過程實施多次的整體進行計時,然后得到每次過程的平均時間,以減少計時誤差帶來的影響。多次實驗,在均值附近正負20%的范圍內的置信度高于95%。為了避免算上首次運行時花在預編譯上的時間,在開始計時前都進行充分的“熱身”運行。)?
函數調用的效率
一個非常簡單的例子,把向量中每個元素加1。(當然這個例子根本不需要調函數,但是,用它主要是為了減少函數執行本身的時間,突出函數解析和調用的過程。)
作為baseline,先看看最直接的實現
% input u: u is a 1000000 x 1 vectorv = u + 1;?
這個過程平均需要0.00105 sec。
而使用長期被要求盡量避免的for-loop?
n = numel(u);% v = zeros(n, 1) has been pre-allocated.for i = 1 : n????v(i) = u(i) + 1;end
所需的平均時間大概是0.00110 sec。從統計意義上說,和vectorization已經沒有顯著差別。無論是for-loop或者vectorization,每秒平均進行約10億次“索引-讀取-加法-寫入”的過程,計算資源應該得到了比較充分的利用。
要是這個過程使用了函數調用呢?MATLAB里面支持很多種函數調用方式,主要的有m-function, function handle, anonymous function, inline, 和feval,而feval的主參數可以是字符串名字,function handle, anonymous function或者inline。
用m-function,就是專門定義一個函數
function y = fm(x)????y = x + 1;
在調用時?
for i = 1 : n????v(i) = fm(u(i));end
function handle就是用@來把一個function賦到一個變量上,類似于C/C++的函數指針,或者C#里面的delegate的作用
fh = @fm;for i = 1 : n????v(i) = fh(u(i));end
anonymous function是一種便捷的語法來構造簡單的函數,類似于LISP, Python的lambda表達式
fa = @(x) x + 1;for i = 1 : n????v(i) = fa(u(i));end
inline function是一種傳統的通過表達式字符串構造函數的過程
fi = inline('x + 1', 'x');for i = 1 : n????v(i) = fi(u(i));end
feval的好處在于可以以字符串方式指定名字來調用函數,當然它也可以接受別的參數。
v(i) = feval('fm', u(i));v(i) = feval(fh, u(i));v(i) = feval(fa, u(i));?
對于100萬次調用(包含for-loop本身的開銷,函數解析(resolution),壓棧,執行加法,退棧,把返回值賦給接收變量),不同的方式的時間差別很大: m-function0.385 secfunction handle0.615 secanonymous function0.635 secinline function166.00 secfeval('fm', u(i))8.328 secfeval(fh, u(i))0.618 secfeval(fa, u(i))0.652 secfeval(@fm, u(i))2.788 secfeval(@fa, u(i))34.689 sec
從這里面,我們可以看到幾個有意思的現象:?
- 首先,調用自定義函數的開銷遠遠高于一個簡單的計算過程。直接寫 u(i) = v(i) + 1 只需要??0.0011 秒左右,而寫u(i) = fm(v(i)) 則需要0.385秒,時間多了幾百倍,而它們干的是同一件事情。這說明了,函數調用的開銷遠遠大于for-loop自己的開銷和簡單計算過程——在不同情況可能有點差別,一般而言,一次普通的函數調用花費的時間相當于進行了幾百次到一兩千次雙精度浮點加法。
- 使用function handle和anonymous function的開銷比使用普通的m-函數要高一些。這可能是由于涉及到句柄解析的時間,而普通的函數在第一次運行前已經在內存中進行預編譯。
- inline function的運行時間則令人難以接受了,竟然需要一百多秒(是普通函數調用的四百多倍,是直接計算的十幾萬倍)。這是因為matlab是在每次運行時臨時對字符串表達式(或者它的某種不太高效的內部表達)進行parse。
- feval(fh, u(i))和fh(u(i)),feval(fa, u(i))和fa(u(i))的運行時間很接近,表明feval在接受句柄為主參數時本身開銷很小。但是,surprising的是 for i = 1 : n????v(i) = feval(@fm, u(i));end比起fh = @fm;for i = 1 : n????v(i) = feval(fh, u(i));end慢了4.5倍 (前者0.618秒,后者2.788秒)。
for i = 1 : n????v(i) = feval(@(x) x + 1, u(i));end比起fa = @(x) x + 1;for i = 1 : n????v(i) = feval(fa, u(i));end竟然慢了53倍(前者0.652秒,后者34.689秒)。
由于在MATLAB的內部實現中,function handle的解析是在賦值過程中進行的,所以預先用一個變量把句柄接下,其效果就是預先完成了句柄解析,而如果直接把@fm或者@(x) x + 1寫在參數列上,雖然寫法簡潔一些,但是解析過程是把參數被賦值到所調函數的局部變量時才進行,每調用一次解析一次,造成了巨大的開銷。 - feval使用字符串作為函數名字時,所耗時間比傳入句柄大,因為這涉及到對函數進行搜索的時間(當然這個搜索是在一個索引好的cache里面進行(除了第一次),而不是在所有path指定的路徑中搜索。)
在2007年以后,MATLAB推出了arrayfun函數,上面的for-loop可以寫為 v = arrayfun(fh, u)
這平均需要4.48 sec,這比起for-loop(需時0.615 sec)還慢了7倍多。這個看上去“消除了for-loop"的函數,由于其內部設計的原因,未必能帶來效率上的正效果。?
元素和域的訪問
除了函數調用,數據的訪問方式對于效率也有很大影響。MATLAB主要支持下面一些形式的訪問:
- array-index??A(i):??
- cell-index:??C{i};??
- struct field:??S.fieldname
- struct field (dynamic):??S.('fieldname')這里主要探索單個元素或者域的訪問(當然,MATLAB也支持對于子數組的非常靈活整體索引)。
對于一百萬次訪問的平均時間
A(i) for a numeric array0.0052 secC{i} for a cell array0.2568 secstruct field0.0045 secstruct field (with dynamic name)1.0394 sec
我們可以看到MATLAB對于單個數組元素或者靜態的struct field的訪問,可以達到不錯的速度,在主流臺式機約每秒2億次(連同for-loop的開銷)。而cell array的訪問則明顯緩慢,約每秒400萬次(慢了50倍)。MATLAB還支持靈活的使用字符串來指定要訪問域的語法(動態名字),但是,是以巨大的開銷為代價的,比起靜態的訪問慢了200倍以上。
關于Object-oriented Programming
MATLAB在新的版本中(尤其是2008版),對于面向對象的編程提供了強大的支持。在2008a中,它對于OO的支持已經不亞于python等的高級腳本語言。但是,我在實驗中看到,雖然在語法上提供了全面的支持,但是matlab里面面向對象的效率很低,開銷巨大。這里僅舉幾個例子。
- object中的property的訪問速度是3500萬次,比struct field慢了6-8倍。MATLAB提供了一種叫做dependent property的屬性,非常靈活,但是,效率更低,平均每秒訪問速度竟然低至2.6萬次(這種速度基本甚至難以用于中小規模的應用中)。
- object中method調用的效率也明顯低于普通函數的調用,對于instance method,每百萬次調用,平均需時5.8秒,而對于static method,每百萬次調用需時25.8秒。這相當于每次調用都需要臨時解析的速度,而matlab的類方法解析的效率目前還明顯偏低。
- MATLAB中可以通過改寫subsref和subsasgn的方法,對于對象的索引和域的訪問進行非常靈活的改造,可以通過它創建類似于數組的對象。但是,一個符合要求的subsref的行為比較復雜。在一個提供了subsref的對象中,大部分行為都需要subsref進行調度,而默認的較優的調度方式將被掩蓋。在一個提供了subsref的類中(使用一種最快捷的實現),object property的平均訪問速度竟然降到1萬5千次每秒。建議[/U]
根據上面的分析,對于撰寫高效MATLAB代碼,我有下面一些建議:?
- 雖然for-loop的速度有了很大改善,vectorization(向量化)仍舊是改善效率的重要途徑,尤其是在能把運算改寫成矩陣乘法的情況下,改善尤為顯著。
- 在不少情況下,for-loop本身已經不構成太大問題,尤其是當循環體本身需要較多的計算的時候。這個時候,改善概率的關鍵在于改善循環體本身而不是去掉for-loop。
- MATLAB的函數調用過程(非built-in function)有顯著開銷,因此,在效率要求較高的代碼中,應該盡可能采用扁平的調用結構,也就是在保持代碼清晰和可維護的情況下,盡量直接寫表達式和利用built-in function,避免不必要的自定義函數調用過程。在次數很多的循環體內(包括在cellfun, arrayfun等實際上蘊含循環的函數)形成長調用鏈,會帶來很大的開銷。
- 在調用函數時,首選built-in function,然后是普通的m-file函數,然后才是function handle或者anonymous function。在使用function handle或者anonymous function作為參數傳遞時,如果該函數被調用多次,最好先用一個變量接住,再傳入該變量。這樣,可以有效避免重復的解析過程。
- 在可能的情況下,使用numeric array或者struct array,它們的效率大幅度高于cell array(幾十倍甚至更多)。對于struct,盡可能使用普通的域(字段,field)訪問方式,在非效率關鍵,執行次數較少,而靈活性要求較高的代碼中,可以考慮使用動態名稱的域訪問。
- 雖然object-oriented從軟件工程的角度更為優勝,而且object的使用很多時候很方便,但是MATLAB目前對于OO的實現效率很低,在效率關鍵的代碼中應該慎用objects。
- 如果需要設計類,應該盡可能采用普通的property,而避免靈活但是效率很低的dependent property。如非確實必要,避免重載subsref和subsasgn函數,因為這會全面接管對于object的接口調用,往往會帶來非常巨大的開銷(成千上萬倍的減慢),甚至使得本來幾乎不是問題的代碼成為性能瓶頸。
總結
以上是生活随笔為你收集整理的【转】提高MATLAB运行效率的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: C++11详解
- 下一篇: 为什么你在互联网上搞不到钱?