上一篇介绍了LDA的理论推导，这一篇进入LDA的实战，大致介绍下我用LDA做过的一些事情，欢迎大家一起探讨~

开源代码

　　我阅读修改过三个版本的LDA开源代码，按时间先后分别是GibbsLDA++/online_twitter_lda/ompi-lda。修改后的代码请移步到我的github https://github.com/nanjunxiao/LDA

1. GibbsLDA++

　　GibbsLDA++是单线程版LDA，由C++编写，代码很干净整洁，很适合入门学习。主要输出以下几个文件：

　　1). tassign文件（一行一个doc，冒号前是wordid,冒号后是topicid）
　　2). theta文件（doc->topic矩阵：行号=doc idx，列号=topic idx，值为文档属于该主题的概率）
　　3). phi文件（topic->word矩阵：行号=topic idx，列号=word idx，值为词在主题上的概率分布）

　　我修改了GibbsLDA++两个bug:

　　1). 内存泄漏。new的数组释放需要delete[]，而不是delete，这个问题比较多。

　　2). 数组访问越界。(double)random() / RAND_MAX产生值范围是[0,1]，int topic = (int)(((double)random() / RAND_MAX) * K); 会有生产K值的可能，这时就越界了。同理double u = ((double)random() / RAND_MAX) * p[K - 1]; u会有等于p[K-1]的可能，这时轮盘赌也越界了。修改办法就是将RAND_MAX改为RAND_MAX+1。

2. ompi-lda

　　ompi-lda是多机多线程版LDA，由C++编写，由王益大神编写。用到了mpi和OpenMP进行多机和多核实现。

　　我修复了infer.cc编写错误等不能编译问题，并且去除了lda.cc中’sampler.UpdateModel(corpus)’，我发现这个没有必要清空模型，采样需要不断累积。

　　我增加了theta和twords文件的输出，并且将依赖的boost库hpp/cpp文件打包到了include文件夹，这样不需要系统安装boost库，可以直接make。这个后面可以用标准C++重写cmd_flag，这样就不需要对boost的依赖了。

3. online_twitter_lda

　　online_twitter_lda是多线程版LDA，由python编写，可以按时间来增量计算。我也对其增加了theta和phi文件的输出。

参数设置

　　主题数K，我的经验一般取sqrt(#docs)/3，可以调研下HDP它不需要提前指定K。当然也可以用kmeans给定步长试探几次，选取sum(类内距离)/sum(类间距离)最小的K。

　　Dirichlet超参数α/β，经验值α=50.0/K，β=0.1

　　迭代次数iter_num一般>500次，当然越大效果越好。

并行实现

　　Gibbs Sampling LDA的并行化，存在着一个主要难题，就是nw、nwsum、nd、ndsum这几个统计量的同步更新问题，如果简单地将文章集拆分成若干份并行计算，A进程和B进程同时启动时，由于每次重新指定主题后，都会修改掉统计量（-1 重新分配topic +1），因此就会造成修改读写冲突，破坏统计量的一致性。

　　为了解决这个问题，David Newman等人提出了AD-LDA算法，该方法将文章切分到不同的机器上（按行拆分），该算法由于可以被转换成map-reduce操作，因此早期的LDA并行化实现都采用AD-LDA，该算法的核心是提出了global update操作，整个步骤如下所示
　　

　　这个算法被看作单机版本的Gibbs Sampling的近似，因为在各个机器上不同进程启动Gibbs Sampling时，nw_p和nwsum_p互不知晓其他机器的存在而进行采样，就会造成前面提到的修改冲突，破坏统计量一致性，所以说这一步其实已经产生了误差，不能完全等同于单机版的Gibbs sampling，但其后的global update在每一轮迭代后都将这个问题尽可能地修复了。

　　并行LDA还有其他一些优化实现方法，比如腾讯的peacock系统，老师木的lightLDA等，what a shame，我还没有深入了解。

我的应用

1. 文本聚类

　　对话题语料进行细粒度子话题聚类，并基于不同时间片子话题间关联关系构建话题的演化历程，是子话题发现与演化的一种思路。

　　算法基本思想是时间片t的子话题通过与其前向时间片t-1和后向时间片t+1的子话题计算相似度，如果相似度超过阈值，两个时间片子话题间就构建一条边关联，可以将子话题演变状态归结为新生、发展、分裂、合并、消亡5种。

　　子话题状态确定伪代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

1).	初始化：t时间片所有子话题的prior和post、fromset、toset都为空，状态为未知。
	t的子话题状态需要在t+1时刻最终确定，因为需要t+1时刻才能确定消亡状态

2).	如果t时间片子话题不为空，遍历所有子话题，计算对t-1时刻的prior；同时计算所有t-1时刻子话题对t时刻的post，超过阈值两个子话题之间建立一条边，分别插入到相应的fromset和toset集合

3).	确定t-1时刻子话题最终状态。
	遍历t-1时刻的所有子话题Tt-1i：
		If  len(Tt-1i.toset)==0：
			消亡状态；

4).	遍历t时刻的所有子话题Tti：
		If  len(Tti.fromset)==0：
				新生状态；
		Elif  len(Tti.fromset)>=2：
				合并状态；
		Elif  Tti ==post(prior(Tti)) and len(prior(Tti).toset)==1：
				发展状态；
		Else：
				分裂状态；

　　大致效果图如下
　　

2. 文档打标签

