格致轩

些许记得杨绛写钱钟书的回忆录中对Bodelian的记忆，一下子我就来到了Bodelian。即使你不进Bodelian，只是在外围如那么多不知道从哪里来要到哪里去的游人一样拍照合影，大概你也能感受到它的三分气势。但是走马观花囫囵吞枣也许真的是不够，当我真正背着我的小黑thinkpad，进去查资料，看书，写程序的时候，我才明白它的好，即使我是在当今这么一个资讯发达的年代下研究热门的计算机视觉。

可乐音(kyle zheng)按：终于可以有那么一两天不用每时每刻关心学术上同行的一举一动，稍微有口喘息的机会来思考一下趣味。

看到人家实验室里面摆着各种有趣的盆景、金鱼还有摇摇脑袋的太阳花，真是有趣味。回来看看我桌上的摆设。什么都没有。除了各种专业书籍，仅有的闲书还是什么决策与判断、Vision之类的哲学书。各种个样的物品都是有他们的功利的目的。我真是一个无趣的人。

朋友要过来看看我之前的照片，翻开看看都是各种学术各种表情严肃僵硬，仿佛我就是一个机器人？一个天天跟电脑接触变得机械麻木的笨蛋。少数的几张去日本、韩国的照片还是表情呆呆的，一点都不能称得上可爱、帅气。完全就是一个没有生气的机器人。无趣的人怎么会有人爱呢？

冬冬说我眼睛总是盯着光芒，这样很不好。想起来，自己其实不是那么天生聪明，记得当年奥数模考100分也才考59，后来各种考试成绩也就能混到前10或者前5边缘线。也许只是因为自己胆子比较大，脸皮比较厚，尝试过一些，坚持了一下，机会就青睐了我一下，于是我装模作样混到现在这等样子。想起来，当年自己横下心把所有的业余爱好，比如足球、比如一个一个全部割舍，花比别人多一倍的时间在学业上，不知道到底值不值得？一个能够为了光芒和虚荣心，舍弃爱好，怎么会爱人？

看一下硅谷教父Paul graham。他自己本身从小就是喜欢编程的黑客。到后来学习什么哲学，再后来在哈佛学习人工智能。发现这些东西他都不爱。最终他在罗德岛找到了画画的快乐。这种经历让他的品味一下子就达到了一定高度。到后来，他折腾Viaweb以及YC看上去都跟他这种特殊的经历和牛逼的品味有关系。有品味的人做到最后肯定不会差到哪去，因为每项看似随机选择往往都会选到对的。对不对？

当然，品味不是等于趣味。各种味也还是有高下之分的。前一段日子受到某个投行高人的感染，对高雅艺术我开始呈现一种假装的喜欢。居然还跟同学装模作样去了趟国家大剧院听希拉里专场，花了不少的银子买了CD。结果发现CD中只有两首曲子听上去那么悦耳。。。本来不爱，假装爱，好痛苦啊。后来进一步跟那个高人交流，原来人家爸妈小时候就让他听。。。不过随着时间的推移，我每天都会把这些音乐折腾到我各种设备中，坚持在背单词，看Times的时候听，渐渐的觉得也还不错。所以，假装爱到最后肯定就真的爱咯，虽然开始有点痛苦。

女生喜欢跳个肚皮桑巴舞、画个画、插个花什么的，貌似选择得比较多。而作为从原始人就承担保护责任的男生，拥有一身肌肉和力量自然是很有必要的。虽然不用去跟猛兽搏斗，但是拿来追女生也许还有加分吧？可惜我从大三进入健身房到现在，还是没有能够坚持把自己弄成型男。中间跟论文实验做斗争那会，居然有了小肚子。因此，内心充满的悔恨和不甘。最近几次去健身房跑到55分钟的时候，已经湿透了的我想起那么多悔恨和不甘，在莫扎特和贝多芬的召唤下，也鼓起勇气继续蹦跑。不过装逼过分也是有代价的，拉伸的时候为了逞能，跟一个能拉80kg的人比赛，拉了40kg三十多下，结果到现在双臂都没完全恢复。。。做计算机视觉和机器学习的人怎么干不过一个做数学的哥们呢？也许大家时间不一样吧。

