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() 
ans = 
8 10 
A() 
ans = 
8 10 12 
A = 
; 
A(2:3,2) 
ans = 
15 
16 
1.2)线性法 
二维矩阵以列优先顺序可以线性展开，可以通过现行展开后的元素序号 
来访问元素。 
A = 
; 
A(6) 
ans = 
16 
A() 
ans = 
13 11 18 
A() 
ans = 
13 
11 
18 
1.3)逻辑法 
用一个和原矩阵具有相同尺寸的0-1矩阵，可以索引元素。在某个位置上为1表示选取元素，否则不选。得到的结果是一个向量。 
A = 6:10; 
A(logical()) 
ans = 
8 10 
A = 
; 
B = ; 
A(logical(B)) 
ans = 
1 4 
----------------------------------------------------- 
2)数组操作和矩阵操作(Array Operations vs. Matrix Operations) 
对矩阵的元素一个一个孤立进行的操作称作数组操作；而把矩阵视为一个整体进行的运算则成为矩阵操作。MATLAB运算符*,/,\,^都是矩阵运算，而相应的数组操作则是.*, ./, .\, .^ 
A=; 
B=; 
A*B % 矩阵乘法 
ans = 
0 1 
1 0 
A.*B % A和B对应项相乘 
ans = 
0 0 
0 0 
------------------------------------------------------ 
3)布朗数组操作(Boolean Array Operations) 
对矩阵的比较运算是数组操作，也就是说，是对每个元素孤立进行的。因此其结果就不是一个“真”或者“假”，而是一堆“真假”。这个结果就是布朗数组。 
D = ; 
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,); 
M = repmat(', ); 
其中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)'; 
= 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 差分运算符 
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(); 
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(); 
count = diff(find()); 
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(,,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(); 
V_A2=cumsum(); 
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
向量累加，每5个元素进行一次： 
工具：CUMSUM , 向量索引 
优化前： 
% Use an arbitrary vector, x 
x = 1:10000; 
y = ; 
for n = 5:5:length(x) 
    y = ; 
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 = ; 
为： 
x = ; 
解(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 = ; 
% The indices associated with the values are summed cumulatively 
indices = ; 
totals = zeros(max(indices),1); 
for n = 1:length(indices) 
    totals(indices(n)) = totals(indices(n)) + values(n); 
end 
优化后： 
indices = ; 
totals = full(sparse(indices,1,values)); 
注意：这一方法开辟了稀疏矩阵的新用途。在使用sparse命令创建稀疏矩阵时，它是对分配到同一个索引的所有值求和，而不是替代已有的数值。这称为"向量累加"，是MATLAB处理稀疏矩阵的方式。
凡是有头像的地方，表示该向量或者矩阵的所有行，在Matlab中的表示，然后呢，在这里发布的时候即变成了头像，呵呵，蛮有意思的！x( : )即表示x的所有行或者列。