Journal of Statistical Physics

, Volume 159, Issue 5, pp 1154–1174

Statistical Mechanics of the Minimum Dominating Set Problem

  • Jin-Hua Zhao
    • 1
  • Yusupjan Habibulla
    • 1
  • Hai-Jun Zhou
    • 1
  1. 1.State Key Laboratory of Theoretical Physics, Institute of Theoretical PhysicsChinese Academy of SciencesBeijingChina
Article

DOI: 10.1007/s10955-015-1220-2

Cite this article as:
Zhao, JH., Habibulla, Y. & Zhou, HJ. J Stat Phys (2015) 159: 1154. doi:10.1007/s10955-015-1220-2

Abstract

The minimum dominating set (MDS) problem has wide applications in network science and related fields. It aims at constructing a node set of smallest size such that any node of the network is either in this set or is adjacent to at least one node of this set. Although this optimization problem is generally very difficult, we show it can be exactly solved by a generalized leaf-removal (GLR) process if the network contains no core. We present a percolation theory to describe the GLR process on random networks, and solve a spin glass model by mean field method to estimate the MDS size. We also implement a message-passing algorithm and a local heuristic algorithm that combines GLR with greedy node-removal to obtain near-optimal solutions for single random networks. Our algorithms also perform well on real-world network instances.

Keywords

Dominating setSpin glassCore percolationLeaf removal Network coarse-grainingBelief propagation

1 Introduction

The minimum dominating set (MDS) problem [1] has fundamental importance in network science. For example, to ensure the proper functioning of a complex networked system such as a nation-wide power grid, it is often necessary to monitor the system’s microscopic dynamics in real-time by placing sensors on the nodes. A sensor may have the capability of observing the instantaneous states of the residing node and all its adjacent nodes in the network [2], so they may not need to occupy all the nodes. We then have the MDS problem: How to place sensors on as few nodes as possible to minimize costs but still ensure that each node is either occupied or adjacent to at least one occupied node? As an example we show in Fig. 1b a minimum dominating set (containing only three nodes) of a small network. A more stringent constraint, which is adopted in lattice glass models [3], is to require an empty node i to be surrounded by at least li occupied nodes, with li being node-dependent. The MDS problem corresponds to the case of li1, while the other limiting case of li=di is just the vertex cover (or independent set) problem [4, 5], where di is node i’s degree (i.e., number of adjacent nodes).
https://static-content.springer.com/image/art%3A10.1007%2Fs10955-015-1220-2/MediaObjects/10955_2015_1220_Fig1_HTML.gif
Fig. 1

An example of minimum dominating set. a a small network with N=11 nodes and M=18 links. bBlue (dark gray) indicates a node being occupied, while cyan (light gray) indicates a node being empty but observed. The three occupied nodes form a MDS Γ0={4,7,10} for this network. c A coarse-grained representation of the network based on the MDS of (b) (Color figure online)

The MDS problem has wide practical applications, such as monitoring large-scale power grids and other transportation systems [2], controlling the spreading of infectious diseases and other network dynamical processes [69], efficient routing in wireless networks [10], and network public goods games (e.g., resource allocation) [11]. Another application is to build a coarse-grained representation for a complex network starting from a MDS. Such an idea has already been applied to multi-document summarization in the field of information extraction [12]. Each node i of the MDS can be regarded as a representative node for a local domain of the network. We can take the subnetwork induced by node i and all its adjacent nodes (except those in the MDS) as a coarse-grained node, and set up an coarse-grained link between two coarse-grained nodes if the two corresponding subnetworks share at least one node or are connected by at least one link in the original network (see Fig. 1c for an example). If such a coarse-graining process is iterated we will then obtain a hierarchical representation for the original network, which may be very useful for understanding the organization of a complex system and for searching and information transmission in such a system.

Exactly solving the MDS problem, however, is extremely difficult in general, since it is a nondeterministic polynomial-complete (NP-complete) optimization problem [1]. Even the task of approximately solving the MDS problem is very hard. For a general network of N nodes, so far the best polynomial algorithms can only guarantee to get dominating sets with sizes not exceeding lnN times of the minimum size [13, 14]. Many local-search algorithms have been proposed to solve the MDS problem heuristically (see review [1] and [2, 6, 7, 9, 15, 16]), but theoretical results on the MDS sizes of random network ensembles are still very rare.

In this work we bring several new theoretical and algorithmic contributions. We show in Sect. 2 that a generalized leaf-removal (GLR) process may cause a core percolation transition, and propose a quantitative theory to describe this percolation. If the network contains no core, GLR reaches an exact MDS; if an extensive core exists, we combine GLR with a local greedy process in Sect. 3 to get an upper bound to the MDS size. We then introduce a spin glass model in Sect. 4 and estimate the MDS size by a replica-symmetric (RS) mean field theory, and implement a message-passing algorithm in Sect. 5 to get near-optimal dominating sets for single random network instances. Our algorithms also perform well on real-world network instances. This work shall be useful both for network scientists who are interested in applying the MDS concept to practical problems, and for applied mathematicians who seek better theoretical understanding on the random MDS problem.

2 Generalized Leaf-Removal and Core Percolation

Consider a simple network W formed by N nodes and M undirected links, each link connecting between two different nodes. Each node with index i{1,2,,N} is either empty (indicated by the occupation state ci=0) or occupied by sensors (ci=1). A node i is regarded as observed if it is occupied or it is empty but adjacent to one or more occupied nodes, otherwise it is regarded as unobserved. We need to occupy a set Γ of nodes to make all the N nodes be observed, and the objective is to make the dominating set Γ as small as possible, i.e., to construct a minimum dominating set. It is easy to verify that the three occupied nodes of Fig. 1b form a MDS for that small network. Notice a network may have more than one MDS.

2.1 The GLR Process

