Vistb's Tiny IT Space

Am I a Geek?

导航

从Matlab的for循环说开……

Posted on 2010-01-19 01:19  Vistb  阅读(37860)  评论(8编辑  收藏  举报

  因为学习和工作的原因,最近又开始使用已经许久没有接触的Matlab。在没有什么特殊考虑的情况下,信手写下了下面的m代码片段:

 1 for i=1:1:(imgHeight-tmpHeight+1)
 2     for j=1:1:(imgWidth-tmpWidth+1)
 3         temp=0;
 4         for m=1:1:tmpHeight
 5             for n=1:1:tmpWidth
 6                 temp=temp+img(i+m-1,j+n-1)*template(m,n);
 7             end
 8         end
 9         if temp>0
10             tmpRst(i+floor(tmpHeight/2),j+floor(tmpWidth/2))=temp;
11         end;
12     end
13 end

  

  外层循环的2个变量长度为300和400,内层的两个为9。出乎我的意料的是,这样一段代码在我的机器上(T5750@2GHz, 2GB DDRII667)竟然要跑1分多钟,而这段代码转换为C++后是准备要在一个实时图像识别系统上跑的。换言之,时间至少必须下降到1/25以内!虽然可以指望C++的效率,但Matlab这样的速度也太离谱了!况且我在Matlab中还要不断的实验,跑一遍就要1、2分钟,实在令人难以接受。下午和师兄们讨论时,无意谈到了这个问题,他们立即对我说,想办法转换为矩阵操作什么的,Matlab中for循环的效率是很低的!

   回寝室后,试验了一下,结果令人乍舌,我把代码改成了下面的样子(注意第三行代码实际上代替了内层for循环,其他的改动我想没什么本质影响):

1 for i=1:1:(imgHeight-tmpHeight+1)
2     for j=1:1:(imgWidth-tmpWidth+1)
3         temp=img(i:(i+tmpHeight-1),j:(j+tmpWidth-1)).*template;
4         temp=sum(sum(temp));
5         tmpRst(i+floor(tmpHeight/2),j+floor(tmpWidth/2))=(temp+abs(temp))/2;
6     end
7 end

  

  结果时间呢,只要了2秒左右!这么一改,效率提升了几十倍!我是学过一些编译原理的,但这种差距实在令我感到很不解。求助Google后,算是得到了满意的答复。


  首先,看看Improving the Speed of MATLAB Calculations这篇文章中怎么说的。

  在前言中,文章中有两段话:

   MATLAB programs are interpretted. This would seem to make it inapproapriate for large scale scientific computing. The power of MATLAB is realized with its extensive set of libraries which are compiled or are carefully coded in MATLAB to utilize "vectorization". The concept of vectorization is central to understanding how to write efficient MATLAB code.
  Vectorized code takes advantage, wherever possible, of operations involving data stored as vectors. This even applies to matrices since a MATLAB matrix is stored (by columns) in contiguous locations in the computer's RAM. The speed of a numerical algorithm in MATLAB is very sensitive to whether or not vectorized operations are used.

  其核心大意就是说为了弥补Matlab程序是解释执行所带来效率问题,我们应该尽量使用“向量化”(vectorization)的命令。Matlab程序执行的效率,对于是否使用了“向量化”命令是非常(very)敏感的!

  其后,文章给出了两条实用的建议。

  第一条,使用向量操作代替循环。 以下举例说明。

1 dx = pi/30;
2 nx = 1 + 2*pi/dx;
3 for i = 1:nx
4   x(i) = (i-1)*dx;
5   y(i) = sin(3*x(i));
6 end


  这段代码是很自然的从C语言的形式转化而来的,但其效率很低!Matlab是实时为变量分配内存的,在第一遍循环时(即i=1时),Matlab为x和y这两个向量(长度均为1)分配内存。以后每执行一次循环,Matlab都会在x和y的末尾附加新的元素。这不仅导致分配内存的调用的增加,也使得x和y的各个元素在内存中的分布不是连续的(就像数据结构中数组和链表的区别)。由此,性能遭到了损失。

  相比之下,下面的代码效率提高不少:

1 = 0:pi/30:2*pi
2 = sin(3*x);


  第一个语句分配了一个连续的内存空间来存储具有多个元素的向量x。类似的,第二个语句在分配内存时,也是分配了一个连续的内存空间来存储具有多个元素的向量y。撇去计算sin的消耗不算,就内存分配命令的执行次数和对向量元素访问的方便程度来说,高下立见。

  第二条,为矩阵和向量预先分配内存。

  文章中指出,虽然Matlab会自动调整变量的大小,我们最好还是预先为变量分配内存空间。因为这样可以使调用内存分配命令的次数降为1,也可以使变量在内存中连续存储(当变量为矩阵时是按列在内存中连续存储)。

  而所谓“预先为变量分配内存空间” ,是指在知道变量的大小的情况下,在变量中的任何一个元素都未被引用之前,创建一个大小和其一致的变量。

  下面是一个例子,代码质量从上至下逐渐提高:

 1 dx = pi/30;
 2 nx = 1 + 2*pi/dx;
 3 nx2 = nx/2;
 4 
 5 for i = 1:nx2
 6   x(i) = (i-1)*dx;
 7   y(i) = sin(3*x(i));
 8 end
 9 
10 for i = nx2+1:nx
11   x(i) = (i-1)*dx;
12   y(i) = sin(5*x(i));
13 end

 

 1 dx = pi/30;
 2 nx = 1 + 2*pi/dx;
 3 nx2 = nx/2;
 4 
 5 = zeros(1,nx);      % 为向量x预分配内存
 6 = zeros(1,nx);      % 为向量y预分配内存
 7 
 8 for i = 1:nx2
 9   x(i) = (i-1)*dx;
10   y(i) = sin(3*x(i));
11 end
12 
13 for i = nx2+1:nx
14   x(i) = (i-1)*dx;
15   y(i) = sin(5*x(i));
16 end


 1 = 0:pi/30:2*pi;     % 计算向量x的值
 2 nx = length(x);
 3 nx2 = nx/2;
 4 
 5 = x;                % 为向量y预分配内存
 6 
 7 for i = 1:nx2
 8   y(i) = sin(3*x(i));
 9 end
10 
11 for i = nx2+1:nx
12   y(i) = sin(5*x(i));
13 end


1 = 0:pi/30:2*pi;                  % 计算向量x的值
2 nx = length(x);
3 nx2 = nx/2;
4 
5 = x;                             % 为向量y预分配内存
6 
7 y(1:nx2) = sin(3*x(1:nx2));        % 计算y的第1部分的值
8 y(nx2+1:nx) = sin(5*x(nx2+1:nx));  % 计算y的第2部分的值



  然后,再来看看这个名为Improving Program Efficiency的PPT。

  除了上篇文章提到的那几点以外,该ppt中还提出了以下几点看法和建议。

  第一,选择合适的数据类型。Matlab有多种数据类型,不同的数据类型可以带来不同的精度,但处理速度也存在差别。double当然可以比int8带来更高的精度,但性能却会下降。不过,我个人对这个建议持保留意见,主要在于有些操作对一些诸如int8类型的非标准类型不支持,而且有时候容易产生误操作(例如相对uint8这样的无符号整形变量)。

  第二,使用tic和toc来测试程序的执行时间。这两个命令配合使用可以测试一段m代码的执行时间。具体的就不多说了,大家可以去查看Matlab的帮助文件。另外,Matlab最近的版本(像R2009b)中出现了类似于性能测试工具的组件,大家可以在Matlab的帮助文件中搜索"Profiling for Improving Performance" 进行进一步了解。

  第三,类似于上一篇文章中提到的使用向量化命令,减少循环。但是,该ppt中还列出了一些常用的可以用来代替循环的向量化命令,列举如下:

  • find        (find values that meet some criteria,寻找符合某些特定条件的矩阵中的元素)
  • sum, prod, diff   (sum 加, product 乘, difference 减)
  • .*, ./        (element by element matrix operations,矩阵间逐元素操作)
  • min, max       (find min or max values,求最小和最大值)
  • zeros, ones     (for initializing arrays,用于初始化变量)

  其中,我觉得find、prod、diff等都是比较少见的(可能由于我才疏学浅,呵呵),大家可以仔细研究一下。尤其是find,非常有用!


  最后,在一篇讨论for循环的效率的帖子中,我也听到了一些不同的意见。

  总的来说,就是对for循环不能一棒子打死,要区别对待(像对goto语句?)。

  在这个帖子中,名为Bruno Luong的作者总结道(本人英文不好,不敢打包票翻译对了,故附上原文~~~):

  • if there is an equivalent vectorized stock function, always use it(如果有等价的向量化命令,毫不犹豫的使用后者)
  • avoid for-loop that call function with non negligible overhead(避免在for循环中调用计算量很大的函数)
  • for loop is desirable when a nested IF condition can be used to save computation time(如果for循环中有if语句,并可以因此而带来时间的节省,那么for语句是值得试试的)
  • for loop is attractive when the result of the preceding iteration(s) can be used to save computation effort of the current calculation(如果上一次的循环得出的结果对本次循环有帮助,可以节省计算量,那么for循环是比较吸引人的)
  • using for loop is not recommended when the large data need to be duplicated inside the loop(当循环中存在大量的数据复制时,for循环是不值得推荐使用的)
  • time it, time it and time it(不停的测试,测试不用for和用for的区别,“唯利是图”就可以了)
  • Read Matt Fig's post!(阅读Matt Fig的帖子!【Matt Fig是哪位大神?我没Google,大家看看他发表了什么高见,然后告诉我一声啊~】)


  哇,写的累死了~~~就此打住吧,时间很晚了,明天还要去看iMax的Avatar(虽然是第三遍,不过依旧很高兴)~~~