1、有很多線性代數問題都需要生成范德蒙多矩陣，對于一個向量X,它的范德蒙多矩陣具 有如下的形式：V=x1八m x1八(m-1)x1,1;x2八m,x2八(m-1),1;.；xnm,xn八(m-1)1如上所示,V的各列對x的元素逐個進行了乘方.下面我們將討論生成此種矩陣的多種方 法.首先想到的第一種方法就是直接使用For循環,具體情況如下面給出的腳本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是逐列構建起來的,在最開始矩陣v是一個空矩陣這種方 法存在許多缺點,最明顯的缺點就是在每次進入循環的時候都要對V重新分配內存.因為向 量化的第一步應該是預先為V分配內存,如下面給出的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被作為一個全1的矩陣進行初始化.然后在For循環中對v的各個列進行賦值. 在For循環中，最后一列未被賦值.這是因為最后一列的元素已經都是1 了，沒有必要在進行 x0的計算.上面給出的代碼可以在Mat lab中的vander函數里找到.但是上面的代碼依舊 存在著兩個問題.首先,v的各列是直接進行計算的,在計算過程中沒有能夠利用對前一列進 行計算的結果.其次，應該盡量避免使用For循環.下面給出的腳本M文件解決了第一

4、個問 題.%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現在是從矩陣V的第二列元素開始逆向計算V的各個列,直到第一列為止.之所以可以這 樣做是因為V的第i列

5、等于第i+1列按元素乘以x.這種方法可以在Mat lab的polyfit函 數中找到.此時如果不消除For循環，則無法進一步對算法進行優化.要消除For循環需要一些 靈活性，并要對Matlab中的函數非常熟悉.通過使用本章前面曾經提到過的數組操作表格, 函數repmat和cumprod可以提供一些有用的功能.下面給出的腳本M文件展示了 repmat 函數的用法.%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 , 一次用于復制x以創建一個每一列都包含x的m+1 列的矩陣;另一次用于創建一個含有應用到包含x的矩陣中的每個元素乘方的矩陣給定這 樣的兩個矩陣，則可以按逐個元素乘幕的方法生成所希望的結果.和一樣，此方法直接計算每 一列,而沒有使用其他列的信息.cumprod函數可以解決這個問題,如下面給出的腳本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復制了乂后，使用了cumpr。d函數來計算V的列.因為cumprod 是從左向右執行的,所以最后的結果需要將V的各列反向.這種方法只使用了一個M文件函數 -repmat.可以通過數組尋址來避免使用此函數,這會加速程序的運行速度.使用數組尋址 的方法如下面給出的腳本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函數來測試一下他們的速度.首先先 要刪除上面6種方法中的前兩行:x(1:6);和m=5;.然后就可以使用下面的腳本M文件測試 他們的執行速度了.%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 在筆者的計算機上運行上面的腳本文件會生成如下的結果. relspeed=1大家一定認為最后一種方法將是最快的，這是因為它使用了最為向量化的解決方法.但 出乎我們的意料，是最快的方法，即使它使用了一個For循環也是如此!這些結果很重要，因 為他們指出了這樣一個事實消除所有的For循環不一定會生成執行速度最快的代碼. 有時預先分配內存并小心使用For循環,從而最大限度的減少內存的使用機會反而會生成 最優的代碼.在消除For循環的過程中，后三種方法都使用了更多的內存.前兩種最快的方法是和, 他們都使用了預先分配內存和For循環.和之間的區別僅

11、僅是使用x(;,ones(1,m)代替了 repmat(x,1,m),從他們的測試結果中可以看出，調用repmat不會對執行速度產生很大的影 響.中的x和m的值可以選擇,可以選擇適當的參數讓方法的運行時間長一些,這樣才可以得 到足夠可靠的時間測試結果.因此,我們選擇了非常大的數組來測試這些方法.同樣,也可以 使用較小的數組來進行測試，這就需要運行多次才能得到可靠的測試結果.下面給出的方法 就使用了For循環來多次運行各種方法.%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現在范德蒙多矩陣共有100行,11歹1比中10000乘以101的矩陣要小很多.在筆者的計 算機上運行此腳本文件會生成下面的結果.relspeed=1我們可以發現結果有了變化.在計算較小的數組的時候，向量化的解決方法成為了執行 速度最快的方法.另外，使用M文件repmat的的執行速度要比慢50%左右.最優使用For循環 的仍舊很有競爭力,它只比最快的慢了 14%.即使是程序結構最差的也只比最快的慢了倍.而 在使用了大型數組的中，最快與最慢的方法在執行速度之間的差距可以達到24倍之多.那么上面給出的6種方法中，到底哪一種方法才是最佳的呢對于非常大型的數組而言, 是最快的,而對于一般的數組而言,則是最快的.