当然，有趣味的人玩什么都是那么有意思。比如林博士居然做了一个“电脑修的好，好人做到老”的电脑恢复U盘。还有各种个样有趣的代码分享。学习学习。这篇博文全是把自己的弱点一一抖落出来，希望看到的人别看不起我！~

待续

有句名言叫做“Rome is not built in a day！”其意思或许是指繁荣昌盛的国家不是那么一天搞定。最近看ICCV2009的录取论文单中有篇论文似乎想卖个关子，它的题目是“Building Rome in a Day”。

事实上，它是用计算机视觉的相关技术来做到这个目标的。即通过计算机视觉方法重建罗马的过程。主要思路就是利用的仅仅是互联网上你我拍摄的许许多多的图像来拼接那个罗马出来。这样看来，它也不是一天建成的，至少需要大家都做贡献。呵呵。

当然，要这样搞，有很多困难。首先这些图像数量虽然数以百万计，甚至更多，利用这些图像的好处在易于获取，不需要专门的设备，不好的地方是图像毫无组织，而且拍摄条件千差万别。要想拼接这些图像怕没有那么难吧。想起当年我和几个同学一起做Imagine Cup2007的时候，也有类似的想法，其中搞定大观园的图像合成就费了死劲，具体可以看这里，有相关的Demo：
With a panoramic view on history: The Knowledge 360 Panorma Systemm

看上去还不错吧？当年也是好歹混到了那个单项Imagine Cup2007世界前4的，当然，有了好idea，其实还要有较强的执行能力才行。人家华盛顿大学可以用集群计算机网络实现了城市级的三维重建，我们没那设备，所以只能玩个小的。
另外，那天看人家玩google map中一个照片自动拼接漫游系统，似乎也很好玩。这个方向似乎很不错啊。期待华人能尽快占据这个领域的最新位置。

附录：
该项目的链接如下：
http://grail.cs.washington.edu/projects/rome/
Abstract：
We present a system that can match and reconstruct 3D
scenes from extremely large collections of photographs such
as those found by searching for a given city (e.g., Rome) on
Internet photo sharing sites. Our system uses a collection
of novel parallel distributed matching and reconstruction
algorithms, designed to maximize parallelism at each stage
in the pipeline and minimize serialization bottlenecks. It is
designed to scale gracefully with both the size of the problem
and the amount of available computation. We have experimented
with a variety of alternative algorithms at each stage
of the pipeline and report on which ones work best in a
parallel computing environment. Our experimental results
demonstrate that it is now possible to reconstruct cities consisting
of 150K images in less than a day on a cluster with
500 compute cores.

现在SNS的产品烧得到处都是，最近发现web2.0的音乐产品也是这样。不过总得来说，创意各有千秋。从这里我感受到自己在学习研究的机器学习和数据挖掘强大的应用前景。这里就我简单的体验感受说说。有的产品使用时间稍微久远，不一定及时。
1亦歌
http://www.1g1g.com
个人感觉这个网站比较炫酷，总会实现这样的效果：一天不经意地点开，突然发现一首天籁之音。酷酷的界面，感觉十分时尚。在实验室网路比较好的时候，基本上很流畅。
问题在于，似乎没有聚类算法或者学习算法分析出users的喜好，只是简单随机罗列歌曲。
2songtaste
http://www.songtaste.com/
这个网站基本上走得是类似douban的歌曲推荐模式。基于假设：大家认为好的那就是好的。不过个人感觉很多事情包括对于音乐的口味，反倒不如douban上那些精英人士推荐的有趣符合我的胃口。
3Midomi
http://www.midomi.com/
这个网站走得是技术路线。不过感觉有些单一。比较有趣的地方在于，通过随意哼唱，能够找到你想要但是忘记名字的歌曲。其中原理比较简单，不过是比 较其中不同users的不同声音。有点当年我们ImagineCup2008第二轮中一个作品通过语音来实现区分关键字的项目创意。
4谷歌音乐
刚刚尝试了谷歌的音乐的测试版本。其中比较好玩的地方在于提供了泡泡聚类不同种类档次的歌曲。不过个人感觉google还是在走精英路线和平民路 线之间做了个trade-off，这个简单的创意把深层次的聚类的算法埋后面。我假象的问题在于：是不是有学习users的算法？如果没有的话，用久了可 能users的忠诚度可能也会下降吧。