Here we extend the leaf-removal idea of [17] (see also more recent work [6, 18, 19]) and consider a generalized leaf-removal process. This dynamics is based on the following two considerations: first, as pointed out in [6, 17], if node i is an unobserved leaf node (which has only a single neighbor, say j), then occupying j but leaving i empty must be an optimal strategy; second, we notice that if i is an empty but observed node and at most one of its adjacent nodes is unobserved, then it must be an optimal strategy not to occupy i. This second point was not considered in the conventional leaf-removal process [6].

The GLR process simplifies the input network W at discrete evolution steps t=0,1,2,. For the convenience of description, let us denote by Wt the simplified network at the start of the t-th evolution step of GLR. W0 at the initial step t=0 is identical to the original network W, and all the nodes of W0 are unobserved. We prove that if GLR makes the whole input network W be observed, then the set of nodes occupied during this process must be a MDS. For this latter purpose, let us denote by Γ0 a MDS of the input network W (there must be at least one such set). The essential idea is to demonstrate that during GLR, we can modify Γ0 in such a way that its size does not change but all the nodes i that are fixed to be occupied (ci=1) are in Γ0 while all the nodes j that are fixed to be unoccupied (cj=0) are not in Γ0. Starting from evolution step t=0, let us perform GLR and modify Γ0 in the following sequential order:
  1. (0)

    As long as there is an isolated node i in network Wt, fix its occupation state to ci=1 and delete it from Wt. All such fixed nodes i must also belong to Γ0.

     
  2. (1)

    As long as there is a leaf node i in network Wt which is not yet observed, fix the occupation state of its unique neighbor j to cj=1 and fix that of i to ci=0 so that j and all its adjacent nodes (including i) are now observed, see Fig. 2 (left panel). We then delete node j and all its connected links from Wt. If j belongs to Γ0 then node i must not belong to it, because otherwise Γ0 could not have been a MDS. On the other hand, if node j does not belong to Γ0 then node i must belong to it, and in this latter case we modify Γ0 by adding j to it and deleting i from it.

     
  3. (2)

    Then as long as there is a node i which is itself observed and which has only a single unobserved neighbor j, delete the link (i,j) from network Wt, see Fig. 2 (right panel). We do not modify Γ0 if node i does not belong to it. If node i does belong to Γ0 then node j must not belong to it, and in this latter case we add j to Γ0 and delete i from it.

     
  4. (3)

    Then as long as there is an observed node i which is not connected to any unobserved node, fix its occupation state to ci=0 and delete it and all its attached links from Wt. Such a node i must not belong to Γ0, for otherwise Γ0 could not have been a MDS.

     
  5. (4)

    If the resulting network Wt is empty or it contains no isolated node nor leaf node, the GLR process stops. If Wt still contains at least one isolated or leaf node, then we increase the evolution step from t to (t+1) and initialize the network Wt+1 as identical to Wt. A node i of Wt+1 is regarded as observed if and only if it is observed in network Wt. We then repeat the above-mentioned operations (1)–(3).

     
If the final simplified network is non-empty, then there must be some nodes that are still unobserved after the GLR process. The subnetwork induced by these unobserved nodes is referred to as the core of the original network W. This core is connected only to observed empty nodes but not to occupied nodes. We denote by ncore the fraction of nodes in this core and by w the fraction of occupied nodes.
https://static-content.springer.com/image/art%3A10.1007%2Fs10955-015-1220-2/MediaObjects/10955_2015_1220_Fig2_HTML.gif
Fig. 2

The two basic operations of the generalized leaf-removal process. White circles denote empty and unobserved nodes, cyan (light gray) circles denote empty but observed nodes, and blue (dark gray) circles denote occupied nodes. Left panel the unique adjacent node j of an unobserved leaf node i is occupied, and all the neighbors of j are observed. Right panel an empty observed node i has only a single unobserved neighbor j, then the link between i and j is deleted (Color figure online)

If the original network W has no core, then the set Γ of occupied nodes by the GLR process must be identical to the final Γ0, which is a MDS modified from the original MDS. We have therefore proven that GLR constructs a MDS for a network W if this network contains no core. (All the above-mentioned modification operations on Γ0 are ignored in the actual implementation of the GLR process. They are introduced here solely for proving that GLR is able to construct a MDS if there is no core.) Furthermore, we notice that if the GLR process finishes with some nodes remaining to be unobserved, the set of nodes occupied during this process must be a MDS for the subnetwork of W induced by all the observed nodes. This is because all these occupied nodes also belong to the modified MDS Γ0, while all those nodes fixed to be unoccupied are outside of Γ0.

We generate many large instances of Erdös–Rényi (ER) and scale-free (SF) random networks and run the GLR process on them (details of the network sampling method are given in Sects. 2.3, 2.4, and 2.5). Some representative results are shown in Fig. 3 for ER networks [20, 21], in Fig. 4 for SF networks generated through the static model [22], and in Fig. 5 for pure SF networks [20, 21]. A major observation is that there is no core in pure SF random networks with minimum node degree dmin=1, therefore a MDS for such a network can be easily constructed by the GLR process. Another major observation is that there is a continuous core percolation transition in ER networks and in SF networks generated through the static model. This core percolation transition occurs at certain threshold value of the mean node degree. For example, for ER networks with N=106 nodes and M=(c/2)N links, when the mean node degree c<2.41 there is no core (ncore=0), and GLR reaches a MDS for the whole network (Fig. 3). The core emerges at c2.41 and its relative size ncore then increases with c continuously from zero. For c>2.41, GLR constructs a MDS only for part of the ER network and it leaves an extensive core of ncoreN unobserved nodes.
https://static-content.springer.com/image/art%3A10.1007%2Fs10955-015-1220-2/MediaObjects/10955_2015_1220_Fig3_HTML.gif
Fig. 3

