matlab代码优化（转来的，没法转，就copy过来了）

复制于http://blog.sciencenet.cn/blog-388372-687651.html

用过Matlab的人都知道，Matlab是一种解释性语言，存在计算速度慢的问题，为了提高程序的运行效率，matlab提供了多种实用工具及编码技巧。

 

1. 循环矢量化
Matlab是为矢量和矩阵操作而设计的，因此，可以通过矢量化方法加速M文件的运行。矢量化是指将for循环和while循环转换为等价的矢量或矩阵操作。下面给出一个循环的例子：
i=0;
for n = 0:0.1:1000
   i=i+1;
   y(i)=cos(n)；
end
那么我们可以矢量化为：
n= 0:0.1:1000;
y=cos(n);
我们可以用tic和toc函数来查看上述各代码运行的时间，采用for循环的程序0.39秒（具体时间和计算机配置有关），而矢量化后几乎耗时为0。

2. 给数组或矩阵预分配内存
   特别是使用大型数组或矩阵时，Matlab进行动态内存分配和取消时，可能会产生内存碎片，这将导致大量闲置内存产生，预分配可通过提前给大型数据结构预约足够空间来避免这个问题。

3. 用函数代替脚本文件
   因为每次调用MATLAB的脚本文件都需要将不必要的中间变量加载到内存中，每执行一次，就加载一次。函数在调用时被编译成了伪代码，只需要加载到内存一次。当多次调用同一个函数时会运行快一些。因此尽量多使用函数文件而少使用脚本文件，也是提高执行效率的一种方法。

4. 用Mex文件编写循环代码
   Matlab提供了与C和C++的接口，那么我们可以在用C或C++语言编写耗时的循环代码，然后通过接口程序在Matlab中转换成dll文件，这就是我们所要的Mex文件，通过这种方法可以极大地提高计算速率。

 

1. 尽量避免使用循环结构
MATLAB变量的基本类型是矩阵，当对矩阵的每个元素循环处理时，运算速度很慢。因此编程时应尽量把数组和矩阵看作一个整体来进行编程，而不是像其他的程序设计语言那样，使用循环结构对矩阵的元素循环进行处理。利用MATLAB提供的用于矢量化操作的函数，把循环矢量化，这样既可以提高编程效率，也可以提高程序的执行效率。下面给出一个循环的例子：
i=0;
for n = 0:0.1:100
         i=i+1;
       y(i)=cos(n)
end
上述程序段把数组中的每个元素都进行函数值计算，这样会耗费大量的运算时间，我们可以把数组看作一个整体来处理，计算函数值，可以修改这个程序段如下。
n = 0:0.1:100;
y = cos(n)
通过使用MATLAB专门提供的测试程序运行时间的函数，可以发现，把数组看作一个整体，进行操作后，执行效率提高约300倍。
另外，在必须使用多重循环的情况下，建议在循环的外环执行循环次数少的，内环执行循环次数多的，这样也可以显著提高程序执行速度。
2. 在使用数组或矩阵之前先定义维数
MATLAB中的变量在使用之前不需要明确地定义和指定维数。但当未预定义数组或矩阵的维数时，当需赋值的元素下标超出现有的维数时，MATLAB 就为该数组或矩阵扩维一次，这样就会大大降低程序的执行效率。因此，在使用数组或矩阵之前，预定义维数可以提高程序的执行效率。
3. 对矩阵元素使用下标或者索引操作
在MATLAB中，矩阵元素的引用可用两个下标来表示。例如：A(i,j) 表示矩阵的第i行第j列的元素；A(1:k,j)表示矩阵A的第j列的前k个元素；A(:,j) 表示矩阵的第j列的所有元素。求矩阵A的第j列元素的平均值的表达式为mean(A(:,j))。
4. 尽量多使用函数文件少使用脚本文件
因为每次调用MATLAB的脚本文件都需要将不必要的中间变量加载到内存中，每执行一次，就加载一次。函数在调用时被编译成了伪代码，只需要加载到内存一次。当多次调用同一个函数时会运行快一些。因此尽量多使用函数文件而少使用脚本文件，也是提高执行效率的一种方法。
5. 在必须使用循环时，可以考虑转换为C-MEX
当必须使用耗时的循环时，可以考虑将循环体中的语句转换为C-MEX。C-MEX是将M文件通过MATLAB的编译器转换为可执行文件，是按照 MEX 技术要求的格式编写相应的程序，通过编译连接，生成扩展名为.dll的动态链接库文件，可以在MATLAB环境下直接执行。这样，循环体中的语句在执行时不必每次都解释(interpret)。一般来说，C-MEX 文件的执行速度是相同功能的M文件执行速率的20~40倍。编写C-MEX不同于M文件，需要了解MATLAB C-MEX规范。幸运的是MATLAB提供了将M文件转换为C-MEX的工具。
6. 内存优化
MATLAB在进行复杂的运算时需要占用大量的内存。合理使用内存和提高内存的使用效率，可以加快运行速度，减少系统资源的占用。
7. 内存管理函数和命令
（1）Clear variablename：从内存中删除名称为variablename的变量。
（2）Clear all：从内存中删除所有的变量。
（3）Save：将指令的变量存入磁盘。
（4）Load：将save命令存入的变量载入内存。
（5）Quit：退出MATLAB，并释放所有分配的内存。
（6）Pack：把内存中的变量存入磁盘，再用内存中的连续空间载回这些变量。考虑到执行效率问题，不能在循环中使用。
8. 节约内存的方法
（1）避免生成大的中间变量，并删除不再需要的临时变量。
（2）当使用大的矩阵变量时，预先指定维数并分配好内存，避免每次临时扩充维数。
（3）当程序需要生成大量变量数据时，可以考虑定期将变量写到磁盘，然后清除这些变量。
当需要这些变量时，再重新从磁盘加载。
（4）当矩阵中数据极少时，将全矩阵转换为稀疏矩阵。

 

