基于混合粒子群优化的贝叶斯网络结构学习方法

尉永清，陈小雪，伊 静，孟媛媛

1(山东师范大学 信息科学与工程学院，济南 250358)2(山东警察学院 公共基础部，济南 250014)3(山东省分布式计算机软件新技术重点实验室，济南 250358)4(山东建筑大学 信息科学与工程学院，济南 250101)

1 引 言

贝叶斯网络以概率论和图论为理论基础，它将概率统计与图论相结合，主要应用于解决复杂系统不确定性数据推理和分析，已经成为处理不确定性问题的有效工具[1，2].目前贝叶斯网络的研究是人工智能[3，4]、基因分析[5]、金融投资与分析[6]、故障诊断[7]、信息融合[8]等领域不确定知识表达和推理技术的主流方法.

贝叶斯网络研究主要集中在贝叶斯网络推理、贝叶斯网络学习和基于贝叶斯网络的应用三个方面.其中，贝叶斯网络学习作为贝叶斯网络的基础，主要是通过分析数据来获得贝叶斯网络的过程.贝叶斯网络学习主要包括参数学习和结构学习两种.参数学习是指已知网络结构，通过不断学习来确定网络参数.贝叶斯网络结构学习指的是通过结合包含专家知识在内的先验信息，以寻找出与样本数据集拟合得最好的网络结构[9].贝叶斯网络结构学习主要包括基于独立性测试的方法[10，11]和基于评分搜索[12-14]的方法.基于评分搜索的方法主要取决于评分函数的确定和搜索策略的选择.与参数学习相比结构学习复杂的多，同时，贝叶斯网络结构学习问题已经证明是一个NP难问题[15，16]，一般采用启发式搜索算法解决评分搜索问题.贝叶斯结构网络结构往往是由专家给出，相对费时、费力.当面对特别复杂的系统时，仅仅依靠专家经验来构建贝叶斯网络基本不可能实现，并且构造正确的贝叶斯网络结构是进行参数学习的前提，因此本文主要研究网络结构学习.

粒子群优化算法[17](Particle Swarm Optimization，PSO)是Kennedy和Eberhart等人通过对鸟群飞行、聚集行为的研究，在1995年提出的一种群体智能优化算法.该算法结构简单、参数较少，在高维空间函数寻优方面具有收敛速度快，解的质量高、鲁棒性好、极易实现等优点.因此在函数优化、神经网络训练 、模式分类、控制工程等领域受到了广泛的关注和学习[18-21].将智能算法与贝叶斯网络结合成为学者的研究热点.2011年冀俊忠等人提出了一种基于蚁群优化的贝叶斯网络结构学习算法.2013年汪春峰等人将无约束优化和遗传算法结合，提出了一种学习贝叶斯网络结构的限制型遗传算法.2014年张平等人将人工蜂群算法应用到贝叶斯网络结构学习中.2014年刘扬等人提出一种信息论结合粒子群优化的贝叶斯网络结构学习算法.2014年高晓光等人提出一种互信息限制的贝叶斯网络结构学习算法.

鉴于此，本文将贝叶斯网络结构学习和粒子群优化算法看成一个组合优化的问题，将粒子群优化算法应用于贝叶斯网络结构学习中，在此基础上，并将遗传算法良好的并行计算能力进行有效的结合，利用混合优化方法寻找最优的贝叶斯网络.同时，结合贝叶斯网络结构特点，设计基于遗传算子的位置更新策略，主要是将粒子群算法中的位置更新策略与遗传算子中的变异算子和交叉算子相结合，以增加学习的精度和效率.最后以标准的Alarm 和Aisa 网络为实例，并与其他算法相比较体现了本文算法在贝叶斯网络结构学习中较强的学习能力.

2 相关介绍

2.1 贝叶斯网络