Generalized leaf-removal on Erdös–Rényi random networks. w and ncore are the fractions of occupied and unobserved nodes, respectively. Cross symbols are results obtained by running GLR on a single ER network of size N=106 and mean degree c; solid lines are the predictions of the percolation theory for N=

https://static-content.springer.com/image/art%3A10.1007%2Fs10955-015-1220-2/MediaObjects/10955_2015_1220_Fig4_HTML.gif
Fig. 4

Generalized leaf-removal on scale-free random networks of decay exponent λ=3.5,3.0,2.667,2.5 (from left to right) generated through the static model [22]. w and ncore are the fractions of occupied and unobserved nodes, respectively. Red dash-dotted lines are results obtained by running GLR on a single network instance of size N=106 and mean degree c, while blue dashed lines are results obtained by the core percolation theory using the degree profile of this network instance as input. Black solid lines are the predictions of the percolation theory for N= (Color figure online)

https://static-content.springer.com/image/art%3A10.1007%2Fs10955-015-1220-2/MediaObjects/10955_2015_1220_Fig5_HTML.gif
Fig. 5

Generalized leaf-removal on pure scale-free random networks with minimum degree d=1. w is the fraction of occupied nodes. The fraction ncore of unobserved nodes is simply ncore=0. Red triangle symbols are results obtained by running GLR on a single pure SF network of size N=106 and decay exponent λ, while blue cross symbols are results obtained by the percolation theory using the degree profile of this network instance as input. The black solid line is obtained by the percolation theory at N= (Color figure online)

Notice the core percolation transition resulting from the GLR optimization process is qualitatively different from the simpler observability transition discussed in [2], which considers the appearance of a giant connected component of observed nodes resulting from an initial set of randomly chosen occupied nodes. We now develop a percolation theory to thoroughly understand the GLR dynamics on random networks.

2.2 The Core Percolation Theory

A random network is characterized by a degree distribution P(d), which gives the fraction of nodes with degree d0 [20]. We assume that there is no correlation between the degrees of adjacent nodes, therefore the degree d of a node reached by following a randomly chosen link obeys the distribution Q(d) of the form
Q(d)P(d)dc,
(1)
where cd0P(d)d is the mean node degree of the network. Consider a link (i,j) of the network W. Let us neglect for the moment the constraint of node i to node j but only consider the other adjacent nodes of j. If the constraint of node i is neglected, then what is the probability αt that node j becomes an unobserved leaf node (i.e., it has no other adjacent node except i) at the start of the t-th GLR evolution step? What is the probability βt that j becomes newly occupied (cj=1) at the tth GLR step? What is the probability γt that j is observed but not occupied at the end of the t-th GLR step? And what is the probability ηt that at the end of the t-th GLR step, node j is an observed and unoccupied node and it has no unobserved adjacent node except i? For an uncorrelated random network these four sets of probability parameters {α0,α1,}, {β0,β1,}, {γ0,γ1,}, and {η0,η1,} can be computed by a set of iterative equations.
The expressions of α0,β0,γ0,η0 for the initial evolution step t=0 are
α0=Q(1),
(2)
β0=1d1Q(d)(1α0)d1,
(3)
γ0=d1Q(d)[(1α0)d1(1α0β0)d1],
(4)
η0=d1Q(d)[(β0+γ0)d1γ0d1].
(5)
Equation (2) is trivial, it simply describes the situation that node j has only a single neighbor (i.e., node i). Equation (3) describes the situation that node j is adjacent to at least one leaf node (except i), which will guarantee j to be occupied at the t=0 GLR step. A random network has only very few short loops and therefore the local network structure around node j is a tree. In the core percolation theory we therefore assume that the adjacent nodes of j are completely independent of each other when j is still unobserved (such an assumption was also exploited in our earlier percolation studies [2325]). Based on this assumption, the probability of all the adjacent nodes (except i) of j not being unobserved leaves is then written in Eq. (3) as the product of the individual probability (1α0) of an adjacent node not being an unobserved leaf. Equation (4) expresses the fact that for node j to be an unoccupied but observed node at the end of the t=0 evolution step, it should not be adjacent to any unobserved leaf node but at least one of its adjacent nodes (except i) should be occupied.

If node j is adjacent to one or more nodes that are occupied at the t=0 evolution step and all its other adjacent nodes (except i) are observed at this evolution step, then at the end of this evolution step j is unoccupied but observed and it is not adjacent to any unobserved node (except i). This then leads to the expression (5) for η0. Notice such an observed but unoccupied node j will be deleted at the end of the t=0 evolution step. After all such nodes are deleted, some unobserved nodes in the remaining network may become isolated or be connected to only a single node. If this is the case, these isolated or leaf nodes will trigger the next (t=1) evolution step.

