MatLab-向量化技巧

1、MatLab 向量化技巧家家 2007-05-17 13:53大家都知道MatLab是一種解釋性語言，它的長處在矩陣運算。因此需要將問題盡量用矩陣表示，并且要避免對單個矩陣元素操作，而用整體矩陣運算。一個重要過程就是所謂將問題矢量化(vectorizing)。一個問題是如果需要對矩陣的不同元素作不同操作，那么怎么辦？舉個例子對矩陣X，作A=5*sin(X)/X，顯然當(dāng)X某元素-0時，A對應(yīng)元素-5。按照通常在其他語言思路編程為：Nx=2000;X=2.0*pi*(rand(Nx,Nx)-0.5);X(1:Nx/2,:)=0.0;%生成NxXNx矩陣X，一半為隨機分布，一半為0。A=zeros(

2、Nx,Nx);form=1:Nxforn=1:Nxif(abs(X(m,n)=eps)A(m,n)=5.0.*sin(X(m,n)./X(m,n);elseA(m,n)=5.0;endendend%(ProgramNo.1)以上代碼通過循環(huán)語句對每一個矩陣元素進行單獨操作，這段代碼在MatLab里效率特別低。為避免循環(huán)，考慮如下方法：方法I：用所謂logicalsubscripting。Nx=2000;X=2.0*pi*(rand(Nx,Nx)-0.5);X(1:Nx/2,:)=0.0;A=zeros(Nx,Nx);I=(abs(X)=eps);A=I.*5.0.*sin(X)./(X+(1-

3、I)+(1-I).*5.0;%(ProgramNo.2)或者上面最后一行用如下代碼代替：A(I)=5.0.*sin(X(I)./X(I);A(I)=5.0;%(ProgramNo.3)一個看起來更有效的辦法是先生成一個所有元素都是5的A矩陣。代碼如下：Nx=2000;X=2.0*pi*(rand(Nx,Nx)-0.5);X(1:Nx/2,:)=0.0;A=5.0;A=A(ones(Nx,Nx);%用Tonystrick生成矩陣A，所有元素=5I=(abs(X)=eps);A(I)=5.0.*sin(X(I)./X(I);%(ProgramNo.4)方法II：用find函數(shù)。Nx=2000;X=

4、2.0*pi*(rand(Nx,Nx)-0.5);X(1:Nx/2,:)=0.0;A=5.0;A=A(ones(Nx,Nx);I=find(abs(X)=eps);A(I)=5.0.*sin(X(I)./X(I);%(ProgramNo.5)我在P41.6G，512MDDR266上作了個測試，每個程序執(zhí)行5遍，選最小執(zhí)行時間：ProgramNo.time(s)171.2023.3132.6442.8053.22顯然ProgramNo.3是最快的。如果X矩陣中0元素很少，那么ProgramNo.2、ProgramNo.3以及ProgramNo.4的執(zhí)行時間會非常接近。大家可以用profiler分

5、析每一個語句所用時間情況。以上是我用過的兩種方法，有沒有更有效的方法呢？還有個問題就是，如何將主對角元素的Index加入到變量I中？例如需要將滿足abs(X)1)if (mod(n,2)=1) n=n*3+1;elseif(mod(n,2)=0) n=n/2;else break;endendf=n;驗證1:, 以前我們可能這樣做：代碼:a=zeros(1,);for k=1:a(k)=jiaogu(k);endall(a)現(xiàn)在我們只需這樣：代碼:all(arrayfun(jiaogu,1:)舉這兩個例子目的是為了讓大家知道很多以前不能優(yōu)化的程序利用這個函數(shù)還可以進一步優(yōu)化。具體了解請doc

6、arrayfun 。另外cellfun和structfun也提供了類似的功能，關(guān)于其用法，幫助文檔里寫的很詳細(xì)。向量化操作的又一重要函數(shù)accumarray的用法總結(jié)向量化編程是MATLAB編程區(qū)別于其他語言的最重要特征之一，MATLAB不斷增強其向量化編程的能力。版本7以后出現(xiàn)的accumarray函數(shù)就是一個很好的例子。accumarray函數(shù)最早出現(xiàn)于7.0版（R14），在隨后的7.1 (R14SP3)， 7.2 (R2006a)版里又對其功能進行增強。大家可以在幫助文檔里了解其基本信息，這里結(jié)合自己的使用經(jīng)驗給出其用法的舉例，大家可以看看其效果。例1：生成一個*1的向量，其元素服從1，

7、之間離散均勻分布，找出都有哪些元素出現(xiàn)了a=unidrnd(,1);方法1：tic;b=union(a,a);toc;Elapsed time is 0. seconds.方法2：tic;c=unique(a);toc;Elapsed time is 0. seconds.方法3tic,d=accumarray(a,1,1);e=find(d);tocElapsed time is 0. seconds.d 給出了每個元素出現(xiàn)的個數(shù)，如果還需要這方面的信息，無疑方法3占絕對優(yōu)勢。例2：在1000*1000的正方形區(qū)域內(nèi)隨機生成個點(坐標(biāo)值是整數(shù))，統(tǒng)計每個坐標(biāo)點上生成的點的個數(shù)。在這個例子下，

8、像例1一樣簡單應(yīng)用union和unique就不行了。通常我們考慮用循環(huán)：p=unidrnd(1000,2);tic;A=zeros(1000);for k=1:A(p(k,1),p(k,2)=A(p(k,1),p(k,2)+1;endtocElapsed time is 0. seconds.而用accumarray,可以這樣： tic;a=accumarray(p,1,1000,1000);tocElapsed time is 0. seconds. isequal(A,a)ans = 1速度整整提高了10倍。例3:1000人，身高分布在170-180，體重在110-200斤，年齡分布在20

9、-50歲，計算身高體重都相等的人的年齡平均值。結(jié)果用矩陣表示：行數(shù)表示身高，列數(shù)表示體重，元素表示年齡的平均值。首先生成數(shù)據(jù)：rand(state,0)height=unidrnd(10,1000,1)+170;rand(state,0)weight=unidrnd(90,1000,1)+110;rand(state,0)old=unidrnd(30,1000,1)+20;利用accumarray計算的語句如下：tic;mo=accumarray(height,weight,old,mean);tocElapsed time is 0. seconds.這個矩陣比較稀疏也可以結(jié)果用稀疏矩陣來表示tic;mo=accumarray(height,weight,old,m