附录：

这次算是google.cn的一个本地化的新举措吧 看了下google.com没有这个东东

http://www.google.cn/music/songscreener

 

互动做得很好的 good visualization

下面是做个产品的google工程师的博客介绍

地址：http://www.google.cn/music/songscreener

我有幸在这个项目中扮演了主力ActionScript开发者和部分交互设计师的角色。

这个项目的交互创意、UI设计主要来自于Eicostudio（2004 年我大三的时候曾经膜拜过工作室创始人之一Rokey大牛的作品，他们都是国内最专业的UI设计师）和Google中国的UX团队。Google中国的 UX团队给予了我极大的帮助，在同Google UX团队的合作中，我能非常清楚地感受到他们每一个人对细节的执着追求和对中文搜索事业的敬业精神，每一个颜色，每一个字体样式，每一个UI组件的使用模 式和位置，都经过反复的斟酌才最终敲定。

同时，项目质量控制上面也要感谢J、M、B、T等谷歌Music团队的工程师的帮助。老实说，我从来没有在一个ActionScript项目中被内存回收折磨得不得不重写一个几千行的类 - 感谢诸位工程师的支持和建议。

具体如何使用我觉得就没必要具体介绍了吧，根据右侧的几个筛选条件进行选择，可以找到你想要的一种特定类型的音乐，十分方便，大家使用起来有什么意见和建议也欢迎在我的Blog留言。

更新：

1. 我只是为Google公司提供挑歌所需的主要技术解决方案，对于计划政策并不了解，也没有任何发言权；
2. 已 知的一个小问题是，在某些配置较落后的机器上，FlashPlayer 10的CPU占用反而比FlashPlayer 9要稍差，但不至于超过80%，这个我经过严格测试了，如果您发现CPU占用过高，请查看是否开启了多个挑歌窗口；一般来说，搜索结果展现出来之 后，CPU占用都在15%以下；
3. 歌曲无法播放很大的可能是因为使用了OpenDNS，这一点我也很无奈，似乎没什么很好的解决办法
4. 歌曲必须在新窗口中播放并非技术原因，如果能够在“挑歌”界面内播放，肯定是最好的。可惜由于一些众所周知的原因，这一点目前还无法做到。

引言：
古语言：三个臭皮匠，顶个诸葛亮。很早就想过通过组合弱分类器来实现一个强分类器的问题。当学习到线性分类器以及著名的决策树之类分类器后，很容易联想到高中时代线性规划里面对不少区域的划分。要是能够找到一种方法，组合这些简单的直线方程，不就可以实现“山寨诸葛亮”的理想？

实际上，在Kearns和Valiant在1989年大作中指出了这种算法的可行性。而后，Freund在1990年以及他和Schapire在1994-1996年提出了boosting整个算法思路，似乎这种算法走到了实际应用的开端。然而直到AdaBoost被viola在其人脸识别系统中运用(2001Viola和Jones)，这种方法才彻底开始暴火。在NIPS会议tutorial中，也可以看到Schapire神采奕奕的样子了。

Adaboost实际具体解决了两个问题：
怎么处理训练样本？
在AdaBoost中，每个样本都被赋予一个权重。如果 某个样本没有被正确分类，它的权重就会被提高， 反之则降低。这样， AdaBoost方法将注意力更多 地放在“难分”的样本上。
怎么合并若分类器成为一个强分类器？
强分类器表示为若干弱分类器的线性加权和形式， 准确率越高的弱学习机权重越高。

