MATLAB 效率再议

关于MATLAB的效率问题，很多文章，包括我之前写的一些，主要集中在使用向量化以及相关的问题上。但是，最近我在实验时对代码进行profile的过程中，发现在新版本的MATLAB下，for-loop已经得到了极大优化，而效率的瓶颈更多是在函数调用和索引访问的过程中。

由于MATLAB特有的解释过程，不同方式的函数调用和元素索引，其效率差别巨大。不恰当的使用方式可能在本来不起眼的地方带来严重的开销，甚至可能使你的代码的运行时间增加上千倍（这就不是多买几台服务器或者增加计算节点能解决的了，呵呵）。

下面通过一些简单例子说明问题。（实验选在装有Windows Vista的一台普通的台式机(Core2 Duo 2.33GHz + 4GB Ram)进行，相比于计算集群, 这可能和大部分朋友的环境更相似一些。实验过程是对某一个过程实施多次的整体进行计时，然后得到每次过程的平均时间，以减少计时误差带来的影响。多次实验，在均值附近正负20%的范围内的置信度高于95%。为了避免算上首次运行时花在预编译上的时间，在开始计时前都进行充分的“热身”运行。）

函数调用的效率

一个非常简单的例子，把向量中每个元素加1。（当然这个例子根本不需要调函数，但是，用它主要是为了减少函数执行本身的时间，突出函数解析和调用的过程。）

作为baseline，先看看最直接的实现

% input u: u is a 1000000 x 1 vector
v = u + 1; 

这个过程平均需要0.00105 sec。

而使用长期被要求尽量避免的for-loop

n = numel(u);
% v = zeros(n, 1) has been pre-allocated.
for i = 1 : n
    v(i) = u(i) + 1;
end

所需的平均时间大概是0.00110 sec。从统计意义上说，和vectorization已经没有显著差别。无论是for-loop或者vectorization，每秒平均进行约10亿次“索引-读取-加法-写入”的过程，计算资源应该得到了比较充分的利用。

要是这个过程使用了函数调用呢？MATLAB里面支持很多种函数调用方式，主要的有m-function, function handle, anonymous function, inline, 和feval，而feval的主参数可以是字符串名字，function handle, anonymous function或者inline。

用m-function，就是专门定义一个函数

function y = fm(x)
    y = x + 1;

在调用时

for i = 1 : n
    v(i) = fm(u(i));
end

function handle就是用@来把一个function赋到一个变量上，类似于C/C++的函数指针，或者C#里面的delegate的作用

fh = @fm;
for i = 1 : n
    v(i) = fh(u(i));
end

anonymous function是一种便捷的语法来构造简单的函数，类似于LISP, Python的lambda表达式

fa = @(x) x + 1;
for i = 1 : n
    v(i) = fa(u(i));
end

inline function是一种传统的通过表达式字符串构造函数的过程

fi = inline('x + 1', 'x');
for i = 1 : n
    v(i) = fi(u(i));
end

feval的好处在于可以以字符串方式指定名字来调用函数，当然它也可以接受别的参数。

v(i) = feval('fm', u(i));
v(i) = feval(fh, u(i));
v(i) = feval(fa, u(i)); 

对于100万次调用（包含for-loop本身的开销，函数解析(resolution)，压栈，执行加法，退栈，把返回值赋给接收变量），不同的方式的时间差别很大：

m-function
0.385 sec

function handle
0.615 sec

anonymous function
0.635 sec

inline function
166.00 sec

feval('fm', u(i))
8.328 sec

feval(fh, u(i))
0.618 sec

feval(fa, u(i))
0.652 sec

feval(@fm, u(i))
2.788 sec

feval(@fa, u(i))
34.689 sec

从这里面，我们可以看到几个有意思的现象：