Following the same line of theoretical considerations, we obtain the expressions of αt, βt, γt, and ηt for the tth GLR evolution step (t1) as
αt={d1Q(d)(η0)d1Q(1),(t=1)d1Q(d)[(l=0t1ηl)d1(l=0t2ηl)d1],(t2)
(6)
βt=d1Q(d)[(1l=0t1αl)d1(1l=0tαl)d1],
(7)
γt=d1Q(d)[(1l=0tαl)d1(1l=0t(αl+βl))d1],
(8)
ηt=d1Q(d)[(l=0tβl+γt)d1(γt)d1]l=0t1ηl.
(9)
Let us denote by γlim the value of γt at the last evolution step t=tlim of the GLR process (notice that the maximal evolution step tlim may approach infinity for a network with N= nodes). Furthermore, we define the accumulated values of αt, βt, and ηt as
αcumt0αt,βcumt0βt,ηcumt0ηt.
There are the following relationships among αcum, βcum, ηcum and γlim:
αcum=d1Q(d)(ηcum)d1,
(10)
βcum=1d1Q(d)(1αcum)d1,
(11)
γlim=d1Q(d)[(1αcum)d1(1αcumβcum)d1],
(12)
ηcum=d1Q(d)[(βcum+γlim)d1(γlim)d1].
(13)
After all the probability parameters αt, βt, γt, ηt (for t=0,1,) for a node j at the end of a link (i,j) are determined by neglecting the constraint associated with node i, we now ask the following two questions: If the constraint of node i to all its adjacent nodes are considered, then what is the probability ncore of i to be unobserved after the whole GLR process? And what is the probability It of node i to be occupied at the t-th GLR evolution step? If node i remains to be unobserved during the whole GLR process, it must not be adjacent to any unobserved leaf node nor to any occupied node, and it must have at least two adjacent nodes after the whole GLR process. Therefore we obtain that
ncore=d2P(d)s=0d2Cds(ηcum)s(1αcumβcumηcum)ds=d1P(d)[(1αcumβcum)d(ηcum)dd(ηcum)d1(1αcumβcumηcum)],
(14)
where Cdsd!/[s!(ds)!] is the binomial coefficient. Notice that if (αcum+βcum+ηcum)=1 then we have ncore=0.
It is easy to see that the probability I0 of a randomly chosen node i to be occupied at the t=0 GLR evolution step is
I0=1P(1)(1α02)d2P(d)(1α0)d.
(15)
The coefficient 1/2 in the second term of the above expression reflects the fact that if node i has only one neighbor j, then i has one-half probability to be occupied if j also has only one neighbor (namely i).
If a randomly chosen node i is not occupied at the t=0 evolution step, then the probability I1 of it being occupied at the t=1 evolution step is
I1=d2P(d)(η0)d+d2P(d)[(1α0)d(1α0α1)ddα1((β0+γ0)d1(γ0)d1)]12d2P(d)dα1(η0)d1.
(16)
All the adjacent nodes of i might haven been deleted at the end of the t=0 evolution step. If this is the case node i becomes isolated at the start of the t=1 evolution step, which leads to the first summation of Eq. (16). The second summation of Eq. (16) corresponds to the other situation of node i not being occupied nor being deleted at the t=0 evolution step but it is adjacent to at least one node that becomes an unobserved leaf at the start of the t=1 evolution step. Notice if node i becomes an unobserved leaf node at the start of the t=1 evolution step with its unique neighbor also being such a leaf node, then i has only one-half probability to be occupied at this evolution step. This last situation leads to the third summation term of Eq. (16), which corrects the over-counted probability of occupation in the second summation term.
Following the same line of theoretical considerations, we obtain the probability It of a randomly chosen node i changing from being unoccupied to being occupied at the tth GLR evolution step (t2):
It=d2P(d)[(l=0t1ηl)d(l=0t2ηl)ddηt1(l=0t2ηl)d1]+d2P(d)[(1l=0t1αl)d(1l=0tαl)d]d2P(d)dαt[(l=0t1βl+γt1)d1(γt1)d1+(l=0t2ηl)d1]12d2P(d)dαt[(l=0t1ηl)d1(l=0t2ηl)d1].
(17)
The probability w of a randomly chosen node i to be occupied during the GLR process is then
w=t0It
(18)
=1P(1)(1α0/2)d2P(k)[(1αcum)d(ηcum)d]t1d2P(d)d[ηt(l=0t1ηl)d1+αt(l=0t1βl+γt1)d1αt(γt1)d1]12t2d2P(d)dαt[(l=0t1ηl)d1+(l=0t2ηl)d1]12d2P(d)dα1(η0)d1.
(19)
Our core percolation theory can be applied both to single finite random network instances and to random network ensembles at the thermodynamic limit N. For each t (starting from t=0), we first compute αt, then use αt as input to compute βt, then use αt and βt as inputs to compute γt, and finally use αt, βt and γt as inputs to compute ηt. For a finite random network of N nodes, the iteration stops if the evolution step t increases to a value tlim such that Itlim<1/N. This is because if NIt<1 the number of newly occupied nodes has a high probability to be zero and then GLR will be unable to continue. For the case of N, the numerical iteration process can be carried out to a sufficiently large evolution step t=tlim until αtlim0.

2.3 Results on Erdös–Rényi Random Networks

We generate an ER random network W of N nodes and M=(c/2)N links by adding links sequentially to an initial network of N isolated nodes. To add a link, we choose two different nodes i and j uniformly at random from the whole node set and set up a link (i,j) between them if this link has not been created before. The mean node degree of the resulting network W is equal to c. When the number N of nodes is sufficiently large the degree distribution P(d) of such a ER network obeys the Poisson distribution [20, 21]
P(d)=cdecd!(d0).
(20)
For this network ensemble, the predicted results of ncore and w by our core percolation theory are in perfect agreement with simulation results (see Fig. 3). Especially, at the thermodynamic limit N, the theory predicts a continuous core percolation phase transition at c2.4102, which is slightly lower than the core percolation phase transition point of c2.7183 caused by the conventional leaf-removal process [17]. Before the GLR-induced core percolation transition occurs, the occupation fraction w obtained by Eq. (19) is equal to the ensemble-averaged MDS size (relative to N), but it is only a lower bound to this size when an extensive core emerges in the random network (ncore>0).

2.4 Results on Scale-Free Random Networks Generated Through the Static Model

Now let us consider GLR-induced core percolation on more heterogeneous random networks. We generate a scale-free network W of N nodes and M=(c/2)N links according to the static model [22]. Each node i{1,2,,N} is first assigned a fitness value θi=iξ/(j=1Njξ), where 0ξ<1 is a control parameter. Then we add links between pairs of these N nodes in a sequential manner. To create a link, two nodes i and j are chosen independently from the set of N nodes, and the probability that i and j being chosen is equal to θiθj; if nodes i and j are different and the link (i,j) has not been created before, this link is added to network. The final network W has a power-law degree distribution P(d)dλ for d1, with degree decay exponent λ=1+1/ξ. In the thermodynamic limit N, an explicit expression for P(d) is obtained as [26]
P(d)=[c(1ξ)]dd!ξ1dxec(1ξ)xxk11/ξ(d0).
(21)
For N=, a continuous core percolation phase transition is observed in such a SF network, and this transition occurs at more and more larger value of the mean node degree c as the decay exponent λ decreases (see Fig. 4 for λ=3.5, 3.0, 8/32.667, and 2.5 and Fig. 6 for 2<λ6). When the decay exponent λ is less then 3.0, theoretical predictions obtained at N= are quantitatively different from theoretical and simulation results obtained on finite (e.g., N=106) network instances, with the deviations become more pronounced as λ is closer to 2.0. Such a finite-size effect is mainly caused by the natural cutoff of maximum node degree in finite networks (it was also observed in our earlier work [23] on another type of percolation transitions). We emphasize that for a give finite value of N, the results of the core percolation theory agree with the simulation results of the actual GLR process very well, especially if we average the theoretical and simulation results over many network instances to reduce fluctuations.
https://static-content.springer.com/image/art%3A10.1007%2Fs10955-015-1220-2/MediaObjects/10955_2015_1220_Fig6_HTML.gif
Fig. 6

Core percolation phase transition in infinite (N=) random scale-free networks generated through the static model [22]. Horizontal axis is the degree decay exponent λ, while vertical axis is the value of the mean node degree c at the phase transition point. Cross symbols are predictions of the core-percolation theory, while the solid line is just a guide for the eye

For random SF networks generated through the static model with N= nodes, the core percolation transition value of mean node degree c is very sensitive to the decay exponent λ in the region of 2<λ<2.5, and it diverges as λ approaches 2.0 from above (Fig. 6). At the other limit of λ, the mean node degree at the phase transition approaches the value of c2.4102, which is just the core percolation phase transition point of an infinite ER random network.

2.5 Results on Pure Scale-Free Random Networks

When N=, a pure scale-free random network has the following degree distribution
P(d)=1k=1kλdλ(d1),
(22)
with λ>2 to ensure a finite value for the mean node degree c. For such a random SF network our core percolation theory predicts that ncore=0, namely there is no core percolation transition and the GLR process will construct a MDS for the whole network. The fraction w of occupied nodes (i.e., the size of a MDS relative to the node number N) decreases with the decreasing of the degree decay exponent λ (see Fig. 5), and it approaches zero as λ approaches 2.0 from above.

We also generate a set of pure SF random networks of finite size N following the same procedure as mentioned in [27] (see also the supplementary information of [23]). The minimum node degree of such a SF network is dmin=1, while the maximum node degree is dmaxN1/(λ1) [27]. When we apply both the GLR process and the core percolation theory on these finite network instances, we find the simulation results on the fraction ncore of nodes in the core and the fraction w of occupied nodes are in perfect agreement with the corresponding theoretical results (see Fig. 5). All these finite SF networks contain no core (ncore=0), and the MDS relative size w is an increasing function of λ.

Figure 5 also demonstrates strong finite-size effect for pure SF random networks of λ<3.0. This finite-size effect is again mainly caused by the cutoff of the maximum node degree of finite networks, which makes the mean node degree of a finite network be smaller than that of an infinite network. For example, at λ=2.1 the mean node degree of an infinite network is c61.49, while that of a finite network of size N=106 is reduced to c5.134.

3 Hybrid Local Algorithm

There is a very simple greedy algorithm in the literature to solve the MDS problem approximately, which is based on the concept of node impact [1, 7, 16]). The impact of an unoccupied node i equals to the number of nodes that will be observed by occupying i. For example, if node i has 3 unobserved neighbors, its impact is 4 if i is itself unobserved and is 3 if i is already adjacent to one or more occupied nodes. Starting from an input network W with all the nodes unobserved, the greedy algorithm selects uniformly at random a node i from the subset of nodes with the highest impact and fix its occupation state to ci=1. All the adjacent nodes of i are then observed. If there are still unobserved nodes in the network, the impact value for each of the unoccupied nodes is updated and the greedy occupying process is repeated until all the nodes are observed. This pure greedy algorithm is very easy to implement and very fast, but we find that it usually fails to reach a true MDS even when the input network contains no core.

