Ensembling, Lagrange
据我观察，名字带“拉”的人一般都很厉害，比如马拉多纳、希拉里、狄波拉、张娜拉、杜拉拉、陈法拉、尼古拉斯·赵四、地铁站的罗拉……
但跟我们今天的大神拉格朗日（Lagrange）一比就都被秒成了渣渣，因为每一个上过高数和物理的我们都被拉大神无情碾压摩擦过啊~
童鞋，请站起来回答拉格朗日中值定律是啥？
童鞋，请站起来回答拉格朗日乘子法是啥？
童鞋，请站起来回答拉格朗日内插公式是啥？
童鞋，请站起来回答拉格朗日点是啥？
童鞋，请站起来回答拉格朗日函数是啥？
童鞋，请站起来回答拉格朗日方程是啥？
……
童鞋，你肿么站不起来了童鞋？
请原谅我召唤出封印在你们心底多年的梦魇，快跟我一起踏碎时空大逃亡，继续闪回到小米数据挖掘大赛那几天。
上文中已经提到我们最终的模型不止一个，而是多个DNN和FM共七八个模型的融合。
所谓模型融合，洋气一点叫Ensembling，在 http://mlwave.com/kaggle-ensembling-guide/ 一文中有详细介绍。简单来说就是综合多个单模型的输出来给出最终的结果，一般会比每个单模型的效果都好，现在已经是各大比赛的常规武器。但对于初涉江湖的小伙伴们还没有太多经验，一开始只是简单的平均，后来尝试了各个模型的输出做为LR的特征，发现没啥卵用，比赛已近尾声就没再折腾更复杂的方法，最终使用线性加权平均，即：
$$
\hat y=w_1\hat y_1+w_2\hat y_2+…+w_M\hat y_M
$$
其中，$w_1+w_2+…+w_M=1$
这里有个问题就是权重系数怎么定？
一开始就是拍脑袋定，单模型效果好的权重就大一点，效果差的就小一点。后来小伙伴使用暴力网格搜索的方法，融合两三个模型还行，再多就已经慢到不可接受。
我看在眼里急在心头，我们是模武双修的种族啊，怎么能只用暴力解决呢？这种目标如此鲜明灿若煌煌皓月的好问题，当然要祭出流光华丽的大模型才相得益彰呢！
第一式：混沌初分模型现！
比如我们有M个单模型分类器，解决K分类问题，测试集包含N条样本，$\hat{y}_{nmk}$ 表示第m个单模型对第n条样本属于第k类的预测概率。
我们采用线性加权平均的方法，融合M个模型得到的预测概率
$$
\hat{y}_{nk}=\sum_{m=1}^{M-1}w_m\hat{y}_{nmk}+(1-\sum_{m=1}^{M-1}w_m)\hat{y}_{nMk}
$$
比赛评价标准为logloss，我们希望融合后的模型在测试集上的logloss尽量小，所以优化目标如下：
$$
\min_{w}f(w)=\min_{w}\{-\frac{1}{N}\sum_{n=1}^{N}\sum_{k=1}^{K}1\{y_n=k\}\ln\hat{y}_{nk}\}\\
=\min_{w}\{-\frac{1}{N}\sum_{n=1}^N\sum_{k=1}^{K}1\{y_n=k\}\ln(\sum_{m=1}^{M-1}w_m\hat{y}_{nmk}+(1-\sum_{m=1}^{M-1}w_m)\hat{y}_{nMk})\}
$$
其中$w=[w_1,..,w_{M-1}]^T$是需要求解的权重参数，我们知道权重之和要归一，所以不需要$w_M$这一变量，用$1-\sum_{m=1}^{M-1}w_m$代替即可。
问题已明确，下面来求解。
第二式：梯度杀！！
这一招早已驾轻就熟，直接求梯度，然后上sgd或adagrad。单条样本的损失函数为
$$
l_n(w)=-\sum_{k=1}^{K}1\{y_n=k\}\ln(\sum_{m=1}^{M-1}w_m\hat{y}_{nmk}+(1-\sum_{m=1}^{M-1}w_m)\hat{y}_{nMk})
$$
梯度为
$$
\frac{\partial l_n}{\partial w_m}=-\sum_{k=1}^{K}1\{y_n=k\}\frac{\hat{y}_{nmk}-\hat{y}_{nMk}}{\hat{y}_{nk}},\quad m=1,…,M-1
$$
adagrad版本的代码如下，注意其中变量名和上面的公式并不完全一致，类别编号是从0到K-1（我就是这么洒脱随性~）1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43import sys, math
classNum = int(sys.argv[1])
modelNum = int(sys.argv[2])
alfa = float(sys.argv[3])
beta = 1.0
print >> sys.stderr, "classNum=" + str(classNum)
print >> sys.stderr, "modelNum=" + str(modelNum)
print >> sys.stderr, "alfa=" + str(alfa)
w = [1.0] * modelNum
n = [0.0] * (modelNum-1)
for i in range(len(w)):
w[i] /= modelNum
print >> sys.stderr, w
for line in sys.stdin:
seg = line.rstrip().split(" ")
label = int(seg[0])
x = []
for tmp in seg[1:]:
seg2 = tmp.split(",")
for i in range(len(seg2)):
seg2[i] = float(seg2[i])
x.append(seg2)
pre = [0.0] * classNum
for m in range(modelNum):
for k in range(classNum):
pre[k] += w[m]*x[m][k]
g = [0.0] * (modelNum-1)
wsum = 0.0
for m in range(modelNum-1):
g[m] = (x[modelNum-1][label]-x[m][label])/pre[label]
n[m] += g[m] * g[m]
lr = alfa/(beta+math.sqrt(n[m]))
w[m] -= lr * g[m]
wsum += w[m]
w[modelNum-1] = 1-wsum
for x in w:
print x
输入数据的格式如下：
label score11,score12,...,score1K score21,score22,...,score2K ... scoreM1,scoreM2,...,scoreMK
举个例子，比如二分类4个模型的数据片段如下：
0 0.912427,0.0875725 0.905673,0.0943273 0.911398,0.088602 0.933166,0.066834
0 0.86061,0.139389 0.865925,0.134075 0.881552,0.118448 0.872196,0.127804
1 0.364299,0.635701 0.32974,0.67026 0.323839,0.676161 0.341047,0.658953
0 0.715563,0.284437 0.713995,0.286005 0.713599,0.286401 0.713383,0.286617
0 0.948186,0.0518137 0.925587,0.0744127 0.929451,0.0705489 0.939952,0.0600485
0 0.58531,0.41469 0.588321,0.411679 0.520456,0.479544 0.671846,0.328154
1 0.0154455,0.984554 0.0150741,0.984926 0.0118866,0.988113 0.0109413,0.989059
1 0.0472074,0.952793 0.0612501,0.93875 0.065446,0.934554 0.061947,0.938053
1 0.430228,0.569772 0.366636,0.633364 0.337206,0.662794 0.5487,0.4513
0 0.97958,0.0204195 0.980348,0.0196517 0.979461,0.0205394 0.985608,0.0143915
我们的测试数据共200多万条，时间复杂度跟单模型数M呈线性关系，跑一次也就分分钟的事。为了验证算法，做了对比，对于3模型融合跑出的结果和暴力求解的结果一致。
问题似乎就这么完美解决了，但人生就像直流充电桩，上一秒还是70A的电流，下一秒就可能跳电。当我们融合更多模型时问题出现了，有的权重算出了负值！就像下面这样：
0.120181463777 gender_cut_ftrl8000_val
0.0929962960391 gender_cut_ftrl_val
0.0245135332374 gender_nn_1005_val
0.269705618425 gender_cut_ftrl2000_val
0.422565675064 gender_cut_ftrl1000_val
0.233316720486 gender_newFM_val
0.0228310140994 gender_xxFM_val
-0.186110321128 gender_result_fm
导致最终模型的预测概率有可能大于1或小于0，这如何能忍，必须再创新招！
第三式：生死符！！！
“这生死符一发作，一日厉害一日，奇痒剧痛递加九九八十一日，然后逐步减退，八十一日之后，又再递增，如此周而复始，永无休止。每年我派人巡行各洞各岛，赐以镇痛止痒之药，这生死符一年之内便可不发。”
好啦，不吓你啦，这是金庸老先生笔下天山童姥的生死符。我这里的生死符简单得很，顾名思义，由权重符号决定单模型生死，正的留，负的弃，留下的单模型重新计算权重，有可能又有新的负值出现，再次使用生死符，直到剩下的单模型权重全部大于等于0。
还以上面那次融合为例，最终把gender_xxFM_val和gender_result_fm都干掉了，剩下的权重为：
0.122788741689 gender_cut_ftrl8000_val
0.0782434412719 gender_cut_ftrl_val
0.0331138038226 gender_nn_1005_val
0.24358375583 gender_cut_ftrl2000_val
0.446465338463 gender_cut_ftrl1000_val
0.0758049189234 gender_newFM_val
在测试集上的logloss = 0.435116，而其中单模型最好的logloss是0.436593。
相应地，这批融合在7分类问题age上效果更加明显，从最好的单模型1.35416降到1.35061。
以上便是比赛期间我们使用的全部招式。
憋走，故事还没结束，忘了塞纳河畔的拉格朗日了吗？
像我这等追求卓越的好青年岂能够留下不完美的话柄，虽然比赛结束了，我还是要把它彻底解决掉——快使出
第四式：铁锁横江！！！！
分析上面的模型，出现负值，是因为少了对w的非负约束，那就加上：
$$
\min_{w}f(w)=\min_{w}\{-\frac{1}{N}\sum_{n=1}^N\sum_{k=1}^{K}1\{y_n=k\}\ln\hat{y}_{nk}\}\\
=\min_{w}\{-\frac{1}{N}\sum_{n=1}^N\sum_{k=1}^{K}1\{y_n=k\}\ln(\sum_{m=1}^{M-1}w_m\hat{y}_{nmk}+(1-\sum_{m=1}^{M-1}w_m)\hat{y}_{nMk})\}
$$
$$
s.t.\quad w_m\ge 0,\quad m=1,2,…,M-1\\
1-\sum_{m=1}^{M-1}w_m\ge 0\qquad\qquad
$$
这就成了带约束的优化问题，该怎么解呢？是时候请出拉格朗日大神了！大神蜜汁微笑，抛出一套对偶宝典，将带约束的极小化问题改造成极大极小问题，正是
第五式：沧海一粟！！！！！
我们先把问题标准化，假定$f(x)$，$c_i(x)$，$h_j(x)$ 是定义在$R^n$上的连续可微函数，将下面约束最优化问题
$$
\min_{x}f(x)\\
s.t.\quad c_i(x)\le 0,\quad i=1,2,…,k\\
\quad \quad h_j(x)=0,\quad j=1,2,…,l
$$
称为原始问题。
然后引入广义拉格朗日函数
$$
L(x,\alpha,\beta)=f(x)+\sum_{i=1}^{k}\alpha_{i}c_i(x)+\sum_{j=1}^{l}\beta_{j}h_j(x)
$$
其中 $\alpha_i$ 和 $\beta_j$ 是拉格朗日乘子，$\alpha_i\ge 0$。然后经过一番推来导去眉来眼去，可以证明
$$
\min_x\max_{\alpha,\beta:\alpha_i\ge 0}L(x,\alpha,\beta)
$$
与原始问题具有相同的x最优解，称作广义拉格朗日函数的极小极大问题。
又有，当$f(x)$，$c_i(x)$，$h_j(x)$ 等满足一定条件时
$$
\max_{\alpha,\beta:\alpha_i\ge 0}\min_xL(x,\alpha,\beta)
$$
也与原始问题具有相同的x最优解，称作广义拉格朗日函数的极大极小问题，可以将此问题也改写成约束的形式：
$$
\max_{\alpha,\beta}\theta_D(\alpha,\beta)=\max_{\alpha,\beta}\min_xL(x,\alpha,\beta)\\
s.t.\qquad \alpha_i\ge 0,\quad i=1,2,…,k
$$
称作原始问题的对偶问题。
这一部分的具体细节懒得写了，估计你们也看不下去，真要有兴趣可以去看书，比如《凸优化》或者李航老师的《统计学习方法》，我上面的公式基本就是抄他的，但请注意书上附录C的KKT条件有误，(C.22)式和(C.23)式不应该包含在里面。
好啦，照此框架，来解决我们的问题：
$$
\max_{\alpha,\alpha_m\ge 0}\min_wL(w,\alpha)\\
=\max_{\alpha,\alpha_m\ge 0}\min_w\{-\frac{1}{N}\sum_{n=1}^N\sum_{k=1}^{K}1\{y_n=k\}\ln(\sum_{m=1}^{M-1}w_m\hat{y}_{nmk}+(1-\sum_{m=1}^{M-1}w_m)\hat{y}_{nMk})-\sum_{m=1}^{M-1}\alpha_mw_m-\alpha_M(1-\sum_{m=1}^{M-1}w_m)\}
$$
极大极小公式列出来很容易，如何求解才头疼。
最早接触拉格朗日对偶是看SVM的推导，最近研究在线分配的shale算法又遇到它，发现基本都是依靠KKT条件推导，然而后面各有各的玩法，没有什么普适的方案。我也试着从KKT出发，但始终没找到什么高效的方法，一赌气，干脆舍弃KKT直接求解极大极小问题。如果哪位高人有更高明的招式，请一定传授小弟。
我自己想的招式便是以不变应万变，依然按sgd的思路，每来一条样本，先固定 $\alpha$ 不动，更新w，这是关于w的极小化问题，使用梯度下降（依然梯度杀）；然后固定w，更新 $\alpha$，这是关于 $\alpha$ 的极大化问题，使用梯度上升（可称梯云纵）。对 $\alpha$ 的非负约束也很简单，更新后如果小于0，则置为0。
如此交错曲折，在解空间的山坡上忽上忽下苦苦寻觅，“路漫漫其修远兮”，正是
第六式：上下求索！！！！！！
具体的梯度计算如下：
$$
\frac{\partial l_n(w|\alpha)}{\partial w_m}=-\sum_{k=1}^{K}1\{y_n=k\}\frac{\hat{y}_{nmk}-\hat{y}_{nMk}}{\hat{y}_{nk}}-\alpha_m+\alpha_M,\quad m=1,…,M-1\\
\frac{\partial l_n(\alpha|w)}{\partial \alpha_m}=-w_m,\quad m=1,…,M-1\qquad\qquad\qquad\qquad\qquad\qquad\\
\frac{\partial l_n(\alpha|w)}{\partial \alpha_M}=-(1-\sum_{m=1}^{M-1}w_m)\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad
$$
adagrad代码如下：
1 | import sys, math |
看一下效果，还是以上面的gender二分类为例，结果如下：
0.123262252313 gender_cut_ftrl8000_val
0.0782768467103 gender_cut_ftrl_val
0.0340907260294 gender_nn_1005_val
0.243746308063 gender_cut_ftrl2000_val
0.447138729994 gender_cut_ftrl1000_val
0.0627800144873 gender_newFM_val
0.00261889445352 gender_xxFM_val
0.00808622794875 gender_result_fm
可以看到gender_xxFM_val和gender_result_fm的权重终于不再是负数，而且很小，接近于0，跟期望的一样，最后的logloss也基本一致，等于0.435125。
说明我们的方法很给力啊！不能骄傲，继续前行。上面的方法对最后一项权重$w_M$ 做了特殊处理，不够优雅，如果我们不专门处理它，而是把权重之和归一放在约束里呢？
第七式：万法归一！！！！！！！
重新设计模型
$$
\hat{y}_{nk}=\sum_{m=1}^{M}w_m\hat{y}_{nmk}
$$
$$
\min_{w}f(w)=
\min_{w}\{-\frac{1}{N}\sum_{n=1}^N\sum_{k=1}^{K}1\{y_n=k\}\ln\hat{y}_{nk}\}\\
=\min_{w}\{-\frac{1}{N}\sum_{n=1}^N\sum_{k=1}^{K}1\{y_n=k\}\ln(\sum_{m=1}^{M}w_m\hat{y}_{nmk})\}\\
s.t.\qquad w_m\ge 0,\quad m=1,2,…,M\\
\sum_{m=1}^{M}w_m=1\qquad
$$
$$
L(w,\alpha,\beta)=-\frac{1}{N}\sum_{n=1}^N\sum_{k=1}^{K}1\{y_n=k\}\ln(\sum_{m=1}^{M}w_m\hat{y}_{nmk})-\sum_{m=1}^{M}\alpha_mw_m+\beta(\sum_{m=1}^{M}w_m-1)
$$
梯度计算
$$
\frac{\partial l_n(w|\alpha,\beta)}{\partial w_m}=-\sum_{k=1}^{K}1\{y_n=k\}\frac{\hat{y}_{nmk}}{\hat{y}_{nk}}-\alpha_m+\beta,\quad m=1,…,M\\
\frac{\partial l_n(\alpha|w,\beta)}{\partial \alpha_m}=-w_m,\quad m=1,…,M\qquad\qquad\qquad\qquad\qquad\\
\frac{\partial l_n(\beta|w,\alpha)}{\partial \beta}=\sum_{m=1}^{M}w_m-1\qquad\qquad\qquad\qquad\qquad\qquad\qquad
$$
代码
1 | import sys, math |
效果
0.117368795386 gender_cut_ftrl8000_val
0.0736857356383 gender_cut_ftrl_val
0.0387213072901 gender_nn_1005_val
0.236047170656 gender_cut_ftrl2000_val
0.463029142921 gender_cut_ftrl1000_val
0.069749655191 gender_newFM_val
0.00104464602712 gender_xxFM_val
0.00131312041348 gender_result_fm
注意，这里权重之和并不是严格的1，而是1.000959573523，毕竟是个迭代算法，没法保证严格的约束啊……所以计算logloss还是要重新归一下，最后logloss = 0.43512，跟前面的结果基本一致。
虽然最后两项的权重已经很接近于0，但看着还是很不爽，何不把很小的权重截成严格的0，就像生死符那样爽。对老司机来说都不是问题，上FTRL（为什么又是我）就好了嘛。
第八式：截权道！！！！！！！！
FTRL的稀疏性就是专职干这的，将w的adagrad改成FTRL，代码：
1 | import sys, math |
注意这里L1正则要设的很大才行，比如100，效果如下：
0.117563109491 gender_cut_ftrl8000_val
0.0734000016105 gender_cut_ftrl_val
0.0378296537773 gender_nn_1005_val
0.234723810866 gender_cut_ftrl2000_val
0.461201320728 gender_cut_ftrl1000_val
0.0747847726281 gender_newFM_val
0 gender_xxFM_val
0 gender_result_fm
哈哈，最后两项终于被彻底干掉，最后的logloss = 0.435117，依然保持的很好。
以上所有方法，在gender的7分类上也做了实验，同样有效，懒得再贴结果。
故事终于讲完了，Ensembling met Lagrange, and they lived happily ever after.
上一篇转到微信公众号时加了段作者简介：
本文作者：张卫，8年多码农生涯，先后在腾讯、搜狗混过些时日，目前在小米负责广告算法。无甚特别，唯好数学物理，正所谓推公式无敌手，推妹子无得手的中二汉子。
结果被人吐槽成了征婚帖，好吧，这次改一改：
本文作者：张卫，8年多码农生涯，先后在腾讯、搜狗混过些时日，目前在小米负责广告算法。曾经年少轻狂，颇好武侠旧学，晨舞葬诗魂的冷月刀，暮弄渡鹤影的寒塘剑；而今终日代码公式为伍，间或寄情文字，重举旧时觞。