International Journal of Statistical Mechanics

Volume 2014 (2014), Article ID 136829, 13 pages

http://dx.doi.org/10.1155/2014/136829

Research Article

The Statistical Mechanics of Random Set Packing and a Generalization of the Karp-Sipser Algorithm

C. Lucibello¹ and F. Ricci-Tersenghi²

¹Dipartimento di Fisica, Università “La Sapienza”, Piazzale Aldo Moro 2, 00185 Rome, Italy

²Dipartimento di Fisica, INFN-Sezione di Roma1, CNR-IPCF UOS Roma Kerberos, Università “La Sapienza”, Piazzale Aldo Moro 2, 00185 Rome, Italy

Received 19 November 2013; Accepted 8 January 2014; Published 10 March 2014

Academic Editor: Hyunggyu Park

Copyright © 2014 C. Lucibello and F. Ricci-Tersenghi. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

We analyse the asymptotic behaviour of random instances of the maximum set packing (MSP) optimization problem, also known as maximum matching or maximum strong independent set on hypergraphs. We give an analytic prediction of the MSPs size using the 1RSB cavity method from statistical mechanics of disordered systems. We also propose a heuristic algorithm, a generalization of the celebrated Karp-Sipser one, which allows us to rigorously prove that the replica symmetric cavity method prediction is exact for certain problem ensembles and breaks down when a core survives the leaf removal process. The -phenomena threshold discovered by Karp and Sipser, marking the onset of core emergence and of replica symmetry breaking, is elegantly generalized to for one of the ensembles considered, where is the size of the sets.

1. Introduction

The maximum set packing is a very much studied problem in combinatorial optimization, one of Karp’s twenty-one NP-complete problems. Given a set and a collection of its subsets labeled by ; a set packing (SP) is a collection of the subsets such that they are pairwise disjoint. The size of a SP is . A maximum set packing (MSP) is an SP of maximum size. The integer programming formulation of the MSP problem reads

The MSP problem, also known in the literature as the matching problem on hypergraphs or the strong independent set problem on hypergraphs, is an NP-Hard problem. This general formulation, however, can be specialized to obtain two other famous optimization problems: the restriction of the MSP problem to sets of size 2 corresponds to the problem of maximum matching on ordinary graphs and can be solved in polynomial time [1]; the restriction where each element of appears exactly 2 times in is the maximum independent set and belongs to the NP-Hard class.

The formulation ((1)–(3)) of the MSP problem, therefore, encodes an ample class of packing problems and, as all packing problems, is related by duality to a covering problem, the minimum set covering problem. Another common specialization of the general MSP problem, known as -set packing, is that in which all sets have size at most . This is one of the most studied specializations in the computer science community, the efforts concentrating on minimal degree conditions to obtain a perfect matching [2], linear relaxations [3, 4], and approximability conditions [5–7]. Motivated by this interest, we choose a -set packing problem ensemble as the principal application of the general analytical framework developed in the following sections. The asymptotic behaviour of random sparse instances of the MSP problem has not been investigated by mathematicians and computer scientists; only in the matching [8] and independent set [9] restrictions some work has been done. Extending some theorems of [8] (on which a part of this work is greatly inspired) to a greater class of problem ensembles is some of the main aims of the present work.

On the other hand also the statistical physics literature is lacking an accurate study of the random MSP problem. One of its specialization though, the matching problem, has been covered since the beginning of the physicists’ interest in optimization problems, with the work of Parisi and Mézard on the weighted and fully connected version of the problem [10, 11]. More recently the matching problem on sparse random graphs has also been accurately studied [12, 13] using the cavity method technique. Also the independent set problem on random graphs [14] and the dual problem to set packing, the set covering problem [15], received some attention by the disordered physics community. The SP problem was investigated with the cavity method formalism in a disguised form, as a glass model on a generalized Bethe lattice, in [16, 17]. This corresponds, as we will see in the next section, to a factor graph ensemble with fixed factor and variable degrees; thus, we will not cover this case in Section 7.

The paper is organized as follows.

(i)In Section 2 we map the MSP problem ((1)–(3)) into a statistical physical model defined on a factor graph and relate the MSP size to the density at infinite chemical potential.(ii)We introduce the replica symmetric (RS) cavity method in Section 3 and give an estimate for the average MSPs size on sparse factor graph ensembles in the thermodynamic limit.(iii)In Section 4 we establish a criterion for the validity of the RS ansatz and introduce the 1RSB formalism.(iv)In Section 5 we propose a generalization of the Karp-Sipser heuristic algorithm [8] to the MSP problem and prove the validity of the RS ansatz for certain ensemble of problems. Moreover we find a relationship between a core emergence phenomena and the breaking of replica symmetry breaking.(v)Section 6 describes the numerical simulations performed.(vi)In Section 7 we apply the analytical tools developed to some problem ensembles. We compare the numerical results obtained from an exact algorithm with the analytical predictions, focusing to greater extent to one ensemble modelling the -set packing.

