[阿Q精品原创]Matlab代码效率篇之原理分析及加速方案--持续更新

最新推荐文章于 2025-05-13 17:15:55 发布

SIPANGYIYOU

最新推荐文章于 2025-05-13 17:15:55 发布

阅读量859

点赞数 3

分类专栏： matlab视觉 matlab加速 matlab 文章标签： matlab 开发语言矩阵

本文链接：https://blog.youkuaiyun.com/SIPANGYIYOU/article/details/121396510

版权

matlab视觉同时被 3 个专栏收录

2 篇文章

订阅专栏

matlab加速

1 篇文章

订阅专栏

matlab

1 篇文章

订阅专栏

Matlab代码效率原理分析及加速方案

序言

大家是否还在挣扎于matlab漫长运行的等待中? 阿Q在这里整理加速matlab代码各种方式, 在这里也分享了阿Q加速matlab代码的经验和技巧,希望对大家有所帮助.Blog后续会时不时更新,撰写不易,如果大家觉得不错的话,希望能点个收藏和关注,阿Q在此谢过啦.

阿Q在这里整理了加速方式的使用表,推荐列和难度列里没有括号的表示简单代码下的星数, 用()表示复杂代码下的星数.

加速方式	适用情况(和局限性)	推荐指数	难度
1.提前分配变量	代码中可能会发生重复且动态分配变量内存的时候,例如for循环中(无局限)	****	*
2.向量化或者矩阵化	使用for循环的时候可以考虑使用(复杂代码难向量化)	*****	* (***)
3.逻辑变量	适用于需要进行条件判定或者选择某些向量/矩阵元素的时候(复杂代码难操作)	*****	同上
4. parfor 的使用	使用parfor(或者其他并行方式)可以方便通过多线程并行来加速代码, 很适合在蒙泰卡诺的重复测试中使用(加速倍率受限严重, 最大不超过多线程的数量)	***	**
5.cell的动态拼接	适合需要向量和矩阵动态拼接的时候(无明显局限)	***	**
6.使用稀疏矩阵	适合算法中存在大量稀疏矩阵,并且在计算上时间消耗过大的情况(稀疏矩阵计算速度虽然很快,但是初始化赋值的操作在Matlab里并不可观. 在初始化的时间消耗上, 编程能力高低会有质的区别)	(***)	* (***)
7.cellfun等加速	适合把for循环中的某些复杂代码部分设计成支持向量化编程函数(加速效果不稳定,加速多少也跟编程能力相关)	***	**
8.1Mex加速	适用于Matlab下编程各类加速都用了之后, 实在提速不上来的时候使用. 将需要加速的部分使用Mex来混合编写. 这种加速效果基本取决于Mex部分开发的算法性能.相比于8.2system方式, 这种方式不需要写文件读写操作, 变量可以直接传递(需要外部语言开发, 加速性能取决于外部开发的程度. 此外, 虽然不用写文件读写接口,但Matlab与外部语言开发的mex库文件之间的数据调用所传输的时间可能远比文件读写要久)	***	****
8.2 system调用	与8.1类似,适用于Matlab下编程各类加速都用了之后, 实在提速不上来的时候使用. 将需要加速的部分直接用外部语言开发. 这种加速效果取决于外部开发的算法性能. 相比于Mex编程, 可以方便后续Matlab上测试好的代码逐渐移植到外部开发上进行系统地提速. 因此, 更推荐一点. (与Matlab交互采用文件接口, 因此需要写文件读写接口)	*****	****
编程习惯	提前分配变量内存 + for循环嵌套时候内外层的设计等	*****	**

[ Hi~ o(￣▽￣)ブ, let’s start]

用过matlab编程的人应该能体会在matlab平台上程序运行效率与代码本身编写水平有很大的关系. 在这篇文章里, 阿Q首先以如下计算为例, 进行提速分析. 根据变量 $x$ 计算并存储变量 $y$ 和 $z$ :
(1) $x$ 取值在0~10000内的整数, 求 $y = c o s (x)$ 的结果;
(2)此外,如果 $y < 0.5$ 的部分, $z = c o s (y)$ ; 反之 $z = s i n (y)$ .

运行效率最慢的代码如下:

%% code.0
clear
time = 0;
num_MC_test = 1000;
min_x = 0;
max_x = 10000;
size_x = max_x-min_x+1;
for i = 1:num_MC_test
    clearvars y z 
    tic;
    for j = 1:size_x
        y(j) = cos(j);
        if y(j) < 0.5
            z(j) = cos(y(j));
        else
            z(j) = sin(y(j));
        end
    end
    time = time + toc;
end
    
disp(['code0 time = ',num2str(time),'s.'])
%code0 time = 2.1776s.

这是一个简单的示例. 大家在实际编程过程中, 编写的算法肯定远比这复杂. 如果一直采用如上编程的方式, 随着代码复杂程度的增加, 运行的时间就会大幅度提高. 那么, 我们应该如何去提升matlab代码效率呢?

第一种提速方式:提前分配变量内存

提速原理分析:在code.0中,可以看到待求变量 $y$ 和 $z$ 并未提前分配内存, 这会降低Matlab程序运行效率. 这是因为,Matlab是一门解释性语言,也就是运行的时候才会告诉机器如何去执行这一句代码.在code.0中执行到y(j) = cos(j) 以及 z(j) = ...的时候, Matlab将会暗自拆解为如下三步进行操作:
(step.1) 生成一个新的长度为j的向量;
(step.2) 将原向量(j-1长度)中的值拷贝给新向量的第0~j-1元素处;
(step.3) 最后执行前面的赋值代码y(j) = cos(j)和z(j) = ....
加速方案: 了解原理之后, 就可以大刀阔斧的改进算法啦~(≧▽≦)/. 以code.0为例,加速代码如下:

%% code.1
clear
num_MC_test = 1000;
min_x = 0;
max_x = 10000;
size_x = max_x-min_x+1;
time = 0;

for i = 1:num_MC_test
    clearvars y z
    tic;
    % code.1: assign memories of the estimated vars in advance: y and z
    y = zeros(1,size_x,'double');
    z = zeros(1,size_x,'double');

    for j = 1:size_x
        y(j) = cos(j);
        if y(j) < 0.5
            z(j) = cos(y(j));
        else
            z(j) = sin(y(j));
        end
    end
    time = time + toc;
end
    
disp(['code1 time = ',num2str(time),'s.'])
%code1 time = 0.34514s. [code0 :2.1776s.]

嘿嘿,果然提速了~效果还不错, 不过我们还能不能更进一步的提速呢?想必大家都跟阿Q一样, 不会满足于此. 在提到新的加速方式前,先纠正一个容易被误解的部分.

变量预分配效率的额外测试, 误区纠正 :
-blog参考:预分配内存提高运行效率以及时间比较

clear
tic;
a(10000,20000) = 0;      %方法一，直接赋值为零
time=toc;
disp(['用直接赋值方式,time=',num2str(time),'秒.'])

clear
tic;
b=zeros(10000,20000);   %方法二，用zeros函数，不指定数据类型
time=toc
disp(['用zeros函数，不指定数据类型,time=',num2str(time),'秒.'])

clear
tic;
c=zeros(10000,20000,'double'); %方法三，用zeros函数，指定数据类型
time=toc
disp(['用zeros函数，指定数据类型,time=',num2str(time),'秒.'])

-缘由: 大家可能会好奇, 为什么阿Q要进行这个额外测试. 起因是如上blog参考文章, 为此阿Q专门进行了测试,可测试结果与博主的结果完全不一致(原文博主效率排序是直接最优,zeros指定类型次之,zeros不指定最低):

用直接赋值方式,time=0.47397秒.
用zeros函数，不指定数据类型,time=0.000129秒.
用zeros函数，指定数据类型,time=3.8e-05秒.

-阿Q的建议: 不同电脑,不同版本的Matlab运行结果可能都有不同. 建议大家最好在自己电脑的Matlab测试相关代码. 阿Q这里的测试与Matlab系列教材书里推荐的一致, 时间效率排序: zeros指定类型效率最高, zeros不定类型次之, 直接法最低.

第二种提速方式:向量化编程取代for循环