• 首先，调用自定义函数的开销远远高于一个简单的计算过程。直接写 u(i) = v(i) + 1 只需要  0.0011 秒左右，而写u(i) = fm(v(i)) 则需要0.385秒，时间多了几百倍，而它们干的是同一件事情。这说明了，函数调用的开销远远大于for-loop自己的开销和简单计算过程——在不同情况可能有点差别，一般而言，一次普通的函数调用花费的时间相当于进行了几百次到一两千次双精度浮点加法。
• 使用function handle和anonymous function的开销比使用普通的m-函数要高一些。这可能是由于涉及到句柄解析的时间，而普通的函数在第一次运行前已经在内存中进行预编译。
• inline function的运行时间则令人难以接受了，竟然需要一百多秒（是普通函数调用的四百多倍，是直接计算的十几万倍）。这是因为matlab是在每次运行时临时对字符串表达式（或者它的某种不太高效的内部表达）进行parse。
• feval(fh, u(i))和fh(u(i))，feval(fa, u(i))和fa(u(i))的运行时间很接近，表明feval在接受句柄为主参数时本身开销很小。但是，surprising的是

  for i = 1 : n
      v(i) = feval(@fm, u(i));
  end
  比起
  fh = @fm;
  for i = 1 : n
      v(i) = feval(fh, u(i));
  end
  慢了4.5倍 (前者0.618秒，后者2.788秒)。
  

  for i = 1 : n
      v(i) = feval(@(x) x + 1, u(i));
  end
  比起
  fa = @(x) x + 1;
  for i = 1 : n
      v(i) = feval(fa, u(i));
  end
  竟然慢了53倍(前者0.652秒，后者34.689秒）。
  

  由于在MATLAB的内部实现中，function handle的解析是在赋值过程中进行的，所以预先用一个变量把句柄接下，其效果就是预先完成了句柄解析，而如果直接把@fm或者@(x) x + 1写在参数列上，虽然写法简洁一些，但是解析过程是把参数被赋值到所调函数的局部变量时才进行，每调用一次解析一次，造成了巨大的开销。

• feval使用字符串作为函数名字时，所耗时间比传入句柄大，因为这涉及到对函数进行搜索的时间（当然这个搜索是在一个索引好的cache里面进行（除了第一次），而不是在所有path指定的路径中搜索。）

在2007年以后，MATLAB推出了arrayfun函数，上面的for-loop可以写为

v = arrayfun(fh, u)

这平均需要4.48 sec，这比起for-loop(需时0.615 sec)还慢了7倍多。这个看上去“消除了for-loop"的函数，由于其内部设计的原因，未必能带来效率上的正效果。

元素和域的访问

除了函数调用，数据的访问方式对于效率也有很大影响。MATLAB主要支持下面一些形式的访问：

• array-index  A(i): 
• cell-index:  C{i}; 
• struct field:  S.fieldname
• struct field (dynamic):  S.('fieldname')

这里主要探索单个元素或者域的访问（当然，MATLAB也支持对于子数组的非常灵活整体索引）。

对于一百万次访问的平均时间

A(i) for a numeric array
0.0052 sec

C{i} for a cell array
0.2568 sec

struct field
0.0045 sec

struct field (with dynamic name)
1.0394 sec

我们可以看到MATLAB对于单个数组元素或者静态的struct field的访问，可以达到不错的速度，在主流台式机约每秒2亿次（连同for-loop的开销）。而cell array的访问则明显缓慢，约每秒400万次（慢了50倍）。MATLAB还支持灵活的使用字符串来指定要访问域的语法（动态名字），但是，是以巨大的开销为代价的，比起静态的访问慢了200倍以上。

关于Object-oriented Programming

MATLAB在新的版本中（尤其是2008版），对于面向对象的编程提供了强大的支持。在2008a中，它对于OO的支持已经不亚于python等的高级脚本语言。但是，我在实验中看到，虽然在语法上提供了全面的支持，但是matlab里面面向对象的效率很低，开销巨大。这里仅举几个例子。

• object中的property的访问速度是3500万次，比struct field慢了6-8倍。MATLAB提供了一种叫做dependent property的属性，非常灵活，但是，效率更低，平均每秒访问速度竟然低至2.6万次（这种速度基本甚至难以用于中小规模的应用中）。
• object中method调用的效率也明显低于普通函数的调用，对于instance method，每百万次调用，平均需时5.8秒，而对于static method，每百万次调用需时25.8秒。这相当于每次调用都需要临时解析的速度，而matlab的类方法解析的效率目前还明显偏低。
• MATLAB中可以通过改写subsref和subsasgn的方法，对于对象的索引和域的访问进行非常灵活的改造，可以通过它创建类似于数组的对象。但是，一个符合要求的subsref的行为比较复杂。在一个提供了subsref的对象中，大部分行为都需要subsref进行调度，而默认的较优的调度方式将被掩盖。在一个提供了subsref的类中（使用一种最快捷的实现），object property的平均访问速度竟然降到1万5千次每秒。

建议

根据上面的分析，对于撰写高效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的接口调用，往往会带来非常巨大的开销（成千上万倍的减慢），甚至使得本来几乎不是问题的代码成为性能瓶颈。

这是一篇转载的文章：http://dahua.spaces.live.com/blog/cns!28AF4251DF30CA42!2459.entry#trackback