機器學習&數據挖掘筆記_23（PGM練習七：CRF中參數的學習）

 

　　前言：

　　本次實驗主要任務是學習CRF模型的參數，實驗例子和PGM練習3中的同樣，用CRF模型來預測多張圖片所組成的單詞，咱們知道在graph model的推理中，使用較多的是factor，而在graph model參數的學習中，則使用較多的是指數線性模型，本實驗的CRF使用的是log-linear模型，實驗內容請參考 coursera課程：Probabilistic Graphical Models 中的assignmnet 7. 實驗code可參考網友的:code實驗對應的模型示意圖以下：

 　　

 

　　CRF參數求解過程：

　　本實驗中CRF模型所表示的條件機率計算公式爲:

 　 

　　其中的分母爲劃分函數，表達式爲：

 　　

　　採用優化方法訓練CRF模型的參數時，主要任務是計算模型的cost和grad表達式。其中cost表達式爲：

 　 

　　grad表達式爲：

 　

　　公式中的2個指望值表示模型對特徵的指望以及數據對特徵的指望，其表達式以下：

 　

　　在計算cost和grad時須要分別計算下面6箇中間量：

 　

　　關於這幾個中間量的計算方法，能夠參考實驗教程中的介紹，或者直接看博文後面貼出的代碼，這裏簡單介紹下其計算方法：

　　log partition function: 需先創建CRF模型對應的clique tree，並對其校訂，校訂過程當中時須要message passing，而最後passing的消息和最後一個clique factor相乘後的val之和取對數就logZ了。

　　weighted feature counts: 當訓練樣本中某個樣本及標籤的值符合CRF模型的某一個特徵時，就將該特徵對應的參數值累加，最後求和便可。

　　regularization cost: 直接計算。

　　model expected feature counts: 計算模型對特徵的指望，一樣須要用到前面校訂好了的clique tree. 當某個特徵的變量所有屬於clique tree中某個clique變量時，求出該clique對應的factor中符合這些特徵變量值的和，注意歸一化。

　　data feature counts: 在計算weighted feature counts的同時，若是某個特徵在樣本中出現，則對相應特徵計數。

　　regularization gradient term: 直接計算。

 

　　matlab知識：

　　ndx = sub2ind(siz,varargin):

　　siz爲一個矩陣的維度向量，varargin輸入的向量表示在矩陣size的位置，返回的是linear index的值。好比sub2ind([3 4],2,4)返回11,表示在3×4大小的矩陣中，第2行第4列爲矩陣的第11個元素。

　　C = horzcat(A1,...,AN):

　　將參數表示的矩陣在水平方向合成一個大矩陣。

 

　　實驗中一些函數簡單說明：

　　[cost, grad] = LRCostSGD(X, y, theta, lambda, i):

　　計算帶L2懲罰項的LR cost. 其中X是輸入矩陣，每一行表明一個樣本，y爲對應的標籤向量。theta爲LR模型的權值參數，lambda爲權值懲罰係數。i表示選擇X矩陣中第i個樣原本計算(循環取,mod實現)。cost和grad分別爲這個樣本的偏差值和輸出對權值的導數值。

　　thetaOpt = StochasticGradientDescent (gradFunc, theta0, maxIter):

　　實驗1的內容。gradFunc是函數句柄，[cost, grad] = gradFunc(theta, i),計算logistic regression在第i個樣本theta處的cost和grad.  theta0爲權值初始值,maxIter爲最大迭代次數。這裏的每次迭代只使用了一個樣本，採用隨機梯度降低法(SGD)更新權值。

　　thetaOpt = LRTrainSGD(X, y, lambda):

　　該函數完成的是用訓練樣本X和標籤y對LR進行參數優化，迭代次數和初始學習率等超參數在函數內部給定，實現該函數時需調用StochasticGradientDescent().

　　pred = LRPredict (X, theta):

　　計算樣本矩陣X在參數theta下的預測標籤pred.

　　acc = LRAccuracy(GroundTruth, Predictions):

　　計算關於真實標籤GroundTruth和預測標籤Predictions的準確率。

　　allAcc = LRSearchLambdaSGD(Xtrain, Ytrain, Xvalidation, Yvalidation, lambdas):

　　實驗2的內容。該函數是計算lambdas向量中每一個lambda在驗證集Xvalidation,Yvalidation上的錯誤率，這些錯誤率保存在輸出變量allAcc中。

　　[P, logZ] = CliqueTreeCalibrate(P, isMax):

