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

Authors: C. Lucibello; F. Ricci-Tersenghi
Journal: International Journal of Statistical Mechanics
 (2014)

---

## 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. Thee-phenomena threshold discovered by Karp and Sipser, marking the onset of core emergence and of replica symmetry breaking, is elegantly generalized to C
s
=
e
/
(d
-
1) for one of the ensembles considered, where d is the size of the sets.


---

## Body

## 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 setF
=
{
1
,
…
,
M
} and a collection of its subsets 𝒮
=
{
S
i
∣
S
i
⊆
F
,
i
∈
V
} labeled by V
=
{
1
,
…
,
N
}; a set packing (SP) is a collection of the subsets S
i 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
(1)
maximize
∑
i
∈
V
‍
n
i
,
(2)
subject  to
∑
i
:
r
∈
S
i
‍
n
i
≤
1
∀
r
∈
F
,
(3)
n
i
=
{
0,1
}
∀
i
∈
V
.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 setsS
i 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 F 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 k-set packing, is that in which all sets S
i have size at most k. 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 k-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 Section2 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 Section3 and give an estimate for the average MSPs size on sparse factor graph ensembles in the thermodynamic limit.
(iii)
In Section4 we establish a criterion for the validity of the RS ansatz and introduce the 1RSB formalism.
(iv)
In Section5 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)
Section6 describes the numerical simulations performed.
(vi)
In Section7 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 k-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 V and to each i
∈
V we associate a variable n
i taking values in {
0,1
} as in (3). F will be our factor nodes set, as its elements acts as hard constrains on the variables n
i through (2). The edge set E is then naturally defined as E
=
{
(
i
,
r
)
∣
i
∈
V
,
r
∈
S
i
⊆
F
}. We call G
=
(
V
,
F
,
E
) the factor graph thus composed and can then rewrite (2) as
(4)
∑
i
∈
∂
r
‍
n
i
≤
1
∀
r
∈
F
.A SP is a configuration{
n
i
} satisfying (4) and its relative size is ρ
(
{
n
i
}
)
=
(
1
/
N
)
∑
i
∈
V
‍
n
i that is simply the fraction of occupied sites.It is now easy to define an appropriate Gibbs measure for the MSPs problem onG through the grand canonical partition function
(5)
Ξ
G
(
μ
)
=
∑
{
n
i
}
‍
∏
i
∈
V
‍
e
μ
n
i
∏
r
∈
F
‍
𝕀
(
∑
i
∈
∂
r
‍
n
i
≤
1
)
.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 𝔾
R
R
(
d
,
c
) defined in Section 7, a d-uniform c-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(6)
ω
G
(
μ
)
=
-
1
μ
N
log
⁡
Ξ
G
(
μ
)
,
and the particle density as
(7)
ρ
G
(
μ
)
=
1
N
〈
∑
i
∈
V
‍
n
i
〉
G
,
μ
=
-
ω
G
(
μ
)
-
μ
∂
μ
ω
G
(
μ
)
.Grand potential and density are related to entropy by the thermodynamic relation(8)
s
G
(
μ
)
=
-
μ
(
ω
G
(
μ
)
+
ρ
G
(
μ
)
)
=
μ
2
∂
μ
ω
G
(
μ
)
.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 ρ
G gives the MSP relative size. Since entropy remains finite in this limit, from (8) we obtain the MSP relative size
(9)
ρ
G
≡
lim
⁡
μ
→
+
∞
ρ
G
(
μ
)
=
lim
⁡
μ
→
+
∞
-
ω
G
(
μ
)
.In this paper we focus on random instances of the MSP problem. As usual in statistical physics we will assume the number of variablesN (the number of subsets S
i) and the number of constrains M to diverge, keeping the ratio N
/
M 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ρ
G is a self-averaging quantity in the thermodynamic limit and we want to compute its asymptotic value
(10)
ρ
=
lim
⁡
N
→
+
∞
𝔼
G
[
ρ
G
]
,
where we denoted with 𝔼
G
[
·
] the expectation over the factor graph ensemble. In the last equation the N 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:
(11)
ω
G
RS
[
ν
^
]
=
1
N
[
∑
r
∈
∂
F
ω
r
[
ν
^
]
+
∑
i
∈
V
(
1
-
|
∂
i
|
)
ω
i
[
ν
^
]
]
,
with the factor and variable contributions given by
(12)
ω
r
[
ν
^
]
=
-
1
μ
log
⁡
[
∑
n
i
∈
∂
r
𝕀
(
∑
i
∈
∂
r
n
i
≤
1
)
∏
j
∈
∂
r
∏
s
∈
∂
j
∖
r
ν
^
s
→
j
(
n
j
)
]
,
ω
i
[
ν
^
]
=
-
1
μ
log
⁡
[
∑
n
i
e
μ
n
i
∏
r
∈
∂
i
ν
^
r
→
i
(
n
i
)
]
.The grand canonical potential is expressed as a function of the factor node to variable node messagesν
^
=
{
ν
^
r
→
i
}. Minimization of ω
G
RS
[
ν
^
] over the messages constrained to be normalized to one yields the fixed point BP equations for the set packing:
(13)
ν
^
r
→
i
(
1
)
=
1
Z
r
→
i
∏
j
∈
∂
r
∖
i
‍
∏
s
∈
∂
j
∖
r
‍
ν
^
s
→
j
(
0
)
,
ν
^
r
→
i
(
0
)
=
1
Z
r
→
i
[
∏
j
∈
∂
r
∖
i
‍
∏
s
∈
∂
j
∖
r
‍
ν
^
s
→
j
(
0
)
=
1
Z
r
→
i
iW
+
e
μ
∑
j
∈
∂
r
∖
i
∏
s
∈
∂
j
∖
r
‍
ν
^
s
→
j
(
1
)
=
1
Z
r
→
i
iW
×
∏
j
′
∈
∂
r
∖
{
i
,
j
}
∏
s
′
∈
∂
j
′
∖
r
‍
ν
^
s
′
→
j
′
(
0
)
]
.The coeeficientsZ
r
→
i are normalization factor. Equations (13) can be simplified introducing the fields {
t
r
→
i
} defined as
(14)
ν
^
r
→
i
(
1
)
ν
^
r
→
i
(
0
)
=
e
-
μ
t
r
→
i
,
yielding
(15)
t
r
→
i
=
1
μ
log
⁡
[
1
+
∑
j
∈
∂
r
∖
i
‍
e
μ
(
1
-
∑
s
∈
∂
j
∖
r
t
s
→
j
)
]
.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
(16)
t
r
→
i
=
max
⁡
{
0
}
∪
{
1
-
∑
s
∈
∂
j
∖
r
‍
t
s
→
j
}
j
∈
∂
r
∖
i
.We note that the messages{
t
r
→
i
} are bounded to take values in the interval [
0,1
] and that if we set the initial value of each t
r
→
i 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 r
→
i if the subtree is connected to the occupied node i. This loss cannot be negative (thus a gain), since we put an additional constrain on the subtree demanding every j neighbour of r to be empty, and cannot be greater than 1 as well, in fact t
r
→
i
=
1 corresponds to the worst case scenario where an otherwise occupied node j
∈
∂
r
∖
i has to be emptied.The Bethe free energy for model (5) on the factor graph G can be expressed as a function of the fixed point messages {
t
r
→
i
} as
(17)
ω
G
RS
(
μ
)
=
-
1
μ
N
[
∑
r
∈
∂
F
log
⁡
(
1
+
∑
j
∈
∂
r
e
μ
(
1
-
∑
s
∈
∂
j
∖
r
t
s
→
j
)
)
=
-
1
μ
N
1
+
∑
i
∈
V
(
1
-
|
∂
i
|
)
log
⁡
(
1
+
e
μ
(
1
-
∑
r
∈
∂
i
t
r
→
i
)
)
(
1
+
∑
j
∈
∂
r
e
μ
(
1
-
∑
s
∈
∂
j
∖
r
‍
t
s
→
j
)
)
]
.We finally arrive to the Bethe estimation of the MSPs relative size, taking the close packing limit of (17) and using ω
G
RS
=
-
ρ
G
RS, which is given by(18)
ρ
G
RS
=
1
N
[
∑
r
∈
F
max
⁡
{
0
}
∪
{
1
-
∑
s
∈
∂
j
∖
r
‍
t
s
→
j
}
j
∈
r
=
1
N
w
+
∑
i
∈
V
(
1
-
|
∂
i
|
)
max
⁡
{
0,1
-
∑
r
∈
i
‍
t
r
→
i
}
[
∑
r
∈
F
max
⁡
{
0
}
∪
{
1
-
∑
s
∈
∂
j
∖
r
‍
t
s
→
j
}
j
∈
r
]
.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
-
|
∂
i
|
)
max
⁡
{
0,1
-
∑
r
∈
i
‍
t
r
→
i
} contributes with 1
-
|
∂
i
| to the sum only if all the incoming t messages are zero. In this case n
i is frozen to 1, that is, the variable i takes part of all the MSPs in G. Obviously all the neighbours of a variable frozen to 1 have to be frozen to 0. To all its |
∂
i
| neighbours r, the frozen to 1 variable i sends a message 1
-
∑
s
∈
∂
i
∖
r
‍
t
s
→
i
=
1, so that we have |
∂
i
| contributions in the first sum of (18) max
⁡
{
0
}
∪
{
1
-
∑
s
∈
∂
j
∖
r
‍
t
s
→
j
}
j
∈
r
=
1 and the total contribution from i correctly sums up to 1. If for a certain i we have a total field τ
i
≡
1
-
∑
r
∈
i
‍
t
r
→
i
<
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 ω
G
RS since it sends a message 1
-
∑
s
∈
i
∖
r
‍
t
s
→
i
≤
0 to each neighbour r; thus, it is not computed in max
⁡
{
0
}
∪
{
1
-
∑
s
∈
∂
j
∖
r
‍
t
s
→
j
}
j
∈
r.The third case is the most interesting. It concerns variablesi which take part to a fraction of the MSPs. We will call them unfrozen variables. The total field on an unfrozen variable i is τ
i
=
0 (thus we have no contribution from the second sum in (18)) and all incoming messages are 0 except for a single t
r
→
i
=
1. To this sole function node r, the node i sends a message 1, so that the contribution of r to the first sum is 1. Actually BP equations impose that r has to have at least another unfrozen neighbour beside i. In other terms the function node r says that whatever MSP we consider, one of my neighbours has to be occupied. The corresponding term max
⁡
{
0
}
∪
{
1
-
∑
s
∈
∂
j
∖
r
‍
t
s
→
j
}
j
∈
r
=
1 in (18) accounts for that.The presence of unfrozen variables is the reason why we cannot express the densityρ
G through the formula
(19)
ρ
G
=
〈
∑
i
∈
V
‍
n
i
〉
.In fact, using the infinite chemical potential formalism, we cannot compute(20)
〈
n
i
〉
=
lim
⁡
μ
→
+
∞
e
μ
(
1
-
∑
r
∈
i
t
r
→
i
)
1
+
e
μ
(
1
-
∑
r
∈
i
t
r
→
i
)
when lim
⁡
μ
→
+
∞
1
-
∑
r
∈
i
‍
t
r
→
i
=
0, and we would have to use the O
(
1
/
μ
) corrections to the fields {
t
r
→
i
}. We bypass the problem using the grand potential ω
G
RS to obtain ρ
G
RS, 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 graphG is sampled from a locally tree-like factor graph ensemble 𝔾
(
N
). We will employ the following notation for the graph ensembles expectations: 𝔼
G
[
·
] for graphs averages; 𝔼
C
0
[
·
] (𝔼
D
0
[
·
]) for expectations over the factor (variable) degree distribution, which we will sometime call root degree distribution; and 𝔼
C
[
·
] (𝔼
D
[
·
]) 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 c and d and the random variables C, C
0, D, and D
0, are related by
(21)
ℙ
[
C
=
k
]
=
(
k
+
1
)
ℙ
[
C
0
=
k
+
1
]
c
,
ℙ
[
D
=
k
]
=
(
k
+
1
)
ℙ
[
D
0
=
k
+
1
]
d
.In Section7 we discuss some specific factor graph ensembles, where C and D are fixed to a deterministic value or Poissonian distributed.With these definitions the distributional equation corresponding to Belief Propagation formula (16) which reads
(22)
T
′
=
d
max
⁡
⁡
{
0
}
∪
{
1
-
∑
s
=
1
D
j
‍
T
s
j
}
j
∈
{
1
,
…
,
C
}
,
where {
D
j
} are i.i.d. random residual variable degrees, C is a random residual factor degree, and {
T
s
j
} are i.i.d. random incoming messages. Since on a given graph each message t
r
→
i takes value only on {
0,1
} we can take the distribution of messages t to be of the form
(23)
P
(
t
)
=
p
δ
(
t
)
+
(
1
-
p
)
δ
(
t
-
1
)
.For this to be a fixed point of (22), the parameter p has to satisfy the self-consistent equation
(24)
p
*
=
𝔼
C
(
1
-
𝔼
D
p
*
D
)
C
.Using (18) and (24), we obtain the replica symmetric approximation to the asymptotic MSP relative size
(25)
ρ
RS
=
d
c
(
1
-
𝔼
C
0
(
1
-
𝔼
D
p
*
D
)
C
0
)
+
𝔼
D
0
(
1
-
D
0
)
p
*
D
0
.It turns out that the RS approximation is exact only when the ratioN
/
M 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 largeN
/
M ratio. In Section 5 we prove that (25) is exact for certain choices of the factor graph ensembles.