　　利用LDA的历史训练结果，做新文章自动打Tag服务同样可以得到不错的结果，该方法通过infer得到文档的theta分布，找到该文档概率最大主题编号，然后根据phi分布获取主题编号的topN词，作为文章标签输出。

　　该方法对新闻文本效果较好，对于微博这种短文本严重依赖分词好坏。

　　代码可参考我的Github https://github.com/nanjunxiao/tag_doc_with_lda

3. 文本分类

　　Kaggle的“Bag of Words Meets Bags of Popcorn”比赛中，我尝试使用LDA的主题分布作为特征进行情感分类，使用linear SVM AUC大概0.93左右。在加入BOW和NGram后，linear SVM可以提高到0.95998。

参考资料

　　马晨. LDA算法漫游指南

一. 几种概率分布

1. Gamma分布

　　
　　他有一个重要性质，)，对于整数可以得出阶乘形式

2. 二项分布(Binomial distribution)

　　二项分布是N重伯努利实验，描述发生k次事件的概率分布，可以理解为抛N次硬币k次正面的概率。
　　

3. beta分布(beta distribution)

　　Beta分布是指一组定义在区间(0,1)的连续概率分布，是对二项分布中p进行建模的概率分布，和二项分布互为共轭，这个后面会给出定义和缘由。有没有觉得二项分布/beta分布形式很像!
　　

　　Beta分布的期望
　　

4. 多项分布(multinomial distribution)

　　多项分布是二项分布的推广扩展，在n次独立试验中每次只输出k种结果中的一个，且每种结果都有一个确定的概率p。多项分布给出了在多种输出状态的情况下，关于成功次数的各种组合的概率。这里直观类比为扔骰子就好了。
　　

5. 狄利克雷分布(dirichlet distribution)

　　dirichlet分布是beta分布在多项情况下的推广，也是多项分布的共轭先验分布，同样的是对多项分布中P进行建模的分布。同样多项分布和dirichlet分布形式也很像。
　　

二. 共轭先验分布(conjugacy prior)

　　所谓共轭，指我们选取一个函数作为似然函数的prior probability distribution，使得后验分布函数(posterior distributions)和先验分布函数形式一致。比如Beta分布是二项式分布的共轭先验概率分布，而狄利克雷分布(Dirichlet分布）是多项式分布的共轭先验概率分布。为什么要这样做呢？归根结底还是为了推导简单漂亮！

　　根据贝叶斯公式，后验分布=似然函数*先验分布，如下面公式：
　　

　　下面以二项分布和Beta分布为例，证明他们互为共轭，多项分布和dirichlet类似。
　　

　　可以看到后验分布(post distribution)又变为beta分布了。这里超参数α,β称之为伪计数(pseudo count)，表示先验计数对似然计数进行修正，这样后验更精确。

三. 参数估计

　　参数估计一直存在两个学派：贝叶斯学派和频率学派。频率学派通过似然函数来选择特定参数值；贝叶斯学派认为参数不是固定值也是变量，给参数赋予先验分布，并使得先验与似然共轭，通过求后验均值来得到参数的估计。不管是哪个学派思想，都要用到似然函数，注意到频率学派所使用的似然函数是N次伯努利实验下的似然函数，但贝叶斯学派所使用的似然函数是二项式分布形式的似然函数

　　当数据量无限大时，贝叶斯方法和频率学派方法所得到的参数估计是一致的。当在有限的数据量下，贝叶斯学派的参数后验均值的大小介于先验均值和频率学派方法得到参数估计之间。

四. LDA推导

1. 概率图模型表示

　　隐含狄利克雷分布（Latent Dirichlet Allocation, LDA）是一个三层贝叶斯生成模型，它把PLSI模型的参数文档-主题/主题-词分布也看做随机变量，并采用Dirichlet分布作为先验。如下图所示，LDA的三层结构被三种颜色表示出来：

　　1). 语料层（红色）：α和β表示语料级别的超参数，每个文档都一样，被整个语料共享。
　　2). 文档层（橙色）：θ是文档级别的变量，每个文档对应一个θ，每个文档产生各个主题z的概率是不同的，生成每个文档需采样一次θ。
　　3). 词层（绿色）：z和w都是单词级别变量，z由θ生成，w由z和β共同生成，一个单词w对应一个主题z。
　　

　　Dirichlet分布是多项分布的共轭先验概率分布，LDA采用服从Dirichlet分布的K维隐含随机变量表示文档的主题混合比例，用一个服从Dirichlet分布的V维隐含随机变量表示主题中词典集的概率分布。在许多应用场合，我们使用对称Dirichlet分布，即各维α相同。其超参数是两个标量：维数K和参数向量各维均值α=(∑α_k )/K。

2. LDA生成文档过程

　　如LDA概率图模型所示，一篇文档生成的过程如下
　　

　　1)．从 Dirichlet分布中取样生成文档i的主题混合比例θi，α为Dirichlet分布的参数，也称为“hyperparamer”。
　　2)．从以θi为参数的主题multinomial分布中取样生成文档i第j个词的主题zi,j
　　3)．从Dirichlet分布中取样生成主题zi,j的词混合比例φzi,j，同样β也是Dirichlet分布的超参数。
　　4)．从以φzi,j为参数的词multinomial分布中采样词语，最终生成词语wi,j

　　整个模型中所有可见变量和隐含变量的联合概率分布为：
　　

　　最终生成一篇文档的概率可以通过将上式的θi及Φ进行积分并对zi进行求和得到。
　　

