Dahua 的个人资料笑对人生,傲立寰宇照片日志列表 工具 帮助
3月30日

MATLAB的思维风格

最近做course project的时候,我发现所写的程序中有几个地方特别能反映MATLAB与别的高级语言很不一样的风格。

例1:一批user去访问文件,比如u1访问了f11, f12, ..., u2访问了f21, f22, ...,每个用户访问的文件数目可能不等,要求输出一个表M, M(i, j)是用户 i 和 j 都访问过的文件的个数。

先用高级语言给出一个比较直接的解决吧,下面是python的代码。(其实C++/C#/Java的也类似,不过python的代码比较简洁)

首先,建立一个访问字典(或者叫map),记录每个文件都被哪些人访问过。

# input is a map named ufiles, 
# ufiles[u] is a set of files accessed by u

acc_dict = {}
for u in users:
    for f in ufiles[u]:
        if f in acc_dict:
            acc_dict[f].add(u)  # add u to an existing set
        else:
            acc_dict[f] = set([u])  # initialize the set

然后,基于这个字典进行计数。对每一个文件,给每一对访问过它的用户的相应计数加1

# initialize the table
M = {}
for u1 in users:
    for u2 in users:
        M[(u1, u2)] = 0

# do counting
for f in acc_dict:
    uset = acc_dict[f]
    for u1 in uset:
        for u2 in uset:
            M[(u1, u2)] += 1

这样就基本完成计数了。

MATLAB没有像python那样丰富的数据结构,比如set(集合)和dict(映射表,字典)。但是,在这个表面上需要依靠经典数据结构完成的问题,我们可以用MATLAB给出一个更加高效简洁的解决方法

给出代码之前,先做点数学分析。令F是一个矩阵,F(u, f) = 1代表u访问过这个文件,否则F(u, f) = 0。那么M其实就可以写成:

M(u1, u2) = sum_f F(u1, f) F(u2, f) --- 很容易看出,这就是u1, u2都访问过的文件的数目。
写成矩阵形式:
M = F * F'

下面的代码就很顺理成章了

% Input: ufiles is a cell array of indices of files accessed by each user
% ufiles(i) is a row vector of the indices accessed by user i

% construct F 
n = length(ufiles)
% to produce usrs like {[1 1 1 1 ...], [2 2 2 ...], ...}
usrs = arrayfun(@(i) i * ones(1, length(ufiles{i})), 1:n, ...
                'UniformOutput', false)

I = [usrs{:}]  % concatenate the indices into one long vector
J = [ufiles{:}]
F = sparse(I, J, 1);

% compute M 
M = F * F'

注意这里用sparse matrix, 那么运算复杂度会显著降低,并且只依赖于实际访问的数量。实际性能上,通过MATLAB的稀疏矩阵乘法比自己用高级语言实现更快。

熟悉kernel的朋友可能已经看出来,MATLAB的做法其实是用file access indicator vector作为每个user的feature,构造出来的M则是相应的Gram matrix。不过,如果没有这方面的思维习惯的话,未必能很快就把一个统计文件访问量的问题联想到矩阵乘法或者内积计算。

例2:给出一组数字,算出每个不同数字出现的次数。

python的实现

counts = {}
for v in values:
    if v in counts:
        counts[v] += 1
    else:
        counts[v] = 1
用高级语言实现,一般需要一个映射表来快速定位某个数字的计数器。

Matlab的实现

svalues = sort(values);
dp = find(diff(svalues));
sp = [1, dp+1];
ep = [dp, length(svalues)];
U = svalues(sp);   % U is the set of distinct values
C = ep - sp + 1;   % C is the corresponding counts of these values

一般情况下,两种实现的时间复杂度都是O(n log(n))。不过,它们体现了不太一样的思维方式。高级语言倾向于通过数据结构加速每个元素的处理,而MATLAB没有复杂的数据结构,它注重于以整体方式进行操作。

3月19日

MATLAB 2008

这几个星期一直忙着作业和考试,这里有些天没有更新了。刚过了mid-term,回到这里和朋友们说说话。

今年3月,Mathworks推出了MATLAB一个重要的新版本,MATLAB 2008a,也叫做MATLAB 7.6。在这个版本里,MATLAB解决了几个长期以来固有的弊端,而且加入了一些重要的能力。这次的更新是非常aggressive的,可能代表了MATLAB的一种历史性的转型。

  1. 完全实现面向对象编程。其实,在MATLAB的早期版本里面,也有class的概念,不过大家如果使用过的话,可能知道那是一种不太好的设计。功能不强,过程繁琐,而且,很多很tricky的地方,尤其是重载numel, subsref这类函数的时候。而新的设计抛开了历史包袱,现在写出来的类和在python里面写的长得差不多,舒服多了。这套设计,吸取(或者“抄袭”)了Python和C#的优点,除了支持封装(encapsulation),继承(inheritance),和多态(polymorphism)这些基本特性以外,还支持了一些新兴的特性,包括属性(property),事件(event),和静态方法(static method)。
  2. 支持Handle类型——用另外一种说法,就是支持函数调用传引用。以前matlab传递参数只有一种方法,copy on write。就是说,当你传一个东西进去,如果它要发生改变,那么,这个东西会整个copy一份,然后修改会在副本上生效。这使得实现动态数据结构变得非常困难。比如一个列表,如果每添加一个元素,都要拷贝整个列表一次,将是什么效果呢?因此,传统上matlab擅长于以矩阵为基础的算法,但是对于以经典动态数据结构为基础的算法,比如动态列表,哈希表,搜索树,图等,就力不从心了。这个新版本终于引入了对引用的支持,这将使MATLAB实现经典数据结构和算法变得前所未有的轻松。现在,数值和统计算法与经典算法越来越多地合流,很多应用都需要同时使用两方面的算法,MATLAB的这个变化正好适应了这种需求。
  3. 引入了名空间的管理。以前,MATLAB所有的函数都在同一个global的名空间下面。比如两个工具包里面出现了同名函数,解决起来很麻烦。比如现在有两个算法叫LDA,一个是Latent Dirichlet Allocation,一个是Linear Discriminant Analysis,在一个应用中需要同时用到两个算法,而写这两个算法的人各自把它们命名为lda.m,那么问题就出来了。一种naive的方法是改名字,不过会直接破坏掉那些toolbox里面对那个函数的依赖。而这个版本,它借鉴其它高级语言的经验,终于引入了namespace,给这个问题一个很好的解决。

从这些特点看来,MATLAB这个版本的重要改变,就是全面吸收其它高级语言的特性,从一个数值运算语言开始迈向一个以数值计算为强项的通用语言,以应对复杂或者更大规模应用的需要。

一直以来,由于matlab缺乏处理高级数据结构和建立复杂应用的能力,它有一个有力的竞争者numpy,这是python里面进行矩阵和数值运算的包,它建基于python这种著名的通用语言,并且提供matlab矩阵的部分能力。这次MATLAB的全面升级,对于numpy无疑提出了严峻的挑战。

除了程序设计结构方面的变化,MATLAB 2008在多个方面也有重要进步。

  1. 它的优化工具箱(optimization toolbox)首次引进了interior point algorithm。interior point algorithm在convex optimization中占有重要地位,并且性能优越。MATLAB optimization toolbox一直以为因为使用老式算法,性能太差,而饱受诟病。这次终于引入interior point,希望它的优化性能能得到显著改善。
  2. 重写核心JIT引擎(运行时编译,可以显著提高运行效率),并且采用了最新的BLAS/Lapack核心,运算速度会有相当程度的提高。另外,还大幅度提高了对sparse matrix的计算速度。
  3. 它的统计工具箱增强了对很多著名的统计算法的有力支持,比如HMM, GMM,还有NMF(Nonnegative Matrix Factorization),并且开始引入对蒙特卡罗采样的支持。

MATLAB 2008刚刚发行,现在学校内部还没开始支持。不过,这个版本确实非常值得期待。

3月7日

课业

这两周经历了这个学期第一个忙碌的高峰。

这里的课业本来就很繁重,而且,很多时候并不均衡。最近,两门课程“不约而同”地把整个学期30%以上的作业量集中到了两个星期——我不记得这两周除了完成各种接续不断的assignment,还做过什么。

学计算机的人都知道循环语句,可是,如果这些人做TA,并且在布置作业的时候使用for循环,那是什么一种情况呢?最近一次作业的一道大题,竟然使用三重循环来布置,实在令人佩服。题目大概是长成这个样子:

1. (a).  (a.1) xxxx, (a.2), xxxx, .... (a.n) xxxx, (a.(n+1)) change some condition, please repeat from (a.1) to (a.n) ...... (b) ..... (c) ..... (d) ........ (e) with a new setting, please repeat (a) to (d) ..... 2 ..... 3 ..... 4 ...... 5 in a new method, please repeat 1 to 4.

除了两门课的作业和实验,还需要完成Network的research proposal。这里的课程project都有一个特点,全部都是不命题的。学生需要自己完成,从survey, proposal, formulation, experiment, 到最后的report的全过程。Network这门课的老师特别强调,project必须探索这个领域的研究前沿,最后的report必须写成可以发表的paper的形式,中间有一次中期检查,最后还有一次答辩。他不只一次的告诉我们,每年都有相当数量的project paper成功发表到top conference上——这应该是他对我们的一种期待——这对于这门课的学生,尤其是不在这个领域的,确实是极具挑战性的。

CS的决策人确实是用心良苦——要求每个学生修的课必须覆盖3大领域(system, theory, AI),而在大部分的课程中,教授都要求学生完成一个从选题到最后报告的完整的课题研究过程。虽然非常辛苦,但是,学生确实可以获得在计算机科学的各大领域开展研究的经验。另外一方面,几年前CS把所属的全部实验室合并成CSAIL——常年有800人工作的超大型实验室,鼓励各个研究小组相互交流。强制性的跨领域修课政策加上大实验室制,使得CSAIL内部的交叉研究特别活跃。很多研究课题都是“杂交”的。每个老师管着不同领域的课题,每个学生参与不同小组的meeting,已经司空见惯。在某些具体的领域,CSAIL的研究未必是最好的,但是,这里确实是世界上交叉学科发展最为蓬勃的地方。