提高MATLAB程序效率的几点原则，这些都是俺在这两年中参加四次数模编写大量m程序总结的经验，加之网上很多英雄也是所见略同。
1.“计算向量、矩阵化，尽量减少for循环。”[/B]
因为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.“循环内大数组预先定义--预先分配空间”[/U]
这一点原则是极其重要的，以至于在编写m程序时编辑器会给出提示“'ver' might be growing inside a loop.Consider prealloacting for speed.”

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内部提供的函数”[/U]
因为matlab内部提供的函数绝对是各种问题的最优算法，那写程序都是他们大师级人物写出来的，程序应该说相当高效，有现成的为什么不用那！     这个原则就不用实际的程序测试了。

 

 

关于MATLAB程序提速的问题，可以参考网上很多智者的文章，都比较经典。也可以看看我的上一篇文章，和网上大部分帖子有点不同，我是以实际的测试程序作为依据对如何提高MATLAB程序速度进行介绍的。      这里我再补充几点大家需要注意的。下面是我在国内一个比较出名的论坛看到的关于m程序提速的帖子，开始还真以为他们谈论的都应该遵循。（尽信书不如无书）
帖子的一部分这样说道：“当要预分配一个非double型变量时使用repmat函数以加速，如将以下代码：
A = int8(zeros(100));
换成：
A = repmat(int8(0), 100, 100);”
凡事不能只凭自己的感觉，没有一点实际的例子，对于权威我们还要有挑战精神那，就不用说现在还不是经典的观点了。下面是我写的测试程序，我本来是想得到这位网友大哥的结果，但是实事不是我们想象的那么简单。

n = 100;
m = 1000;
for k=1:n
   tic
   A = int8(ones(m));
   t1(k) = toc;
   tic
   B = repmat(int8(1),m,m);
   t2(k) = toc;
end
plot(1:n,t1,'r',1:n,t2)
isequal(A,B)
可以看出下面的红线是t1，而且最后的一句返回1，说明两种形式返回的值完全一样。

由此我想说的是，不管是在我们做论文，还是写博客的时候，别直接从网上或者别人文章那儿找点知识定理之类的补充自己那苍白无力的文章。最好是自己动手编一下，“实践是检验真理的唯一标准”。

 

经过这一测试，我感觉有必要，也有责任对这个论坛上的一大批经典谈论加以测试。尽管这个结论是错误的但这还不足以证明论坛上的帖子都不是经典。
还有一点关于m程序提速的这样说到：“在必须使用多重循环时下，如果两个循环执行的次数不同，则在循环的外环执行循环次数少的，内环执行循环次数多的。这样可以显著提高速度。”

n=1000;
A = ones(1000)*13;
for k=1:n
   tic
   for i=1:10
       for j=1:1000
           A(i,j)=A(i,j)*15;
       end
   end
   t1(k)=toc;
   tic
   for i=1:1000
       for j=1:10
           A(i,j)=A(i,j)*16;
       end
   end
   t2(k)=toc;
end
t2(t1>10^9)=[];
t1(t1>10^9)=[];
t1(t2>10^9)=[];
t2(t2>10^9)=[];%去除外界因素影响所产生的寄点
plot(1:size(t1,2),t1,'r',1:size(t1,2),t2)

 

由这个时间图可以看出for循环的嵌套顺序对于速度是有影响的，虽然相对来说差别不是很大，但是毕竟论坛上的观点是正确的。至于他所说的“显著”二字就没必要加上了。
此论坛还有一些提速的观点列举如下：
“遵守Performance Acceleration的规则
    关于什么是“Performance Acceleration”请参阅matlab的帮助文件。我只简要的将
