MATLAB的效率杂谈(2007-7)
在我们的领域，MATLAB是进行实验的最主要的工具之一。我听到过很多抱怨说MATLAB很慢，可是这往往不是MATLAB的问题，而是因为自己程序没有写好。MATLAB是我非常欣赏的一种语言，堪称灵活易用和高效运算的结合典范。
作为解释语言，MATLAB进行高效运算的秘诀有几个方面：
(1) JIT即时编译技术。Matlab在第一次加载运行一个函数的时候，会在后台对它进行某种程度的编译和优化，使得后续运行更快。因此，除了第一次运行需要进行语句解释之外，后面运行的其实是已经放在内存中的经过优化的中间码。所以，很多时候我们可能会看到第一次运行一个新函数，比它后面几次的运行明显慢一些。不过目前的JIT技术还不是非常成熟，和标准的编译语言相比还有相当差距，仅凭这个MATLAB还不能称为高效。
(2) 向量化(Vectorization)。这是MATLAB最著名的特色。向量化配合经过高度优化的数值运算引擎，是其高效运算的最重要的基石——很多MATLAB的使用者都明白这一点。能够转化成矩阵操作的规则运算过程，使用MATLAB计算远比自己手工用C/C++实现高效。我自己做过很多次对比profile，MATLAB在关键的核心运算上的实现可以比自己用C/C++按照标准的routine进行实现快几十倍。
其实这不完全是MATLAB的功劳，其实MATLAB是建基于BLAS和LAPACK等的基础数学运算包的基础上的。Intel和AMD都发布了这些东西的vendor-implementation，并且针对各自的CPU进行了大量的优化，这非个人之力所能企及。因此，我从来都强烈不建议个人用C/C++去实现数值运算的关键代码（学习数值分析课者除外）。
不过，对于向量化这个事情，也不宜极端化，下面的一些例子说明在哪些时候，for-loop比vectorization更适合：
- (a) 粗粒度的算法流程控制。比如，你要实现一个迭代算法，循环做一个事情直到收敛。只要循环几次或者几十次，但是每个循环内部要进行大量的复杂运算，那么你就没有必要费心把这层循环给vectorize掉了。除非收敛结果有某种close-form的解析解。
- (b) 如果向量化可能导致产生巨型矩阵，则使用前要三思！很多情况下，向量化是一种用空间换时间的行为，就是通过把数据组织成某种方式，从而使得内建的高效引擎能得以应用。但是，有些时候要处理大量的数据，可能导致其组织过程需要耗费额外的数百兆乃至千兆内存空间，那么这有可能造成效率的不升反降。原因有几个方面：
第一，数据组织过程也是需要时间的，最起码它也需要大量的操作进行密集的内存读写和用于定位的偏移量计算。这方面增加的时间 vs. 向量化节省的运算时间，孰轻孰重，需要衡量。
第二，分配额外的大块内存是一件非常耗时的事情，它可能导致虚拟内存的使用，那么当对这块矩阵进行读写和计算时可能涉及频繁的内存与外存交换区的I/O，这会造成效率的严重下降。我一直不赞同在Out of Memory的情况下，通过增大虚存来解决问题，这样，即使你勉强让程序能继续运行下去，运行时间也会变得极为缓慢——这时应该做的是对程序进行重新思考和设计，降低其对内存的耗用。
第三，vectorization有些时候还可能增加实际运算次数，这往往出现在不适合的向量化过程中。这样，即使你通过生搬硬套的向量化过程让操作变成矩阵运算，但是增加的无用计算使得即使是更高效的引擎也无法挽回你的损失了。
说了这么些，不是想劝说从向量化的道路退回去，而是提醒在做一些事情的时候，要考虑得周全一些。
(3) Inplace Operation。这是一个很有趣的事情，MATLAB的对象管理策略是Copy-on-Write，就是说你把一个矩阵传给一个函数的时候，会先传引用，而不产生副本，而当函数要对这个矩阵修改的时候，才会制造出它的副本出来，让函数去修改。这样听上去很完美，既保护了输入矩阵不被改变，又避免了大量无谓的复制。不过，这个策略真的没有缺陷么？可以看看下面的例子
A = rand(2000, 2000);
for i = 1 : 1000
A = f(A);
end
function A = f(A)
% a temporary copy of A is produced for modification
A(1) = A(1) + 1;
% the modified temporary copy is assigned to A
return;
在上面的代码中，只是想对A的第一个元素做1000次加法，结果却导致了对整个有四百万元素的大矩阵做了1000次副本复制！而且这些副本其实没有用，只要f(A)直接对A的原矩阵进行修改，这些巨大的浪费就能避免。你可能跟我argue说，这里完全可以写成A(1) = A(1) + 1000，不用这么搞。我想说明的是，我只是想举一个简单例子说明，Copy-on-write的策略为什么可能造成巨大的效率浪费。其实，很多小修改没法像上面那么容易合并。
以前，为了解决这个问题，f函数的作者需要自己写一个mex function来避开MATLAB的保护机制，直接改写A。后来matlab自己也意识到这种策略在实际中的问题，于是他们改写了解释器，采用了一种办法：A = f(A, ...) 当它发现这样的定义和调用的时候，它会很智能地了解到你其实是想改写A这个输入，于是它就把操作改成inplace，直接对A进行，避免拷贝。这种修改的策略，在2006b后才比较正式的运用，在旧版Matlab中仍旧不同程度地存在上述问题。我在2007a中，测试了两种f的写法：
function y = f0(x) y = x + 1; end
function x = f1(x) x = x + 1; end
确实发现大量调用f1比f0快了好多倍。所以，如果确实是想改变input的话，函数定义中要让input parameter和output parameter同名，以触发这种智能的解释机制。并且建议升级到最新的2007a版本，这个版本在这方面进行了重要改善。
MATLAB 效率再议(2008-6)
关于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。
1)用m-function，就是专门定义一个函数
function y = fm(x)
y = x + 1;
在调用时
for i = 1 : n
v(i) = fm(u(i));
end
2)function handle就是用@来把一个function赋到一个变量上，类似于C/C++的函数指针，或者C#里面的delegate的作用
fh = @fm;
for i = 1 : n
v(i) = fh(u(i));
end
3)anonymous function是一种便捷的语法来构造简单的函数，类似于LISP, Python的lambda表达式
fa = @(x) x + 1;
for i = 1 : n
v(i) = fa(u(i));
end
4)inline function是一种传统的通过表达式字符串构造函数的过程
fi = inline('x + 1', 'x');
for i = 1 : n
v(i) = fi(u(i));
end
5)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代码，我有下面一些建议：
- 虽然for-loop的速度有了很大改善，vectorization（向量化）仍旧是改善效率的重要途径，尤其是在能把运算改写成矩阵乘法的情况下，改善尤为显著。
- 在不少情况下，for-loop本身已经不构成太大问题，尤其是当循环体本身需要较多的计算的时候。这个时候，改善概率的关键在于改善循环体本身而不是去掉for-loop。
- MATLAB的函数调用过程（非built-in function）有显著开销，因此，在效率要求较高的代码中，应该尽可能采用扁平的调用结构，也就是在保持代码清晰和可维护的情况下，尽量直接写表达式和利用built-in function，避免不必要的自定义函数调用过程。在次数很多的循环体内（包括在cellfun, arrayfun等实际上蕴含循环的函数）形成长调用链，会带来很大的开销。
- 在调用函数时，首选built-in function，然后是普通的m-file函数，然后才是function handle或者anonymous function。在使用function handle或者anonymous function作为参数传递时，如果该函数被调用多次，最好先用一个变量接住，再传入该变量。这样，可以有效避免重复的解析过程。
- 在可能的情况下，使用numeric array或者struct array，它们的效率大幅度高于cell array（几十倍甚至更多）。对于struct，尽可能使用普通的域（字段，field）访问方式，在非效率关键，执行次数较少，而灵活性要求较高的代码中，可以考虑使用动态名称的域访问。
- 虽然object-oriented从软件工程的角度更为优胜，而且object的使用很多时候很方便，但是MATLAB目前对于OO的实现效率很低，在效率关键的代码中应该慎用objects。
- 如果需要设计类，应该尽可能采用普通的property，而避免灵活但是效率很低的dependent property。如非确实必要，避免重载subsref和subsasgn函数，因为这会全面接管对于object的接口调用，往往会带来非常巨大的开销（成千上万倍的减慢），甚至使得本来几乎不是问题的代码成为性能瓶颈。
注：(Functions that are frequently used or that can take more time to execute are often implemented as executable files. These functions are called built-ins.
Identifying Built-In Functions：
As with M-file functions, you can identify which functions are built-ins using the exist function. This function identifies built-ins by returning the number 5:
exist reshape
ans =
5
)
让MATLAB更快(2008-5)
MATLAB是一门在数值运算方面非常高效的语言，对于提高MATLAB的效率，有一些大家熟知的方法，比如避免使用for循环，向量化，使用C-mex写核心，等等。这些原则大体是没错的，但是教条地运用有时反而会降低效率。
根据问题的规模，适当选择向量化的策略
向量化的基本思想是以空间换时间，通过把要处理的对象转化成一个矩阵，集中处理。这比起非向量化的处理方式往往需要消耗更多的内存。在内存紧张的时候，内存分配的开销可能是很大的。
比如，给定一个列向量v1，和一个行向量v2, 计算矩阵M，使得M(i, j) = v1(i) + v2(j)。一种常用的向量化实现：
m = length(v1);
n = length(v2);
M = repmat(v1, [1, n]) + repmat(v2, [m, 1]);
这个过程中，除了M以外，还要创建两个大小为m x n的临时矩阵。repmat创建临时矩阵有一定的overhead，主要花费在元素索引上。由于元素加法的向量化收益本身并不显著，扣除repmat的额外开销后，得益就更小了。另外，当这些矩阵很大，接近物理内存容纳能力的时候，必然导致和虚拟内存进行频繁的交换，耗费大量的IO时间，向量化可能得不偿失。这种情况下，下面这种向量化程度略低的方式可能是更好的选择：
M = zeros(m, n);
if m < n
for i = 1 : m
M(i, :) = v1(i) + v2;
end
else
for i = 1 : n
M(:, i) = v1 + v2(i);
end
end
这个例子是要说明在设计向量化的时候，要对其它潜在的影响因素也进行考虑，权衡得失。对于这个例子本身，在2007版以后的matlab，有一个通用的简洁的写法。
M = bsxfun(@plus, v1, v2);
MATLAB内部对于bsxfun有很专业的优化，这个写法在时间和空间上都比较节省。
对一组d x d矩阵求逆，假设这些矩阵放在大小为d x d x n的数组里。在一般条件下，多个矩阵求逆没有特别的向量化方法，通常就是逐个求
n = size(A, 3);
B = zeros(size(A));
for i = 1 : n
B(:, :, i) = inv(A(:, :, i));
end
不过，当d很小而n很大的时候，比如对大量2x2矩阵求逆（这在设计几何的问题中很常见），那么直接利用2x2矩阵的求逆共识，进行向量化了：
a = reshape(A(1, 1, :), [1, n]);
b = reshape(A(1, 2, :), [1, n]);
c = reshape(A(2, 1, :), [1, n]);
d = reshape(A(2, 2, :), [1, n]);
g = 1 ./ (a .* d - b .* c);
r11 = d .* g;
r12 = -b .* g;
r21 = -c .* g;
r22 = a .* g;
B = reshape([r11; r21; r12; r22], [2, 2, n]);
在正式的代码中，还需要处理可能为奇异(a*d-b*c=0)的情况。
某些非数值计算问题的向量化
有些问题，表面上不像是数值计算，但是通过适当的转化，也可以通过向量化解决。
非等量复制: 比如有一组数[10, 20, 30], 各自要复制成[3, 2, 5]份，最终产生出 [10 10 10 20 20 30 30 30 30 30]。这个问题，虽然简单，但是因为不是那么“整齐”，因此向量化不是特别直接。下面给出一种思路：
首先产生出起始位置标志向量[1 0 0 1 0 1 0 0 0 0]，然后通过cumsum，产生出[1 1 1 2 2 3 3 3 3 3]，最后以此为索引，可以得到结果。
% vs: the values to be copied
% ns: the numbers of copies
vs = vs(ns > 0);
ns = ns(ns > 0);
sp = [1, cumsum(ns(1:end-1)) + 1];
result = vs(cumsum(sp));
对于图像处理或者二维信号问题，善用filter。
比如，要对每个像素，基于其局部邻域进行一些计算。即使整个计算本身不是线性的，但是只要能分解成线性局部运算的组合，就可以利用filter。
比如，对于图像I, 产生V，使得V(i, j)是像素I(i, j)周围w x w邻域的像素的方差。
I = im2double(I);
h = ones(w, w) / (w * w);
M = imfilter(I, h, 'symmetric');
M2 = imfilter(I.^2, h, 'symmetric');
V = M2 - M.^2;
