一种基于置信传播的算法求解随机约束满足问题

A Belief Propagation Based Algorithm for Solving Stochastic Constraint Satisfaction Problem

摘要: 为了求解具有增长域的随机约束满足问题(CSP)，提出一种基于置信传播的算法即NBP* (new-selected belief propagation*, NBP*)。在置信传播算法中，当BP方程不收敛时，算法就会终止。然而算法在经过多次迭代之后，虽然约束发送给变量的信息没有达到收敛条件，但是仍有部分信息是准确的，所以当算法的BP方程不收敛时，提出利用最后一次迭代得到的约束发送给变量的信息来计算变量的边际概率，当赋值不满足约束时，根据边际概率确定的变量顺序挑选下一个变量进行赋值，得到NBP*算法。数值实验表明：这种算法可以在可满足性相变区域找到解，并且有效提高了置信传播算法的求解效率。

Abstract: In order to solve the stochastic constraint satisfaction problem (CSP) with a growing domain, a belief propagation based algorithm called NBP* (new selected belief propagation*, NBP*) is proposed. In the belief propagation algorithm, when the BP equation does not converge, the algorithm will terminate. However, after multiple iterations, although the information sent by the constraint to the variable did not meet the convergence condition, some information was still accurate. Therefore, when the BP equation of the algorithm does not converge, it is proposed to use the information sent by the constraint to the variable from the last iteration to calculate the marginal probability of the variable. When the assignment does not meet the constraint, the next variable is selected according to the variable order determined by the marginal probability for assignment, resulting in the NBP* algorithm. Numerical experiments have shown that this algorithm can find solutions in the satisfiable phase transition region and effectively improve the efficiency of the confidence propagation algorithm.

1. 绪论

约束满足问题(Constraint Satisfaction Problem, CSP)是人工智能和计算机科学领域一个非常重要的研究对象。在计算复杂性理论中，NP-完全问题都可以表述为CSP，例如子集合加总问题、布尔可满足性问题、顶点涵盖问题、图着色问题、N-puzzle问题等，因此CSP具有非常高的研究价值。

对于CSP的研究，早期主要是在A，B，C，D四种的模型[1] [2] [3]的基础上展开的。但之后Achlioptas等人指出随着问题规模的增大，这四种模型会表现出平凡渐近无解性[4]。为解决这一问题，学者们先后提出了许多改进的模型[4] [5] [6] [7] [8]。2000年，Xu等对B模型进行了修改，提出了RB (Revised B)模型[5]，该模型克服了B模型的平凡渐近无解性，证明了RB模型存在精确的可满足性相变现象，并且RB模型可以在相变区域产生大量难解的实例[9] [10]。

由于RB模型具有精确相变，并且它的取值域会随着问题规模增大而增大，所以受到了国内外学者的广泛关注，成为国际上应用最为广泛的难解实例的产生模型，被广泛应用于ACM等国际算法竞赛。自从相变现象被发现后，统计物理学中的自旋玻璃理论就提供了一个有力的研究工具，在接近相变点时，基于物理学中的无序系统的置信传播(belief propagation, BP)算法展现出惊人的结果[11]，越来越多的学者利用消息传递算法中的置信传播算法来解决CSP。文献[12]中提出两种基于置信传播的信息传递算法求解RB模型产生的随机实例。文献[13]赵春艳等人提出了一种基于变量熵的BP算法。文献[14]中提出在BP方程加入惩罚值的置信传播算法。2019年，文献[15]提出置信传播和模拟退火相结合的算法求解RB模型。文献[16]中，提出一种基于异步更新的置信传播算法求解RB模型。文献[17]中提出一种基于残差的消息传递算法。以上算法都不同程度地提高了置信传播算法(belief propagation algorithm, BP)的性能。

本文基于文献[12]的算法，在置信传播方程迭代不收敛时，提出一种基于置信传播的算法即NBP* (new-selected belief propagation*, NBP*)。若算法迭代收敛，则将具有最大边际概率的变量固定，并令该变量取边际概率最大的分量的值，同时按照边际概率由大到小排列，确定变量的顺序，若迭代不收敛，则利用最后一次迭代信息的结果计算边际概率，当赋值不满足约束时，根据边际概率确定的变量顺序挑选下一个变量进行赋值，并将其固定到边际概率最大的分量上，从而保证算法继续进行。数值实验结果表明，与BP算法相比，NBP*算法可以找到更多的有解实例。