其规则总结如下7条：
1、只有使用以下数据类型，matlab才会对其加速：
logical,char,int8,uint8,int16,uint16,int32,uint32,double 而语句中如果使用了非以上的数据类型则不会加速，如：numeric,cell,structure,single,function handle,java classes,user classes,int64,uint64
2、matlab不会对超过三维的数组进行加速。
3、当使用for循环时，只有遵守以下规则才会被加速：a、for循环的范围只用标量值来表示；
b、for循环内部的每一条语句都要满足上面的两条规则，即只使用支持加速的数据类型，只使用三维以下的数组；c、循环内只调用了内建函数（build-in function）。
4、当使用if、elseif、while和switch时，其条件测试语句中只使用了标量值时，将加速运行。
5、不要在一行中写入多条操作，这样会减慢运行速度。即不要有这样的语句：
x = a.name; for k=1:10000, sin(A(k)), end;
6、当某条操作改变了原来变量的数据类型或形状（大小，维数）时将会减慢运行速度。
7、应该这样使用复常量x = 7 + 2i，而不应该这样使用：x = 7 + 2*i,后者会降低运行速度。”
“尽量用向量化的运算来代替循环操作。如将下面的程序：
i=0;
for t = 0:.01:10
   i = i+1;
   y(i) = sin(t);
end
替换为：
t = 0:.01:10;
y = sin(t);
速度将会大大加快。最常用的使用vectorizing技术的函数有：All、diff、ipermute、permute、reshape、ueeze、y、find、logical、prod、shiftdim、sub2ind、cumsum、ind2sub、ndgrid、repmat、sort、sum 等。”
“优先使用matlab内建函数，将耗时的循环编写进MEX-File中以获得加速。
b、使用Functions而不是Scripts 。”
“ 绝招：你也许觉得下面两条是屁话，但有时候它真的是解决问题的最好方法。
1、改用更有效的算法
2、采用Mex技术，或者利用matlab提供的工具将程序转化为C语言、Fortran语言。关于如何将M文件转化为C语言程序运行，可以参阅本版帖子：“总结：m文件转化为c/c++语言文件，VC编译”。 ”

 

 

除了m程序提速的问题，这里还列出了《MATLAB代码矢量化指南（译）》
一、基本技术 
----------------------------------------------------- 
1)MATLAB索引或引用(MATLAB Indexing or Referencing) 
在MATLAB中有三种基本方法可以选取一个矩阵的子阵。它们分别是 下标法，线性法和逻辑法(subscripted,linear, and logical)。 
如果你已经熟悉这个内容，请跳过本节

1.1)下标法 
非常简单，看几个例子就好。 
A = 6:12; 
A([3,5]) 
ans = 
8 10 
A([3:2:end]) 
ans = 
8 10 12

A = 
[11 14 17; 
12 15 18;
13 16 19]; 
A(2:3,2) 
ans = 
15 
16

1.2)线性法 
二维矩阵以列优先顺序可以线性展开，可以通过现行展开后的元素序号 
来访问元素。 
A = 
[11 14 17;
12 15 18;
13 16 19]; 
A(6) 
ans = 
16

A([3,1,8]) 
ans = 
13 11 18 
A([3;1;8]) 
ans = 
13 
11 
18

1.3)逻辑法 
用一个和原矩阵具有相同尺寸的0-1矩阵，可以索引元素。在某个位置上为1表示选取元素，否则不选。得到的结果是一个向量。 
A = 6:10; 
A(logical([0 0 1 0 1])) 
ans = 
8 10 
A = 
[1 2 
3 4]; 
B = [1 0 0 1]; 
A(logical(B)) 
ans = 
1 4 
----------------------------------------------------- 
2)数组操作和矩阵操作(Array Operations vs. Matrix Operations) 
对矩阵的元素一个一个孤立进行的操作称作数组操作；而把矩阵视为一个整体进行的运算则成为矩阵操作。MATLAB运算符*,/,\,^都是矩阵运算，而相应的数组操作则是.*, ./, .\, .^ 
A=[1 0 ;0 1]; 
B=[0 1 ;1 0]; 
A*B % 矩阵乘法 
ans = 
0 1 
1 0 
A.*B % A和B对应项相乘 
ans = 
0 0 
0 0 
------------------------------------------------------ 
3)布朗数组操作(Boolean Array Operations) 
对矩阵的比较运算是数组操作，也就是说，是对每个元素孤立进行的。因此其结果就不是一个“真”或者“假”，而是一堆“真假”。这个结果就是布朗数组。 
D = [-0.2 1.0 1.5 3.0 -1.0 4.2 3.14]; 
D >= 0 
ans = 
0 1 1 1 0 1 1 
如果想选出D中的正元素： 
D = D(D>0) 
D = 
1.0000 1.5000 3.0000 4.2000 3.1400 
除此之外，MATLAB运算中会出现NaN,Inf,-Inf。对它们的比较参见下例 
Inf==Inf返回真 
Inf<1返回假 
NaN==NaN返回假 
同时，可以用isinf，isnan判断，用法可以顾名思义。在比较两个矩阵大小时，矩阵必须具有相同的尺寸，否则会报错。这是 你用的上size和isequal,isequalwithequalnans(R13及以后)。 
------------------------------------------------------ 
4)从向量构建矩阵(Constructing Matrices from Vectors) 
在MATLAB中创建常数矩阵非常简单，大家经常使用的是： 
A = ones(5,5)*10 
但你是否知道，这个乘法是不必要的？ 
A = 10; 
A = A(ones(5,5)) 
A = 
10 10 10 10 10 
10 10 10 10 10 
10 10 10 10 10 
10 10 10 10 10 
10 10 10 10 10 
类似的例子还有： 
v = (1:5)'; 
n = 3; 
M = v(:,ones(n,1)) 
M =