为了进一步提速,是时候祭出Matlab强烈推荐的大杀器了:向量化编程. (不管阻挡代码运行效率的是何方妖魔鬼怪, 在向量化面前一切,都是渣渣,吼吼吼)

提速原理分析: Matlab号称矩阵实验室,本身就是为矩阵和向量运算服务的. Matlab本身是c开发的一门高级语言, 在运算矩阵和向量上人家效率绝对是杠杠的, 这相当于为我们高效率编程开了一个特别通道. 如果我们把代码的运算部分, 尽可能的向量化或者矩阵化运算, 那么代码运行效率将会大幅度提升.
加速方案: 那么以code.0为例,加速代码如下:

%% code.2
clear
num_MC_test = 1000;
min_x = 0;
max_x = 10000;
size_x = max_x-min_x+1;
time = 0;

for i = 1:num_MC_test
    clearvars x y z
    tic;
    % code.1: assign memories of the estimated vars in advance: z
    z = zeros(1,size_x,'double');
    
    % code.2: step.1 vectorize the variable x
    x = min_x:max_x;
    
    % code.2: step.2 For variable y, replace the `for cycle` by vectorized operation
    y = cos(x);
    for j = 1:length(y)
        if y(j) < 0.5
            z(j) = cos(y(j));
        else
            z(j) = sin(y(j));
        end
    end
    time = time + toc;
end
    
disp(['code2 time = ',num2str(time),'s.'])
%code2 time = 0.32192s.[code1: 0.34514s.  code0 :2.1776s.]

em…提速好像不明显诶. 让我们来回顾下, 我们先是通过step.1提前向量化变量x, 然后step.2直接向量化计算变量y的值 (注意:这样一来,y变量计算不用for循环了,所以我们就不用提前在code.1处给y变量预分配内存了). 阿Q表示很郁闷,说好的Matlab大杀器,为什么提速这么少呢?

其实, 这提速效果是由规模和问题复杂程度来决定的, code0所处理的问题,在for循环中需要处理的计算很简单, 而且在for规模不大的情况下,向量化提速的效果其实并不会很明显. 因此,接下来的一个规模扩大的额外测试就有必要给大家展示下了(复杂计算就不测试了):

-扩大规模测试:我们将x的最大取值从1e4增加到1e5和1e6来看下测试结果

[max_x = 1e4;]
code0 time = 2.1776s.
code1 time = 0.34514s.
code2 time = 0.32192s.

[max_x = 1e5;]
code0 time = 23.6586s.
code1 time = 3.8735s.
code2 time = 2.1596s.

[max_x = 1e6;]
code0 time = 244.0673s.
code1 time = 42.0152s.
code2 time = 20.9031s.

可以看出,随着规模扩大,向量化编程提速的效果会越来越明显:从1e4提速一丢丢,到1e6提速超一倍.那么试想,对于1e4 x 1e4的矩阵, 如果类似向量化或矩阵化操作, 程序效率将会大大提升.

看到这里, 不知道老铁们会不会觉得变量z部分的处理,看起来太不美观了.而且,我们是否也能将z变量的部分进行向量化提速呢? 嘿嘿,答案是可以的~阿Q后续在第三种提速方式中进行改进.

矩阵化提速: 矩阵化提速与向量化提速的方式类似,就是把原本的问题预先组合成矩阵(类似code2 step.1).此外,矩阵化的操作方式比向量化会复杂很多.这里就不详细阐述了.
向量化技巧示例: Matlab中有很多内置函数是支持向量化和矩阵化处理的,比如norm, cross等.这里就不一一列举了, 主要说几个常用的:
(1) 给定一个由 $n$ 个3维向量构成的矩阵 $X=(x_1,...,x_n)$ ,如何计算每一个3维向量的 $L_2$ 范数,并构成一个新的 $n$ 维度向量 $y=(||x_1||_2,...,||x_2||_2)$ . 可以从如下示例代码看到向量化编程下,效率提速了三个数量级.

%%example.1
n = 100;
X = rand(3,n);
tic;
y = sqrt(sum(X.^2,1));
time = toc;

disp(['example.1 speedy_time = ',num2str(time),'s.'])
%----------------------------------------------------------------------
tic;
y = zeros(1,n,'double');
for i = 1:n
    y(i) = norm(X(:,i));