3. 贝叶斯推断-Gibbs采样

　　计算LDA中未知的隐含变量z¬m,n, θ¬m和φk的后验概率是概率图模型的推断（Inference）问题。主要的算法分为精确推断和近似推断两类。LDA用精确推断解起来很困难，所以一般采用近似推断来学习隐含变量。比如LDA原始论文[Blei 2003]中使用的mean-field variational EM算法和[Griffiths 2002]使用的Gibbs采样（Gibbs Sampling），其中Gibbs采样更为简单精确。

　　Gibbs采样（Gibbs Sampling）是一种Markov Chain Monte Carlo算法，通过生成一个观察样本序列来逼近多变量的联合概率分布，其基本思想是使用一个马尔可夫链条，sample出一系列的状态点，使其最终的平稳分布状态就是我们给定的那个联合概率分布（该联合概率就是LDA里的文档集被生成的概率）。该算法每次选取联合概率的一个变量维度，根据其他变量值采样当前维度变量值(完全条件概率，full conditionals)。不断迭代该过程，直到收敛。

　　整个文本训练集生成的联合概率为，有了联合概率分布，紧接着就可以推导full conditionals
　　
　　
　　
　　有没有觉得full conditionals公式很漂亮？而且含义很直观，其实就是p(topic|doc)*p(word|topic)，这个就是doc->topic->word产生式路径，amazing!

　　gamma函数的出现有两个作用：换掉积分；在分子分母同时出现gamma函数时利用gamma是阶乘的特性，而约去。

　　做完多次主题采样后，假设这是已经进入细致平稳状态，根据后验期望作为参数估计，就可以得到θmat(doc->topic)和φmat(topic->word)两个重要的矩阵，至此隐变量推断就完成了！
　　
　　

　　根据上面推导，Gibbs采样推断LDA隐变量分布的算法伪代码如下：
　　

　　1). 初始化：首先对整个语料所有词遍历一遍，为其随机分配一个主题，即zm,n=k~Mult(1/K),其中m表示第m篇文档，n表示文档中的第n个词，k表示主题，K表示主题的总数。然后将对应的m文档中k主题出现的词个数、m文档的词总数、k主题对应的词t个数、k主题对应的总词数（ n(k)m, nm, n(t)k, nk）分别加1。
　　2). 采样过程：遍历整个语料的所有词，假如当前文档m的词t对应主题为k，则n(k)m, nm, n(t)k, nk分别减1，即先拿出当前词t。然后根据计算出的Gibbs updating rule，采样出新的主题k’，在对应的n(k’)m, nm, n(t)k’和nk’上分别加1。P(zi|z-i,d,w)称为Gibbs updating rule，即排除当前词的主题，根据其他词的主题估计当前词分配各个主题的概率。
　　
　　3). 重复迭代2 采样过程。
　　4). 迭代终止后（主题分布和词分布收敛或到达最大迭代次数）输出文档-主题分布θ和主题-词分布φ：
　　

　　采样过程对应的C++代码如下：
　　

　　其中821-831行对应轮盘赌法，求得新采样的主题编号k’。注意812和837行其实是可以省略的，因为前后p[k]的相对比例不变，经过825归一化概率相同。825行有out of bound风险，应改为RAND_MAX+1，这个下一篇会详细说。

4. 复杂度分析

　　1). 时间复杂度：Collapsed Gibbs Sampling LDA每次迭代都要把语料所有词属于K个主题概率采样一遍，因此每轮迭代循环数是doc_num * per_word_num * K = corpus_word_count *K，所以最终的时间复杂度=O(iter_num * corpus_word_count * K)。

　　2). 空间复杂度：算法运行中需要保存几个数组，耗费空间最多的有两个nw和nd，其中nw占用内存空间最大=主题数K * vocabulary_size，所以算法的空间复杂度=O(K * V)。

五. 参考资料

　　1). David M.Blei & Andrew Y.Ng. Latent Dirichlet Allocation
　　2). Gregor Heinrich. Parameter estimation for text analysis
　　3). 靳志辉(Rickjin). LDA数学八卦 2013.2.8
　　4). 马晨. LDA算法漫游指南

　　今天整理资料发现之前收集的比较好的机器学习cheat sheet和tips，这里备忘分享下~
　　
　　
　　

　　GBM(gradient boosting machine)是一种ensemble机器学习框架，gradient boosting是boosting方法的一种，通过迭代拟合损失函数的负梯度值，我们常见的adaboost也是一种boosting方法，它通过迭代改变样本weight。

GBRT传统推导

　　GBM基本思想就是让损失函数持续下降，每次迭代先求得关于累加函数F(x)的负梯度(- gradient)，然后通过弱学习器f(x)去拟合该负梯度，然后将弱学习器f累加到F得到新的F。这里将函数F类比于参数theta就好理解了，平时给定模型我们如何迭代求解参数theta？对嘛，我们先求关于theta的目标函数负梯度，然后再和原来的theta累加更新为新的theta。这样两者就统一起来了，只不过GBM模型开始不是给定的，需要对f相加求得F，这也是为啥是ensemble的原因。

　　GBM可以处理分类/回归/排序问题，统一优化Object目标，区别仅在于损失函数不同而已，这个后面会推导。上面提到的弱学习器是回归模型，可以用各种回归模型，常见的是采用Tree based的Regression Tree，因此本文着重介绍GBRT(gradient boosting regression tree)。注:这里不用GBDT特指分类，统一使用GBRT，因为Object目标分类/回归没有本质区别，弱学习器统一做回归。