1 1 1 
2 2 2 
3 3 3 
4 4 4 
5 5 5 
事实上，上述过程还有一种更加容易理解的实现方法： 
A = repmat(10,[5 5]); 
M = repmat([1:5]', [1,3]); 
其中repmat的含义是把一个矩阵重复平铺，生成较大矩阵。更多详细情况，参见函数repmat和meshgrid。 
----------------------------------------------------- 
5)相关函数列表(Utility Functions) 
ones 全1矩阵 
zeros 全0矩阵 
reshape 修改矩阵形状 
repmat 矩阵平铺 
meshgrid 3维plot需要用到的X－Y网格矩阵 
ndgrid n维plot需要用到的X－Y－Z...网格矩阵 
filter 一维数字滤波器，当数组元素前后相关时特别有用。 
cumsum 数组元素的逐步累计 
cumprod 数组元素的逐步累计 
eye 单位矩阵 
diag 生成对角矩阵或者求矩阵对角线 
spdiags 稀疏对角矩阵 
gallery 不同类型矩阵库 
pascal Pascal 矩阵 
hankel Hankel 矩阵 
toeplitz Toeplitz 矩阵

==========================================================
二、扩充的例子 
------------------------------------------------------ 
6)作用于两个向量的矩阵函数 
假设我们要计算两个变量的函数F 
F(x,y) = x*exp(-x^2 - y^2) 
我们有一系列x值，保存在x向量中，同时我们还有一系列y值。我们要对向量x上的每个点和向量y上的每个点计算F值。换句话说，我们要计算对于给定向量x和y的所确定的网格上的F值。 
使用meshgrid，我们可以复制x和y来建立合适的输入向量。然后 
可以使用第2节中的方法来计算这个函数。 
x = (-2:.2:2); 
y = (-1.5:.2:1.5)'; 
[X,Y] = meshgrid(x, y); 
F = X .* exp(-X.^2 - Y.^2); 
如果函数F具有某些性质，你甚至可以不用meshgrid，比如 
F(x,y) = x*y ，则可以直接用向量外积 
x = (-2:2); 
y = (-1.5:.5:1.5); 
x'*y 
在用两个向量建立矩阵时，在有些情况下，稀疏矩阵可以更加有 
效地利用存储空间，并实现有效的算法。我们将在第8节中以一个 
实例来进行更详细地讨论. 
-------------------------------------------------------- 
7)排序、设置和计数(Ordering, Setting, and Counting Operations) 
在迄今为止讨论过的例子中，对向量中一个元素的计算都是独立于同一向量的其他元素的。但是，在许多应用中，你要做的计算则可能与其它元素密切相关。例如，假设你用一个向量x来表示一 个集合。不观察向量的其他元素，你并不知道某个元素是不是一个冗余元素，并应该被去掉。如何在不使用循环语句的情况下删除冗余元素，至少在现在，并不是一个明显可以解决的问题。 
解决这类问题需要相当的智巧。以下介绍一些可用的基本工具

max 最大元素 
min 最小元素 
sort 递增排序 
unique 寻找集合中互异元素（去掉相同元素） 
diff 差分运算符[X(2) - X(1), X(3) - X(2), ... X(n) - X(n-1)] 
find 查找非零、非NaN元素的索引值 
union 集合并 
intersect 集合交 
setdiff 集合差 
setxor 集合异或

继续我们的实例，消除向量中的多余元素。注意：一旦向量排序后，任何多余的元素就是相邻的了。同时，在任何相等的相邻元素在向量diff运算时变为零。这是我们能够应用以下策略达到目的。我们现在在已排序向量中，选取那些差分非零的元素。