end
time = toc;

disp(['example.1 time = ',num2str(time),'s.'])
%----------------------------------------------------------------------
%example.1 speedy_time = 0.000171s.
%example.1 time = 0.10153s.

(2) 给定两个 $3n \times 3$ 矩阵 $A$ 和 $B$ , 其中 $A=(A_1^T,...,A_n^T)^T,B=(B_1^T,...,B_n^T)^T$ .
求矩阵 $C=((A_1B_1)^T,...,(A_nB_n)^T)^T$ .

%%example.2
n = 10000;
A = rand(3*n,3);
B = rand(3*n,3);
%----------------------------------------------------------------------
%% speedy version
tic;
C = zeros(3*n,3,'double');
%compute the 1st column of C
tmp_B = B(:,1);
tmp_B = reshape(tmp_B,3,[])';
tmp_B = kron(tmp_B,ones(3,1));
C(:,1) = sum(tmp_B.*A,2);

%compute the 2nd column of C
tmp_B = B(:,2);
tmp_B = reshape(tmp_B,3,[])';
tmp_B = kron(tmp_B,ones(3,1));
C(:,2) = sum(tmp_B.*A,2);

%compute the 3rd column of C
tmp_B = B(:,3);
tmp_B = reshape(tmp_B,3,[])';
tmp_B = kron(tmp_B,ones(3,1));
C(:,3) = sum(tmp_B.*A,2);
time = toc;

disp(['example.2 speedy_time = ',num2str(time),'s.'])
%----------------------------------------------------------------------
%% for cycle version
tic;
C = zeros(3*n,3,'double');
for i = 1:n
    C(3*i-2:3*i,:) = A(3*i-2:3*i,:)*B(3*i-2:3*i,:);
end
time = toc;

disp(['example.2 for_time = ',num2str(time),'s.'])
%----------------------------------------------------------------------
%example.2 speedy_time = 0.002039s.
%example.2 for_time = 0.024967s.

第三种提速方式:逻辑向量或矩阵来取代if,find

这个方式的好处是可以将原本if和find等条件判别的代码变成向量化或者矩阵化代码.这是一个非常有用的技巧,在类似的加速blog中几乎没有见过有人介绍这个方式.阿Q在以前编写或者帮人家修改的代码中常常用到,用的好的话,配合向量化编程(第二种提速方式),可以让原本的程序提速数十倍到上百倍, 而且程序代码上看起来也会很美观.

原理分析: 逻辑向量(或矩阵)就是直接对当前需要判断的向量进行逻辑运算, 例如示例code中求解 $z$ 变量前需要判断 $y < 0.5$ 的条件, 我们就可以用简单的逻辑变量logical_var = y<0.5;来替换. 那么logical_var的变量是针对y向量, 凡是小于0.5的将会标记为1, 反之标记为0; 这样一来,我们就可以快速将符合条件的元素按照logical_var来向量化或者矩阵化提取.
加速方案:那么话不多说, 直接看加速示例就会对逻辑变量的操作和用途一目了然

%% code.3
clear
num_MC_test = 1000;
min_x = 0;
max_x = 1e4;
size_x = max_x-min_x+1;
time = 0;

for i = 1:num_MC_test
    clearvars x y z
    tic;
    % code.2: step.1 vectorize the variable x
    x = min_x:max_x;
    
    % code.2: step.2 For variable y, replace the `for cycle` by vectorized operation
    y = cos(x);
    
    % code.3: use logical variable to replace the `if` 
    z = cos(y);
    logical_var = y > 0.5;
    z(logical_var) = sin(y(logical_var));

    time = time + toc;
end

disp(['code3 time = ',num2str(time),'s.'])
%code3 time = 0.22806s.[code2: 0.32192s. code1: 0.34514s. code0: 2.1776s. ]
%注意,这里时间跟最一开始code0,code1代码中时间不完全一致了,因为在扩大规模测试那的时候,阿Q重新统计了下时间.

