1、有很多线性代数问题都需要生成范德蒙多矩阵，对于一个向量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
2、 %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 i
3、n 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 循环.下
4、面给出的脚本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的各个列,
5、直到第一列为止.之所以可以这样做是因为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; %highes
6、t 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文件所示.%vander
7、5.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文件函
8、数-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和t
9、oc函数来测试一下他们的速度.首先先要删除上面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;t
10、icvander4times(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 循环不一定会生成执行速度最快的代码.有时预先分配内
11、存并小心使用For 循环,从而最大限度的减少内存的使用机会反而会生成最优的代码.在消除For 循环的过程中,后三种方法都使用了更多的内存.前两种最快的方法是vander2.m 和vander3.m,他们都使用了预先分配内存和For循环.vander5.m和vander6.m之间的区别仅仅是使用x(;,ones(1,m)代替了repmat(x,1,m),从他们的测试结果中可以看出,调用repmat不会对执行速度产生很大的影响.testvander.m中的x和m的值可以选择,可以选择适当的参数让方法的运行时间长一些,这样才可以得到足够可靠的时间测试结果.因此,我们选择了非常大的数组来测试这些方法.
12、同样,也可以使用较小的数组来进行测试,这就需要运行多次才能得到可靠的测试结果.下面给出的方法就使用了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;t
13、icfor 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
14、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则是最快的.在Matlab的polyfit函数中使用了vander3.m.如果向量化对内存的要求没有对计算机平台的内存管理能力带来更大需求,那么完全向量化将是最优的解决方法.而另一方面,如果问题的尺寸非常巨大,那么为变量预先分配内存并小心地使用For 循环以减少对内存的需求将是最好的办法.随着范德蒙多矩阵的维数的增长,它的计算将变得越来越复杂,但范德蒙多矩阵通常都用于矩阵尺寸相对较小的应用中.综上所述,在通常情况下,vander6.m是最佳的算法.
