撰写快速的 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行上使用了这个方法。
升级会员 