可以看到,当阿Q使用了逻辑变量后,在也不用if判断,而且也淘汰了for循环的使用,代码上明显简化和美观了不少. 同时,在时间上又有了明显的提速. 所以说, 如果说向量化是关公,那么逻辑变量就是青龙偃月刀,提速问题就如颜良文丑, 分分钟被关公耍大刀K.O. 主要的提速技巧,已经说明了,那么还有几种常用的加速方案,在这里也分享下.

第四种提速方式:parfor替换for循环

在示例code的测试中,阿Q用了num_MC_test来进行重复的Monte Carlo测试, 在用Matlab做实验测试算法性能的时候,这是大家经常用的手段. 可一旦有了这样的重复测试,就避免不了使用for循环. 而这样的循环有的时候没办法写成向量化编程的形式. 遇到这样的情况那就到了并行for循环表演时刻了.

更详细的并行编程,参考:matlab官网或https://blog.youkuaiyun.com/enjoyyl/article/details/41929033(个人觉得比较详细)
加速示例:

%% code.4
clear
num_MC_test = 1000;
min_x = 0;
max_x = 1e5;
size_x = max_x-min_x+1;
time = zeros(1,num_MC_test,'double');

tic;
for i = 1:num_MC_test
    % code.2: step.1 vectorize the variable x
    x = min_x:max_x;
    
    % code.2: step.2 For variable y, replace the `for cycle` by vectorized operation
    y = cos(x);
    
    % code.3: use logical variable to replace the `if` 
    z = cos(y);
    logical_var = y > 0.5;
    z(logical_var) = sin(y(logical_var));
end
time = toc;

disp(['code4 time = ',num2str(time),'s.'])

%% code.5
clear
num_MC_test = 1000;
min_x = 0;
max_x = 1e5;
size_x = max_x-min_x+1;

tic;
parfor i = 1:num_MC_test
    
    % code.2: step.1 vectorize the variable x
    x = min_x:max_x;
    
    % code.2: step.2 For variable y, replace the `for cycle` by vectorized operation
    y = cos(x);
    
    % code.3: use logical variable to replace the `if` 
    z = cos(y);
    logical_var = y > 0.5;
    z(logical_var) = sin(y(logical_var));

end
time = toc;

disp(['code5 time = ',num2str(time),'s.'])

%code4 time = 1.4536s.
%code5 time = 1.0836s.

可以看到,通过使用parfor是可以利用并行计算对MCtest进行提速.

parfor的骚操作 之save&load: 如果有用过parfor的老铁们,肯定有遇到过parfor的各种警告,比如slice之类的. for循环里的代码越复杂, parfor就越容易爆出各种错误, 让人非常头疼. 其实很多类似的问题,有一个很骚操作的解决办法: 就是将报错的代码部分,重新封装成一个函数.这里就以并行保存save和load来举例, parfor是不允许语句中出现save和load函数的, 可如果我们想在每次重复测试的时候,还涉及到生成,保存和读取当前测试下的数据,那这个骚操作就可以完美避免parfor无脑报错.

 %% code.6
clear
num_MC_test = 10;

parfor i = 1:num_MC_test
    % code 6: initialize the var x with variable length
    min_x = 0;
    max_x = unidrnd(100);
    size_x = max_x-min_x+1;

    % code.2: step.1 vectorize the variable x
    x = min_x:max_x;
    
    % code.2: step.2 For variable y, replace the `for cycle` by vectorized operation
    y = cos(x);
    
    % code.3: use logical variable to replace the `if` 
    z = cos(y);
    logical_var = y > 0.5;
    z(logical_var) = sin(y(logical_var));
    
    % code 6: save current variable x into a mat file
    mat_file = sprintf('tmp_%d.mat',i);
    AQ_save(mat_file,x);
end

% code 6: add a function to replace the 'save' function
function AQ_save(filename,x)
save(filename,'x')
end

如此一来,save和load函数的功能就可以在parfor中实现了. 如果大家遇到了parfor中某些奇怪问题又无从下手的时候,或许可以尝试这种方式来解决.

第五种提速方式:用cell代替向量/矩阵的动态拼接

在向量的拼接上, Matlab允许我们直接用[]来进行拼接向量和矩阵(如[a,b],其中a和b是向量).这是一种静态拼接的方式, 那么当我们在for循环中需要根据具体情况来进行动态拼接时, 我们也可以用类似的方式进行拼接,例如y = [y;[a,b]].这种方式虽然简单,实际上效率却并不高.