2. Statistical Physics Description

In order to turn the MSP combinatorial optimization problem into a useful statistical physical model let us recast ((1)–(3)) into a graphical model using the factor graph formalism [18, 19]. We define our variable nodes set to be and to each we associate a variable taking values in as in (3). will be our factor nodes set, as its elements acts as hard constrains on the variables through (2). The edge set is then naturally defined as . We call the factor graph thus composed and can then rewrite (2) as

A SP is a configuration satisfying (4) and its relative size is that is simply the fraction of occupied sites.

It is now easy to define an appropriate Gibbs measure for the MSPs problem on through the grand canonical partition function

Only SPs contribute to the partition function, and in the close packing limit, as we will call the limit , the measure is dominated by MSPs. Equation (5) is also a model for a particle gas with hard core repulsion and chemical potential located on a hypergraph and as such has been studied mainly on lattice structures and more in general on ordinary graphs. Model (5) has been studied on a generalized Bethe lattice (i.e., the ensemble defined in Section 7, a -uniform -regular factor graph) in [16, 17] as a prototype of a system with finite connectivity showing a glassy behaviour. This has been the only approach, although disguised as a hard spheres model, from the statistical physics community to a general MSP problem.

The grand canonical potential is defined as

and the particle density as

Grand potential and density are related to entropy by the thermodynamic relation