## 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:
(11)
ω
G
RS
[
ν
^
]
=
1
N
[
∑
r
∈
∂
F
ω
r
[
ν
^
]
+
∑
i
∈
V
(
1
-
|
∂
i
|
)
ω
i
[
ν
^
]
]
,
with the factor and variable contributions given by
(12)
ω
r
[
ν
^
]
=
-
1
μ
log
⁡
[
∑
n
i
∈
∂
r
𝕀
(
∑
i
∈
∂
r
n
i
≤
1
)
∏
j
∈
∂
r
∏
s
∈
∂
j
∖
r
ν
^
s
→
j
(
n
j
)
]
,
ω
i
[
ν
^
]
=
-
1
μ
log
⁡
[
∑
n
i
e
μ
n
i
∏
r
∈
∂
i
ν
^
r
→
i
(
n
i
)
]
.The grand canonical potential is expressed as a function of the factor node to variable node messagesν
^
=
{
ν
^
r
→
i
}. Minimization of ω
G
RS
[
ν
^
] over the messages constrained to be normalized to one yields the fixed point BP equations for the set packing:
(13)
ν
^
r
→
i
(
1
)
=
1
Z
r
→
i
∏
j
∈
∂
r
∖
i
‍
∏
s
∈
∂
j
∖
r
‍
ν
^
s
→
j
(
0
)
,
ν
^
r
→
i
(
0
)
=
1
Z
r
→
i
[
∏
j
∈
∂
r
∖
i
‍
∏
s
∈
∂
j
∖
r
‍
ν
^
s
→
j
(
0
)
=
1
Z
r
→
i
iW
+
e
μ
∑
j
∈
∂
r
∖
i
∏
s
∈
∂
j
∖
r
‍
ν
^
s
→
j
(
1
)
=
1
Z
r
→
i
iW
×
∏
j
′
∈
∂
r
∖
{
i
,
j
}
∏
s
′
∈
∂
j
′
∖
r
‍
ν
^
s
′
→
j
′
(
0
)
]
.The coeeficientsZ
r
→
i are normalization factor. Equations (13) can be simplified introducing the fields {
t
r
→
i
} defined as
(14)
ν
^
r
→
i
(
1
)
ν
^
r
→
i
(
0
)
=
e
-
μ
t
r
→
i
,
yielding
(15)
t
r
→
i
=
1
μ
log
⁡
[
1
+
∑
j
∈
∂
r
∖
i
‍
e
μ
(
1
-
∑
s
∈
∂
j
∖
r
t
s
→
j
)
]
.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
(16)
t
r
→
i
=
max
⁡
{
0
}
∪
{
1
-
∑
s
∈
∂
j
∖
r
‍
t
s
→
j
}
j
∈
∂
r
∖
i
.We note that the messages{
t
r
→
i
} are bounded to take values in the interval [
0,1
] and that if we set the initial value of each t
r
→
i 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 r
→
i if the subtree is connected to the occupied node i. This loss cannot be negative (thus a gain), since we put an additional constrain on the subtree demanding every j neighbour of r to be empty, and cannot be greater than 1 as well, in fact t
r
→
i
=
1 corresponds to the worst case scenario where an otherwise occupied node j
∈
∂
r
∖
i has to be emptied.The Bethe free energy for model (5) on the factor graph G can be expressed as a function of the fixed point messages {
t
r
→
i
} as
(17)
ω
G
RS
(
μ
)
=
-
1
μ
N
[
∑
r
∈
∂
F
log
⁡
(
1
+
∑
j
∈
∂
r
e
μ
(
1
-
∑
s
∈
∂
j
∖
r
t
s
→
j
)
)
=
-
1
μ
N
1
+
∑
i
∈
V
(
1
-
|
∂
i
|
)
log
⁡
(
1
+
e
μ
(
1
-
∑
r
∈
∂
i
t
r
→
i
)
)
(
1
+
∑
j
∈
∂
r
e
μ
(
1
-
∑
s
∈
∂
j
∖
r
‍
t
s
→
j
)
)
]
.We finally arrive to the Bethe estimation of the MSPs relative size, taking the close packing limit of (17) and using ω
G
RS
=
-
ρ
G
RS, which is given by(18)
ρ
G
RS
=
1
N
[
∑
r
∈
F
max
⁡
{
0
}
∪
{
1
-
∑
s
∈
∂
j
∖
r
‍
t
s
→
j
}
j
∈
r
=
1
N
w
+
∑
i
∈
V
(
1
-
|
∂
i
|
)
max
⁡
{
0,1
-
∑
r
∈
i
‍
t
r
→
i
}
[
∑
r
∈
F
max
⁡
{
0
}
∪
{
1
-
∑
s
∈
∂
j
∖
r
‍
t
s
→
j
}
j
∈
r
]
.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
-
|
∂
i
|
)
max
⁡
{
0,1
-
∑
r
∈
i
‍
t
r
→
i
} contributes with 1
-
|
∂
i
| to the sum only if all the incoming t messages are zero. In this case n
i is frozen to 1, that is, the variable i takes part of all the MSPs in G. Obviously all the neighbours of a variable frozen to 1 have to be frozen to 0. To all its |
∂
i
| neighbours r, the frozen to 1 variable i sends a message 1
-
∑
s
∈
∂
i
∖
r
‍
t
s
→
i
=
1, so that we have |
∂
i
| contributions in the first sum of (18) max
⁡
{
0
}
∪
{
1
-
∑
s
∈
∂
j
∖
r
‍
t
s
→
j
}
j
∈
r
=
1 and the total contribution from i correctly sums up to 1. If for a certain i we have a total field τ
i
≡
1
-
∑
r
∈
i
‍
t
r
→
i
<
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 ω
G
RS since it sends a message 1
-
∑
s
∈
i
∖
r
‍
t
s
→
i
≤
0 to each neighbour r; thus, it is not computed in max
⁡
{
0
}
∪
{
1
-
∑
s
∈
∂
j
∖
r
‍
t
s
→
j
}
j
∈
r.The third case is the most interesting. It concerns variablesi which take part to a fraction of the MSPs. We will call them unfrozen variables. The total field on an unfrozen variable i is τ
i
=
0 (thus we have no contribution from the second sum in (18)) and all incoming messages are 0 except for a single t
r
→
i
=
1. To this sole function node r, the node i sends a message 1, so that the contribution of r to the first sum is 1. Actually BP equations impose that r has to have at least another unfrozen neighbour beside i. In other terms the function node r says that whatever MSP we consider, one of my neighbours has to be occupied. The corresponding term max
⁡
{
0
}
∪
{
1
-
∑
s
∈
∂
j
∖
r
‍
t
s
→
j
}
j
∈
r
=
1 in (18) accounts for that.The presence of unfrozen variables is the reason why we cannot express the densityρ
G through the formula
(19)
ρ
G
=
〈
∑
i
∈
V
‍
n
i
〉
.In fact, using the infinite chemical potential formalism, we cannot compute(20)
〈
n
i
〉
=
lim
⁡
μ
→
+
∞
e
μ
(
1
-
∑
r
∈
i
t
r
→
i
)
1
+
e
μ
(
1
-
∑
r
∈
i
t
r
→
i
)
when lim
⁡
μ
→
+
∞
1
-
∑
r
∈
i
‍
t
r
→
i
=
0, and we would have to use the O
(
1
/
μ
) corrections to the fields {
t
r
→
i
}. We bypass the problem using the grand potential ω
G
RS to obtain ρ
G
RS, 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 graphG is sampled from a locally tree-like factor graph ensemble 𝔾
(
N
). We will employ the following notation for the graph ensembles expectations: 𝔼
G
[
·
] for graphs averages; 𝔼
C
0
[
·
] (𝔼
D
0
[
·
]) for expectations over the factor (variable) degree distribution, which we will sometime call root degree distribution; and 𝔼
C
[
·
] (𝔼
D
[
·
]) 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 c and d and the random variables C, C
0, D, and D
0, are related by
(21)
ℙ
[
C
=
k
]
=
(
k
+
1
)
ℙ
[
C
0
=
k
+
1
]
c
,
ℙ
[
D
=
k
]
=
(
k
+
1
)
ℙ
[
D
0
=
k
+
1
]
d
.In Section7 we discuss some specific factor graph ensembles, where C and D are fixed to a deterministic value or Poissonian distributed.With these definitions the distributional equation corresponding to Belief Propagation formula (16) which reads
(22)
T
′
=
d
max
⁡
⁡
{
0
}
∪
{
1
-
∑
s
=
1
D
j
‍
T
s
j
}
j
∈
{
1
,
…
,
C
}
,
where {
D
j
} are i.i.d. random residual variable degrees, C is a random residual factor degree, and {
T
s
j
} are i.i.d. random incoming messages. Since on a given graph each message t
r
→
i takes value only on {
0,1
} we can take the distribution of messages t to be of the form
(23)
P
(
t
)
=
p
δ
(
t
)
+
(
1
-
p
)
δ
(
t
-
1
)
.For this to be a fixed point of (22), the parameter p has to satisfy the self-consistent equation
(24)
p
*
=
𝔼
C
(
1
-
𝔼
D
p
*
D
)
C
.Using (18) and (24), we obtain the replica symmetric approximation to the asymptotic MSP relative size
(25)
ρ
RS
=
d
c
(
1
-
𝔼
C
0
(
1
-
𝔼
D
p
*
D
)
C
0
)
+
𝔼
D
0
(
1
-
D
0
)
p
*
D
0
.It turns out that the RS approximation is exact only when the ratioN
/
M 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 largeN
/
M 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 {
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 p
∈
[
0,1
], from (22) we obtain the dynamical system
(26)
p
′
=
𝔼
C
(
1
-
𝔼
D
p
D
)
C
≡
f
(
p
)
.The stability criterion(27)
|
f
′
(
p
*
)
|
<
1
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 changingt messages induced by a change in a single message t
r
→
i (1
→
0 or 0
→
1). This is given by
(28)
N
ch
=
𝔼
[
∑
(
s
,
j
)
∑
a
0
,
a
1
b
0
,
b
1
𝕀
(
t
s
→
j
=
b
0
⟶
b
1
∣
t
r
→
i
=
a
0
⟶
a
1
)
]
,
where a
0, a
1, b
0, and b
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
(29)
N
ch
=
∑
s
=
0
+
∞
‍
(
C
¯
D
¯
)
s
∑
a
0
,
a
1
b
0
,
b
1
‍
P
(
t
s
=
b
0
⟶
b
1
∣
t
o
=
a
0
⟶
a
1
)
,
where t
s is a message at distance s from the tree root o, and we defined the average residual degrees C
¯
=
𝔼
C
[
C
] and D
¯
=
𝔼
D
[
D
]. The stability condition N
ch
<
+
∞ yields a constraint on the greatest eigenvalue λ
M of the transfer matrix P
(
b
0
→
b
1
∣
a
0
→
a
1
):
(30)
C
¯
D
¯
λ
M
<
1
.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 Section7.1 to the ensemble 𝔾
R
P. 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 graphG from a locally tree-like ensemble 𝔾. We call Q
r
→
i
(
t
r
→
i
) the distribution of messages on the directed edge r
→
i over the states of the system. We still expect the messages t
r
→
i to take values 0 or 1, so that Q
r
→
i can be parametrized as
(31)
Q
r
→
i
(
t
r
→
i
)
=
q
r
→
i
δ
(
t
r
→
i
)
+
(
1
-
q
r
→
i
)
δ
(
t
r
→
i
-
1
)
.The 1RSB Parisi parameterx
∈
[
0,1
] has to be properly rescaled in order to correctly take the limit μ
↑
∞. Therefore we introduce the new 1RSB parameter y
=
μ
x which stays finite in the close packing limit and takes a value in [
0
,
+
∞
). The reweighting factor e
-
y
ω
iter
r
→
i is defined as
(32)
e
-
y
ω
iter
r
→
i
=
Z
r
→
i
∏
j
∈
∂
r
∖
i
∏
s
∈
∂
j
∖
r
Z
s
→
j
.Last equation combined with (13) and (14) gives ω
iter
r
→
i
=
-
t
r
→
i. Averaging over the whole ensemble we can then write the zero temperature 1RSB message passing rules (also called Survey Propagation equations):
(33)
q
′
=
d
∏
j
=
1
C
(
1
-
∏
s
=
1
D
j
q
s
j
)
e
y
+
(
1
-
e
y
)
∏
j
=
1
C
(
1
-
∏
s
=
1
D
j
q
s
j
)
.In preceding equation{
q
s
j
} are i.i.d.r.v. on [
0,1
] and, as usual, C is the random variable residual degree and {
D
j
} are random independent factors residual degrees. Fixed points of (33) take the form
(34)
P
(
q
)
=
p
0
δ
(
q
)
+
p
1
δ
(
q
-
1
)
+
p
2
P
2
(
q
)
,
where P
2
(
q
) is a continuous distribution on [
0,1
] and p
2
=
1
-
p
0
-
p
1. Parameters p
0 and p
1 have to satisfy the closed equations
(35)
p
1
=
𝔼
C
(
1
-
𝔼
D
(
1
-
p
0
)
D
)
C
,
p
0
=
1
-
𝔼
C
(
1
-
𝔼
D
p
1
D
)
C
.Solutions of (35) with p
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].From the stable fixed point of (33) we can calculate the 1RSB free energy functional ϕ
(
y
) as
(36)
-
y
ϕ
(
y
)
=
d
c
𝔼
log
⁡
[
(
1
-
e
y
)
∏
j
=
1
C
0
‍
(
1
-
∏
s
=
1
D
j
‍
q
s
j
)
+
e
y
]
+
𝔼
(
1
-
D
0
)
log
⁡
[
(
1
-
e
y
)
(
1
-
∏
r
=
1
D
‍
q
r
)
+
e
y
]
,
and 1RSB density, ρ
1
RSB
(
y
)
=
-
(
∂
y
ϕ
(
y
)
/
∂
y
), as
(37)
ρ
1
RSB
(
y
)
=
d
c
𝔼
e
y
(
1
-
∏
j
=
1
C
0
(
1
-
∏
s
=
1
D
j
q
s
j
)
)
(
1
-
e
y
)
∏
j
=
1
C
0
(
1
-
∏
s
=
1
D
j
q
s
j
)
+
e
y
+
𝔼
(
1
-
D
0
)
e
y
∏
r
=
1
D
0
q
r
(
1
-
e
y
)
(
1
-
∏
r
=
1
D
0
q
r
)
+
e
y
,
with expectations intended over 𝔾 and over fixed point messages {
q
s
} and {
q
s
j
}. Since the free energy functional ϕ
(
y
) and the complexity Σ
(
ρ
) are related by the Legendre transform
(38)
Σ
(
ρ
)
=
-
y
ρ
-
y
ϕ
(
y
)
,
with ∂
Σ
/
∂
ρ
=
-
y, through (36) we can compute the complexity taking the inverse transform. Equilibrium states, that is, MSPs, are selected by
(39)
ρ
1
RSB
=
arg max
ρ
⁡
{
ρ
:
Σ
(
ρ
)
≥
0
}
or equivalently taking the 1RSB parameter y to be
(40)
y
s
=
arg max
y
∈
[
0
,
+
∞
]
⁡
ϕ
(
y
)
.In the static 1RSB phase we expectΣ
(
ρ
s
)
=
0 so that from (38) we have ϕ
(
y
s
)
=
-
ρ
s. We will see that this is generally true except for the ensemble 𝔾
R
P
(
2
,
c
) of Section 7, corresponding to maximum matchings on ordinary graphs, where the equilibrium state have maximal complexity and y
s
=
+
∞. The relation
(41)
ρ
1
RSB
=
-
ϕ
(
y
s
)
is always valid though, since for y
s
↑
∞ complexity stays finite.


