【转载】让你的MATLAB运行效率更快一些吧！

转自 https://www.digquant.com.cn/forum.php?mod=viewthread&tid=258

一、速度

1.1 For-循环

1、改变算法，多用矩阵运算（尤其是矩阵乘法），尽量减少for循环；

2、减少for循环中的函数调用；

传统观点认为for-loop是影响性能的致命环节，让我们来对此验证：

tic
toc

Elapsed time is 0.000239 seconds.

tic
toc
for i=1:1000000
end

Elapsed time is 0.000050 seconds.

从上面的实验结果可以得出以下结论：

1、tic/toc语句的时间开销可以忽略不计
2、for-loop语句本身的时间开销也非常小，关键的影响效率的地方不在于循环本身，而是在于循环的内部。
3、tic/toc不一定要成对出现，一个tic后面可以有多个toc，但需要需要重新计时的时候，要再次执行tic。
4、toc的结果可以用变量接收下来，如：

T=toc
T(k)=toc;

接下来我们就借助for循环，分析一下其他的各个影响效率的因素。

内建函数

tic
for i=1:1000000
cos(0);
end
toc

Mean elapsed time is 0.032866 seconds.

m-函数

tic
for i=1:1000000
func(i);
end
toc

Mean elapsed time is 0.185556 seconds.

function func( ~ )
end

匿名函数

tic
for i=1:1000000
funca(i);
end
toc

Mean elapsed time is 0.561228 seconds.

funca=@(x)'';

内联函数

tic
for i=1:1000000
funci(i);
end
toc

Mean elapsed time is 19.5606 seconds.

funci=inline('','x');

image

从上面的实验结果可以得出以下结论：
1、内联函数的调用时间开销最小，约为for-loop本身的10倍
2、m-函数的调用时间开销约为内联函数的6倍，约为for-loop本身的60倍
3、匿名函数的调用时间开销约为m-函数的3倍，约为for-loop本身的187倍
4、内联函数的调用时间开销过大，尽量不要在循环中使用
5、另外MEX-函数的调用时间开销，理应介于内联函数和m-函数之间

矩阵索引

tic
A=zeros(1000000,1);
for i=1:1000000
A(i)=i;
end
toc

Mean elapsed time is 0.007592 seconds.

tic
A=zeros(1000000,1);
for i=1:1000000
A(i,1)=i;
end
toc

Mean elapsed time is 0.007954 seconds.

tic
A=zeros(1000000,1);
for i=1:1000000
A(i:i,1)=i;
end
toc

Mean elapsed time is 0.663598 seconds.

tic
A=zeros(1000000,1);
for i=1:1000000
A(i,:)=i;
end
toc

Mean elapsed time is 0.273345 seconds.

tic
A=zeros(1000000,1);
for i=1:1000000
A(i,1:1)=i;
end
toc

Mean elapsed time is 0.730042 seconds.

tic
A=zeros(1000000,1);
for i=1:1000000
A(i:i,1:1)=i;
end
toc

Mean elapsed time is 1.00852 seconds.

image

二、内存

2.1 内存分配

tic
A=zeros(1000000,1);
for i=1:1000000
A(i)=i;
end
toc

Mean elapsed time is 0.009025 seconds.

tic
% A=zeros(1000000,1);
for i=1:1000000
A(i)=i;
end
toc

Mean elapsed time > 20 minutes.

因此，如果不预先分配好内存，将会大大增加仿真时间，拖慢执行效率。

所幸的是，由于这个现象的重要性，Matlab的编辑器能够发现并提示这个问题，会用红的波浪线~标记出来。

三、向量化

N=0:0.1:1000;
for i=0:10000
y(i)=cos(N(i))；
end

向量化：

N=0:0.1:1000;
y=cos(N);

MATLAB向量化函数
accumarray函数
arrayfun函数
bsxfun函数
cellfun函数
spfun函数

3.1 accumarray函数

A = accumarray(subs,val)
A = accumarray(subs,val,sz)
A = accumarray(subs,val,sz,fun)
A = accumarray(subs,val,sz,fun,fillval)
A = ccumarray(subs,val,sz,fun,fillval,issparse)

val = 101:105;
subs = [1; 2; 4; 2; 4]
A = accumarray(subs, val)

A =

101
206
0
208

subs =
1 1 1
2 1 2
2 3 2
2 1 2
2 3 2

val =
101
102
103
104
105

1、val的元素个数与subs的行数是一致的。

2、A = accumarray(subs, val)的实现过程分成2步。

第一步
是把val中的元素，按照subs对应行所给出的下标放到一个新的cell矩阵B中（cell是为了方便解释，也就是说B矩阵中的每个位置可以放入多个数值），注意，subs的值是B的下标，不是val的。举例来说，subs第一行[ 1 1 1],意思就是把val中第一个元素（val（1））放入到B（1，1，1）的位置，依次类推，val（2）放入到B(2 1 2)，val（3）放入到B(2 3 2)，val（4）放入到B(2 1 2)，val（5）放入到B(2 3 2)。此时，可以看到B（1，1，1）中有1个数（val（1））；B(2 1 2)有2个数（val（2），val（4））；B(2 3 2)也有2个数（val（3），val（5））。

第二步
把B中每个单元中的数分别累加，并放入到A的对应位置。