% 初次尝试。不太正确! 
x = sort(x(:)); 
difference = diff(x); 
y = x(difference~=0);

这离正确结果很近了，但是我们忘了diff函数返回向量的元素个数比输入向量少1。在我们的初次尝试中，没有考虑到最后一个元素也可能是相异的。为了解决这个问题，我们可以在进行差分之前给向量x加入一个元素，并且使得它与以前的元素一定不同。一种实现的方法是增加一个NaN。

% 最终的版本。 
x = sort(x(:)); 
difference = diff([x;NaN]); 
y = x(difference~=0);

我们使用(:)运算来保证x是一个向量。我们使用~=0运算，而不用find函数，因为find函数不返回NaN元素的索引值，而我们操作中差分的最后元素一定是NaN。这一实例还有另一种实现方式：

y=unique(x);

后者当然很简单，但是前者作为一个练习并非无用，它是为了练习使用矢量化技术，并示范如何编写你自己的高效代码。此外，前者还有一个作用：Unique函数提供了一些超出我们要求的额外功能，这可能降低代码的执行速度。 
     假设我们不只是要返回集合x，而且要知道在原始的矩阵里每个相异元素出现了多少个“复本”。一旦我们对x排序并进行了差分，我们可以用find来确定差分变化的位置。再将这个变化位置进行差分，就可以得到复本的数目。这就是"diff of find of diff"的技巧。基于以上的讨论，我们有：

% Find the redundancy in a vector x 
x = sort(x(:)); 
difference = diff([x;max(x)+1]); 
count = diff(find([1;difference])); 
y = x(find(difference)); 
plot(y,count)

这个图画出了x中每个相异元素出现的复本数。注意，在这里我们避开了NaN，因为find不返回NaN元素的索引值。但是，作为特例，NaN和Inf 的复本数可以容易地计算出来：

count_nans = sum(isnan(x(:))); 
count_infs = sum(isinf(x(:)));

另一个用于求和或者计数运算的矢量化技巧是用类似建立稀疏矩阵的方法实现的。这还将在第9节中作更加详细的讨论. 
------------------------------------------------------- 
8)稀疏矩阵结构(Sparse Matrix Structures)

在某些情况下，你可以使用稀疏矩阵来增加计算的效率。如果你构造一个大的中间矩阵，通常矢量化更加容易。在某些情况下，你可以充分利用稀疏矩阵结构来矢量化代码，而对于这个中间矩阵不需要大的存储空间。 
假设在上一个例子中，你事先知道集合y的域是整数的子集， 
{k+1,k+2,...k+n}；即， 
y = (1:n) + k

例如，这样的数据可能代表一个调色板的索引值。然后，你就可以对集合中每个元素的出现进行计数（构建色彩直方图？译者）。这是对上一节中"diff of find of diff"技巧的一种变形。现在让我们来构造一个大的m x n矩阵A，这里m是原始x向量中的元素数， n是集合y中的元素数。 
A(i,j) = 1 if x(i) = y(j) 
0 otherwise

回想一下第3节和第4节，你可能认为我们需要从x和y来构造矩阵A。如果当然可以，但要消耗许多存储空间。我们可以做得更好，因为我们知道，矩阵A中的多数元素为0，x中的每个元素对应的行上只有一个值为1。 
以下就是构造矩阵的方法（注意到y(j) = k+j，根据以上的公式）： 
x = sort(x(:)); 
A = sparse(1:length(x), x+k, 1, length(x), n);

现在我们对A的列进行求和，得到出现次数。 
count = sum(A); 
在这种情况下，我们不必明确地形成排序向量y，因为我们事先知道 
y = 1:n + k.

这里的关键是使用数据，（也就是说，用x控制矩阵A的结构）。由于x在一个已知范围内取整数值，我们可以更加有效地构造矩阵。 假设你要给一个很大矩阵的每一列乘以相同的向量。使用稀疏矩阵，不仅可以节省空间，并且要比在第5节介绍的方法更加快速. 下面是它的工作方式： 
F = rand(1024,1024); 
x = rand(1024,1); 
% 对F的所有行进行点型乘法. 
Y = F * diag(sparse(x)); 
% 对F的所有列进行点型乘法. 
Y = diag(sparse(x)) * F;

我们充分利用矩阵乘法算符来执行大规模运算，并使用稀疏矩阵以防止临时变量变得太大。 
-------------------------------------------------------- 
9)附加的例子(Additional Examples) 
下面的例子使用一些在本技术手册中讨论过的方法，以及其它一些相关方法。请尝试使用tic 和toc（或t=cputime和cputime-t），看一下速度加快的效果。 
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
用于计算数组的双重for循环。 
使用的工具：数组乘法 
优化前： 
A = magic(100); 
B = pascal(100); 
for j = 1:100 
  for k = 1:100; 
      X(j,k) = sqrt(A(j,k)) * (B(j,k) - 1); 
  end 
