由於博客園對markdown支持不完善（或者我不太會用），一些公式和引用展示不正常或不好看，建議下載pdf版閱讀
https://github.com/AlexPengW/resource/tree/master/pdf
Abstract
本文首先介紹CRF[1]（條件隨機場）模型及用來提高CRF求解和預測效率的Forward-Backward 算法，然后介紹經典的無約束優化算法Line Search[2]，DFP[3]，BFGS，L-BFGS[4], IIS[5]，給出一些必要性的證明，最后結合CRF++代碼講解CRF和L-BFGS的代碼實現。
CRF模型
CRF模型概述
CRFs are a type of discriminative undirected probabilistic graphical model. It is used to encode known relationships between observations and construct consistent interpretations. It is often used for labeling or parsing of sequential data, such as natural language text or biological sequences and in computer vision. —— Wikipedia
本文主要描述Linear-chain CRF，
記\(P(y|x,\lambda )\)為CRF模型定義的條件概率， 即給定觀測序列\(x\)時，任意標注序列\(y\)的條件概率,則\(P(y|x,\lambda )\)定義如下：
\begin{equation}
\begin{cases}
P\left(y|x,\lambda \right)&=\dfrac{1}{Z(x)}exp\left(\sum_{i,k}\lambda_kt_k\left(y_{i-1},y_i,x,i\right)+\sum_{i,l}\mu_ls_l\left(y_i,x,i\right)\right) \
Z\left(x\right)&=\sum_yexp\left(\sum_{i,k}\lambda_kt_k\left(y_{i-1},y_i,x,i\right)+\sum_{i,l}\mu_ls_l\left(y_i,x,i\right)\right)
\end{cases}
\end{equation}
其中\(s_l(y_i,x,i)\)為狀態特征函數，\(t_k(y_{i-1},y_i,x,i)\)為轉移特征函數。特征函數由用戶定義，CRF++代碼里面用Unigram模板提取狀態特征，用Bigram模板提取轉移特征，對應的特征函數返回0或1。\(\mu_l\)和\(\lambda_k\)分別是對應的特征權值，即待求解的模型參數。
上面的概率公式具體來歷可以參考統計學習方法[6]書中的最大熵模型，或者這篇博客最大熵模型[7]，這里不再敘述。
為了簡化公式，記\(s(y_i,x,i)=s(y_{i-1},y_i,x,i)\), \(f(y_{i-1},y_i,x,i)\)或者是個狀態函數\(s(y_{i-1},y_i,x,i)\)或者是個轉移函數\(t(y_{i-1},y_i,x,i)\)，記\(N\)為某個樣本序列長度，不同樣本序列長度可能不同，記\(L\)為標記序列\(y\)在任意位置可能取到的標記個數，記\(T\)為特征數量。則\(P(y|x,\lambda )\)也可以定義如下：
\begin{equation}\begin{cases}
P\left(y|x,\lambda \right)&=\dfrac{1}{Z(x)}exp\left(\sum_j\lambda_jF_j\left(y,x\right)\right)\
F_j\left(y,x\right)&=\sum_{i=1}^Nf_j\left(y_{i-1}, y_i,x,i\right)\
Z(x)&=\sum_yexp\left(\sum_j\lambda_jF_j\left(y,x\right)\right)
\end{cases}\end{equation}
CRF模型求解
給定訓練樣本\(\{(x^m,y^m)\}, m =\{1,...,M\}\), 則CRF模型的對數似然函數為：
\begin{equation}
\mathcal{L}(\lambda)=\sum_mlogP\left(y^m|x^m,\lambda\right)=\sum_m\left(\sum_j\lambda_jF_j\left(x^m,y^m\right)-logZ\left(x^m\right)\right)
\end{equation}
可以證明上面的函數是關於\(\lambda\)的凹函數，我只試了一個不優美的證明方法，所以就不列出來了，簡單來說就是利用[定理A function is convex if and only if it is convex when restricted to any line that intersects its domain，這個定理也容易證明]將高維上的證明轉化為一維上的證明，然后證明一維函數二階導數小於0即可。通常我們求解時都會在后面減去\(L_2\)正則化項\(\dfrac{\sum_j\lambda_j^2}{2C}\)，減去\(L_2\)正則化項后就是嚴格凹函數(strictly concave)，由嚴格凹函數性質知\(\mathcal{L}(\lambda)\)有唯一一個全局最大值，且最大值處導出為\(0\)。
將對數似然對\(\lambda_j\)求導得：
\begin{equation}
\left{\begin{aligned}
\dfrac{\partial\mathcal{L}(\lambda)}{\partial\lambda_j}&=\sum_m\left(F_j\left(x^m,y^m\right)-\dfrac{\sum_y\left(F_j\left(y,x^m\right)exp\left(\sum_j\lambda_jF_j\left(y,x^m\right)\right)\right)}{Z(x^m)}\right)\
&=\sum_m\left(F_j\left(x^m,y^m\right)-\sum_yP(y|x^m,\lambda)F_j\left(y,x^m\right)\right)\
&=\sum_m\left(F_j\left(x^m,y^m\right)-E_{P(Y|x^m,\lambda)}[F_j\left(Y,x^m\right)]\right)\
&=M\sum_m\left(\widetilde{P}(y^m,x^m)F_j(y^m,x^m)-\widetilde{P}(x^m)\sum_yP(y|x^m,\lambda)F_j(y,x^m)\right)\
&=M\left(E_{\widetilde{P}(Y,X)}[F_j(y^m,x^m)]-E_{P(Y,X)}[F_j(Y,x^m)]\right)
\end{aligned}
\right.
\end{equation}
注意上面的\(\widetilde{P}(Y,X)\)是X，Y的聯合經驗概率分布， 概率之和為1，所以左邊有乘以\(M\)，\(P(Y,X)\)是模型分布（\(P(y,x)=\widetilde{P}(x)P(y|x,\lambda)\)），CRF++代碼是按照\(\sum_m\left(F_j(x^m,y^m)-E_{P(Y|x^m,\lambda)}[F_j(Y,x^m)]\right)\)來計算梯度的(准確的說CRF++會對似然函數和導數取負號，轉化為最小化問題)。
上面的梯度式子中，對於每個樣本要計算特征\(F_j(Y,x^m)\)在模型條件分布\(P(Y|x^m,\lambda)\)上的期望\(E_{P(Y|x^m,\lambda)}[F_j\left(Y,x^m\right)]\)。記\(g\)為梯度數組。
本文采用一種與論文不同的思路來推引出Forward Backward算法，寫的比較繞。
暴力求解
假設用最直接的方法求梯度\(\sum_m\left(F_j\left(x^m,y^m\right)-\sum_yP(y|x^m,\lambda)F_j\left(y,x^m\right)\right)\)，對於后面一項，\(y\) 有\(L^N\)種不同取值，對於每個取值\(y\), 計算\(P(y|x^m,\lambda)\)，然后將計算結果累加到梯度中，即：
\(F_j(y,x^m)\)即第 \(j\)個特征在 \((y, x^m)\)中出現的次數。計算 \(P(y|x^m,\lambda)\)概率復雜度為 \(O(NT)\)(這里假設提前分別統計存儲每個位置上的特征，所以要遍歷每個位置，共 \(N\)個位置，每個位置最多可能 \(T\)個特征。只考慮非規范化概率，即不考慮 \(Z(x)\)（這項所有 \(y\)都一樣，可以在非規范化概率計算完后，累加即得到，不影響整體計算復雜度）；后面將概率累加到數組 \(g\)，復雜度也是 \(O(NT)\)(如前面所說每個位置特征單獨存儲，所以單獨相加，最差情況下這個復雜度是 \(O(NT)\), 當然實際上，由於 \(P(y|x^m,\lambda)\)計算時涉及計算指數和對數，實際計算時間會是 \((g[])\)的倍數）。所以對於每個樣本總體計算復雜度是
\(O(Y集合大小)*(O(P(y|x^m,\lambda)+O(g[]))=O（L^N*N^T）\)
為了節省CPU，考慮下怎么加速，加速方法一般就是減少重復計算。很容易可以發現對於兩個序列 \(y^1,y^2\), 如果他們只有最后一個位置的標記不同，那么 \(F_j(y,x)=\sum_{i=1}^Nf_j(y_{i-1}, y_i,x,i)\)的前 \(N-1\)部分累加和是相等的, 如果枚舉 \(Y\)時按照一定的順序，使每次 \(y\)只在一個位置發生變化，那 \(O(P(y|x^m,\lambda)）\)的計算時間一般情況也能得到明顯提高，不過即使能降低 \(O(P(y|x^m,\lambda）)\)的計算復雜度，由於 \(O(g[]+=)\)並未變化， \(O(Y集合大小)\)也未變化，所以整體計算復雜度未變。 所以只要采取枚舉\(y\)的計算策略，單一靠降低單個\(y\)的\(O(P(y|x^m,\lambda）)\)的計算復雜度是不能降低整體計算復雜度的。
Forward-Backward算法
由於枚舉\(y\)的暴力計算方法無法降低總體計算復雜度，我們可以尋求將\(y\)打包計算的方法。比如如果能將\(Y\)打包成一些子集\(\{Y^q|q=\{1,...,Q\}\}\), 其中不同子集可能相交，將\(\dfrac{\partial\mathcal{L}(\lambda)}{\partial\lambda_j}\)進行適當變換，使計算復雜度可以表示成
\begin{equation}
O(Q)*(O(P(Y^q|x^m,\lambda)+O(g[]))
\end{equation}
並且使\(Q\) 遠小於\(L^N\)， 並且計算\(O(P(Y^q|x^m,\lambda)\)時不能采取之前枚舉\(y\)累加的辦法，如果枚舉累加，那么總體復雜度還跟之前一樣了，所以各個子集的選取要方便一起打包高效計算。
對\(\dfrac{\partial\mathcal{L}(\lambda)}{\partial\lambda_j}\)做如下變化
\begin{equation}
\left{\begin{aligned}
\dfrac{\partial\mathcal{L}(\lambda)}{\partial\lambda_j}&=\sum_m\left(F_j\left(x^m,y^m\right)-\sum_yP(y|x^m,\lambda)F_j\left(y,x^m\right)\right)\
&=\sum_m\left(F_j\left(x^m,y^m\right)-\sum_{i=1}^N\left(\sum_yf_j\left(y_{i-1}, y_i,x^m,i\right)P\left(y|x^m,\lambda\right)\right)\right)\
&=\sum_m\left(F_j\left(x^m,y^m\right)-\sum_{i=1}^N\left(\sum_{l,l'}\left(\sum_{{y|y_{i-1}=l \ and\ y_i=l'}}f_j\left(y_{i-1}, y_i,x^m,i\right)P\left(y|x^m,\lambda\right)\right)\right)\right)
\end{aligned}
\right.
\end{equation}
如果\(f_j\)是單特征，那么上式可以進一步化簡為
\(\dfrac{\partial\mathcal{L}(\lambda)}{\partial\lambda_j}=\sum_m\left(F_j\left(x^m,y^m\right)-\sum_{i=1}^N\left(\sum_l\left(\sum_{\{y|y_i=l\}}s_j\left(y_i,x^m,i\right)P\left(y|x^m,\lambda\right)\right)\right)\right)\)
上式中的\(l,l'\)是相鄰標記可能出現的各種組合,一共\(L^2\)種，右邊的\(\{y|y_{i-1}=l \ and\ y_i=l'\}\)是為轉移特征找到的一種\(y\)集合, 集合里面每個\(y\)滿足\(y_{i-1}=l \ and\ y_i=l'\)，一共\((N-1)L^2\)種\(y\)集合。\(\{y|y_i=l\}\)是我們為狀態特征找到的集合，一共\(NL\)種, 記\(Y^q\)為任一找到的集合，\(q\in{\{1,...,NL\}}\)為狀態特征集合，\(q\in{\{NL+1,...,(N-1)L^2+NL\}}\)為轉移特征集合，記
\(P(Y^q|x^m,\lambda)=\sum_{\{y|y\in{Y^q}\}}P(y|x^m,\lambda)\)
對於狀態特征
\(\dfrac{\partial\mathcal{L}(\lambda)}{\partial\lambda_j}=\sum_m\left(F_j(x^m,y^m)-\sum_{\{q|q=1,...,NL\}}s_j(y_i,x^m,i)P(Y^q|x^m,\lambda)\right)\)
對於轉移特征
\(\dfrac{\partial\mathcal{L}(\lambda)}{\partial\lambda_j}=\sum_m\left(F_j(x^m,y^m)-\sum_{\{q|q\in{(NL+1,...,(N-1)L^2+NL)\}}}t_j(y_{i-1},y_i,x^m,i)P(Y^q|x^m,\lambda)\right)\)
現在的問題就是尋找\(P(Y^q|x^m,\lambda)\)的高效算法。
記\(S_{y_n=l}\)為某個狀態特征\(Y\)子集，\(S_{y_{n-1}=l, y_n=l'}\)為某個轉移特征\(Y\)子集。
對於兩個狀態特征子集\(S_{y_n=a},S_{y_n=b}\), 存在一個一一映射，對於任意\(y^a\in{S_{y_n=a}}\), 存在一個唯一的\(y^b\in{S_{y_n=b}}\), 使\(y^a,y^b\)只在第\(n\)個位置上不同，在其他位置都是一樣的。考慮\(S_{y_n=a},S_{y_n=b}\)的兩個子集\(S_{y_{n-1}=c\ and\ y_n=a\ and\ y_{n+1}=d},S_{y_{n-1}=c\ and\ y_n=b\ and\ y_{n+1}=d}\), 兩個子集的元素也存在一一對應關系，使他們只在位置\(n\)不同。
記
\(\alpha_{n-1}(c)=\sum_{y\in\alpha_{n-1,c}}exp\left(\sum_j\left(\lambda_j\sum_{i=1}^{n-1}f_j\left(y_{i-1}, y_i,x,i\right)\right)\right)\)
\(\beta_{n+1}(d)=\sum_{y\in\beta_{n+1,c}}exp\left(\sum_j\left(\lambda_j\sum_{i=2}^{N-n}f_j\left(y_{i-1}, y_i,x,i\right)\right)\right)\)
\(M_n(c,a)=exp\left(\sum_j\lambda_jf_j\left(c,a,x,n\right)\right)\)
其中\(\alpha_{n-1,c}\)為所有長度為\(n-1\),終點為\(c\)的標記序列集合，\(\beta_{n+1,d}\)為所有長度為\(N-n\),起點為\(d\)的標記序列集合，\(\alpha_{n-1}(c)\)表示集合\(\alpha_{n-1,c}\)前向非規范化概率和，\(\beta_{n+1}(d)\)表示集合\(\beta_{n+1,d}\)的后項非規范化概率和。則
\(P(S_{y_{n-1}=c\ and\ y_n=a\ and\ y_{n+1}=d}|x,\lambda)=\dfrac{1}{Z(x)}\alpha_{n-1}(c)M_n(c,a)M_{n+1}(a,d)\beta_{n+1}(d)\)
上面的公式可以這樣理解：
記\(\alpha_{n,c,a}\)為長度為\(n\)且最后兩個標記分別為\(c,a\)的標記序列集合，則\(\alpha_{n,c,a}\)與\(\alpha_{n-1,c}\)兩個集合元素個數一樣，將集合\(\alpha_{n-1,c}\)里面每個元素尾部加個\(a\)就得到集合\(\alpha_{n,c,a}\)。記\(\alpha_n(c,a)\)表示集合\(\alpha_{n,c,a}\)的前向概率和，則
\(\alpha_n(c,a)=\sum_{y\in\alpha_{n,c,a}}exp\left(\sum_j\lambda_j\sum_{i=1}^nf_j\left(y_{i-1}, y_i,x,i\right)\right)\)
\(=\sum_{y\in\alpha_{n,c,a}}exp\left(\sum_j\lambda_j\left(f_j\left(c,a,x,n\right)+\sum_{i=1}^{n-1}f_j\left(y_{i-1}, y_i,x,i\right)\right)\right)\)
\(=exp\left(\sum_j\lambda_jf_j\left(c,a,x,n\right)\right)\sum_{y\in\alpha_{n,c,a}}exp\left(\sum_j\sum_{i=1}^{n-1}f_j\left(y_{i-1}, y_i,x,i\right)\right)\)
\(=M_n(c,a)\sum_{y\in\alpha_{n,c,a}}exp\left(\sum_j\sum_{i=1}^{n-1}f_j\left(y_{i-1}, y_i,x,i\right)\right)\)
\(=M_n(c,a)\alpha_{n-1}(c)\)
同理，記\(\alpha_{n+1,c,a,d}\)為長度為\(n+1\)且最后三個標記分別為\(c,a,d\)的標記序列集合,能得出
\(\alpha_{n+1}(c,a,d)=\alpha_{n-1}(c)M_n(c,a)M_{n+1}(a,d)\)
將集合\(\alpha_{n+1,c,a,d}\)每個元素擴充一下（每個元素尾部分別追加集合\(\beta_{n+1,d}\)里面的元素）便得到集合\(S_{y_{n-1}=c\ and\ y_n=a\ and\ y_{n+1}=d}\)，用類似上面推出\(\alpha_n(c,a)\)的方法，也可以推出（有些繞就不推了）
\(P(S_{y_{n-1}=c\ and\ y_n=a\ and\ y_{n+1}=d}|x,\lambda)=\dfrac{1}{Z(x)}\alpha_{n-1}(c)M_n(c,a)M_{n+1}(a,d)\beta_{n+1}(d)\)
同理可以推出
\(P(S_{y_{n-1}=c\ and\ y_n=b\ and\ y_{n+1}=d}|x,\lambda)=\dfrac{1}{Z(x)}\alpha_{n-1}(c)M_n(c,b)M_{n+1}(b,d)\beta_{n+1}(d)\)
\(P(S_{y_{n-1}=c\ and\ y_n=a\ and\ y_{n+1}=d}|x,\lambda)\)和\(P(S_{y_{n-1}=c\ and\ y_n=b\ and\ y_{n+1}=d}|x,\lambda)\)具有公共子結構\(\alpha_{n-1}(c)和\beta_{n+1}(d)\)，也即存在可以復用的子結構來提高求解速度。並且由\(\alpha,\beta\)定義知道子結構之間存在遞推關系：
\begin{equation}
\begin{cases}
\alpha_{n}(a)=\sum_l\alpha_{n-1}(l)M_n(l,a) \
\beta_{n}(a)=\sum_l\beta_{n+1}(l)M_{n+1}(a,l)
\end{cases}
\end{equation}
這樣便可以用經典的動態規划算法來提高求解效率，避免子問題重復求解。
直接由\(\alpha_{n}(a),\beta_{n}(a)\)的定義能得出
\(P(S_{y_n=a}|x,\lambda)=\dfrac{1}{Z(x)}\alpha_{n}(a)\beta_{n}(a)\)
也可以將\(S_{y_n=a}\)分解成互不相交子集合的形式累加出來
\(P(S_{y_n=a}|x,\lambda)=\sum_{l,l'}P(S_{y_{n-1}=l\ and\ y_n=a\ and\ y_{n+1}=l'}|x,\lambda)\)
\(=\dfrac{1}{Z(x)}\sum_{l,l'}\alpha_{n-1}(l)M_n(l,a)M_{n+1}(a,l')\beta_{n+1}(l')\)
\(=\dfrac{1}{Z(x)}\sum_{l'}\{\beta_{n+1}(l')M_{n+1}(a,l')(\sum_l\alpha_{n-1}(l)M_n(l,a))\}\)
\(=\dfrac{1}{Z(x)}\sum_{l'}\{\beta_{n+1}(l')M_{n+1}(a,l')\alpha_{n}(a)\}\)
\(=\dfrac{1}{Z(x)}\alpha_{n}(a)\sum_{l'}\beta_{n+1}(l')M_{n+1}(a,l')\)
\(=\dfrac{1}{Z(x)}\alpha_{n}(a)\beta_{n}(a)\)
同理可以推出
\(P(S_{y_{n-1}=l, y_n=l'}|x,\lambda)=\dfrac{1}{Z(x)}\alpha_{n}(a)M_n(l,l')\beta_{n}(a)\)
這樣便可以將各個\(Y\)子集\(S_{y_n=l}\)和\(S_{y_{n-1}=l, y_n=l'}\)的概率用\(\alpha_{n}(a),\beta_{n}(a)\)表示出來。並且由公式(7)可以用動態規划法高效求出\(\alpha_{n}(a),\beta_{n}(a)\)。
現在來根據公式(5)計算下Forward-Backward算法的總體時間復雜度。
由於多了求\(\alpha_{n}(a),\beta_{n}(a)\)的開銷，所以總體復雜度要加上\(O(\alpha_{n}(a)+\beta_{n}(a))\)
第一步先計算所有\(M_n(l,l')\),再計算\(\alpha_{n}(a),\beta_{n}(a)\)，總體計算復雜度是\(O(NL^2)\)
第二步計算\(Z(x)=\sum_l\alpha_N(l)\) 為\(O(L)\)復雜度
\(P(Y^s|x^m,\lambda)\)的計算復雜度為\(O(1)\)
集合大小\(S=(N-1)L^2+NL\)
\(O(g[]+=)\)未變，仍為\(O(T)\)
所以總體時間復雜度是\(O(TNL^2)+O(NL^2)=O(TNL^2)\)， 比暴力求解\(O(TNL^N)\)復雜度明顯優越。
后面做預測時用的Viterbi算法和Forward-Backward算法很類似，只是將遞推公式中的\(\sum\)換成\(Max\)，時間復雜度是\(O(NL^2)\)，如果用暴力方法預測，則時間復雜度為\(O(TNL^N)\)。
為了便於后續對着代碼解釋，這里將導數公式整理下，
\begin{equation}
\left{\begin{aligned}
\dfrac{\partial\mathcal{L}(\lambda)}{\partial\lambda_j}&=\sum_m\left(\sum_is_j(y_i^m,x^m,i)\right)+\sum_m\left(\sum_i\left(\sum_{l}s_j(l,x^m,i)P(S_{y_i=l}|x^m,\lambda)\right)\right)\
&=\sum_m\left(\sum_is_j(y_i^m,x^m,i)\right)+\sum_m\left(\sum_i\left(\sum_{l}s_j(l,x^m,i)\dfrac{1}{Z(x^m)}\alpha_{n}(a)\beta_{n}(a)\right)\right)\ \ 狀態特征\
\dfrac{\partial\mathcal{L}(\lambda)}{\partial\lambda_j}&=\sum_m\left(\sum_i\left(t_j(y_{i-1}^m,y_i^m,x^m,i)\right)\right)+\sum_m\left(\sum_i\left(\sum_{l,l'}\left(t_j(l,l',x^m,i)P(S_{y_{i-1}=l,y_i=l'}|x^m,\lambda)\right)\right)\right)\
&=\sum_m\left(\sum_it_j(y_{i-1}^m,y_i^m,x^m,i)\right)+\sum_m\left(\sum_i\left(\sum_{l,l'}t_j(l,l',x^m,i)\dfrac{1}{Z(x^m)}\alpha_{n}(a)M_n(l,l')\beta_{n}(a)\right)\right)\ \ 轉移特征
\end{aligned}
\right.
\end{equation}
\(S_{y_i=l}\)，\(S_{y_{i-1}=l, y_i=l'}\)的定義見上文。
Unconstrained Optimization
無約束優化問題數學描述如下：
其中 \(x\in{R^n}\)即 \(n\)維實數， \(f:R^n\rightarrow R\)是一個光滑函數, 由於 \(max_x\ f(x)=min_x-f(x)\), 所以我們只討論最小化問題。
一些無約束優化問題有 Closed-Form Solution，這樣就可以直接套公式求出解，比如 \(min_x\ \ (x-a)^2\)的解為 \(x=a\)。對於沒有閉包形式解的函數，或者閉包形式的解計算和存儲代價比較大的時候，就要尋求其他方法。
一般可以通過給定一個初始 \(x_0\), 然后不斷迭代，搜索一系列的 \(x_0,x_1,...,x_k\)，直到我們認為 \(x_k\)已經近似接近最優解，或者 \(f(x_k)\)已經足夠接近最優值，或者達到我們心里預期值，或者我們的算法在 \(x_k\)時由於迭代條件無法滿足，無法繼續迭代了（當然這時也可以換個初始點開始新一輪迭代）。
如果 \(f(x)\)是個強凸函數，由凸函數性質知道 \(f(x)\)存在唯一最小值，並且我們可以通過一定准則判斷是否達到或者接近最小值或者知道接下來進一步優化性價比比較低了，比如可以判斷導數是否足夠接近0。對於非凸函數，如果我們能推出 \(f(x)\)在定義域內有下界，記為 \(Bound\)，也可以根據 \(f(x_k)-Bound\)大小來判斷 \(x_k\)是否已經接近全局最優值，或者根據梯度判讀是否接近局部最優值。如果非凸函數沒有其他已知信息，就很難找到一個合理的判定准則，但可以根據已知問題領域的知識，估計出一個可以接受的值 \(Accept\), 根據 \(\left|f(x_k)-Accept\right|\)來判定是否已經達到預期優化目標。
Line Search
Line Search[^Line Search]顧名思義就是沿着一個線尋找下一個\(x_{k+1}\)， 即
從上面可以看出Line-Search算法分為兩步
Bracket Stage
如圖一所示，Bracket Stage或者找到一個滿足Strong Wolfe Conditions的\(\alpha_*\) ,那整個算法就迭代結束了，或者找到一個區間\((\alpha_{lo},\alpha_{hi})\)使這個區間存在滿足Strong Wolfe Conditions的步長，圖中寫的算法如果上面兩種case都找不到就繼續循環，實際寫代碼時，一般會設置個循環次數限制，避免這里耗費太多時間，或者進入死循環，畢竟我們上面只是證明存在滿足條件的步長，不代表我們可以在有限步驟內能找到這個滿足條件的步長。
Zoom Stage
Zoom Stage會不斷縮小搜索區間\((\alpha_{lo},\alpha_{hi})\)，直到找到滿足Strong Wolfe Conditions的步長。
用定理1一樣的證明方法，可以證明如果\(\alpha_{max}\)選的足夠大，則在調用\(Zoom\)之前，初始化搜索區間\((0,\alpha_{max})\)和后續搜索區間\((\alpha_i,\alpha_{max})\)都存在滿足Strong Wolfe Conditions的步長。
下面證明
定理2: Zoom Stage初始搜索區間\((\alpha_{lo},\alpha_{hi})\)和后續迭代更新后的區間都存在滿足Strong Wolfe Conditions的步長
顯而易見兩個步驟的區間都是隨着各自循環逐漸減小的，如果能夠證明定理2，則起碼說明我們沒有盲目在根本不可能存在滿足Strong Wolfe Conditions的區間搜索，並且我們搜索區間逐漸減小，離目標也越來越近了。並且假設統一記上面區間為\((a,b)\)，注意這里\(a\)可能大於\(b\)，也可能小於\(b\),則可以證明\(min(h(a),h(b))\)在兩個步驟循環里面是單調遞減的，即每迭代一次找到一個更小值，這樣即使最后沒有找到滿足Strong Wolfe Conditions的步長，我們也可以返回搜索到的最小值對應步長。
容易證明Zoom函數第一次傳入的搜索區間和后續更新后的搜索區間都同時滿足下面兩個條件(這里不證了)：
（1）\(h'(a)(b-a)<0\) 且\(a\)滿足sufficient decrease condition，即\(F(a)<=0\)
（2）如果\(b\)滿足sufficient decrease condition則\(h(a)<h(b)\)
當搜索區間同時滿足上面兩個條件時，那么存在\(\alpha\in(a,b)\)滿足Strong Wolfe Conditions條件。分兩種情況來證明
case 1:
假設\(b\)滿足sufficient decrease condition， 此時有\(h(a)<h(b)\)
由\(h(a)<h(b)\)知存在\(c\in{(a,b)}\)使得\(h'(c)(b-a)>0\)，聯合\(h'(a)(b-a)<0\) 得到存在\(d\in{(a,c)}\)使\(h'(d)=0\)， 進而知道沿着點\(a\)像\(d\)的方向走，能找到\(h'(\alpha)=0\)的點，記第一次遇見的\(h'(\alpha)=0\)的點為\(\alpha\)，則\(\alpha\)滿足curvature condition。又由\(h'(a)(b-a)<0\) 知，從\(a\)走向\(\alpha\)的過程，\(h(x)\)函數值是單調遞減的，所以\(h(\alpha)<h(a)<h(b)\)， 在聯合\(a,b\)都滿足sufficient decrease condition能推出\(\alpha\)也滿足sufficient decrease condition。綜合以上得出\(\alpha\)滿足sufficient decrease condition。
case 2:
假設\(b\)不滿足sufficient decrease condition，那么由\(F(b)>0\)和\(F(a)<=0\)知道存在\(c\in{(a,b)}\)使得\(F'(c)(b-a)>0\)，又由\(h'(a)(b-a)<0\) 和\(a\)不滿足curvature condition能推出\(F'(a)(b-a)<0\)，由導數連續性知道從\(a\)朝着\(c\)的方向走，能找到第一個\(F'(\alpha)=0\)的點，並且這個走的方向，\(F(x)\)值是單調減少的，所以\(F(\alpha)<F(a)\)，即\(\alpha\)滿足sufficient decrease condition條件，並且由\(F'(\alpha)=0\)得出\(\alpha\)滿足curvature condition，從而得證。
Convergence Of Line Search
前面一小節已經提到可以證明，Line Search算法搜索到的函數值是單調遞減的，但函數值單調遞減不代表\(x_k\)能夠收斂到駐點（如果是凸函數就是全局最小點），單調遞減也可能收斂到其他點，或者如果函數無下界，則\(x_k\)也可能發散。
如果目標函數\(f(x)\)滿足Lipschitz continuity，記\(cos^2(\theta_k)=\dfrac{-f'(x_k)p_k}{|f(x_k)|\ |p_k|}\)即負搜索方向與導數夾角，則容易證明(證明可以參考Numerical optimization[2:1])
上式暗含着 \(cos(\theta_k)||f'(x_k)||\rightarrow0\)，如果每次能保證每次搜索方向與導數夾角余弦大於某一個正整數（Steepest Descent就是每次夾角余弦為1），則 \(||f'(x_k)||\)收斂到0，即證明了在目標函數滿足Lipschitz continuity且搜索方向滿足一定條件下，Line Search算法會收斂到一個駐點，如果是凸函數，就收斂到全局最小點。
Gradient Gased Method
基於梯度的優化算法，搜索方向形式為
其中 \(B\)是正定矩陣，這樣就保證了搜索方向是函數值下降方向
Steepest Descent方法取 \(B=I\)
Newton法取 \(B\)為Hessian矩陣的逆
Quasi-Newton方法根據最近n次搜索的函數信息來近似一個Hessian矩陣的逆作為 \(B\)
下面將詳細分析各個方法，以及在Line Search方法下的收斂速度
記 \(x^*\)是搜索方向前方駐點
Steepest Descent
令搜索方向\(p_k=-If'(x_k)=-f'(x_k)\)就得到Steepest Descent算法，正如其名，\(-f'(x_k)\)方向是函數局部極小領域內函數值下降最快的方向，假設每次Line Search搜索到的\(\alpha\)都是局部最優值即Exact Line Search，則\(h'(\alpha)=0\)，即
，可以看出相鄰搜索方向相互垂直。
如果目標函數的二階導數Hessian矩陣滿足$lI\preceq\ G(x)\preceq LI \(，\)l,L \(分別是Hessian矩陣\)G \(的最小最大特征值，那么可以證明在Exact Line Search下<center>\)\left | x_{k+1}-x^* \right |_Q^2<=\left(\dfrac{L-l}{L+l}\right)^2\left | x_k-x^* \right |_Q^2$
也即可以證明Line Search線性收斂，這個具體怎么證明我還沒精力去看，只看了 pluskid [8]在步長固定的情況下對線性收斂的證明。
Newton Method
記\(G\)為目標函數Hessian矩陣，令\(p_k=-G^{-1}f'(x_k)\)，則得到Newton Method。顯然如果\(G\)是正定矩陣，則\(f'(x_k)p_k<0\)，即\(p_k\)是函數下降方向，記\(x^s\)為函數二階泰勒近似的全局最小值，則\(p_k=x^s-x_k\)，即搜索方向指向\(x^s\)，所以對於二次函數一次搜索就找到全局最小值。
假設初始點\(x_0\)離\(x^*\)足夠近，並且\(G\)滿足Lipschitz continuity，每次步長固定為1，可以證明搜索到的點序列\({x_k}\)收斂到\(x^*\)， 並且收斂速度是二次的，證明詳見Numerical optimization[2:2]。
Quasi-Newton Methods
Newton Method要計算和存儲Hessian矩陣，並求Hessian矩陣的逆，所以計算和存儲代價比較大，Quasi-Newton Methods通過近似Hessian矩陣的逆來減少計算量
DFP
DFP[3:1]算法通過最近\(k\)次迭代的導數及步長信息來近似計算一個Hessian矩陣的逆\(H=G^{-1}\)，通過迭代增量更新減小計算代價。迭代更新公式如下：
\begin{equation}
H_{k+1}=H_k-\dfrac{H_ky_ky_k^TH_k}{y_k^TH_ky_k}+\dfrac{\delta_k\delta_k^T}{y_k^T\delta_k}
\end{equation}
其中\(\delta_k=x_{k+1}-x_k\)，\(y_k=g_{k+1}-g_k\)，\(g_k=f'(x_k)\)，上面公式可以參考文獻[3:2]看看怎么推出的。也可以參考Numerical optimization[2:3]的推導方法，即先通過帶等式約束的最小化求出\(G\)的近似\(B\)的迭代更新公式，然后通過Sherman–Morrison formula變化得到\(H\)的迭代公式。可以采用兩次應用rank-one形式變化公式得到\(H\)（可以參考邏輯回歸模型和算法文中的變換方法）[9]，也可以采取一次rank-two形式變換公式（可以參考Quasi-Newton methods for minimization里面的變化方法）。
推導都要求\(H_{k+1}\)滿足\(H_{k+1}y_k=\delta_k\)且\(H_{k+1}=H_{k+1}^T\)，因為\(G\delta_k=y_k\)且，\(G^T=G\)，\(H_{k+1}\)想要近似\(G^{-1}\)自然要滿足這些性質。
文獻[3:3]證明了在二次正定函數下即\(f(x)=f(x_0)+a^Tx+x^TGx\)，使用公式（9）生成的\(H\)，在Exact Line Search下，經過\(N\)次迭代后收斂到\(G^{-1}\)，\(N\)是數據的緯度，並且搜索到全局最小值。證明過程如下：
如果\(H_k\)是對稱正定矩陣，容易證明出\(H_{k+1}\)也是對稱正定矩陣，文獻[3:4]給出了遞推證明，但證明有點問題，正確證明應該是\(<p|p><q|q>-<p|q>^2\) 和\(<x|\delta_k>^2\)不能同時為0，因為如果\(<p|p><q|q>-<p|q>^2=0\)，則\(x=\lambda{y_k}\)，進而\(<x|\delta_k>=<\lambda{y_k}|\delta_k>=<\delta_k|G|\delta_k>\neq0\)，所以得證。
上面是在\(f(x)\)是正定函數和Exact Line Search下證明出\(H_{k+1}\)保持正定性，其實對於\(f(x)\)是任意普通凸函數的情況，也可以證明在Line Search下（不要求Exact），\(H_{k+1}\)保持正定性，從上面推導可以看出只要\(y_k^T\delta_k>0\)即可推出\(H_{k+1}\)保持正定性，由sufficient decrease condition知\(f'(x_{k+1})^Tp_k>c_2f'(x_k)^Tp_k\)，由\(\delta_k=\alpha_kp_k\)知\(f'(x_{k+1})^T\delta_k>c_2f'(x_k)^T\delta_k>f'(x_k)^T\delta_k\)，進而推出\((f'(x_{k+1})-f'(x_k))^T\delta_k>0\)即\(y_k^T\delta_k>0\)。
\(H_k\)正定保證了搜索方向下降。
如果能夠證明\(H^NG\delta_i=\delta_i，0\leq{i} <N\)，則易知\(H^N=G^{-1}\)，即\(H\)經過\(N\)次迭代收斂到\(G^{-1}\)
接下來證明下面一組公式
\begin{equation}
\delta_i^TG\delta_j=0，0\leq{i} < j <k
\end{equation}
\begin{equation}
H_kG\delta_i=\delta_i，0\leq{i} <k
\end{equation}
當\(k=1\)，顯然成立，下面遞推證明，假設已知\(k\)時，上面一組公式成立，則
\(g_k=a+Gx_k=a+G(x_{i+1}+\delta_{i+1}+...+\delta_{k-1})=g_{i+1}+G(\delta_{i+1}+...+\delta_{k-1})\)
\(\delta_i^Tg_k=\delta_i^Tg_{i+1}+\delta_i^TG(\delta_{i+1}+...+\delta_{k-1})=\delta_ig_{i+1}=0，0\leq{i} <k\)
聯合上面公式和公式(11)可以得到
\(0=\delta_i^Tg_k=-\alpha_k\delta_{i}^TGH_kg_k=-\alpha_k\delta_{i}^TG\delta_k，0\leq{i} <k\)
上式和\(k\)公式(10)合並起來即知\(k+1\)時，公式(10)仍然成立。
\(H_{k+1}G\delta_i=H_{k}G\delta_i-\dfrac{H_ky_ky_k^TH_kG\delta_i}{y_k^TH_ky_k}+\dfrac{\delta_k\delta_k^TG\delta_i}{y_k^T\delta_k}=H_{k}G\delta_i-\dfrac{H_ky_ky_k^TH_kG\delta_i}{y_k^TH_ky_k}+\dfrac{\delta_k\delta_k^TG\delta_i}{y_k^T\delta_k}，0\leq{i} <k+1\)
由\(k+1\)時的公式(10)和\(k\)時公式(11)可以推出
\(H_{k+1}G\delta_i=H_{k}G\delta_i-\dfrac{H_ky_ky_k^T\delta_i}{y_k^TH_ky_k}+\dfrac{\delta_k\delta_k^Ty_i}{y_k^T\delta_k}=H_{k}G\delta_i-\dfrac{H_ky_k(G\delta_k)^T\delta_i}{y_k^TH_ky_k}+\dfrac{\delta_k\delta_k^TG\delta_i}{y_k^T\delta_k}=H_{k}G\delta_i-0-0=\delta_i，0\leq{i} <k+1\)
即得證\(k+1\)時公式(11)成立。從而證明了在二次正定函數下，經過\(N\)次迭代后收斂到Hessian矩陣\(G\)的逆。
BFGS
BFGS與DFP方法類似，迭代公式推導方法也類似。
其中 \(s_k\)即為DFP中的 \(\delta_k\)，通過最小化上面式子即求出 \(H_k\)的迭代更新公式
據說特定條件下DFP和BFGS都能保證全局收斂性，並且都具有超線性收斂速度，但in practice BFGS比DFP更魯棒高效，所以使用更普遍。
L-BFGS
DFP和BFGS要計算和存儲\(H\)，存儲和迭代計算\(H\)的代價都是\(O(N^2)\)。在大型廣告和推薦系統中，\(N\)可能達到上億級別，這樣存儲和計算代價太高。根據\(H_k\)的遞推公式，將\(H_k\)做\(m\)次展開，得
\begin{equation}
\left{\begin{aligned}
H_k&=V_{k-1}^TV_{k-2}^T...V_{k-m}^TH_{k-m}V_{k-m}...V_{k-2}V_{k-1}\
&+V_{k-1}^TV_{k-2}^T...V_{k-m+1}^Ts_{k-m}\rho_{k-m}s_{k-m}^TV_{k-m+1}...V_{k-2}V_{k-1}\
&+V_{k-1}^TV_{k-2}^T...V_{k-m+2}^Ts_{k-m+1}\rho_{k-m+1}s_{k-m+1}^TV_{k-m+2}...V_{k-2}V_{k-1}\
&+...\
&+V_{k-1}^Ts_{k-2}\rho_{k-2}s_{k-2}^TV_{k-1}\
&+s_{k-1}\rho_{k-1}s_{k-1}^T
\end{aligned}
\right.
\end{equation}
其中\(V_k=I-\rho_{k}y_ks_k^T\)
如果取\(H_{k-m}\)為對角陣，則\(H_{k-m}\)也可以當做一個\(O(N)\)的向量存儲，或者可以將\(H_{k-m}\)分解為兩個向量的乘積形式。相當於將\(H_k\)全部用向量表示出來，\(V_k\)也是用向量表示的。用最近\(m\)次搜索的函數信息來近似\(H_k\)。這樣存儲代價變為\(O(Nm)\)。
我們最終目的是通過\(p_k=-H_kg_k\)求出迭代搜索方向。如果將\(g_k\)從右向左，右乘以\(H_k\)展開式，由於都是向量點擊運算，（因為上文已知\(H_k\)展開式里面每個元素都是向量，\(g_k\)右乘\(H_k\)遇見的每個向量或者是航向量，這是就是向量點擊；或者遇見列向量，這時已經乘出來的部分就是個常數)，所以計算代價是\(O(Nm^2)\)，觀察\(H_k\)展開式可以看出，里面存在可以復用的子運算結構。比如\(H_k\)展開式每行都是對稱結構，相鄰兩行間大部分運算也是一樣的。
記\(q_{k-m+1}=V_{k-m+2}...V_{k-2}V_{k-1}g_k\)表示\(H_k\)展開式某一行（其實是第三行）的右邊，則易知
\(\alpha_{k-m+1}=\rho_{k-m+1}s_{k-m+1}^Tq_{k-m+1}\)
\(q_{k-m}=q_{k-m+1}-\alpha_{k-m+1}y_{k-m+1}=V_{k-m+1}q_{k-m+1}（第二行右邊）\)
即\(\alpha,q\)存在遞推關系，從下往上，可以將\(H_kg_k\)展開式每一行右邊遞推出來，存到臨時變量\(\alpha,q\)。
觀察\(H_k\)展開式左邊，可以發現前\(r\)行的左邊前\(m-r+1\)項是完全一樣的，即可以將前\(r\)行的左邊部分提取出作為公共乘子。另一方面我們需要將所有行累加起來，所以可以采取從上到下迭代累加方式。記\(r_r\)表示前\(r\)行扣除左邊最長公共乘子的累加和，\(r_1=H_{k-m}q_{k-m-1}=H_{k-m}V_{k-m}...V_{k-2}V_{k-1}\)表示前1行扣除公共乘子的累加和，則
記 \(\beta_{r+1}=\rho_{k-m+r-1}y_{k-m+r-1}^Tr_r\)，則
易知 \(r_{m+1}=H_kg_k\)，且上面總的 計算時間復雜度是 \(O(Nm)\)。將上述計算過程描述出來即如下圖所示（圖片來自Numerical optimization [2:4]， \(\bigtriangledown\ f_k\)即 \(g_k\)）
IIS
IIS即 Improved Iterative Scaling Algorithm的思想如下，對於求解凹函數\(f(x)\)的最大值問題，如果可以找到一個新的更簡單的凹函數\(g(x)\)使得\begin{equation}
\begin{cases}
g(x)&\geq f(x)\
g(x_k)&=f(x_k)\
g'(x_k)&=f'(x_k)
\end{cases}
\end{equation}
那么對\(f(x)\)的迭代求最小值可以轉化為對\(g(x)\)迭代求最小值（\(g(x)\)可能需要在每次迭代后重新構造）。因為
case 1：如果\(g(x)\)能找到一個\(g(x_{k+1})>g(x_k)\)，那么由公式（13）知道
即 \(g(x)\)的下一個搜索到的更大值位置 \(x_{k+1}\)也使 \(f(x)\)的值增大。
case 2：如果 \(g(x)\)無法找到一個比 \(g(x_k)\)更大的值，那么由 \(g(x)\)是凹函數，可知 \(g'(x_k)=0\)，又由公式（13）知 \(f'(x_k)=0\)所以 \(f(x_k)\)也是 \(f(x)\)的最大值，求解可以結束了。
問題主要是怎么構造滿足公式（13）的凹函數 \(g(x)\)，使 \(g(x)\)或者 \(g'(x)\)足夠簡單方便求解。比如 \(x\)是多維向量，如果可以找到一個新的 \(g(x)\)，使 \(g(x)\)可以分解成各個獨立無關的部分，每個部分只含其中一個維度的變量。或者使 \(\dfrac{\partial g(x)}{\partial x_j}\)只含單一維度變量。那找 \(g(x)\)的更小值就簡單很多。
對於類似CRF這種形式的目標函數 \(\mathcal{L}(\lambda)\)（參見公式（3）），文獻 [5:1]闡述了一種構造求解方法。
記 \(\lambda^k\)是當前搜索到的最優點， \(f(\delta)=\mathcal{L}(\lambda^k+\delta)-\mathcal{L}(\lambda^k)\)，則最大化 \(f(\delta)\)與最大化 \(\mathcal{L}(\lambda)\)是等價的。 \(\lambda^*=\lambda^k+\delta^*\)，其中 \(\lambda^k,\delta^*\)分別是 \(\mathcal{L}(\lambda),f(\delta)\)的最大值。
下面證明一個不等式
記 \(LogF(x)=-log(x)-1+x\)，則 \(LogF(1)=0\)且 \(LogF'(x)=1-\dfrac{1}{x}<=0,\ x\in{(0,1)}\), 由 \(LogF(x)\)在 \((0,1)\)區間單調下降，且 \(LogF(1)=0\)，推出 \(LogF(x)>0,\ x\in{(0,1)}\)，上面不等式得證。
利用上面不等式得
\begin{equation}
\begin{cases}
f(\delta)&\geq \sum_m\left(\sum_j\delta_jF_j\left(x^m,y^m\right)+1-\dfrac{Z_{\lambda^k+\delta}\left(x^m\right)}{Z_{\lambda^k}\left(x^m\right)}\right)\
&=\sum_m\left(\sum_j\delta_jF_j\left(x^m,y^m\right)+1-\sum_yp_{\lambda^k(y|x^m)}exp\left(\sum_j\delta_jF_j(x^m,y)\right)\right)
\end{cases}
\end{equation}
記 \(F^\#(x,y)=\sum_jF_j(x,y)\)，知 \(\sum_j\dfrac{F_j(x^m,y)}{F^\#(x^m,y)}=1\)，應用 Jensen's inequality，得
\begin{equation}
\begin{cases}
f(\delta)&\geq \sum_m\left(\sum_j\delta_jF_j\left(x^m,y^m\right)+1-\sum_yp_{\lambda^k(y|x^m)}exp\left(\sum_j\delta_jF^#(x^m,y)\dfrac{F_j(x^m,y)}{F^#(x^m,y)}\right)\right)\
&\geq \sum_m\left(\sum_j\delta_jF_j\left(x^m,y^m\right)+1-\sum_yp_{\lambda^k(y|x^m)}\sum_j\left(\dfrac{F_j(x^m,y)}{F^#(x^m,y)}exp\left(\delta_jF^#(x^m,y)\right)\right)\right)
\end{cases}
\end{equation}
記上式右邊為 \(g(\delta)\)，則$ g(\delta) \(是凹函數且滿足公式（13）的3個條件。可以驗證\)g(0)=f(0) \(且\)g'(0)=f'(0) \(,有上面不等式知\)f(\delta)\geq\ g(\delta) \(。 現在只要找出一個\)g(\delta)>g(0) \(即可， \)\dfrac{\partial\ g(\delta)}{\partial\delta_j}= \sum_m\left(\sum_jF_j\left(x^m,y^m\right)+1-\sum_yp_{\lambda^k(y|x^m)}F_j(x^m,y)exp\left(\delta_jF^#(x^m,y)\right)\right)$
導數為 \(0\)的點即為 \(g(\delta)\)最大值的點，上面梯度只含單一維度變量，令 \(\dfrac{\partial\ g(\delta)}{\partial\delta_j}=0\)可以直接用牛頓法求出各個 \(\delta_j\)，並且求解速度很快，二次收斂速度。
代碼分析
CRF代碼分析
關鍵函數執行流程如下，附注釋
Encoder::learn(const char *templfile,//模板文件，用於生成特征，包括Unigram（狀態特征）和Bigram（轉移特征）
const char *trainfile,//訓練文件，里面有(x,y)標注樣本
size_t maxitr,//最大迭代次數，這里相當於求函數梯度次數，當迭代次數超過maxitr時就結束繼續優化
size_t freq,//特征最小出現頻次，小於此值則過濾掉
double eta, //連續3次相鄰兩次迭代似然函數值之差小於eta的話，就結束繼續優化
double C,//L2正則項權重縮放因子
unsigned short thread_num//求梯度並發線程數
)
{
EncoderFeatureIndex feature_index;
feature_index.open(templfile, trainfile))//讀入模板文件， 根據訓練文件提取標注集合，即看看y可能取哪些標注
std::ifstream ifs(trainfile);
while (ifs) {
TaggerImpl *_x = new TaggerImpl();
_x->open(&feature_index, &allocator);
if (!_x->read(&ifs) || !_x->shrink()) {//讀入訓練樣本， 每個TaggerImpl對應一個訓練樣本；先TaggerImpl::read讀入一個樣本及對應的標注序列，
//並調用TaggerImpl::shrink()生成特征，分配特征id（特征字符串一樣，則id一樣，同時記錄每個特征出現的頻次）；每個樣本每個位置的特征id集合打包成一個vector保存在 FeatureCache中；這里一個id是只根據x，未用y（標注序列）生成的，所以每個id如果是Unigram生成的，實際上對應L個連續的id，L是標注集合大小，即標注取值個數，每個id如果是Bigram生成的則對應L*L個連續的id，即當相鄰的y_(i-1),y_i可能取值個數。
}
_x->set_thread_id(line % thread_num);//這個樣本分配給第(line % thread_num)個線程處理
}
feature_index.shrink(freq, &allocator); //所有訓練數據讀完后，根據頻次統計，將低頻特征去掉， 並重新分配特征id，以使特征空間緊湊；即一些特征id丟掉后，特征id會出現不連續情況，所以重新分配id，以使max(id)以內的所有id都有用。
//因為后面優化過程中分配的空間跟max(id)有關系,FeatureIndex::size() { return maxid_; }
//樣本數據(x,y)都讀進來了，狀態和轉移特征都生成了，下面可以跑優化，求解每個特征對應的權重了
std::vector <double> alpha(feature_index.size()); //alpha即待求解的特征權重數組
runCRF(x, &feature_index, &alpha[0],
maxitr, C, eta, shrinking_size, thread_num, false);
}
bool runCRF(const std::vector<TaggerImpl* > &x,
EncoderFeatureIndex *feature_index,
double *alpha,
size_t maxitr,
float C,
double eta,
unsigned short thread_num) {
double old_obj = 1e+37;
int converge = 0;
LBFGS lbfgs;
std::vector<CRFEncoderThread> thread(thread_num);
//這里線程初始化，代碼省略
for (size_t itr = 0; itr < maxitr; ++itr) {
//這里啟動CRFEncoderThread線程，代碼省略,啟動后，CRFEncoderThread::run()(詳見下一個代碼注釋)得到執行，
//這里等待所有CRFEncoderThread線程計算完畢，多個線程正忙碌着跑forward-backward算法計算梯度thread[0].expected，也同時計算似然函數值thread[0].obj，並用當前特征權重跑viterbi算法預測每個樣本的標注序列，然后與樣本的真實標注序列比對，記錄預測錯誤數到thread[0].err和thread[0].zeroone
//這里所有線程已經計算完畢，下面將所有線程的結果合並
for (size_t i = 1; i < thread_num; ++i) {
thread[0].obj += thread[i].obj;//將函數值累加
thread[0].err += thread[i].err;//預測錯誤的的標注數，一個樣本的一個位置預測錯則計數一次
thread[0].zeroone += thread[i].zeroone;//預測錯誤的樣本（序列）數，整個樣本只要有一個標注預測錯則計數一次，多個標注錯也只計數一次
}
for (size_t i = 1; i < thread_num; ++i) {
for (size_t k = 0; k < feature_index->size(); ++k) {
thread[0].expected[k] += thread[i].expected[k];//所有線程梯度值累加起來
}
}
for (size_t k = 0; k < feature_index->size(); ++k) {
thread[0].obj += (alpha[k] * alpha[k] //(2.0 * C));//函數值加上L2項
thread[0].expected[k] += alpha[k] //梯度值加上L2項
}
//判斷相鄰兩次迭代的目標函數值之差
double diff = (itr == 0 ? 1.0 : std::abs(old_obj -
thread[0].obj)/old_obj);
old_obj = thread[0].obj;
converge += (diff < eta) ? 1 : 0;
//如果連續3次目標函數值diff小於閥值，或者達到最大迭代次數，則停止迭代
if (itr > maxitr || converge == 3) {
break;
}
//lbfgs里面iflag=0（gnorm / xnorm <= eps 時）會導致這里返回值<=0，被認為失敗， 感覺應該認為是已經足夠接近最優值才對。
if (lbfgs.optimize(feature_index->size(),
&alpha[0],
thread[0].obj,
&thread[0].expected[0], orthant, C) <= 0) {
return false;
}
}//itr < maxitr
return true;
}
void CRFEncoderThread::run() {
obj = 0.0; err = zeroone = 0;
std::fill(expected.begin(), expected.end(), 0.0);//梯度初始化為0
for (size_t i = start_i; i < size; i += thread_num) {
obj += x[i]->gradient(&expected[0]);//obj為似然函數值，x即樣本，對每個樣本調用TaggerImpl::gradient(double *expected)算梯度，並跑viterbi算法預測標注序列
int error_num = x[i]->eval();//評估預測值與真實標注的誤查
err += error_num;//預測錯誤的的標注數，一個樣本的一個位置預測錯則計數一次
if (error_num) {
++zeroone;//預測錯誤的樣本（序列）數，整個樣本只要有一個標注預測錯則計數一次，多個標注錯也只計數一次
}
}
}
接下來就是比較重要的計算梯度的代碼，這里可以對照公式（8）來看，進入這個函數已經確定樣本了，即最外面的累加m下標已經確定。
先看公式的前半部分，即經驗分布期望。對於狀態特征，node_[i][answer_[i]]等於\(y_i^m\)即第\(m\)個樣本在第\(i\)個位置的真實標注值，node_[i][answer_[i]]->fvector存放着這個位置標注確定后，對應的所有狀態特征id，因為一個位置可能有多個Unigram模板提取出來的不同特征。轉移特征經驗分布類似。
再看公式的后半部分，即模型分布期望。下面代碼的i對應着公式里面的i，即序列下標。node_[i][j]對應着l(狀態特征)，Node::calcExpectation里面的lpath對應着l,l'轉移特征。公式（8）只是對其中一個特征做計算，這里要對所有特征計算，因為一旦\(l\)確定后，這個序列位置\(i\)的所有特征的\(\dfrac{1}{Z(x^m)}\alpha_{n}(a)\beta_{n}(a)\)值都是一樣的，所以把這個模型概率加到所有特征id即可（因為特征函數返回0或1）。轉移特征模型分布也類似。
TaggerImpl::gradient()函數返回似然函數值Z_ - s ，可以對照公式（3）的來看怎么計算的。
上面的梯度和函數值計算時是把公式前面加個負號，將最大化問題轉化為最小化問題。
buildLattice()可以參考hankcs的CRF++代碼分析里面的圖
注意計算前向后向概率時為避免浮點溢出，采用將值取對數保存方式，即logsumexp方式，詳見hankcs的計算指數函數的和的對數
double TaggerImpl::gradient(double *expected) {
buildLattice();//構建Node 和 Path ； 並計算Node（對於狀態特征）和Path權重（對應轉移特征）
forwardbackward();//根據公式（7）遞推計算前向概率alpha和后向概率beta
double s = 0.0;
for (size_t i = 0; i < x_.size(); ++i)
for (size_t j = 0; j < ysize_; ++j)
node_[i][j]->calcExpectation(expected, Z_, ysize_); //計算梯度中的狀態特征模型分布期望，和轉移特征模型分布期望。
for (size_t i = 0; i < x_.size(); ++i) {
for (const int *f = node_[i][answer_[i]]->fvector; *f != -1; ++f) {
--expected[*f + answer_[i]];//梯度中的狀態特征經驗分布期望，對應狀態特征，一個f實際上對應L個特征，即f是根據樣本x值生成的，未考慮y，對應不同的y就有不同的特征，f只是這一組特征的起始值，所以這里加上answer_[i]，即實際樣本真實標注的特征。
}
s += node_[i][answer_[i]]->cost; // UNIGRAM cost公式（3）似然函數值里面的狀態特征權重
const std::vector<Path *> &lpath = node_[i][answer_[i]]->lpath;
for (const_Path_iterator it = lpath.begin(); it != lpath.end(); ++it) {
if ((*it)->lnode->y == answer_[(*it)->lnode->x]) {//判斷當前path是否是樣本真實標注值對。
for (const int *f = (*it)->fvector; *f != -1; ++f) {
--expected[*f +(*it)->lnode->y * ysize_ +(*it)->rnode->y];//梯度中的轉移特征經驗分布期望，這里加(*it)->lnode->y * ysize_ +(*it)->rnode->y的含義與狀態特征類似，只不過這里一個f對應L^2個特征。
}
s += (*it)->cost; // BIGRAM COST公式（3）似然函數值里面的轉移特征權重
break;
}
}
}
viterbi(); // call for eval() 跑viterbi算法，用當前特征權重預測標注序列y，后面會那預測值於真實標注值比較，以評估當前特征權重的模型好壞。
return Z_ - s ;//返回似然函數值，即公式（3）的值，取負即可
}
Node::calcAlpha()和Node::calcBeta()可以參考公式（7）遞推公式，
注意Node::calcBeta()與公式（7）有點不同，最后多加了一項beta += cost; 等於比公式（7）中的值多了第一個節點的狀態特征權重。所以后續使用時也稍微變化了下
void Node::calcAlpha() {
alpha = 0.0;
for (const_Path_iterator it = lpath.begin(); it != lpath.end(); ++it) {
alpha = logsumexp(alpha,
(*it)->cost +(*it)->lnode->alpha,
(it == lpath.begin()));
}
alpha += cost;//前向后向概率保存的是取對數后的值，防止保存原始值時中間計算溢出； log(....)
}
void Node::calcBeta() {
beta = 0.0;
for (const_Path_iterator it = rpath.begin(); it != rpath.end(); ++it) {
beta = logsumexp(beta,
(*it)->cost +(*it)->rnode->beta,
(it == rpath.begin()));
}
beta += cost;//這里比公式（7）多加了第一個節點的狀態特征權重
}
void Node::calcExpectation(double *expected, double Z, size_t size) const {
const double c = std::exp(alpha + beta - cost - Z);//上面beta多加了一個cost，所以這里減掉
for (const int *f = fvector; *f != -1; ++f) {
expected[*f + y] += c; //狀態特征期望
}
for (const_Path_iterator it = lpath.begin(); it != lpath.end(); ++it) {
(*it)->calcExpectation(expected, Z, size);//計算轉移特征期望
}
}
L-BFGS代碼分析
void LBFGS::lbfgs_optimize(int size, //數組x和g的長度
int msize,//即lbfgs里面的m，使用最近m次迭代的函數信息
double *x,//當前x值，對於crf++就是特征權重數組
double f,//當前函數值，對於crf++就是似然函數值
const double *g,//當前梯度
double *diag,//存放H(-m)即初始對角陣
double *w,//存放最近m次迭代的ρ，y，s，以及計算搜索方向時用來存放臨時變量α，q，r，等
bool orthant, //是否采取L1正則化
double C, //正則化權重因子
double *v,//對於L2正則化v=g，對於L1正則化，v是偽梯度
double *xi,//用於L1正則化時存放x的符號，用來限定Line Search在一個象限里面
int *iflag//見下文解釋
) {
double yy = 0.0;
double ys = 0.0;
int bound = 0;
int cp = 0;
if (orthant) {
--xi;
pseudo_gradient(size, v, x, g, C);//計算偽梯度
}
if (*iflag == 1) goto L172;//*iflag為1意味着上次執行了Line search(mcsrch)，並且未找到滿足Wolfe condition的step，但計算出一個新的step，並同時更新了x數組，返回到CRF主程序，讓CRF計算step步長時的函數值和梯度值，執行到這里時，說明CRF計算完畢，重新進來直接跳到L172，讓mcsrch評估當前步長是否滿足Wolfe condition
if (*iflag == 2) goto L100;//這個分支代碼里未用，感覺是由用戶自己生成diag（H(-m)）來用
//執行到這里意味着iflag_ = 0，iflag_=0只會出現一次，只有第一次調用lbfgs時才為0
//w_.resize(size * (2 * msize + 1) + 2 * msize);
//w[0]未用
//w[ (1, size) ] 存放lbfgs迭代時的臨時變量q或r; 在計算出line search搜索方向后，w[ (1, size) ]也用來存放梯度f'(x_k)，以便后面確定好步長后，來計算梯度變化（w[iypt + npt + i] = g[i] - w[i];）注意每次傳進來的g可能不是f'(x_k)，也有可能是f'(x_k+step*p_k)
//w[(size + 1, size + msize)] 存放ρ
//w[(size + msize + 1, size + 2*msize)] 存放α
//w[(size + 2*msize + 1, size + 2*msize+msize*size)] 存放s; 初始化s存儲區域基址；ispt = size + (msize << 1);
//w[(size + 2*msize + msize*size +1, size + 2*msize+2*msize*size)] 存放y ;初始化y存儲區基址；iypt = ispt + size * msize;
//第一次調用lbfgs，做初始化
if (*iflag == 0) {
point = 0;
for (int i = 1; i <= size; ++i) {
diag[i] = 1.0;//H（-m）初始化為單位矩陣
}
ispt = size + (msize << 1);//s下標基址
iypt = ispt + size * msize;//y下標基址
for (int i = 1; i <= size; ++i) {
w[ispt + i] = -v[i] * diag[i];//初始化第一次迭代搜索方向
}
stp1 = 1.0 / std::sqrt(ddot_(size, &v[1], &v[1]));
}
// MAIN ITERATION LOOP
while (true) {
//執行到這里，或者是從第一次調用lbfgs進來的（iflag=0），或者就是從while循環下面mcsrch之后跳轉過來的，mcsrch找到滿足Wolfe condition的step后，會繼續下一輪循環，跑到這里以重新用lbfgs計算新的迭代方向。所以也意味着運行到這里后，下次傳給mcsrch的是新的搜索方向。
++iter;
info = 0;//每次計算出新的搜索方向，info要置為1，這樣mcsrch知道傳進來的參數是f(x_k),f'(x_k),mcsrch會把這個值保存下來，后續跟其他f(x_k+step*p_k),f'(x_k+step*p_k)比較，判斷是否滿足Wolfe condition
if (orthant) {//L1正則化初始化xi，如前所述，執行到這里說明開啟了新一輪Line Search，即x_k已經找到，開始尋找x_(k+1)，所以初始化符號向量，以將后續Line Search限制在一個象限
for (int i = 1; i <= size; ++i) {
xi[i] = (x[i] != 0 ? sigma(x[i]) : sigma(-v[i]));
}
}
if (iter == 1) goto L165;//第一次迭代，沒有歷史s，y等搜索信息，所以直接用上面iflag=0計算的初始搜索方向
ys = ddot_(size, &w[iypt + npt + 1], &w[ispt + npt + 1]);
yy = ddot_(size, &w[iypt + npt + 1], &w[iypt + npt + 1]);
//初始化diag=H(-m)
for (int i = 1; i <= size; ++i) {
diag[i] = ys / yy;
}
L100:
cp = point;
if (point == 0) cp = msize;
w[size + cp] = 1.0 / ys; //存放ρ
//初始化w[(1,size)] 為-g
for (int i = 1; i <= size; ++i) {
w[i] = -v[i];
}
bound = min(iter - 1, msize); //用最近bound次的迭代信息（ρ，y，s）來求近似-H(-m)*g
//下面步驟迭代計算-H(-m)*g，可以參考L-BFGS小節LBFGS two loop recursion那個圖
//迭代計算q和α
//α= ρ * s_t * q
//q= q - a*y
cp = point;//point指向最新的搜索信息s，y的下一個位置，即未來存放s，y的下標，所以下面cp先減一。
for (int i = 1; i <= bound; ++i) {
--cp;
if (cp == -1) cp = msize - 1;
double sq = ddot_(size, &w[ispt + cp * size + 1], &w[1]);
int inmc = size + msize + cp + 1;
iycn = iypt + cp * size;
w[inmc] = w[size + cp + 1] * sq;
double d = -w[inmc];
daxpy_(size, d, &w[iycn + 1], &w[1]);
}
//初始化r=-H(-m)*q
for (int i = 1; i <= size; ++i) {
w[i] = diag[i] * w[i];
}
//迭代更新r
//beta = ρ_i * y_t * r
//r = r + (alpha_i - beta) * s
for (int i = 1; i <= bound; ++i) {
double yr = ddot_(size, &w[iypt + cp * size + 1], &w[1]);
double beta = w[size + cp + 1] * yr;
int inmc = size + msize + cp + 1;
beta = w[inmc] - beta;
iscn = ispt + cp * size;
daxpy_(size, beta, &w[iscn + 1], &w[1]);
++cp;
if (cp == msize) cp = 0;
}
if (orthant) {//L1正則化，則校正搜索方向，以使搜索方向不與偽梯度相反
for (int i = 1; i <= size; ++i) {
w[i] = (sigma(w[i]) == sigma(-v[i]) ? w[i] : 0);
}
}
//將搜索方向r=-H*g存放到w[ispt + point * size]
for (int i = 1; i <= size; ++i) {
w[ispt + point * size + i] = w[i];
}
L165:
// OBTAIN THE ONE-DIMENSIONAL MINIMIZER OF THE FUNCTION
// BY USING THE LINE SEARCH ROUTINE MCSRCH
nfev = 0;//nfev表示mcsrch評估函數值和梯度次數，置為0，每評估一次mcsrch里面會加一，達到最大次數后就認為搜索失敗，整個算法求解就失敗了，CRF主程序也退出。
stp = 1.0;
if (iter == 1) {
stp = stp1;
}
for (int i = 1; i <= size; ++i) {
w[i] = g[i];//存儲梯度f'(x_k)，供后面找到滿足Wolfe condition步長后用來計算梯度變化y=f'(x_k+step*p_k)-f'(x_k)
}
L172:
//每次計算出新搜索方向后，第一次調用mcsrch時， info=0， 導數g存於 w[1，size]中，
//當后續不斷迭代mcsrch知直到找出符合Wolfe condition的stp，即當info=1時，這時便可以將導數變化g-w[1,size]存到w[iypt + npt]中，並設置
//ρ, 存於w[size + cp] = 1.0 / ys (詳見L100處代碼)
//在找到Wolfe condition步長前，iflag=1，會直接跳轉鰲L172. 或者達到最大maxfev = 20還未找到合適步長，就認為失敗
mcsrch_->mcsrch(size, &x[1], f, &v[1], &w[ispt + point * size + 1],
&stp, &info, &nfev, &diag[1]);
if (info == -1) {
if (orthant) {//限制一次Line Search搜索的一系列點都在指定象限
for (int i = 1; i <= size; ++i) {
x[i] = (sigma(x[i]) == sigma(xi[i]) ? x[i] : 0);
}
}
*iflag = 1; //未找到合適步長，就用插值方法計算出新的步長，並更新x數組，返回CRF程序，讓CRF計算新的x處的函數值和梯度
return;
}
if (info != 1) {//尋找失敗，可能是達到maxfev = 20，也可能是其他情況
std::cerr << "The line search routine mcsrch failed: error code:"
<< info << std::endl;
*iflag = -1;
return;
}
//找到了滿足Wolfe condition的stp，
npt = point * size;
for (int i = 1; i <= size; ++i) {
w[ispt + npt + i] = stp * w[ispt + npt + i];//計算存儲s=stp*p_k=stp*H(-m)*g
w[iypt + npt + i] = g[i] - w[i];//計算存儲y=f'(x_k+step*p_k)-f'(x_k)
}
++point;
if (point == msize) point = 0;//保留最近m次搜索的s，y等信息,循環復用數組
double gnorm = std::sqrt(ddot_(size, &v[1], &v[1]));
double xnorm = max(1.0, std::sqrt(ddot_(size, &x[1], &x[1])));
if (gnorm / xnorm <= eps) {//梯度的模足夠小就停止繼續
*iflag = 0; // OK terminated
return;
}
//繼續循環，用l-bfgs求下一次的搜索方向。
}
return;
}
[^Line Search]:Moré J J, Thuente D J. Line search algorithms with guaranteed sufficient decrease[J]. ACM Transactions on Mathematical Software (TOMS), 1994, 20(3): 286-307.
Lafferty, J., McCallum, A., Pereira, F. Conditional random fields: Probabilistic models for segmenting and labeling sequence data. Proc. 18th International Conf. on Machine Learning. Morgan Kaufmann: 282–289. 2001. ↩︎
Gilbert, J. C., & Lemarechal, C. (2006). Numerical optimization: theoretical and practical aspects. Springer. ↩︎ ↩︎ ↩︎ ↩︎ ↩︎
R. Fletcher and M. J. D. Powell, A Rapidly Convergent Descent Method for Minimization , Computer Journal 6 (1963), 163–168. ↩︎ ↩︎ ↩︎ ↩︎ ↩︎
J.NOCEDAL, Updating quasi-Newton matrices with limited storage,MathematicsofComputation, 35 (1980), pp. 773–782. ↩︎
Berger A. The improved iterative scaling algorithm: A gentle introduction[J]. Unpublished manuscript, 1997. ↩︎ ↩︎
統計學習方法，李航 ↩︎
最大熵模型http://blog.csdn.net/zgwhugr0216/article/details/53240577 ↩︎
pluskid. Gradient Descent, Wolfe's Condition and Logistic Regression http://freemind.pluskid.org/machine-learning/gradient-descent-wolfe-s-condition-and-logistic-regression/ ↩︎
邏輯回歸模型和算法http://blog.csdn.net/zgwhugr0216/article/details/50145075 ↩︎