算法介绍：
通过研究在Schapire的大作中提到了一个Toy Game的例子，这里给出了一个类似的Matlab代码。终于感到了这个算法强大。
首先是程序需要产生一些随机的样本数据，然后分别调用其他的matlab函数实现分类结果输出。代码如下：
-----------------
clear all
clc
tr_n=200;   %the population of the train set
te_n=200;   %the population of the test set
weak_learner_n=20;  %the population of the weak_learner
tr_set=[1,5;2,3;3,2;4,6;4,7;5,9;6,5;6,7;8,5;8,8];
te_set=[1,5;2,3;3,2;4,6;4,7;5,9;6,5;6,7;8,5;8,8];
tr_labels=[2,2,1,1,2,2,1,2,1,1]';
te_labels=[2,2,1,1,2,2,1,2,1,1]';

figure;
subplot(2,2,1);
hold on;axis square;
indices=tr_labels==1;
plot(tr_set(indices,1),tr_set(indices,2),'b*');
indices=~indices;
plot(tr_set(indices,1),tr_set(indices,2),'r*');
title('Training set');
subplot(2,2,2);
hold on;axis square;
indices=te_labels==1;
plot(te_set(indices,1),te_set(indices,2),'b*');
indices=~indices;
plot(te_set(indices,1),te_set(indices,2),'r*');
title('Training set');
% Training and testing error rates
tr_error=zeros(1,weak_learner_n);
te_error=zeros(1,weak_learner_n);
for i=1:weak_learner_n
    adaboost_model=adaboost_tr(@threshold_tr,@threshold_te,tr_set,tr_labels,i);
    [L_tr,hits_tr]=adaboost_te(adaboost_model,@threshold_te,tr_set,tr_labels);
    tr_error(i)=(tr_n-hits_tr)/tr_n;
    [L_te,hits_te]=adaboost_te(adaboost_model,@threshold_te,te_set,te_labels);
    te_error(i)=(te_n-hits_te)/te_n;
end

subplot(2,2,3);
plot(1:weak_learner_n,tr_error);
axis([1,weak_learner_n,0,1]);
title('Training Error');
xlabel('weak classifier number');
ylabel('error rate');
grid on;

subplot(2,2,4);axis square;
plot(1:weak_learner_n,te_error);
axis([1,weak_learner_n,0,1]);
title('Testing Error');
xlabel('weak classifier number');
ylabel('error rate');
grid on;
-----------------
这里需要另外分别撰写两个函数，其中一个为生成adaboost模型的训练函数，另外为测试测试样本的测试函数。代码如下：
------------
function adaboost_model=adaboost_tr(tr_func_handle,te_func_handle,train_set,labels,no_of_hypothesis)
% 训练函数

adaboost_model = struct('weights',zeros(1,no_of_hypothesis),...
                        'parameters',[]); %cell(1,no_of_hypothesis));

sample_n = size(train_set,1);
samples_weight = ones(sample_n,1)/sample_n;

for turn=1:no_of_hypothesis
    adaboost_model.parameters{turn} = tr_func_handle(train_set,samples_weight,labels);
    [L,hits,error_rate] = te_func_handle(adaboost_model.parameters{turn},...
                                         train_set,samples_weight,labels);
    if(error_rate==1)
        error_rate=1-eps;
    elseif(error_rate==0)
        error_rate=eps;
    end

    % The weight of the turn-th weak classifier
    adaboost_model.weights(turn) = log10((1-error_rate)/error_rate);
    C=likelihood2class(L);
    t_labeled=(C==labels);    % true labeled samples

    % Importance of the true classified samples is decreased for the next weak classifier
    samples_weight(t_labeled) = samples_weight(t_labeled)*...
                    ((error_rate)/(1-error_rate));                

    % Normalization
    samples_weight = samples_weight/sum(samples_weight);