## 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 {
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 p
∈
[
0,1
], from (22) we obtain the dynamical system
(26)
p
′
=
𝔼
C
(
1
-
𝔼
D
p
D
)
C
≡
f
(
p
)
.The stability criterion(27)
|
f
′
(
p
*
)
|
<
1
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 changingt messages induced by a change in a single message t
r
→
i (1
→
0 or 0
→
1). This is given by
(28)
N
ch
=
𝔼
[
∑
(
s
,
j
)
∑
a
0
,
a
1
b
0
,
b
1
𝕀
(
t
s
→
j
=
b
0
⟶
b
1
∣
t
r
→
i
=
a
0
⟶
a
1
)
]
,
where a
0, a
1, b
0, and b
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
(29)
N
ch
=
∑
s
=
0
+
∞
‍
(
C
¯
D
¯
)
s
∑
a
0
,
a
1
b
0
,
b
1
‍
P
(
t
s
=
b
0
⟶
b
1
∣
t
o
=
a
0
⟶
a
1
)
,
where t
s is a message at distance s from the tree root o, and we defined the average residual degrees C
¯
=
𝔼
C
[
C
] and D
¯
=
𝔼
D
[
D
]. The stability condition N
ch
<
+
∞ yields a constraint on the greatest eigenvalue λ
M of the transfer matrix P
(
b
0
→
b
1
∣
a
0
→
a
1
):
(30)
C
¯
D
¯
λ
M
<
1
.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 Section7.1 to the ensemble 𝔾
R
P. 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 graphG from a locally tree-like ensemble 𝔾. We call Q
r
→
i
(
t
r
→
i
) the distribution of messages on the directed edge r
→
i over the states of the system. We still expect the messages t
r
→
i to take values 0 or 1, so that Q
r
→
i can be parametrized as
(31)
Q
r
→
i
(
t
r
→
i
)
=
q
r
→
i
δ
(
t
r
→
i
)
+
(
1
-
q
r
→
i
)
δ
(
t
r
→
i
-
1
)
.The 1RSB Parisi parameterx
∈
[
0,1
] has to be properly rescaled in order to correctly take the limit μ
↑
∞. Therefore we introduce the new 1RSB parameter y
=
μ
x which stays finite in the close packing limit and takes a value in [
0
,
+
∞
). The reweighting factor e
-
y
ω
iter
r
→
i is defined as
(32)
e
-
y
ω
iter
r
→
i
=
Z
r
→
i
∏
j
∈
∂
r
∖
i
∏
s
∈
∂
j
∖
r
Z
s
→
j
.Last equation combined with (13) and (14) gives ω
iter
r
→
i
=
-
t
r
→
i. Averaging over the whole ensemble we can then write the zero temperature 1RSB message passing rules (also called Survey Propagation equations):
(33)
q
′
=
d
∏
j
=
1
C
(
1
-
∏
s
=
1
D
j
q
s
j
)
e
y
+
(
1
-
e
y
)
∏
j
=
1
C
(
1
-
∏
s
=
1
D
j
q
s
j
)
.In preceding equation{
q
s
j
} are i.i.d.r.v. on [
0,1
] and, as usual, C is the random variable residual degree and {
D
j
} are random independent factors residual degrees. Fixed points of (33) take the form
(34)
P
(
q
)
=
p
0
δ
(
q
)
+
p
1
δ
(
q
-
1
)
+
p
2
P
2
(
q
)
,
where P
2
(
q
) is a continuous distribution on [
0,1
] and p
2
=
1
-
p
0
-
p
1. Parameters p
0 and p
1 have to satisfy the closed equations
(35)
p
1
=
𝔼
C
(
1
-
𝔼
D
(
1
-
p
0
)
D
)
C
,
p
0
=
1
-
𝔼
C
(
1
-
𝔼
D
p
1
D
)
C
.Solutions of (35) with p
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].From the stable fixed point of (33) we can calculate the 1RSB free energy functional ϕ
(
y
) as
(36)
-
y
ϕ
(
y
)
=
d
c
𝔼
log
⁡
[
(
1
-
e
y
)
∏
j
=
1
C
0
‍
(
1
-
∏
s
=
1
D
j
‍
q
s
j
)
+
e
y
]
+
𝔼
(
1
-
D
0
)
log
⁡
[
(
1
-
e
y
)
(
1
-
∏
r
=
1
D
‍
q
r
)
+
e
y
]
,
and 1RSB density, ρ
1
RSB
(
y
)
=
-
(
∂
y
ϕ
(
y
)
/
∂
y
), as
(37)
ρ
1
RSB
(
y
)
=
d
c
𝔼
e
y
(
1
-
∏
j
=
1
C
0
(
1
-
∏
s
=
1
D
j
q
s
j
)
)
(
1
-
e
y
)
∏
j
=
1
C
0
(
1
-
∏
s
=
1
D
j
q
s
j
)
+
e
y
+
𝔼
(
1
-
D
0
)
e
y
∏
r
=
1
D
0
q
r
(
1
-
e
y
)
(
1
-
∏
r
=
1
D
0
q
r
)
+
e
y
,
with expectations intended over 𝔾 and over fixed point messages {
q
s
} and {
q
s
j
}. Since the free energy functional ϕ
(
y
) and the complexity Σ
(
ρ
) are related by the Legendre transform
(38)
Σ
(
ρ
)
=
-
y
ρ
-
y
ϕ
(
y
)
,
with ∂
Σ
/
∂
ρ
=
-
y, through (36) we can compute the complexity taking the inverse transform. Equilibrium states, that is, MSPs, are selected by
(39)
ρ
1
RSB
=
arg max
ρ
⁡
{
ρ
:
Σ
(
ρ
)
≥
0
}
or equivalently taking the 1RSB parameter y to be
(40)
y
s
=
arg max
y
∈
[
0
,
+
∞
]
⁡
ϕ
(
y
)
.In the static 1RSB phase we expectΣ
(
ρ
s
)
=
0 so that from (38) we have ϕ
(
y
s
)
=
-
ρ
s. We will see that this is generally true except for the ensemble 𝔾
R
P
(
2
,
c
) of Section 7, corresponding to maximum matchings on ordinary graphs, where the equilibrium state have maximal complexity and y
s
=
+
∞. The relation
(41)
ρ
1
RSB
=
-
ϕ
(
y
s
)
is always valid though, since for y
s
↑
∞ 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 ano
(
n
) 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 Figure1 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).Algorithm 1:Algorithm  1 generalized Karp-Sipser (GKS).
Require: a factor graph G
= (V, F, E)
Ensure: a set packing V
′
V
′
=
∅
add to V
′ all isolated variable nodes and remove them
from G
remove from G any isolated factor node
while  V is not empty do
if  G has any pendant then
choose a pendant i uniformly at random
add i to V
′
remove i from G, then remove its factor neighbours
and their variable neighbours
else
pick uniformly at random a variable node and add
it to V
′, remove it from G, then remove its factor
neighbours and their variable neighbours
end if
remove from G any isolated factor node
end while
Return  V
′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 Algorithm1.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 callphase 1 the execution of the algorithm up to the point where the first nonpendant variable is added to V
′. It is trivial to show thatphase 1 is enough to find an MSP on a tree factor graph. The interesting thing though is thatphase 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 inphase 1. We will show how the emergence of a core is directly related to replica symmetry breaking. Let us define as usual
(42)
f
(
x
)
≡
𝔼
C
(
1
-
𝔼
D
x
D
)
C
.The functionf is continuous, nonincreasing, and satisfies the relation 1
≥
f
(
0
)
≥
f
(
1
)
=
ℙ
[
C
=
0
]. It turns out that, in the large graph limit,phase 1 of GKS is characterized by the solutions of the system of equations
(43)
p
=
f
(
r
)
,
r
=
f
(
p
)
.We notice that system (43) is equivalent to 1RSB equation (35) once the substitutions p
→
p
0 and r
→
1
-
p
1 are made.Lemma 1.
The system of (43) admits always a (unique) solution p
*
=
r
*.Proof.
By monotony and continuity the functiong
(
p
)
=
f
(
p
)
-
p has a single zero in the interval [
0,1
].If other solutions are present the relevant one is the one with smallestp, as the following theorems certify.Theorem 2.
Let(
p
,
r
) be the solution with smallest p of (43). Then the density of factor nodes surviving p
h
a
s
e
1 in the large graphs limit is given by
(44)
ψ
1
=
1
+
p
~
-
r
~
+
c
p
𝔼
D
[
p
D
-
r
D
]
with
(45)
r
~
=
𝔼
C
0
(
1
-
𝔼
D
p
D
)
C
0
,
p
~
=
𝔼
C
0
(
1
-
𝔼
D
r
D
)
C
0
.
In particular if the smallest solution is p
*
=
r
* the graph is depleted with high probability in p
h
a
s
e
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 survivesphase 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(
p
,
r
) be the solution with smallest p of (43). Then the density of variables assigned in p
h
a
s
e
1 in the large graphs limit is
(46)
ρ
1
=
d
c
(
1
-
r
~
)
+
𝔼
D
0
(
1
-
D
0
)
p
D
0
,
where p
~ has been defined in previous theorem. In particular if the smallest solution is p
*
=
r
* the replica symmetric cavity method prediction (25) is exact.These two theorems imply that the RS prediction forρ, (25), holds true as long asphase 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 p
≠
p
* 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 equationp
=
f
2
(
p
) instead of the two of (43), since the value of r is uniquely determined by the value of p. Then previous theorems state that the RS results hold as long as f
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 f
2 has a unique fixed point. Since those are special cases of the weighted MSP problem, we conjecture that the f
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., O
(
N
)). In presence of a core the GKS algorithm makes some suboptimal random choices (afterphase 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 duringphase 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 structureF, is equivalent to an independent set problem on an appropriate ordinary graph G. The node set of G will be the variable node set of F and we add an edge to G between each pair of nodes having a common factor node neighbour in F. Therefore each neighborhood of a factor node in F forms a clique in G. We then solve the maximum set problem for G 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:𝔾
P
R
(
N
,
d
,
c
), 𝔾
R
P
(
N
,
d
,
c
), 𝔾
P
P
(
N
,
d
,
c
), and 𝔾
R
R
(
N
,
d
,
c
). Subscript R or P 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, d
=
𝔼
D
0
D
0 and c
=
𝔼
C
0
C
0. They are constituted by factor graphs having N variable nodes and M
=
⌊
N
d
/
c
⌋ 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)
𝔾
R
P
(
N
,
d
,
c
): each element G has N variables and M
=
⌊
N
d
/
c
⌋ factors. Every variable node has fixed degree d. G is obtained linking each variable with d factors chosen uniformly and independently at random. The factor nodes degree distribution obtained is Poissonian of mean c with high probability. This ensemble is a model for the k-set packing and will be the main focus of our attention.
(ii)
𝔾
P
R
(
N
,
d
,
c
): each element G has N variables and M
=
⌊
N
d
/
c
⌋ factor. Every factor node has fixed degree c. G is built linking each factors with c variables chosen uniformly and independently at random. The variable nodes degree distribution obtained is Poissonian of mean d with high probability.
(iii)
𝔾
P
P
(
N
,
d
,
c
): each element G has N variables and M
=
⌊
N
d
/
c
⌋ factors. G is built adding an edge (
i
,
r
) with probability c
/
N independently for each choice of a variable i and a factor r. The factor graph obtained has w.h.p Poissonian variable nodes degree distribution of mean d and Poissonian factor nodes degree distribution of mean c.
(iv)
𝔾
R
R
(
N
,
d
,
c
): it is constituted of all factor graphs of N variable nodes of degree d and M
=
N
d
/
c factor nodes of degree c (N
d has to be multiple of c). 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 k-set packing.We will omit the argumentN when we refer to an ensemble in the N
↑
∞ limit.