Here we introduce an improved local algorithm by combining the GLR process with the impact-based greedy process. We call this new algorithm the GLR-Impact hybrid algorithm. Given an input network W with all the nodes unobserved, we first carry out the GLR process to simplify W as far as possible. If all the nodes are observed during this initial GLR process, a MDS of network W is then constructed. For the nontrivial case of some nodes being left unobserved after this GLR process, we first occupy a randomly chosen node from the subset of highest-impact nodes and then perform the GLR process again to further simplify the network as far as possible. We keep repeating such a occupying-followed-by-GLR process until there is no unobserved node left in the network.

The GLR-Impact hybrid algorithm is also very easy to implement and very fast. Its performance is demonstrated in Fig. 7 for single ER networks and regular random (RR) networks. (All the nodes of a RR network have the same integer degree c but the network is otherwise completely random [23]). This hybrid local algorithm outperforms the pure greedy algorithm considerably for c10, but it is still inferior to the belief-propagation-guided decimation (BPD) algorithm of Sect. 5.
https://static-content.springer.com/image/art%3A10.1007%2Fs10955-015-1220-2/MediaObjects/10955_2015_1220_Fig7_HTML.gif
Fig. 7

Constructing dominating sets for Erdös–Rényi networks (left panel) and regular random networks (right panel). The relative sizes w of dominating sets obtained by a single running of the pure greedy, the hybrid, and the BPD algorithm with x=10 on 96 ER or RR network instances of N=105 and (mean) degree c are compared (fluctuations are of order 104 and are not shown). The ensemble-averaged MDS relative sizes obtained by the replica-symmetric mean field theory are also shown (RS)

