机器学习&数据挖掘笔记_21（PGM练习五：图模型的近似推理）

 

　　前言：

　　此次练习完成的是图模型的近似推理，参考的内容是coursera课程：Probabilistic Graphical Models . 上次实验PGM练习四：图模型的精确推理 中介绍的是图模型的精确推理，但在大多数graph上，其精确推理是NP-hard的，因此有必要采用计算上可行的近似推理。本实验中的近似推理分为2个部分，LBP(loop belief propagation算法)和MCMC采样。实验code可参考：实验code可参考网友的：code.

 

　　算法流程：

　　LBP(loop belief propagation)算法近似推理：

　　(a): 输入factor list F和Evidence E.

　　(b): 用evidence E掉F中的每一个factor, 获得新的F.

　　(c): 利用F中的factor构造bethe cluster graph，该graph中的factor要么是single的要么是pairwise的.

　　(d): 按照某种策略选择一条message来发射，并更新message passing先后的MESSAGE矩阵。

　　(e): 继续步骤(d)直到收敛，收敛是指先后MESSAGE矩阵中的message的val值都很是接近.

　　(f): 步骤(e)完成后，说明bethe cluster graph已是calibration的. 求某个变量的边缘几率时，找到该graph中含有该变量的某个cluster，并对其求和积分掉其它变量，最后归一化便可。

　　Gibbs算法图模型的近似推理过程以下：

 　　

　　由此可知，图摸型的状态指的是全部节点的一个assignment,并非指单独分析某一个节点。

　　Metropolis-Hastings(简称MH算法)算法流程以下：

 　　

　　由上面的流程可知，MH算法与Gibbs不一样之处在于，它每次更新state时是将图中全部节点都更新，而不是一个个来。其最关键之处是建议分布的设计。

 　　

　　matlab知识：

　　在matlab中尽可能少使用for-loop,由于速度很是慢。

　　[I,J] = ind2sub(siz,IND):

　　其中siz是一个线性递增矩阵(从1开始按列递增)的长和宽，IND为须要进行index的向量值。返回的I和J就是这些值所在的行和列。

　　ind = find(X, k):

　　找到矩阵X中元素不为0的最多前k个元素，ind为这k个元素的下标(linear indices ).记住find函数找到的是下标index便可。

　　S = sparse(i,j,s,m,n):

　　生成一个稀疏矩阵，矩阵中的下标分别存在向量I,j中，对应的值存在向量s中。生成的稀疏矩阵大小为m*n.

　　n = nnz(X):

　　算出矩阵X中非0元素的个数。

 

　　实验中一些函数简单说明：

　　V2F = VariableToFactorCorrespondence(V, F):

　　该函数的做用是找出变量所对应的factor集。其中V就是变量构成的向量，F是全部的factor构成的向量。

　　[toy_network, toy_factors] = ConstructToyNetwork(on_diag_weight, off_diag_weight):

　　该函数的做用是构造一个pair-wise的小网络，网络大小4×4,每一个节点都只能取0和1. 网格中每条edge上的2个节点若是取值相同的话，则边的权值大小为on_diag_weight, 反之为off_diag_weight. 返回值toy_network就是所构成的网络结构体，里面包含的成员有names、card、edges、var2factors. 下面是一个toy network的例子：

 　　

　　[i, j] = NaiveGetNextClusters(P, m):

　　实验1的内容。参数P为cluster graph, 里面包含的变量有clusterList, edges. 参数m表示已经传递了m条消息。I→j表示须要传递的第m+1条消息。按照j下标的升序来选edge.

　　[i j] = SmartGetNextClusters(P,Messages,oldMessages,m):

　　另外一种选择须要传送message的方法。不过在函数内部参数Messages和oldMessages主要用于smoothing messages和informed message scheduling两种方法，这两种方法都是为了提升BP收敛速度的。这里没有去实现它。课件中有对应的理论。

　　[i j] = GetNextClusters(P,Messages,oldMessages,m,useSmart):

　　用参数useSmart来选择上面的两种message passing方法。

　　P = CreateClusterGraph(F, Evidence):

　　实验2的内容。用factorlist F来构造cluster graph，这比PGM练习4中构造clique tree要简单不少，由于不要求是tree，且只需edge是两相邻节点交集的子集。Evidence为观察到的变量。返回的P包含clusterList和edges两部分。该函数内部构造的是bethe cluster tree.

　　converged = CheckConvergence(mNew, mOld):

　　实验3的内容。该函数是用来判断message passing过程是否收敛，若是message passing先后的message矩阵中全部元素对应的value差值都很小(e-6)时，则表明已差很少收敛。

