基于边际概率分布重新进行单变量选取的置信传播算法求解约束满足问题

Belief Propagation Algorithm for Solving Constraint Satisfaction Problem Based on Marginal Probability Distribution with Re-Selection of Univariate

摘要: 针对RB模型这一类具有增长取值域的随机约束满足问题，提出一种基于边际概率分布重新进行单变量选取的置信传播算法。该算法在置信传播方程不收敛时，通过边际概率分布顺序由大到小找到下一个变量进行重新赋值，从而消去变量的过程。实验结果表明：这种重新挑选变量进行赋值的置信传播算法能在可满足相变区域找到问题的解，有效地提高了置信传播地求解效率。

Abstract: Aiming at the RB model which has a growing value range of stochastic constraint satisfaction problem, a belief propagation algorithm based on marginal probability distribution is proposed. When the belief propagation equation does not converge, the algorithm finds the next variable to be reassigned by the order of marginal probability distribution from large to small, so as to eliminate the process of variables. Experimental results show that the belief propagation algorithm based on reselecting variables for assignment can find the solution of the problem in the region that can satisfy the phase change, and effectively improves the efficiency of belief propagation.

1. 引言

约束满足问题(Constraint Satisfaction Problem, CSP)是人工智能领域的一个重要研究方向。约束满足问题的求解涉及到计算复杂性理论、图论、组合优化等多个学科领域，不仅在计算机科学、统计物理学、离散数学、信息论等领域中都具有非常重要的理论研究价值；而且在视觉处理、资源分配、时序推理、自然语言处理、机器学习等实际问题中也得到了广泛关注和研究。

CSP涉及一组变量和这些变量需要满足的一组约束条件。每个变量都有一个非空的可能取值域，每个约束描述了一个变量子集以及子集之间各变量的不相容赋值集合，CSP的目标是在满足所有约束条件的情况下，找到一组变量值(解)。N皇后、图着色、合取式范式的可满足性(satisfiability, SAT)和最小顶问题点覆盖等都是CSP。一般情况下，随机CSP是NP完全问题。为了研究CSP，早期的许多理论研究和算法实验基本都是围绕四个标准的二元CSP模型A、B、C、D展开的 [1] [2] [3] ，Achlioptas [4] 等人指出这四个模型具有平凡渐近不可满足性的缺点，为了克服这一缺点，Molloy [5] [6] 等人先后提出了一些改进模型，这些模型在理论上能产生非平凡的难解实例，但是它们约束的产生并不是简单而自然的。2000年，许可和李未提出了RB模型 [7] (revised B)，RB模型是一个典型的具有增长取值域的随机约束满足问题，在RB模型中，变量的取值域随问题规模的增大呈多项式级增长，这与以往具有固定增长域的随机约束满足问题存在不同的特点，克服了B模型不能产生难解实例的缺点，引起了国内外学者的广泛关注。许可等人 [7] 利用二阶矩方法严格证明RB模型不仅存在可满足性相变现象，而且可以得到精确的相变点。文献 [8] [9] 分别从理论和实验上说明了求解RB模型随机实例的难度随着问题规模的增加呈多项式级增长，证明了RB模型在相变点附近可以产生大量的难解实例。因此，RB模型是一个非常值得研究的具有精确相变现象并且能从相变区域产生大量难解实例的NP完全问题。

RB模型可满足相变现象已经被证明，但求解RB模型的随机实例仍然具有挑战性。2011年，文献赵春艳等人引入统计物理学理论中的空腔场方法，提出了由置信传播算法引导的信息传播(Belief Propagation，简称BP)算法求解具有增长取值域的约束满足问题 [10] ，这些算法在接近理论相变点时都能有效找到实例的解。赵春艳等人提出了一种基于变量熵来挑选变量的置信算法，且进一步探讨了RB模型的解空间 [11] 。王晓峰等人从理论上分析了置信传播算法在RB模型实例上的收敛性 [12] 。原志强等人提出了两种改进的模拟退火(Revised Simulated Annealing，简称RSA)算法求解增长值域的约束满足问题 [13] 。吴拨荣等人将置信传播和模拟退火相结合的算法求解RB模型 [14] 。李飞龙等人提出了基于禁忌搜索(Tabu Search，简称TS)算法求解RB模型并将禁忌搜索算法与置信传播算法相结合 [15] 。范如梦等人提出了基于度启发式和最少约束值启发式的回溯算法(Heuristic Backtracking Algorithm, HBT)，通过度启发式来确定待赋值变量顺序，然后利用最少约束值启发式对选择的变量进行赋值，来求解RB模型 [16] 。赵春艳等人提出了不同约束紧度下p-RB模型存在精确的可满足性相变现象 [17] 。2022年，张学才等人提出两种基于动态度的回溯算法 [18] 求解RB模型，与经典回溯算法相比，具有显著的优越性。2023年，林童提出基于回溯的置信满足算法 [19] ，实验结果表示该算法能够更有效地找到问题的解。大量算法研究表明，如何挑选变量和如何对变量赋值是解决问题的关键。