2. 预备知识

2.1. RB模型

CSP定义为一个三元组$\langle X,D,C\rangle$ ，其中X是一组变量，记为$X=\left\{x_1,x_2,x_3,\cdots ,x_N\right\}$ ，D表示值域，记为$D=\left\{d_1,d_2,d_3,\cdots ,d_N\right\}$ ，C表示$t=rN\ln N$ 个约束组成一组约束，记为$C=\left\{C_1,C_2,C_3,\cdots ,C_t\right\}$ ，每个变量$x_i$ 都从它们对应的取值域$d_i$ 中取值，其中$\left|d_i\right|=d\left(i=1,2,3,\cdots ,N\right)$ ，$d=N^{\alpha }$ ，$\alpha >0$ 为常数。RB模型的随机实例$\text{RB}\left(N,k,\alpha ,r,p\right)$ 由以下两个步骤产生：

a) 从$C_N^k$ 中随机可重复地选取t个约束，每个约束$C_a\left(1\le a\le t\right)$ 包含了从变量集合中随机挑选的$k\left(k\ge 2\right)$ 个不同的变量构成的集合$X_a$ 。

b) 对每个约束$C_a\left(1\le a\le t\right)$ ，选定一个相应的不相容赋值集合来限制$X_a$ 中k个变量的某些赋值为不可满足的取值。其中$X_a$ 是从变量集合X中随机挑选的不重复的k个变量组成的；$Q_a$ 是从$X_a$ 中所有可能的$d^k$ 个取值中随机挑选不重复的$p{d}^k$ 个k元赋值，$0<p<1$ 表示约束紧度。

这t个约束就构成了RB模型的一个随机实例。求解这个随机实例就是找到同时满足t个约束的变量的一组赋值。

2.2. RB模型可满足相变现象

文献[5]中利用二阶矩方法严格证明了RB模型存在精确的可满足性相变现象。用$Pr\left(SAT\right)$ 表示RB模型生成的随机实例可满足的概率，则有以下结论成立：

定理1 [5] 设$P_{cr}=1-{\text{e}}^{-\alpha /r}$ ，若$\alpha >1/k$ ，$r>0$ 是两个常数，且$k{\text{e}}^{-\alpha /r}\ge 1$ ，则