In the close packing limit (i.e., we recover the MSP problem, since in this limit the Gibbs measure is uniformly concentrated on MSPs and gives the MSP relative size. Since entropy remains finite in this limit, from (8) we obtain the MSP relative size

In this paper we focus on random instances of the MSP problem. As usual in statistical physics we will assume the number of variables (the number of subsets ) and the number of constrains to diverge, keeping the ratio finite. We will refer to this limit as the thermodynamic limit. Instances of the MSP problem will be encoded in factor graph ensembles which we assume to be locally tree-like in the thermodynamic limit.

The MSP relative size is a self-averaging quantity in the thermodynamic limit and we want to compute its asymptotic value

where we denoted with the expectation over the factor graph ensemble. In the last equation the dependence is encoded in the graph ensemble considered. Computing (10) is not an easy task and some approximation have to be taken. We will employ the cavity method from the statistical physics of disordered systems [19, 20], using both the replica symmetric (RS) and the one-step replica symmetry breaking (1RSB) ansatz. We will prove in Section 5 that the RS ansatz is exact in a certain region of the phase space, while in Section 7 we will give numerical evidence that the 1RSB approximation gives very good results outside the RS region.

3. Replica Symmetry

3.1. Bethe Approximation on a Single Instance

The RS cavity method has been known for many decades outside the statistical physics community as the Belief Propagation (BP) algorithm and only in recent years the two approaches have been bridged [18, 19]. We start with a variational approximation to the grand potential equation (6) of an instance of the problem, the Bethe free energy approximation:

with the factor and variable contributions given by

The grand canonical potential is expressed as a function of the factor node to variable node messages . Minimization of over the messages constrained to be normalized to one yields the fixed point BP equations for the set packing:

The coeeficients are normalization factor. Equations (13) can be simplified introducing the fields defined as

yielding

Since we are interested in the close packing limit to solve the problem we will straightforwardly apply the zero temperature cavity method [21]. The related BP equations which can be found as the limit of (15) read

We note that the messages are bounded to take values in the interval and that if we set the initial value of each in the discrete set , at each BP iteration, all messages will take value either 0 or 1. These values can be directly interpreted as the occupational loss occurring in the subtree if the subtree is connected to the occupied node . This loss cannot be negative (thus a gain), since we put an additional constrain on the subtree demanding every neighbour of to be empty, and cannot be greater than 1 as well, in fact corresponds to the worst case scenario where an otherwise occupied node has to be emptied.

The Bethe free energy for model (5) on the factor graph can be expressed as a function of the fixed point messages as

We finally arrive to the Bethe estimation of the MSPs relative size, taking the close packing limit of (17) and using , which is given by

Let us examine the various contributions to (18) since we want to convince ourselves that it exactly counts the MSP size, at least on tree factor graphs. The term contributes with to the sum only if all the incoming messages are zero. In this case is frozen to 1, that is, the variable takes part of all the MSPs in . Obviously all the neighbours of a variable frozen to 1 have to be frozen to 0. To all its neighbours , the frozen to 1 variable sends a message , so that we have contributions in the first sum of (18) and the total contribution from correctly sums up to 1. If for a certain we have a total field (two or more incoming messages are equal to one) the variable is frozen to 0; that is, it does not take part of any MSPs. It correctly does not contribute to since it sends a message to each neighbour ; thus, it is not computed in .

The third case is the most interesting. It concerns variables which take part to a fraction of the MSPs. We will call them unfrozen variables. The total field on an unfrozen variable is (thus we have no contribution from the second sum in (18)) and all incoming messages are 0 except for a single . To this sole function node , the node sends a message 1, so that the contribution of to the first sum is . Actually BP equations impose that has to have at least another unfrozen neighbour beside . In other terms the function node says that whatever MSP we consider, one of my neighbours has to be occupied. The corresponding term in (18) accounts for that.

The presence of unfrozen variables is the reason why we cannot express the density through the formula

In fact, using the infinite chemical potential formalism, we cannot compute

when , and we would have to use the corrections to the fields . We bypass the problem using the grand potential to obtain , also addressing a problem reported in [22] of extending an analysis suited for weighted matchings and independent sets to the unweighted case.

3.2. Ensemble Averages

To proceed further in the analysis and since one is often concerned with the average properties of a class of related factor graphs, let us consider the case where the factor graph is sampled from a locally tree-like factor graph ensemble . We will employ the following notation for the graph ensembles expectations: for graphs averages; () for expectations over the factor (variable) degree distribution, which we will sometime call root degree distribution; and () for expectations over the excess degree distribution of factor (variable) nodes conditioned to have at least one adjacent edge, which we will sometime call residual degree distribution. The quantities and and the random variables , , , and , are related by

In Section 7 we discuss some specific factor graph ensembles, where and are fixed to a deterministic value or Poissonian distributed.

With these definitions the distributional equation corresponding to Belief Propagation formula (16) which reads

where are i.i.d. random residual variable degrees, is a random residual factor degree, and are i.i.d. random incoming messages. Since on a given graph each message takes value only on we can take the distribution of messages to be of the form

For this to be a fixed point of (22), the parameter has to satisfy the self-consistent equation

Using (18) and (24), we obtain the replica symmetric approximation to the asymptotic MSP relative size

It turns out that the RS approximation is exact only when the ratio is sufficiently low. In the SP language (25) holds true only when the number of the subsets among which we choose our MSP is not too big compared to the number of elements of which they are composed.

In the next section we will quantitatively establish the limits of validity of the RS ansatz and introduce the one-step replica symmetry breaking formalism which provides a better approximation to the exact results in the regime of large ratio. In Section 5 we prove that (25) is exact for certain choices of the factor graph ensembles.

4. Replica Symmetry Breaking

4.1. RS Consistency and Bugs Propagation

Here we propose two criterions in order to check the consistency of the RS cavity method, iIf any of those fails.

The first criterion is the assumption of unicity of the fixed point of (24) and its dynamical stability under iteration. We restrict ourselves to the subspace of distributions with support on , although it is possible to extend the following analysis to the whole space of distributions over with an argument based on stochastic dominance following [22]. Characterizing the distributions over as in (23) with a real parameter , from (22) we obtain the dynamical system

The stability criterion

suggests that the RS approximation to the MSP size (25) is exact as long as (27) is satisfied. This statement will be made rigorous in Section 5.

The second method we use to check the RS stability, called bugs proliferation, is the zero temperature analogous of spin glass susceptibility. We will compute the average number of changing messages induced by a change in a single message ( or ). This is given by

where , , , and take value in . Since a random factor graph is locally a tree, and assuming correlations decay fast enough, last equation can be expressed as

where is a message at distance from the tree root , and we defined the average residual degrees and . The stability condition yields a constraint on the greatest eigenvalue of the transfer matrix :

The two methods presented above give equivalent conditions for the RS ansatz to hold true and they simply express the independence for a finite subgraph from the tail boundary conditions.

4.2. The 1RSB Formalism

We are going to develop the 1RSB formalism for the MSP problem and then apply it in Section 7.1 to the ensemble . We will not check the coherence of the 1RSB ansatz through the interstate and intrastate susceptibilities [23]; we are then not guaranteed against the need of further steps of replica symmetry breaking in order to recover the exact solution. Even in the worst case scenario though, when the 1RSB solution is trivially exact only in the RS region and a full RSB ansatz is needed otherwise, the 1RSB prediction for MSP relative size should be everywhere more accurate than the RS one and possibly very close to the real value. We will refer to the textbook of Montanari and Mézard [19] for a detailed exposition of the 1RSB cavity method.

Let us fix a factor graph from a locally tree-like ensemble . We call the distribution of messages on the directed edge over the states of the system. We still expect the messages to take values 0 or 1, so that can be parametrized as

The 1RSB Parisi parameter has to be properly rescaled in order to correctly take the limit . Therefore we introduce the new 1RSB parameter which stays finite in the close packing limit and takes a value in . The reweighting factor is defined as

Last equation combined with (13) and (14) gives . Averaging over the whole ensemble we can then write the zero temperature 1RSB message passing rules (also called Survey Propagation equations):

In preceding equation are i.i.d.r.v. on and, as usual, is the random variable residual degree and are random independent factors residual degrees. Fixed points of (33) take the form

where is a continuous distribution on and . Parameters and have to satisfy the closed equations

Solutions of (35) with correspond to replica symmetric solutions and their instability marks the onset of a spin glass phase. In this new phase the MSPs are clustered according to the general scenario displayed by constraint satisfaction problems [24].

From the stable fixed point of (33) we can calculate the 1RSB free energy functional as

and 1RSB density, , as

with expectations intended over and over fixed point messages and . Since the free energy functional and the complexity are related by the Legendre transform

with , through (36) we can compute the complexity taking the inverse transform. Equilibrium states, that is, MSPs, are selected by

or equivalently taking the 1RSB parameter to be

In the static 1RSB phase we expect so that from (38) we have . We will see that this is generally true except for the ensemble of Section 7, corresponding to maximum matchings on ordinary graphs, where the equilibrium state have maximal complexity and . The relation

is always valid though, since for complexity stays finite.

5. A Heuristic Algorithm and Exact Results

In this section we propose a heuristic greedy algorithm to address the problem of MSP. It is a natural generalization of the algorithm that Karp and Sipser proposed to solve the maximum matching problem on Erdös-Rényi random graphs [8]; therefore, we will call it generalized Karp-Sipser (GKS). Extending their derivation concerning the leaf removal part of the algorithm we are able to prove that the RS prediction for MSP density is exact as long as the stability criterion (27) is satisfied. We will not give the proofs of the following theorems as they are lengthy but effortless extension of those given in [8]. In order to find the maximum matching on a graph Karp and Sipser noticed that as long as the graph contains a node of degree one (a leaf), its unique edge has to belong to one of the perfect matchings.

They considered the simplest randomized algorithm one can imagine: as long as there is any leaf remove it from the graph, otherwise remove a random edge; then iterate until the graph is depleted. They studied the behaviour of this leaf-removal algorithm on random graphs and were able to prove that it grants w.h.p a maximum matching (within an error).

To generalize some of their results we need to extend the definition of leaf to that of pendant. We call pendant a variable node whose factor neighbours all have degree one, except for one at most. Stating the same concept in different words, all of the neighbours of a pendant have the pendant itself as their sole neighbour, except for one of them at most. See Figure 1 for a pictorial representation of a pendant (in red). The GKS algorithm is articulated in two phases: a pendant removal and a random occupation phase. We give the pseudocode for the generalized Karp-Sipser algorithm (see Algorithm 1).

Require: a factor graph (, , )

Ensure: a set packing

 add to all isolated variable nodes and remove them

 from

 remove from any isolated factor node

 while   is not empty do

      if   has any pendant then

         choose a pendant uniformly at random

         add to

         remove from , then remove its factor neighbours

         and their variable neighbours

      else

         pick uniformly at random a variable node and add

         it to , remove it from , then remove its factor

         neighbours and their variable neighbours

      end if

      remove from any isolated factor node

end while

Return  

Algorithm 1: Algorithm  1 generalized Karp-Sipser (GKS).

Figure 1: (left) A representation of a factor graph containing a pendant, depicted in red. (right) The factor graph after the removal of the pendant as occurs in the inner part of the while loop in Algorithm 1.

At each step the algorithm prioritizes the removal of pendants over that of random variable nodes. We notice that the removal of a pendant is always an optimal choice in order to achieve a MSP; we have no guarantees though on the effect of the occupation of a random node. We call phase 1 the execution of the algorithm up to the point where the first nonpendant variable is added to . It is trivial to show that phase 1 is enough to find an MSP on a tree factor graph. The interesting thing though is that phase 1 is also able to deplete nontree factor graphs and find an MSP as long as the factor graphs are sufficiently sparse and large enough.

We call core the subset of the variable nodes which has not been assigned to the MSP in phase 1. We will show how the emergence of a core is directly related to replica symmetry breaking. Let us define as usual

The function is continuous, nonincreasing, and satisfies the relation . It turns out that, in the large graph limit, phase 1 of GKS is characterized by the solutions of the system of equations

We notice that system (43) is equivalent to 1RSB equation (35) once the substitutions and are made.

Lemma 1. The system of (43) admits always a (unique) solution .

Proof. By monotony and continuity the function has a single zero in the interval .

If other solutions are present the relevant one is the one with smallest , as the following theorems certify.

Theorem 2. Let be the solution with smallest of (43). Then the density of factor nodes surviving in the large graphs limit is given by

with

In particular if the smallest solution is the graph is depleted with high probability in .

As already said we will not give the proof of this and of the following theorem, as they are lengthy and can be obtained from the derivation given in Karp and Sipser’s article [8] with little effort even if not in a completely trivial way. Theorem 2 affirms that as soon as (43) develops another solution a core survives phase 1. This phenomena coincides with the need for replica symmetry breaking in the cavity formalism. Let us now establish the exactness of the cavity prediction for in the RS phase.

Theorem 3. Let be the solution with smallest of (43). Then the density of variables assigned in in the large graphs limit is

where has been defined in previous theorem. In particular if the smallest solution is the replica symmetric cavity method prediction (25) is exact.

These two theorems imply that the RS prediction for , (25), holds true as long as phase 1 manages to deplete w.h.p. the whole factor graph. A similar behaviour has been observed in other combinatorial optimization problem, for example, the random XORSAT [25]. Conversely it is easy to prove, given the equivalence between (43) and (35), that if a solution with exists the RS fixed point is unstable in the 1RSB distributional space. Therefore the system is in the RS phase if and only if (35) admits a unique solution.

We notice that we could have also looked for the solutions of the single equation instead of the two of (43), since the value of is uniquely determined by the value of . Then previous theorems state that the RS results hold as long as has a single fixed point. In [22] it has been proven that the RS solution of the weighted maximum independent set and maximum weighted matching holds true as long as the corresponding squared cavity operator has a unique fixed point. Since those are special cases of the weighted MSP problem, we conjecture that the condition holds for the general case too, as it does for the unweighted case. Probably it would be possible to further generalize the Karp-Sipser algorithm to cover analytically the weighted case as well.

The scenario arising from the analysis of the algorithm and from 1RSB cavity method is the following: in the RS phase maximum set packings form a single cluster and it is possible to connect any two of them with a path involving MSPs separated by a rearrangement of a finite number of variables [26]. In the RSB phase instead MSPs can be grouped into many connected (in the sense we mentioned before) clusters, each one of them defined by the assigning of the variables in the core. Two MSPs who differ on the core are always separated by a global rearrangement of the variables (i.e., ). In presence of a core the GKS algorithm makes some suboptimal random choices (after phase 1).

We did not make an analytical study of the random variable removal part of the GKS algorithm. This has been done by Wormald for the matching problem, associating to the graph process a set of differential equations amenable to analysis [9] (something similar has been done for random XORSAT as well [25]). In that case it turns out that the algorithm still achieves optimality and yields an almost perfect matching on the core. An appropriate analysis of the second phase of the GKS algorithm and the reason of its failure for most of the SP ensembles we covered deserve further studies.

6. Numerical Simulations

While the RS cavity (24) and (25) can be easily computed, the numerical solution of the 1RSB density evolution (33) is much more involved (although simplified by the zero temperature limit) and has been obtained with a standard population dynamics algorithm [27].

The implementation of the GKS algorithm is straightforward. While it is very fast during phase 1, we noticed a huge slowing down in the random removal part. We were able to find set packings for factor graphs with hundred of thousands of nodes.

In order to test the cavity and GKS predictions, we also computed the exact MSP size on factor graphs of small size. First we notice that a set packing problem, coded in a factor graph structure , is equivalent to an independent set problem on an appropriate ordinary graph . The node set of will be the variable node set of and we add an edge to between each pair of nodes having a common factor node neighbour in . Therefore each neighborhood of a factor node in forms a clique in . We then solve the maximum set problem for using an exact algorithm [28] implemented in the igraph library [29]. Since the time complexity is exponential in the size of the graph we performed our simulations on graphs containing up to only one hundred nodes.

7. Applications to Problem Ensembles

We will now apply the methods developed in the previous sections to some factor graph ensembles, each modelling a class of MSP problem instances. We consider graph ensembles containing nodes with Poissonian random degree or regular degree: , , , and . Subscript or indicates whether the type of nodes to which they refer (variable nodes for the first subscript, factor nodes for the second) have regular or Poissonian random degree, respectively. We parametrize these ensemble by their average variable and factor degree; that is, and . They are constituted by factor graphs having variable nodes and factor nodes but they differ both for their elements and for their probability law. We define the ensembles giving the probability of sampling one of their elements.

(i): each element has variables and factors. Every variable node has fixed degree . is obtained linking each variable with factors chosen uniformly and independently at random. The factor nodes degree distribution obtained is Poissonian of mean with high probability. This ensemble is a model for the -set packing and will be the main focus of our attention.(ii): each element has variables and factor. Every factor node has fixed degree . is built linking each factors with variables chosen uniformly and independently at random. The variable nodes degree distribution obtained is Poissonian of mean with high probability.(iii): each element has variables and factors. is built adding an edge with probability independently for each choice of a variable and a factor . The factor graph obtained has w.h.p Poissonian variable nodes degree distribution of mean and Poissonian factor nodes degree distribution of mean .(iv): it is constituted of all factor graphs of variable nodes of degree and factor nodes of degree ( has to be multiple of ). Every factor graph of the ensemble is equiprobable and can be sampled using a generalization of the configuration model for random regular graph [30]. This ensemble is a model for the -set packing.

We will omit the argument when we refer to an ensemble in the limit.

7.1.

This is the ensemble with variable node degrees fixed to and Poissonian factor node degree that is . The case corresponds to the maximum matching problem on Erdös-Rényi random graph [8, 12].

Density evolution (22) for reduces to

Considering distributions of the form fixed points of (47) reads

The last equation admits only one solution for each value of and , as it is easily seen through a monotony argument considering the left and right hand side of the equation. The values of as a function of satisfying are the critical points delimiting the RS phase () and are given by

For a core survives phase 1 of the GKS algorithm as stated by Theorem 2 and shown in Figure 2. We notice that the same threshold value is found in the minimum set covering problem on the dual factor graph ensemble [31], where the application of a leaf removal algorithm, that is, phase 1 of GKS, unveils an analogous core emergence phenomena. In the RS phase, the MSP density is equal to the minimum set covering density, but in the RSB phase we cannot compare the cavity predictions for these two dual problems, since in [31] the 1RSB solution to the minimum set covering is not provided.

Figure 2: Results from GKS algorithm applied to . The set packing size after phase 1 (blue line) and after the algorithm stops (red line) are reported. The vertical line is and is drawn as a visual aid to determine the point where the RS assumptions stop to hold true. Above a nonempty core emerges continuously (green line).

In the matching case, that is, equation (49) expresses the notorious -phenomena discovered by Karp and Sipser, while for higher values of it provides an extension of the critical threshold.

We can recover the same critical condition equation (49) through the bug propagation method, as the transfer matrix non-zero elements are only:

which give . The average branching factor is so that (30) and (48) yield (49). The analytical value for the relative size of MSPs, that is, the particle density , is

In Figure 3 we compared from (51) as a function of for some values with an exact algorithm applied to finite factor graphs (as explained in Section 6), both above and below . Clearly for the RS approximation is increasingly more inaccurate.

Figure 3: RS cavity method analytical prediction equation (51) for MSP size in compared with the average size of MSPs obtained from 1000 samples of factor graphs with 100 variable nodes (dots). Dashed parts of the lines are the RS estimations in the RSB phase, that is, for , where it is no longer exact.

We continue our analysis of above the critical value through the 1RSB cavity method as outlined in Section 4.2. Fixed point messages of (33) are distributed as

where

and has to be determined through (33). Equation (53) admits always an RS solution which is stable up to , as already noticed. Above a new stable fixed point, with , continuously arises and we study it numerically with a population dynamics algorithm.

The 1RSB free energy, as a function of the Parisi parameter , takes the form

As prescribed by the cavity method, the value which maximizes over yields the correct free energy; therefore, we have .

Unsurprisingly, as they belong to different computational classes, the cases and show qualitatively different pictures. In the case of maximum matching on the Poissonian graph ensemble, numerical estimates suggest that complexity is an increasing function of on the whole real positive axis. Correct choice for parameter is then , as already conjectured in [12], and we find that maximum matching size prediction from 1RSB cavity method fully agrees with rigorous results from Karp and Sipser [8] and with the size of the matchings given by their algorithm (see Figure 4). The 1RSB ansatz is therefore exact for .

Figure 4: Maximum set packing size for as a function of for . Blue line corresponds to the RS analytical value (51), red line is from the Karp-Sipser algorithm and dotted black line is the 1RSB solution (54) with . Karp-Sipser and 1RSB cavity method yield the same and exact result.

The case analysis does not yield such a definite result. The complexity of states is no more a strictly increasing function of . It reaches its maximum in , the choice of that selects the most numerous states, which could be those where local search greedy algorithms are more likely to be trapped. Then it decreases up to the finite value where complexity changes sign and takes negative values. Therefore is the correct choice for the Parisi parameter which maximizes . Plotted as a function of , complexity has a convex nonphysical part, with extrema the RS solution (on the right) and the point corresponding to the dynamic 1RSB solution (on the left) and a concave physically relevant for (see Figure 5). The 1RSB seems to be in very good agreement with the exact algorithm and we are inclined to believe that no further steps of replica symmetry breaking are needed in this ensemble. The GKS algorithm instead falls short of the exact value (see Figure 6); therefore, it constitutes a lower bound which is not strict but at least it could probably be made rigorous carrying on the analysis of the GKS algorithm beyond phase 1.

Figure 5: Complexity in , with and , as a function of the relative MSP size .

Figure 6: MSP size in as a function of . RS and 1RSB cavity method is confronted with the GKS algorithm and with an exact algorithm on factor graphs of 100 variable nodes.

7.2.

The ensemble is constitute factor graphs containing factor nodes of degree fixed to and variable nodes of Poissonian random degree of mean . It has statistical properties with respect to the MSP problem quite different from those encountered in , as we will readily show. The MSP problem on with is equivalent to the well known problem of maximum independent sets on Poissonian graphs [14, 32]. The real parameter characterizing discrete support distributions of messages has to satisfy the fixed point (24) that reads

At fixed value of the RS ansatz holds up to the critical value , which is implicitly given by first derivative condition (27):

Although critical condition equation (56) is not as elegant as the one we obtained for the ensemble , it can be easily solved numerically for as a function of . For the threshold value is exactly . For instead is an increasing function of the factor degree . Thanks to (25) we readily compute the MSP size in the RS phase:

We can see in Figure 7 that the MSP size is a decreasing function in both arguments as expected. Equation (57) can be taken as the RS estimate for MSP size for . The RS estimate is strictly greater than the average size of SPs given by the GKS algorithm at all values of (see Figure 7).

Figure 7: RS cavity method analytical prediction (57) for MSP size in (dot-dashed) compared with the set packing size from phase 1 of the GKS algorithm (dashed lines) and with a complete run of the algorithm (continuous lines).

7.3.

We will now briefly examine our MSP model on the ensemble where both factor nodes and variable nodes have Poissonian random degrees of mean and , respectively. From (23) and (24) we obtain the fixed point condition for the probability distribution of messages:

As usual is the parameter characterizing the distribution of messages, . Equation (58) admits one and only one fixed point solution for each value of and . In fact is continuous, strictly decreasing, and . The first derivative condition defines the critical line through

The curve separates the RS phase from the RSB phase in the parametric space (see Figure 8). The unbounded RS region shares some resemblance with the corresponding (although -discretized) region of and is at variance with the compact area of the RS phase in .

Figure 8: MSP relative size for the ensemble as given by (60). Above the boundary separating RS and RSB phase values are only approximation to the exact value of .

We can compute the MSP relative size through (25) and obtain

which holds only as an approximation in the RSB phase. We can see from Figure 8 that is decreasing both in and as was observed in the other ensembles as well.

7.4.

The MSP problem on the ensemble , the straightforward generalization to factor graphs of the random regular graphs ensemble, poses some simplification to the cavity formalism thanks to his homogeneity. It has already been the object of preliminary studies by Weigt and Hartmann [16] and then a much more deep work of Hansen-Goos and Weigt [17] who disguised it as a hard spheres model on a generalized Bethe lattice. The authors studied through the cavity method this hard spheres model on both at finite chemical potential and in the close packing limit and found out that the 1RSB solution is unstable in the close packing limit, therefore suggesting the need of a full RSB treatment of the problem.

8. Conclusions

We studied the average asymptotic behaviour of random instances of the maximum set packing problem, both from a mathematical and a physical viewpoint. We contributed to the known list of models where the replica symmetric cavity method can be proven to give exact results, thanks to the generalization of an algorithm (and of its analysis) first proposed by Karp and Sipser [8]. Moreover, our analysis address a problem reported in recent work on weighted maximum matchings and independent sets on random graphs [22], where the authors could not extend their results to the unweighted cases. We achieve here the desired result making use the grand canonical potential instead of the direct computation of single variable expectations. We also extend their condition for the system to be in what physicists call a replica symmetric phase, namely, the uniqueness of the fixed point of the square of a certain operator (which is the analogue of the one defined in (42)), to the more general setting of maximum set packing (although without weights). On some problem ensembles, where the assumptions of Theorems 2 and 3 no longer hold and the RS cavity method fails, we used the 1RSB cavity method machinery to obtain an analytical estimation of the MSP size. Numerical simulations show very good agreement of the 1RSB estimation with the exact values, although comparisons have been done only with small random problems due to the exact algorithm being of exponential time complexity. The GKS algorithm instead generally fails to find any MSPs except for some special instances of the problem.

Some questions remain open for further investigation. To validate the 1RSB approach the stability of the 1RSB solution has to be checked against more steps of replica symmetry breaking. Moreover a thorough analysis of the second phase of the GKS algorithm could shed some light on the mechanism of replica symmetry breaking and give a rigorous lower bound to the average maximum packing size.

Regarding the problem of looking for optimal solutions in single samples, we believe that an efficient heuristic algorithm, able to obtain near-to-optimal configurations also in the RSB phase, could be constructed following the ideas of [33].

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This research has received financial support from the European Research Council (ERC) through Grant agreement no. 247328 and from the Italian Research Minister through the FIRB project no. RBFR086NN1.

References

S. Micali and V. V. Vazirani, “An O algoithm for finding maximum matching in general graphs,” in Proceedings of the 21st Annual IEEE Symposium on Foundations of Computer Science, pp. 17–27, 1980.

N. Alon, J.-H. Kim, and J. Spencer, “Nearly perfect matchings in regular simple hypergraphs,” Israel Journal of Mathematics, vol. 100, pp. 171–187, 1997.

Z. Füredi, J. Kahn, and P. D. Seymu, “On the fractional matching polytope of a hypergraph,” Combinatorica, vol. 13, no. 2, pp. 167–180, 1993.

Y. H. Chan and L. C. Lau, “On linear and semidefinite programming relaxations for hypergraph matching,” Mathematical Programming, vol. 135, pp. 123–148, 2011.

V. Rödl and L. Thoma, “Asymptotic packing and the random greedy algorithm,” Random Structures & Algorithms, vol. 8, no. 3, pp. 161–177, 1996.

M. M. Halldórsson, “Approximations of weighted independent set and hereditary subset problems,” Journal of Graph Algorithms and Applications, vol. 4, no. 1, pp. 1–6, 2000.

E. Hazan, S. Safra, and O. Schwartz, “On the complexity of approximating K-set packing,” Computational Complexity, vol. 15, no. 1, pp. 20–39, 2006.

R. M. Karp and M. Sipser, “Maximum matchings in sparse random graphs,” in Proceedings of the 22nd Annual Symposium on Foundations of Computer Science, vol. 12, pp. 364–375, IEEE Computer Society, 1981.

N. C. Wormald, “The differential equation method for random graph processes and greedy algorithms,” in Lectures on Approximation and Randomized Algorithms, pp. 73–155, 1999.

G. Parisi and M. Mézard, “Replicas and optimization,” Journal de Physique Lettres, vol. 46, no. 17, pp. 771–778, 1985.

G. Parisi and M. Mézard, “A replica analysis of the travelling salesman problem,” Journal de Physique, vol. 47, pp. 1285–1296, 1986.

H. Zhou and Z.-C. Ou-Yang, “Maximum matching on random graphs,” In press, http://arxiv.org/abs/cond-mat/0309348

L. Zdeborová and M. Mézard, “The number of matchings in random graphs,” Journal of Statistical Mechanics, vol. 2006, Article ID P05003, 2006.

L. Dallasta, A. Ramezanpour, and P. Pin, “Statistical mechanics of maximal independent sets,” Physical Review E, vol. 80, Article ID 061136, 2009.

M. Mézard and M. Tarzia, “Statistical mechanics of the hitting set problem,” Physical Review E, vol. 76, no. 4, Article ID 041124, 10 pages, 2007.

M. Weigt and A. K. Hartmann, “Glassy behavior induced by geometrical frustration in a hard-core lattice-gas model,” Europhysics Letters (EPL), vol. 62, Article ID 021005, p. 533, 2003.

H. Hansen-Goos and M. Weigt, “A hard-sphere model on generalized bethe lattices: statics,” Journal of Statistical Mechanics, vol. 2005, Article ID P04006, 2005.

F. R. Kschischang, B. J. Frey, and H.-A. Loeliger, “Factor graphs and the sum-product algorithm,” IEEE Transactions on Information Theory, vol. 47, no. 2, pp. 498–519, 2001.

A. Montanari and M. Mézard, Information, Physics and Computation, Oxford University Press, 2009.

G. Parisi, M. Mézard, and M. A. Virasoro, Spin Glass Theory and Beyond, World Scientific, Singapore, 1987.

M. Mézard and G. Parisi, “The cavity method at zero temperature,” Journal of Statistical Physics, vol. 22, no. 1-2, pp. 1–34, 2003.

D. Gamarnik, T. Nowicki, and G. Swirszcz, “Maximum weight independent sets and matchings in sparse random graphs. Exact results using the local weak convergence method,” Random Structures and Algorithms, vol. 28, no. 1, pp. 76–106, 2006.

A. Montanari, G. Parisi, and F. Ricci-Tersenghi, “Instability of one-step replica-symmetry-broken phase in satisfiability problems,” Journal of Physics A, vol. 37, no. 6, p. 2073, 2004.

F. Krzakala, A. Montanari, F. Ricci-Tersenghi, G. Semerjian, and L. Zdeborová, “Gibbs states and the set of solutions of random constraint satisfaction problems,” Proceedings of the National Academy of Sciences of the United States of America, vol. 104, no. 25, pp. 10318–10323, 2007.

M. Mézard, F. Ricci-Tersenghi, and R. Zecchina, “Two solutions to diluted p-spin models and XORSAT problems,” Journal of Statistical Physics, vol. 111, no. 3-4, pp. 505–533, 2002.

G. Semerjian, “On the freezing of variables in random constraint satisfaction problems,” Journal of Statistical Physics, vol. 130, no. 2, pp. 251–293, 2008.

M. Mézard and G. Parisi, “The bethe lattice spin glass revisited,” The European Physical Journal B, vol. 20, no. 2, pp. 217–233, 2001.

S. Tsukiyama, M. Ide, H. Ariyoshi, and I. Shirakawa, “A new algorithm for generating all the maximal independent sets,” SIAM Journal on Computing, vol. 6, no. 3, pp. 505–517, 1977.

G. Csárdi and T. Nepusz, “The igraph software package for complex network research,” InterJournal, Complex Systems, p. 1695, 2006.

N. C. Wormald, “Models of random regular graphs,” in Surveys in Combinatorics, London Mathematical Society Lecture Note, pp. 239–298, 1999.

S. Takabe and K. Hukushima, “Minimum vertex cover problems on random hypergraphs: replica symmetric solution and a leaf removal algorithm,” In press, http://arxiv.org/abs/1301.5769

S. Sanghavi, D. Shah, and A. S. Willsky, “Message passing for maximum weight independent set,” IEEE Transactions on Information Theory, vol. 55, no. 11, pp. 4822–4834, 2009.

M. Weigt and H. Zhou, “Message passing for vertex covers,” Physical Review E, vol. 74, Article ID 046110, 19 pages, 2006.