matlab有效生成范德蒙多矩阵

有很多线性代数问题都需要生成范德蒙多矩阵，对于一个向量x,它的范德蒙多矩阵具有如下的形式:V=x1m x1(m-1)x1,1;x2m,x2(m-1),1;.;xnm,xn(m-1)1如上所示,V的各列对x的元素逐个进行了乘方.下面我们将讨论生成此种矩阵的多种方法.首先想到的第一种方法就是直接使用For 循环,具体情况如下面给出的脚本M文件所示.%vander1.m%construct a Vandermonde matrix x=(1:6); %conlumn vector for input data m=5; %highest power to compute V=;for I=1:m+1 %build V conlumn by column V=V x.(m+1-i);end在上面给出的方法中,矩阵V是逐列构建起来的,在最开始矩阵是一个空矩阵.这种方法存在许多缺点,最明显的缺点就是在每次进入循环的时候都要对V重新分配内存.因为向量化的第一步应该是预先为V分配内存,如下面给出的M文件所示.%vander2.m%construct a Vandermonde matrix.x=(1:6); %conlumn vector for input data m=5; %highest power to compute n=length(x); %number of elements in xV=ones(n,m+1); %preallocate memory for result for i=1:m %build V column by columnV(:,i)=x.(m+1-i);end在这里V被作为一个全1的矩阵进行初始化.然后在For 循环中对的各个列进行赋值.在For 循环中,最后一列未被赋值.这是因为最后一列的元素已经都是1了,没有必要在进行x.0的计算.上面给出的代码可以在Matlab中的vander函数里找到.但是上面的代码依旧存在着两个问题.首先,v的各列是直接进行计算的,在计算过程中没有能够利用对前一列进行计算的结果.其次,应该尽量避免使用For 循环.下面给出的脚本M文件解决了第一个问题.%vander3.m%construct a Vandermonde matrix.x=(1:6); %conlumn vector for input data m=5; %highest power to compute n=length(x); %number of elements in xV=ones(n,m+1); %preallocate memory for result for i=m:-1:1 %build V column by column V(:,i)=x.*V(:,i+1);end现在是从矩阵V的第二列元素开始逆向计算V的各个列,直到第一列为止.之所以可以这样做是因为V的第i列等于第i+1列按元素乘以x.这种方法可以在Matlab 的polyfit 函数中找到. 此时如果不消除For循环,则无法进一步对算法进行优化.要消除For循环需要一些灵活性,并要对Matlab中的函数非常熟悉.通过使用本章前面曾经提到过的数组操作表格,函数repmat 和cumprod 可以提供一些有用的功能.下面给出的脚本M文件展示了repmat 函数的用法.%vander4.m%construct a Vandermonde matrix.x=(1:6); %conlumn vector for input data m=5; %highest power to compute n=length(x); %number of elements in xp=m:-1:0; %column powersV=repmat (x,1,m+1).repmat (p,n,1);在这种方法中两次使用了repmat ,一次用于复制x以创建一个每一列都包含x的m+1列的矩阵;另一次用于创建一个含有应用到包含x的矩阵中的每个元素乘方的矩阵.给定这样的两个矩阵,则可以按逐个元素乘幂的方法生成所希望的结果.和vander2.m一样,此方法直接计算每一列,而没有使用其他列的信息.cumprod函数可以解决这个问题,如下面给出的脚本M文件所示.%vander5.m%construct a Vandermonde matrix.x=(1:6); %conlumn vector for input data m=5; %highest power to compute n=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 复制了x 后,使用了cumprod函数来计算V的列.因为cumprod 是从左向右执行的,所以最后的结果需要将V的各列反向.这种方法只使用了一个M文件函数-repmat.可以通过数组寻址来避免使用此函数,这会加速程序的运行速度.使用数组寻址的方法如下面给出的脚本M文件所示.%vander6.m%construct a Vandermonde matrix.x=(1:6); %conlumn vector for input data m=5; %highest power to compute n=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文件测试他们的执行速度了.%testvander.m%script file to test vandermonde implementationsx=randan(10000,1); %column vector for input datam=100; %highest power to computetimes=zeros(1,5);ticvander1times(1)=toc;ticvander2times(2)=toc;ticvander3times(3)=toc;ticvander4times(4)=toc;ticvander5times(5)=toc;relspeed=times/min(times) %relative speed results在笔者的计算机上运行上面的脚本文件会生成如下的结果.relspeed= 24.194 1.632 1 7.2114 2.057 2.0563大家一定认为最后一种方法将是最快的,这是因为它使用了最为向量化的解决方法.但出乎我们的意料,vander3.m是最快的方法,即使它使用了一个For循环也是如此!这些结果很重要,因为他们指出了这样一个事实-消除所有的For 循环不一定会生成执行速度最快的代码.有时预先分配内存并小心使用For 循环,从而最大限度的减少内存的使用机会反而会生成最优的代码.在消除For 循环的过程中,后三种方法都使用了更多的内存.前两种最快的方法是vander2.m 和vander3.m,他们都使用了预先分配内存和For循环.vander5.m和vander6.m之间的区别仅仅是使用x(;,ones(1,m)代替了repmat(x,1,m),从他们的测试结果中可以看出,调用repmat不会对执行速度产生很大的影响.testvander.m中的x和m的值可以选择,可以选择适当的参数让方法的运行时间长一些,这样才可以得到足够可靠的时间测试结果.因此,我们选择了非常大的数组来测试这些方法.同样,也可以使用较小的数组来进行测试,这就需要运行多次才能得到可靠的测试结果.下面给出的方法就使用了For循环来多次运行各种方法.%testvander2.m%scipt file to test vandermonde implementations x=randan(100,1); %column vector for input datam=10; %highest power to computeN=1:500;times=zeros(1,6);ticfor i=N,vander1,endtimes(1)=toc;ticfor i=N,vander2,endtimes(2)=toc;ticfor i=N,vander3,endtimes(3)=toc;ticfor i=N,vander4,endtimes(4)=toc;ticfor i=N,vander5,endtimes(5)=toc;ticfor i=N,vander6,endtimes(6)=toc;relspeed=times/min(times) %relative speed results现在范德蒙多矩阵共有100行,11列,比testvander.m中10000乘以101的矩阵要小很多.在笔者的计算机上运行此脚本文件会生成下面的结果.relspeed= 2.5123 1.209 1.1444 3.8932 1.5075 1我们可以发现结果有了变化.在计算较小的数组的时候,向量化的解决方法vander6.m成为了执行速度最快的方法.另外,使用M文件repmat的vander5.m的执行速度要比vander6.m慢50%左右.最优使用For循环的vander3.m仍旧很有竞争力,它只比最快的vander6.m慢了14%.即使是程序结构最差的vander1.m也只比最快的vander6.m慢了2.5倍.而在使用了大型数组的testvander.m中,最快与最慢的方法在执行速度之间的差距可以达到24倍之多.那么上面给出的6种方法中,到底哪一种方法才是最佳的呢?对于非常大型的数组而言,vander3.m是最快的,而对于一般的数组而言,vaner6.m则是最快