上篇的CRF++源碼閱讀中， 我們看到CRF++如何處理樣本以及如何構造特征。本篇文章將繼續探討CRF++的源碼，並且本篇文章將是整個系列的重點，會介紹條件隨機場中如何構造無向圖、前向后向算法、如何計算條件概率、如何計算特征函數的期望以及如何求似然函數的梯度。本篇將結合條件隨機場公式推導和CRF++源碼實現來講解以上問題。
開啟多線程
我們接着上一篇encoder.cpp文件中的learn函數繼續看，該函數的下半部分將會調用具體的學習算法做訓練。目前CRF++支持兩種訓練算法，一種是擬牛頓算法中的LBFGS算法，另一種是MIRA算法， 本篇文章主要探討LBFGS算法的實現過程。在learn函數中，訓練算法的入口代碼如下：
switch (algorithm) { case MIRA: //MIRA算法的入口 if (!runMIRA(x, &feature_index, &alpha[0], maxitr, C, eta, shrinking_size, thread_num)) { WHAT_ERROR("MIRA execute error"); } break; case CRF_L2: //LBFGS-L2正則化的入口函數 if (!runCRF(x, &feature_index, &alpha[0], maxitr, C, eta, shrinking_size, thread_num, false)) { WHAT_ERROR("CRF_L2 execute error"); } break; case CRF_L1: //LBFGS-L1正則化的入口函數 if (!runCRF(x, &feature_index, &alpha[0], maxitr, C, eta, shrinking_size, thread_num, true)) { WHAT_ERROR("CRF_L1 execute error"); } break; }
runCRF函數中會初始化CRFEncoderThread數組，並啟動每個線程，源碼如下：
bool runCRF(const std::vector<TaggerImpl* > &x, EncoderFeatureIndex *feature_index, double *alpha, size_t maxitr, float C, double eta, unsigned short shrinking_size, unsigned short thread_num, bool orthant) { ... //省略代碼
for (size_t itr = 0; itr < maxitr; ++itr) { //開始迭代， 最大迭代次數為maxitr，即命令行參數-m for (size_t i = 0; i < thread_num; ++i) { thread[i].start(); //啟動每個線程，start函數中會調用CRFEncoderThread類中的run函數 } for (size_t i = 0; i < thread_num; ++i) { thread[i].join(); //等待所有線程結束 }
... //省略代碼
CRFEncoderThread類中的run函數調用gradient函數，完成一系列的核心計算。源碼如下：
void run() { obj = 0.0; err = zeroone = 0; std::fill(expected.begin(), expected.end(), 0.0); //excepted變量存放期望 for (size_t i = start_i; i < size; i += thread_num) {//每個線程並行處理多個句子， 並且每個線程處理的句子不相同， size是句子的個數 obj += x[i]->gradient(&expected[0]); //x[i]是TaggerImpl對象，代表一個句子， gradient函數主要功能： 1. 構建無向圖 2. 調用前向后向算法 3. 計算期望 int error_num = x[i]->eval(); err += error_num; if (error_num) { ++zeroone; } } }
構造無向圖
我們知道條件隨機場是概率圖模型，幾乎所有的概率計算都是在無向圖上進行的。那么這個圖是如果構造的呢？答案就在gradient函數第一個調用 —— buildLattice函數中。該函數完成2個核心功能，1. 構建無向圖 2. 計算節點以及邊上的代價，先看一下無向圖的構造過程：
void TaggerImpl::buildLattice() { if (x_.empty()) { return; } feature_index_->rebuildFeatures(this); //調用該方法初始化節點（Node）和邊（Path），並連接 ... //省略代碼 }
void FeatureIndex::rebuildFeatures(TaggerImpl *tagger) const { size_t fid = tagger->feature_id(); //取出當前句子的feature_id，上篇介紹構造特征的時候，在buildFeatures函數中會set feature_id const size_t thread_id = tagger->thread_id(); Allocator *allocator = tagger->allocator(); allocator->clear_freelist(thread_id); FeatureCache *feature_cache = allocator->feature_cache();
//每個詞以及對應的所有可能的label，構造節點 for (size_t cur = 0; cur < tagger->size(); ++cur) { //遍歷每個詞， const int *f = (*feature_cache)[fid++]; //取出每個詞的特征列表，詞的特征列表對應特征模板里的Unigram特征 for (size_t i = 0; i < y_.size(); ++i) { //每個詞都對應不同的label， 每個label用數組的下標表示，每個特征+當前的label就是特征函數 Node *n = allocator->newNode(thread_id); //初始化新的節點，即Node對象 n->clear(); n->x = cur; //當前詞 n->y = i; //當前詞的label n->fvector = f; //特征列表 tagger->set_node(n, cur, i); //有一個二維數組node_存放每個節點 } }
//從第二個詞開始構造節點之間的邊，兩個詞之間有y_.size()*y_.size()條邊 for (size_t cur = 1; cur < tagger->size(); ++cur) { const int *f = (*feature_cache)[fid++]; //取出每個邊的特征列表，邊的特征列表對應特征模板里的Bigram特征 for (size_t j = 0; j < y_.size(); ++j) {//前一個詞的label有y_.size()種情況，即y_.size()個節點 for (size_t i = 0; i < y_.size(); ++i) {//當前詞label也有y_.size()種情況，即y_.size()個節點 Path *p = allocator->newPath(thread_id);//初始化新的節點，即Path對象 p->clear();
//add函數會設置當前邊的左右節點，同時會把當前邊加入到左右節點的邊集合中 p->add(tagger->node(cur - 1, j), //前一個節點 tagger->node(cur, i)); //當前節點 p->fvector = f; } } } }
圖構造完成后， 接下來看看節點和邊上的代價是如何計算的。那么代價是什么？我的理解就是特征函數值乘以特征的權重。這部分源碼在buildLattice函數中，具體如下：
for (size_t i = 0; i < x_.size(); ++i) { for (size_t j = 0; j < ysize_; ++j) { feature_index_->calcCost(node_[i][j]); //計算節點的代價 const std::vector<Path *> &lpath = node_[i][j]->lpath; for (const_Path_iterator it = lpath.begin(); it != lpath.end(); ++it) { feature_index_->calcCost(*it); //計算邊的代價 } } } //節點的代價計算函數 void FeatureIndex::calcCost(Node *n) const { n->cost = 0.0; #define ADD_COST(T, A) \ do { T c = 0; \ for (const int *f = n->fvector; *f != -1; ++f) { c += (A)[*f + n->y]; } \ //取每個特征以及當前節點的label，即為特征函數，且值為1，特征函數乘以權重（alpha_[*f + n->y]）是代價，特征函數為1所以代價=alpha_[*f + n->y]*1，對所有代價求和 n->cost =cost_factor_ *(T)c; } while (0) //cost_factor_是代價因子 if (alpha_float_) { ADD_COST(float, alpha_float_); } else { ADD_COST(double, alpha_); //將會在這里調用， 上一篇內容可以看到，CRF++初始化的是alpha_變量 } #undef ADD_COST }
//邊的代價計算函數與節點類似，不再贅述
看完源碼，我們舉個例子來可視化一下無向圖，仍然用上一篇中構造特征的那個例子。如果忘記了，出門左轉回顧一下。上一個例子中有三個詞，假設這三個詞分別是“我”、“愛”、“你”。構建的無向圖如圖一所示。
圖一
這個例子中，有三個詞和三個label，每個label用0,1，2表示，之前我們說過用數組下標代替label。每個詞有3個節點，且這三個節點的特征列表f是一樣的，由於label不一樣，所以他們的特征函數值不一樣。由於沒有bigram特征，所有邊上的特征列表都是f=[-1]。大部分資料的無向圖前后會加一個start節點和stop節點，加上后可以便於理解和公式推導。CRF++源碼中沒加，所以我們這里就沒有表示。在這里node_[0][0]對應就是最左上角的節點，代表“我”這個詞label為0的節點。我們再看一下node_[0][0]這個節點的代價如何計算的，node_[0][0]的cost = alpha_[0 + 0] + alpha_[3 + 0] = alpha_[0] + alpha_[3]，由於alpha_第一次節點初始化為0，所以cost=0。其余節點和邊計算方法類似。
前向-后向算法
有了無向圖，我們就可以在圖上進行前向-后向算法。利用前向-后向算法，很容易計算標記序列在位置i（詞）的label是yi的條件概率，以及在位置i-1（前一個詞）與位置i（當前詞）的label是yi-1與yi的條件概率。進行CRF++源碼閱讀之前先看一下條件隨機場矩陣的表示形式。對一個句子的每一個位置（單詞） i=1,2,…，n+1，定義一個 m 階矩陣（m 是標記 yi 取值的個數），i=0代表start節點， i=n+1代表stop節點。
\begin{aligned} M_i(x) &= \left \{ M_i(y_{i-1},y_i|x)\right \} \\ M_i(y_{i-1},y_i|x)&= \exp \left \{ W_i(y_{i-1} ,y_i|x)\right \}\\ W_i(y_{i-1},y_i|x)&= \sum_{k=1}^Kw_k \cdot f_k(y_{i-1},y_i,x,i) \end{aligned}
\begin{align} f_k(y_{i-1},y_i,x,i) = \left \{ \begin{aligned} &t_k(y_{i-1},y_i,x,i), \ \ k = 1,2,...,K_1 \\ &s_t(y_i,x,i), \ \ \ \ \ \ \ \ \ \ k = K_1 + l ; l = 1,2,...,K_2 \end{aligned}\right. \end{align}
Wi 的解釋：當前節點代價 + 與該節點相連的一條邊的代價。
節點之間的轉移概率，用矩陣的形式表現如下：
\begin{aligned} M_1(x) &= \begin{bmatrix} M_1(0,0|x) & M_1(0,1|x) &M_1(0,2|x) \\ 0 & 0 &0 \\ 0 & 0 &0 \end{bmatrix} \\ \\ M_2(x) &=\begin{bmatrix} M_2(0,0|x) & M_2(0,1|x) & M_2(0,2|x)\\ M_2(1,0|x) & M_2(1,1|x) & M_2(1,2|x)\\ M_2(2,0|x) & M_2(2,1|x) & M_2(2,2|x) \end{bmatrix} \\ \\ M_i(x) \ &\mathbf{has \ the \ same \ form \ with} \ M_2(X), \ i = 3,...,n\\ \\ M_{n+1}(x) &=\begin{bmatrix} 1 & 0 & 0 \\ 1 & 0 & 0 \\ 1&0&0 \end{bmatrix} \\ \end{aligned}
Mi 的解釋：以 \begin{aligned} M_2(2,1|x) \end{aligned} 為例，代表第2個位置（第2個詞）label是1，前一個詞label是2，計算Wi，再取exp后的值。接下來，我們看一下用矩陣表示的前向-后向算法。
對i = 0, 1, 2, ... n+1， 定義前向向量αi(x)，對於起始狀態i = 0：
\begin{align} \alpha_0(y|x) = \left \{ \begin{aligned} &1, \ \ y = start \\ &0, \ \ else \end{aligned}\right. \end{align}
對於之后的狀態 i=1,2,...,n+1，遞推公式為：
\begin{aligned} a_i^T(y_i|x) = a^T_{i-1}(y_{i-1}|x)M_i(y_{i-1},y_i|x) \end{aligned}
假設label個數是m，α是m*1的列向量，Mi(yi-1,yi|x) 是m*m的矩陣，α解釋：前一個單詞每個節點的α分別乘以（與當前節點相連的邊的代價 + 當前節點的代價），再求和 。
同樣，后向算法β計算， 對於i = 0, 1, 2, ..., n+1，定義后向向量βi(x)：
\begin{align} \beta_{n+1}(y_{n+1}|x) = \left \{ \begin{aligned} &1, \ \ y_{n+1} = stop \\ &0, \ \ else \end{aligned}\right. \end{align}
向前遞推公式如下：
\begin{aligned} \beta_i(y_i|x) = M_i(y_i,y_{i+1}|x)\beta_{i+1}(y_{i+1}|x) \end{aligned}
βi是m*1的列向量， Mi(yi,yi+1|x)是m*m的矩陣。β解釋：（當前詞與下一個詞連接的邊的代價 + 下一個詞的代價） 分別乘以下一個詞的β，再相加。
由前向-后向向量定義不難得到：
\begin{aligned} Z(x) = a_n^T(x) \cdot \mathbf{1} = \mathbf{1}^T \cdot \beta_1(x) \end{aligned}
需要注意一下，矩陣表示形式的代價是對特征函數乘以權重加和后再取exp的值， 而上面的CRF++ calcCost函數中並沒有取exp值。
接下來繼續看下α和β在CRF++中是如何計算的。在gradient函數中調用的forwardbackward函數即是這部分的核心代碼，具體如下：
void TaggerImpl::forwardbackward() { if (x_.empty()) { return; } for (int i = 0; i < static_cast<int>(x_.size()); ++i) { //前向算法 for (size_t j = 0; j < ysize_; ++j) { node_[i][j]->calcAlpha(); } } for (int i = static_cast<int>(x_.size() - 1); i >= 0; --i) { //后向算法 for (size_t j = 0; j < ysize_; ++j) { node_[i][j]->calcBeta(); } } Z_ = 0.0; for (size_t j = 0; j < ysize_; ++j) { //計算Z(x) Z_ = logsumexp(Z_, node_[0][j]->beta, j == 0); } return; } void Node::calcAlpha() { alpha = 0.0; for (const_Path_iterator it = lpath.begin(); it != lpath.end(); ++it) { //這里遍歷當前節點的左邊(path)的集合， 對應的就是Mi(yi-1,yi|x)矩陣中的某一列 alpha = logsumexp(alpha, (*it)->cost +(*it)->lnode->alpha, (it == lpath.begin())); //函數里面回取exp，因此邊的代價 + 上一個節點的α，會轉化成相乘，取完exp還會再取log，取log為了方式直接exp導致的溢出 } alpha += cost; //統一加上當前節點的代價， Mi(yi-1,yi|x)每列中每個元素都加了當前節點的代價， 只不過CRF++是在后面統一加上 } 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; //這里需要注意，在矩陣的推導過程中，沒有加當前節點的代價，但是CRF++里面加了， 后續我們會看到有一個減當前節點代價的一段代碼 } // log(exp(x) + exp(y)); // this can be used recursivly // e.g., log(exp(log(exp(x) + exp(y))) + exp(z)) = // log(exp (x) + exp(y) + exp(z))
// 這部分取log的操作是為了防止直接取exp溢出，具體的解釋以及推導參考 計算指數函數的和的對數
inline double logsumexp(double x, double y, bool flg) { if (flg) return y; // init mode const double vmin = std::min(x, y); const double vmax = std::max(x, y); if (vmax > vmin + MINUS_LOG_EPSILON) { return vmax; } else { return vmax + std::log(std::exp(vmin - vmax) + 1.0); } }
閱讀完上述代碼會發現，這里的α計算除了沒有對最終結果取exp以外，跟上面矩陣推導的α計算是一樣的。可以利用矩陣方法和CRF++的算法具體算一下α或β的值，對比一下理解的會更深， 這個過程並不復雜。
概率計算
有了α和β，就可以進行條件概率和期望的計算。一個句子在位置i的label是yi的條件概率，以及在位置i-1與位置i標記為yi-1與yi的概率：
\begin{aligned} P(Y_i= y_i|x) &= \frac{a_i^T(y_i|x) \beta_i(y_i|x)}{Z(x)} \\ P(Y_{i-1} = y_{i-1} ,Y_i= y_i|x) &=\frac{a_{i-1}^T(y_{i-1}|x)M_i(y_{i-1},y_i|x)\beta_i(y_i|x)}{Z(x)} \end{aligned}
第一個式子可以說是節點的概率，第二個式子是節點之間邊的概率。有了條件概率，就可以計算特征函數fk 關於條件分布 P(Y|X) 的數學期望是：
\begin{aligned} E_{p(Y|X)}[f_k] &= \sum_yP(y|x)f_k(y,x) \\ &=\sum_{i=1}^{n+1}\sum_{y_{i-1}\ y_i}f_k(y_{i-1},y_i,x,i) \frac{a_{i-1}^TM_i(y_{i-1},y_i|x)\beta_i(y_i|x)}{Z(x)} \end{aligned}
計算特征函數的期望是因為后續計算梯度的時候會用到。這里，如果fk是unigram特征（狀態特征），對應的條件概率是節點的概率， 如果是bigram特征（轉移特征），條件概率就是邊的概率。繼續看下CRF++中是如何計算條件概率和特征函數的期望的，代碼在gradient函數中：
for (size_t i = 0; i < x_.size(); ++i) { //遍歷每一個節點的，遍歷計算每個節點和每條邊上的特征函數，計算每個特征函數的期望 for (size_t j = 0; j < ysize_; ++j) { node_[i][j]->calcExpectation(expected, Z_, ysize_); } } void Node::calcExpectation(double *expected, double Z, size_t size) const { //狀態特征的期望 const double c = std::exp(alpha + beta - cost - Z); //這里減去一個多余的cost，剩下的就是上面提到的節點的概率值 P(Yi=yi | x)，這里已經取了exp，跟矩陣形式的計算結果一致 for (const int *f = fvector; *f != -1; ++f) { expected[*f + y] += c; //這里會把所有節點的相同狀態特征函數對應的節點概率相加，特征函數值*概率再加和便是期望。由於特征函數值為1，所以直接加概率值 } for (const_Path_iterator it = lpath.begin(); it != lpath.end(); ++it) { //轉移特征的期望 (*it)->calcExpectation(expected, Z, size); } } void Path::calcExpectation(double *expected, double Z, size_t size) const { const double c = std::exp(lnode->alpha + cost + rnode->beta - Z); //這里計算的是上面提到的邊的條件概率P(Yi-1=yi-1,Yi=yi|x)，這里取了exp，跟矩陣形式的計算結果一致 for (const int *f = fvector; *f != -1; ++f) { expected[*f + lnode->y * size + rnode->y] += c; //這里把所有邊上相同的轉移特征函數對應的概率相加 } }
至此，CRF++中前后-后向算法、條件概率計算以及特征函數的期望便介紹完畢，接下來看看如何計算似然函數值和梯度。
計算梯度
條件隨機場的訓練，我們這里主要看CRF++中應用的LBFGS算法。先做簡單的推導， 再結合實際的CRF++源碼去理解。條件隨機場模型如下：
\begin{aligned} P_w(y|x) = \frac{\exp \left \{ \sum_{k=1}^K w_kf_k(x,y)\right \}}{ \sum_y \left \{ \exp \sum_{i=1}^n w_if_i(x,y)\right \}} \end{aligned}
\begin{aligned} f_k(y,x) = \sum_{i=1}^nf_k(y_{i-1},y_i,x,i), k=1,2,...,K \end{aligned}
訓練函數的對數似然如下：
\begin{aligned} L(w) &= \log \prod_{t}P_w(y^t|x^t) \\ &= \sum_{t} \log P_w(y^t|x^t) \\ &= \sum_{t} \left \{ \sum_{k=1}^Kw_kf_k(y^t,x^t)-\log Z_w(x) \right \} \end{aligned}
t代表所有的訓練樣本， 一般使用m來表示，但是上面已經把m給用了， 為了避免歧義， 我們用t來表示訓練樣本。我們求似然函數最大值來求解最優參數w，同時也可以對似然函數加負號，通過求解最小值來求最優的w。這里我們與CRF++保持一致，將似然函數取負號，再對wj求導，推導如下：
\begin{aligned} \frac{\partial L(w)}{\partial w_j} &= \sum_{t} \left \{ \frac{\sum_y \left \{ f_i(x^t,y^t)\exp \sum_{i=1}^K w_if_i(x^t,y)\right \}}{Z_w(x)} - f_j(y^t,x^t) \right \} \\ &= \sum_{t} \left \{ \sum_y P(y|x^t)f_j(y, x^t) - f_j(y^t,x^t) \right \} \\ &= \sum_{t} \left \{ E_{P(y|x)}[f_j(y,x)] - f_j(y^t,x^t) \right \} \end{aligned}
對於一個句子來說，特征函數的期望減去特征函數真實值就是我們要計算的梯度，Σt 代表對所有句子求和得到最終的梯度。接下來看下CRF++中是如何實現的，代碼還是在gradient函數中：
for (size_t i = 0; i < x_.size(); ++i) { //遍歷每一個位置（詞） for (const int *f = node_[i][answer_[i]]->fvector; *f != -1; ++f) { //answer_[i]代表當前樣本的label，遍歷每個詞當前樣本label的特征，進行減1操作，遍歷所有節點減1就相當於公式中fj(y,x) --expected[*f + answer_[i]]; //狀態特征函數期望減去真實的狀態特征函數值 } s += node_[i][answer_[i]]->cost; // UNIGRAM cost 節點的損失求和 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]) { for (const int *f = (*it)->fvector; *f != -1; ++f) { --expected[*f +(*it)->lnode->y * ysize_ +(*it)->rnode->y]; //轉移特征函數期望減去真實轉移特征函數值 } s += (*it)->cost; // BIGRAM COST 邊損失求和 break; } } } viterbi(); // call for eval() 調用維特比算法做預測，為了計算分類錯誤的次數，算法詳細內容下篇介紹 return Z_ - s ; //返回似然函數值，看L(w)推導的最后一步，大括號內有兩項，其中一項是logZw(x)，我們知道變量Z_是沒有取exp的結果，我們要求這一項需要先對Z_取exp，取exp再取log相當於還是Z_，因此 logZw(x) = Z_ //再看另一項，是對當前樣本代價求和，正好這一項是沒有取exp的因此該求和項就等於s， 之前說過CRF++是對似然函數取負號，因此返回Z_ - s
至此，一個句子的似然函數值和梯度就計算完成了。公式的Σt 是對所有句子求和，CRF++的求和過程是在run函數調用gradient函數結束后由線程內匯總，然后所有線程結束后再匯總。runCRF函數剩下的代碼便是所有線程完成一輪計算后的匯總邏輯，如下：
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]; //梯度值求和 } } size_t num_nonzero = 0; if (orthant) { // L1 根據L1或L2正則化，更新似然函數值 for (size_t k = 0; k < feature_index->size(); ++k) { thread[0].obj += std::abs(alpha[k] / C); if (alpha[k] != 0.0) { ++num_nonzero; } } } else { //L2 num_nonzero = feature_index->size(); for (size_t k = 0; k < feature_index->size(); ++k) { thread[0].obj += (alpha[k] * alpha[k] /(2.0 * C)); thread[0].expected[k] += alpha[k] / C; } } ...省略代碼 if (lbfgs.optimize(feature_index->size(), &alpha[0], thread[0].obj, &thread[0].expected[0], orthant, C) <= 0) { //傳入似然函數值和梯度等參數，調用LBFGS算法 return false; }
最終調用LBFGS算法更新w，CRF++中的LBFGS算法最終是調用的Fortran語言編譯后的C代碼，可讀性比較差，本篇文章暫時不深入介紹。至此，一次迭代的計算過程便介紹完畢。
總結
通過這篇文章的介紹，已經了解到了CRF++如何構建無向圖、如何計算代價、如何進行前向-后向算法、如何計算特征函數的期望以及如何計算梯度。寫這篇文章耗時最長，花了整整一天的時間。力求這篇文章通俗易懂，理論結合實踐。希望能夠把條件隨機場這個比較枯燥的算法詮釋好。文中有可能仍然有表達不通順或者表達不通俗的地方，甚至可能會有表達錯誤的地方，如果存在上述問題歡迎評論區留言，我將第一時間更新。