end 
优化后： 
A = magic(100); 
B = pascal(100); 
X = sqrt(A).*(B-1); 
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
用一个循环建立一个向量，其元素依赖于前一个元素使用的工具：FILTER, CUMSUM, CUMPROD 
优化前： 
A = 1; 
L = 1000; 
for i = 1:L 
  A(i+1) = 2*A(i)+1; 
end 
优化后： 
L = 1000; 
A = filter([1],[1 -2],ones(1,L+1)); 
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
如果你的向量构造只使用加法或乘法，你可使用cumsum或cumprod函数。 
优化前： 
n=10000; 
V_B=100*ones(1,n); 
V_B2=100*ones(1,n); 
ScaleFactor=rand(1,n-1); 
for i = 2:n 
   V_B(i) = V_B(i-1)*(1+ScaleFactor(i-1)); 
end 
for i=2:n 
   V_B2(i) = V_B2(i-1)+3; 
end 
优化后： 
n=10000; 
V_A=100*ones(1,n); 
V_A2 = 100*ones(1,n); 
ScaleFactor=rand(1,n-1); 
V_A=cumprod([100 1+ScaleFactor]); 
V_A2=cumsum([100 3*ones(1,n-1)]); 
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
向量累加，每5个元素进行一次： 
工具：CUMSUM , 向量索引 
优化前： 
% Use an arbitrary vector, x 
x = 1:10000; 
y = []; 
for n = 5:5:length(x) 
   y = [y sum(x(1:n))]; 
end 
优化后(使用预分配)： 
x = 1:10000; 
ylength = (length(x) - mod(length(x),5))/5; 
% Avoid using ZEROS command during preallocation 
y(1:ylength) = 0; 
for n = 5:5:length(x) 
   y(n/5) = sum(x(1:n)); 
end 
优化后(使用矢量化,不再需要预分配)： 
x = 1:10000; 
cums = cumsum(x); 
y = cums(5:5:length(x)); 
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
操作一个向量，当某个元素的后继元素为0时，重复这个元素： 
工具：FIND, CUMSUM, DIFF 
任务：我们要操作一个由非零数值和零组成的向量，要求把零替换成为它前面的非零数值。例如，我们要转换下面的向量： 
a=2; b=3; c=5; d=15; e=11; 
x = [a 0 0 0 b 0 0 c 0 0 0 0 d 0 e 0 0 0 0 0]; 
为： 
x = [a a a a b b b c c c c c d d e e e e e e]; 
解(diff和cumsum是反函数)： 
valind = find(x); 
x(valind(2:end)) = diff(x(valind)); 
x = cumsum(x); 
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
将向量的元素累加到特定位置上 
工具：SPARSE 
优化前： 
% The values we are summing at designated indices 
values = [20 15 45 50 75 10 15 15 35 40 10]; 
% The indices associated with the values are summed cumulatively 
indices = [2 4 4 1 3 4 2 1 3 3 1]; 
totals = zeros(max(indices),1); 
for n = 1:length(indices) 
   totals(indices(n)) = totals(indices(n)) + values(n); 
end 
优化后： 
indices = [2 4 4 1 3 4 2 1 3 3 1]; 
totals = full(sparse(indices,1,values));

注意：这一方法开辟了稀疏矩阵的新用途。在使用sparse命令创建稀疏矩阵时，它是对分配到同一个索引的所有值求和，而不是替代已有的数值。这称为"向量累加"，是MATLAB处理稀疏矩阵的方式。

 