　　實驗3的內容。對clique tree P進行校訂，在校訂過程當中同時求得劃分函數Z的log值:logZ. 求logZ時只能用sum-product不能使用max-product,因此此時的isMax=0. 其方法是將最後一次傳送的message和最後一個clique相乘獲得的factor，而後將factor中的val求和便可。

　　featureSet = GenerateAllFeatures(X, modelParams):

　　X是訓練樣本矩陣，由於在CRF中須要同時輸入多張圖片(本實驗中多張圖片構成一個單詞)，因此這裏X中的每一行表明一張圖片。結構體modelParams有3個成員：numHiddenStates，表示CRF中隱含節點的個數，這裏爲26(26個字母); numObservedStates,表示CRF中觀察節點的個數，這裏爲2(每一個像素要麼爲0,要麼爲1); lambda,權值懲罰係數。返回值featureSet包括2個成員：numParams, CRF中參數的個數，需考慮權值共享狀況，便可能有多個特徵共用一個權值; features, 裝有多個feature的向量，且每一個feature又是一個結構體。該feature結構體中有3個成員，以下：var,特徵所包含的變量;assignment,由於特徵通常爲指示函數，因此表示特徵只在assignment處的值爲1,其它處爲0; paramIdx,特徵所對應的參數在theta中的索引。

　　features = ComputeConditionedSingletonFeatures (X, modelParams):

　　計算輸入圖像中單個像素的特徵，若是輸入X爲3*32大小，由於有26個字母，因此總共的特徵數爲3*32*26=2496. 又假設每一個像素可取0或1兩個值，因此總共的參數個數爲2*32*26=1664. 很明顯有些特徵是共用相同參數的。獲得的features.var爲圖片序列的編號，features.assignment爲對應字母的編號，features.paramIdx由像素值決定參數的位置。

　　features = ComputeUnconditionedSingletonFeatures (len, modelParams):

　　len爲圖片序列中圖片的個數。若是len=3,則該函數實現後的特徵向量長度爲3*26=78，參數個數爲26. features.var爲圖片編號，features.assignment爲字母編號的，features.paramIdx位置和字母編號同樣。

　　features = ComputeUnconditionedPairFeatures (len, modelParams):

　　len爲圖片序列中圖片的個數。若是len=3, 則該函數實現後的特徵向量長度爲2*26*26=1352，參數個數爲26×26=676.  其中的features.var爲圖片序列編號的組合，features.assignment爲字母序號編號的組合，features.paramIdx爲所在字母序號中對應的位置。

　　由上面學到的3種特徵可知，特徵的var都是與輸入圖片序列的標號有關，特徵的assignment都是與字母的序號有關，paramIdx可能與字母序列以及圖片序列編號都  有關。實驗教程中給出的3種特徵以下：

 　　

　　上面第2種特徵爲應該爲conditionedPairFeatures，但和實驗code沒有對應起來，實驗code中該特徵被替換成了conditionedSigleFeatures.

　　VisualizeCharacters (X):

　　可視化樣本X，由於X中一個樣本序列，可能包含多個字母，該函數將X中所含的字母顯示在一張圖上。

　　[nll, grad] = InstanceNegLogLikelihood(X, y, theta, modelParams):

　　實驗4和5的內容。其中參數X,y,modelParams和前面介紹的同樣，注意X矩陣對CRF來講只算一個樣本。參數theta爲列向量，大小numParams x 1，是整個CRF模型中共享的參數。這2個實驗的實現主要按照博文前面介紹的算法來計算，代碼以下：

% function [nll, grad] = InstanceNegLogLikelihood(X, y, theta, modelParams)
% returns the negative log-likelihood and its gradient, given a CRF with parameters theta,
% on data (X, y). 
%
% Inputs:
% X            Data.                           (numCharacters x numImageFeatures matrix)
%              X(:,1) is all ones, i.e., it encodes the intercept/bias term.
% y            Data labels.                    (numCharacters x 1 vector)
% theta        CRF weights/parameters.         (numParams x 1 vector)
%              These are shared among the various singleton / pairwise features.
% modelParams  Struct with three fields:
%   .numHiddenStates     in our case, set to 26 (26 possible characters)
%   .numObservedStates   in our case, set to 2  (each pixel is either on or off)
%   .lambda              the regularization parameter lambda
%
% Outputs:
% nll          Negative log-likelihood of the data.    (scalar)
% grad         Gradient of nll with respect to theta   (numParams x 1 vector)
%
% Copyright (C) Daphne Koller, Stanford Univerity, 2012

