近日来，抽空再读了一遍点集拓扑(Point Set Topology)，这是我第三次重新学习这个理论了。我看电视剧和小说，极少能有兴致看第二遍，但是，对于数学，每看一次都有新的启发和收获。

代数，分析，和拓扑，被称为是现代数学的三大柱石。最初读拓扑，是在两三年前，由于学习流形理论的需要。可是，随着知识的积累，发现它是很多理论的根基。可以说，没有拓扑，就没有现代意义的分析与几何。我们在各种数学分支中接触到的最基本的概念，比如，极限，连续，距离（度量），边界，路径，在现代数学中，都源于拓扑。

拓扑学是一门非常奇妙的学科，它把最直观的现象和最抽象的概念联系在一起了。拓扑描述的是普遍使用的概念（比如开集，闭集，连续），我们对这些概念习以为常，理所当然地使用着，可是，真要定义它，则需要对它们本质的最深刻的洞察。数学家们经过长时间的努力，得到了这些概念的现代定义。这里面很多第一眼看上去，会感觉惊奇——怎么会定义成这个样子。

首先是开集。在学习初等数学时，我们都学习开区间 (a, b)。可是，这只是在一条线上的，怎么推广到二维空间，或者更高维空间，或者别的形体上呢？最直观的想法，就是“一个不包含边界的集合”。可是，问题来了，给一个集合，何谓“边界”？在拓扑学里面，开集(Open Set)是最根本的概念，它是定义在集合运算的基础上的。它要求开集符合这样的条件：开集的任意并集和有限交集仍为开集。

我最初的时候，对于这样的定义方式，确实百思不解。不过，读下去，看了和做了很多证明后，发现，这样的定义一个很重要的意义在于：它保证了开集中每个点都有一个邻域包含在这个集合内——所有点都和外界（补集）保持距离。这样的理解应该比使用集合运算的定义有更明晰的几何意义。但是，直观的东西不容易直接形成严谨的定义，使用集合运算则更为严格。而集合运算定义中，任意并集的封闭性是对这个几何特点的内在保证。

另外一个例子就是“连续函数”(Continuous Function)。在学微积分时，一个耳熟能详的定义是“对任意的epsilon > 0，存在delta > 0，使得 。。。。”，背后最直观的意思就是“足够近的点保证映射到任意小的范围内”。可是，epsilon, delta都依赖于实空间，不在实空间的映射又怎么办呢？拓扑的定义是“如果一个映射的值域中任何开集的原像都是开集，那么它连续。”这里就没有epsilon什么事了。

这里的关键在于，在拓扑学中，开集的最重要意义就是要传递“邻域”的意思——开集本身就是所含点的邻域。这样连续定义成这样就顺理成章了。稍微把说法调节一下，上面的定义就变成了“对于f(x)的任意领域U，都有x的一个邻域V，使得V里面的点都映射到U中。”

这里面，我们可以感受到为什么开集在拓扑学中有根本性的意义。既然开集传达“邻域”的意思，那么，它最重要的作用就是要表达哪些点靠得比较近。给出一个拓扑结构，就是要指出哪些是开集，从而指出哪些点靠得比较近，这样就形成了一个聚集结构——这就是拓扑。

可是这也可以通过距离来描述，为什么要用开集呢，反而不直观了。某种意义上说，拓扑是“定性”的，距离度量是“定量”的。随着连续变形，距离会不断变化，但是靠近的点还是靠近，因此本身固有的拓扑特性不会改变。拓扑学研究的就是这种本质特性——连续变化中的不变性。

在拓扑的基本概念中，最令人费解的，莫过于“紧性”(Compactness)。它描述一个空间或者一个集合“紧不紧”。正式的定义是“如果一个集合的任意开覆盖都有有限子覆盖，那么它是紧的”。乍一看，实在有点莫名其妙。它究竟想描述一个什么东西呢？和“紧”这个形容词又怎么扯上关系呢？

一个直观一点的理解，几个集合是“紧”的，就是说，无限个点撒进去，不可能充分散开。无论邻域多么小，必然有一些邻域里面有无限个点。上面关于compactness的这个定义的玄机就在有限和无限的转换中。一个紧的集合，被无限多的小邻域覆盖着，但是，总能找到其中的有限个就能盖全。那么，后果是什么呢？无限个点撒进去，总有一个邻域包着无数个点。邻域们再怎么小都是这样——这就保证了无限序列中存在极限点。

Compact这个概念虽然有点不那么直观，可是在分析中有着无比重要的作用。因为它关系到极限的存在性——这是数学分析的基础。了解泛函分析的朋友都知道，序列是否收敛，很多时候就看它了。微积分中，一个重要的定理——有界数列必然包含收敛子列，就是根源于此。

在学习拓扑，或者其它现代数学理论之前，我们的数学一直都在有限维欧氏空间之中，那是一个完美的世界，具有一切良好的属性，Hausdorff, Locally compact, Simply connected，Completed，还有一套线性代数结构，还有良好定义的度量，范数，与内积。可是，随着研究的加深，终究还是要走出这个圈子。这个时候，本来理所当然的东西，变得不那么必然了。

• 两个点必然能分开？你要证明空间是Hausdorff的。
• 有界数列必然存在极限点？这只在locally compact的空间如此。
• 一个连续体内任意两点必然有路径连接？这可未必。

一切看上去有悖常理，而又确实存在。从线性代数到一般的群，从有限维到无限维，从度量空间到拓扑空间，整个认识都需要重新清理。而且，这些绝非仅是数学家的概念游戏，因为我们的世界不是有限维向量能充分表达的。当我们研究一些不是向量能表达的东西的时候，度量，代数，以及分析的概念，都要重新建立，而起点就在拓扑。

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