本文基于文献 [10] 中的算法，提出了根据边际概率分布重新进行单变量选取的置信传播算法(New-selected belief propagation, NBP)来求解RB模型这一类具有可变定义域的随机CSP。在NBP算法中，若BP方程收敛，则利用变量的边际概率分布选取最大边际概率的变量作为待赋值的变量，并令该变量取其最大边际概率分量的值。同时将待赋值变量根据边际概率分布进行降序排列，以此作为挑选变量的顺序；若迭代不收敛，则重新挑选变量，选择当前变量的后一位待赋值变量，并将其值固定到该变量对应的最大边际概率的分量上，再重新进行迭代。数值实验结果表明：与BP算法相比，NBP算法通过重新挑选待赋值变量可以找到更多实例的解，提高了随机实例的求解概率。

2. RB模型

RB模型是由N个变量组成的变量集合 $X=\left\{x_1,x_2,x_3,\cdots ,x_N\right\}$ 和t个约束组成的约束集合 $C=\left\{C_1,C_2,C_3,\cdots ,C_t\right\}$ 构成的，每个变量都从它们对应的定义域 $D=\left\{d_1,d_2,d_3,\cdots ,d_N\right\}$ 中取值，其中 $\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)$ 由以下步骤产生：

从 $C_N^k$ 中随机可重复地选取t个约束，其中 $t=rN\ln N$ ， $r>0$ 为常数。每个约束 $C_a\left(1\le a\le t\right)$ 包含了从变量集合中随机挑选的 $k\left(k\ge 2\right)$ 个不同的变量构成的集合集合 $X_a$ 和一个相应的不相容赋值集合来限制 $X_a$ 中k个变量的某些赋值为不可满足的取值。其中 $X_a$ 是从变量集合X中随机挑选的不重复的k个变量组成的； $Q_a$ 是从 $X_a$ 中所有可能的 $d^k$ 个取值中随机挑选不重复的 $p{d}^k$ 个k元赋值， $0<p<1$ 表示约束紧度。这t个约束就构成了RB模型的一个随机实例。求解这个随机实例就是至少找到一个解，也就是找到同时满足t个约束的N个变量的一组赋值。

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

定理1 [7] 设 $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}{c}1p<p_{cr}\\ 0p>p_{cr}\end{array}$ (1)

定理2 [7] 设 $r_{cr}=-\alpha /\ln \left(1-p\right)$ ，若 $\alpha >1/k$ ， $0<p<1$ 是两个常数，且 $k\ge 1/\left(1-p\right)$ ，则

