CRF
许多随机变量组成一个无向图G = {V, E}，V代表顶点，E代表顶点间相连的边，
每个顶点代表一个随机变量，边代表两个随机变量间存在相互影响关系(变量非独立)，
如果随机变量根据图的结构而具有对应的条件独立性，
具体来说，两个没有边连接随机变量V1、V2，在其它随机变量O都确定的情况下，是独立的。
即 P(V1, V2 | O) = P(V1 | O) * P(V2 | O)
那么这被称为【成对马尔科夫性】，另有不同定义的【局部马尔科夫性】、【全局马尔科夫性】，它们互为充要条件（此处无证明）
对于满足成对马尔科夫性无向图，可以对最大团作定义
一个最大团是一组随机变量（顶点）的集合，且要满足两个条件
(1).这一组顶点之间，两两都有边相连
(2).任意一个不在这组内的顶点，不能和该组顶点的每一个都有边相连
那么一个无向图G，可以唯一地表示为一个最大团的集合C = {C1, C2, ...}
我们可能会对这组随机变量的联合分布感兴趣，比如计算P(V1=a,V2=b,V3=c...)
可以证明，无向图G的联合分布，可以被最大团表示为 phi1(C1)phi1(C2)....phin(Cn) / Z
其中，phi1~phin称为最大团C1~Cn上的势函数，Z是所有可能的随机变量取值组合下，phi1(C1)phi1(C2)....phin(Cn)的和
可以看出来Z实际上就是一个归一项
因为势函数一般要求是严格正的，所以会用一种指数函数的形式来表示
即phi1(C1)phi1(C2)....phin(Cn) = exp(E1(C1))exp(E2(C2))...exp(En(Cn)) = exp(En(Cn) + E1(Cn) + ... + En(Cn))
其中这个E1~En，可以看做是对某个最大团的随机变量的当前取值的打分
总结起来，如果要求P(Y1Y2...Yn)的值，则应该计算当前每个最大团的分数，求和，并执行softmax，softtmax的底是所有可能的变量值组合。
linear-CRF
线性的CRF，每个Y都只和前一个或后一个随机变量相连，
如 Y1——Y2——Y3——...Yn
每个最大团都是邻近的两个随机变量，如(Y1, Y2)、(Y2、Y3)等
线性CRF通常用于序列预测中，
比如，假设输入为X，输出为序列Y，每一个随机变量的取值都在T = {'B', 'E', 'M', 'O'}之中，
那么计算某个特定序列的条件概率。
更具体地，比如计算
P('BMEOBMEO' | X) / P(Y | X) = exp(score('BMEOBMEO', X)) / ∑Y exp(score(Y), X)
∑Y exp(Y | X)这项作为归一项，实际是求同样长度的所有序列组合的exp分数之和，
对于目标序列'BMEOBMEO'，它长度为8，那么同样长度的序列包括'BBBBBBBB'、'BBBBBBBM'等等....
linear-CRF的前向计算
这个同样长度的所有序列，数量是非常巨大的，
确切地说，应该是O(|T|^n)的量级，|T|是状态集的大小，n是序列长度。
如果每一个序列都要计算一次分数，那稍长一点的序列计算时间都会长到无法接受。
此时，根据linear-CRF图结构的特性，可以采用动态规划的方式，减少重复计算量，降低时间复杂度。
考虑Score(Y1:n-1, X)和Score(Y1:n, X)的区别，在两个无向图中，后者比前者增加了两条边，
(1).边Yn-1——Yn (2).边X——Yn
所以Score(Y1:n, X) = Score(Y1:n-1, X) + E-tran(Yn-1, Yn) + E-emmi(X, Yn)
其中E(Yn-1, Yn)又可称为转移分数，E(X, Yn)又可称为发射分数。
发现了该递推关系以后，再思考一个非常重要的递推中间变量：
记 g(k, tag) = exp(Score(Y1:k-1Yk = tag, X))
也就是，在时间点k上，Yk等于tag，但是前面的Y1:k-1状态随意的分数，
因为Y1:k-1状态随意，它实质上是Y1:k-1字序列所有可能组合，且Yk以tag为结尾的分数和
具体地，g(3, 'B') = ∑(Y1:k-1) exp(Score(Y1:Y2Y3='B', X)) = exp(Score('BBB', X)) + exp(Score('BEB', X)) + exp(Score('BMB', X)) + ... + exp(Score('OOB', X))
最后一个等式，很容易计算到，有 4 * 4 = 16个加和项
定义了g(k, tag)，很重要一个点就是计算递推关系，
即，假如知道了g(k, tag)，那对于计算g(k+1, Yk+1)又有什么帮助呢？首先展开g(k+1, Yk+1)
g(k+1, Yk+1) = ∑(Y1:k) exp(Score(Y1:kYk+1, X))
= ∑(Y1:k) exp[Score(Y1:k, X) + E-tran(Yk, Yk+1) + E-emis(X, Yk+1)]
= ∑(Y1:k) exp[Score(Y1:k, X)]exp[E-tran(Yk, Yk+1)]exp[E-emis(X, Yk+1)]
= ∑(Yk)∑(Y1:k-1) exp[Score(Y1:k-1Yk, X)]exp[E-tran(Yk, Yk+1)]exp[E-emis(X, Yk+1)]
= ∑(Yk) g(k, Yk)exp[E-tran(Yk, Yk+1)]exp[E-emis(X, Yk+1)]
这样，就产生了g(k, tag)的递推公式，也举一个简单的例子
g(3, 'B') = g(2, 'B')exp(E-tran('B', 'B'))exp(E-emis(X, 'B')) + g(2, 'E')exp(E-tran('E', 'B'))exp(E-emis(X, 'B')) + g(2, 'M')exp(E-tran('M', 'B'))exp(E-emis(X, 'B')) + g(2, 'O')exp(E-tran('O', 'B'))exp(E-emis(X, 'B'))
在这种情况下，要计算全序列的exp分数和，就要先计算每个时间点上，以某个状态结束的全序列exp分数和，并向后递推
每次递推都要用前一个时间的全状态，组合后一个时间的每个状态，故递推一次时间复杂度为O(|T|^2)
总时间复杂度为O(n|T|^2)，比起强行计算的O(|T|^n)大大减少
linear-CRF的数值优化
我们得到了重要的递推公式
g(k+1, Yk+1) = ∑(Yk) g(k, Yk)exp[E-tran(Yk, Yk+1)]exp[E-emis(X, Yk+1)]
但是，大量exp的相乘，可能导致浮点运算的数值不稳定，故在代码中，会采取一些方法优化数值稳定性
首先求和可以写成对向量的求和
sum(g(k, Yk+1)exp(E-tran('B', Yk+1))exp(E-emis(X, Yk+1)) + g(k, 'E')exp(E-tran('E', Yk+1))exp(E-emis(X, Yk+1)) + g(k, 'M')exp(E-tran('M', Yk+1))exp(E-emis(X, Yk+1)) + g(k, 'O')exp(E-tran('O', Yk+1))exp(E-emis(X, Yk+1)))
进一步，调研向量中的每个分量，以第二个为例
g(k, Yk+1)exp(E-tran('B', Yk+1))exp(E-emis(X, Yk+1))
= exp(log(g(k, 'E')exp(E-tran('E', Yk+1))exp(E-emis(X, Yk+1)))) #先log后exp，数值还是一样的
= exp(log(g(k, 'E')) + E-tran('E', Yk+1) + E-emis(X, Yk+1)) #log掉后两项自带的exp，直接可以用分数相加
这个时候，可以发现原来乘积的第一项，已经由g函数变成了log g函数，但等式左边还是g函数
如果想要保持两边一致，不妨左侧也加上log，就会变成
log g(k+1, Yk+1) = log(sum(exp(vector)))，其中vector = (log g(k, Yk+1)+E-tran('E', Yk+1) + E-emis(X, Yk+1))
右侧这个log(sum(exp(xxx)))的计算，在pytorch里有API，直接可以看torch.logsumexp的文档
实际上，这里都还没有讲到真正数值稳定计算的部分，只是做了恒等变换，变成方便观察的形式,
对于logsumexp的数值稳定计算，假设其内部的vector为(v1, v2, ...vn)，且max(v1, v2, ...vn) = vmax，可进行如下变形
logsumexp([v1, v2, ..., vn])
= log(sum(exp([v1-vmax+vmax, v2-vmax+vmax, ..., vn-vmax+vmax])))
= log(sum([exp(vmax)exp(v1-vmax), exp(vmax)exp(v2-vmax), ... , exp(vmax)exp(vn-vmax)]))
= log(exp(vmax) * sum([exp(v1-vmax), exp(v2-vmax), ..., exp(vn-vmax)]))
= vmax + log(sum(exp[v1-vmax, v2-vmax, ..., vn-vmax]))
这样子，(v1, v2, ..., vn)中的最大分量vmax将会被直接提出来，不进行一遍exp后再log的计算，
其余的较小值进行exp的计算，使得数值偏差尽可能小，达到数值稳定的效果
linear-CRF的学习
在此前linear-CRF的计算过程中，我们可以使用递推公式(状态转移方程)，使得计算指定序列P(Y | X)的时间复杂度降为O(nT^2)
但是，使用递推公式要求E-tran和E-emis函数是已知的
如果仅有观测数据集(X, Y)对，却不值得E-tran和E-emis，应如何计算E-tran和E-emis呢
实际上，当知道确切E-tran、E-emis时，计算出来的exp(score('BMEOBMEO', X)) / ∑Y exp(score(Y), X)可以被认为是条件概率
而不知道E-tran、E-emis时，根据初始化参数计算得的exp(score('BMEOBMEO', X)) / ∑Y exp(score(Y), X)则可以被认为是似然
求解E-tran、E-emis的过程，即对具体的数据集(X, Y)求使得exp(score(Y, X)) / ∑Y' exp(score(Y'), X)最大的E-train、E-emis参数
其实就是最大似然，可以使用梯度下降等优化方法求解。实际操作中，会去最小化负对数似然，即，
minimize log(∑Y' exp(score(Y'), X)) - score(Y, X)
那具体的参数应该如何设置呢，按照流行的LSTM + linear-CRF的方法，
E-tran被直接设置为一个 T * T 的参数矩阵，矩阵代表了序列状态间的转移分数
E-emis则不直接设置为参数，而是由LSTM + dense_layer把X映射成一个feature，
feature的维度是 n * T，第一个维度代表时间(序列长度)，第二维度代表每个时间状态上都有T个发射分数，
这T个发射分数则对应每个时间段的E-emis
即E-emis(Yk, X) = E-emis(Yk, feature) = feature[k, tags.index(Yk)]
这样，E-emis的参数便被暗含在LSTM + dense_layer的参数中
LSTM + linear-CRF的动机
既然如此，那么可能出现疑问，为什么E-tran可以单独设置，E-emis却不单独设置呢？
换句话说，既然单纯的linear-CRF就能完整解决一个序列预测的问题，为什么还要在前面加一个LSTM呢？
这个问题其实又可以反过来问，既然单纯的LSTM已经可以完整解决一个序列预测的问题，为什么还要在后面加一个CRF呢？
实际上，这两种方法的能力各有侧重，有刚刚好可以互补，所以在序列预测中被放在一起，成为了非常流行的结构。
LSTM的优点是具有非常强大的特征提取能力，因为它可以使用许多非线性函数、上下文关系等，构建出高维又稠密的非线性特征，
如果不使用LSTM的特征提取，像E-trans一样，设置一个纯粹参数化的E-emis矩阵，
那这种矩阵相乘的特征提取是线性的，且同时还是稀疏的，也无法结合X序列的上下文关系，效果差距非常非常大。
LSTM虽然具有以上特点，但是如果仅仅使用它进行序列预测，它则缺少一个显式地对转移关系的建模。
很明显地，LSTM只能模拟发射函数，却不能模拟转移函数，它可以提取复杂的特征，但是有时会对不合理的局部过于宽容。
例如，BMEO分词标注中(可以是分词标注、实体标注，以分词为例)，
B代表词语开始，E代表词语结束，M代表词语中部，O代表单字词，
逻辑上而言，一组距离最近的B和E的中间只能存在M，而像BB、MB这种局部结构，是绝对错误的，
但LSTM更关注总体的loss最小化，不显式地对邻近转移关系进行建模，所以总是免不了出现这样的局部小错
linear-CRF则带来了一个转移矩阵参数，能够很有效地解决这个问题。
这就是LSTM与linear-CRF配合使用的原因，一句话总结，LSTM负责特征提取，linear-CRF负责邻近转移关系管理。
但是写到这里还有一个有趣的思考：虽然linear-CRF要注重转移关系，
但是E-trans是否可以像E-emis一样，使用一个复杂网络来模拟，而不直接使用一个线性矩阵呢？
这个似乎是一个很有趣的话题，到时要去算算间复杂度是否能允许这种情况下的前向计算，说不定是一个很有趣的点？
linear-CRF的预测
在linear-CRF的计算问题中，刚刚已经讨论了两种，一个是概率计算问题，一个是参数学习问题
概率计算问题知道E-tran和E-emis，要计算P(Y | X)，
参数学习问题知道一批(X, Y)，希望参数估计E-tran，E-emis，
最后剩下一个预测问题，也称解码问题，模型训练好以后，得到了估计的E-tran，E-emis，此时来一个新的X，最可能的Y是什么？
即，预测问题，知道E-tran，E-emis，X，要计算argmaxY P(Y | X)
这个式子可以作一些简化，
argmaxY P(Y | X)
= argmaxY exp(score(Y, X)) / ∑Y exp(score(Y), X)
= argmaxY exp(score(Y, X)) #normalize项对Y的取值不影响
= argmaxY score(Y, X) #exp是单调函数
这三个问题，实际上跟HMM的三个问题也是完全对应的。
linear-CRF的预测问题，实际上和概率计算问题是极其类似的，
只不过前者是求最大的概率路线，后者是前全路线的概率和，
同样地，如果不使用动态规划方法，枚举的时间复杂度是O(T^n)，使用后时间复杂度是O(nT^2)
然而最大概率路线的算法还有个独有的名字，也就是平常非常常见的viterbi算法233333
在前向概率计算时，我们定义的递推公式单元是
g(k, tag) = exp(Score(Y1:k-1Yk = tag, X))
思想是，只要知道了k时间点Yk为特定状态的全路线概率，就可以知道k+1时间点Yk+1位特定状态的全路线概率
同理，定义一个 h(k, tag) = argmax Y1:k-1 score(Y1:k-1Yk=tag, X)
思想是，只要知道了k时间点Yk为特定状态的最大概率路线，希望能够递推出k+1时间点Yk+1位特定状态的最大概率路线
假设，后者这个路线，在k时间点Yk要经过状态t，且k时间点Yk经过状态t的最大概率路线是Yt1:k-1，
那么后者的路线必定是Yt1:k-1,Yk=t,Yk+1，否则，不存在更大概率的路线。
h(k+1, Yk+1)
= argmax Y:k score(Y1:kYk+1, X)
= argmax Y:k score(Y1:kYkYk+1, X)
= argmax Y:k-1Yk score(Y1:k-1Yk, X) + E-tran(Yk, Yk+1) + E-emis(Yk, X)
= Y1:k-1argmaxYk (argmax Y:k-1 score(Y1:k-1Yk, X) + E-tran(Yk, Yk+1) + E-emis(Yk, X))
= Y1:k-1argmaxYk (h(k, Yk) + E-tran(Yk, Yk+1) + E-emis(Yk, X))
这样就可以求得最大的概率路线