一点MATLAB程序加速技巧

 
 
转自：http://www.ilovematlab.cn/thread-264312-1-1.html

很多人都嚷嚷着避免循环来提高程序速度，的确，某些时候，循环是影响程序速度的因素，但同时还别忘了还有判断流（if），以及各种复杂的数据结构（cell,struct,class）。除此之外还有没有其他因素呢？对了，还有程序调用、递归。

偶然看到http://undocumentedmatlab.com/上的几篇关于编程的文后，让我有了新的发现——上面都是一些常见，并且容易被大家发现的地方，然而在编程过程中还有很多被大家忽视的细节，这些细节往往使得你的程序开销成倍增长，而你却毫无察觉。这就是——变量的操作。

下面，我们来说说变量存亡、变量赋值、size变化对程序时间开销的影响。

1. 变量存亡与个数对于时间的影响
比如，a=10;b=a;与a=10;b=10; 显然是前者高效。
又如某段程序
     
     

      
      
 
       
(calculate value=...)
if value>7
cal=value+8;
(some operations on cal);
end
...
 
      
      
      复制代码
     
     
更高效的方法是这样的：
     
     

      
      
 
       
(calculate value=...)
if value>7
value=value+8;
(some operations on value);
end
...
 
      
      
      复制代码
     
     
对于变量，尽量使用同名或者等价操作(如果没有前后历时数据需要参与计算的话)，因为matlab内部的“copy-on-write”机制免去了重复不必要的开辟空间与重新读写操作，节约了时间。

关于clear，其实matlab内部有种机制即计算后很多被编译的结果是暂存在内存的，如果使用了clear(清空普通变量)或clear
 all(清空所有包括内存中的global,persistent变量，以及一些被记住的m函数,随机发生器状态)，有些时候，时间是大大增加了的。确保确实不用某些变量以及函数之后，再进行清空操作，这样才会节省时间，否则适得其反。

2. 数组size变化对于时间的影响
一般来说，size动态变化本身就是一个影响程序性的因素，应该极力避免。因为如果程序出现一些意想不到的错误，调试的时候极难发现。但是有时没法避免，那么可以考虑使用预分配的方法（关于这个，下文详细有论述）。当然，有些时候根本没法知道上限，因此无法预分配。那么这个时候动态扩大数组就很容易增加程序开销（mlint中经常会有提示），如何才能尽最大可能减少这种开销呢？下面看例子——
     
     

      
      
 
       
% 下列语句运行于 Matlab 7.12 (R2011a):
% #1: 直接赋值到特定的超界下标
data=[]; tic, for idx=1:100000; data(idx)=1; end, toc
   => Elapsed time is 0.075440 seconds. 

% #2: 单次超界赋值
data=[]; tic, for idx=1:100000; data(end+1)=1; end, toc
   => Elapsed time is 0.241466 seconds.    % 3 times slower

% #3: 用cat方法自动扩大数组
data=[]; tic, for idx=1:100000; data=[data,1]; end, toc
   => Elapsed time is 22.897688 seconds.   % 300 times slower!!!
 
      
      
      复制代码
     
     
明显，作为我们常常习惯使用的第三种方式，恰恰是最耗时的方式。最佳的做法是第一种（如果知道上限的话）——直接用下标赋值；如果无法估计上限，那么就用第二种。

3.预分配/动态增加
预分配/动态增加是一种c/c++中经常使用的有效的策略。对于变量可变大小的情况，这会节约运行时间，减少出错的可能性。
对于可估计size上限的时候，使用预分配比如a=cell(5,8);a=zeros(3,1);等等，同时在程序计算完毕时候，将多余空间还回去。但是对于不可预知空间上限，那么可以使用固定增量法，比如一次增加500或者50等等。

4.初始化/赋值
这是最容易忽视的地方了。很多人以为这里并不会出现多少时间损耗。诚然，在计算量不太大的时候，这种方法产生的时间损耗可以不计，但是，一旦计算规模大了，这种损耗是非常可观的。比如下面的例子：
     
     

      
      
 
       
% 运行环境： Matlab 7.12 (R2011a):
matrix = magic(3);

% #1: 直接赋值-快速并且是线性时间代价
data=[]; tic, for idx=1:10000; data(:,(idx*3-2):(idx*3))=matrix; end, toc
   => Elapsed time is 0.969262 seconds.