　　[P MESSAGES] = ClusterGraphCalibrate(P,useSmartMP):

　　实验4的内容。和PGM练习4里面的同样，也是找到一条message，而后发送，不一样的是，PGM练习4中message passing的次数固定了，而这里其次数需达到cluster graph收敛的次数才行。

　　M = ComputeApproxMarginalsBP(F,E):

　　实验5的内容。利用LBP算法进行对变量的条件几率进行近似推理。具体过程见前面的LBP算法流程。

　　LogBS = BlockLogDistribution(V, G, F, A):

　　实验6的内容。V为须要计算条件几率的变量集(若是在Gibbs状况下，则通常每次只取一个变量，此时V的长度就为1; 若是长度不为1,则需它们之间的值是相同的)，每一个变量的取值范围必须同样，G为cluster graph，F为factorlist，A为当前graph的一个assignment. 该函数的做用是计算V中每一个变量的log条件几率值. 计算依据以下：

 　　

　　由上面的公式可知，采样Xi时只将与包含Xi的那些factor相乘(不考虑归一化时).同时，程序中为了防止数值下溢(numerical underflow)问题,上面值的计算将在log空间进行。

　　[v] = randsample(vals,numSamp,replace,weightIncrements):

　　该函数实际上完成的是一个多项式分布采样。其中参数vals通常表示采样的最大值，numSamp表示须要采样的样本数。weightIncrements表示多项式分布的区间(由于是用0～1之间的随机数来投票).但GibbsTrans()一直经过不了测试，课程老师说matlab 2012后的版本都不行，是个bug。但从函数内部的代码也没看出有和matlab版本相关的代码，因此只要用到这个函数的，提交答案时都经过不了。下面是PGM课程团队在新的实验教程中对bug的提示:

  

　　A = GibbsTrans(A, G, F):

　　实验7的内容。转移矩阵采用Gibbs方法。A为图G的一个assignment, F为factorlist.该函数的做用是从当前的A获得下一个assignment,也保持在A中。内部调用函数BlockLogDistribution()获得了其分布，而后使用randi()函数实现多项式的随机采样。

　　M = ExtractMarginalsFromSamples(G, samples, collection_indx):

　　参数G仍为graph, sample为采样到的全部样本(包括mix-time前面的样本), collection_indx为samples中用于inference样本的下标。输出M为每一个样本的margin值。内部实现很简单，其实就是求估计变量的几率，用符合要求的样本个数除以样本总数。

　　E2F = EdgeToFactorCorrespondence(V, F):

　　算出图中每条edge所对应的factor编号。

　　[M, all_samples] = MCMCInference(G, F, E, TransName, mix_time, num_samples,  sampling_interval, A0):

　　实验8和实验12的内容。MCMC的接口函数。其中的G、F、E与前面的同样;TransName表明MCMC转移矩阵的类型; mix_time表示分布达到稳定前的时间。num_samples为所需sampling的样本。sampling_interval为采样间隔，通常都设置为1. A为Markov Chain的初始状态，即初始的assignment. 该函数输出值M求的就是每个变量的条件几率，是调用函数ExtractMarginalsFromSamples()来实现的。

　　ogp = LogProbOfJointAssignment(F, A):

　　求出assignment A与F中全部相关factor对应var处的乘积，就是联合几率。这个乘积与MH算法公式:

　　

　　里那个分子分母有关，若是A为当前时刻的，则logp值为公式中的分母;若是A为状态转移获得的下一时刻，则lopg为公式中的分子。为何这里2个会与实验教程公式1中相对应呢?主要是由于公式中的转移矩阵T此时为1,直接能够约掉。另外由于公式中是相除，因此几率值不用归一化。

　　A = MHUniformTrans(A, G, F):

　　实验9的内容。转移矩阵采用MH方法。一样是从当前的assignment获得下一个assignment，只不过这里采用的是均匀分布取下一状态. 内部可采用相似的LogProbOfJointAssignment()思想实现。

　　[sci sizes] = scomponents(A):

　　生成的sci,sizes是2个列向量。其中的sci是返回的每一个顶点所对应连通区域的标号, sizes是对应连通区域中元素的个数。

　　A = MHSWTrans(A, G, F, variant):

　　实验10和实验11的内容，实现MH转移几率。采用Swendsen-Wang建议分布的MH算法。该算法有2种建议分布类型, 二者的主要不一样之处体如今状态转移几率那里。下面是Swendsen-Wang状态转移的示意图：

 　　

　　 由图能够看出，其转移过程主要有下面几个步骤：