贝叶斯网络(Bayesian Networks，BNs)作为概率信息的载体，主要由联合概率分布的图形来表示.通常一个贝叶斯网络由有向无环图(简称DAG)和条件概率表(CPT)两部分组成：DAG中的每个节点代表一个随机变量，根据节点之间的概率依赖关系由有向边连接；CPT主要表示节点间的依赖强度.

对于一个n元有限随机变量V={v1，v2，v3…vn}，Θ为变量X的联合概率分布.则贝叶斯网络可由二元组BN=(G，Θ)来表示，其中G=(V，E)是一个有向无环图，V={v1，v2，v3…vn}是一个节点集合，V中的变量与图中的节点相对应.E为有向边集合，表示随机变量间的依赖关系.Θ是一个网络参数向量，为节点的条件概率分布集合，表示节点之间的依赖程度向量Θ由若干分量组成，用Θ={P[vi|π(vi)]，vi∈V}表示贝叶斯再给定节点vi的父母节点时的条件概率，且每个分量表示网络结构中节点vi的条件概率，为一个网络参数.π(vi)表示节点vi的父母节点，根据概率论的链式规则，变量vi的联合概率分布为：

(1)

若没有父节点，则π(vi)=Φ.

2.2 评分标准

贝叶斯网络结构学习主要思想是首先结合先验知识找到一个和样本数据集合拟合最好的网络结构，然后通过评分函数对网络结构与样本集合匹配程度进行测试选取最优的网络.对于基于评分搜索的方法主要解决的是评分函数的选取和搜索方法的选择两个问题.其中评分函数是衡量生成的贝叶斯网络结构与数据样本匹配程度的度量标准.本文选用使用广泛的贝叶斯信息标准(Bayesian information criterion，BIC)作为适应度函数，适应度值越大则表明匹配程度越高.

对于一组随机变量V={v1，v2，v3…vn}，D={D1，D2，D3…Dn}是与这组变量对应的独立分布的样本集.G是v1，v2，v3…vn为节点的贝叶斯网络.贝叶斯网络结构学习的主要目的是首先结合先验知识找到一个和样本数据集合D拟合最好的网络结构G.然后通过评分函数对网络结构与样本集合匹配程度进行测试选取最优的网络.BIC评分函数是在样本满足独立同分布的假设前提下，用对数似然度量结构与数据的拟合程度.BIC评分函数为[24，25]：

(2)

评分函数确定之后，接下来的任务就是选择一个合适的搜索算法应用于贝叶斯网络结构中以寻找出分数最高的网络结构.对于一个包含n个节点的有向无环图，f(n)是组成无环图所有可能的个数，Robinson给出的f(n)的计算公式如下[26]：

f(1)=1

(3)

从上式中可以看出，随着n的增大贝叶斯网络的空间不断增大.贝叶斯网络结构学习问题已经证明是一个NP-Hard问题，故采一般地确定性的精确算法寻找最优的网络结构，通常很难达到目标.

3 贝叶斯网络结构的学习过程

3.1 粒子群优化算法

粒子群优化算法(Particle Swarm Optimization，PSO)起源于对鸟群捕食过程中的迁徙和聚集的模拟，是一种基于全局搜索的群体智能算法[27].在粒子群算法中，粒子可表示成一个二元组(X，V)，其中，X、V分别表示粒子的位置、速度.每个粒子都代表着搜索空间中的一个解，所有的粒子都对应一个由优化函数所决定的适应度值，解的质量通过适应度值的大小来衡量.粒子群优化算法首先初始化一群随机粒子，然后在每一次迭代中，每个粒子根据自身的最优值(Pbest)和群体的最优值(Gbest)不断地更新自己的位置找到最优解.

假设在d维空间中，粒子的位置定义为Xi=(Xi1，Xi2…Xid)，Vi=(Vi1，Vi2…Vid)，粒子的速度和位置更新公式如下[28]：

(4)

(5)

文献[29]中证明了粒子的进化过程与速度无关，从而提出了一种简化的粒子群算法(Simple Particle Swarm Optimization，SPSO)如下：

(6)