关于MATLAB的效率问题，很多文章，包括我之前写的一些，主要集中在使用向量化以及相关的问题上。但是，最近我在实验时对代码进行profile的过程中，发现在新版本的MATLAB下，for-loop已经得到了极大优化，而效率的瓶颈更多是在函数调用和索引访问的过程中。
由于MATLAB特有的解释过程，不同方式的函数调用和元素索引，其效率差别巨大。不恰当的使用方式可能在本来不起眼的地方带来严重的开销，甚至可能使你的代码的运行时间增加上千倍（这就不是多买几台服务器或者增加计算节点能解决的了，呵呵）。
下面通过一些简单例子说明问题。（实验选在装有Windows Vista的一台普通的台式机(Core2 Duo 2.33GHz + 4GB Ram)进行，相比于计算集群, 这可能和大部分朋友的环境更相似一些。实验过程是对某一个过程实施多次的整体进行计时，然后得到每次过程的平均时间，以减少计时误差带来的影响。多次实验，在均值附近正负20%的范围内的置信度高于95%。为了避免算上首次运行时花在预编译上的时间，在开始计时前都进行充分的“热身”运行。） 
函数调用的效率
一个非常简单的例子，把向量中每个元素加1。（当然这个例子根本不需要调函数，但是，用它主要是为了减少函数执行本身的时间，突出函数解析和调用的过程。）
作为baseline，先看看最直接的实现
% input u: u is a 1000000 x 1 vectorv = u + 1; 
这个过程平均需要0.00105 sec。
而使用长期被要求尽量避免的for-loop 
n = numel(u);% v = zeros(n, 1) has been pre-allocated.for i = 1 : n    v(i) = u(i) + 1;end
所需的平均时间大概是0.00110 sec。从统计意义上说，和vectorization已经没有显著差别。无论是for-loop或者vectorization，每秒平均进行约10亿次“索引-读取-加法-写入”的过程，计算资源应该得到了比较充分的利用。
要是这个过程使用了函数调用呢？MATLAB里面支持很多种函数调用方式，主要的有m-function, function handle, anonymous function, inline, 和feval，而feval的主参数可以是字符串名字，function handle, anonymous function或者inline。
用m-function，就是专门定义一个函数
function y = fm(x)    y = x + 1;
在调用时 
for i = 1 : n    v(i) = fm(u(i));end
function handle就是用@来把一个function赋到一个变量上，类似于C/C++的函数指针，或者C#里面的delegate的作用
fh = @fm;for i = 1 : n    v(i) = fh(u(i));end
anonymous function是一种便捷的语法来构造简单的函数，类似于LISP, Python的lambda表达式
fa = @(x) x + 1;for i = 1 : n    v(i) = fa(u(i));end
inline function是一种传统的通过表达式字符串构造函数的过程
fi = inline('x + 1', 'x');for i = 1 : n    v(i) = fi(u(i));end
feval的好处在于可以以字符串方式指定名字来调用函数，当然它也可以接受别的参数。
v(i) = feval('fm', u(i));v(i) = feval(fh, u(i));v(i) = feval(fa, u(i)); 
对于100万次调用（包含for-loop本身的开销，函数解析(resolution)，压栈，执行加法，退栈，把返回值赋给接收变量），不同的方式的时间差别很大： m-function0.385 secfunction handle0.615 secanonymous function0.635 secinline function166.00 secfeval('fm', u(i))8.328 secfeval(fh, u(i))0.618 secfeval(fa, u(i))0.652 secfeval(@fm, u(i))2.788 secfeval(@fa, u(i))34.689 sec
从这里面，我们可以看到几个有意思的现象： 

• 首先，调用自定义函数的开销远远高于一个简单的计算过程。直接写 u(i) = v(i) + 1 只需要  0.0011 秒左右，而写u(i) = fm(v(i)) 则需要0.385秒，时间多了几百倍，而它们干的是同一件事情。这说明了，函数调用的开销远远大于for-loop自己的开销和简单计算过程——在不同情况可能有点差别，一般而言，一次普通的函数调用花费的时间相当于进行了几百次到一两千次双精度浮点加法。

• 使用function handle和anonymous function的开销比使用普通的m-函数要高一些。这可能是由于涉及到句柄解析的时间，而普通的函数在第一次运行前已经在内存中进行预编译。

• inline function的运行时间则令人难以接受了，竟然需要一百多秒（是普通函数调用的四百多倍，是直接计算的十几万倍）。这是因为matlab是在每次运行时临时对字符串表达式（或者它的某种不太高效的内部表达）进行parse。