We also test the performance of the hybrid algorithm on single SF random networks generated through the static model [22] (see Fig. 8). The GLR-Impact algorithm still outperforms the pure greedy algorithm on these heterogeneous networks, and its performance approaches that of the BPD algorithm as the network becomes more and more heterogeneous (i.e., as the decay exponent λ approaches 2.0 from above).
https://static-content.springer.com/image/art%3A10.1007%2Fs10955-015-1220-2/MediaObjects/10955_2015_1220_Fig8_HTML.gif
Fig. 8

Constructing dominating sets for scale-free random networks generated through the static model [22]. The relative sizes w of dominating sets obtained by a single running of the pure greedy, the hybrid, and the BPD algorithm with x=10 on 96 SF network instances of N=105 and (mean) degree c are compared (fluctuations are of order 104 and are not shown). The degree decay exponent is λ=3.5, 3.0, 2.667, and 2.5, respectively

Real-world networks are often very heterogenous, with a small fraction of highly connected nodes [21]. As a test of the algorithms introduced in this work, we apply the GLR process, the pure greedy algorithm, the hybrid algorithm, and the BPD algorithm to a set of twelve real-world networks. Among these network instances, five are infrastructure networks: European express road network (RoadEU [28]), road network of Texas (RoadTX [29]), power grid of western US states (Grid [30]), and two Internet networks at the autonomous systems level (IntNet1 and IntNet2 [31]); three are information networks: Google webpage network (WebPage [29]), European email network (Email [32]), and research citation network (Citation [31]); three are social contact networks: collaboration network of condensed-matter authors (Author [32]), peer-to-peer interaction network (P2P [33]),and on-line friendship network (Friend [34]); the remaining one is the biological network of protein-protein interactions (PPI [35]).

The numerical results are summarized in Table 1. The GLR process is able to simplify these networks considerably. After GLR, the remaining number of unobserved nodes is often much smaller than the total number N of nodes in the original network. The BPD algorithm performs slightly better than the GLR-Impact hybrid algorithm, and both BPD and the hybrid algorithm outperform the pure greedy algorithm in all the twelve network instances.
Table 1

Results on twelve real-world network instances

Network

N

M

dmax

Core

Greedy

Hybrid

BPD

RoadEU

1177

1417

10

306

428

389

387

PPI

2361

6646

64

17

550

539

539

Grid

4941

6594

19

603

1564

1485

1481

IntNet1

6474

12,572

1458

8

660

656

656

Author

23,133

93,439

279

9052

3686

3612

3604

Citation

34,546

420,877

846

11,178

3335

3168

3095

P2P

62,586

147,892

95

35

12,710

12,582

12,582

Friend

196,591

950,327

14,730

6097

42,536

41,633

41,672

Email

265,214

364,481

7636

470

18,183

18,181

18,181

WebPage

875,713

4,322,051

6332

162,439

81,288

79,928

80,769

RoadTX

1,379,917

1,921,660

12

560,582

477,729

437,503

425,774

IntNet2

1,696,415

11,095,298

35,455

211,244

187,592

183,516

183,248

N and M are, respectively, the total number of nodes and links in the network; dmax is the maximum node degree of the network; the column marked by ‘Core’ records the number of nodes that are left unobserved after the GLR process; the columns marked by ‘Greedy’, ‘Hybrid’, and ‘BPD’ record the sizes of the dominating sets constructed by a single running of the pure greedy, the hybrid, and the BPD algorithm, respectively

4 Spin Glass Model and Replica-Symmetric Mean Field Theory

If a given network instance W contains an extensive core, the GLR process can only give a lower bound to the MDS size. We now discuss the issue of estimating the MDS size by way of a mean field theory. We introduce a partition function Z as
Z=c_iW{exci[1(1ci)ji(1cj)]},
(23)
where c_(c1,c2,,cN) denotes one of the 2N possible occupation configurations, x>0 is a re-weighting parameter, and i denotes node i’s set of adjacent nodes. The constraint of each node i leads to a multiplication term [1(1ci)ji(1cj)], which equals to 0 if i and all its adjacent nodes are empty and equals to 1 if otherwise. The partition function therefore only takes into account all the dominating sets, and at x it is contributed exclusively by the MDS configurations.

4.1 Replica-Symmetric Mean Field Theory

We solve the spin glass model (23) by a RS mean field theory, which can be understood from the angle of Bethe-Peierls approximation [36] or derived alternatively through partition function expansion [37, 38]. The marginal probability qic of node i’s occupation state being c ({0,1}) is expressed as
qic=excjicjqji(cj,c)δ0cjiqji(0,0)ciexcijicjqji(cj,ci)jiqji(0,0),
(24)
where the Kronecker symbol δmn=1 if m=n and δmn=0 if otherwise. The quantity qji(cj,ci) is defined as the joint probability that node i is in occupation state ci and its adjacent node j is in occupation state cj when the constraint of node i is not considered. This probability can be evaluated through the following belief-propagation (BP) equation:
qji(cj,ci)=excjkjickqkj(ck,cj)δ0ci+cjkjiqkj(0,0)ci,cjexcjkjickqkj(ck,cj)kjiqkj(0,0),
(25)
where ji denotes the subset obtained by deleting node i from set j.
The total free energy F is related to the partition function by F(1/x)lnZ. According to the RS mean field theory, its expression is
F=iWfi(i,j)Wf(i,j),
(26)
where fi and f(i,j) are the free energy contributions of a node i and a link (i,j) between nodes i and j:
fi=1xln[ciexcijicjqji(cj,ci)jiqji(0,0)],
(27)
f(i,j)=1xln[ci,cjqij(ci,cj)qji(cj,ci)].
(28)
From Eqs. (26) and (24) we can compute the free energy density fF/N and the mean occupation fraction w=(1/N)iWqi+1. The entropy density of the system is then evaluated as s=(wf)x.