end

% Normalization
adaboost_model.weights=adaboost_model.weights/sum(adaboost_model.weights);
-------------

function [L,hits]=adaboost_te(adaboost_model,te_func_handle,test_set,...
                              true_labels)
%测试函数
hypothesis_n=length(adaboost_model.weights);
sample_n=size(test_set,1);
class_n=length(unique(true_labels));
temp_L=zeros(sample_n,class_n,hypothesis_n);

for i=1:hypothesis_n
    [temp_L(:,:,i),hits,error_rate]=te_func_handle(adaboost_model.parameters{i},...
                                                    test_set,ones(sample_n,1),true_labels);
    temp_L(:,:,i)=temp_L(:,:,i)*adaboost_model.weights(i);
end
L=sum(temp_L,3);
hits=sum(likelihood2class(L)==true_labels);

-------
其中上面函数由于体积太大，另外还需要分别撰写两个阈值函数和一个隶属分配函数。
-------------
function model=threshold_tr(train_set,sample_weights,labels)
% 训练阈值函数

model=struct('min_error',[],'min_error_thr',[],'pos_neg',[],'dim',[]);

sample_n=size(train_set,1);
min_error=sum(sample_weights);
min_error_thr=0;
pos_neg='pos';

% for each dimension
for dim=1:size(train_set,2)
    sorted=sort(train_set(:,dim),1,'ascend');
    
    % for each interval in the specified dimension
    for i=1:(sample_n+1)
        if(i==1)
            thr=sorted(1)-0.5;
        elseif(i==sample_n+1)
            thr=sorted(sample_n)+0.5;
        else
            thr=(sorted(i-1)+sorted(i))/2;
        end
       
        ind1=train_set(:,dim)<thr;
        ind2=~ind1;
        tmp_err=sum(sample_weights((labels.*ind1)==2))+sum(sample_weights((labels.*ind2)==1));
       
        if(tmp_err<min_error)
            min_error=tmp_err;
            min_error_thr=thr;
            pos_neg='pos';
            model.dim=dim;
        end
       
        ind1=train_set(:,dim)<thr;
        ind2=~ind1;
        tmp_err=sum(sample_weights((labels.*ind1)==1))+sum(sample_weights((labels.*ind2)==2));
       
        if(tmp_err<min_error)
            min_error=tmp_err;
            min_error_thr=thr;
            pos_neg='neg';
            model.dim=dim;
        end
    end
end
model.min_error=min_error;
model.min_error_thr=min_error_thr;
model.pos_neg=pos_neg;
---------------------
function [L,hits,error_rate]=threshold_te(model,test_set,sample_weights,true_labels)
% 测试阈值函数
feat=test_set(:,model.dim);
if(strcmp(model.pos_neg,'pos'))
    ind=(feat>model.min_error_thr)+1;
else
    ind=(feat<model.min_error_thr)+1;
end

hits=sum(ind==true_labels);
error_rate=sum(sample_weights(ind~=true_labels));

L=zeros(length(feat),2);

L(ind==1,1)=1;
L(ind==2,2)=1;

--------------
function classes=likelihood2class(likelihoods)
% 隶属分配函数
[sample_n,class_n] = size(likelihoods);
maxs = (likelihoods==repmat(max(likelihoods,[],2),[1,class_n]));

classes=zeros(sample_n,1);
for i=1:sample_n
    classes(i) = find(maxs(i,:),1);
end
------------
算法其他资源：
在OpenCV.org.cn论坛上，还有这个算法关于人脸识别运用的代码讲解，这里广告一下。

总结：
Adaboost的思想主要是将一系列粗略的规则加权组合起来得到高度精确的规则。其中，
算法优点在于1）易于执行；2）分类精度较高。其缺点在于：1）容易受到噪声干扰，
这也是大部分算法的缺点；2）执行效果依赖于弱分类器的选择。


