Matlab高性能编程——代码优化和并行计算

Jeremy Lin  @HQU
Update: 2014/4/29 

Matlab代码优化

Matlab是一种高级计算机语言，同时也是一个用于算法开发，数据可视化，数据分析和数值计算的交互式工作环境。尽管Matlab软件提供了大量专业化的工具箱，使用户避免了很多编程工作，但是在实际工作中仍不免需要自行编写Matlab代码以应对各种纷繁复杂的应用。我们需要明白Matlab是一种专门为数组运算而设计的语言，因此在程序设计中要注意充分利用这一优点来加快运算速度。

本文从以下三个方面展开：

程序语法分析 ：Code Analyzer（M-Lint）
程序性能分析 ：Profiler
运算性能提升 ：向量化、预分配

1.程序语法分析

关于Matlab的程序语法分析，我们要充分利用Code Analyzer的各种提示，修改出现的各种warning，来改善程序。所谓的Code Analyzer就是在m文件编辑器右边的一栏warning提示，当程序出现warning时，它附带地提供了一些可参考的方法。
在下图中，我们可见到两个warning，一个是没有分号结束符，这样会使代码结果直接输出，有过编程经验的人都知道，在代码的运行过程中不断输出中间结果，会大大增加程序的运行时间；另一个warning是没有对其中的变量进行初始化，这样会让程序在运行过程中不断重复初始化变量，同时也会大大消耗内存。



一个简单的例子如下所示：



当我们没有去掉warning时，程序的运行时间是：5.570289s，而当我们添加了代码的分号结束符后，程序的运行时间是：0.008179s，对比发现程序提速非常明显。

2.程序性能分析

对程序的性能分析我们主要利用Profiler来深度了解代码的性能，找出瓶颈代码，然后做出修改以提升性能。
在Matlab中打开Profiler只要点击Run and Time即可，位置如下图所示：



我们从Profiler可以得到

函数调用的总次数
每次函数调用耗时状况

然后根据各个函数的调用情况，结合实际情况，把可以替换的函数直接替换，不能替换的高调用函数尽量优化。

3.运算性能提升

运算性能提升主要利用两种方法：内存预分配和向量化。
首先我们来看看内存预分配的必要性：
需要明白，当你未进行内存预分配时，比如执行如下的代码：

>> x=4;
>> x(2)=7;
>> x(3)=12;

它在内存中的过程是这样的：



从上图我们可以发现，改变数组大小是很耗费内存的，因此最好进行内存预分配，同时它能够减少存储器碎片。

接下来，我们来看看向量化，所谓向量化即是将for循环和while循环转换为等价的向量或矩阵运算，它可以大大加速运算，还增强了代码的可读性。

我们从一个简单的例子开始：

原程序：

A=1;
for x=1:10000000
f(x)=A*sin((x-1)/(2*pi));
end
运行时间：12.714 s

向量化后：

<span style="font-size:14px;">A=1;
x=0:10000000-1;
f=A*sin(x/(2*pi));
运行时间：0.544 s</span>

可以发现，向量化后，程序的运行速度提高了23倍左右。

虽然一维的向量化很简单，但是，当被评估的函数有两个变量时，优化的方法可能就会复杂一些。

MATLAB使用meshgrid来实现对二维函数的评估，该函数的语法为

[C, R] = meshgrid(c, r)

该函数将由行向量c和r指定的域变换成数组C和R，这两个数组能用来评估有着两个变量的函数和三维表面图（注意，在meshgrid的输入和输出中，列总是首先列出）。

输出数组C是向量c的副本，R是r的副本。例如，假设我们想形成一个二维函数，该函数的元素是坐标向量x和y的值的平方和，其中x=0, 1, 2; y=0, 1。向量r由坐标的行向量构成：r=[0, 1, 2]；类似的，c由坐标的列向量构成：c=[0 1]（注意，此处的r和c均为行向量）。将这两个向量代入meshgrid可得到如下数组：

>> [C,R]=meshgrid(c,r)
C =
0     1
0     1
0     1
R =
0     0
1     1
2     2

我们感兴趣这个函数的实现：

>> h=R.^2+C.^2

h =
0     1
1     2
4     5

例子：