4.2 Belief-Propagation Iterations

According to Eq. (25) each probability distribution qji(cj,ci) has the property that qji(1,1)=qji(1,0). Therefore in the numerical computations qji(cj,ci) can be represented by three non-negative real numbers qji(0,0), qji(0,1), and qji(1,0), which satisfy in addition the normalization condition
qji(0,0)+qji(0,1)+2qji(1,0)=1.
(29)
We initialize qji(cj,ci) and qij(ci,cj) for each link (i,j) of the network between two nodes i and j, for example setting qji(0,0)=qji(0,1)=qji(1,0)=1/4. We then perform BP iteration a number T of times at a given value of the re-weighting parameter x, until a fixed-point solution of Eq. (25) is reached or T exceeds a pre-specified number (e.g., 1000). In each BP iteration step we treat all the nodes of the network in a random order. When node j is examined, the output messages qji(cj,ci) to all its adjacent nodes ij are updated according to Eq. (25). The difference Δji(t) between an updated message qji(t) at the t-th BP step and the old message qji(t1) at the (t1)-th BP step is defined as
Δji(t)|qji(0,0)(t)qji(0,0)(t1)|+|qji(0,1)(t)qji(0,1)(t1)|+2|qji(1,0)(t)qji(1,0)(t1)|.
(30)
If the maximal value among the set of 2M difference values {Δji(t)} is less than certain pre-specified threshold value (e.g., 103 or even smaller), then BP iteration is regarded as being converged. At a fixed point of Eq. (25) we then compute the free energy density f, the mean occupation fraction w, and the entropy density s through the RS mean field theory. As an example, we show in Fig. 9 the results obtained on a single ER random network of size N=106 and mean degree c=10.
https://static-content.springer.com/image/art%3A10.1007%2Fs10955-015-1220-2/MediaObjects/10955_2015_1220_Fig9_HTML.gif
Fig. 9

Replica-symmetric (RS) mean field theory and belief-propagation (BP) results for ER random networks of mean degree c=10. The RS results are obtained by population dynamics simulations, while the BP results are obtained on a single ER network instance of N=106 nodes. The BP iteration converges to a fixed point only for x<8.22. a Occupation fraction w. b Free energy density f. c Entropy density s. d Entropy density s as a function of occupation fraction w

For ER networks with mean degree c>2.41 and regular random networks with integer degree c3, we find that when the re-weighting parameter x is larger than certain threshold value, BP iteration is unable to converge to a fixed point. Such a non-convergence phenomenon indicates that, when the random network system has an extensive core, it will be in a spin glass phase at sufficiently large values of x. Systematic theoretical investigations on this spin glass phase will be reported in another publication.

4.3 Ensemble-Averaged Properties

A random network ensemble is characterized by a degree distribution P(d). We perform population dynamics simulations using Eqs. (25), (24) and (26) to obtain ensemble-averaged results. First, we create a long array A of N (e.g., 105) elements to store a set of messages, each of which represents a probability distribution qji(cj,ci) in the form of a three-dimensional vector satisfying Eq. (29): qji(qji(0,0),qji(0,1),qji(1,0)). We then repeatedly update elements of this array by the following procedure: (1) generate a random integer d1 according to the degree probability distribution Q(d); (2) draw (d1) elements qkj from array A uniformly at random, and then use these (d1) elements as input messages to Eq. (25) to compute a new message qji; (3) replace a randomly chosen element of array A with this new message. The message array A is expected to reach a steady state after it is updated a sufficient number of times (e.g., after each element of this array is updated 10,000 times on average).

We then keep updating the message array A and at the same time compute the thermodynamic quantities f, w, and s. For example, the free energy density f is obtained by
f=fi¯c2f(i,j)¯,
(31)
where fi¯ is the average of the free energy node contribution fi over all the nodes, and f(i,j)¯ is the average of the free energy link contribution f(i,j) over all the links. We generate many samples of fi and f(i,j) to compute their averages fi¯ and f(i,j)¯. The procedure of obtaining a sample of fi is the same as that of updating an element of the message A, the only difference being that the degree di of node i should be generated according to the distribution P(d) instead of Q(d). A sample of f(i,j) is obtained very easily through Eq. (28) by picking two messages qji and qij uniformly at random from the message array A.

For ER random networks with mean degree c=10, we compare in Fig. 9 the results obtained by this RS population dynamics with the results obtained by BP iteration on a single network instance. The ensemble-averaged results are in perfect agreement with the BP iteration results (provided the BP iteration is able to converge).

The entropy density s as a function of the mean occupation fraction w can be obtained from these RS population dynamics results (see for example Fig. 9d). In some random network systems, the entropy density s become negative if w decreases below certain threshold value w0, indicating that there is no dominating set with relative size below w0. We therefore take the value w0 as the ensemble-averaged MDS relative size. For ER networks of c=10, we obtain from Fig. 9 that w00.120 (the corresponding value of x is x8.637). In some other random network systems (e.g., ER random networks with c<2.41, before the core percolation transition), the entropy density s approaches a non-negative limiting value as w approaches a limiting value w0 from above. For these latter cases, we simply take w0 as the ensemble-averaged MDS relative size.

The ensemble-averaged results on the MDS sizes of ER and RR networks are shown in Fig. 7. For ER networks with mean node degree c<2.41 (before the core-percolation transition), the RS mean field results coincide with the results predicted by the core percolation theory. When the random network contains an extensive core, the results obtained by the pure greedy algorithm and the GLR-Impact algorithm are higher than the RS mean field predictions, but the results obtained by the BPD algorithm of the next section are very close to the RS mean field predictions.

5 Belief-Propagation-Guided Decimation Algorithm