${\lim }_{N\to \infty }Pr\left(SAT\right)=\{\begin{array}{ll}1\hfill & p<p_{cr}\hfill \\ 0\hfill & p>p_{cr}\hfill \end{array}$ (1)

定理表明：当N充分大时，随着参数p的增加，当$p>p_{cr}$ 时，RB模型随机实例可满足的概率趋近于；当$p<p_{cr}$ 时，RB模型随机实例可满足的概率趋近于1。由此可见，RB模型随机实例可满足的概率发生了从1降到0，这就是可满足性相变现象，$p_{cr}$ 是相变点。

2.3. RB模型的因子图

由于RB模型在$k\ge 2$ 时都是NP完全的，一般地，我们在生成二元RB模型随机实例时，取$k=2$ 即可。将二元RB模型随机表示成因子图形式，如图1所示。用圆形表示变量节点，记作$i,j,k,\cdots$ ；用正方形表示约束节点，记作$a,b,c,\cdots$ ；若约束$\alpha$ 包含变量i中，则用边连接变量节点i与约束节点$\alpha$ ，记作$\left(\alpha ,i\right)$ 。

每条边上包含了约束a发送给变量i的信息${\eta }_{a\to i}\left(s_i\right)$ 和变量i传递给约束a的信息$u_{i\to a}\left(s_i\right)$ ，${\eta }_{a\to i}\left(s_i\right)$ 表示约束a传递信息让变量i取值为$s_i\left(s_i\in D\right)$ 的概率；$u_{i\to a}\left(s_i\right)$ 表示变量i在没有约束a的信息下取值为$s_i\left(s_i\in D\right)$ 的概率。

2.4. BP迭代方程

根据统计物理学中的空腔场理论，可得到BP迭代方程：

$u_{i\to a}^{\left(t\right)}\left(s_i\right)=\frac{1}{Z^{i\to a}}{\displaystyle \prod _{b\in \delta \left(i\right)\backslash a}{\eta }_{b\to i}^{\left(t\right)}\left(s_i\right)}$ , (2)

${\eta }_{a\to i}^{\left(t+1\right)}\left(s_i\right)=\frac{1}{Z^{a\to i}}{\displaystyle \sum _{s_j\in D,j\in \delta \left(a\right)\backslash i}\delta \left(s_i,s_j\right)u_{j\to a}^{\left(t\right)}\left(s_j\right)}$ , (3)

其中：

$\delta \left(s_i,s_j\right)=\{\begin{array}{l}0,\left(s_i,s_j\right)\in Q_a\\ 1,\left(s_i,s_j\right)\notin Q_a\end{array}$ , (4)

$Z^{i\to a}$ 和$Z^{a\to i}$ 是归一化因子，$\delta \left(i\right)$ 表示与变量i相关的所有约束，$\delta \left(i\right)\backslash a$ 表示与变量i相关的不包括约束a的其余相关的约束；$\delta \left(a\right)$ 表示约束a中包含的所有变量，$\delta \left(a\right)\backslash i$ 表示约束a除去变量包含i的其他变量。

Figure 1. Factor graph of binary RB model

图1. 二元RB模型的因子图

3. NBP*算法

文献[12]中的算法，在BP方程收敛后，计算每个变量取值的边际概率分布，通过先按照最大边际概率挑选变量再利用最大概率分量进行赋值方法执行算法，并且提出的算法表明BP方程会在接近相变点时出现不收敛现象，从而使得算法失效。然而在BP方程迭代到最大迭代次数时，虽然约束传递给变量的信息没有达到收敛条件，但是有一部分信息是准确的，此时根据最后一次迭代的信息计算边际概率分布，对变量进行赋值，再判断赋值是否满足所有约束；当赋值不满足约束时，根据边际概率确定的变量顺序挑选下一个变量进行赋值，从而增加算法收敛的可能性，最后找到满足约束的一组赋值。

在执行算法时，将变量分为三种类型，A类型：已赋值的变量；B类型：与A类型的变量在同一个约束中的变量；C类型：其它变量。

NBP*算法步骤：

输入：二元RB模型一个随机实例的因子图，BP方程的最大迭代次数$t_{\max }$ 和精度$\epsilon$ 。 输出：实例的解或者算法不收敛。 1) $n=1$ ，$t=0$ ，随机初始化每条边上约束发送给变量的信息${\eta }_{a\to i}^{\left(0\right)}\left(s_i\right)$ 。 2) $t=t+1$ ，将${\eta }_{a\to i}^{\left(t-1\right)}\left(s_i\right)$ 代入公式(3)计算得到每条边上变量传递给约束的信息$u_{i\to a}^{\left(t-1\right)}\left(s_i\right)$ (如果$\delta \left(i\right)\backslash a=\varphi$ ，则令$u_{i\to a}^{\left(t-1\right)}\left(s_i\right)=1/d$ )，然后将$u_{i\to a}^{\left(t-1\right)}\left(s_i\right)$ 代入公式(4)更新得到${\eta }_{a\to i}^{\left(t\right)}\left(s_i\right)$ 。 3) 如果$\left|{\eta }_{a\to i}^{\left(t\right)}\left(s_i\right)-{\eta }_{a\to i}^{\left(t-1\right)}\left(s_i\right)\right|<\epsilon$ ，则令${\eta }_{a\to i}^{\ast }\left(s_i\right)={\eta }_{a\to i}^{\left(t\right)}\left(s_i\right)$ ，并执行步骤5)；否则，执行步骤2)。 4) 如果$t=t_{\max }$ ，则令${\eta }_{a\to i}^{\ast }\left(s_i\right)={\eta }_{a\to i}^{\left(t\right)}\left(s_i\right)$ 。 5) 计算每个变量i的取值的边际概率 ${\eta }_i\left(s_i\right)=\frac{{\displaystyle \prod _{a\in \delta i}{\eta }_{a\to i}^{\ast }\left(s_i\right)}}{{\displaystyle \sum _{s_i\in D}{\displaystyle \prod _{a\in \delta i}{\eta }_{a\to i}^{\ast }\left(s_i\right)}}}$ (5) 6) 挑选所有${\eta }_i\left(s_i\right)$ 中的最大值，将其所对应的变量i作为第一个赋值的变量$x_1$ ，对应的取值分量$s_i$ 赋给$x_1$ ，并将每个变量按照${\eta }_i\left(s_i\right)$ 中的边际概率降序排列，再取每个变量的最大边际概率降序排列，得到所有变量的一个顺序$V_1$ 。 7) $n=n+1$ ，$t=0$ ，初始化每条边上约束发送给变量的信息${\eta }_{a\to i}^{\left(0\right)}\left(s_i\right)$ ： 对于A类型的变量：跳过； 对于B类型的变量：如果变量i的取值$s_i$ 与约束中A类型的变量j的取值满足$\left(s_i,s_j\right)\in Q_a$ ，则${\eta }_{a\to i}^{\left(0\right)}\left(s_i\right)=0$ ；否则，${\eta }_{a\to i}^{\left(0\right)}\left(s_i\right)=1$ ； 对于C类型的变量：随机初始化${\eta }_{a\to i}^{\left(0\right)}\left(s_i\right)$ 。 8) $t=t+1$ 对于A类型的变量：跳过； 对于B类型的变量：${\eta }_{a\to i}^{\left(t\right)}\left(s_i\right)={\eta }_{a\to i}^{\left(t-1\right)}\left(s_i\right)$ ； 对于C类型的变量：将${\eta }_{a\to i}^{\left(t-1\right)}\left(s_i\right)$ 代入公式(3)计算得到每条边上变量发送给约束的信息$u_{i\to a}^{\left(t-1\right)}\left(s_i\right)$ (如果$\delta \left(i\right)\backslash a=\varphi$ ，则令$u_{i\to a}^{\left(t-1\right)}\left(s_i\right)=1/d$ )，然后将$u_{i\to a}^{\left(t-1\right)}\left(s_i\right)$ 代入公式(4)更新得到${\eta }_{a\to i}^{\left(t\right)}\left(s_i\right)$ 。 9) 如果$\left|{\eta }_{a\to i}^{\left(t\right)}\left(s_i\right)-{\eta }_{a\to i}^{\left(t-1\right)}\left(s_i\right)\right|<\epsilon$ ，则令${\eta }_{a\to i}^{\ast }\left(s_i\right)={\eta }_{a\to i}^{\left(t\right)}\left(s_i\right)$ ，并执行步骤11)；否则，执行步骤8)。 10) 如果$t=t_{\max }$ ，则令${\eta }_{a\to i}^{\ast }\left(s_i\right)={\eta }_{a\to i}^{\left(t\right)}\left(s_i\right)$ 。 11) 对于B类型和C类型的变量，利用公式(5)计算变量取值的边际概率。 12) 挑选所有${\eta }_i\left(s_i\right)$ 中的最大值，将其所对应的变量i作为第n个赋值的变量$x_n$ ，对应的取值分量$s_i$ 赋给$x_n$ ，并将每个变量按照${\eta }_i\left(s_i\right)$ 中的边际概率降序排列，再取每个变量的最大边际概率降序排列，得到所有变量的一个顺序$V_n$ ，且如果该变量i的赋值不满足约束则执行步骤n)。 13) 如果$n<N$ ，执行步骤g)；如果$n=N$ ，则找到上满足约束的一组赋值，并输出这组赋值作为实例的解。 14) 根据确定变量$x_n$ 时的变量顺序$V_n$ ，如果$V_n$ 有变量未赋值，则将$V_n$ 中当前变量的后一位变量进行赋值，并且将该变量作为第n个待赋值的变量，如果该赋值不满足约束则执行步骤14)；如果该赋值满足约束则执行步骤7)；如果$V_n$ 中无变量可取，则输出没找到实例的解。