3.2 粒子编码

由于PSO算法的搜索空间是由贝叶斯网络结构组成.假设一个贝叶斯网络结构有n个节点组成，则可以用一个n×n维的邻接矩阵来表示.本文采用文献[28，30]矩阵编码方式进行编码，如果把“粒子的位置”也定义为一个有向无环图，则粒子的当前位置也可以用矩阵G表示，G={gij}，其中gij定义如下：

(7)

本文算法中的贝叶斯网络结构的个体可表示为：g11g21…gn1g12g22…gn2…g1ng2n…gnn

图1 粒子的位置表示Fig.1 Particle position

根据上述编码方式，图1的贝叶斯网络结构编码为0011001000010000

3.3 基于遗传算子的位置更新策略

通过对粒子的变异操作，提高了全局搜索能力，避免陷入局部最优解；将粒子分别与局部最优值和个体最优值交叉，保护了粒子的最优状态值也同时增强了种群的多样性.在提高搜索效率的同时减小了粒子更新的随机性，通过对贝叶斯网络逻辑关系的推理，获得准确的贝叶斯网络结构.采用标准的Alarm[22]和Asia[23]网络验证了本文算法用于贝叶斯网络结构学习的可行性及优越性，与其他算法相比取得了较好的实验效果.

3.3.1 变异算子操作

(8)

式中，r1是[0，1)的随机数，A代表变异操作，ω代表变异概率.

本文算法对贝叶斯网络结构设计了3种变异操作如图2所示.

图2 变异操作后的扑拓图Fig.2 A graph of variation operations

图2中(a)表示当前模型，图2(b)表示从1→4减边操作，图2(c)表示从2→3反转边操作，图2(d)表示从2→4加边操作，图2(e)表示从2→4多边操作，因为导致回路，记为不合法的多边操作.通过加边、删边或反转边三种操作之后，再计算网络得分.由图2可以看出，粒子在交叉变异之后可能产生非法的拓扑图，如图2(e)所示.对此本文采用文献[34]所提出的修复算子对非法的贝叶斯网络结构进行去环操作.具体修复步骤如下：

1)求出与网络拓扑图相对应矩阵的传递闭包；

2)通过判断闭包矩阵对角线上的元素是否全为0，判断是否为合法的网络结构.如果全为0，代表合法的网络结构；否则，保留主对角线上非0元素对应的节点(这些节点均位于环内).

3)任取环内的一点，求出它位于闭环内的所有父节点；

4)对于这些父节点指向该节点的任一边，进行删除或者反转操作，使得网络结构中不存在有向环.

3.3.2 交叉算子操作

交叉是产生新个体的核心，因为交叉后的子代是父代的继承和重组，所以种群适应度较高，交叉的结果直接影响后代的质量和收敛的结果.为了丰富种群的多样性本文采用两点交叉的方式.在粒子与个体最优粒子交叉的过程中，将会产生新粒子，主要是由父母粒子的公共部分和随机部分组成，公式如下：

(9)

(10)

式中，Cp代表新粒子与个体最优粒子的交叉操作，Cg表示个体最优粒子与全局最优粒子的交叉过程.c1，c2是学习因子，c1代表个体最优粒子的交叉概率，c2全局最优粒子的交叉概率，r2，r3是[0，1)的随机数.

3.3.3 粒子位置更新

粒子群算法通过遗传算子的操作，粒子的位置更新公式为：

(11)

惯性权重可以衡量算法的全局搜索能力和局部搜索能力，当ω过大会时有利于全局搜素，当ω过小时有利于局部搜索.学习因子是调节粒子自身经验和群体经验的权重，主要反应粒子之间的交流，影响粒子的运动轨迹.c1，c2过大过小都会对结果造成不利影响.通过以上分析，通过公式(12～14)线性的调节ω，c1，c2.

(12)

(13)

(14)

式中，t表示当前迭代次数，Itmax代表最大的迭代次数.ωs，c1s，c2s代表初始值，ωe，c1e，c2e代表最终值.