　　GBM算法流程图如下：
　　　　　　

　　首先计算负梯度方向，使用CART进行回归拟合gm，然后优化最优步长ρm，最后累加到F更新F。

GBRT xgboost版本推导

　　利用梯度能够很直观理解，也画龙点睛到了GBRT中的gradient。但陈天奇大牛给出了更一般的推导，对loss做了二阶泰勒展开，引入了二阶倒数海森矩阵，并加入了正则项整体求最优，更加精确通用的描述了GBRT的来龙去脉。(陈天奇PPT http://homes.cs.washington.edu/~tqchen/pdf/BoostedTree.pdf)

　　机器学习目标可以描述为Obj = Loss + regularization。

　　回归的Loss可以选择square loss=(y-F)^2，y和F取值实数；分类的loss可以选择hinge loss=max(0,1-yF)或者logistic loss=ln(1+exp(-yF) )，y取值+1/-1，F取值实数；排序的loss可以选择pointwize(转为二分类问题)的hinge loss或者logistic loss，pairwize(xi-xj作为样本新x，转为二分类问题)的hinge loss，listwize的NDCG loss。下面给出几种loss function的效果对比，这里假设真实y=1情况，横坐标m表示预测值，纵坐标表示loss。
　　　　　　　　　　　　　　

　　蓝色的是0-1loss，用于分类表示错误个数，往往作为和其他loss对比的baseline，红色的表示hinge loss，黄色表示logistic loss，绿色表示adaboost loss=exp(-yF)，黑色表示square loss。从上图可以看出:Hinge/logistic对于噪音函数不敏感，因为当m<0时，他们的反应不大，而square loss与adaboost loss可能更爱憎分明，尤其是square loss，因此对于分类问题square loss不太常用，更适合回归问题。

　　可见回归/分类/排序并未有本质上的区别，都是去最小化Obj，唯一不同的就是哪种loss function更适合而已。

　　因为现在的参数可以认为是在一个函数空间里面，我们不能采用传统的如SGD之类的算法来学习我们的模型，因此我们会采用一种叫做additive training的方式(boosting就是指additive training的意思)。每一次保留原来的模型不变，加入一个新的函数f到我们的模型中。
　　

　　现在还剩下一个问题，我们如何选择每一轮加入什么f呢？答案是非常直接的，选取一个f来使得我们的目标函数尽量最大地降低。

　　
　　
　　将目标Obj做二阶泰勒展开，除去常数项(包括l(y,y^(t-1) ) )，求得每个样本的一阶导g和二阶导h，将目标函数按叶子节点规约分组，得到下图。
　　

　　如果树结构是固定的时候，上式中Obj有闭式最小值解-叶子节点score Wj，如上图。

　　然而不幸的是，这时的树还是未知的，不过可以按照Gain最大化去构造。如果暴力的枚举所有CART树分裂情况，计算太复杂了，这里可以采用贪心算法：
　　

1
2
3

　　遍历所有特征：
	　　对每个特征所有取值排序
	　　线性扫描1遍，确定该特征最好分裂值

　　我们可以发现，当推导目标的时候，像计算叶子节点分数和split分支这样的策略会自然地出现，而不再是像我们从书上学到的利用启发式规则，比如基于gini进行分支，叶子节点采用平均值。一切都有理可循，make sense!

　　GBRT采用CART根据value值split成二叉树，因此适合数值特征，针对类别特征，需要进行one-hot-encoder编码。这也解释了我在kaggle实战(一)中提到的，GBRT对高维稀疏特征效果不好，以及对于年月日这种特征，不进行one-hot编码直接采用数值效果好的原因。

xgboost

　　陈天奇大牛开源了xgboost工具包(https://github.com/dmlc/xgboost)， 这是一个GBRT大规模并行开源框架。xgboost是各种比赛的利器，我参加的kaggle比赛基本都要用xgboost跑一组结果，同时也可应用到工业界。

　　推荐大家有时间阅读学习下代码，可以参考陈天奇的PPT和网上的”xgboost导读和实战”，个人建议画出代码的UML类图，不清楚的细节用gbd断点打印调试。
　　

总结

　　gradient版本将f类比于参数，通过f对负梯度进行回归，通过负梯度逐渐最小化Object目标；xgboost版本通过使得当前Object目标最小化，构造出回归树f，更直接。两者都是求得f对历史累积F进行修正。

　　对Obj进行二阶泰勒展开，是否可以像我上一篇无约束优化算法总结中写的，将传统的对负梯度的回归，改为对牛顿方向-H^-1*g的回归？GBRT是否有类似的NBRT名字？我感觉是可以的，至于为什么没有这么求解的还不清楚，求指导 :-)

　　机器学习的优化目标一般可以表示成Obj = Loss + regularization，有了这样一个优化目标剩下就是如何求解最小化问题了。机器学习绝大多数Obj都是非线性的，很难有形式解，这时需要各种近似优化算法来求Obj的极小化。如果Obj是凸函数，那么极小值等价于全局最小值，否则极小值未必是全局最优值。这里将Obj抽象为f(x)，并假设f为凸函数，且二阶连续可导。

梯度下降法

　　这个很常见也很简单，参数修正方向就是负梯度方向，步长可以定长也可以单独优化，即Xk+1 = Xk – lambda*grad.

牛顿法

　　牛顿法最初就是用来迭代求解方程根，想必大家高中时候都接触过。原理是利用泰勒公式在x0处一阶展开，即f(x)=f(x0)+f’(x0)(x-x0)。求解方程f(x)=0得，x=x1=x0-f(x0)/f’(x0)，由于泰勒展开只是近似等价，所以这里的x1也只是比x0更接近f(x)=0的近似解，所以需要不断的迭代来逼近真实x。牛顿法通过下图一目了然：
　　　　　　　　　　　　　　　　

　　牛顿法不是用来求根的吗，那怎么用来求最优化？直观的，求f(x)的最小值，必要条件是f’(x)=0，so，可以用牛顿法来求f’(x)=0的解！

　　下面还是来公式推导下，还是利用泰勒公式,二阶展开，f(x+deltax)=f(x)+f’(x)deltax+1/2f’’(x)deltax^2。当deltax无限趋近于0，上面式子近似为f’(x)deltax+f’’(x)deltax^2=0，即deltax = - f’(x)/f’’(x)。所以迭代公式为xk+1 = xk-f’(x)/f’’(x)。

　　牛顿法比梯度下降更容易收敛，速度更快，一般认为牛顿法利用了曲线更多的信息（二阶导数）。根据wiki上的解释，从几何上说，牛顿法就是用一个二次曲面去拟合你当前所处位置的局部曲面，而梯度下降法是用一个平面去拟合当前的局部曲面，通常情况下，二次曲面的拟合会比平面更好，所以牛顿法选择的下降路径会更符合真实的最优下降路径（from 知乎）。如下形象图，红色曲线是利用牛顿法迭代求解，绿色曲线是利用梯度下降法求解。
　　　　　　　　　　　　　　　　　　

　　扩展下对于多元x，梯度变为梯度向量，二阶导变为海森矩阵(Hessian matrix)，分别用字母g和H表示，迭代方向变为-H^-1g称为牛顿方向。

　　阻尼牛顿法相比牛顿法增加了步长lambda的精确搜索。

拟牛顿法

　　牛顿法存在两个主要缺点：

　　1. 需要计算海森矩阵和它的逆，计算复杂度高而且空间复杂度为O(N^2)
　　2. 有时海森矩阵无法保证正定，算法失效。

　　针对上面问题，研究者提出拟牛顿法(Quasi-Newton methond)，基本思想是不直接求海森阵或其逆，而是通过梯度构造出海森矩阵(或其逆)的近似，而且是正定对称的。

　　符号上用B表示对海森矩阵H的近似，D表示对逆矩阵H^-1的近似，Sk=Xk+1-Xk，yk=gk+1-gk。

　　各种拟牛顿法都需要满足拟牛顿条件：yk=Bk+1Sk或者Sk=Dk+1yk。下面大致介绍下DFP/BFGS/L-BFGS拟牛顿算法。

1. DFP算法

　　在步骤3还需要计算gk+1，下同。

2. BFGS算法

　　这时从海森阵计算逆计算量也不小，通过Sherman-Morrison公式，可以直接得到B的逆的递推关系，将B的逆替换为D，这样我们得到另一种BFGS。

3.L-BFGS算法

　　上面的拟牛顿法都没有解决空间复杂度O(N^2)的问题，为了解决该问题，我们就不能存储海森或其逆矩阵了。想想时间换空间，我们就只能存储生成海森的s和y向量了，需要海森时计算得出。而且，向量s和y也不是所有都存，而是固定存最新的m组。每次计算D时，只利用最新的m组s和y，这样空间复杂度降到了O(mN).

　　若记ρk=1/yk^T*sk，Vk=I-ρkyksk^T，算法2.3的步骤6可以近似为
　　　　　　

　　由于V和ρ标号增加方向相反所以需要前后分别遍历，总计两次遍历才能计算D。

　　由BFGS算法流程可知，求Dk*gk获取搜索方向才是最终目的，因此论文(updating quasi-newton matrices with limited storage)设计出了一种Dk*gk快速算法：

　　注意倒数第二行beta应该是i下标，最后算出的rL即为Dk*gk的值。

后记

　　本文算法截图来自博客(http://blog.csdn.net/itplus/article/details/21897715) ,感谢原作者的分享。

　　上面讨论的都是f(x)可导的情况，针对不可导还有一些变种算法，比如OWL-QN，有时间再研究总结。

回顾了下2011年开发的元搜索自动抽取组件，当时的定位是全自动通用。算法适用于所有搜索引擎和新加的微博检索，或者说类似结构的网页，因为业务比较专一吧，效果还是不错的，准确率召回率均超过95%。大概总结下:-)

代码请移步到我的github https://github.com/nanjunxiao/Arise

算法介绍

如上图大致是一个检索结果页面的DOM树结构，算法大致流程如下：

1. 找到所有锚文本链接

2. 进行分组（Xpath相同的锚文本链接）

3. 求出每组的最小父块(最小父节点子树)，找出字符数超过body子树字符总数一半的最小父块，这时可能有多个最小父块，找子树最小的（即Xpath路径最长的），即html-…-div-…-div。这之后可能还有多个相同的最小父块，这时的最小父块就为数据区域块，如图阴影区域。需要记录下最小父块为数据区域块的链接组和其在链接组中的位置的对应关系。

4. 对最小父块是数据区域块的链接组求出平均锚文本链接长度（文本链接比），最长的即为记录锚文本链接，比如绿色文本链接比大，绿色的就为记录锚文本链接

5. 求记录锚文本链接相对数据区域块的最大父节点，即数据区域的第一层孩子节点，作为数据记录块。

6. 找到数据记录块后，此处利用第一个和最后一个传统数据记录进行扩展，可以识别出其中的新型数据记录。

   （想用标签名、class属性识别新型数据记录，很难对所有搜索引擎都适用，从而做成通用的）

7. 数据记录块生成子树，抽取所需元信息

这里强调下3/4的顺序，由于右侧广告推广的文本链接比可能也很大，通过3可以去掉右侧广告噪音的影响。

性能指标

1. 准确率/召回率都超过95%

2. 平均处理速度10ms/page

后记

由于当时的定位是全自动通用,所以本算法采用了启发式规则。

如果现在重做，我可能会直接使用类似scrapy的xpath配置文件，毕竟短平快；或者复杂点利用记录的父/子/左右兄弟节点位置/标签/字数等属性构造特征，直接二分类来搞。

上一篇都是针对小数据集的，入门不建议从大数据集开始，可以不用考虑机器内存，不用out-of-core的online learning，不用考虑分布式，可以专注模型本身。

接下来我做了两个广告CTR预估相关的比赛，不过比赛当时都已经closed了，还好，我们还可以提交结果看看close时能排到的位置。

比赛实战

6. Display Advertising Challenge

Predict click-through rates on display ads. https://www.kaggle.com/c/criteo-display-ad-challenge

这是一个广告CTR预估的比赛，由知名广告公司Criteo赞助举办。数据包括4千万训练样本，500万测试样本，特征包括13个数值特征，26个类别特征，评价指标为logloss。

CTR工业界做法一般都是LR，只是特征会各种组合/transform，可以到上亿维。这里我也首选LR，特征缺失值我用的众数，对于26个类别特征采用one-hot编码，数值特征我用pandas画出来发现不符合正态分布，有很大偏移，就没有scale到[0,1]，采用的是根据五分位点（min,25%,中位数,75%,max）切分为6个区间(负值/过大值分别分到了1和6区间作为异常值处理)，然后一并one-hot编码，最终特征100万左右，训练文件20+G。

强调下可能遇到的坑：1.one-hot最好自己实现，除非你机器内存足够大(需全load到numpy，而且非sparse);2.LR最好用SGD或者mini-batch，而且out-of-core模式(http://scikit-learn.org/stable/auto_examples/applications/plot_out_of_core_classification.html#example-applications-plot-out-of-core-classification-py), 除非还是你的内存足够大;3.Think twice before code.由于数据量大，中间出错重跑的话时间成品比较高。

我发现sklearn的LR和liblinear的LR有着截然不同的表现，sklearn的L2正则化结果好于L1，liblinear的L1好于L2，我理解是他们优化方法不同导致的。最终结果liblinear的LR的L1最优，logloss=0.46601，LB为227th/718，这也正符合lasso产生sparse的直觉。

我也单独尝试了xgboost，logloss=0.46946，可能还是和GBRT对高维度sparse特征效果不好有关。Facebook有一篇论文把GBRT输出作为transformed feature喂给下游的线性分类器，取得了不错的效果，可以参考下。（Practical Lessons from Predicting Clicks on Ads at Facebook）

我只是简单试验了LR作为baseline，后面其实还有很多搞法，可以参考forum获胜者给出的solution，比如：1. Vowpal Wabbit工具不用区分类别和数值特征；2.libFFM工具做特征交叉组合；3.feature hash trick；4.每个特征的评价点击率作为新特征加入；5.多模型ensemble等。

7. Avito Context Ad Clicks

Predict if context ads will earn a user’s click. https://www.kaggle.com/c/avito-context-ad-clicks

跟上一个CTR比赛不同的是，这个数据没有脱敏，特征有明确含义，userinfo/adinfo/searchinfo等特征需要和searchstream文件 join起来构成完整的训练/测试样本。数据包含392356948条训练样本，15961515条测试样本，特征基本都是id类别特征和query/title等raw text特征。评价指标还是logloss。

由于数据量太大，跑一组结果太过耗时，根据比赛6的参考，目前我只选择liblinear lasso LR做了一组结果。最终目标是预测contextual ad，为了减小数据量，*searchstream都过滤了非contextual的，visitstream和phonerequeststream及params目前我都没有使用，但其实都是很有价值的特征（比如query和title各种similarity），后面可以尝试。

对于这种大数据，在小内存机器上sklearn和pandas处理起来已经非常吃力了，这时就需要自己定制实现left join和one-hot-encoder了，采用按行这种out-of-core方式，不过真心是慢啊。类似比赛6，price数值特征还是三分位映射成了类别特征和其他类别特征一起one-hot，最终特征大概600万左右，当然要用sparse矩阵存储了，train文件大小40G。

Libliear貌似不支持mini-batch,为了省事没办法只好找一台大内存服务器专门跑lasso LR了。由于上面过滤了不少有价值信息，也没有类似libFM或libFFM做特征交叉组合，效果不好，logloss只有0.05028，LB排名248th/414。

对于该比赛需要好好调研下大牛们的做法，看看相关paper了，自己瞎搞跑一遍太耗时间了，加油吧~

总结与感悟

通过参加kaggle提高了自己的机器学习实战能力，对问题和数据有了一些感觉，大致了解了各模型的适用场景。当然还有很多需要提高，比如特征组合/transform/hash trick，模型ensemble方法等，实现的scalable(比如采用pipeline)。

Ps:一定要挑选几个适合自己的高效工具包，并对其中2-3个看过源码，最好能做到定制优化。希望大家都加入到kaggle，欢迎一起探讨提高~

第一次接触kaggle比赛，是在听完台大林轩田老师的机器学习基石和技法课程之后。都说实践出真知，为了系统的巩固下机器学习实战技能，成为一名合格的数据挖掘工程师，我踏入了kaggle大门。

参加比赛很耗费时间和精力，由于本人已经工作，只能利用业余时间有选择的参加，希望能从中学到东西。我选择的比赛都是有监督的学习（当然用到了非监督方法，比如Bag of Words Meets Bags of Popcorn就用到了Collapsed Gibbs LDA），概括起来就是分类/排序/回归。下面就跟大家分享下近两个月我参加的比赛，欢迎一起探讨。

工欲善其事，必先利其器

首先介绍下我经常使用的机器学习工具:

1. scikit-learn. 涵盖了基本能想到的各种机器学习算法，由于本人python党，我把它当作matlab和R的替代品。
2. xgboost. 华盛顿大学机器学习大牛陈天奇出品的GBRT框架，果然是刷比赛利器。
3. liblinear/libsvm. 台大林智仁团队的佳作，工业界很多也在用。
4. pandas. 处理数据python包，DataFrame那叫一个好用。

机器学习流程

拿到一个比赛，我的一般套路是：

1. 读懂比赛介绍，明确是哪类问题:分类/排序/回归。
2. 数据特征处理。这个是最耗时也是最重要的，正所谓“数据和特征决定了效果上限，模型和算法决定了逼近这个上限的程度”。其实这点我还有很大欠缺，汗！
3. Cross validation数据集切分。数据集很大完全可以hold out一份作为测试集（不是待提交结果的测试集，此处是用来CV的），数据集偏小就需要K-fold或者Leave-one-out了，如果训练集有时序关系，还要注意测试集选取最后时间片的。这点我自我批评，有时为了省事，直接就提交结果做CV了。咳咳，这有点像imagenet比赛作弊了，只是我没用小号增加提交次数。
4. 常用算法/默认参数跑结果作为baseline。这个需要一些经验和直觉，一般来说Tree Based的Random Forest和GBRT效果都不会太烂，如果特征维度很大很稀疏这时就需要试试线性SVM和LR了。
5. 接下来就是调参了，这个我也没用太多经验，一般就是GridSearchCV或者RandomizedSearchCV。有人推荐Hyperopt库，接下来调研下。
6. 迭代。为了取得比较好的结果，下面就是2/3/4/5不断迭代了。
7. Blending.上面说的都是单模型，最后让你结果更general/low variance，提升一个档次的就是结果ensemble了（不是指gbrt/rf的ensemble，是多种模型的融合）。这里我一般就是简单的多种模型结果的averaging（weighted）or voting，这里推荐一篇ensemble selection paper(http://www.cs.cornell.edu/~alexn/papers/shotgun.icml04.revised.rev2.pdf)。

比赛实战

1. Bike Sharing Demand

Forecast use of a city bikeshare system. https://www.kaggle.com/c/bike-sharing-demand

这是一个回归问题，最后预测租车数量。这里需要注意一点，最后总数实际等于casual+registered。原始共10个特征，包括datetime特征，season/holiday等类别特征，temp/atemp等数值特征，没有特征缺失值。评价指标为RMSLE，其实就是RMSE原来的p和a加1取ln。

当时正在研究GBRT，所以使用了xgboost。由于使用RMSLE，xgboost自带的loss是square loss，eval_metric是RMSE，这时两种选择1.修改xgboost代码，派生新的优化objective，求新objective的gradient（一阶导）/hessian（二阶导），派生新的eval_metric；2.训练数据的y做ln(y+1)转化，最后预测时再做exp(y^)-1就转回来了。当然2简单了，我也是这么实施的。

关于数据特征处理，datetime转成y/m/d/h/dayofweek，y/m等类别特征由于有连续性，这里没有做one-hot编码。经过cv最后cut掉了日/season。

Xgboost参数其实没有怎么去调，shrinkage=0.1，tree_num=1000，depth=6，其他默认。

效果三次提升拐点分别是：1.RMSE转换为RMLSE(square loss转为square log loss)，说明预测值的范围很大，log转化后bound更tight了；2.cut了日/season特征；3.转换为对casual和registered的分别回归问题，在加和。最后RMLSE结果为0.36512，public LB最好为30位，最终private LB为28，还好说明没有overfit。

2. Bag of Words Meets Bags of Popcorn

Use Google’s Word2Vec for movie reviews. https://www.kaggle.com/c/word2vec-nlp-tutorial

这是一个文本情感二分类问题。25000的labeled训练样本，只有一个raw text 特征”review“。评价指标为AUC，所以这里提交结果需要用概率，我开始就掉坑里了，结果一直上不来。

比赛里有教程如何使用word2vec进行二分类，可以作为入门学习材料。我没有使用word embeddinng，直接采用BOW及ngram作为特征训练，效果还凑合，后面其实可以融合embedding特征试试。对于raw text我采用TfidfVectorizer(stop_words=’english’, ngram_range=(1,3), sublinear_tf=True, min_df=2)，并采用卡方检验进行特征选择，经过CV，最终确定特征数为200000。

单模型我选取了GBRT/NB/LR/linear SVC。GBRT一般对于维度较大比较稀疏效果不是很好，但对于该数据表现不是很差。NB采用MultinomialNB效果也没有想象的那么惊艳。几个模型按效果排序为linear SVC(0.95601)>LR(0.94823)>GBRT(0.94173)>NB(0.93693)，看来线性SVM在文本上还是很强悍的。

后续我又采用LDA生成主题特征，本来抱着很大期望，现实还是那么骨感，采用上述单模型AUC最好也只有0.93024。既然单独使用主题特征没有提高，那和BOW融合呢？果然work了!后面试验证实特征融合还是linear SVC效果最好，LDA主题定为500，而且不去除停用词效果更好，AUC为0.95998。

既然没有时间搞单模型了，还有最后一招，多模型融合。这里有一个原则就是模型尽量多样，不一定要求指标最好。最终我选取5组不是很差的多模型结果进行average stacking，AUC为0.96115，63位。最终private LB跌倒了71st，应该融合word enbedding试试，没时间细搞了。

3. Titanic: Machine Learning from Disaster

Predict survival on the Titanic. https://www.kaggle.com/c/titanic

硬二分类问题，不需要预测概率，给出0/1即可，评价指标为accuracy。说句题外话，网上貌似有遇难者名单，LB上好几个score 1.0的。有坊间说，score超过90%就怀疑作弊了，不知真假，不过top300绝大多数都集中在0.808-0.818。这个题目我后面没有太多的改进想法了，求指导啊~

数据包括数值和类别特征，并存在缺失值。类别特征这里我做了one-hot-encode，缺失值是采用均值/中位数/众数需要根据数据来定，我的做法是根据pandas打印出列数据分布来定。

模型我采用了DT/RF/GBDT/SVC，由于xgboost输出是概率，需要指定阈值确定0/1，可能我指定不恰当，效果不好0.78847。效果最好的是RF，0.81340。这里经过筛选我使用的特征包括’Pclass’,’Gender’, ‘Cabin’,’Ticket’,’Embarked’,’Title’进行onehot编码，’Age’,’SibSp’,’Parch’,’Fare’,’class_age’,’Family’ 归一化。我也尝试进行构建一些新特征和特征组合，比如title分割为Mr/Mrs/Miss/Master四类或者split提取第一个词，添加fare_per_person等，pipeline中也加入feature selection，但是效果都没有提高，求指导~

4. San Francisco Crime Classification

Predict the category of crimes that occurred in the city by the bay. https://www.kaggle.com/c/sf-crime

这是一个多分类问题，一般三种处理方法：one vs all, one vs one, softmax，信息损失逐渐递减。87809条训练数据，数据包括datetime/类别/数值特征，没有缺失值，label共39种。评价指标为logloss，这里要说下和AUC的区别，AUC更强调相对排序。

我抽取后特征包括year,m,d,h,m,dow,district,address,x,y，模型选择softmax objective的LR和xgboost。这两个模型对特征挑食，有不同的偏好，LR喜好0/1类别或者locale到0-1的数值特征，而xgboost更喜好原始的数值特征，而且对缺失值也能很好的处理。所以对于LR就是2个归一化的数值特征和8个待one-hot编码的特征，对于xgboost是8个原始数值特征（包括year/m/d等，具有连续性）和2个待one-hot编码的特征。

LR效果要略好于xgboost效果，logloss分别为2.28728/2.28869，最好位置为3rd，目前跌到4th，后面找时间再搞一搞。

5. Caterpillar Tube Pricing

Model quoted prices for industrial tube assemblies. https://www.kaggle.com/c/caterpillar-tube-pricing

这也是一个回归问题，预测tube报价。30213条训练样本，特征分散在N个文件中，需要你left join起来。评价指标为RMLSE，哈哈，是不是很熟悉？对，跟bike sharing的一样，所以怎么转换知道了吧？看，做多了套路trick也就了然了。感觉这个需要领域知识，但其实有些特征我是不知道含义的，anyway，先merge所有特征不加domain特征直接搞起。

这是我见过小样本里特征处理最麻烦的(后面的CTR大数据处理更耗时)，它特征分散在多个文件中，还好我们有神器pandas，直接left join搞定。这里有些trick需要注意，比如comp_*文件要用append不能join，这样正好是一个全集，否则就会多个weight特征了。特征存在缺失值，这里我全部采用0值，不知是否恰当？

模型我主要试了RF和xgboost，RF tree_num=1000，其他默认值，RMLSE=0.255201，主要精力放在了xgboost上，调了几次参数(depth=58,col_sample=0.75,sample=0.85,shrinkage=0.01,tree_num=2000)，最好RMLSE=0.231220，最好位置120th，目前跌倒206th了，看来需要好好搞搞特征了！

Join Kaggle

你看我搞的还凑合吧，so，快来加入kaggle吧，求组队:-)

Welcome to Hexo! This is your very first post. Check documentation for more info. If you get any problems when using Hexo, you can find the answer in troubleshooting or you can ask me on GitHub.

Quick Start

Create a new post

1

$ hexo new "My New Post"

More info: Writing

Run server

1

$ hexo server

More info: Server

Generate static files

1

$ hexo generate

More info: Generating

Deploy to remote sites

1

$ hexo deploy

More info: Deployment