1. Introduction
The difficulty in understanding complex systems manifests not only in finding a proper description of their interaction topology, but also in quantifying the collective and cooperative behaviors among their constituents at multiple scales [1]. Message passing algorithms prove to be an efficient tool in tackling hard computational problems formulated in graphical models [2], which is a simple description of complex interacted systems. In a typical algorithmic procedure, one defines messages on the connections of basic units in a graphical representation, derives iterative equations of messages to account for the propagation of correlations among variables, and extracts information on problem solutions through messages after a sufficiently large number of updating steps. Among the main roots of the message passing algorithms are coding theory in information theory research [3] and statistical physics of disordered systems [4]. Message passing algorithms have been heavily explored in optimization problems [5, 6], inference problems [7], learning problems [8, 9], and so on. Different frameworks of message passing algorithms are based on mean-field theory at a certain level of approximation, typical examples of which include the belief propagation (BP) algorithm [10], the survey propagation algorithm [11], the cluster variational algorithm [12, 13], and so on.
Constraint satisfaction problems (CSPs) [5] are among the most typical computational tasks in complex systems research. A typical instance of CSPs can be formulated as a bipartite graph with two types of nodes denoting m constraints and n variables respectively and edges only between constraints and variables. In a random setting, edges can be established from a null graph with only nodes by connecting a given constraint node with multiple variable nodes chosen among the variable node set with a uniform probability. A basic computational question for CSPs is to find the existence and the number of assignments to variables which satisfy all the constraints in a given instance. CSPs are heavily discussed in computational complexity theory in theoretical computer science [5], as many prototype NP-hard problems can be formulated as CSPs. In the language of statistical physics, satisfiability and optimization problems can be understood as the energy landscapes of configurations of a physical system with continuous or discrete states, while their global optima correspond to the ground-state energy of landscapes. Thus we can adopt statistical mechanical frameworks, such as the mean-field theory of spin glasses [4], to estimate the energy and the entropy of their ground-state solutions. Typical CSPs studied with the spin glass theory are the random k-satisfiability (k-SAT) problem [14, 15], the minimum vertex cover problem [16, 17], and the q-coloring problem [18]. These problems focus on random CSPs with fixed domains, whose size of variables is independent of the variable number n (that is, 2, 2, and q as given integers, respectively). However, many practical problems, such as the Latin square problem, the n-queens problem and so on, can be transformed into random CSPs with growing domains. Among the models for these CSPs, model RB (revised B) proposed in [19] is a popular one, whose domain size grows polynomially with the variable number n.
Spin glass mean-field theory on satisfiability problems shows that fundamental structural transitions exist in the geometric organization of ground-state configurations. Taking the random k-SAT (k ≥ 3) problem [14] as an example, with an increase in constraint densities (ratio of constraint number to variable number), there are clustering, condensation, and SAT/UNSAT phase transitions separated by multiple thresholds, where drastic changes in the numbers and sizes of the clusters of ground-state configurations take place. Furthermore, there is a deep connection between these solution space phase transitions and the algorithmic performance in finding solutions. Specifically, around these thresholds, messages in mean-field theory often show large fluctuations, rendering impossible the convergence of messages from which a typical message passing algorithm can extract information. Meanwhile, each mean-field theory has an approximation at some level in the organization of the solution space. For example, for the BP algorithm [10], solutions are assumed to be organized in a single cluster, or a pure state. When a level of approximation fails, there are often two choices: whether we generalize the current mean-field assumption to a theoretical framework with a more complicated and computationally demanding form, or we introduce modifications at an algorithmic level to achieve better performance by leveraging the existing form of message passing algorithm.
Our focus in this paper is to suppress the fluctuation of iterative messages around thresholds, a critical problem for message-passing-based algorithms, thus to improve the performance of existing BP algorithm for solving CSPs at the algorithmic level. Our main contribution here is to define residuals in the message updating process, based on which the existing BP algorithm is further modified. We show that, for a specific CSP model RB, compared with its usual BP algorithm, our modification of the algorithm can significantly improve the convergence of the updated messages and increase the probability of finding ground-state solutions around the satisfiability thresholds at a lower computational cost.
The layout of the paper is as follows. In Section 2 , we briefly introduce the concept and some known results of model RB, which is a typical CSP with growing domains. Section 3 presents the maximal residual BP (MRBP) algorithm, which modifies the current BP algorithm based on residuals in the process of updating messages. In Section 4 , we illustrate and analyze the numerical results of our algorithm. In Section 5 , we conclude the paper with some discussions.
2. Model RB: an example of constraint satisfaction problems
Over the past few decades, studies on CSPs have mainly focused on the initial standard CSP models, named A, B, C and D [20–22]. However, it has been proven by Achlioptas et al that the random instances generated by the four standard models have trivial asymptotic insolubility with an increase in problem size [23]. Therefore, in order to overcome this defect, various improved models from different perspectives in numerous efforts have been proposed successively [19, 24–27]. Model RB, proposed by Xu and Li to avoid the trivial unsatisfiability of model B [19], is a modification on the domain size of variables and the number of constraints of model B without special restrictions on the constraint structure.
Model RB consists of a set of n variables X = {x1, x2, … , xn} and a set of m constraints C = {C1, C2, … , Cm} with $m={rn}\mathrm{ln}n$, where r > 0 is a constant. Each variable xi (i = 1, 2, … , n) has a nonempty domain D = {1, 2, … , dn} of dn = nα (α > 0 is a constant) possible values. Each constraint Ca (a = 1, 2, … , m) involves k ( ≥ 2) different variables, and there is a set of incompatible assignments Qa ⊂ Dk which specifies the disallowable combinations of values for the k variables. More precisely, we generate a random instance of model RB with the following two steps.
A random instance of model RB is simply the conjunction of the above m randomly selected constraints.
| • | Step 1. To construct a single constraint, we randomly select k ( ≥ 2) distinct variables out of the n variables, and then randomly select exactly pdnk (p is the constraint tightness) different incompatible assignments out of dnk possible ones as the set Qa to restrict the values of these k variables. |
| • | Step 2. We repeat the above step to obtain $m={rn}\mathrm{ln}n$ independent constraints and their corresponding sets Qa (a = 1, 2, ⋯ ,m). |
Let Pr(SAT) be the probability that a random instance is satisfied. It is shown in [19] that
Let ${p}_{s}=1-{{\rm{e}}}^{-\tfrac{\alpha }{r}}$. If $\alpha \gt \tfrac{1}{k}$, $r\gt 0$ are two constants and $k{{\rm{e}}}^{-\tfrac{\alpha }{r}}\geqslant 1$, then
$\begin{eqnarray}\mathop{\mathrm{lim}}\limits_{n\to \infty }\Pr (\mathrm{SAT})=\left\{\begin{array}{lll}1, & {\rm{if}} & p\lt {p}_{s},\\ 0, & {\rm{if}} & p\gt {p}_{s}.\end{array}\right.\end{eqnarray}$
Let ${r}_{s}=-\tfrac{\alpha }{\mathrm{ln}(1-p)}$. If $\alpha \gt \tfrac{1}{k}$, $0\lt p\lt 1$ are two constants and $k\geqslant \tfrac{1}{1-p}$, then
$\begin{eqnarray}\mathop{\mathrm{lim}}\limits_{n\to \infty }\Pr (\mathrm{SAT})=\left\{\begin{array}{lll}1, & {\rm{if}} & r\lt {r}_{s},\\ 0, & {\rm{if}} & r\gt {r}_{s}.\end{array}\right.\end{eqnarray}$
The two theorems show that model RB exhibits a sharp threshold of satisfiability at the critical values ps and rs, as p and r are two control parameters of the problem. Nevertheless, we cannot rigorously determine the exact value of the satisfiability threshold for some random CSPs with fixed domains [28, 29]. It has been theoretically proven that almost all instances of model RB have no tree-like resolution proofs of less than exponential size [30]. In other words, these random instances of model RB in the transition region are extremely difficult to solve. Numerical results, using complete and incomplete search algorithms, have confirmed the exact satisfiability thresholds and the hardness of forced and unforced satisfiable instances [31]. Although model RB shows an asymptotic phase transition phenomenon, it is still very challenging to find solutions of a random instance, and the size of the problem has been restricted to n ∼ 100 [32]. Therefore, benchmarks based on model RB are widely used in various kinds of algorithm competitions such as CSP, SAT, pseudo-Boolean and so on, and these results have verified the intrinsic computational hardness of these benchmarks.
Both analytical and algorithmic studies have been carried out on model RB. On the analytical side, Fan and Shen defined a large class of random CSP models d-k-CSP which unified several related models including model RB [33, 34]. Later, it was proven that the restriction on k can be weaker, which can simplify the generation of random instances [35]. In order to study the transition behaviors of model RB due to finite-size effects, Zhao et al used finite-size scaling analysis to bound the width of the transition region [36]. On the algorithmic side, statistical physics concepts and methods, especially the replica method and the cavity method in spin glass mean-field theory, have played a very significant role in constructing solutions for random CSPs [37, 38]. Zhao et al proposed two different types of message passing algorithms guided by BP to solve random CSPs with large domains like model RB [39]. Subsequently, Zhao et al proposed a reinforced BP algorithm, which can effectively force the BP equations to converge to a solution of model RB by introducing an external field with a certain probability [40]. Moreover, it was shown that the solutions of model RB are grouped in disconnected clusters before the theoretical satisfiability phase transition threshold [40]. Recently, Zhao et al proposed three kinds of BP decimation algorithms, which improved the performance of BP algorithms by improving the updating method of the BP iterative equation [41].
3. Maximal residual belief propagation algorithm
Here we propose the MRBP algorithm, which combines modifications in the updating process based on residuals with the usual BP algorithm. We first lay down the general procedure of the BP algorithm and the relevant equations, and then we define the residuals of messages in the updating process and present our algorithm in detail.
For k ≥ 2, model RB is NP-complete. Therefore, we take the binary case (k = 2) of model RB as the tested problem. Model RB with k ≥ 3 can be reduced to the case of k = 2 in polynomial time.
3.1. Belief propagation algorithm
An instance of binary model RB admits a natural factor graph representation, as shown in figure 1. A factor graph is a bipartite undirected graph, in which one type of nodes represents the variables (denoted by i, j, ⋯ ) and the other types of nodes represent the constraints (denoted by a, b, ⋯ ). An edge (a,i) connects a constraint node a and a variable node i if and only if a contains i. The average degree of each variable node in the graph is $r\mathrm{ln}n$. Locally such a graph is tree-like; the typical size of a loop in the graph scales like $\mathrm{ln}n/\mathrm{lnln}n$ for large n.
Figure 1. Part of a factor graph representing a random instance of binary model RB. Circles and squares denote the variable nodes and the constraint nodes, respectively. In the message passing algorithm of model RB, μa→i(σi) and ηi→a(σi) are the messages passed between a and i along the edge (a,i) in opposite directions. |
An instance of model RB can be treated as a statistical mechanical system with interacting variables with discrete states and an interaction topology in the form of a factor graph. In a random instance of binary model RB, we let $\vec{\sigma }=({\sigma }_{1},{\sigma }_{2},\ldots ,{\sigma }_{n})$ (σi ∈ D, i = 1, 2, ⋯ ,n) be an assignment to the n variables. We further let $\vec{{\sigma }_{a}}=({\sigma }_{a1},{\sigma }_{a2})$ be an assignment to the two variables involved in the constraint a ( ∈ {1, 2, ⋯ ,m}). We define the energy of the constraint a as ${f}_{a}(\vec{{\sigma }_{a}})=1$ if $\vec{{\sigma }_{a}}\in {Q}_{a}$ and ${f}_{a}(\vec{{\sigma }_{a}})=0$ if $\vec{{\sigma }_{a}}\notin {Q}_{a}$. The total energy of an assignment $\vec{\sigma }$, denoted as $E(\vec{\sigma })$, is the sum of energy terms on all the m constraints, which is the number of unsatisfied constraints computed following
$\begin{eqnarray}E(\vec{\sigma })=\sum _{a=1}^{m}{f}_{a}(\vec{{\sigma }_{a}}).\end{eqnarray}$
Based on the locally tree-like property of the factor graph representation of model RB, we derive the framework of the BP algorithm [4]. Associated with each edge (a, i) there are two cavity messages μa→i(σi) and ηi→a(σi) passing in opposite directions. See figure 1 for an example. Specifically, μa→i(σi) is the probability that the constraint a is satisfied if the variable i takes the value σi, and ηi→a(σi) is the probability that the variable i takes the value σi in the absence of the constraint a. Let $\{{\mu }_{a\to i}^{t}({\sigma }_{i})\}$ and $\{{\eta }_{i\to a}^{t}({\sigma }_{i})\}$ denote the set of messages passed along all the edges in the factor graph at time t = ( 0, 1, ⋯ ). According to the cavity method in disordered systems of statistical physics [4, 39], the self-consistent BP equations of cavity messages can be derived as
$\begin{eqnarray}{\eta }_{i\to a}^{t}({\sigma }_{i})=\displaystyle \frac{1}{{Z}_{i\to a}}\prod _{b\in V(i)\setminus a}{\mu }_{b\to i}^{t}({\sigma }_{i}),\end{eqnarray}$
$\begin{eqnarray}{\mu }_{a\to i}^{t+1}({\sigma }_{i})=\displaystyle \frac{1}{{Z}_{a\to i}}\sum _{j\in V(a)\setminus i,{\sigma }_{j}\in D}{f}_{a}(\vec{{\sigma }_{a}}){\eta }_{j\to a}^{t}({\sigma }_{j}).\end{eqnarray}$
Here Zi→a and Za→i are both normalization factors, V(i)\a the set of constraints adjacent to the variable i removing a, and V(a)\i the set of variables adjacent to the constraint a removing i.In a typical procedure of the BP algorithm, messages on the edges in the factor graph of a given problem instance are uniformly initialized at random in [0, 1]. A message iteration process is carried out, in which messages are updated following the self-consistent equations (4 ) and (5 ). After a sufficiently long iterative process, the messages may converge to a fixed point. The convergence criterion can be measured by a precision parameter ϵ ( ≪ 1): when the maximal absolute difference between two consecutive steps of the messages on all edges of a factor graph is smaller than ϵ, we can say that the message updating process converges. With the fixed point of messages, the statistical mechanical properties of model RB in an average sense can be extracted, such as the energy and entropy of ground-state solutions, and so on.
To find a solution configuration of a given problem instance, we adopt the BP decimation (BPD) algorithm, which assigns values to variables one by one based on the marginal probability after message updating. In order to find a solution configuration of the problem, we can assign the value of each variable based on its marginal probability distribution. This is the basic procedure of the BPD algorithm. With the fixed point of messages $\{{\mu }_{a\to i}^{* }({\sigma }_{i})\}$ for all edges (a, i), the marginal probability of a variable σi (i = 1, 2, ⋯ ,n) is computed following
$\begin{eqnarray}{\mu }_{i}({\sigma }_{i})=\displaystyle \frac{\prod _{a\in V(i)}{\mu }_{a\to i}^{* }({\sigma }_{i})}{\sum _{{\sigma }_{i}\in D}\prod _{a\in V(i)}{\mu }_{a\to i}^{* }({\sigma }_{i})}.\end{eqnarray}$
In a decimation step, the variable with the largest marginal probability is selected to be assigned to its corresponding most biased component. In the decimation procedure, we iteratively switch between a variable fixing process and a message updating process after a graph simplification process upon fixing variables.3.2. Maximal residual belief propagation algorithm
The convergence of messages in the BP iterative equations affects the performance of the algorithm. However, the BP equations usually do not converge when the control parameter approaches the theoretical satisfiability threshold. The core idea of this work is the introduction of the residuals for messages. We modify the message updating process, as we dynamically update the passing messages between the constraints and the variables based on the maximal residual in BP iterative equations, to improve the convergence and the performance of BP algorithm. The residuals of messages measure the relative difference sent from constraints to variables between two consecutive steps during the iteration of messages. If we consider all the messages in the BP algorithm as a time series are $\{{\mu }_{a\to i}^{t}\}$ with t = 0, 1, ⋯ , the residual on an edge (a,i) at time t is defined as
$\begin{eqnarray}{\gamma }_{a\to i}^{t}=\displaystyle \frac{| {\mu }_{a\to i}^{t}-{\mu }_{a\to i}^{t-1}| }{{\mu }_{a\to i}^{t}}.\end{eqnarray}$
A large residual indicates a large fluctuation of the message ${\mu }_{a\to i}^{t}$, thus the messages at these steps are far from the value after convergence of the whole message passing algorithm, if there is one.Based on the residual on each edge at a time step, the MRBP algorithm dynamically adjusts the order of updating messages according to the real-time status of the messages, in which the messages connected to the edge with the largest residual are prioritized to update. See figure 2 for an example. At time t, we select the edge (aγ, iγ) with the largest residual. We first adopt equation (4 ) to compute the messages sent from the variables iγ and jγ associated with the constraint node aγ to the neighboring constraint nodes connected to iγ and jγ except aγ (denoted by V(iγ) ∪ V(jγ)\aγ), as indexed by the number 1 in figure 2. We need to reduce the greediness of the algorithm to achieve better performance. Hence, the updating messages are considered from the point of view of node aγ rather than edge (aγ, iγ). Therefore, we use equation (4 ) to calculate the messages sent from the nodes iγ and jγ connected to aγ to V(iγ) ∪ V(jγ)\aγ simultaneously. Then, we follow equation (5 ) to update the messages that the constraint nodes V(iγ) ∪ V(jγ)\aγ send to the variables connected to the constraints except iγ and jγ, as indexed by the number 2 in figure 2. In order to avoid falling into the local optimum in the iteration process, we mark the edge with the largest residual after the message updating, and the residuals on unmarked edges are recalculated. The selection of the current edge with the largest residual and the recalculation of residuals on unmarked edges are carried out iteratively until all edges are marked to ensure that messages on all edges are updated during the BP iteration process. In this way, the messages between constraints and variables are constantly updated till they converge to a fixed point or the maximal number of iterations is reached.
Figure 2. Local procedure of message updating based on maximal residuals on a factor graph. The edge (aγ, iγ) with the largest residual is selected for updating and then marked with 0, and the numbers 1 and 2 represent the order of message updating. |
During the decimation process of the BPD algorithm, we iteratively fix variables with values until we have an assignment for all variables. Once a variable in a factor graph is assigned a certain value, variables in the factor graph are divided into three different types, as illustrated in figure 3.
| • | Type A: The variables that have been assigned with some values. |
| • | Type B: The variables not assigned with values, yet connected to the same constraints as the assigned ones. |
| • | Type C: Others. |
Figure 3. In the decimation process of the BPD algorithm, we classify the variables into three different types, indicated by A, B and C. |
The MRBP algorithm is detailed in the following pseudocode.
MRBP algorithm
| Input: A factor graph of a random instance generated by model RB, maximal iterations ${t}_{\max }$ used in the subroutines MRBP-UPDATE and MRBP-UPDATE*, and a precision parameter ϵ. |
| Output: A solution or ‘UNCONVERGED’. |
| Step 1 At T = 0: |
| 1.1 Run the subroutine MRBP-UPDATE. |
| 1.2 For each variable i, compute the marginal probability μi(σi) by equation ( |
| 1.3 Select the variable i* with the highest marginal probability from the n variables, and assign it to its most biased value ${\sigma }_{{i}^{* }}^{* }=\mathrm{argmax}{\mu }_{{i}^{* }}({\sigma }_{{i}^{* }})$. |
| Step 2 From T = 1 to n − 1: |
| 2.1 Run the subroutine MRBP-UPDATE*. |
| 2.2For each free (unfixed) variable i, compute the marginal probability μi(σi) by equation ( |
| 2.3 Select the variable i* with the highest marginal probability from the (n − T) free variables, and assign it to its most biased value ${\sigma }_{{i}^{* }}^{* }=\mathrm{argmax}{\mu }_{{i}^{* }}({\sigma }_{{i}^{* }})$. |
| 2.4 Check the energy function $E(\vec{\sigma })$ for the variables with assigned values according to equation ( |
| Step 3 If the assignment ${\vec{\sigma }}^{* }=({\sigma }_{1}^{* },{\sigma }_{2}^{* },\ldots ,{\sigma }_{n}^{* })$ of the n variables satisfies the m constraints simultaneously, output ${\vec{\sigma }}^{* }$ as a solution of the instance, otherwise output ‘UNCONVERGED’. |
The subroutine MRBP-UPDATE is to update the messages on edges of a factor graph until a fixed point or a maximal number of iteration steps.
Subroutine MRBP-UPDATE:
| Step 1 At t = 0: for each edge (a, i), randomly initialize ${\mu }_{a\to i}^{0}({\sigma }_{i})\in [0,1]$. |
| Step 2 For the constraints from a = 1 to m: |
| 2.1 For each edge (a,i) with i ∈ V(a), use equation ( |
| 2.2 Use equation ( |
| Step 3 From t = 1 to ${t}_{\max }$: |
| 3.1 For each edge (a,i), calculate the residual ${\gamma }_{a\to i}^{t}$ using equation ( |
| 3.2 For i ∈ V(aγ), use equation ( |
| 3.3 For each edge (a,i), if $| {\mu }_{a\to i}^{t+1}({\sigma }_{i})-{\mu }_{a\to i}^{t}({\sigma }_{i})| \lt \varepsilon $ holds, break the loop and set ${\mu }_{a\to i}^{* }({\sigma }_{i})={\mu }_{a\to i}^{t}({\sigma }_{i})$. Otherwise go to 3.1 until all edges are marked. |
| 3.4 If $t={t}_{\max }$, output ‘UNCONVERGED’. |
The subroutine MRBP-UPDATE* is to obtain the fixed point of message updating in the decimation process, during which some variables have already been assigned with certain values.
Subroutine MRBP-UPDATE*:
| Step 1 At t = 0: |
| • For variables of Type A: skip; |
| • For variables of Type B: if the assignment of the variable i and the value of the fixed variable satisfy the corresponding constraint a, set ${\mu }_{a\to i}^{0}({\sigma }_{i})=1;$ otherwise, set ${\mu }_{a\to i}^{0}({\sigma }_{i})=0;$ |
| • For variables of Type C: uniformly initialize the values of ${\mu }_{a\to i}^{0}({\sigma }_{i})\in [0,1]$ at random. |
| Step 2 For the constraints from a = 1 to m: |
| 2.1 For each edge (a,i), where i ∈ V(a). |
| • For variables of Type A and Type B: skip; |
| • For variables of Type C: use equation ( |
| 2.2 Update ${\mu }_{a\to i}^{1}({\sigma }_{i})$. |
| • For variables of Type A: skip; |
| • For variables of Type B: let ${\mu }_{a\to i}^{1}({\sigma }_{i})={\mu }_{a\to i}^{0}({\sigma }_{i});$ |
| •For variables of Type C: use equation ( |
| Step 3 From t = 1 to ${t}_{\max }$: |
| 3.1 Compute the residual ${\gamma }_{a\to i}^{t}$. |
| • For variables of Type A and Type B: skip; |
| • For variables of Type C: for each edge (a,i), calculate the residual ${\gamma }_{a\to i}^{t}$ using equation ( |
| 3.2 For i ∈ V(aγ), use equation ( |
| 3.3 Determine whether the iterative equations converge or not. |
| • For variables of Type A and Type B: skip; |
| • For variables of Type C: for each edge (a,i), if $| {\mu }_{a\to i}^{t+1}({\sigma }_{i})-{\mu }_{a\to i}^{t}({\sigma }_{i})| \lt \varepsilon $ holds, break the loop and set ${\mu }_{a\to i}^{* }({\sigma }_{i})={\mu }_{a\to i}^{t}({\sigma }_{i})$. Otherwise go to 3.1 until all edges are marked. |
| 3.4 If $t={t}_{\max }$, output ‘UNCONVERGED’. |
4. Results
Here we test our MRBP algorithm based on residuals and show that it not only reduces the fluctuation of the messages, but also significantly improves the convergence of the BP iterative equations. Hence, it can effectively construct a solution of a random instance generated by model RB when approaching the satisfiability threshold.
As a representative parameter set to generate random instances of binary model RB, we take α = 0.8 and r = 3. We should mention that, in random instances with other values of parameters of α and r, we can obtain similar results on the algorithmic performance around thresholds. In Table 1, for different problem sizes n, corresponding quantities and thresholds are shown. In the MRBP algorithm, we take ${t}_{\max }=400$ and ϵ = 10−4.
Table 1. For different problem sizes n, the domain size dn, the constraint number m, and the theoretical satisfiability threshold ps obtained by Theorem |
| n | α | r | dn | m | ps |
|---|---|---|---|---|---|
| 20 | 0.8 | 3 | 11 | 180 | ≃ 0.23 |
| 40 | 0.8 | 3 | 19 | 443 | ≃ 0.23 |
| 60 | 0.8 | 3 | 26 | 737 | ≃ 0.23 |
| 80 | 0.8 | 3 | 33 | 1052 | ≃ 0.23 |
As the first part of our results, we consider the performance of the MRBP algorithm. The fraction of successful runs over 50 random instances, which are generated by binary model RB for n ∈ {20, 40, 60, 80}, is shown in figure 4. It is observed that almost all instances are satisfiable when the constraint tightness p < 0.19. The MRBP algorithm can construct a solution efficiently when p < 0.20. However, the algorithm fails with probability tending to 1 when p > 0.22. Therefore, the algorithm shows a good performance when p approaches the theoretical satisfiability threshold ps ≃ 0.23. Unfortunately, the algorithm fails in the extreme hard region, which may be closely related to the structural transition of the solution space. It is worth noting that the performance of the MRBP algorithm is almost independent of the problem size n.
Figure 4. Fraction of satisfiable instances as a function of p obtained by the MRBP algorithm in solving 50 random instances generated by model RB. |
In the decimation process of the MRBP algorithm, we measure the entropy of the selected variable at each time step T for different problem sizes n. We define the entropy of the decimated variable at T as
$\begin{eqnarray}S({x}_{i}^{T})=-\sum _{{\sigma }_{i}\in D}{\mu }_{i}({\sigma }_{i})\mathrm{ln}{\mu }_{i}({\sigma }_{i}).\end{eqnarray}$
In figure 5, the results of the entropy in function of T/n for p = 0.1 and n = {20, 40} are presented. For a given n, the entropy curve shows a concave shape with a minimum at T/n ≈ 0.5, where half of the variables are fixed with certain values. For n = 20 and 40, the entropy $S({x}_{i}^{T})$ is always less than the maximum in the case in which each variable is evaluated with equal probability from its domain, which guarantees that the MRBP algorithm can select a polarized variable among the remaining free variables to fix.
Figure 5. Entropy $S({x}_{i}^{T})$ of the decimated variables at each step T in the MRBP algorithm for binary model RB as a function of T/n at p = 0.1. |
We then consider the mean degree of freedom of the n variables in the decimation procedure of the BPD algorithm, which is defined as
$\begin{eqnarray}s(p)=\displaystyle \frac{1}{n}\sum _{i=1}^{n}\sum _{{\sigma }_{i}\in D}[-{\mu }_{i}({\sigma }_{i})\mathrm{ln}{\mu }_{i}({\sigma }_{i})].\end{eqnarray}$
In figure 6, the relationship between the mean degree of freedom s(p) and the constraint tightness p for n = 40 is reported. It is shown that, with an increase in constraint tightness p, s(p) monotonously decreases to 0. Therefore, we suspect that as p increases, some variables are frozen at certain assignments, which prevents the MRBP algorithm from escaping from the locally optimal solution. Further improving the performance of the algorithm involves studying in more detail the structural evolution of the solution space of model RB.
Figure 6. Mean degree of freedom s(p) on n variables in function of p during the execution of MRBP algorithm on a random instance generated by binary model RB for n = 40. |
As the second part of our results, we present a comparison between the usual BP algorithm (algorithm 3 in [39]) and the MRBP algorithm. Figure 7 shows their performance comparison in terms of message convergence and solving power. It can be seen from figure 7 (a), for most instances, that if the BP algorithm converges, it can usually construct a solution of the instance. However, this is not the case for the MRBP algorithm. As shown in figure 7 (b), when p < 0.21, if the MRBP algorithm converges, it almost always finds the solution of the instance. However, when p ≥ 0.21, the assignment obtained after the algorithm has converged may not be the solution of the instance. In other words, compared with algorithm 3 in [39], although the convergence of the MRBP algorithm has been greatly improved, the constructed assignment may not be a solution of the instance. As we can see from figure 7 (c), the convergence of the MRBP algorithm is significantly better than algorithm 3 in [39]. For the MRBP algorithm, the number of convergent instances decreases first and then increases with an increase in p. Moreover, the convergence is the worst at p = 0.21, where the probability of the MRBP algorithm finding a solution of a random instance is extremely low. In figure 7 (d), the MRBP algorithm shows a consistently higher success probability in finding ground-state solutions. It is obvious that the MRBP algorithm not only improves the convergence of BP iterative equations but also greatly improves the efficiency of finding solutions.
Figure 7. (a) Number of convergent instances and number of instances for which a solution can be found by algorithm 3 in [39] when solving 50 random instances of binary model RB for n = 40. (b) Results of the same content with (a) yet obtained from the MRBP algorithm. (c) and (d) Comparison of the two algorithms on the number of convergent instances and the number of successful instances, respectively. |
We also compare the solving efficiency of the MRBP algorithm with three related algorithms based on BP (i.e. ABP algorithm with γ = 0.5, SBP algorithm with γ = 0.5, and VABP algorithm) in [41] and algorithm 3 in [39]. The results of 50 random instances generated by binary model RB for n = 40 are shown in figure 8, which illustrate that the MPBP algorithm can find a solution of the problem with a high probability when approaching the satisfiability phase transition region.
Furthermore, from the perspective of computational cost, we compare the average number of iterations required at each time step for the decimated variable if the BP equations converge. In figure 9, the diagram of average iterations is shown for n = 60 at p=0.15, where the five algorithms can construct a solution of a random instance with probability 1. We can see that the number of iterations required by the MRBP algorithm in the early decimation process is much lower than that of the ABP algorithm with γ = 0.5, SBP algorithm with γ = 0.5, VABP algorithm in [41], and algorithm 3 in [39]. After approximately 90% of the variables are fixed, the iterations required by the five algorithms are drastically reduced. Therefore, compared with other novel related algorithms, the MRBP algorithm greatly improves the convergence of message passing algorithms.
Finally, we analyze the time complexity of the MRBP algorithm. The typical time complexity of BP algorithms is O(M) as M is the size of the edges in a graph instance. In our modification on BP algorithms, there are an extra O(1) message updating steps after each iteration of messages on all edges. Thus the complexity of the running time scales as $O({n}^{2}\mathrm{ln}n)$, in which ${n}^{2}\mathrm{ln}n$ is simply the size of the edges in model RB here. In figure 10, we show the running time of the MRBP algorithm as a function of p for n = 20. We can see that the running time increases sharply when p is close to the satisfiability threshold, a common phenomenon for search algorithms approaching thresholds in CSPs.
Figure 10. Running time of the MRBP algorithm as a function of p. We consider here n = 20. |
Summing the above results together, for the problem of model RB, the MRBP algorithm greatly improves the convergence of BP iterative equations and shows significantly higher probabilities in finding solutions approaching the satisfiability threshold with a large reduction of iterations in computational cost.
5. Conclusion
In this paper, we propose an improved message passing algorithm based on residuals of messages to solve CSPs. The residual of messages is to quantify the fluctuation of messages as a time series between two consecutive time steps. As a modification to the usual message updating process, those messages with maximal residuals are given higher priority in the updating process. We test our MRBP algorithm on model RB, a random CSP with growing domains, which has an exact satisfiability threshold by previous analytical study. Numerical results show that our algorithm outperforms the BP algorithm with a better convergence of message updating, a higher probability in constructing a solution for the problem, and a lower computational cost.
To explore the potential implication of residuals on optimization problems and message passing algorithms, several lines of research can be carried out as future work. The first one is to combine our algorithm with a detailed description of the landscape of hard optimization problems to develop better problem-solving algorithms. The second one is to extend our algorithm to combinatorial optimization and CSPs with fixed domains, such as the minimum vertex cover problem and the k-SAT problem, which have a more complex and yet a more refined picture of solution space structure of ground-state configurations. At the same time, the idea of residuals of messages can be introduced into message passing algorithms in contexts other than sparse graphs, such as the approximated message passing algorithms on dense graphs, cluster variational message passing algorithms on lattice structures, and so on.