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))
這樣就可以求得最大的概率路線