原因分析: 因为这种情况下,y向量需要预分配的内存大小我们预先是不知道的,只有在for循环处理完成后才清楚.因此,在每次for循环中都会遇到第1种情况遇到的问题,即y变量需要重复开辟内存.
解决办法: 解决的办法常用的有两种:
(1)通过cell变量来动态存储,最后进行拼接;
(2)灵活运用向量化编程(见前面的第二和三部分).
具体示例: 对于给定的 $x$ 向量(服从[-1,1]区间的均匀分布),我们需要提取其中小于0.5的部分,并动态拼接成一个字符串向量 $y=[y_1,...,y_n]$ , 其中 $y_i$ 正则表达为'x_%d = %.3f;\n'.

%% Method.5 test
size_x = 10000;
x = randn(1,size_x);
%----------------------------------------------------------------------
% approach.1: append variable using vec/array
tic;
y = '';
for i = 1:size_x
    if x(i) < 0.5
        tmp_str = sprintf('x_%d = %.3f;\n', i, x(i));
        y = [y,tmp_str];
    end
end
time = toc;
y1 = y;
disp(['approach.1 time = ',num2str(time),'s.'])
%----------------------------------------------------------------------
% approach.2: append variable using cell
clearvars y 
tic;
id_current = 1;
for i = 1:size_x
    if x(i) < 0.5
        y{id_current} = sprintf('x_%d = %.3f;\n', i, x(i));
        id_current = id_current + 1;
    end
end
time = toc;
y2 = cell2mat(y);
disp(['approach.2 time = ',num2str(time),'s.'])
%----------------------------------------------------------------------
% approach.3: append variable using vectorization
clearvars y 
tic;
logical_x = x < 0.5;
id = 1:size_x;
y = sprintf('x_%d = %.3f;\n', [id(logical_x);x(logical_x)]);
time = toc;

y3 = y;
disp(['approach.3 time = ',num2str(time),'s.'])
%----------------------------------------------------------------------
%approach.1 time = 0.14093s.
%approach.2 time = 0.061806s.
%approach.3 time = 0.005425s.

结论:
(1)可以看到, 如果能使用向量化编程加速的方式来加速,当然是最快的(见approach.3);
(2)如果实在无法向量化,则可以考虑用cell变量来进行动态拼接(见approach.2).

第六种提速方式:用稀疏矩阵替换稠密矩阵

使用Matlab处理大规模数据时候,经常遇到矩阵内部元素为0的情况.在这种情况下,可以考虑使用将原本的矩阵改为用稀疏矩阵来存储,然后使用稀疏计算来提速.在稀疏矩阵的使用过程中,需要注意两点:
(1)稀疏矩阵的赋值效率;
(2)稀疏矩阵的优点在于内存和计算效率的提升.
实际上,稀疏矩阵的赋值效率有时不尽人意,甚至不如稠密矩阵的初始化. 因此,在稀疏矩阵的构建上需要花些心思. 一旦构建了稀疏矩阵后,相关的计算效率就会大幅度提升.至于赋值效率这块,在这里就暂时不讨论了,后面阿Q会发布关于"稀疏矩阵的赋值方式和效率篇",后续发布后会链接于此处,敬请期待.

稀疏矩阵参考web推荐:官方web(yyds)
https://ww2.mathworks.cn/help/matlab/sparse-matrices.html?s_tid=CRUX_lftnav

第七种提速方式:用cellfun, boxfun, arrayfun来代替for循环

有的Blog很推荐cellfun,boxfun和arrayfun来向量化编程. 我的理解是,cellfun这种是针对for循环中复杂问题的代码来进行向量化编程的方式, 可以赋予用户自定义函数具有支持向量化编程的能力. 一般来说,使用这种方式来替换for循环,的确是可以提高代码性能的. 可用的不好,同样也会明显地降低代码性能.

问题: 在向量化编程的时候,最难受的就是必须要调用的函数如果不支持向量化的输入参数呢?
原理分析: Matlab自带的很多函数是支持输入参数的向量形式的, 而用户自定义的函数却没有这个功能. 一旦自定义函数内部过于复杂,我们也很难去进行向量化, 于是不得不用for循环来调用自定义函数.