data=[]; tic, for idx=1:100000; data(:,(idx*3-2):(idx*3))=matrix; end, toc
   => Elapsed time is 9.558555 seconds.

% #2: 用cat方法自增 – 很慢, 时间代价是二次的
data=[]; tic, for idx=1:10000; data=[data,matrix]; end, toc
   => Elapsed time is 2.666223 seconds.
data=[]; tic, for idx=1:100000; data=[data,matrix]; end, toc
   => Elapsed time is 356.567582 seconds.
 
      
      
      复制代码
     
     
赋值最好是整体赋值（即“一次性到位”），尽量避免单个元素的操作。

对于变量初始化：
     
     

      
      
 
       
% MathWorks 建议的编程方式
clear data1, tic, data1 = zeros(1000,3000); toc
   => Elapsed time is 0.016907 seconds.

% 一种更快的方案 - 500 times faster!
clear data1, tic, data1(1000,3000) = 0; toc
   => Elapsed time is 0.000034 seconds.
 
      
      
      复制代码
     
     
由于zeros函数的实现机制其实并不是等价于for i=1:1000*3000,
 data1(i)=0; end，而是一个内存块整体分配，且共享一个值(0)。所以比单个元素赋值要快。但为何data1(1000,3000) = 0;这种方案更加快呢？或许这相当于预先告知最大分配空间，然后告知最后一个值为0，剩余的元素被程序默认为0。

结合上面关于初始化和赋值的加速技巧，我们可以扩展演变得到下面的技巧：
     
     

      
      
 
       
scalar = pi;  % for example...

data = scalar(ones(1000,3000));           % 方案A: 87.680 ms
data(1:1000,1:3000) = scalar;             % 方案 B: 28.646 ms
data = repmat(scalar,1000,3000);          % 方案 C: 17.250 ms
data = scalar + zeros(1000,3000);         % 方案 D: 17.168 ms
data(1000,3000) = 0; data = data+scalar;  % 方案 E: 16.334 ms
 
      
      
      复制代码
     
     
从上面的例子可以看出，初始化或者赋值的最佳方案是D和E，而不是大多数人所习惯的A和B。D和E中所使用的便是标量和数组的运算（自动扩张大小）。


结论：
+如果变量size上界可估计，则使用预分配；如果变量size按固定大小或者百分比递增，则考虑使用固定增量法。
+直接超界赋值往往更加高效。
+使用cell数组储存多个变动大小的数据，然后使用cell2mat函数将结果转换为普通数组。
+复用已存在的数组。
+尽量使用标量与数组（矩阵）的运算，来代替数组运算。

更新：
1.结构体变量的初始化与赋值
     
     

      
      
 
       
tic;S(250,500) = struct( 'a', [], 'b', [], 'c',[]);toc
tic;W(250,500).a = [];   W(250,500).b = [];   W(250,500).c = [];toc
Elapsed time is 0.000379 seconds.
Elapsed time is 0.005867 seconds.
 
      
      
      复制代码
     
     
可见，单个元素赋值与整体赋值相比，效率是非常低下的。普通数组、cell数组也是一样。

对于元素重复的情况，有下面的例子：
     
     

      
      
 
       
tic;S(1:1250,1:1500) = struct('a', 20, 'b', 30, 'c', 40);toc
tic;W = repmat(struct('a', 20, 'b', 30, 'c', 40), [1250 1500]);toc
tic;X=struct('a', 20, 'b', 30, 'c', 40);X1=X(ones(1250,1500));toc
Elapsed time is 0.264818 seconds.
Elapsed time is 0.106745 seconds.
Elapsed time is 0.166025 seconds.
 
      
      
      复制代码
     
     
因此，对于构造具有重复元素的结构体，显然repmat方法更快，其次是使用ones索引方法。

但是对于普通数组，情况则有所不同，比如：
     
     

      
      
 
       
tic;S(1:1250,1:1500) = 9;toc
tic;W = repmat(9, [1250 1500]);toc
tic;X=9;X1=X(ones(1250,1500));toc
tic;Y=9;Y1=Y*ones(1250,1500);toc
Elapsed time is 0.016590 seconds.
Elapsed time is 0.015504 seconds.
Elapsed time is 0.097858 seconds.
Elapsed time is 0.015138 seconds.
 
      
      
      复制代码
     
     