For a given network W, the RS mean field theory gives an estimate for the occupation probability qi+1 of each node i, see Eq. (24). Such information is exploited in a BPD algorithm to construct a near-optimal dominating set. (Such an algorithm and its extensions have already been successfully applied to many other combinatorial optimization problems, e.g., the K-satisfiability problem [39, 40] and the vertex-cover problem [5]). At each round of the BPD process, unoccupied nodes with the highest estimated occupation probabilities are added to the dominating set, and the occupation probabilities for the remaining unoccupied nodes are then updated.

If a node j is unobserved (it is empty and has no adjacent occupied node), the output message qji(cj,ci) on the link (j,i) between j and node i is updated according to Eq. (25). On the other hand, if node j is empty but observed (it has at least one adjacent occupied node), this node then presents no restriction to the occupation states of all its unoccupied neighbors. For such a node j, the output message qji(cj,ci) on the link (j,i) is then updated according to the following equation:
qji(cj,ci)=excjkjickqkj(ck,cj)cj,ciexcjkjickqkj(ck,cj).
(32)
Similar to Eq. (32), the marginal probability distribution qici for an observed empty node i is evaluated according to
qici=excijicjqji(cj,ci)ciexcijicjqji(cj,ci).
(33)
It is easy to verify from Eq. (32) that qji(0,0)=qji(0,1) and qji(1,0)=qji(1,1). Notice that if all the nodes in the set ji are observed, then we derive from Eq. (32) that qji(0,0)=qji(1,0)=qji(0,1)=qji(1,1)=1/4. Because of this property, we need only to consider the links between unobserved nodes and the links between unobserved and observed nodes. All the other links (which are between observed nodes) do not need to be considered in the BP iteration Eqs. (25) and (32).
We implement the BPD algorithm as follows:
  1. (0)

    Input the network W, set all the nodes to be empty and unobserved and set all the probability distributions qji(cj,ci) to be the uniform distribution. Set the re-weighting parameter x to a sufficiently large value (e.g., x=10). Then perform the BP iteration a number T0 of rounds (e.g., T0=200). After these T0 iterations we compute the occupation probability qi+1 of each node i using Eq. (24).

     
  2. (1)

    Then occupy a small fraction r (e.g., r=0.01) of the unoccupied nodes that having the highest estimated occupation probabilities.

     
  3. (2)

    Then simplify network W by first deleting all the links between observed nodes, and then deleting all the isolated observed nodes.

     
  4. (3)

    If the resulting network W still contains unobserved nodes, we perform BP iteration for a number of T1 rounds (e.g., T1=10). The output message of an node i is updated either according to Eq. (24) or according to Eq. (33), depending on whether i is unobserved or observed. We then repeat operations (1)–(3) until all the nodes are observed.

     
In addition, we may first carry out the GLR process to simplify the network W as far as possible before running the BPD process. For real-world networks with some nodes being highly connected, we find that such a GLR simplifying step reduces the BPD running time considerably and also slightly reduces the size of the constructed dominating set.

The results of the BPD algorithm for random networks and for real-world networks are compared with the results obtained by the local heuristic algorithms in Figs. 7,  8, and Table 1. For ER and RR random networks, the BPD algorithm considerably beats both the pure greedy algorithm and the GLR-Impact hybrid algorithm; for very heterogeneous (e.g., scale-free) networks, the BPD algorithm only slightly outperforms the GLR-Impact algorithm.

6 Discussions

In this work, we proposed two heuristic algorithms (a GLR-Impact local algorithm and a BPD message-passing algorithm) and presented a core percolation theory and a replica-symmetric mean field theory for solving the network dominating set problem algorithmically and theoretically. We found that the GLR process may lead to a core percolation transition in the network (see Figs. 3 and  4). Our numerical results shown in Figs. 7,  8 and Table 1 suggested that the GLR-Impact algorithm and the BPD algorithm can construct near-optimal dominating sets for random networks and real-world networks.

There are many theoretical issues remaining to be investigated. An easy extension of the core percolation theory is to consider GLR with a subset of initially occupied nodes. By optimizing this initial subset (e.g., following the methods of [4143]), we may reach an improved lower-bound to the MDS size. Core percolation on degree-correlated random networks [44] and in the more general lattice glass problem [3] are also very interesting. When the random network has an extensive core, we observed that the belief-propagation Eq. (25) fails to converge at large values of the re-weighting parameter x (see Fig. 9), indicating a spin glass phase transition. A systematic study of the spin glass phase will be carried out using the first-step replica-symmetry-breaking mean field theory [40, 45, 46], which may in addition offer an improved estimate on the ensemble-averaged MDS size. The possible deep connections between core percolation and the complexity of the random MDS problem will also be addressed by adapting the long-range frustration theory [24, 25].

The methods of this work can be readily extended to the MDS problem of directed networks. Our theoretical and algorithmic results on the directed MDS problem will soon be reported in an accompanying paper [47]. A more challenging problem is the connected dominating set problem [48] which has the additional constraint that the nodes in the dominating set should induce a connected subnetwork. Our present work may stimulate further theoretical studies on this hard problem.

Acknowledgments

Part of this work was done when H.-J. Zhou was participating in the “Collective Dynamics in Information Systems 2014” Program of the Kavli Institute for Theoretical Physics China (KITPC). H.-J. Zhou thanks Chuang Wang for a helpful discussion, and Alfredo Braunstein, Yang-Yu Liu, Federico Ricci-Tersenghi, and Yi-Fan Sun for helpful comments on the manuscript; J.-H. Zhao and H.-J. Zhou thank Prof. Zhong-Can Ou-Yang for support. Research partially supported by the National Basic Research Program of China (grant number 2013CB932804) and by the National Natural Science Foundation of China (grand numbers 11121403 and 11225526).

Copyright information

© Springer Science+Business Media New York 2015