tic
u0=1/(4*pi);
v0=1/(4*pi);
for r=1:1000
u0x=u0*(r-1);
for c=1:1000
v0y=v0*(c-1);
f(r,c)=sin(u0x+v0y);
end
end
t1=toc

运行时间： 2.4337 s

向量化后：

tic
u0=1/(4*pi);
v0=1/(4*pi);
r=0:1000-1;
c=0:1000-1;
[C,R] = meshgrid(c,r);
g=sin(u0*R+v0*C);
t2=toc
运行时间： 0.0590s

向量后快了41倍。

OTHER TIPS:

以列向量存储：以列作为双重for循环的外循环会比以行作为外循环运算速度更快
逻辑索引运算性能更好：例子如下

%%
N=2000;
A=magic(N);
A1=magic(N);
A2=magic(N);
myRef=1e6;
%%
tic
ix=1;
vals=zeros(size(A(:)));
for jj=1:N
for ii=1:N
if A(ii,jj)>myRef
vals(ix) = A(ii,jj);
ix=ix+1;
end
end
end

% vals(ix:end)=[];
toc

%%

运行时间：Elapsed time is 5.466597 seconds.

tic
vals=A2(A2>myRef);
toc

运行时间：Elapsed time is 0.084450 seconds.

In-place操作 ：减少临时变量的使用

N=3e3;
x=rand(N);
tic;
y=x*1.2;
toc;
%Elapsed time is 0.073008 seconds.
%% In-place 操作
tic;
x=x*1.2;
toc;
%Elapsed time is 0.037856 seconds.

MATLAB并行计算

关于Matlab并行计算，这一块我实际的经验较少，只做一点浅显的介绍，更多相关的资料见下文的参考资料。

1、Matlab并行计算原理

Matlab的并行计算实质还是主从结构的分布式计算。当你初始化Matlab并行计算环境时，你最初的Matlab进程自动成为主节点，同时初始化多个Matlab计算子节点。Parfor的作用就是让这些子节点同时运行Parfor语句段中的代码。Parfor运行之初，主节点会将Parfor循环程序之外变量传递给计算子节点。子节点运算过程时互不干扰，运算完毕，则应该有相应代码将各子节点得到的结果组合到同一个数组变量中，并返回到Matlab主节点。当然，最终计算完毕应该手动关闭计算子节点。

2、初始化Matlab并行计算环境

这里讲述的方法仅针对多核机器做并行计算的情况。设机器的CPU核心数量是CoreNum双核机器的CoreNum2，依次类推。CoreNum以不等于核心数量，但是如果CoreNum小于核心数量则核心利用率没有最大化，如果CoreNum大于核心数量则效率反而可能下降。因此单核机器就不要折腾并行计算了，否则速度还更慢。

下面一段代码初始化Matlab并行计算环境：

%Initialize Matlab Parallel Computing Enviornment by Xaero | Macro2.cn
CoreNum=2; %设定机器CPU核心数量，我的机器是双核，所以CoreNum=2
if matlabpool('size')<=0 %判断并行计算环境是否已然启动
matlabpool('open','local',CoreNum); %若尚未启动，则启动并行环境
else
disp('Already initialized'); %说明并行环境已经启动。
end

运行成功后会出现如下语句：

Starting matlabpool using the 'local' configuration ... connected to 2 labs.

如果运行出错，按照下面的办法检测：

首先运行：

matlabpool size

如果出错，说明你没有安装Matlab并行工具箱。确认安装了此工具箱后，运行：

matlabpool open local 2;

如果出错，证明你的机器在开启并行计算时设置有问题。

3、终止Matlab并行计算环境

用上述语句启动Matlab并行计算环境的话，在你的内存里面有CoreNum个Matlab进程存在，每个占用内存都在百兆以上。（可以用Windows任务管理器查看），故完成运行计算后可以将其关闭。关闭的命令很简单：

matlabpool close

Reference：

[1]  Rafael C. Gonzalez. Digital Image Processing Using MATLAB.

[2] mathworks  webinars

[3] Matlab Help

[4] Matlab并行编程方法 Rachel Zhang的专栏

本文地址：http://blog.csdn.net/linj_m/article/details/9730717
更多资源请关注 博客：LinJM-机器视觉
 微博：林建民-机器视觉