• feval(fh, u(i))和fh(u(i))，feval(fa, u(i))和fa(u(i))的运行时间很接近，表明feval在接受句柄为主参数时本身开销很小。但是，surprising的是 for i = 1 : n    v(i) = feval(@fm, u(i));end比起fh = @fm;for i = 1 : n    v(i) = feval(fh, u(i));end慢了4.5倍 (前者0.618秒，后者2.788秒)。
  for i = 1 : n    v(i) = feval(@(x) x + 1, u(i));end比起fa = @(x) x + 1;for i = 1 : n    v(i) = feval(fa, u(i));end竟然慢了53倍(前者0.652秒，后者34.689秒）。
  由于在MATLAB的内部实现中，function handle的解析是在赋值过程中进行的，所以预先用一个变量把句柄接下，其效果就是预先完成了句柄解析，而如果直接把@fm或者@(x) x + 1写在参数列上，虽然写法简洁一些，但是解析过程是把参数被赋值到所调函数的局部变量时才进行，每调用一次解析一次，造成了巨大的开销。

• feval使用字符串作为函数名字时，所耗时间比传入句柄大，因为这涉及到对函数进行搜索的时间（当然这个搜索是在一个索引好的cache里面进行（除了第一次），而不是在所有path指定的路径中搜索。）

  在2007年以后，MATLAB推出了arrayfun函数，上面的for-loop可以写为 v = arrayfun(fh, u)
  这平均需要4.48 sec，这比起for-loop(需时0.615 sec)还慢了7倍多。这个看上去“消除了for-loop"的函数，由于其内部设计的原因，未必能带来效率上的正效果。 
  元素和域的访问
  除了函数调用，数据的访问方式对于效率也有很大影响。MATLAB主要支持下面一些形式的访问：
  

• array-index  A(i):  

• cell-index:  C{i};  

• struct field:  S.fieldname

• struct field (dynamic):  S.('fieldname')

这里主要探索单个元素或者域的访问（当然，MATLAB也支持对于子数组的非常灵活整体索引）。
对于一百万次访问的平均时间
A(i) for a numeric array0.0052 secC{i} for a cell array0.2568 secstruct field0.0045 secstruct field (with dynamic name)1.0394 sec
我们可以看到MATLAB对于单个数组元素或者静态的struct field的访问，可以达到不错的速度，在主流台式机约每秒2亿次（连同for-loop的开销）。而cell array的访问则明显缓慢，约每秒400万次（慢了50倍）。MATLAB还支持灵活的使用字符串来指定要访问域的语法（动态名字），但是，是以巨大的开销为代价的，比起静态的访问慢了200倍以上。
关于Object-oriented Programming
MATLAB在新的版本中（尤其是2008版），对于面向对象的编程提供了强大的支持。在2008a中，它对于OO的支持已经不亚于python等的高级脚本语言。但是，我在实验中看到，虽然在语法上提供了全面的支持，但是matlab里面面向对象的效率很低，开销巨大。这里仅举几个例子。

• object中的property的访问速度是3500万次，比struct field慢了6-8倍。MATLAB提供了一种叫做dependent property的属性，非常灵活，但是，效率更低，平均每秒访问速度竟然低至2.6万次（这种速度基本甚至难以用于中小规模的应用中）。

• object中method调用的效率也明显低于普通函数的调用，对于instance method，每百万次调用，平均需时5.8秒，而对于static method，每百万次调用需时25.8秒。这相当于每次调用都需要临时解析的速度，而matlab的类方法解析的效率目前还明显偏低。

• MATLAB中可以通过改写subsref和subsasgn的方法，对于对象的索引和域的访问进行非常灵活的改造，可以通过它创建类似于数组的对象。但是，一个符合要求的subsref的行为比较复杂。在一个提供了subsref的对象中，大部分行为都需要subsref进行调度，而默认的较优的调度方式将被掩盖。在一个提供了subsref的类中（使用一种最快捷的实现），object property的平均访问速度竟然降到1万5千次每秒。

建议[/U]
根据上面的分析，对于撰写高效MATLAB代码，我有下面一些建议：

 

1. 虽然for-loop的速度有了很大改善，vectorization（向量化）仍旧是改善效率的重要途径，尤其是在能把运算改写成矩阵乘法的情况下，改善尤为显著。

2. 在不少情况下，for-loop本身已经不构成太大问题，尤其是当循环体本身需要较多的计算的时候。这个时候，改善概率的关键在于改善循环体本身而不是去掉for-loop。

3. MATLAB的函数调用过程（非built-in function）有显著开销，因此，在效率要求较高的代码中，应该尽可能采用扁平的调用结构，也就是在保持代码清晰和可维护的情况下，尽量直接写表达式和利用built-in function，避免不必要的自定义函数调用过程。在次数很多的循环体内（包括在cellfun, arrayfun等实际上蕴含循环的函数）形成长调用链，会带来很大的开销。

4. 在调用函数时，首选built-in function，然后是普通的m-file函数，然后才是function handle或者anonymous function。在使用function handle或者anonymous function作为参数传递时，如果该函数被调用多次，最好先用一个变量接住，再传入该变量。这样，可以有效避免重复的解析过程。

5. 在可能的情况下，使用numeric array或者struct array，它们的效率大幅度高于cell array（几十倍甚至更多）。对于struct，尽可能使用普通的域（字段，field）访问方式，在非效率关键，执行次数较少，而灵活性要求较高的代码中，可以考虑使用动态名称的域访问。

6. 虽然object-oriented从软件工程的角度更为优胜，而且object的使用很多时候很方便，但是MATLAB目前对于OO的实现效率很低，在效率关键的代码中应该慎用objects。

7. 如果需要设计类，应该尽可能采用普通的property，而避免灵活但是效率很低的dependent property。如非确实必要，避免重载subsref和subsasgn函数，因为这会全面接管对于object的接口调用，往往会带来非常巨大的开销（成千上万倍的减慢），甚至使得本来几乎不是问题的代码成为性能瓶颈。