提高MATLAB程序效率的几点原则，这些都是俺在这两年中参加四次数模编写大量m程序总结的经验，加之网上很多英雄也是所见略同。这里先介绍最最重要三个原则。
1.“计算向量、矩阵化，尽量减少for循环。”
因为MATLAB本来就是矩阵实验室的意思，他提供了极其强大而灵活的矩阵运算能力，你就没必要自己再用自己编写的for循环去实现矩阵运算的功能了。另外由于matlab是一种解释性语言，所以最忌讳直接使用循环语句。但在有些情况下，使用for循环可以提高程序的易读性，在效率提高不是很明显的情况下可以选择使用for循环。
口说无凭，下面是利用tic与toc命令计算运算所用时间的方法，测试两种编程的效率。需要说明的是没有办法精 确计算程序执行时间，matlab帮助这样写到“Keep in mind that tic and toc measure overall elapsed time. Make sure that no other applications are running in the background on your system that could affect the timing of your MATLAB programs.”意思是说在程序执行的背后很可能有其他程序在执行，这里涉及到程序进程的的问题，m程序执行的过程中很可能有其他进程中断m程序来利 用cup，所以计算出来的时间就不仅是m程序的了，在这种情况下我把那些寄点去掉，进行多次计算求他的平均时间。
n = 100;
A(1:1000,1:1000) = 13;
C(1:1000,1) = 15;
D(1:1000,1) = 0;
for k = 1:n
D(:) = 0;
tic
for i = 1:1000
for j = 1:1000
D(i) = D(i) + A(i,j)*C(j);
end
end
t1(k) = toc;
%------------------
D(:) = 0;
tic
D = A*C;
t2(k) = toc;
end
u = t1./t2;
u(u<0) = [];
plot(u)
p = mean(u)
t1、t2分别代表用for循环编程和矩阵化编程计算矩阵乘向量所用时间，u代表时间的比值。u(u<0) = [];是认为t1不可能小于t2，所以去掉不可能出现的情况。然后画出图形计算平均值。
经多次试验结果大致相同，其中一次结果如下：
p =
9.6196
------------t1时间是t2的十倍左右。
2.“循环内大数组预先定义--预先分配空间”
这一点原则是极其重要的，以至于在编写m程序时编辑器会给出提示“'ver' might be growing inside a loop.Consider prealloacting for speed.”——这里插一句，最新的MATLAB R2008a版本编辑程序时对很多错误的编程习惯会给出警告，功能很完美。
clear
n = 50;
m = 1000;
for k = 1:n
A = [];
tic
A(1:m,1:m) = 3;
for i = 1:m
A(i,i) = i;
end
t1(k) = toc;
%------------
A = [];
tic
for j = 1:m
A(j,j) = j;
end
t2(k) = toc;
end
t2(t1>10^9) = [];
t1(t1>10^9) = [];
plot([t1;t2]')
t1、t2分别表示预先分配空间和循环中分配空间的时间，下图上面一条t2、下面t1
3.“尽可能利用matlab内部提供的函数”
因为matlab内部提供的函数绝对是各种问题的最优算法，那写程序都是他们大师级人物写出来的，程序应该说相当高效，有现成的为什么不用那！
这个原则就不用实际的程序测试了。
实际上还有很多具体的编程技巧，在接下来的一段时间我会相继写出或者由英语文章翻译过来。
skip to main |
skip to sidebar
2 条评论:
下面是摘自振动论坛eight的文章：
《提高matlab运行速度的一点心得》
首先推荐使用matlab 2006a版本，该版本优点很多（不过有一个小bug，就是通过GUI自动生成的m文件居然一大堆warning，希望在已经发布了的2006b版本中有改善），其中对于编程人员来说比较突出的一个就是编辑窗口的自动语法检查功能。这可以在一定程度上避免使用没有被定义或赋值的变量，另外，也可以帮助你优化代码，【例1】的【方案3】就是因为我看到matlab编辑窗口的warning而得到的启发。顺便提一下，虽然matlab不像其他语言那样，对变量采用“先定义，后使用”的规则，但是，我的经验是，在使用一个变量之前，最好先对它进行“定义”，这里的“定义”是指为它分配空间，这样不但可以提高运行的速度（这在matlab的帮助中也提到，详见Preallocating Arrays一节），而且还可以减少出错的几率，特别是在循环赋值、且变量大小不固定的时候（对此可参阅这个帖子：http://www.chinavib.com/forum/viewthread.php?tid=23732]）。
下面说说如何对matlab提速的问题，我会使用两个例子来说明。
【例1】任务描述：根据A的取值使用imshow函数显示矩阵B
A = randn(100, 100);
B = zeros(size(A));
【方案1】
[X,Y] = find(A > 0.6);
For i = 1:length(X)
B(X(i),Y(i)) = 1;
End
【方案2】
B = zeros(size(A));
X = find(A > 0.6);
B(X) = 1;
【方案3】：
B = zeros(size(A));
X = logical(A > 0.6);
B(X) = 1;
事实上，【方案1】到【方案2】的改进在“再谈Matlab的多维数组问题”一文中已经提及过，但是没想到自己在矩阵输入的时候注意到了，但是在矩阵输出的时候就忘记了，看来程序是需要不断修改、优化的，技巧也是需要不断巩固的。至于【方案3】，是得益于matlab的warning提示。
然而，这并不是表示所有类似的地方都可以用logical代替find，当遇到循环次数与X有关时，用find会更有效，对此可以参考【例2】的【方案2】，这是用logical实现不到的。
【例2】任务描述：有一个四维矩阵A（大小为61*73*61*210），其中前三维表示一个包含大脑结构的立方体，最后一维表示大脑中每个点对应的一个长度为210的时间序列。另外有一个三维矩阵Mask（大小是61*73*61），B是二值的，其中1表示该点是前景点（大脑），0表示该点是背景点。任务是对A中属于前景点的时间序列进行EMD处理，从而判断该前景点是否属于激活区。
显然，这个问题要利用循环来完成。在本人写的“再谈Matlab的多维数组问题”一文中已经提到可以把多维数组转化为一维数组来处理，在这里，也要利用这个思想。而且，从下文可以看到，正是因为 “多维变一维”的出现，才令程序得到更进一步的提速。
首先利用reshape函数把四维矩阵A变成二维矩阵B，把三维矩阵Mask变成一维矩阵C：
B = reshape(A, 61*73*61, 210);
C = Mask(:);
t = 1:210;
【方案1】
iTotalVoxel = 61*73*61;
for k = 1:iTotalVoxel
if C(k) == 1
temp = B(k,: );
imf = emd(t, temp);
…
end
end
【方案2】
D = find(C);
iTotalBrainVoxel = length(D);
for k = 1: iTotalBrainVoxel
temp = B(D(k),: );
imf = emd(t, temp);
…
end
【方案1】明显是基于C语言的套路，而【方案2】则充分避免了matlab的弱点――循环，经过改进以后（“多维转一维”为此提供了保证，多维的话，恐怕要使用形如A(X(k),Y(k),Z(k),:)的形式了），由于循环次数的降低（大约降低为原来的1/3），故运行时间大致上减少了一半。
由此可见，在matlab中，想加快运行速度，不但要减少循环的层数，而且，还要减少循环的次数。
以下是对【例2】（实现大脑激活区检测）的不同实现的结果比较：（核心工作是对73000个前景点相应的时间序列进行EMD处理，生成10个左右的imf，其中抽取每个imf时平均迭代大概10次左右）
版本1：完全用matlab编写——运行时间大约是200分钟，也就是说，要在matlab做73000*100次循环，最头痛的是迭代时的需要产生样条包络，默认的spline函数耗时相当多；
版本2：用matcom转换成cpp文件后在Borland C++ Builder中运行，完全脱离matlab——由于spline函数耗时太多，因此转换前改用了折线包络，而非样条包络，运行时间为15分钟左右，不过这个结果是对仿真数据，非实际数据而言的，因为我仍未解决matrix.h和matlib.h不能共存的问题，故无法对实际数据进行测试，不过一般来说实际数据比仿真数据运行速度更慢；
版本3：把核心代码（基于样条包络的EMD算法）做成dll文件后在matlab中调用——整个程序，从数据输入到数据输出，只在一个地方使用了一层for循环（就是【例2】中【方案2】的循环），且结合上述优化方案，对于实际数据大概5分钟就能得出结果。
小结：要学好matlab，有效地使用matlab，一定要摆脱C＋＋的思想。能不用循环的地方，尽量不要使用（例如求极值点等这些算法基本上可以不使用循环便可实现），逐渐抛弃C＋＋的逐点运算的思想，多从矩阵的整体（或分块）上考虑。虽然matlab 2006a版对循环已经改进不少，但是，循环的确是造成程序运行速度降低的主要原因。matlab提供的远远不止在调用函数上的方便（例如在C++中编写fft、dwt等函数可能需要几十甚至几百行，而在matlab中只需要一两个语句），运行速度慢或许是没有使用好它，让它发挥出所长所致的。想matlab更高效地为你服务，那就需要不断修改、优化你的代码吧（我的程序编写大概用了一个星期，而修改、优化的时间就用了两个多星期，呵呵）。最后，套用某人的一句话来作结：与其抱怨matlab运行速度慢，不如先改进一下你的算法吧。
提高matlab运行速度和节省空间的一点心得（之三）
首先说说Matlab与其他语言的差异：例如对于C或者C++来说，只要算法的思想不变、采用的数据结构相同，不同人写出来的语句在效率上一般不会产生太大的差别。所以，对于C来说，程序的好坏一般由算法来决定。但是，在matlab中，同样的算法、同样的结构、同样的流程，如果采用的语句不一样，在效率上就会大大不同。所以，我认为，使用matlab比使用其他语言更加困难，也显得matlab更难以掌握。另外，由于matlab在存储管理上的不便，使得在同时提高时空两域的效率变得更加困难，特别是在空间上（因为在时间上matlab提供了profiler这个非常有用的工具，但是在空间上就没有）。当需要处理大量的数据时，精简时空两域的程序语句就尤为重要了。
空间上：
1. 建议使用A = logical(sparse(m,n))，不建议使用 A = sparse(false(m,n))，两者结果一样，但是后者生成m×n的临时矩阵，浪费空间，且当m、n很大时，后者不一定能申请成功；
2. 使用sparse几点注意：
a) 只能用在二维以下的矩阵上；
b) 由于matlab按照“先行后列”的方式读取数据（即先把第一列所有行读取完以后再读取第二列的各行），因此定义稀疏矩阵时，最好“行数>列数”，这样有利于寻址和空间的节省（自己试试a=sparse(10,5); whos a和b= sparse(5,10);whos b就知道了）；
c) 对大型矩阵用sparse非常有效（不但节省空间，而且加快速度，强烈推荐！这在动态申请数组空间的时候尤其方便，当然了，数组不是太大的时候也可以使用eval即字符串的方法），但对小型矩阵使用反而增加存储量（自己试试a=false(5,1); whos a和b=logical(sparse(5,1));whos b就知道了），相信这是由稀疏矩阵需要存储额外的信息引起的。
3. 尽量按照精度来选择数据的类型，例如，如果只需用到0－255之间的整数，则定义矩阵为uint8型就ok了，定义方式：A = zeros(10,10,’uint8’)；可以用intmin(‘uint8’)和intmax(‘uint8’)返回该种类型的最值。
时间上：
1. 在for循环中，清零操作用赋值语句 A = B，其中B是在for循环外的一个同A大小一样的全0阵，不要使用A(:) = 0；但这样会大大影响后面的逐点运算速度。这个问题要请教高手，那就是“个别语句的改动会引起其他语句的执行速度”。例如分别运行3万多次下面代码，但执行时间有较大差别：
iLen = length(find(alRegion)); % 0.58 s
if iLen >= iThreshold_2 % 0.05 s
…
end
和
iLen = sum(alRegion); % 0.37 s
if iLen >= iThreshold_2 % 0.40 s
…
end
2. Find函数较慢，可用logical函数代替，但是，当需要取得满足条件的下标时，就无法使用logical函数，这已经在我之前的帖子中提及过。不过，大家有没有想过，连find函数都可以进行优化的，方法是用“基本矩阵进行显式逻辑引用”来代替find。例如，假设矩阵A是一个61*73*20的三维逻辑矩阵，如果用下面的语句循环3万多次，需要的时间是13 s ：
B = find(A == true);
如果采用下面的方法，则只需要不到0.7 s：
首先定义一个索引矩阵：
a3iCubicIdx = uint32(1:iTotalVoxel); % uint32可以根据需要调整，这里省略了条件判断
a3iCubicIdx = reshape(a3iCubicIdx, [iVolLen, iVolWdh, iVolHgh]);
然后在循环中写以下代码：
a3iTemp = a3iCubicIdx(iXmin:iXMax,iYMin:iYMax,iZMin:iZMax);
B = aiTemp(A(iXmin:iXMax,iYMin:iYMax,iZMin:iZMax));
当然了，改进的前提是知道矩阵A的非零元（即值为true的元素）大致的分布，也就是能够求出iXmin:iXMax,iYMin:iYMax,iZMin:iZMax这个范围。现在终于明白并体会到cwit所说的“连num2str都优化过”的含义了。
3. 不断优化代码，例如corrcoef函数，matlab自带的corrcoef函数求整个矩阵所有列的相关系数，因为我只需要求出某一列跟其他各列的相关系数，所以参照corrcoef函数自己写了一个，不但把速度提了上去，而且还发现了：repmat(5，100，1)的速度并不比ones(100,1)*5 快，另外，别小看一个小矩阵的转置操作，当循环次数很大的时候，有没有转置就相差很远了。
4. 使用逻辑运算符&、| 时，两个操作对象最好是logical类型，否则速度会减慢。
5. 二维矩阵转置操作可以用以下三种方法进行，三者的效率基本一样（时空），如果遇到三维以上的矩阵要转置，用permute命令较为方便：
a) A = A’;
b) A = permute(A,[2,1]);
c) A = shiftdim(A,1);
6. “使用eval方式动态存储多个一维数组”比“使用二维数组动态存储多个一维数组”要快，即：eval(['A_', num2str(k),' = B;']);比 A(k,:) = B; 快，其中B是一个一维数组，k表示循环次数。注：并非所有B都进行存储，只存储满足某个条件的B，另外，对预申请空间A不成功，这是对结论的补充说明。值得注意的是，如果对B是一个稀疏的一维数组，则eval方式的优势荡然无存，当k增大时反而增加系统开销。
7. 当矩阵很大时，利用A(:,k+1:end)=[];去掉多余元素操作时会减慢程序的运行，因此，如果后续处理中没有用到这些多余元素，则没有必要使用这个语句，即不管就是了。
8. 当需要对很大的一个矩阵进行操作时，可以考虑使用循环来完成。例如corrcoef函数，如果处理的对象矩阵A是100*180时（即对100个列向量求它们两两之间的相关程度，假设需要的只是前面99个与第100个向量的相关系数，其他不需要用到），直接用corrcoef(A)会比较慢，这时候可以考虑把矩阵A分为5个部分，每个子块与第100个向量进行相关，这样速度会更快。
9. 局部比较、赋值比全局比较、赋值要快（呵呵，这是废话），假设A、B都是三维逻辑矩阵，如果只想对某个局部（例如X_1:X_2,Y_1:Y_2,Z_1:Z_2这个立方体）进行比较和赋值，则推荐使用B(X_1:X_2,Y_1:Y_2,Z_1:Z_2) = B(X_1:X_2,Y_1:Y_2,Z_1:Z_2) & ~A(X_1:X_2,Y_1:Y_2,Z_1:Z_2)，这比B(A) = false或B = B&~A速度上都要快不少。
后记：
1. 此心得原来打算在cwit给了我回复，指正文中错误之后再整理并发布的，不过cwit很忙，所以不知道要等到什么时候，因此先献丑了，有错误大家也帮忙指正一下。特别是happy、bainhome、jimin，还有Tla等高手。（注：之前写的版本与这个版本略有不同，主要是删除了“循环总次数固定无论使用多少个for循环速度不变”这句，因为考虑到矩阵维数会影响速度）
2. 上述技巧都是出自我在处理大型矩阵时候自己总结出来的个人心得，转载时请注明。
3. 实验室中由于只有我一个人比较关心matlab的运行速度和存储空间问题（尤其是速度问题），所以不像cwit那样，有一个团队可以互相讨论互相提高，因此，我的心得难免有错误或理解不当的地方，请大家见谅。
4. 错误的地方待cwit给我完整的回复后我会跟贴更正或补充。
发表评论