3.4 算法实现步骤

本文提出一种混合粒子群优化的贝叶斯网络结构学习算法，简称MMPSO算法.在该算法中，粒子所处的位置代表一种网络结构图，用邻接图来表示每一个候选网络结构，如何产生初始解以及对每一个解的邻域进行搜索是解决问题的关键.

由于初始种群的选择对PSO算法的寻优性能影响较大，对此，本文通过采用文献[35]提出的最大权生成树(MWST)算法，产生网络结构的初始边集，然后生成与贝叶斯网络结构拟合度最好的树结构.以此作为初始结构图模型，通过对其任意添加边、删除边或反转边的操作以产生该图的所有邻接图.在相应邻接矩阵中选用一定数量的矩阵作为初始种群，然后在利用MMPSO搜索出最优的网络结构.具体的步骤如下：

实验步骤：

Step1. 设置初始参数：种群规模N，交叉概率，最大迭代次数等.

Step2. 产生初始种群：通过MWST算法得到初始结构图，并由此生成该图的所有邻接图，选取N个结构图作为初始粒子群；

图3 本文算法的流程图Fig.3 Flow chart of algorithm

Step3. 根据公式(2)计算每个粒子的BIC评分值作为适应度值；并且根据个体极值找出全局最大的评分值和对应的初始全局粒子，同时更新每个粒子的最优解；

Step4. 根据公式(8)，调节权重系数对粒子进行变异操作，产生新粒子，以保证良好的全局搜索能力和局部搜索能力.

Step5. 根据公式(9)，将新粒子与个体最优粒子进行交叉操作，同时更新个体极值与全局极值；

Step6. 根据公式(10)，将粒子的个体极值与全局极值进行交叉操作，更新新的个体极值与全局极值.

Step7. 通过公式(2)重新计算每个粒子的BIC值，并更新种群的全局最优解.

Step8. 判断粒子是否早熟收敛，如果早熟收敛则执行混沌优化搜索策略，否则继续执行粒子群算法.

Step9. 检查终止条件(达到足够好的位置或达到最大迭代次数)；如果达到条件，运行终止并输出全局最优解.否则，返回Step 4.

4 实验结果与分析

4.1 数据集

本文仿真实验的硬件环境是内存4G，CPU为Intel(R)CoreTM，2 Duo 3.30GHz.实验平台Windows7，程序实现采用Matlab2014b，以及软件工具包FullBNT-1.0.4.

表1 标准测试网络数据组成Table 1 Standard test network data

为了测试MMPSO的性能，本文采用通用的Benchmark数据集Alarm网络、Asia网络、网络来完成实验.具体数据组成如表1所示.

4.2 参数设置及评价指标

根据标准的Alarm、Asia网络，首先利用BNT工具依照标准概率表生成样本量为500、1000、2000、3000、5000的数据集，然后使用本文算法分别对产生的多组数据进行结构学习，并将运行的结果与MMHC算法[36]、GA算法[31]、PSO算法[28]、MMACO算法[37]进行对比，以此来验证本文算法的效果.在仿真实验中，设置初始粒子种群为20，算法的参数设置如下：(注：初始化种群的规模为N，蚂蚁数量为m，信息挥发因子为ρ，权重参数α，β，交叉概率记为pc，变异概率记为pm)最大迭代次数Tmax=100.

表2 参数设置Table 2 Parameter setting

实验采用的评价指标是BIC评分值和结构汉明距离(SHD)，评价各种算法所学到的最优结构与标准网络结构之

间的差异.其中

1)标准BIC评分准则，表示算法学习到最优结构时的BIC得分值，分值越高结果越好.

图4 标准的Aisa网Fig.4 Aisa network图5 学习后的Aisa网Fig.5 Learned Aisa network

2)结构汉明距离(SHD)：用来度量学习所得网络结构与真实网络结构之间的差异程度，汉明距离用网络结构中多余边、反转边和丢失边三者之和来表示.取值越小代表结果越好.

