The statistical mechanics of random set packing and a generalization of the Karp-Sipser algorithm

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 analytical 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 $e$-phenomena threshold discovered by Karp and Sipser, marking the onset of core emergence and of replica symmetry breaking, is elegantly generalized to $c_s = \frac{e}{d-1}$ for one of the ensembles considered, where $d$ is the size of the sets.


Introduction
The maximum set packing is a very much studied problem in combinatorial optimization, one of Karp's twenty-one NPcomplete problems. Given a set = {1, . . . , } and a collection of its subsets S = { | ⊆ , ∈ } labeled by = {1, . . . , }; a set packing (SP) is a collection of the subsets such that they are pairwise disjoint. The size of a SP S ⊆ S is |S |. A maximum set packing (MSP) is an SP of maximum size. The integer programming formulation of the MSP problem reads subject to ∑ : ∈ ≤ 1 ∀ ∈ , = {0, 1} ∀ ∈ .
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 S 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 asset 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][6][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 2 International Journal of Statistical Mechanics 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.

Statistical Physics Description
In order to turn the MSP combinatorial optimization problem into a useful statistical physical model let us recast ((1)- ) 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 {0, 1} 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 ({ }) = (1/ ) ∑ ∈ 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 G ( , ) defined in Section 7, a -uniformregular 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 International Journal of Statistical Mechanics 3 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 E [⋅] 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.

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 messageŝ= {] → }. Minimization of RS [̂] 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 aŝ 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 [0, 1] and that if we set the initial value of each → in the discrete set {0, 1}, 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 → = 1 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 4

International Journal of Statistical Mechanics
We finally arrive to the Bethe estimation of the MSPs relative size, taking the close packing limit of (17) and using RS = − RS , 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 (1 − | |) max{0, 1 − ∑ ∈ → } contributes with 1 − | | 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 1−∑ ∈ \ → = 1, so that we have | | contributions in the first sum of (18) max{0} ∪ {1 − ∑ ∈ \ → } ∈ = 1 and the total contribution from correctly sums up to 1. If for a certain we have a total field ≡ 1 − ∑ ∈ → < 0 (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 RS since it sends a message 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 = 0 (thus we have no contribution from the second sum in (18)) and all incoming messages are 0 except for a single → = 1. To this sole function node , the node sends a message 1, so that the contribution of to the first sum is 1. 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 max{0} ∪ {1 − ∑ ∈ \ → } ∈ = 1 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 lim → +∞ 1 − ∑ ∈ → = 0, and we would have to use the (1/ ) corrections to the fields { → }. We bypass the problem using the grand potential RS to obtain RS , also addressing a problem reported in [22] of extending an analysis suited for weighted matchings and independent sets to the unweighted case.

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 G( ). We will employ the following notation for the graph ensembles expectations: for expectations over the factor (variable) degree distribution, which we will sometime call root degree distribution; and E [⋅] (E [⋅]) 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 , 0 , , and 0 , 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) 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 {0, 1} 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 onestep 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.
International Journal of Statistical Mechanics

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 {0, 1}, although it is possible to extend the following analysis to the whole space of distributions over [0, 1] with an argument based on stochastic dominance following [22]. Characterizing the distributions over {0, 1} as in (23) with a real parameter ∈ [0, 1], 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 → (1 → 0 or 0 → 1). This is given by where 0 , 1 , 0 , and 1 take value in {0, 1}. 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 = E [ ] and = E [ ]. The stability condition ch < +∞ yields a constraint on the greatest eigenvalue of the transfer matrix ( 0 → 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.

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 G . 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 G. 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 ∈ [0,1] 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 [0, +∞). The reweighting factor − → iter is defined as Last equation combined with (13) and (14) gives → iter = − → . 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 [0, 1] and, as usual, is the random variable residual degree and { } are random independent factors residual degrees. Fixed points of (33) take the form where 2 ( ) is a continuous distribution on [0, 1] and 2 = 1 − 0 − 1 . Parameters 0 and 1 have to satisfy the closed equations Solutions of (35) with 2 = 0 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].

6
International Journal of Statistical Mechanics From the stable fixed point of (33) we can calculate the 1RSB free energy functional ( ) as and 1RSB density, 1RSB ( ) = −( ( )/ ), as with expectations intended over G 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 = arg max In the static 1RSB phase we expect Σ( ) = 0 so that from (38) we have ( ) = − . We will see that this is generally true except for the ensemble G (2, ) 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.

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 leafremoval 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).
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 1 ≥ (0) ≥ (1) = P[ = 0]. It turns out that, in the large graph limit, phase 1 of GKS is characterized by the solutions of the system of equations Require: a factor graph = ( , , ) Ensure: a set packing = 0 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  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 phase 1 in the large graphs limit is given by In particular if the smallest solution is * = * the graph is depleted with high probability in phase 1.
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 phase 1 in the large graphs limit is wherẽ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 = 2 ( ) 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 2 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 2 has a unique fixed point. Since those are special cases of the weighted MSP 8 International Journal of Statistical Mechanics problem, we conjecture that the 2 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.

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.

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 : G ( , , ), G ( , , ), G ( , , ), and G ( , , ). 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, = E 0 0 and = E 0 0 . 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) G ( , , ): 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. 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) G ( , , ): 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.

G ( , )
. This is the ensemble with variable node degrees fixed to and Poissonian factor node degree that is ∼ 0 ∼ Poisson( ). The case = 2 corresponds to the maximum matching problem on Erdös-Rényi random graph [8,12].
International Journal of Statistical Mechanics 9 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 | ( * )| = 1 are the critical points ( ) delimiting the RS phase ( < ( )) and are given by For > /( − 1) 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 /( − 1) 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.
In the matching case, that is, = 2 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 ( 0 → 1 | 0 → 1 ) non-zero elements are only: which give = −1 . The average branching factor is = ( − 1) 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.
We continue our analysis of G above the critical value through the 1RSB cavity method as outlined in Section 4.2. Fixed point messages of (33) are distributed as and 2 ( ) has to be determined through (33). Equation  which is stable up to , as already noticed. Above a new stable fixed point, with 2 > 0, 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 [0, +∞] yields the correct free energy; therefore, we have ( ) = − 1RSB .
Unsurprisingly, as they belong to different computational classes, the cases = 2 and ≥ 3 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 = 2.
The ≥ 3 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. G ( , ). The ensemble G ( , ) 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 G , as we will readily show. The MSP problem on G ( , ) with = 2 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)

7.2.
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 G , it can be easily solved numerically for as a function of . For = 2 the threshold value is exactly (2) = . For > 2 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).

G ( , ).
We will now briefly examine our MSP model on the ensemble G ( , ) 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: * = − ( * −1) .
As usual is the parameter characterizing the distribution of messages, ( ) = ( ) + (1 − ) (1 − ). Equation (58) admits one and only one fixed point solution * for each value of and . In fact is continuous, strictly decreasing, and (0) > 0, (1) < 1. The first derivative condition | ( )| = 1 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 G ( , ) and is at variance with the compact area of the RS phase in G . We can compute the MSP relative size through (25) and obtain = (1 − * ) + (1 − * ) ( −1) for < ( ) , (60) 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. G ( , ). The MSP problem on the ensemble G ( , ), 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 G ( , ) 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.

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].