### 7.1.𝔾
R
P
(
d
,
c
)

This is the ensemble with variable node degrees fixed tod and Poissonian factor node degree that is C
~
C
0
~
Poisson
(
c
). The case d
=
2 corresponds to the maximum matching problem on Erdös-Rényi random graph [8, 12].Density evolution (22) for 𝔾
R
P reduces to
(47)
T
′
=
d
max
⁡
{
0
}
∪
{
1
-
∑
s
=
1
d
-
1
‍
T
s
j
}
j
∈
{
1
,
…
,
C
}
.Considering distributions of the formP
(
t
)
=
p
δ
(
t
)
+
(
1
-
p
)
δ
(
1
-
t
) fixed points of (47) reads
(48)
p
*
=
e
-
c
p
*
d
-
1
≡
f
(
p
*
)
.The last equation admits only one solution for each value ofd and c, as it is easily seen through a monotony argument considering the left and right hand side of the equation. The values of c as a function of d satisfying |
f
′
(
p
*
)
|
=
1 are the critical points c
s
(
d
) delimiting the RS phase (c
<
c
s
(
d
)) and are given by
(49)
c
s
(
d
)
=
e
d
-
1
.Forc
>
e
/
(
d
-
1
) a core survivesphase 1 of the GKS algorithm as stated by Theorem 2 and shown in Figure 2. We notice that the same threshold value e
/
(
d
-
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.Figure 2
Results from GKS algorithm applied to𝔾
R
P
(
3
,
c
). The set packing size after phase 1 (blue line) and after the algorithm stops (red line) are reported. The vertical line is c
s
=
e
/
2 and is drawn as a visual aid to determine the point where the RS assumptions stop to hold true. Above c
s a nonempty core emerges continuously (green line).In the matching case, that is,d
=
2 equation (49) expresses the notorious e-phenomena discovered by Karp and Sipser, while for higher values of d 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 P
(
b
0
→
b
1
∣
a
0
→
a
1
) non-zero elements are only:
(50)
P
(
0
⟶
1
∣
1
⟶
0
)
=
p
d
-
1
,
P
(
1
⟶
0
∣
0
⟶
1
)
=
p
d
-
1
,
which give λ
M
=
p
d
-
1. The average branching factor is C
¯
D
¯
=
c
(
d
-
1
) so that (30) and (48) yield (49). The analytical value for the relative size of MSPs, that is, the particle density ρ, is
(51)
ρ
(
d
,
c
)
=
d
c
(
1
-
p
*
)
+
(
1
-
d
)
p
*
d
for
c
<
c
s
(
d
)
.In Figure3 we compared ρ from (51) as a function of c for some values d with an exact algorithm applied to finite factor graphs (as explained in Section 6), both above and below c
s. Clearly for c
>
c
s the RS approximation is increasingly more inaccurate.Figure 3
RS cavity method analytical prediction equation (51) for MSP size in 𝔾
R
P
(
d
,
c
) 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 c
>
e
/
(
d
-
1
), where it is no longer exact.We continue our analysis of𝔾
R
P above the critical value c
s through the 1RSB cavity method as outlined in Section 4.2. Fixed point messages of (33) are distributed as
(52)
P
(
q
)
=
p
0
δ
(
q
)
+
p
1
δ
(
q
-
1
)
+
p
2
P
2
(
q
)
,
where
(53)
p
1
=
e
-
c
(
1
-
p
0
)
d
-
1
,
p
0
=
1
-
e
-
c
p
1
d
,
p
2
=
1
-
p
0
-
p
1
,
and P
2
(
q
) has to be determined through (33). Equation (53) admits always an RS solution p
1
=
1
-
p
0
=
p
* which is stable up to c
s, as already noticed. Above c
s a new stable fixed point, with p
2
>
0, continuously arises and we study it numerically with a population dynamics algorithm.The 1RSB free energy, as a function of the Parisi parametery, takes the form
(54)
ϕ
(
y
)
=
-
d
y
c
𝔼
log
⁡
[
(
1
-
e
y
)
∏
j
=
1
C
‍
(
1
-
∏
s
=
1
d
-
1
‍
q
s
j
)
+
e
y
]
+
d
-
1
y
𝔼
log
⁡
[
(
1
-
e
y
)
(
1
-
∏
r
=
1
d
‍
q
r
)
+
e
y
]
.As prescribed by the cavity method, the valuey
s which maximizes ϕ
(
y
) over [
0
,
+
∞
] yields the correct free energy; therefore, we have ϕ
(
y
s
)
=
-
ρ
1
RSB.Unsurprisingly, as they belong to different computational classes, the casesd
=
2 and d
≥
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 y on the whole real positive axis. Correct choice for parameter y is then y
s
=
+
∞, 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 d
=
2.Figure 4
Maximum set packing size for𝔾
R
P
(
d
,
c
) as a function of c for d
=
2. 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 y
≫
1. Karp-Sipser and 1RSB cavity method yield the same and exact result.Thed
≥
3 case analysis does not yield such a definite result. The complexity of states Σ is no more a strictly increasing function of y. It reaches its maximum in y
d, the choice of y 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 y
s where complexity changes sign and takes negative values. Therefore y
s is the correct choice for the Parisi parameter which maximizes ϕ
(
y
). 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 y
∈
[
y
d
,
y
s
] (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 beyondphase 1.Figure 5
Complexity in𝔾
R
P
(
d
,
c
), with d
=
3 and c
=
4.0, as a function of the relative MSP size ρ.Figure 6
MSP size in𝔾
R
P
(
3
,
c
) as a function of c. 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.𝔾
P
R
(
d
,
c
)

The ensemble𝔾
P
R
(
d
,
c
) is constitute factor graphs containing factor nodes of degree fixed to c and variable nodes of Poissonian random degree of mean d. It has statistical properties with respect to the MSP problem quite different from those encountered in 𝔾
R
P, as we will readily show. The MSP problem on 𝔾
P
R
(
d
,
c
) with d
=
2 is equivalent to the well known problem of maximum independent sets on Poissonian graphs [14, 32]. The real parameter p characterizing discrete support distributions of messages has to satisfy the fixed point (24) that reads
(55)
p
*
=
(
1
-
e
d
(
p
*
-
1
)
)
c
-
1
.At fixed value ofc the RS ansatz holds up to the critical value d
s
(
c
), which is implicitly given by first derivative condition (27):
(56)
(
c
-
1
)
p
*
(
c
-
2
)
/
(
c
-
1
)
e
(
p
*
-
1
)
d
s
d
s
=
1
.Although critical condition equation (56) is not as elegant as the one we obtained for the ensemble 𝔾
R
P, it can be easily solved numerically for d
s as a function of c. For c
=
2 the threshold value is exactly d
s
(
2
)
=
e. For c
>
2 instead d
s is an increasing function of the factor degree c. Thanks to (25) we readily compute the MSP size in the RS phase:
(57)
ρ
=
d
c
(
1
-
p
*
c
/
(
c
-
1
)
)
+
(
1
-
p
*
d
)
e
d
(
p
*
-
1
)
for
c
<
c
s
(
d
)
.We can see in Figure7 that the MSP size ρ
(
d
,
c
) is a decreasing function in both arguments as expected. Equation (57) can be taken as the RS estimate for MSP size for d
>
d
s
(
c
). The RS estimate is strictly greater than the average size of SPs given by the GKS algorithm at all values of d
>
d
s
(
c
) (see Figure 7).Figure 7
RS cavity method analytical prediction (57) for MSP size in 𝔾
P
R
(
d
,
c
) (dot-dashed) compared with the set packing size fromphase 1 of the GKS algorithm (dashed lines) and with a complete run of the algorithm (continuous lines).


### 7.3.𝔾
P
P
(
d
,
c
)

We will now briefly examine our MSP model on the ensemble𝔾
P
P
(
d
,
c
) where both factor nodes and variable nodes have Poissonian random degrees of mean c and d, respectively. From (23) and (24) we obtain the fixed point condition for the probability distribution of messages:
(58)
p
*
=
e
-
c
e
d
(
p
*
-
1
)
.As usualp is the parameter characterizing the distribution of messages, P
(
t
)
=
p
δ
(
t
)
+
(
1
-
p
)
δ
(
1
-
p
). Equation (58) admits one and only one fixed point solution p
* for each value of d and c. In fact f is continuous, strictly decreasing, and f
(
0
)
>
0
,
f
(
1
)
<
1. The first derivative condition |
f
′
(
p
)
|
=
1 defines the critical line d
s
(
c
) through
(59)
d
s
(
c
)
=
-
1
p
*
log
⁡
(
p
*
)
.The curved
s
(
c
) separates the RS phase from the RSB phase in the c
-
d parametric space (see Figure 8). The unbounded RS region shares some resemblance with the corresponding (although d-discretized) region of 𝔾
P
R
(
d
,
c
) and is at variance with the compact area of the RS phase in 𝔾
R
P.Figure 8
MSP relative size for the ensemble𝔾
P
P
(
d
,
c
) 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
(60)
ρ
=
d
c
(
1
-
p
*
)
+
(
1
-
d
p
*
)
e
d
(
p
-
1
)
for
d
<
d
s
(
c
)
,
which holds only as an approximation in the RSB phase. We can see from Figure 8 that ρ is decreasing both in c and d as was observed in the other ensembles as well.


### 7.4.𝔾
R
R
(
d
,
c
)

The MSP problem on the ensemble𝔾
R
R
(
d
,
c
), 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 𝔾
R
R
(
d
,
c
) 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.


## 7.1.𝔾
R
P
(
d
,
c
)

This is the ensemble with variable node degrees fixed tod and Poissonian factor node degree that is C
~
C
0
~
Poisson
(
c
). The case d
=
2 corresponds to the maximum matching problem on Erdös-Rényi random graph [8, 12].Density evolution (22) for 𝔾
R
P reduces to
(47)
T
′
=
d
max
⁡
{
0
}
∪
{
1
-
∑
s
=
1
d
-
1
‍
T
s
j
}
j
∈
{
1
,
…
,
C
}
.Considering distributions of the formP
(
t
)
=
p
δ
(
t
)
+
(
1
-
p
)
δ
(
1
-
t
) fixed points of (47) reads
(48)
p
*
=
e
-
c
p
*
d
-
1
≡
f
(
p
*
)
.The last equation admits only one solution for each value ofd and c, as it is easily seen through a monotony argument considering the left and right hand side of the equation. The values of c as a function of d satisfying |
f
′
(
p
*
)
|
=
1 are the critical points c
s
(
d
) delimiting the RS phase (c
<
c
s
(
d
)) and are given by
(49)
c
s
(
d
)
=
e
d
-
1
.Forc
>
e
/
(
d
-
1
) a core survivesphase 1 of the GKS algorithm as stated by Theorem 2 and shown in Figure 2. We notice that the same threshold value e
/
(
d
-
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.Figure 2
Results from GKS algorithm applied to𝔾
R
P
(
3
,
c
). The set packing size after phase 1 (blue line) and after the algorithm stops (red line) are reported. The vertical line is c
s
=
e
/
2 and is drawn as a visual aid to determine the point where the RS assumptions stop to hold true. Above c
s a nonempty core emerges continuously (green line).In the matching case, that is,d
=
2 equation (49) expresses the notorious e-phenomena discovered by Karp and Sipser, while for higher values of d 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 P
(
b
0
→
b
1
∣
a
0
→
a
1
) non-zero elements are only:
(50)
P
(
0
⟶
1
∣
1
⟶
0
)
=
p
d
-
1
,
P
(
1
⟶
0
∣
0
⟶
1
)
=
p
d
-
1
,
which give λ
M
=
p
d
-
1. The average branching factor is C
¯
D
¯
=
c
(
d
-
1
) so that (30) and (48) yield (49). The analytical value for the relative size of MSPs, that is, the particle density ρ, is
(51)
ρ
(
d
,
c
)
=
d
c
(
1
-
p
*
)
+
(
1
-
d
)
p
*
d
for
c
<
c
s
(
d
)
.In Figure3 we compared ρ from (51) as a function of c for some values d with an exact algorithm applied to finite factor graphs (as explained in Section 6), both above and below c
s. Clearly for c
>
c
s the RS approximation is increasingly more inaccurate.Figure 3
RS cavity method analytical prediction equation (51) for MSP size in 𝔾
R
P
(
d
,
c
) 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 c
>
e
/
(
d
-
1
), where it is no longer exact.We continue our analysis of𝔾
R
P above the critical value c
s through the 1RSB cavity method as outlined in Section 4.2. Fixed point messages of (33) are distributed as
(52)
P
(
q
)
=
p
0
δ
(
q
)
+
p
1
δ
(
q
-
1
)
+
p
2
P
2
(
q
)
,
where
(53)
p
1
=
e
-
c
(
1
-
p
0
)
d
-
1
,
p
0
=
1
-
e
-
c
p
1
d
,
p
2
=
1
-
p
0
-
p
1
,
and P
2
(
q
) has to be determined through (33). Equation (53) admits always an RS solution p
1
=
1
-
p
0
=
p
* which is stable up to c
s, as already noticed. Above c
s a new stable fixed point, with p
2
>
0, continuously arises and we study it numerically with a population dynamics algorithm.The 1RSB free energy, as a function of the Parisi parametery, takes the form
(54)
ϕ
(
y
)
=
-
d
y
c
𝔼
log
⁡
[
(
1
-
e
y
)
∏
j
=
1
C
‍
(
1
-
∏
s
=
1
d
-
1
‍
q
s
j
)
+
e
y
]
+
d
-
1
y
𝔼
log
⁡
[
(
1
-
e
y
)
(
1
-
∏
r
=
1
d
‍
q
r
)
+
e
y
]
.As prescribed by the cavity method, the valuey
s which maximizes ϕ
(
y
) over [
0
,
+
∞
] yields the correct free energy; therefore, we have ϕ
(
y
s
)
=
-
ρ
1
RSB.Unsurprisingly, as they belong to different computational classes, the casesd
=
2 and d
≥
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 y on the whole real positive axis. Correct choice for parameter y is then y
s
=
+
∞, 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 d
=
2.Figure 4
Maximum set packing size for𝔾
R
P
(
d
,
c
) as a function of c for d
=
2. 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 y
≫
1. Karp-Sipser and 1RSB cavity method yield the same and exact result.Thed
≥
3 case analysis does not yield such a definite result. The complexity of states Σ is no more a strictly increasing function of y. It reaches its maximum in y
d, the choice of y 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 y
s where complexity changes sign and takes negative values. Therefore y
s is the correct choice for the Parisi parameter which maximizes ϕ
(
y
). 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 y
∈
[
y
d
,
y
s
] (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 beyondphase 1.Figure 5
Complexity in𝔾
R
P
(
d
,
c
), with d
=
3 and c
=
4.0, as a function of the relative MSP size ρ.Figure 6
MSP size in𝔾
R
P
(
3
,
c
) as a function of c. 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.𝔾
P
R
(
d
,
c
)

The ensemble𝔾
P
R
(
d
,
c
) is constitute factor graphs containing factor nodes of degree fixed to c and variable nodes of Poissonian random degree of mean d. It has statistical properties with respect to the MSP problem quite different from those encountered in 𝔾
R
P, as we will readily show. The MSP problem on 𝔾
P
R
(
d
,
c
) with d
=
2 is equivalent to the well known problem of maximum independent sets on Poissonian graphs [14, 32]. The real parameter p characterizing discrete support distributions of messages has to satisfy the fixed point (24) that reads
(55)
p
*
=
(
1
-
e
d
(
p
*
-
1
)
)
c
-
1
.At fixed value ofc the RS ansatz holds up to the critical value d
s
(
c
), which is implicitly given by first derivative condition (27):
(56)
(
c
-
1
)
p
*
(
c
-
2
)
/
(
c
-
1
)
e
(
p
*
-
1
)
d
s
d
s
=
1
.Although critical condition equation (56) is not as elegant as the one we obtained for the ensemble 𝔾
R
P, it can be easily solved numerically for d
s as a function of c. For c
=
2 the threshold value is exactly d
s
(
2
)
=
e. For c
>
2 instead d
s is an increasing function of the factor degree c. Thanks to (25) we readily compute the MSP size in the RS phase:
(57)
ρ
=
d
c
(
1
-
p
*
c
/
(
c
-
1
)
)
+
(
1
-
p
*
d
)
e
d
(
p
*
-
1
)
for
c
<
c
s
(
d
)
.We can see in Figure7 that the MSP size ρ
(
d
,
c
) is a decreasing function in both arguments as expected. Equation (57) can be taken as the RS estimate for MSP size for d
>
d
s
(
c
). The RS estimate is strictly greater than the average size of SPs given by the GKS algorithm at all values of d
>
d
s
(
c
) (see Figure 7).Figure 7
RS cavity method analytical prediction (57) for MSP size in 𝔾
P
R
(
d
,
c
) (dot-dashed) compared with the set packing size fromphase 1 of the GKS algorithm (dashed lines) and with a complete run of the algorithm (continuous lines).


## 7.3.𝔾
P
P
(
d
,
c
)

We will now briefly examine our MSP model on the ensemble𝔾
P
P
(
d
,
c
) where both factor nodes and variable nodes have Poissonian random degrees of mean c and d, respectively. From (23) and (24) we obtain the fixed point condition for the probability distribution of messages:
(58)
p
*
=
e
-
c
e
d
(
p
*
-
1
)
.As usualp is the parameter characterizing the distribution of messages, P
(
t
)
=
p
δ
(
t
)
+
(
1
-
p
)
δ
(
1
-
p
). Equation (58) admits one and only one fixed point solution p
* for each value of d and c. In fact f is continuous, strictly decreasing, and f
(
0
)
>
0
,
f
(
1
)
<
1. The first derivative condition |
f
′
(
p
)
|
=
1 defines the critical line d
s
(
c
) through
(59)
d
s
(
c
)
=
-
1
p
*
log
⁡
(
p
*
)
.The curved
s
(
c
) separates the RS phase from the RSB phase in the c
-
d parametric space (see Figure 8). The unbounded RS region shares some resemblance with the corresponding (although d-discretized) region of 𝔾
P
R
(
d
,
c
) and is at variance with the compact area of the RS phase in 𝔾
R
P.Figure 8
MSP relative size for the ensemble𝔾
P
P
(
d
,
c
) 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
(60)
ρ
=
d
c
(
1
-
p
*
)
+
(
1
-
d
p
*
)
e
d
(
p
-
1
)
for
d
<
d
s
(
c
)
,
which holds only as an approximation in the RSB phase. We can see from Figure 8 that ρ is decreasing both in c and d as was observed in the other ensembles as well.


## 7.4.𝔾
R
R
(
d
,
c
)

The MSP problem on the ensemble𝔾
R
R
(
d
,
c
), 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 𝔾
R
R
(
d
,
c
) 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].



---

*Source: 136829-2014-03-10.xml*