撰写快速的 MATLAB 程序代码.docx
- 文档编号:9272562
- 上传时间:2023-02-03
- 格式:DOCX
- 页数:19
- 大小:212.11KB
撰写快速的 MATLAB 程序代码.docx
《撰写快速的 MATLAB 程序代码.docx》由会员分享,可在线阅读,更多相关《撰写快速的 MATLAB 程序代码.docx(19页珍藏版)》请在冰豆网上搜索。
撰写快速的MATLAB程序代码
撰写快速的MATLAB程序代码
MATLAB程序语言是剖析式(parsed)的,原始码是实时直译的。
像C++和Fortran的程序语言比较快,因为它们在变成计算机原生(native)语言前先被编译。
实时剖析的优点是较与平台不相关、具强固性(robustness)和更容易除错。
缺点是在速度上会变慢,负荷会变多,及受限的低阶控制。
要弥补速度上的损失,MATLAB提供了能帮助加速的程序代码。
此文章讨论这些和其它的策略来改善MATLAB原始码的速度。
*数组预先配置
*向量化的计算和逻辑
*向量化的参考(Vectorizedreferencing)
警告!
*先学习程序语言:
最佳化(不仅是MATLAB)需要熟悉程序语言的语法和函式。
本文并不是MATLAB的教学。
*使用批注:
最佳化的程序代码是精练和神秘的。
为了帮助他人和你自己,记得要加批注。
*不要在它需要被最佳化程序代码之前最佳化它:
在最佳化程序代码之前,考虑它是否值得这些心力。
如果程序代码很快就会被修改或扩充,无论如何它会被重写。
*只要最佳化必要的部份:
确定程序确实是速度的瓶颈。
如果它不是,最佳化只会把程序代码变得难以理解。
1.速度剖析器(TheProfiler)
MATLAB5.0和较新的版本都包含了一个叫作profiler的工具,它帮助你决定程序里哪里是瓶颈。
考虑下面的程序:
functionresult=example1(Count)
fork=1:
Count
result(k)=sin(k/50);
ifresult(k)<-0.9
result(k)=gammaln(k);
end
end
要分析程序的效率,首先启动profiler,并清除任何旧的profiler数据:
>>profileon
>>profileclear
现在执行程序吧。
将输入的自变量变得大些或小些,让它跑个几秒。
>>example1(5000);
然后输入
>>profreport(’example1’)
profiler会在函式上产生一个HTML报表,然后开启一个浏览器窗口。
根据系统,profiler的结果可能会与此例有些不同。
在蓝色的「example1」连结上按一下会提供更多的细节:
显示了最花时间的行数,伴随着时间、时间的百分比和行数。
代价最高的行数是第4和第7行。
2.预先配置数组
MATLAB的矩阵变量有动态的自变量行和列的能力。
例如:
>>a=2
a=
2
>>a(2,6)=1
a=
200000
000001
MATLAB会自动地改变矩阵的尺寸大小。
在内部,矩阵数据存储器必须被预先配置予较大的尺寸。
如果矩阵可重复地改变大小--像是在for循环里--此负荷会变得很显著。
要避免经常性的重新配置,用zero命令来「预先配置」矩阵。
思考下列的程序代码:
a
(1)=1;
b
(1)=0;
fork=2:
8000
a(k)=0.99803*a(k-1)-0.06279*b(k-1);
b(k)=0.06279*a(k-1)+0.99803*b(k-1);
end
profiler计算此程序代码的时间花费了0.47秒。
在执行for循环之后,二个数组都是长度8000的列向量,因此要预先配置,建立8000个元素的空的a和b列向量。
a=zeros(1,8000);%预先配置
b=zeros(1,8000);
a
(1)=1;
b
(1)=0;
fork=2:
8000
a(k)=0.99803*a(k-1)-0.06279*b(k-1);
b(k)=0.06279*a(k-1)+0.99803*b(k-1);
end
经此修改,程序代码只花费0.14秒(快超过三倍)。
预先配置通常很容易,在本例中它只需要决定正确的预先配置大小及新增二行。
如果最终的数组尺寸大小会变动那么会怎么样?
一个可能的方法,是在循环执行完之后使用数组尺寸大小的上界(upperbound),并且裁剪掉超出的部份:
a=zeros(1,10000);%预先配置
count=0;
fork=1:
10000
v=exp(rand
(1)*rand
(1));
ifv>0.5%conditionallyaddtoarray
count=count+1;
a(count)=v;
end
end
a=a(1:
count);%裁剪结果
此程序的平均的执行时间是0.42秒(没有预先配置),及0.18秒(有预先配置)。
巢状数组(cellarray)也能被预先配置,使用cell命令来建立想要的尺寸大小。
对常常变化大小的巢状数组预先配置内存,甚至比对双精度数组更有效益。
3.向量化(Vectorization)
要「向量化」一个计算代表说要用向量运算取代平行运算。
此策略通常增进十倍的速度。
好的向量化是一种必须被发展的技巧,它需要熟悉MATLAB程序语言和有创造力。
3.1向量化的运算
许多标准的MATLAB函式是「向量化的」,它们能在数组上操作,就如同函式已经各别地被套用于每一个元素之上。
>>sqrt([1,4;9,16])
ans=
12
34
>>abs([0,1,2,-5,-6,-7])
ans=
012567
考虑下列的函式:
functiond=minDistance(x,y,z)
%找出介于一个点集合和原点间旳最小距离
nPoints=length(x);
d=zeros(nPoints,1);%预先配置
fork=1:
nPoints%计算每一点的距离
d(k)=sqrt(x(k)^2+y(k)^2+z(k)^2);
end
d=min(d);%取得最小距离
对于每一个点,介于它和原点间的距离被计算且储存在d。
然后最小距离可用min来找到。
要向量化距离的运算,用向量运算代换for循环:
functiond=minDistance(x,y,z)
%找到介于一个点集合和原点之间的最小距离
d=sqrt(x.^2+y.^2+z.^2);%计算每一点的距离
d=min(d);%取得最小的距离
修正后的程序代码以向量运算来执行距离运算。
x,y和z数组首先使用逐元素的power运算(.^,用来做乘法和除法的逐元素运算子是.*和./)来计算平方。
平方的部份用向量加法相加。
最后向量总和的平方根被逐元素计算,产生了一个距离数组。
minDistance程序的第一版50000个点花了0.73秒。
向量化的版本花费的时间少于0.04秒,快了超过18倍。
一些对于向量化计算有用的函式:
min,max,repmat,meshgrid,sum,cumsum,diff,prod,cumprod,filter。
3.2向量化的逻辑(VectorizedLogic)
前一节告诉你如何去向量化纯粹的计算。
瓶颈程序代码通常包含条件式逻辑。
就像计算,MATLAB的逻辑运算子被向量化了:
>>[1,5,3]<[2,2,4]
ans=
101
二个数组被逐元素比较。
逻辑运算以二元的元素传回「logical」数组。
这如何使用?
对于在逻辑数组运算上,MATLAB有一些强力的函式:
*find:
找出非零元素的索引值(index)。
*any:
如果向量中有任何一个元素为非零(或是一个矩阵中的每一行),则为True。
*all:
如果向量中所有的元素为非零(或是一个矩阵中的每一行),则为True。
>>find([1,5,3]<[2,2,4])
ans=
13
>>find(eye(3)==1)
ans=
1
5
9
find函式传回「向量逻辑运算传回True」的索引(index)。
在第一个命令中,1<2为真,5<2为伪,而3<4为真,所以find回报说第一个和第三个比较为真。
在第二个命令中find传回identity矩阵等于1的索引。
索引1,5和9对应到3x3矩阵的对角线。
any和all函式很简单,不过偶尔非常有用。
如果任何元素为真,any会传回true。
如果所有的元素为真,all会传回true。
范例1:
移除元素
当在某些逐元素的情况下,数组元素必须被移除时,此情况通常会产生。
例如,这段程序代码从数组x中,移除所有的NaN和inf元素:
i=find(isnan(x)|isinf(x));%找出x中的坏元素
x(i)=[];%删除它们
另一个方式是,
i=find(~isnan(x)&~isinf(x));%找出不是Nan也不是inf的元素
x=x(i);%保留那些元素
范例2:
片段的(peicewise)函式
在许多应用程序中,需要以片段的(piecewise)定义来计算函式(例如,splineinterpolants)。
sinc函式特别有趣,因为在0处有可移除的singularity。
对于x~=0,sinc能以y=sin(x)./x来计算。
此程序使用find和向量化计算来各别处理二个案例:
functiony=sinc(x)
%对于一堆x值,逐元素计算sinc函式。
y=ones(size(x));%将y全设为0,因为sinc(0)=1
i=find(x~=0);%找出所有非零的x值
y(i)=sin(x(i))./x(i);%计算sinc,其中x~=0
范例3:
以meshgrid来绘制图像
meshgrid函式要二个输入向量并且藉由复制第一个向量于各列上,和第二个向量于各行上,来将它们转换它们为矩阵。
>>[x,y]=meshgrid(1:
5,1:
3)
x=
12345
12345
12345
y=
11111
22222
33333
上面的矩阵运作的方式像是一个宽度为5,高度为3影像的地图。
对于已知的像素,位置x可以由x读取,而位置y可由y读取。
这似乎像是一个内存的免费(gratuitous)用法,x和y简单地纪录了行和列的位置,不过这很有用。
例如,要画出椭圆(ellipse)。
%对一个宽150,高100的影像来建立x和y
[x,y]=meshgrid(1:
150,1:
100);
%原点为(60,50),尺寸为15x40的椭圆
Img=sqrt(((x-60).^2/15^2)+((y-50).^2/40^2))>1;
%画出影像
imagesc(Img);colormap(copper);
axisequal;axisoff;
画线几乎是一样,只要在方程式中改变即可
[x,y]=meshgrid(1:
150,1:
100);
%线y=x*0.8+20
Img=(abs((x*0.8+20)-y)>1);
imagesc(Img);colormap(copper);
axisequal;axisoff;
Polar函式可以被藉由先将x和y变量以cart2pol函式来转换,然后被绘出。
[x,y]=meshgrid(1:
150,1:
100);
[th,r]=cart2pol(x-75,y-50);%转成极坐标
%对中心(75,50)旋转
Img=sin(r/3+th);
imagesc(Img);colormap(hot);
axisequal;axisoff;
范例4:
多项式内插
有时候不仅仅会需要执行运算于二个向量的每一对元素间,也要对每个元素的组合。
此范例会在polynomialfitting算法中告诉你如何做这件事。
己知n个点x1,x2,x3,...xn和n个对应函式值y1,y2,y3,...yn,n-1维度内插多项式系数c0,c1,c2,...,cn-1可以藉由解矩阵方程式而找到。
MATLAB早就准备好去解决这样的矩阵方程式,不过困难的部份要被设定。
functionc=polyint(x,y)
%给定一个点集合和函式值x和y,计算内插函式。
x=x(:
);%确定x和y都是行向量
y=y(:
);
n=length(x);%n=点数
%%%建构矩阵的左手边%%%
xMatrix=repmat(x,1,n);%建立一个nxn矩阵,每一行都是x
powMatrix=repmat(n-1:
-1:
0,n,1);%建立另一个nxn指数矩阵
A=xMatrix.^powMatrix;%计算指数(power)
c=Ay;%解决系数的矩阵方程式
建构左手边矩阵的策略首先是建立二个底数(base)和指数(exponent)的nxn矩阵,然后将它们用逐元素的power操作数.^来放在一起。
repmat函式(「复制矩阵」)被用来建立base矩阵xMatrix和指数矩阵powMatrix。
xMatrix藉由在行上面重复行向量xn次来形成。
相似地,powMatrix是一个有n-1,n-2,n-3,...,0重复递减列向量n次。
MATLAB4.0和较早的版本没有repmat函式,请参考6.5节里的另一个方式。
二个矩阵也可以用meshgrid来建立:
[powMatrix,xMatrix]=meshgrid(n-1:
-1:
0,x);
此函式只有一个范例;对重要的(serious)多项式内插法使用标准的polyfit函式。
左手边矩阵用像号角(Horner-like)的建构(polyfit.m45-48行)方式更通用且更具算法效率。
4参考运算(ReferencingOperations)
MATLAB里的参考(Referencing)是多变且强到足以让我们保留一节来讨论的。
透彻地了解参考(Referencing)让你对编程中各种情况能加以向量化。
4.1下标vs.索引
下标是用来参考矩阵元素最常见的方法,例如,A(3,9)指到第三列,第九行。
索引是另一个参考方法。
考虑一个10x10矩阵A。
在内部,MATLAB线性地储存矩阵数据为一个一维的、100个元素的数据数组。
一个索引指到一个元素在此一维数组中的位置。
例如,A(83)指到在第三列,第九行上的元素。
介于下标和索引间的转换可以sub2ind和ind2sub函式来达成。
然而,因为这些是m-file函式,而非快速的内建运算,它更有效率地直接计算转换。
对一个二维的、尺寸大小为MxN的矩阵A,介于下标(i,j)和索引(index)间的转换为:
A(i,j)<-->A(i+(j-1)*M)
A(index)<-->A(rem(index-1,M)+1,floor(index/M)+1)
索引也能延伸到N维矩阵,以先增加行、然后是列,然后是第三个维度,等等....的索引。
下标则是直观地扩充了,A(...,dim4,dim3,row,col)。
4.2向量化的下标
用子矩阵(submatrices)而非各别的元素是很有用的。
这用索引或下标的向量来达成。
如果A是一个二维矩阵,向量下标参考有下列语法:
A(rowv,colv)
其中rowv是一个有M个元素的列向量,而colv是一个有N个元素的行向量。
二者可能是任意长度且它们的元素可能是依任何顺序排成的。
如果是一个矩阵的话,它会被重新塑型(reshaped)成一个向量。
在向量下标中使用列向量或行向量是没有差别的。
在运算的右手边(rhs),向量下标参考会传回一个元素尺寸大小为MxN的子矩阵。
如果向量下标矩阵在左手边(lhs,left-handside),rhs的结果必须是MxN或是纯量大小的。
如果有任何元素在目的地的参考被重复了,例如,A(1,2)和A(2,2)这种不明确、令人混淆的指派:
A([1,2],[2,2])=[1,2;3,4]
这是在其较大索引的来源数组中指定的(dominate)值。
>>A=zeros
(2);A([1,2],[2,2])=[1,2;3,4]
A=
02
04
向量下标参考在较高维度中直观地延伸扩充了。
4.3向量索引(VectorIndices)
多重元素也可被以向量索引参考。
A(indexv)
其中indexv是一个索引数组。
在右手边(rhs),向量索引参考传回和indexv尺寸大小相同的矩阵,如果indexv是一个3x4矩阵,A(indexv)是3x4矩阵:
然而向量下标参考到块状(block-shaped)子矩阵是有限制的,向量索引可以参考到任意形状。
如果向量索引参考是在左手边(lhs),右手边(rhs)一定会传回一个和indexv同样大小的矩阵或纯量。
使用向量下标,混淆不清重复的指派会使用较后面所指派的值。
>>A=zeros
(2);A([3,4,3,4])=[1,2,3,4]
A=
03
04
4.4参考的万用字符
在一个指到整列或整行的下标中使用万用字符,:
。
例如,A(:
1)指到第一行中的每一列--整个第一行。
这可以结合向量下标,A([2,4],:
)指到第二和第四列。
当万用字符被用于向量索引时,整个矩阵会被参考到。
在右手边(rhs),这总是会传回一个行向量。
A(:
)=行向量
这常常有用:
例如,如果一个函式输入必须是一个列向量,使用者的输入可以被以A(:
).'快速地重新塑形成列向量。
A(:
).'是转成一个行向量然后转置成列向量(译注:
A(:
)'也行)。
在左手边(lhs)的A(:
)指派了所有A的元素,不过不会改变其尺寸大小。
例如,A(:
)=8将所有A矩阵的元素变为8。
4.5以[]删除子矩阵
在一个矩阵中的元素可以藉由指派空矩阵来删除之。
例如,A([3,5])=[]从A删除了第三和第五个元素。
如果这以索引参考来达成,矩阵会被重新塑形(reshape)成列向量。
如果除了一个以外,所有的下标是万用字符,也能以下标来删除。
A(2,:
)=[]删除了第二列。
像A(2,1)=[]或甚至是A(2,1:
end)=[]是不被允许的。
5.杂七杂八的技巧
5.1转换任何数组成一个行向量
强制一个数组成行向量通常是有用的,例如,当撰写一个期望行向量作为输入的函式时。
此简易的技巧会将输入数组转换成一个行向量(如果数组是一个列向量、矩阵、N维数组或已经是一个行向量)。
x=x(:
);%将x转换成一个行向量
藉由在此运算后加上转置操作数「.'」,便可以将数组转成列向量了。
5.2取得数组中的元素个数
size函式传回一个内含数组尺寸大小的向量。
用prod,数组中元素的个数可以被简要地计算,用
nElements=prod(size(x))%计算x里的元素个素
用ndims(x)来取得数组中维度的数目。
注意:
在MATLAB6.1时,MATLAB引入了一个numel函式来计算元素的个数,不过此函式不回溯相容。
5.3不使用if叙述来限制(bound)一个值
如果问题是「x值必须至少是a,至多为b」,直观的方法是用if和elseif叙述式来写这段程序代码。
不幸地是,MATLAB的if叙述很慢。
使用min和max函式来将x限制在lowerBound和upperBound范围之内的方法是:
x=max(x,lowerBound);%令x>=lowerBound
x=min(x,upperBound);%令x<=upperBound
如果x是一个任意大小的矩阵,这一样也会作用于每个元素之上。
5.4找出矩阵或N维数组的的min/max
给定一个矩阵输入(没有其它的输入),min和max函式沿着行作运算,它会找出每行内的极值。
通常找出整个矩阵元素的极值会更有用。
重复min或max函式于行极值上是合理的;
min(min(A))用来找出一个二维矩阵A的最小值。
此方法只对min/max呼叫一次(使用转换到行向量的技巧)且决定了极值元素的位置:
[MinValue,MinIndex]=min(A(:
));%找到A中最小的元素
%极小值为MinValue,它的索引为MinIndex
MinSub=ind2sub(size(A),MinIndex);%将MinIndex转换为下标
最小的元素是A(MinIndex)或以下标参考的方式是A(MinSub
(1),MinSub
(2),...)。
(类似地,用max来取代min以求得极大值。
)
5.5不以repmat来「重复向量/将向量如瓦片般铺设」(Repeating/tilingvectorswithoutrepmat)
在MATLAB5.3版时,repmat函式以m-file而非内建函式来实作。
虽然它非常快,它仍有因使用m-file而带来的间接成本。
因此最标准的MATLABm-file函式内嵌地执行repmat,而非呼叫它(例如,请见meshgrid.m41,42,55行),这并不令人惊讶。
然而repmat有为了运算于矩阵上的泛用性,而通常我们只需要排列一个向量或纯量而已。
要重复一个行向量y于行上n次,
A=y(:
ones(1,n));%相当于A=repmat(y,1,n);
要在列上重复一个列向量xm次,
A=x(ones(1,m),:
);%相当于A=repmat(x,m,1);
要重复一个纯量s到一个mxn的矩阵内
A=s(ones(m,n));%相当于A=repmat(s,m,n);
此方法避免了呼叫m-file函式的负担。
它决不会比repmat慢。
鉴定家应该注意到repmat.m本身在第52-53行上使用了这个方法。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 撰写快速的 MATLAB 程序代码 撰写 快速
![提示](https://static.bdocx.com/images/bang_tan.gif)