显然是标量乘法和repmat方法更加快。

2.欠维度运算
如果数组有同样运算分别作用各行/列，那么最好用矩阵运算（比如矩阵乘法/除法/幂等），或者采用bsxfun，而不是数组运算。以乘法举例如下：
     
     

      
      
 
       
a=rand(2000,1);b=rand(1,4000);
tic;c=a*b;toc
tic;d=repmat(a,1,4000).*repmat(b,2000,1);toc
tic;f=bsxfun(@times,a,b);toc
Elapsed time is 0.111471 seconds.
Elapsed time is 0.280428 seconds.
Elapsed time is 0.072298 seconds.
 
      
      
      复制代码
     
     

3.稀疏矩阵
对于大量存在0元素的矩阵来说，使用稀疏矩阵(sparse)能更有效节约储存空间与计算时间。MATLAB中的稀疏矩阵与普通矩阵一样，支持所有内建的算术、索引、逻辑操作，还有数组运算、矩阵运算。这些操作运算能作用于稀疏矩阵之间（返回稀疏矩阵），稀疏矩阵与普通矩阵之间（大多数情况返回普通矩阵）。

注意：
1）稀疏矩阵不支持下列函数：
特征值和奇异值相关函数、矩阵分解相关函数、生成随机数和复数的函数、特殊数学函数（比如贝塞尔、伽马函数等）、按位运算函数(biit-wise)。
2）稀疏矩阵支持的最大元素个数为2^31-1。
3）尽量使用(i,j)索引。

稀疏矩阵的构造/初始化
如果按元素（elementwise）赋值初始化一个稀疏矩阵，会造成索引操作的巨大开销（以下例子均来自MATLAB帮助文件）：
     
     

      
      
 
       
S1 = spalloc(1000,1000,100000);
tic;
for n = 1:100000
    i = ceil(1000*rand(1,1));
    j = ceil(1000*rand(1,1));
    S1(i,j) = rand(1,1);
end
toc
Elapsed time is 26.281000 seconds.
 
      
      
      复制代码
     
     
但是，构造索引向量和相应的值进行初始化，那会非常快：
     
     

      
      
 
       
i = ceil(1000*rand(100000,1)); %i-subindex
j = ceil(1000*rand(100000,1)); %j-subindex
v = zeros(size(i));
for n = 1:100000
    v(n) = rand(1,1); %value
end

tic;
S2 = sparse(i,j,v,1000,1000);
toc
Elapsed time is 0.078000 seconds.
 
      
      
      复制代码
     
     
稀疏矩阵的操作
由于稀疏矩阵也是按列储存，因此列存取比行要快：
     
     

      
      
 
       
% 按行添加
S = sparse(10000,10000,1);
tic; 
for n = 1:1000
   A = S(100,:) + S(200,:);
end; 
toc
Elapsed time is 1.208162 seconds.

% 按列添加
S = sparse(10000,10000,1);
tic;
for n = 1:1000
   B = S(:,100) + S(:,200);
end;
toc
Elapsed time is 0.088747 seconds.
 
      
      
      复制代码
     
     
因此，如果可能，尽量使用列运算，然后转置结果即可：
     
     

      
      
 
       
S = sparse(10000,10000,1);
tic; 
for n = 1:1000
   A = S(100,:)' + S(200,:)';
   A = A';
end; 
toc
Elapsed time is 0.597142 seconds.
 
      
      
      复制代码
     
     
补充：转置所需的时间可以忽略不计，但是要注意的是，行数巨大的稀疏矩阵在转置时会消耗相当大的内存，即使其中非零元素非常少。

补充内容 (2013-11-13 14:17):
关于循环以及判断的补充：
在matlab里面，for循环的速度是远远快于while的(c/c++中相差不大)。因此，尽量用for而不是while。如果空间不足，就用while吧。

补充内容 (2013-11-13 14:20):
if的速度与switch相比，到底哪个更快，暂时还没有定论。我是倾向于if的，switch更像是一个函数，而不是命令过程。如果有人有独到看法可以提出来。

补充内容 (2013-11-13 14:25):
关于内存空间：
如果存在大量临时变量（用的次数不多），会使得连续内存空间变少，分配大规模矩阵(或reshape运算)不能保证成功。因此应该避免。