1. 对graph随机取一个初始状态;
2. 将graph中不一样值两节点之间的边断开，这样相同值的节点就构成了一个个小连通图;
3. 将上面获得的每一个小连通图中的每条边以几率(1-qij)随机断开，这样就获得了更多小小连通图;
4. 以均等几率选取图中的一个小连通图，并将该小图记为Y;
5. 以R几率分布对Y中节点值进行改变，不论是改变前仍是改变后，Y中全部节点的值都是同样的;

　　套用公式1中的MH理论，此时建议分布为：

    

　　两个的比为：

 　

　　其中记：

 　

　　而练习中的2个Swendsen-Wang算法不一样主要体如今qij选取不一样以及R分布的不一样。算法2中的qij计算公式为：

 　

　　其R采用的是BlockLogDistribution.m中的方法。

 

　　相关理论知识点：

　　LBP算法不只限于clique trees，它对任意的cluster graph都适应。只不过是不一样的cluster graph的收敛速度和精度不一样。因为cluster graph中有feedback环，因此它只能用于近似推理。在LBP算法中，message passing order并无严格要求，不用规定root节点。

　　若是给定样本(IID样本)后，则能够利用这些样原本估计一些统计量，且由Hoeffding Bound和Chernoff Bound理论保证，当样本数越多时，其估计越准确。Hoeffding Bound也称为additive bound, Chernoff Bound被称为Chernoff Bound.

　　若是给定统计量估计的错误率，则由这2个bound均可以计算出小于该错误率所需的最小样本数。

　　两个bound虽然存在，可是做用不大，由于有些统计量须要的样本数呈指数增长。

　　对BNs进行forward sampling比较简单(网络参数给定的前提下)，先sampling父节点，而后sampling子节点。若是已知某些变量的观察值，则直接将那些不符合观测值的样本点去掉便可。

　　Forward sampling不适用于Mns.

　　Regular Markov Chains是指任意两个状态之间都有连线，而且每个状态都有本身的一个环(self-transition).

　　无论初始化状态如何，Regular Markov Chains都收敛于一个惟一的稳定分布(离散的话则是多项式分布)，对这个分布进行采样比较简单。

　　若是某个分布P的样本很差直接采样，则咱们能够构造一个Regular Markov Chains T，使得T的稳定分布就是分布P，这样从T上采样就可获得P的样本。

　　mixing: Regular Markov Chains达到稳定状态时则称为mixed. 可是很难证实Regular Markov Chains是否已mixed，顶多能验证它有没有mixed. 采用MC sampling方法的难点就在于对mixed状态的判断。通常可从多个不一样初始状态值出发，对chain采用窗口计算一些统计量，经过可视化这些统计量来判断。

　　即便Regular Markov Chains已经mixed，稳定分布下连续采样获得的样本也是相关的，并不属于IID，因此这些样本不能直接用于统计推断。按照经验来讲，若是Regular Markov Chains收敛速度越快，则这些相邻样本的相关性越弱，也就越趋于IID，越有用。

　　前面的sampling都是假设咱们已经构造好了P所对应的Regular Markov Chain. 那么到底该怎样来构造呢？对应有不少的理论。不过对于PGM中的两种模型BNs和MNs而言，已有理论证实：对Gibbs chian的采样过程是可行的。

　　若是graph中全部的factor都是positive的，则Gibbs chain是regular的。而且Gibbs chain采样某个元素值时只与该元素所在的factor有关，与其它factor无关。固然，前提是咱们知道怎样从这些factor中来采样。

　　精确推理在general networks上是NP-hard的。

　　Reversible Chains: 知足细致平衡的chain. 细致平衡是指当前状态为x，且下一状态为x'的几率与当前状态为x', 且下一状态为x的几率相等。

　　构造分布P对应的MCMC链的一个通用的理论框架是MH(Metropolis-Hasngs)算法(Gibbs也属于MH算法的一种)。

　　MH框架下须要一个建议分布Q，以及一个接受率A。该建议分布必须是reversible的。为了缩短mix time，Q应该越宽越好，可是若是太宽则接受率又会太低。

　　通常来讲，当graph越稀疏就越适合用message passing的方法来inference，若是越dense，则越适合用sampling的方法。

　　为了使inference更准确，message passing方法中应该使cluster graph 中的cluster更大;在sampling方法中应该使建议分布每次转移覆盖更多变量。

 

　　参考资料：

       Daphne Koller，Probabilistic Graphical Models Principles and Techniques书籍第12章

       coursera课程：Probabilistic Graphical Models 

　　PGM练习四：图模型的精确推理 

       实验code可参考网友的：code.