4. 数值结果及分析

在数值实验中，本文取$N\in \left\{20,40,60,80,100\right\}$ ，$k=2$ ，$\alpha =0.8$ ，$r=3$ 生成50个二元RB模型的实例。由定理1可以推出二元RB模型的可满足性相变点是$p_s\approx 0.234$ 。对于不同的变量数目N，变量定义域的取值域d和约束的数目t如表1所示。在算法中，取最大迭代次数$t_{\max }={10}^3$ ，精度$ϵ={10}^{-4}$ 。

Table 1. The parameters of binary RB model under different problem scales N

表1. 二元RB模型在不同问题规模N下的参数值

N	α	r	d	t
20	0.8	3	11	180
40	0.8	3	19	443
60	0.8	3	26	737
80	0.8	3	33	1052
100	0.8	3	40	1382

4.1. NBP*算法的结果

在数值实验中，取$N\in \left\{20,40,60,80,100\right\}$ 时，本文在随机生成的50个二元RB模型随机实例上运行NBP*算法，算法找到解的概率如图2所示：

Figure 2. The probability of NBP* algorithm solving on 50 random instances

图2. NBP*算法在50个随机实例上的求解概率

从图2可以看出，在N = 20时，NBP*算法在$p\le 0.17$ 时所有实例都可以找到解，当$p\ge 0.25$ 时算法失效；在N = 100时，NBP*算法在$p\le 0.19$ 时所有实例都可以找到解，当$p\ge 0.22$ 时算法失效。