为什么说是cellfun等可以处理复杂问题的向量化编程? 对比前面的加速方式三和四的就能理解了一点了: 在加速方式三和四中, for循环内部的处理只需要使用矩阵和向量的基本操作, 然后结合支持向量化处理的函数(比如cross())来进行向量化编程, 就可以替换掉for.
解决方案: 或许可以利用cellfun这类函数赋予了用户自定义函数具有向量化处理的功能.

clear;
m = 10000;
A = cell(m,1);
B = cell(m,1);
for i = 1:m
    n = randperm(10,1);
    tmp_data = [];
    tmp_data(:,1) = rand(n,1);
    tmp_data(:,2) = randperm(n,n);
    tmp_data(:,3) = i*ones(n,1);
    A{i,1} = tmp_data;
    
    tmp_data = [];
    tmp_data(:,1) = rand(n,1);
    tmp_data(:,2) = randperm(n,n);
    tmp_data(:,3) = i*ones(n,1);
    B{i,1} = tmp_data;
end

B_mat = cell2mat(B);

%% for cycle test
tic;

size_A = size(A,1);
test0 = cell(size_A,1);

for i = 1:size_A
    tmp_A = A{i};
    tmp_B = B{i};
    [tmp_min,tmp_min_id] = min(tmp_A(:,1));
    tmp_id = tmp_A(tmp_min_id,2);
    test0{i,1} = tmp_B(tmp_id,1);
    
    tmp_num = size(tmp_A,1);
    
    for j = 1:tmp_num
        tmp_id_A = tmp_A(:,2) == j;
        tmp_id_B = tmp_B(:,2) == j;
        
        tmp_s = tmp_A(tmp_id_A,1) + tmp_B(tmp_id_B,1);
        test0{i,1} = test0{i,1} + tmp_s;
    end
end
toc;


%%
tic;
test4 = cellfun(@(X,Y)test_operator2(X,Y),A,B);
toc;

%%
tic;
test3 = cellfun(@(X)test_operator(X,B_mat),A);
toc;

function output = test_operator(A,B_mat)
[~,id_min] = min(A(:,1));
id = ismember(B_mat(:,[2,3]), A(id_min,[2,3]),'rows');
output = B_mat(id,1);

end

function output = test_operator2(A,B)
[~,id_min] = min(A(:,1));
id = B(:,2) == A(id_min,2);
output = B(id,1);

tmp_num = size(A,1);
    
for j = 1:tmp_num
    tmp_id_A = A(:,2) == j;
    tmp_id_B = B(:,2) == j;
    
    tmp_s = A(tmp_id_A,1) + B(tmp_id_B,1);
    output = output + tmp_s;
end

end

%for test 时间已过 0.172028 秒。
%cellfun test1 时间已过 0.147736 秒。
%cellfun test2 时间已过 166.449052 秒。

总结: 可以看到使用cellfun如果用的不好,可以使得程序运行速度变得极慢. 这里的例子只是进行了简单的测试, 实际上for循环内容如果越发复杂,cellfun的提速效果会更加明显.

第八种提速方式:用外部语言编程来加速核心算法

8.1 Mex 混合编程

Matlab支持外部语言来进行混合编程, 也就是Mex混合编程. 支持C/C++, Java, Python, .NET 和COM对象的混合编程.

采用Mex混合编程的Blog教程不少. 从我个人使用过得Mex混编经验看, 强烈建议官网上的使用手册,总结的很系统而且也可以一步步上手. 如果大家想要系统地掌握Mex编程,最好还是详细去阅读Mex官网文档.

具体如何使用Mex加速,这里我就不详细评估和测试了,如果大家希望我帮忙整理系统的Mex使用说明和加上自己的经验以及注意事项,可以反馈给我. 我会看情况来决定是否写一篇Matlab效率之Mex混合编程的文章.

8.2 system调用外部命令

Matlab的system命令可以直接调用系统中的cmd命令. 因此,如果为了加速代码运行速度, 可以在各自的系统平台里进行更底层的编程开发,比如利用C++进行编程.