参考资料：
[1] Cuneyt Mertayak的Matlab代码。
[2] R. Schapire, A Short Introduction to BoostingIn Proceedings of the Sixteenth International Joint Conference on Artificial Intelligence, pp. 1401- 1405 1999.
[3] P. Viola, M. Jones, Robust Real-Time Face Detection. International Journal of Computer Vision 57(2): 137-154, 2004.

Matlab2008b号称Matlab7.7，广告说是对于向量化和非向量化的区分已经不是那么明显了。直接使用for循环而不用向量化似乎效率也是可以的。但是事实胜于雄辩。这里实验了测试了一下三种方法的效率，效果明显，看来差距确实不是那么大了，但是还存在差距。
机器配置：
Intel T2080 1.73GHz 双核
内存1.5G
操作系统Windows XP sp3
Matlab 2008b （Matlab 7.7.0）
测试程序：
% Purpose :
% Caculate the squares from 1 to 10 000 in the region of int
clear all
clc

maxcount=1;
% runing period
tic;
% using tic/toc to record the time cost.
for jj=1:maxcount
        clear square
        for ii=1:10000
             square(ii)=ii^2;
        end
end
ave1=(toc)/maxcount; %direct caculating way without pre defining array

maxcount=50;
tic;
for jj=1:maxcount
        clear square
        square=zeros(1,10000);
        for ii=1:10000
                square(ii)=ii^2;
        end
end
ave2=(toc)/maxcount; %direct caculating way with pre defining array

maxcount=100;
tic;
for jj=1:maxcount
        clear square
        ii=1:10000;
            square=ii.^2;
end
ave3=(toc)/maxcount; %vectorilized caculating way
%display the results
fprintf('Loop/uninitialized array= %9.5f\n',ave1);
fprintf('Loop/initialized arrray= %9.5f\n',ave2);
fprintf('vectorized= %9.5f\n',ave3);


实验结果：
每次测得的时间会略有差距，但是对比其他环境Matlab2007a，for循环的计算速率明显有了提高。
Matlab 2008b测试结果
Loop/uninitialized array=   0.14229
Loop/initialized arrray=   0.00019
vectorized=   0.00008

Matlab2007a测试结果
Loop/uninitialized array=   0.47000
Loop/initialized arrray=   0.00021
vectorized=   0.00010

有的时候总有那么点胡思乱想，可能是坐的时间足够久了，“屁股也能决定大脑”。于是，响应同志们的号召，借着year of ox或者是year of bull到来的机会，参加了饭局。发现原来长时间不说话，长时间看不说人话的书，开始进入所谓的geek路线了。 也就有那么几个比较差劲的笑话作为铺垫吧。
那个笑话原来是说： 在堪萨斯城有幸听到 80 多高龄的于牧师讲道。他说一年轻人要和一女孩处对象，不知如何和人接洽。于是牧师建议，说圣诞节快到了，你不如买一礼物送她吧。年轻人买了礼物送去，结果被轰走，原来他送了人家一百斤米。女孩说：你当我是饭桶？年轻人改口说：这米是送你哥的。女孩怒不可遏：骂我哥饭桶！你哥才饭桶！你给我把米背回去！！！年轻人只好把米背走，边走边骂自己饭桶。
其实说白了，还是那个“如果你只有一把锤子，你看见的所有将是钉子”。我最近老是中着魔咒般的日子，似乎也开始走上这个路线。不过也许走多了也就好了，就像当年刚刚待北京一年，回家发现开始不适应方言全部的环境一样，总想和人说普通话，回北京又想和人说方言。 末了，为了打发自己在看论文写码累了能有东西陪伴，顺便买了只冰河世纪的抱了只大榛子的松鼠为伴。

motivation
为啥呢？
Introduction
无他，看看是否体温曲线在24小时内有高有低。
discussion and results
按powerful sleep的讲法，白天要有日光照、要有体育锻炼，晚上要有黑屋。注意什么时候入睡，而不是什么时候起床。控制好每天睡觉睡下和起来的时间符合4个周期长度。当然这个周期人人不同。