${\lim }_{N\to \infty }Pr\left(SAT\right)=\{\begin{array}{c}1r<r_{cr}\\ 0r>r_{cr}\end{array}$ (2)

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

3. 根据边际概率分布重新进行单变量选取的置信传播算法

3.1. RB模型的因子图和BP迭代方程

RB模型在 $k\ge 2$ 时都是NP完全的，因此我们取 $k=2$ 。将二元RB模型随机表示成因子图形式，如图1所示。圆圈表示变量节点，记作 $i,j,k,\cdots$ ；方框表示约束节点，记作 $a,b,c,\cdots$ ；若变量i在约束α中，则变量节点i与约束节点α之间用边连接，记作 $\left(\alpha ,i\right)$ 。由于RB模型约束数为 $O\left(N\ln N\right)$ ，所以因子图具有局部树状结构。

Figure 1. Factor graph of binary RB model

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

每条边上定义两种信息 ${\eta }_{a\to i}\left(s_i\right)$ 和 $u_{i\to a}\left(s_i\right)$ ， ${\eta }_{a\to i}\left(s_i\right)$ 是约束a发送给变量i的信息，表示约束a传递信息给变量i令其取值为 $s_i\left(s_i\in D\right)$ 的概率； $u_{i\to a}\left(s_i\right)$ 是变量i传递给约束a的信息，表示变量i在没有约束a的信息下取值为 $s_i\left(s_i\in D\right)$ 的概率。根据统计物理学中的空腔场理论，可得到BP迭代方程：

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

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

其中：

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

这里 $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除去变量包含的其他变量。

3.2. NBP算法

事实上，文献 [10] 中提出的第二种算法是利用BP方程的结果，每次挑选一个变量进行赋值，逐步消去变量，但是BP方程会在接近相变点时出现不收敛的现象，从而使得算法失效，为了提高算法的求解效率，NBP算法根据迭代方程计算边际概率，将概率最大的边际概率所对应的变量固定，令该变量取其最大边际概率分量的值，并按照每个变量的最大边际概率大小对变量进行排序，若迭代收敛，算法运行；若算法迭代不收敛或赋值不满足约束，将按照变量顺序对下一个变量进行赋值，将其固定到最大边际概率分量上，若变量顺序中无变量可挑选，则结束算法，通过边际概率重新挑选单变量进行赋值，来提高算法收敛的可能性，从而找到满足约束的一组赋值。

置信传播算法是在BP方程收敛后，计算每个变量取值的边际概率分布，通过先按照最大边际概率挑选变量再利用最大概率分量进行赋值方法执行算法，在执行算法中，将变量分为三类，A类型：已赋值的变量；B类型：与A类型的变量在同一个约束中的变量；C类型：其它变量。

NBP算法步骤：

输入：二元RB模型一个随机实例的因子图，BP方程的最大迭代次数 $t_{\max }$ 和精度 $\epsilon$ 。

输出：实例的解或者算法不收敛。

a) $n=1$ ， $t=0$ ，随机初始化每条边上约束发送给变量的信息 ${\eta }_{a\to i}^{(0)}\left(s_i\right)$ 。

b) $t=t+1$ ，将 ${\eta }_{a\to i}^{(t-1)}\left(s_i\right)$ 代入公式(3)计算得到每条边上变量传递给约束的信息 $u_{i\to a}^{(t-1)}\left(s_i\right)$ (如果 $\delta \left(i\right)\backslash a=\varphi$ ，则令 $u_{i\to a}^{(t-1)}\left(s_i\right)=1/d$ )，然后将 $u_{i\to a}^{(t-1)}\left(s_i\right)$ 代入公式(4)更新得到 ${\eta }_{a\to i}^{(t)}\left(s_i\right)$ 。

c) 如果 $\left|{\eta }_{a\to i}^{(t)}\left(s_i\right)-{\eta }_{a\to i}^{(t-1)}\left(s_i\right)\right|<\epsilon$ ，则令 ${\eta }_{a\to i}^{\ast }\left(s_i\right)={\eta }_{a\to i}^{(t)}\left(s_i\right)$ ，并执行步骤e)；否则，执行步骤b)。

d) 如果 $t=t_{\max }$ ，则输出不收敛。

e) 计算每个变量i的取值的边际概率

${\eta }_i\left(s_i\right)=\frac{{\displaystyle \underset{a\in \partial i}{\prod }{\eta }_{a\to i}^{\ast }\left(s_i\right)}}{{\displaystyle \underset{s_i\in D}{\sum }{\displaystyle \underset{a\in \partial i}{\prod }{\eta }_{a\to i}^{\ast }\left(s_i\right)}}},$ (6)

f) 挑选所有 ${\eta }_i\left(s_i\right)$ 中的最大值，将其所对应的变量i作为第一个赋值的变量 $x_1$ ，对应的取值分量 $s_i$ 赋给 $x_1$ ，并将变量按照 ${\eta }_i\left(s_i\right)$ 中最大边际概率降序排列，得到所有变量的一个变量顺序V₁。

g) $n=n+1$ ， $t=0$ ，初始化每条边上约束发送给变量的信息 ${\eta }_{a\to i}^{(0)}\left(s_i\right)$ ：

对于A类型的变量：跳过；

对于B类型的变量：如果变量i的取值 $s_i$ 与约束中A类型的变量j的取值满足 $\left(s_i,s_j\right)\in Q_a$ ，则 ${\eta }_{a\to i}^{(0)}\left(s_i\right)=0$ ；否则， ${\eta }_{a\to i}^{(0)}\left(s_i\right)=1$ ；

对于C类型的变量：随机初始化 ${\eta }_{a\to i}^{(0)}\left(s_i\right)$ 。

h) $t=t+1$

对于A类型的变量：跳过；

对于B类型的变量： ${\eta }_{a\to i}^{(t)}\left(s_i\right)={\eta }_{a\to i}^{(t-1)}\left(s_i\right)$ ；