4.3 实验结果分析

为了验证本文算法用于贝叶斯网络结构的效果，主要从三个方面进行实验：算法的可行性、算法的准确性、算法的收敛性.

图6 标准的Alarm网Fig.6 Alarm network图7 学习之后的Alarm网Fig.7 Learned Alarm network

4.3.1 可行性分析

本文以经典的Alarm网和Asia网作为仿真模型来评价算法的性能.图4，图5分别标准和学习之后的Aisa网络，图6，图7分别标准和学习之后的Alarm网络.从上图中可以看出，MMPSO在Aisa网的结构中，学习的最优效果是仅出现从long到smoking的反向边，准确率高于其他算法.在复杂度更

表3 不同算法对Alarm网络学习的结果Table 3 Results of different algorithms for learning Alarm networks

高的Alarm网的结构学习中，与标准的网络对比，分别多3条边，少一条边，反向一条边.(多11至15，27至37，33至15；少：17-26，27-36；反：30至29)最有效果是仅有一条反向边，较其他算法精确率较高.

表4 不同算法对Aisa网络学习的结果Table 4 Results of different algorithms for learning Aisa networks

4.3.2 准确性分析

为了验证MMPSO算法学习网络结构的准确性，将本文算法与与MMHC算法[36]、GA算法[31]、PSO算法[28]、MMACO算法[37]分别对产生的多组数据进行结构学习，并将运行结果进行比较.表3，表4给出了各种算法独立运行10次得到的BIC评分的均值和SHD值，其中表示独立运行10次BIC评分的均值和方差.

从表3至表4的数据中可以看出，在样本容量相同的情况下，MMPSO算法学习的网络结构更接近真实网络，说明了本文算法的有效性.在Alarm数据集上，MMHC学习效果最差，本文算法所学习的效果最好，MMACO次之，虽然随着样本数量的逐渐增大，MMACO与MMPSO的得分值逐渐接近，但是本文算法在保证得分较高的情况下，SHD取值最小.在Alarm数据集上，与其他算法相比，MMPSO的SHD取值与MMABC取值相同，但MMPSO的得分值高于MMACO，得到了比较好的结果.

4.3.3 收敛性分析

为了验证MMPSO算法的收敛性能，以Aisa为例，用2000组数据做实验并将本文算法与PSO和MMACO算法的结果进行对比.用2000组数据做实验并将MMPSO与PSO和MMACO算法的结果进行对比.

表5 3种算法学习Aisa网络的结果比较Table 5 Results of learning Aisa networks

设置每种算法的最大迭代次数为50次，初始粒子为20，表5列出了三种算法分别运行10次，学习Aisa网络的BIC得分和相应有效时间的均值.由表8和收敛结果图可以看出，

图8 学习Aisa网络的结果收敛图Fig.8 Convergence of Aisa network

本文算法虽然在时间上略高于PSO算法，但是MMPSO算法学习到的网络结构得分最高，收敛效果最好，得到了最好的结构.由此证明了MMPSO算法的有效性以及可行性.

5 总结与展望

本文将粒子群优化法和遗传算法的优点进行有效的结合，并将混沌优化策略引入混合算法中形成了一种混合优化算法，同时将这种算法应用到贝叶斯网络结构学习方法中.通过对粒子的变异操作，提高了全局搜索能力，避免陷入局部最优解；将粒子分别与局部最优值和个体最优值交叉，保护了粒子的最优状态值也同时增强了种群的多样性.在提高搜索效率的同时减小了粒子更新的随机性.通过实验结果可以看出，本文在样本容量相同的情况下，学习的网络结构更接近真实网络，并且在最短的有效时间内，学习到分值较高的网络结构，收敛效果最好.但是，在实验中发现，某些明显错误边，特别是在多余边问题上的改善不是很明显，这是算法的不足之处.算法仍需改进优化，将是下一步研究的重点.