4.2. 算法分析

NBP*算法和BP算法在50个实例上的求解概率对比分别如图3~7所示。

从图3~7中可以清楚地看出，与BP算法相比，NBP*算法的求解效率有明显效提高。当$0.19\le p\le 0.22$ 内，NBP*算法可以找到更多有解实例，这就说明了当置信传播方程不收敛，仍有部分信息是准确的，可以利用这些信息来提高算法的求解效率。当变量数目$N=\left\{20,40,60,80\right\}$ ，可以通过图3~6看出，当BP算法在接近可满足性相变点时失效，而NBP*算法却仍然可以找到一部分实例的解。

Figure 3. Comparison of NBP* algorithm and BP algorithm results at N = 20

图3. N = 20时NBP*算法和BP算法结果对比

Figure 4. Comparison of NBP* algorithm and BP algorithm results at N = 40

图4. N = 40时NBP*算法和BP算法结果对比

Figure 5. Comparison of NBP* algorithm and BP algorithm results at N = 60

图5. N = 60时NBP*算法和BP算法结果对比

Figure 6. Comparison of NBP* algorithm and BP algorithm results at N = 80

图6. N = 80时NBP*算法和BP算法结果对比

Figure 7. Comparison of NBP* algorithm and BP algorithm results at N = 100

图7. N = 100时NBP*算法和BP算法结果对比

N = 60，NBP*算法在策略影响下求解50个实例时，在$0.18\le p\le 0.23$ 处利用最后一次迭代信息策略找到解和未利用最后一次迭代信息找到解的情况如图8所示。

Figure 8. The situation of finding solutions under the influence of NBP* algorithm strategy when N = 60

图8. N = 60时NBP*算法策略影响下找到解的情况

从图8可以看出，当N = 60时，随着p的增大，在找到解的实例中，NBP*算法利用最后一次迭代的结果找到解的实例与未利用最后一次信息找到解的比值逐渐增大，$p=0.22$ 时，BP算法失效时，NBP*算法利用最后一次迭代结果仍能找到解，这说明了在置信方程不收敛时，最后一次迭代的结果其中有以部分信息是准确的，可以用来寻找实例的解。$p\ge 0.23$ 时，NBP*算法利用策略也无法找到解。

在p = 0.19处，NBP*算法和BP算法在求解二元RB模型的随机实例的运行时间结果对如图9所示。

Figure 9. Comparison of runtime between NBP* algorithm and BP algorithm when p = 0.19

图9. p = 0.19时NBP*算法和BP算法求解的运行时间对比

图9结果表明，算法运行时间随问题规模的增加呈指数型增长。NBP*算法运行时间明显多于BP算法的运行时间。

5. 总结

本文基于置信传播算法，在置信传播方程迭代不收敛时，利用最后一次迭代信息的结果计算边际概率，当赋值不满足约束时，根据边际概率确定的变量顺序挑选下一个变量进行赋值，保证算法继续进行。数值实验结果表明，与BP算法相比，NBP*算法确实提高了随机实例的求解概率，但增加了运行时间。将来的工作可以结合这一思路，进一步优化算法的结果。