matlab有效生成范德蒙多矩陣

1、有很多線性代數(shù)問題都需要生成范德蒙多矩陣，對(duì)于一個(gè)向量X,它的范德蒙多矩陣具 有如下的形式：V=x1八m x1八(m-1)x1,1;x2八m,x2八(m-1),1;.；xnm,xn八(m-1)1如上所示,V的各列對(duì)x的元素逐個(gè)進(jìn)行了乘方.下面我們將討論生成此種矩陣的多種方 法.首先想到的第一種方法就是直接使用For循環(huán),具體情況如下面給出的腳本M文件所 示.%construct a Vandermonde matrixx=(1:6),; %conlumn vector for input datam=5; %highest power to computeV=;for I=1:m+1 %bui

2、ld V conlumn by columnV=V x.(m+1-i);end在上面給出的方法中,矩陣V是逐列構(gòu)建起來的,在最開始矩陣v是一個(gè)空矩陣這種方 法存在許多缺點(diǎn),最明顯的缺點(diǎn)就是在每次進(jìn)入循環(huán)的時(shí)候都要對(duì)V重新分配內(nèi)存.因?yàn)橄?量化的第一步應(yīng)該是預(yù)先為V分配內(nèi)存,如下面給出的M文件所示.%construct a Vandermonde matrix.x=(1:6),; %conlumn vector for input datam=5;%highest power to computen=length(x); %number of elements in xV=ones(n,m+1)

3、; %preallocate memory for resultfor i=1:m %build V column by columnV(:, i)=x(m+1-i);end在這里V被作為一個(gè)全1的矩陣進(jìn)行初始化.然后在For循環(huán)中對(duì)v的各個(gè)列進(jìn)行賦值. 在For循環(huán)中，最后一列未被賦值.這是因?yàn)樽詈笠涣械脑匾呀?jīng)都是1 了，沒有必要在進(jìn)行 x0的計(jì)算.上面給出的代碼可以在Mat lab中的vander函數(shù)里找到.但是上面的代碼依舊 存在著兩個(gè)問題.首先,v的各列是直接進(jìn)行計(jì)算的,在計(jì)算過程中沒有能夠利用對(duì)前一列進(jìn) 行計(jì)算的結(jié)果.其次，應(yīng)該盡量避免使用For循環(huán).下面給出的腳本M文件解決了第一

4、個(gè)問 題.%construct a Vandermonde matrix.x=(1:6),; %conlumn vector for input datam=5;%highest power to computen=length(x); %number of elements in xV=ones(n,m+1); %preallocate memory for resultfor i=m：-1:1 %build V column by columnV(:, i)=x.*V(:, i+1);end現(xiàn)在是從矩陣V的第二列元素開始逆向計(jì)算V的各個(gè)列,直到第一列為止.之所以可以這 樣做是因?yàn)閂的第i列

5、等于第i+1列按元素乘以x.這種方法可以在Mat lab的polyfit函 數(shù)中找到.此時(shí)如果不消除For循環(huán)，則無法進(jìn)一步對(duì)算法進(jìn)行優(yōu)化.要消除For循環(huán)需要一些 靈活性，并要對(duì)Matlab中的函數(shù)非常熟悉.通過使用本章前面曾經(jīng)提到過的數(shù)組操作表格, 函數(shù)repmat和cumprod可以提供一些有用的功能.下面給出的腳本M文件展示了 repmat 函數(shù)的用法.%construct a Vandermonde matrix.x=(1:6);%conlumnvector forinputdatam=5;%highest power to computen=length(x);%numberof

6、elementsin xp=m:-1:0;%columnpowersV=repmat (x,1,m+1)repmat (p,n,1);在這種方法中兩次使用了 repmat , 一次用于復(fù)制x以創(chuàng)建一個(gè)每一列都包含x的m+1 列的矩陣;另一次用于創(chuàng)建一個(gè)含有應(yīng)用到包含x的矩陣中的每個(gè)元素乘方的矩陣給定這 樣的兩個(gè)矩陣，則可以按逐個(gè)元素乘幕的方法生成所希望的結(jié)果.和一樣，此方法直接計(jì)算每 一列,而沒有使用其他列的信息.cumprod函數(shù)可以解決這個(gè)問題,如下面給出的腳本M文件 所示.%construct a Vandermonde matrix.x=(1:6); %conlumn vector f

7、or input datam=5;%highest power to computen=length(x); %number of elements in xV=ones(n,m+1);V(;,2:end)二cumprod(repmat(x,1,m),2);V=V(:,m+1：-1:1);在這里，在使用了repmat復(fù)制了乂后，使用了cumpr。d函數(shù)來計(jì)算V的列.因?yàn)閏umprod 是從左向右執(zhí)行的,所以最后的結(jié)果需要將V的各列反向.這種方法只使用了一個(gè)M文件函數(shù) -repmat.可以通過數(shù)組尋址來避免使用此函數(shù),這會(huì)加速程序的運(yùn)行速度.使用數(shù)組尋址 的方法如下面給出的腳本M文件所示.%co