function [nll, grad] = InstanceNegLogLikelihood(X, y, theta, modelParams)
    % featureSet is a struct with two fields:
    %    .numParams - the number of parameters in the CRF (this is not numImageFeatures
    %                 nor numFeatures, because of parameter sharing)
    %    .features  - an array comprising the features in the CRF.
    %
    % Each feature is a binary indicator variable, represented by a struct 
    % with three fields:
    %    .var          - a vector containing the variables in the scope of this feature
    %    .assignment   - the assignment that this indicator variable corresponds to
    %    .paramIdx     - the index in theta that this feature corresponds to
    %
    % For example, if we have:
    %   
    %   feature = struct('var', [2 3], 'assignment', [5 6], 'paramIdx', 8);
    %
    % then feature is an indicator function over X_2 and X_3, which takes on a value of 1
    % if X_2 = 5 and X_3 = 6 (which would be 'e' and 'f'), and 0 otherwise. 
    % Its contribution to the log-likelihood would be theta(8) if it's 1, and 0 otherwise.
    %
    % If you're interested in the implementation details of CRFs, 
    % feel free to read through GenerateAllFeatures.m and the functions it calls!
    % For the purposes of this assignment, though, you don't
    % have to understand how this code works. (It's complicated.)
    
    featureSet = GenerateAllFeatures(X, modelParams); %由於擬合的是條件機率，因此須要使用X

    % Use the featureSet to calculate nll and grad.
    % This is the main part of the assignment, and it is very tricky - be careful!
    % You might want to code up your own numerical gradient checker to make sure
    % your answers are correct.
    %
    % Hint: you can use CliqueTreeCalibrate to calculate logZ effectively. 
    %       We have halfway-modified CliqueTreeCalibrate; complete our implementation 
    %       if you want to use it to compute logZ.
    
    nll = 0;
    grad = zeros(size(theta));
    %%%
    % Your code here:
    % 計算cost
    ctree = CliqueTreeFromFeatrue(featureSet.features, theta, modelParams); %對整個展開的CRF對應的graph而言
    [ctree,logZ] = CliqueTreeCalibrate(ctree, 0); %對tree進行校訂，並求出劃分函數的對數
    [featureCnt,weightCnt] = WeightFeatureCnt(y, featureSet.features, theta);
    weightedFeatureCnt = sum(weightCnt);
    regCost = (modelParams.lambda/2)*(theta * theta'); 
    nll = logZ-weightedFeatureCnt+regCost;
    
    % 計算grad
    mFeatureCnt = ModelFeatureCount (ctree, featureSet.features, theta);%求模型指望時，不能使用y信息
    regGrad = modelParams.lambda* theta;
    grad = mFeatureCnt-featureCnt+regGrad;
end

%% 該函數實現的功能是對每一個特徵都創建一個factor
function ctree = CliqueTreeFromFeatrue(features, theta, modelParams)
    n = length(features);
    factors = repmat(EmptyFactorStruct(),n,1);
    for i=1:n
        factors(i).var = features(i).var; 
        factors(i).card = ones(1,length(features(i).var))*modelParams.numHiddenStates; %難道都是y變量？
        factors(i).val = ones(1, prod(factors(i).card));
        % 給該factor賦特徵值
        factors(i) = SetValueOfAssignment(factors(i), features(i).assignment, exp(theta(features(i).paramIdx)));
    end
    ctree = CreateCliqueTree(factors);
end

%% 該函數是求輸入樣本是否知足各個特徵，若是知足特徵i，則counts(i)=1,且weighted(i)填入相應的權值。
function [counts, weighted] = WeightFeatureCnt(y, features, theta) %這裏要使用y值，由於要用y來計算指示特徵
    %注意特徵向量的長度和參數向量的長度並不相同,由於多個特徵能夠共用一個參數，因此通常參數向量要短些
    counts = zeros(1,length(theta));  
    weighted = zeros(length(theta), 1);
    for i = 1:length(features)
        feature = features(i);
        if all(y(feature.var)==feature.assignment) %判斷所給的y是否知足特徵所描述的
            counts(feature.paramIdx) = 1;
            weighted(i) = theta(feature.paramIdx);
        end
    end
end

%% 該函數是計算模型的特徵指望值，利用模型對應校訂好的clique tree來計算，每一個特徵由其對應的clique中歸一化的val構成
function mFeatureCnt = ModelFeatureCount (ctree, features, theta)
    mFeatureCnt = zeros(1,length(theta)); %提早開闢空間有利於matlab運算速度
    for i = 1:length(features)
        mIdx = features(i).paramIdx;
        cliqueIdx = 0;
        for j = 1:length(ctree.cliqueList)
            if all(ismember(features(i).var,ctree.cliqueList(j).var))
                cliqueIdx = j; %在clique tree上找到包含第i個特徵全部元素的clique
                break;
            end
        end
        eval = setdiff(ctree.cliqueList(cliqueIdx).var, features(i).var);
        featureFactor = FactorMarginalization(ctree.cliqueList(cliqueIdx),eval); %獲得只包含該特徵變量的factor
        idx = AssignmentToIndex(features(i).assignment,featureFactor.card);
        mFeatureCnt(mIdx) = mFeatureCnt(mIdx) + featureFactor.val(idx) / sum(featureFactor.val);%歸一化
    end
end

 

　　相關理論知識點：

　　learning按照模型結構是否已知，數據是否徹底能夠分爲4類。好比HMM屬於結構已知但數據不徹底那一類(由於模型中的狀態變量不能觀測)。

　　PGM learning的任務有：probabilistic queries(for new instance);Specific prediction task(for new instance);Knowledge discovery(for distribution).

　　overfitting分爲參數的overfitting和結構的overfitting.

　　PGM中比較容易整合先驗知識。

　　MLE(最大似然估計)與充分統計量密切相關。

　　在BN中，若是有disjoint set的參數，則可將似然函數分解成局部的似然函數乘積。若是是table CPD的話，則局部似然函數又能進一步按照變量中每一維分解。

　　數據越少需使用越簡單的模型，這樣泛化性能纔好，不然很容易過擬合。

　　MLE的缺點是不能很好的判斷其參數估計的可信度。好比在下列兩種狀況下用拋硬幣估計硬幣朝上的機率時，使用MLE有結果：(a). 10次有7次朝上，這時估計硬幣朝上的機率爲0.7;(b). 10000次有7000次朝上,硬幣朝上的機率也被判爲0.7. 雖然估計的結果都爲0.7,但很明顯第二種情形的估計結果更可信，第一種情形過擬合。

　　爲了克服MLE的缺點，可採用Bayesian估計(也叫最大後驗估計)，Bayesian估計的抗噪能力更強，相似於權值懲罰。貝葉斯估計是把模型中的參數也當作是一個隨機變量，這樣在估計該參數時會引入該參數的先驗。爲了達到共軛分佈的目的，該先驗分佈通常取Dirichlet分佈。

　　2個變量的Dirichlet分佈曲線能夠在2維平面上畫出，是由於Dirichlet變量之間有和爲1的約束，至關於減小了一個自由量。

　　在貝葉斯網絡中，若是參數在先驗中是相互獨立的，則這些參數在後驗中也是相互獨立的。

　　劃分函數的對數對參數求導，直接把劃分函數按照定義代入公式，可得求導結果爲模型對特徵的指望。在採用MLE估計時，能夠推出loss對參數的導數爲數據對特徵的指望減去模型對特徵的指望, MRF下的證實以下：

 　　

　　因此最終MRF的梯度公式爲：

 　　

　　相對應CRF是用的log條件似然函數，其梯度計算公式爲：

 　　

　　MRF和CRF的loss函數的優化都屬於凸優化。而且二者計算梯度時都須要圖模型的inference,其中MRF每次迭代需inference一次，而CRF每次迭代的每一個樣本處需inference一次。從這點看貌似CRF的計算量要大些。不過因爲MRF須要擬合聯合機率，而CRF不須要，因此總的來講，CRF計算量要小些。

　　在MRF和CRF中也能夠採用MAP估計，這裏的先驗函數通常爲高斯先驗和拉普拉斯先驗，二者的公式分別爲：

 　　

　　

　　其中高斯先驗相似L2懲罰，使得不少參數在0附近(由於接近0時導數變小，懲罰力度變小)，但不必定爲0,因此對應的模型很dense，而拉普拉斯先驗相似L1懲罰，使得不少參數都爲0(由於其導數不變，即懲罰力度沒變)，對應模型比較sparse.

 

　　參考資料：

       實驗code可參考網友的:code

       PGM練習三：馬爾科夫網絡在OCR上的簡單應用

       coursera課程：Probabilistic Graphical Models