说到这里, 大家可能会很好奇, Mex和system命令调用的速度比较情况. 从阿Q我个人的经验来看, 可以比较负责任的建议采用system命令会更好.

以C++ 编程为例, 原因如下: 1) 直接使用IDL集成开发平台进行C++代码下的算法开发,可以直接先开发成库和可执行文件, 方便后续使用C++下进一步基于原先写好的代码来一步步搭建出相关的应用.
使用方式:
(1) 使用system命令在Matlab中混合调用, 需要在C++上给测试算法写好文件读写的接口. 以文件操作为媒介,来让Matlab来调用外部开发的功能模块.
(2) 如果需要统计算法时间统计,可以直接在实际的C++算法部分进行时间统计.

编程习惯的问题

在matlab代码编写的时候经常会遇到多层for循环嵌套的情况,在编写算法之初,可以先不考虑加速的问题,因为一旦加速之后代码的可读性会急剧下降.

建议: 那么在一开始编写的时候遇到for循环多层嵌套, 可以好好思考下循环的设计.

尽量每一层for循环的迭代次数上, 从内向外看, 要尽可能的由大到小. 即最内的循环次数可以多, 要处理部分放在内层循环.
示例代码:


m = 10;
n = 1000;
size_A = [ 100,100];

num_test = 10; 
%-----------------------------------------------------
tic;
for k = 1:num_test
    for i = 1:m
        for j = 1:n
            A = rand(size_A);
            min_svd_value = svds(A,1,'smallest');
        end
    end
end

time_1 = toc / num_test;
tmpstr = sprintf('time cost for test 1: %.3f',time_1);
disp(tmpstr);

%-----------------------------------------------------
tic;
for k = 1:num_test
    for i = 1:n
        for j = 1:m
            A = rand(size_A);
            min_svd_value = svds(A,1,'smallest');
        end
    end
end
time_2 = toc / num_test;
tmpstr = sprintf('time cost for test 2: %.3f',time_2);
disp(tmpstr);
%-----------------------------------------------------
%time cost for test 1: 51.169
%time cost for test 2: 61.210

总结:虽然测试的规模较小,而且for循环内部代码比较简单, 不过从结果也可以看出来从内到外,好的编程设计应该是循环次数依次递减的. 换言之, 最外层的迭代次数最小,最内层迭代次数最大.

看到这里,给大家一个总结吧. 总而言之, Matlab代码效率上加速主要还是得靠向量化, 各类加速方式可以参考Blog最前面的列表. 至于或许还有其他加速策略或者是利于加速的编程习惯,我暂时还没有想起来. 后面持续更新的时候,如果某天突然想起来了,会跟大家进行分享.

在这里,对所有阅读关注和点收藏的朋友,阿Q表示衷心的感谢~ 你们的支持是我更新的动力.暂时先发布这部分内容,后续会进行持续更新. 如有书写错误的地方也请指正,我会更正的.

这里说个题外话,打个广告. 在代码加速这块, 导师和课题组的代码都是我负责加速的. 以我们公司为例, 在我们发表的论文里,
[TAES]Yuanxin Wu*, Qi Cai and Tnieu-Kien Truong, “Fast RodFIter for Attitude Reconstruction from Inertial Measurements”, IEEE Trans. on Aerospace and Electronic Systems(TAES), vol. 55, no. 1, pp. 419-428, 2019.
以及
[IJCV]. Qi Cai, Yuanxin Wu*, Lilian Zhang and Peike Zhang, “Equivalent Constraints for Two-View Geometry: Pose Solution/Pure Rotation IYuanxin Wu*dentification and 3D Reconstruction”, International Journal of Computer Vision(IJCV), vol. 127, no. 2, pp. 163-180, 2019.

代码找了其他人(公司)加速不起来,后来让我去加速,从原来的版本加速了超过200多倍, 那么大家在写论文的时候就可以在时间性能上进行更好的展示.我们自己刚建立的公司支持为高校师生提供Matlab代码加速服务, 后续考虑c/c++和Python. 暂时在咸鱼上试拓展业务,后面反馈不错的话,就开始转到公司的拓展业务上. 报表和保密协议在咸鱼上有说明.