8、nstruct a Vandermonde matrix.x=(1:6); %conlumn vector for input datam=5;%highest power to computen=length(x); %number of elements in xV=ones(n,m+1);V(;,2:end)二cumprod(x(:,ones(1,m),2);V=V(;,m+1：-1:1);上面一共給出了 6種方法.下面我們使用tic和toc函數(shù)來測(cè)試一下他們的速度.首先先 要?jiǎng)h除上面6種方法中的前兩行:x(1:6);和m=5;.然后就可以使用下面的腳本M文件測(cè)試 他們的執(zhí)行速度了.%s

9、cript file to test vandermonde implementationsx=randan(10000,1);%column vector for input datam=100;times=zeros(1,5);m=100;times=zeros(1,5);%highest power to computeticvanderltimes(1)=toc;ticvander2times(2)=toc;ticvander3times(3)=toc;ticvander4times(4)=toc;ticvander5times(5)=toc;relspeed=times/min(ti

10、mes) %relative speed results 在筆者的計(jì)算機(jī)上運(yùn)行上面的腳本文件會(huì)生成如下的結(jié)果. relspeed=1大家一定認(rèn)為最后一種方法將是最快的，這是因?yàn)樗褂昧俗顬橄蛄炕慕鉀Q方法.但 出乎我們的意料，是最快的方法，即使它使用了一個(gè)For循環(huán)也是如此!這些結(jié)果很重要，因 為他們指出了這樣一個(gè)事實(shí)消除所有的For循環(huán)不一定會(huì)生成執(zhí)行速度最快的代碼. 有時(shí)預(yù)先分配內(nèi)存并小心使用For循環(huán),從而最大限度的減少內(nèi)存的使用機(jī)會(huì)反而會(huì)生成 最優(yōu)的代碼.在消除For循環(huán)的過程中，后三種方法都使用了更多的內(nèi)存.前兩種最快的方法是和, 他們都使用了預(yù)先分配內(nèi)存和For循環(huán).和之間的區(qū)別僅

11、僅是使用x(;,ones(1,m)代替了 repmat(x,1,m),從他們的測(cè)試結(jié)果中可以看出，調(diào)用repmat不會(huì)對(duì)執(zhí)行速度產(chǎn)生很大的影 響.中的x和m的值可以選擇,可以選擇適當(dāng)?shù)膮?shù)讓方法的運(yùn)行時(shí)間長一些,這樣才可以得 到足夠可靠的時(shí)間測(cè)試結(jié)果.因此,我們選擇了非常大的數(shù)組來測(cè)試這些方法.同樣,也可以 使用較小的數(shù)組來進(jìn)行測(cè)試，這就需要運(yùn)行多次才能得到可靠的測(cè)試結(jié)果.下面給出的方法 就使用了For循環(huán)來多次運(yùn)行各種方法.%scipt file to test vandermonde implementationsx=randan(100,1); %column vector for in

12、put datam=10;%highest power to computeN=1:500;times=zeros(1,6);ticfor i=N,vander1,endtimes(1)=toc;ticfor i=N,vander2,end times(2)=toc;ticfor i=N,vander3,end times(3)=toc;ticfor i=N,vander4,end times(4)=toc;ticfor i=N,vander5,end times(5)=toc;ticfor i=N,vander6,endtimes(6)=toc;relspeed=times/min(time

13、s) %relative speed results現(xiàn)在范德蒙多矩陣共有100行,11歹1比中10000乘以101的矩陣要小很多.在筆者的計(jì) 算機(jī)上運(yùn)行此腳本文件會(huì)生成下面的結(jié)果.relspeed=1我們可以發(fā)現(xiàn)結(jié)果有了變化.在計(jì)算較小的數(shù)組的時(shí)候，向量化的解決方法成為了執(zhí)行 速度最快的方法.另外，使用M文件repmat的的執(zhí)行速度要比慢50%左右.最優(yōu)使用For循環(huán) 的仍舊很有競(jìng)爭(zhēng)力,它只比最快的慢了 14%.即使是程序結(jié)構(gòu)最差的也只比最快的慢了倍.而 在使用了大型數(shù)組的中，最快與最慢的方法在執(zhí)行速度之間的差距可以達(dá)到24倍之多.那么上面給出的6種方法中，到底哪一種方法才是最佳的呢對(duì)于非常大型的數(shù)組而言, 是最快的,而對(duì)于一般的數(shù)組而言,則是最快的.