提高matlab代码运行效率
Matlab是一种解释性语言，追求的是方便性、灵活性以及交互性，因此在快速性上要比C语言这种性能强劲著称的稍逊一筹。然而，通过一些手段，我们也能让MATLAB语言快起来，甚至和C差不多了!
1可行方法
1.1循环矢量化
MATLAB变量的基本类型是矩阵，当对矩阵的每个元素循环处理时，运算速度很慢。利用MATLAB提供的用于矢量化操作的函数，把循环矢量化，这样既可以提高编程效率，也可以提高程序的执行效率。下面给出一个循环的例子：
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函数来查看上述各代码运行的时间，可以发现，把数组看作一个整体，进行操作后，执行效率提高约300倍。
另外，在必须使用多重循环的情况下，建议在循环的外环执行循环次数少的，内环执行循环次数多的，这样也可以显著提高程序执行速度。
下面再举一个例子
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的十倍左右。
1.2在使用数组或矩阵之前先定义维数
MATLAB中的变量在使用之前不需要明确地定义和指定维数。但当未预定义数组或矩阵的维数时，当需赋值的元素下标超出现有的维数时，MATLAB 就为该数组或矩阵扩维一次，这样就会大大降低程序的执行效率。因此，在使用数组或矩阵之前，预定义维数可以提高程序的执行效率。
1.3 对矩阵元素使用下标或者索引操作
在MATLAB中，矩阵元素的引用可用两个下标来表示。例如：A(i,j) 表示矩阵的第i行第j列的元素；A(1:k,j)表示矩阵A的第j列的前k个元素；A(:,j) 表示矩阵的第j列的所有元素。求矩阵A的第j列元素的平均值的表达式为mean(A(:,j))。
1.4用函数代替脚本文件
因为每次调用MATLAB的脚本文件都需要将不必要的中间变量加载到内存中，每执行一次，就加载一次。函数在调用时被编译成了伪代码，只需要加载到内存一次。当多次调用同一个函数时会运行快一些。因此尽量多使用函数文件而少使用脚本文件，也是提高执行效率的一种方法。
1.5用Mex文件编写循环代码
Matlab提供了与C和C++的接口，那么我们可以在用C或C++语言编写耗时的循环代码，然后通过接口程序在Matlab中转换成dll文件，这就是我们所要的Mex文件，可以在MATLAB环境下直接执行。通过这种方法可以极大地提高计算速率。一般来说，C-MEX 文件的执行速度是相同功能的M文件执行速率的20~40倍。
1.6尽可能利用matlab内部提供的函数
因为matlab内部提供的函数绝对是各种问题的最优算法，那写程序都是他们大师级人物写出来的，程序应该说相当高效，有现成的为什么不用那！这个原则就不用实际的程序测试了。
1.7给数组或矩阵预分配内存
特别是使用大型数组或矩阵时，Matlab进行动态内存分配和取消时，可能会产生内存碎片，这将导致大量闲置内存产生，预分配可通过提前给大型数据结构预约足够空间来避免这个问题。
这一点原则是极其重要的，以至于在编写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
1.8内存管理函数和命令
1)Clear variablename：从内存中删除名称为variablename的变量。
2)Clear all：从内存中删除所有的变量。
3)Save：将指令的变量存入磁盘。
4)Load：将save命令存入的变量载入内存。
5)Quit：退出MATLAB，并释放所有分配的内存。
6)Pack：把内存中的变量存入磁盘，再用内存中的连续空间载回这些变量。考虑到执行效率问题，不能在循环中使用。
1.9节约内存的方法
1)避免生成大的中间变量，并删除不再需要的临时变量。
2)当使用大的矩阵变量时，预先指定维数并分配好内存，避免每次临时扩充维数。
3)当程序需要生成大量变量数据时，可以考虑定期将变量写到磁盘，然后清除这些变量。当需要这些变量时，再重新从磁盘加载。
4)当矩阵中数据极少时，将全矩阵转换为稀疏矩阵。
1.10其他经验
1)多重循环让内层循环次数多
在必须使用多重循环时下，如果两个循环执行的次数不同，则在循环的外环执行循环次数少的，内环执行循环次数多的，这样可以提高速度。
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循环的嵌套顺序对于速度是有影响的，虽然相对来说差别不是很大。
2.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 =; for k=1:10000, sin(A(k)),
end;
6)当某条操作改变了原来变量的数据类型或形状(大小，维数)时将会减慢运行速度。
7)应该这样使用复常量x = 7 + 2i，而不应该这样使用：x = 7 +
2*i,后者会降低运行速度。”
8)“尽量用向量化的运算来代替循环操作。
最常用的使用vectorizing技术的函数有：All、diff、ipermute、permute、reshape、ueeze、y、find、logical、prod、shiftdim、sub2ind、cumsum、ind2sub、ndgrid、repmat、sort、sum 等。”优先使用matlab内建函数，将耗时的循环编写进MEX-File中以获得加速。
3.MATLAB代码矢量化指南
这里列出了《MATLAB代码矢量化指南(译)》
3.1基本技术
1)MATLAB索引或引用
在MATLAB中有三种基本方法可以选取一个矩阵的子阵。它们分别是下标法，线性法和逻辑法(subscripted,linear,
and logical)。
如果你已经熟悉这个内容，请跳过本节
1、下标法非常简单，看几个例子就好。A = 6:12; A([3,5])ans =8 10A([3:2:end])ans =8 10 12
A = [11 14 17; 12 15 18;
13 16 19]; A(2:3,2)ans =1516
2、线性法
二维矩阵以列优先顺序可以线性展开，可以通过现行展开后的元素序号来访问元素。
A = [11 14 17;
12 15 18;
13 16 19]; A(6)ans =16
A([3,1,8])ans =13 11 18A([3;1;8])ans =131118
3、逻辑法
用一个和原矩阵具有相同尺寸的0-1矩阵，可以索引元素。在某个位置上为1表示选取元素，否则不选。得到的结果是一个向量。
A = 6:10; A(logical([0 0 1 0 1]))ans =8 10A = [1 2 3 4]; B = [1 0 0 1]; A(logical(B))ans =1
4
2)数组操作和矩阵操作
对矩阵的元素一个一个孤立进行的操作称作数组操作；而把矩阵视为一个整体进行的运算则成为矩阵操作。MATLAB运算符*,/,\,^都是矩阵运算，而相应的数组操作则是.*, ./, .\,
.^A=[1
0 ;0 1]; B=[0
1 ;1 0]; A*B
% 矩阵乘法ans
=0
11
0A.*B
% A和B对应项相乘ans
=0
00
0
3)布朗数组操作
对矩阵的比较运算是数组操作，也就是说，是对每个元素孤立进行的。因此其结果就不是一个“真”或者“假”，而是一堆“真假”。这个结果就是布朗数组。D
= [-0.2 1.0 1.5 3.0 -1.0 4.2 3.14]; D
>= 0ans
=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)从向量构建矩阵
在MATLAB中创建常数矩阵非常简单，大家经常使用的是：
A = ones(5,5)*10
但你是否知道，这个乘法是不必要的？A = 10; A = A(ones(5,5)) A =10 10 10 10
1010 10 10 10
1010 10 10 10
1010 10 10 10
1010 10 10 10 10
类似的例子还有：v = (1:5)'; n = 3; M = v(:,ones(n,1)) M =
1 1 12 2 23 3 34 4 45 5 5
事实上，上述过程还有一种更加容易理解的实现方法：
A = repmat(10,[5
5]); M = repmat([1:5]', [1,3]);
其中repmat的含义是把一个矩阵重复平铺，生成较大矩阵。更多详细情况，参见函数repmat和meshgrid。
5)相关函数列表
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 矩阵
3.2扩充的例子
6)作用于两个向量的矩阵函数
假设我们要计算两个变量的函数FF(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)排序、设置和计数
在迄今为止讨论过的例子中，对向量中一个元素的计算都是独立于同一向量的其他元素的。但是，在许多应用中，你要做的计算则可能与其它元素密切相关。例如，假设你用一个向量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)稀疏矩阵结构
在某些情况下，你可以使用稀疏矩阵来增加计算的效率。如果你构造一个大的中间矩阵，通常矢量化更加容易。在某些情况下，你可以充分利用稀疏矩阵结构来矢量化代码，而对于这个中间矩阵不需要大的存储空间。
假设在上一个例子中，你事先知道集合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)附加的例子
下面的例子使用一些在本技术手册中讨论过的方法，以及其它一些相关方法。请尝试使用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处理稀疏矩阵的方式。
4.建议
对于撰写高效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的接口调用，往往会带来非常巨大的开销(成千上万倍的减慢)，甚至使得本来几乎不是问题的代码成为性能瓶颈。