对于C类型的变量：将 ${\eta }_{a\to i}^{(t-1)}\left(s_i\right)$ 代入公式(3)计算得到每条边上变量发送给约束的信息 $u_{i\to a}^{(t-1)}\left(s_i\right)$ (如果 $\delta \left(i\right)\backslash a=\varphi$ ，则令 $u_{i\to a}^{(t-1)}\left(s_i\right)=1/d$ )，然后将 $u_{i\to a}^{(t-1)}\left(s_i\right)$ 代入公式(4)更新得到 ${\eta }_{a\to i}^{(t)}\left(s_i\right)$ 。

i) 如果 $\left|{\eta }_{a\to i}^{(t)}\left(s_i\right)-{\eta }_{a\to i}^{(t-1)}\left(s_i\right)\right|<\epsilon$ ，则令 ${\eta }_{a\to i}^{\ast }\left(s_i\right)={\eta }_{a\to i}^{(t)}\left(s_i\right)$ ，并执行步骤k)；否则，执行步骤h)。

j) 如果 $t=t_{\max }$ ，则BP方程不收敛，执行步骤n)。

k) 对于B类型和C类型的变量，利用公式(6)计算变量取值的边际概率。

l) 挑选所有 ${\eta }_i\left(s_i\right)$ 中的最大值，将其所对应的变量i作为第n个赋值的变量 $x_n$ ，对应的取值分量 $s_i$ 赋给 $x_n$ ，并将 ${\eta }_i\left(s_i\right)$ 中所有边际概率降序排列，再按照每个变量的最大边际概率排列，得到此时的一个变量顺序 $V_n$ ，且如果该变量的赋值不满足约束则执行步骤n)；

m) 如果 $n<N$ ，执行步骤g)；如果 $n=N$ ，则找到上满足约束的一组赋值，并输出这组赋值作为实例的解。

n) 根据变量 $x_n$ 的变量顺序 $V_n$ ，如果 $V_n$ 中在变量 $x_n$ 后还有变量可取，则将 $V_n$ 中当前赋值变量的后一位变量进行赋值，若赋值满足约束，则该变量成为第n个赋值变量，并执行步骤g)；如果该赋值不满足约束则执行步骤n)；如果 $V_n$ 无变量可取，则输出未找到实例的解。

4. 数值结果及分析

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

Table 1. Parameters of the binary RB model corresponding to different numbers of variables N

表1. 二元RB模型对应不同的变量数目N时的参数

4.1. 算法结果

本文在随机生成的50个二元RB模型随机实例上运行NBP算法，算法找到解的概率如图2所示：

Figure 2. Probability of NBP algorithm solving on 50 random instances

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

从图2可以看出，在N = 20时，NBP算法在 $p\le 0.18$ 时找到解的概率是100%，当 $p\ge 0.25$ 时找不到任何解；在N = 80时，NBP算法在 $p\le 0.19$ 时找到解的概率是100%，当 $p\ge 0.22$ 时找不到任何解。

4.2. 算法分析

$N=\left\{20,40,60,80\right\}$ 时，NBP算法和BP算法在50个实例上的求解概率对比分别如图3~6所示。从图中可以看出，与BP算法相比，NBP算法可以明显提高对实例的求解概率，尤其是在区域 $0.19\le p\le 0.22$ 内，NBP算法可以有效地找到更多有解实例。

在p = 0.19和处，NBP算法和BP算法在求解二元RB模型的随机实例的运行结果对如图7所示。结果表明，算法运行时间随变量数目的增加呈指数型增长。由于NBP算法重新挑选变量进行赋值，则其运行时间明显多于BP算法的运行时间。

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

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

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

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

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

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

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

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

Figure 7. Comparison of runtime between NBP algorithm and BP algorithm for solving at p = 0.19

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

5. 总结

本文提出了根据边际概率分布重新进行单变量选取的置信传播算法来求解具有可变定义域的随机约束满足问题。在NBP算法中，若BP方程迭代不收敛，则重新挑选变量，选择当前变量的后一位待赋值变量，并将其值固定到该变量对应的最大边际概率的分量上，再重新进行迭代，保证算法继续进行。数值实验结果表明，与BP算法相比，NBP算法通过重新挑选待赋值变量可以找到更多实例的解，提高了随机实例的求解概率。但因为重新挑选变量赋值，增加了算法的运行时间，在以后的研究中可以通过进一步改进算法，提高算法的运行速度。