注：accumarray默认的是把每个单元中的数累加，因为对每个单元中的数的默认处理函数是sum。可以通过A = accumarray(subs,val,[],[@fun](https://github.com/fun "@fun"))的调用格式来指定其他的处理函数，比如说mean。对指定的fun函数的要求是，接受列向量输入，输出单个的数值型，字符型或逻辑型变量。A的维数与B相同，A中的元素默认为零。A的大小为max(subs(1))×max(subs(2))×max(subs(3))…

很显然，A的维数与subs的列数相等。

A = accumarray(subs, val)
A = accumarray(subs,val,sz)
sz 可以用来指定A大小，但是不能小于A = accumarray(subs, val)得到的A的大小。比如A = accumarray(subs, val)的到A是一个3×4的二维矩阵，那么sz应当为一个包含2个元素的向量sz=[m1,m2] （sz向量的长度和A的维数相等），其中，m1大于等于3，m2大于等于4. 但是，当得到的A是一个p×1的一维向量时，sz=[m，1]，m大于等于p。另外，sz可以赋值为空，表示由函数自动决定A的大小。
A = accumarray(subs,val,sz,fun)
fun可以指定专门的处理函数，默认的处理函数为sum
A = accumarray(subs,val,sz,fun,fillval)
fillval指定A中元素的默认值。可以等于NaN
A = ccumarray(subs,val,sz,fun,fillval,issparse)
isspares选择A是否使用稀疏矩阵的格式
A = accumarray({subs1, subs2, ...}, val,...)
{subs1, subs2, …}，等同于A = accumarray(subs, val,…)，此时，subs=[subs1, subs2, …]或者=[subs1；subs2； …]

例子：
1000人，身高分布在170_{180cm，体重在110}100斤，年龄分布在20~50岁，计算身高体重都相等的人的年龄平均值。结果用矩阵来表示：行数表示身高，列数表示体重，矩阵元素表示年龄的平均值。

height=unidrnd(10,1000,1)+170; %身高的数据
weight=unidrnd(90,1000,1)+110; %体重的数据
old=unidrnd(30,1000,1)+20;
mo=accumarray([height,weight],old,[],@mean);
%或
mo=accumarray([height,weight],old,[],@mean,0,true);

3.2 arrayfun函数

arrayfun函数实现的是将指定的函数应用到给定数组在内的所有元素。这样以前不可避免的循环现在可以向量化了。

生成一个这样的n×n矩阵

a:a(i,j)=dblquad(@(u,v)sin(u)*sqrt(v),0,i,0,j)，以n=10为例。
a=zeros(10);
for ii=1:10
for jj=1:10
a(ii,jj)=dblquad(@(u,v) sin(u)+sqrt(v),0,ii,0,jj);
end
End
%现在只需要如下调用
[J,I]=meshgrid(1:10);
a1=arrayfun(@(ii,jj) dblquad(@(u,v) sin(u)+sqrt(v),0,ii,0,jj),I,J);

3.3 bsxfun函数

以前，当我们想对一个矩阵A的每一列或每一行与同一个向量a进行某些操作（比较大小、乘除等）时，只能用循环方法或者利用repmat函数将要操作的向量a复制成和A一样尺寸的矩阵，进而进行操作。从Matlab R2007a开始，有了更有效的方法，那就是bsxfun函数。
有如下矩阵：

image

向量为b=[1 2 3]T，请找出b在A矩阵列中的位置loc=[1,4]。

方法1：

A=[1 2 3 1;2 3 4 2;3 3 8 3]
b=[1;2;3]
loc = find(all(bsxfun(@eq,A,b)))
%把A的每一列和b用==（@eq）来判断，找出全1的列；

方法2：

loc = find(arrayfun(@(n)all(A(:,n)==b),1:4))
%用arrayfun对n进行1:4的遍历；

方法3：

loc =find(all(~bsxfun(@minus,A,b)))
%把A的每一列和b来相减( @minus )，求反后找出全1的列；

方法4：

loc=find( arrayfun(@(n) isequal(A(:,n),b),1:4))
%用arrayfun对n进行1:4的遍历后用isequal函数来判断A的每列和b是否相等；

方法5：

loc=find(b'*A==sum(b.^2))
%的转置和A相乘，然后和b.^2每列的和进行比较，找到相等的;

3.4 cellfun函数

A={'Hello', 'MATLAB', 'I love MATLAB', 'MATLAB is powerful', 'MATLAB is the language of technical computer'}

A={‘Hello’, ‘MATLAB’, ‘I love MATLAB’, ‘MATLAB is powerful’, ‘MATLAB is the language of technical computer’};

cellfun(@length,A)

ans =
5 6 13 18 44

3.5 spfun函数

a=sparse([1 3 20 60 100],[2 20 30 60 80],1:5)

</pre>

a =

(1,2) 1
(3,20) 2
(20,30) 3
(60,60) 4
(100,80) 5

sa=spfun(@(x) x.^2+1,a)

</pre>

sa =
(1,2) 2
(3,20) 5
(20,30) 10
(60,60) 17
(100,80) 26

四、函数化

尽量使用内建函数，内建函数的速度是最快的。
m-函数的执行效率也很高。
MEX-函数的执行效率仅次于内建函数，将耗时的代码写成MEX-函数,将大大提高运行速度。
匿名函数，内联函数，以及一些面向对象方法，尽量不要在执行次数多的循环体内使用。

五、预分配内存

A= zeros(1000, 1);
A = int8(zeros(100, 1));
A = zeros(1000, 1, 'int8');

常用的预分配内存函数：

zeros
ones
eye