Universality of the SAT-UNSAT (jamming) threshold in non-convex continuous constraint satisfaction problems

I Introduction

The exact description of the glassy phases of high-dimensional sphere systems and the discovery that universal predictions at jamming match finite dimensional observations Charbonneau et al. (2017) has renewed the interest in the statistical physics of random constraint satisfaction problems (CSP) with continuous variables. In CSPs, one seeks assignments of a set of $N$ variables that satisfy a system of constraints. Sphere systems clearly belong to this class: one needs to find the positions of $N$ spheres inside a box with the conditions that the spheres do not overlap. Particularly interesting in a statistical physics perspective is the case where the constraints are taken at random from some ensemble Monasson et al. (1999); Mézard et al. (2002); Krzakala et al. (2007). In that case, quite generically in the limit of large systems one observes a phase transition as the density of constraints is increased passing from a satisfiable (SAT) phase where admissible configurations exist to an unsatisfiable (UNSAT) phase, where the minimum number of unsatisfied constraints is larger than zero. This SAT/UNSAT threshold is clearly analogue to the jamming transition in soft-spheres O’Hern et al. (2003); Liu et al. (2011), that separates the low density region where spheres do not overlap from the high density one where no-overlap configurations do not exist. The jamming transition in spheres has highly universal features, with exponents that appear to be independent of the dimension and of the protocol used to produce the jammed configurations O’Hern et al. (2003); Liu et al. (2011); Wyart (2012); Lerner et al. (2013); Charbonneau et al. (2012, 2015, 2017); Goodrich et al. (2016); Kallus (2016). While some exponents, like e.g. the ones relating the number of contacts or the pressure to the packing fraction in the UNSAT phase, take simple semi-integer values Goodrich et al. (2016), other exponents, e.g. the ones describing the distributions of forces and the interparticle distances have non trivial, presumably non-rational values Charbonneau et al. (2014a).

In order to understand the origin of this universality, it is important to study the SAT/UNSAT transition in different CSP models. A crucial ingredient for jamming is the continuous nature of variables. Jamming is the point where the volume of the set of solutions to the problem continuously shrinks to zero, and in its vicinity scaling laws can emerge. To this aim, in Franz and Parisi (2016) it was suggested to study the random perceptron problem Gardner and Derrida (1988) as a prototype of a CSP with continuous variables, generalized to a region of non-convex optimization. It was found that the nontrivial criticality and universality of the SAT/UNSAT transition point was associated to non-convexity. In the convex regime jamming is reached from a liquid, ergodic phase and it is hypostatic and non-critical. In the non-convex regime jamming is reached from a glassy phase, it is critical and in the same universality class of spheres. This led to the conjecture that it exists a large set of continuous CSP that belong to the same universality class. The perceptron emerges therefore as the simplest continuous CSP where glassy phenomena and jamming can be studied. This paradigm has been fruitfully applied to study the vibrational spectrum of glasses at low temperatures. In Franz et al. (2015) the spectrum of the Hessian matrix of the energy minima in the UNSAT phase of the perceptron has been computed showing that it captures essential features of the vibrations of low temperature glasses. Furthermore, in Altieri et al. (2016) it has been shown how to study systematically the free energy landscape of the model using a Thouless-Anderson-Palmer approach Thouless et al. (1977) to obtain, in particular, the vibrational spectrum in the SAT phase. In Franz and Spigler (2016), the avalanches characterizing the glassy phase around jamming have also been studied. Thanks to these studies, the non-convex perceptron now emerges as the simplest model that captures, at the mean field level, all the most important features of the glass and jamming transitions.

The scope of this paper is to give a detailed account of the space of solutions of the random perceptron model. In particular we carefully study the scaling behavior close to jamming, coming both from the SAT and UNSAT phase. The paper is organized as follows. In Sec. II we give a general formulation of continuous CSPs, we discuss the properties of the SAT-UNSAT transition, and we briefly discuss the case of sphere packings to motivate some denominations that are used throughout the paper. In Sec. III we define the random perceptron model, we introduce the replica method to solve it, and we give the main equations needed for its study. In Sec. IV we discuss the zero-temperature phase diagram of the model, and in Sec. V we completely characterize the SAT/UNSAT transition (jamming) line and its critical properties. Finally, we present concluding remarks and perspectives for future work.

II Continuous constraint satisfaction problems

In an abstract form continuous CSPs can be formulated in the following way: find a $N$-dimensional vector $\vec{X}=\{X_{i}\}_{i=1\cdots N}\in\hbox{\msytw R}^{N}$ that satisfies the set of $M$ constraints

The constraints are specified by some real functions $h_{\mu}(\vec{X}):\hbox{\msytw R}^{N}\rightarrow\hbox{\msytw R}$, which can be either deterministic or random (i.e. they may contain some quenched disorder). One can associate this problem to an optimization one, defining a Hamiltonian function taking positive values if at least one constraint is not satified and zero if all the constraints are satisfied. There is a large choice for such a Hamiltonian; here in analogy with harmonic soft spheres O’Hern et al. (2003) we choose

Other choices of $v(h)$ can be considered, provided $v(h>0)=0$ and $v(h<0)>0$. The analysis of the space of solutions can be performed, following Derrida and Gardner Gardner and Derrida (1988), from the study of the partition function

where the measure ${\cal D}\vec{X}$ may include some additional normalization constraint.

One is typically interested in the limit $N\rightarrow\infty$, with the number of constraints $M$ scaled appropriately to have a non-trivial limit. In this limit, a sharp SAT-UNSAT phase transition emerges Monasson et al. (1999); Mézard et al. (2002); Krzakala et al. (2007), and can be characterized by looking at the zero temperature limit of the partition function. In the SAT phase, the ground state energy is equal to zero, and the partition function reduces to the volume of the space of satisfying assignments (whose logarithm is the entropy of solutions). The corresponding homogeneous measure on the space of the solutions gives identical weights to all configurations that satisfy the constraints. In the UNSAT phase, the ground state energy is non-zero and the partition function is dominated by the ground state configurations.

One is usually interested in the free energy per particle,

where the overline indicates an average over the quenched disorder, if it is present in the constraints. This quantity has a finite limit for $N\rightarrow\infty$ and allows one to extract easily all the thermodynamic information; moreover, in presence of quenched disorder, the free energy per particle is usually self-averaging Mézard et al. (1987), i.e. its fluctuations due to the disorder vanish for $N\rightarrow\infty$. In general, ${\rm f}=e-Ts$, where $e$ is the thermodynamic energy and $s$ the thermodynamic entropy. In particular, in the SAT phase, when $T\rightarrow 0$, $e\rightarrow 0$ and $s$ has a finite limit; as a result $-\beta{\rm f}\rightarrow s$. On the contrary, in the UNSAT phase, when $T\rightarrow 0$, $e$ has a finite limit, and $Ts\rightarrow 0$, and as a result ${\rm f}=e$.

Note that besides the thermodynamic (or “equilibrium”) SAT-UNSAT transition, defined as above from the partition function in Eq. (3), one can study many other similar transitions that happen out of equilibrium. Indeed, in the SAT phase at high enough density of constraints, the equilibrium state of the system is a collection of distinct thermodynamic states Krzakala et al. (2007). One can restrict the study to one of these states (also called the “state following” formalism Franz and Parisi (1995); Krzakala and Zdeborová (2010); Rainone et al. (2015); Rainone and Urbani (2016)) and study the SAT-UNSAT transition of the Boltzmann-Gibbs measure restricted to that particular state. It has been shown for sphere systems that this does not change the critical properties of the transition Rainone and Urbani (2016), hence in the following we restrict our study to the equilibrium setting.

II.2 Distribution of gaps

Besides the free energy and its derived quantities, another interesting observable, in particular in the context of jamming, is the probability distribution of “gaps”. In fact, it has been shown that this quantity encodes important information about the marginal stability of jammed configurations that we are going to discuss below Wyart (2012).

A gap variable is just a constraint function $h_{\mu}(\vec{X})$ – the name “gap” originates from the fact that $h_{\mu}=0$ corresponds to a constraint being on the verge of becoming unsatisfied (or a “contact”), and the value of $h_{\mu}$ is thus the distance (“gap”) to this configuration. The gap probability distribution $\rho(h)$ is defined as

where the brackets denote a thermodynamical average, while the overline denotes an average over quenched disorder, if present. From $\rho(h)$ one can derive

which is the average fraction of unsatisfied constraints, or fraction of contacts. To compute the distribution of gaps, we note that the partition function can be written as

Hence, because $-\beta{\rm f}=\overline{\log Z}/N$, we get

In the SAT phase, there are no unsatisfied constraints, and $\rho(h)=0$ for $h<0$, so that $z=0$, while in the UNSAT phase $\rho(h)$ is non-zero for $h<0$ and $z>0$. A particularly important quantity is the limit value of $z$ when one approaches the SAT-UNSAT transition coming from the UNSAT phase where $z>0$. This limit is usually strictly positive, hence $z$ jumps discontinuously at the transition. If this limiting value is such that the number of unsatisfied constraints is exactly equal to the number of degrees of freedom, i.e. $Mz=N$, then the system is said to be isostatic Liu et al. (2011). We can therefore define an isostaticity index $c=(M/N)z$, which is equal to 1 if the system is isostatic. More generally, the system is said to be hypostatic, isostatic or hyperstatic whenever the total number of violated constraints is less ($c<1$), equal ($c=1$) or higher ($c>1$) than the total number of degrees of freedom.

We will be particularly interested in observables of the form ${\cal O}[\vec{X}]=M^{-1}\sum_{\mu=1}^{M}{\cal O}(h_{\mu})\theta(-h_{\mu})$, that are functions of the negative gaps. We define the thermodynamic and disorder average

Special cases of this class of observables are the thermodynamic energy $e=\overline{\left\langle H[\vec{X}]\right\rangle}/N$, the “pressure” $p$, and the fraction of contacts:

II.3 Some consequences of isostaticity: force distribution and soft modes

We now formulate in our abstract continuous CSP setting some well-known consequences of isostaticity Lerner et al. (2013); Charbonneau et al. (2015), i.e. the condition that the number of unsatisfied constraints is equal to the number of degrees of freedom, and $c=1$. We omit for notational simplicity the $\vec{X}$-dependence of the gaps $h_{\mu}$ and we define, with respect to the Hamiltonian given in Eq. (2):

where the set of “contacts” ${\cal C}=\{\mu:h_{\mu}<0\}$ is the set of unsatisfied constraints, of size $|{\cal C}|=Mz=Nc$, the matrix $\SS$ has dimension $Mz\times N$, and the “contact force” vector $\vec{f}$ has dimension $Mz$ because we consider only the non-zero components $f_{\mu}$. We therefore have $\vec{F}=\SS^{T}\vec{f}$. At zero temperature, we are especially interested in minima of the Hamiltonian, which correspond to vanishing “total forces” $\vec{F}=\SS^{T}\vec{f}=\vec{0}$. Furthermore, one is often interested in the small harmonic vibrations around a minimum, which are described by the Hessian matrix

It is particularly interesting to consider the above structure in the jamming limit, i.e. at the SAT-UNSAT transition. Note that “contacts”, i.e. constraints for which $h_{\mu}<0$ in the UNSAT phase close to the transition, must then become marginally satisfied, i.e. $h_{\mu}=0$, right at the transition. Moreover, the average contact force is proportional to the pressure, $|{\cal C}|^{-1}\sum_{\mu\in{\cal C}}f_{\mu}=-[h]/=p/$, which indeed vanishes at jamming. For contacts, one can thus define scaled forces $f^{s}_{\mu}=f_{\mu}/p$, that remain finite at the jamming transition. These scaled forces still satisfy the condition $\SS^{T}\vec{f}^{s}=0$, where the matrix $\SS$ is calculated at the jamming point. Although by definition both $\SS$ and $\vec{f}$ are fully determined by the configuration $\vec{X}$, one can ask in general how many solutions $\vec{f}$ to $\SS^{T}\vec{f}=\vec{0}$ can be found, for fixed $\SS$. Because this is a linear equation, for hypostatic systems ($c<1$) there are in general no solutions, while for hyperstatic systems ($c>1$) there are in general $N(c-1)$ solutions, and for isostatic systems ($c=1$) there is a single solution. Therefore, for an isostatic system, the vector of scaled forces must necessarily coincide with the unique solution of $\SS^{T}\vec{f}^{s}=\vec{0}$. This is particularly useful for numerical calculations, because the force vector vanishes at jamming and it is therefore difficult to evaluate the scaled forces with high precision, while the matrix $\SS$ remains finite and determining its unique zero mode requires much lower precision Charbonneau et al. (2015). Moreover, while in the SAT phase in principle the contact forces vanish, one can still define effective contact forces by following the procedure outlined in Brito and Wyart (2006). Then, one finds that the scaled contact forces also converge to the zero mode of $\SS$ when approaching jamming from the SAT phase. Finally, one can consider the Hessian matrix at jamming. The second term in Eq. (12) vanishes because $h_{\mu}=0$ for contacts, and therefore ${\cal H}=\SS^{T}\SS$. From this one can deduce that in general ${\cal H}$ has $N(1-c)$ zero modes for hypostatic systems, while it has no zero modes for hyperstatic systems; therefore at jamming, which separates the two situations, one finds a large number of very small eigenvalues of the Hessian (“soft modes”) Lerner et al. (2013).

We have seen in this section that, for general continuous CSPs, isostaticity at jamming implies that the contact forces are fully determined by the matrix $\SS$, that also determines the soft modes of the Hessian. This mathematical structure has many other interesting consequences. We will not provide further details on these aspects. The interested reader can consult Ref. (Charbonneau et al., 2015, Supplementary Information) for a detailed review of the properties of $\SS$ in the context of spheres, Ref. Franz et al. (2015) for a study of the Hessian in the perceptron, and Ref. Wyart (2012); Lerner et al. (2013) for a discussion of soft modes in sphere packings.

II.4 Sphere packing as a constraint satisfaction problem

In order to motivate the interest for the observables discussed above (and their names), we discuss here more explicitly the sphere packing problem and its formulation as a continuous CSP. The $d$-dimensional sphere packing problem is indeed a special case of the general setting discussed above. One considers $n$ points in a $d$ dimensions volume, $\vec{x}_{i}\in V\subset\hbox{\msytw R}^{d}$, $i=1,\cdots,n$, hence the total number of degrees of freedom is $N=dn$. There are $M=n(n-1)/2$ constraints corresponding to all possible distinct particle pairs $\mu=\left\langle i,j\right\rangle$ (e.g. with $i<j$), of the form

In this case, $\sigma_{i}$ can be interpreted as the diameter of particle $i$, and then $h_{\left\langle i,j\right\rangle}$ is precisely the physical gap between particles $i,j$. In particular, if $h_{\left\langle i,j\right\rangle}<0$, then the two particles $i,j$ overlap and thus feel a repulsive interaction. This justifies the name “gap” given to $h_{\mu}$, and the name “fraction of contacts” given to $z$. Note that the fraction of contacts defined in Eq. (6) is normalized to the total number of constraints; in particle systems, it is customary to consider instead the average number of contacts per particle,

Hence, an isostaticity index $c=1$ corresponds to $z_{p}=2d$ contacts per particle. The Hamiltonian $H[\vec{X}]$ in Eq. (2), with the choice $v(h)=h^{2}\theta(-h)/2$, is precisely the Hamiltonian of a system of soft harmonic repulsive spheres O’Hern et al. (2003); Liu et al. (2011); the partition function in Eq. (3) is the thermodynamical canonical partition function of the model, and $\rm f$ the associated free energy per particle. The distribution of gaps is simply related to the “radial distribution function” of the particle system, and the average gap is related to the pressure. Finally, $f_{\mu}$ defined in Eq. (11) is precisely the modulus of the contact force associated to the particle pair $\left\langle ij\right\rangle$, and the matrix $\SS$ encodes the network of contacts in the particle packing Charbonneau et al. (2015). In this case, one is then interested in the limit $n,|V|\rightarrow\infty$ with fixed density $n/|V|$.

Note that in this case, in the SAT phase at $T\rightarrow 0$, the Boltzmann-Gibbs measure becomes a uniform measure over all the configurations satisfying the non-overlapping constraint, which coincides with the equilibrium Boltzmann-Gibbs distribution of a system of hard spheres. In the UNSAT phase (which cannot exist for hard spheres) there are overlaps and, at zero temperature, one has instead a mechanically stable assembly of soft repulsive spheres.

In the sphere packing problem, at least when $\sigma_{ij}=\sigma$ (all particles are identical), there is an additional complication because, even within the SAT phase, the system can form a crystal, a phase in which the particles are arranged over a regular lattice in $\hbox{\msytw R}^{d}$. Crystals are usually denser than disordered arrangements, and therefore the SAT-UNSAT transition usually happens within the crystal phase, although the situation is somewhat uncertain in large dimensions Torquato and Stillinger (2010). In this case, the SAT-UNSAT transition (also called “close packing” point) has very different properties with respect to the same transition in disordered assemblies. Yet, one can restrict to the study of disordered configurations and in that case the SAT-UNSAT transition (also called “jamming transition” or “random close packing” in this case) has robustly universal critical properties Parisi and Zamponi (2010); Liu et al. (2011).

III Definition of the model and basic equations

Following Franz and Parisi (2016), in this paper we want to study the properties of the simplest continuous CSP that exhibits a phenomenology similar to that of spheres (once one restricts to their disordered phase). While spheres do not have quenched disorder (the constraints are deterministic functions), the perceptron model has quenched disorder. The presence of quenched disorder eliminates any possible crystal phase and one can then study the equilibrium properties of the model in the disordered phase. Before getting into the study of the phase diagram, and its critical properties at the SAT-UNSAT transition, in this section we give the basic definition of the model and the main equations needed for its study.

The perceptron model, on which we concentrate in the rest of this paper, is defined by the linear functions

with the normalization $\vec{X}\cdot\vec{X}=N$. The vectors $\vec{\xi}^{\mu}$, that following the neural network literature Gardner and Derrida (1988) we call “patterns”, have components $\xi_{i}^{\mu}$ which are independent Gaussian variables of zero average and unit variance. The partition function reads

where the measure ${\cal D}\vec{X}$ contains the constraint $\vec{X}\cdot\vec{X}=N$, i.e. it is the uniform measure on the $N$-dimensional sphere of radius $\sqrt{N}$. For positive $\sigma$ the model can be interpreted as a classifier of the random patterns $\vec{\xi}^{\mu}$, and defines a convex optimization problem. For negative $\sigma$ the model cannot be interpreted as a classifier. It still is a legitimate CSP, but it is non-convex. It has been observed in Franz and Parisi (2016) that this is the interesting regime to describe jamming and glassy phases. In both cases one considers the thermodynamic limit $N,M\rightarrow\infty$ for values of $\alpha=M/N$ fixed.

The disorder average of the free energy can be computed via the replica method:

The free energy $\rm f$ is obtained through standard manipulations Mézard et al. (1987) from the expression of $Z^{n}$ for integer $n$. At the end of a computation sketched in Appendix A, the free energy can be written as a saddle point over the $n\times n$ “replica overlap matrix”

where $\vec{X}_{a}$ are replicas of the configurations of the system. Notice that the spherical contraint on the $\vec{X}$ implies $Q_{aa}=1$. The explicit expression for the free-energy is

where $s.p.$ denotes the saddle point over the values $Q_{ab}$ and

is the replicated free energy. The saddle point equation for $Q_{ab}$ is thus

Finding the solution of such kind of equations is an extremely difficult task. Furthermore, once the solution is found a sensible analytic continuation down to $n\rightarrow 0$ must be taken, which is an additional complication. The solution to both difficulties is provided by the use of hierarchical matrices. Here we do not review in detail this construction, which can be found in several reviews, e.g. Mézard et al. (1987).

III.2 Hierarchical ansatz for the saddle point solution

The general solution of Eq. (22), in the limit $n\rightarrow 0$, is given in terms of a continuous Parisi hierarchical matrix. In this case the matrix $Q_{ab}$ is parametrized by a function $q(x)$ in the interval $x\in$ whose general form is plotted in Fig. 1. For a hierarchical matrix, the first term of $S(Q)$ reads Mézard and Parisi (1991)

where the dot denotes differentiation with respect to $x$, and

The function $f(x,h)$ verifies the Parisi equation (the dot represents a $x$-differentiation, while the prime a $h$-differentiation):

while for $x\notin[x_{m},x_{M}]$ one has $\dot{q}(x)=0$ and $\dot{f}(x,h)=0$, hence $f(x,h)$ is independent of $x$. The boundary condition for $f(x,h)$ at $x=1$ (or, equivalently, at $x=x_{M}$) is:

It is useful to introduce the inverse function of $q(x)$, called $x(q)$, which is defined for $q\in[q_{m},q_{M}]$ where $q_{m}$ and $q_{M}$ are defined in Fig. 1. Defining $f(q,h)\equiv f(x(q),h)$, Eqs. (27,28) become

where now the dot represents $q$-differentiation. Thus, for a given function $x(q)$ and $q_{m}$ and $q_{M}$ we can solve Eq. (29) to compute $f(q,h)$. Then the replicated free energy computed on this particular function is

where we have introduced $\lambda(q)=\lambda(x(q))$, given by

The next step is to find the variational equations for the function $x(q)$, or equivalently $q(x)$.

III.3 Variational equations

In order to obtain the equations that determine the function $x(q)$ we need to impose that the function $f(q,h)$ that appears in Eq. (30) satisfies Eq. (29). A simple way to impose Eq. (29) is to add a Lagrange multiplier $P(q,h)$ to $s[x(q)]$. The new variational free energy thus becomes

Taking the variational equation with respect to $P(q_{M},h)$ and $P(q,h)$ gives back Eq. (29). Now we can take the variational equations with respect to $f(q,h)$, $f(q_{m},h)$, and $x(q)$ Mézard et al. (1987). The resulting equations are:

Note that $P(q,h)$ is normalised to 1 for all $q$. Whenever $x(q)$ has a continuous part (i.e. $\dot{x}(q)\neq 0$), one can differentiate Eq. (34) w.r.t. $q$; the first and second derivatives lead, respectively, to explicit expressions of $\lambda(q)$ and $x(q)$ as functions of $f(q,h)$ and $P(q,h)$:

We stress once again that Eqs. (35) and (36) only hold whenever $\dot{x}(q)\neq 0$.

III.4 Iterative solution of the saddle point equations

The thermodynamic value of the free-energy at given values of the parameters $(\sigma,\alpha)$ can be obtained from the numerical solution of the variational equations. This can be obtained by iteration according to the following procedure:

Start with a guess for $x(q)$, $q_{m}\leq q\leq q_{M}$;

Use Eqs. (29) and (33) to obtain an estimate of $f(q,h)$ and $P(q,h)$;

Get $q_{m}$ from the ratio of Eqs. (34) and (35), computed in $q=q_{m}$;

Obtain a new guess for $q_{M}$ from Eqs. (35) and (31) computed in $q=q_{M}$;

Use Eq. (31) to obtain $\lambda(q)$ and Eq. (36) to obtain $x(q)$;

The numerical solution of the partial differential equations (29) and (33) can be obtained by discretizing the profile $x(q)$ through a stepwise function with $K$ steps (this is known as the $K$RSB approximation) and then increasing $K$ until convergence. The details of this procedure are described in Appendix B.

III.5 Distribution of gaps and contacts

From the general relation in Eq. (8), and using Eq. (III.3) we get

In particular the fraction of contacts is given by

IV The zero temperature phase diagram

The formulae derived in Sec. III hold for any value of the control parameters: the density of constraints $\alpha$, the parameter $\sigma$ that enters into the constraints, and the inverse temperature $\beta$. As a function of these control parameters, the order parameter function $q(x)$ has different forms, giving rise to several distinct phases and phase transitions. To simplify the study of this phase diagram, here we specialise to the zero temperature limit where a sharp SAT-UNSAT (or, equivalently, jamming) phase transition is found. Note that at finite temperature there is always a finite probability of violating some constraint, and the transition is smoothed out. The zero temperature phase diagram of the perceptron has been discussed for $\sigma\geq 0$ in Gardner and Derrida (1988), and for $\sigma<0$ (but $|\sigma|$ not too large) in Franz et al. (2015). Here, we discuss the complete phase diagram for all values of $\sigma$ and $\alpha$, with the result given in Fig. 3, and characterize all the phases and phase transitions that appear.

The solution of replica equations like the ones derived in Sec. III is usually obtained through a series of steps. One starts by the simplest possible solution, called the “replica symmetric” (RS) solution. It corresponds to a constant $q(x)=q_{M}$, in which case the matrix $Q_{ab}=q_{M}$ for all $a\neq b$. Because $\dot{q}(x)=0$, one also has $f(q,h)=f(q_{M},h)=\log\gamma_{1-q_{M}}\star e^{-\beta v(h)}$ and $P(q,h)=P(q_{M},h)=\gamma_{q_{M}}(h+\sigma)$. Thus, at the RS level we have

When we take the zero temperature limit of this expression, we need to specify whether we are in a SAT or UNSAT phase. We should note that the RS ansatz amounts to assume that the space of solutions form a unique connected component, or equivalently that the free energy has a single minimum Mézard et al. (1987); Krzakala et al. (2007). In the SAT phase the value of $q_{M}$ remains finite for $T\rightarrow 0$: because $q_{M}$ measures the similarity between two solutions, when the volume of the space of solutions is finite, two typical solutions differ and thus $q_{M}<1$. Instead, in the UNSAT phase, under the assumption that there is a unique minimum, all the replicas converge towards the same state when $T\rightarrow 0$ and therefore one has $q\rightarrow 1$ in that limit.

Taking the zero temperature limit with constant $q_{M}<1$, we get

Note that in this case the function $-\beta{\rm f}_{RS}(q_{M})$ in Eq. (40) has a finite limit for $T\rightarrow 0$, which gives the entropy of the system, $s(q_{M})=\lim_{T\rightarrow 0}[-\beta{\rm f}_{RS}(q_{M})]$, i.e. the logarithm of the volume of the space of solutions.

The saddle point equation for $q_{M}$ can be obtained either by taking the derivative of Eq. (40) with respect to $q_{M}$ or by considering explicitly the RS ansatz in Eq. (34). In the last case we obtain

One can solve Eq. (43) numerically, and in the SAT phase it is found, as expected, that $q_{M}<1$. Within this solution, the SAT-UNSAT transition point is reached when $q_{M}\rightarrow 1$. Taking this limit in Eq. (43), using the asymptotic properties of the error function, we get the equation for the critical satisfiability threshold $\alpha_{J}(\sigma)$ within the replica symmetric ansatz:

For $\sigma>0$, our optimization problem is convex, the RS solution is always stable (see Sec. IV.1.3), and Eq. (44) gives the correct result for the SAT-UNSAT transition line Gardner and Derrida (1988), as shown in Fig. 3.

IV.1.2 The UNSAT phase

Within the RS ansatz, in the UNSAT phase, there is a unique energy minimum and the free energy converges to its energy. Correspondingly, $q_{M}\rightarrow 1$. At very low temperature, the system performs harmonic vibrations around that minimum, and in that case one can show that $q_{M}=1-\chi T+O(T^{2})$. We thus take the $T\rightarrow 0$ limit with $1-q_{M}=\chi T\rightarrow 0$ at the same time. Plugging this scaling in $f(q_{M},h)=\log\gamma_{1-q_{M}}\star e^{-\beta v(h)}$, with $v(h)=h^{2}\theta(-h)/2$, we get at leading order in $\beta$ (see Appendix C)

and the saddle point equation for $\chi$ is given by

where $\alpha_{J}(\sigma)$ is given in Eq. (44). Eq. (47) has a solution for $\chi$ only for $\alpha>\alpha_{J}(\sigma)$, which indeed defines the UNSAT phase. In the SAT phase, $q_{M}$ remains less than one for $T\rightarrow 0$, while in the UNSAT phase we have $q_{M}=1-\chi T$. To match the two scalings, when $\alpha\rightarrow\alpha_{J}(\sigma)$ from the UNSAT phase, we should have that $\chi\rightarrow\infty$, which indeed follows from Eq. (47). Furthermore, inserting the saddle point equation (47) in the expression of the energy Eq. (46), we obtain

Plugging the RS ansatz in Eq. (39), and using the asymptotic scaling of $f(q_{M},h)$ in Eq. (45), one gets the fraction of contacts in the UNSAT phase (at the replica symmetric level):

The isostaticity condition is $c=\alpha z=1$, and the isostaticity index $c$ is reported in Fig. 2 along the jamming line. For $\sigma>0$ the system is hypostatic on the jamming line, while it becomes isostatic at $\sigma=0$ Franz et al. (2015).

IV.1.3 Stability of the RS solution

After the RS solution, and the associated phase diagram, has been discussed, one should investigate the stability of this solution towards replica symmetry breaking (RSB). RSB can be associated to a continuous de Almeida-Thouless (dAT) instability Mézard et al. (1987), or to the discontinuous appearance of a RSB solution, usually called a “Random First Order Transition” (RFOT) Castellani and Cavagna (2005). In this section we concentrate on the first mechanism, which is the relevant one for the transition in the UNSAT region and in the SAT region at moderately negative values of $\sigma$.

The dAT continuous instability of the RS solution can be discussed by computing the eigenvalues of the Hessian matrix, defined as Mézard et al. (1987)

on the RS saddle point solution $q_{a\neq b}=q_{M}$, where $q_{M}$ satisfies the RS saddle point equation. The continuous breaking of replica symmetry is associated to the vanishing of an eigenvalue of the Hessian matrix.

Here, instead of computing explicitly the eigenvalues of $H_{ab;cd}$, we follow an equivalent procedure that makes use of the fullRSB equations derived in Sec. III. Replica symmetry breaking means that the function $q(x)$ is not a constant; continuous RSB means that $q(x)$ becomes non-constant in a continuous way, and therefore close to the instability $q(x)$ is very close to a constant. Usually, the deviation from a constant is localized around a particular point $x$. We thus assume that there is a single value of $x$ where $\dot{q}(x)$ is continuously becoming different from zero. Slightly in the unstable phase, around this point $x$, Eq. (35) holds. Upon approaching the instability point from the unstable phase, the function $q(x)$ tends to a constant and Eq. (35) reduces to its expression computed on the RS solution:

In the SAT phase ($T\rightarrow 0$ at finite $q_{M}$), Eq. (51) reduces to

In the UNSAT phase ($T\rightarrow 0$ with $q_{M}=1-\chi T$), it instead becomes

IV.2 The nature of the RSB phase

We have seen in Sec. IV.1.3 that replica symmetry must be spontaneously broken in the region delimited by the dAT instability, reported in the phase diagram of Fig. 3. We now characterize the nature of the RSB transition and of the broken symmetry phase. We follow a recipe based on experience with this kind of transitions, along the following logical steps:

The next step is to evaluate whether $m<1$ or $m>1$. In fact, because the function $q(x)$ is defined for $x\in$, a consistent RSB solution requires $m<1$. If this is not the case, then the dAT line cannot be a transition line, it must be preceded by a discontinuous transition of the RFOT kind. In fact one can show that the case $m=1$ separates the two regimes. When $m=1$, the dAT instability splits into two RFOT-like transition lines: a “dynamical transition” where a 1RSB solution appears discontinuously but the free energy remains analytical, and a Kauzmann (or condensation) transition which corresponds to a true phase transition to a spin glass phase Castellani and Cavagna (2005); Krzakala et al. (2007). We will discuss further this situation in Sec. IV.3.

If $m<1$, the dAT instability corresponds to a true continuous phase transition between a “paramagnetic” RS and a “spin glass” RSB phase. However, the dAT instability can give rise either to a fullRSB phase, with a continuous $q(x)$ for $x\in[x_{m},x_{M}]$, or to a 1RSB phase, with $x_{m}=x_{M}=m$ and $q(x)=q_{m}$ for $0<x<m$ and $q(x)=q_{M}$ for $m<x<1$. To determine which one is the case, the final step is to investigate the value of $\dot{q}(m)$. In fact, a fullRSB solution requires $\dot{q}(m)>0$, because $q(x)$ must be an increasing functionSee Baviera and Virasoro (2015) for an attempt of interpreting decreasing solutions in the replica formalism. of $x$. If this is not the case, then the transition is a continuous transition to a 1RSB solution. To compute $\dot{q}(m)$, we derive Eq. (36), expressed as a function of $x$, with respect to $x$. We get

By evaluating Eq. (54) on the RS solution, with $x=m$, we obtain the desired result.

These steps can be performed for both dAT instabilities, in the SAT and UNSAT phases, leading to a full characterization of the RSB transition.

IV.2.2 The UNSAT phase

in the zero temperature limit in the UNSAT phase the breaking point is

The breaking point thus tends to zero proportionally to $\sqrt{T}$, and therefore in the UNSAT phase the dAT instability always leads to a consistent continuous phase transition (a discontinuous transition is never present in this case). We remark that the scaling of the breaking point along the dAT line, $m\propto\sqrt{T}$, is different from the one observed in the Sherrington-Kirkpatrick model in the zero temperature limit, where instead the breaking point remains finite. The origin of this difference will be clarified in Sec. V.

Moreover, the zero temperature limit of the slope $\dot{q}(m)$ at the breaking point, in the UNSAT phase, is given by

which is of order $\sqrt{T}$ and always positive. This implies that the dAT instability line in the UNSAT phase is a transition from a RS solution to a fullRSB one Franz et al. (2015), see Fig. 3. We will see in Sec. V that the scaling with $T$ of both $m$ and $\dot{q}(m)$ is in perfect agreement with the complete scaling form of $q(x)$ in the UNSAT phase.

IV.3 The 1RSB free-energy in the SAT phase

To compute $\alpha_{\rm dyn}(\sigma)$ and $\alpha_{\rm K}(\sigma)$, we are therefore led to consider the 1RSB entropy for $m\approx 1$, where we can write to the leading order

When $m\rightarrow 1$, $q_{0}$ can be obtained by setting to zero the derivative of Eq. (40), which gives the saddle point equation (43). The dynamical transition point $\alpha_{\rm dyn}(\sigma)$ corresponds to the lowest value of $\alpha$ where the variational equation for $q_{1}$, obtained by setting the derivative of Eq. (61) to zero, admits a solution $q_{1}>q_{0}$. At the Kauzmann transition point $\alpha_{\rm K}(\sigma)$, the breaking point $m$, obtained by setting the derivative with respect to $m$ of Eq. (59) to zero, first becomes smaller than one. Upon further increasing $\alpha$, it can be shown that the 1RSB solution becomes unstable and undergoes a Gardner transition towards a fullRSB phase Gardner (1985). The equation that controls the Gardner instability of the 1RSB solution can also be obtained using the fullRSB equations, similarly to what has been done in Sec. IV.1.3 for the RS solution. In the following we assume that the continuous part of $q(x)$ appears in correspondence to the value $q_{1}$, as it usually happens at a Gardner transition Gardner (1985) (the stability towards a continuous breaking at $q_{0}$ can be discussed along similar lines). If $q_{1}$, $q_{0}$ and $m$ satisfy the saddle point conditions, the instability point $\alpha_{\rm G}(\sigma)$ of the 1RSB transition is given by the following condition:

The dynamical, the Kauzmann and Gardner transition lines are plotted in Fig. 3.

V The SAT-UNSAT transition and its critical properties

Having established the phase diagram of the zero temperature random perceptron model (Fig. 3), we discuss here the critical properties of the jamming (or SAT-UNSAT) transition line. The jamming line at $\sigma>0$ has been discussed already in Gardner and Derrida (1988); Franz and Parisi (2016) and it is known to be non-critical, in the sense that the system is hypostatic (see Fig. 2), which according to the analysis of Wyart (2012) does not lead to marginal stability. In this section we want to describe the jamming transition for $\sigma<0$, which falls in the fullRSB phase (Fig. 3): in this case, the system is isostatic at jamming and a non-trivial critical behavior appears Franz and Parisi (2016). The jamming point can be approached both from the SAT and UNSAT phase. From the SAT (unjammed) phase, upon approaching the transition the volume of the space of solutions shrinks to zero and the self-overlap in a cluster of solutions is asymptotically close to one. From the UNSAT (jammed) phase, the energy goes to zero upon approaching the transition. In both cases, universality naturally emerges with a set of nontrivial critical exponents that characterize the scaling of physical quantities on both sides of the transition. We did not attempt to solve numerically the fullRSB equations (see Charbonneau et al. (2014a) for a numerical study of the equations in the hard sphere case); in fact, the values of the critical exponents can be extracted analytically via a scaling analysis of the equations, that we discuss in the rest of this section.

For later use let us discuss the asymptotic behavior of $f(q,h)$ and $P(q,h)$ for $h\rightarrow\pm\infty$ which are generically valid, in particular in the scaling solutions of our interest. For $h\rightarrow\pm\infty$ the boundary condition for $f(q,h)$ reduces to $f(q_{M},h\rightarrow\infty)=0$ and $f(q_{M},h\rightarrow-\infty)=-\frac{h^{2}}{2}\frac{\beta}{1+\beta(1-q_{M})}$. Using Eq. (29) one readily finds

For $h\rightarrow\infty$, Eq. (33) becomes simply $\dot{P}=P^{\prime\prime}/2$, which has a unique solution compatible with the boundary condition

For $h\rightarrow-\infty$ the equation is slightly more complicated; but still, a Gaussian form of the kind $P(q,h)\sim\sqrt{D(q)}e^{-D(q)h^{2}}$ satisfies Eq. (33) for large negative $h$, where $D(q)$ verifies

V.2 Approaching jamming from the SAT phase

In the SAT phase we can take the limit $T=1/\beta\rightarrow 0$ by simply replacing $e^{-\beta v(h)}\rightarrow\theta(h)$. The volume of the space of solutions is finite and the resulting equations are well defined and give a value $q_{M}<1$. This is due to the fact that zero-energy configurations are not isolated: they form clusters where two typical solutions have overlap $q_{M}<1$. Solutions that belong to different clusters have overlap $q<q_{M}$ whose statistical properties are described by the function $x(q)$: it represents the (average) probability that two configurations have an overlap smaller than $q$ Mézard et al. (1987). In the jamming limit, the cluster volumes go to zero and therefore the self-overlap of a cluster, $q_{M}$, goes to one.

Close to jamming the fullRSB equations develop a scaling regime, as can be deduced from a numerical analysis Charbonneau et al. (2014a) (see Sec. III.4 and Appendix B for details on how to solve the fullRSB equations numerically). In the scaling regime, it is convenient to make the following change of variables:

where $\epsilon$ is the linear distance from the jamming line (either in $\alpha$ or $\sigma$). In addition we assume that $q_{M}=1-\epsilon^{\kappa}$ where $\kappa$ is an exponent to be determined by the equations. In the jamming limit $\epsilon\rightarrow 0$, one has $y\in[0,1/\epsilon]\rightarrow[0,\infty)$ and $q\in[q_{m},q_{M}]\rightarrow[q_{m},1]$.

We wish to show that a scaling solution exists in this limit, when $y$ is large and $q\sim 1$. It has the following form:

While the functions $p_{-}$ and $p_{+}$ and the equations that they verify are peculiar of the perceptron model and depend on the parameters $\sigma$ and $\alpha$, on the jamming line the function $p_{0}$ as well as the exponents $a$ and $\kappa$ are universal (i.e. independent of the precise location on the jamming line). In Sec. V.2.2 we show that the universal equations determining $p_{0}$ and the critical exponents also coincide with the ones obtained for the jamming transition of hard spheres in high dimension. Note that for finite $\epsilon$, the scaling solution (68) is cutoff when $q\sim q_{M}=1-\epsilon^{\kappa}$ and $y\sim y_{M}=Y\epsilon^{-1}$. Finally, note that while we will be able to prove the existence of the scaling solution, and compute the values of $a$ and $\kappa$, we will not be able to prove that $\epsilon$ is proportional to the distance from the jamming line: at present, this follows from the numerical solution of the equations Charbonneau et al. (2014a).

V.2.2 Proof of the scaling form

The scaling analysis is carried out along the lines of Charbonneau et al. (2014b, a); Franz and Parisi (2016).

Scaling of $f(q,h)$ – The function $m(q,h)$ introduced in Eq. (68) satisfies the equation

where the differential equation comes from Eq. (29) and the asymptotes from Eq. (64). Let us inspect the value of $\frac{y(q)}{\widehat{\lambda}(q)}$ for $\epsilon\rightarrow 0$ and $q\rightarrow 1$, according to Eq. (68):

This expression has a crossover for $1-q\sim\epsilon^{\kappa}$, when the two terms in the denominator are of the same order. The scaling form is obtained for $1-q\gg\epsilon^{\kappa}$, in which case we have

Plugging this result and the scaling form (68) in Eq. (69), we obtain a scaling equation for the function ${\cal M}(t)$:

This non-linear equation, with boundary conditions at $t\rightarrow\pm\infty$, admits a unique solution for each value of $\kappa$.

Scaling of $P(q,h)$ – The existence of the functions $p_{-}$ and $p_{+}$ and the corresponding scaling variables can be obtained by the analysis of the equation for $P$ at large negative and positive arguments, respectively. From Eq. (65) it follows that for $h\rightarrow\infty$, $P(q,h)$ remains a finite Gaussian, $P(q,h\rightarrow\infty)\sim\gamma_{q}(h+\sigma)$. For finite $h>0$, the Gaussian will be deformed to a finite function $p_{+}(h)$ as it appears in Eq. (68). For $h\rightarrow-\infty$, on the other hand, we have from Eq. (V.1) that for $\beta\rightarrow\infty$ and in the scaling regime of Eq. (71):

If $\kappa<2$ (which, we will see, is the case), this equation admits a scaling solution $D(q)\sim D_{J}\,(1-q)^{-2(\kappa-1)/\kappa}$ where the term $2D(q)^{2}$ is negligible. One concludes therefore that in the scaling regime, $P(q,h\rightarrow-\infty)\sim\sqrt{D(q)}e^{-D(q)h^{2}}$ has the form

which has been proven asymptotically for $h\rightarrow-\infty$ but can be extended to the whole regime where $q\sim 1$ and $|h|\sim(1-q)^{(\kappa-1)/\kappa}$, as in Eq. (68). Note that even if it has a scaling form, the function $p_{-}$ is not uniquely determined by the scaling regime, and remains non-universal, see Charbonneau et al. (2014b).

The “matching” regime with $p_{0}(t)$ must be introduced in Eq. (68) to smoothly match the two regimes for negative and positive $h$. Here is where universality appears. Matching $p_{0}$ and $p_{+}$ requires that

Matching $p_{-}$ and $p_{0}$ requires that

The scaling variable $t=h/\sqrt{1-q}$ appearing in the matching regime is naturally the same as for $f(q,h)$. This is in fact the only choice that leads, once plugged in Eq. (33), to a non-trivial equation for $p_{0}(t)$, namely:

Recall that in Eq. (77), ${\cal M}(t)$ depends on $\kappa$; it turns out that Eq. (77) also admits a unique solution for $p_{0}(t)$ satisfying the correct asymptotic conditions, but only for a given choice of $a=a(\kappa)$. For a given $\kappa$, Eqs. (72) and (77) thus determine ${\cal M}(t)$, $p_{0}(t)$ and $a$.

Determination of the exponent $\kappa$ – The exponent $\kappa$ can be fixed using Eq. (36) which can be equivalently written as

In the scaling regime, Eq. (71) gives the left hand side. In the right hand side, we can note that $m^{\prime\prime}(q,h)$ and $m^{\prime}(q,h)^{2}[1+m^{\prime}(q,h)]$ both vanish outside the scaling regime, because of the asymptotic behavior of $m(q,h)$. Therefore, the right hand side only receives contribution from the regime $h\sim\sqrt{1-q}$. We obtain

Because both ${\cal M}(t)$ and $p_{0}(t)$ depend on $\kappa$, this is an equation for $\kappa$. The numerical solution of the system of Eqs. (72), (77) and (79), gives $\kappa=1.41574\ldots$ while within the numerical precision one finds the relation Charbonneau et al. (2014b)

The importance of these “scaling” relations will be further discussed in Sec. V.5 and Sec. V.6.

V.3 Zero temperature fullRSB solution in the UNSAT phase

We now turn to the analysis of the approach to the jamming transition from the UNSAT phase. Before doing that, however, we need to study the behavior of the fullRSB solution for $T\rightarrow 0$ in the UNSAT phase. In Sec. IV.2 we have shown that in the UNSAT phase, close to the dAT instability line $(\sigma=0,\alpha>2)$ of the RS solution, the solution $q(x)$ has the following properties:

The Edwards-Anderson order parameter is $q_{M}=1-\chi T$

The breaking point at the instability transition line is $m=\hat{m}\sqrt{T}$

The slope of $q(x)$ at the breaking point on the instability line is $\dot{q}(m)\sim\sqrt{T}$

Additionally, here we want to show that, like in the Sherrington-Kirkpatrick model Parisi and Toulouse (1980); Mézard et al. (1987):

$P(q,h)$ is smooth and regular at zero temperature

The function $\beta x(q)$ admits a finite zero temperature limit

For small temperature in the UNSAT phase, the initial condition for $f(q_{M},h)$ is given by (Appendix C):

with ${\cal F}(x)$ defined in Eq. (56). In the fullRSB region, as in the RS case, we define $q_{M}=1-\chi T$ and we introduce

As in Sec. V.2.2, we introduce $m(q,h)=\widehat{\lambda}(q)\widehat{f}^{\prime}(q,h)=\lambda(q)f^{\prime}(q,h)$, which satisfies Eq. (69) with the modification $m(q,h\rightarrow-\infty)=-h\chi\widehat{\lambda}(q)/(1+\chi\widehat{\lambda}(q))$ and initial condition

The equation for $P(q,h)$ is still Eq. (33). Finally, Eq. (35) becomes

and Eq. (36) becomes identical to Eq. (78). The distribution of gaps, given in Eq. (37), becomes

Because the breaking point $m=x(q_{M})\propto\sqrt{T}$, one has $y(q_{M})\sim 1/\sqrt{T}$ and thus $y(q)$ extends up to infinity in the zero temperature limit. The scaling behavior of $y(q)$ for $q\rightarrow 1$ is expected to be Parisi and Toulouse (1980):

We now check that this scaling is consistent. First of all, this implies that $q(x)\sim 1-AT^{2}/x^{2}$ for some constant $A$. For $x=\hat{m}\sqrt{T}$, we get $q_{M}=1-\chi T$ which matches the behavior on the instability line. Also, one has $\dot{q}(x)\propto T^{2}/x^{3}$ and then $\dot{q}(m)\propto\sqrt{T}$, as it is the case along the instability line. The final check can be obtained from Eq. (78). We assume that the scaling behavior of $m(q,h)$ for $q\rightarrow 1$ is

which agrees with the initial condition (83). Plugging this ansatz inside Eq. (78) we get

which confirms the consistency of Eq. (86) and provides an explicit expression of $y_{\chi}$. The only point left to verify is that $P(1,0)$, and more generally $P(1,h)$ are finite. First, we note that for $h\rightarrow\infty$ we have $P(q,h)\rightarrow\gamma_{q}(h+\sigma)$, which suggests that for $h\rightarrow\infty$, $P(1,h)$ is finite and has smooth corrections in $q$. Next, we can observe that Eq. (V.1) becomes for $q\rightarrow 1$:

which therefore indicates that for $h\rightarrow-\infty$, $P(1,h)$ is finite and has corrections proportional to $\sqrt{1-q}$. Plugging $P(q,h)\sim P(1,h)+\sqrt{1-q}\,\delta P(h)$ in Eq. (33), and using Eq. (87), we obtain

V.4 The jamming limit from the UNSAT phase

The jamming limit from the UNSAT phase is obtained by considering the $T=0$ equations of Sec. V.3 in the limit $\chi\rightarrow\infty$, as in the RS case. This is because $q_{M}$ is finite in the SAT phase when $T=0$, while in the UNSAT phase $q_{M}=1-\chi T$ for $T\rightarrow 0$: matching the two regimes for $T\sim 0$ requires the divergence of $\chi$ at the jamming point.

We know that in the UNSAT phase for $q\rightarrow 1$ we have, from Eq. (86) and (82):

In this regime, we must have a different scaling solution in which $y/\widehat{\lambda}\propto 1/(1-q)$: but this is exactly the jamming scaling solution that was already derived in Sec. V.2, which indeed must emerge from the regular zero temperature UNSAT solution upon approaching the jamming point.

To summarize, when $\chi\rightarrow\infty$ and in the region of $q\rightarrow 1$, we expect two different scaling solutions: when $1-q\ll 1-q_{*}$, we have the “regular” UNSAT scaling of Sec. V.3; for $1-q\gg 1-q_{*}$ we have instead the “jamming” scaling solution of Sec. V.2. The matching point $q_{*}$ is determined by the condition that $\chi P(1,0)\sqrt{1-q_{*}}\sim 1$, and the value of $y(q_{*})$ is then simply

Note that in the jamming solution, Eq. (68), we have $P(q_{*},0)\sim(1-q_{*})^{-a/\kappa}$, and we know from Eq. (90) that in the regular solution $P(q,h)$ changes by a very small amount, $\sim\sqrt{1-q}$, when $q\rightarrow 1$. We conclude that $P(1,0)\sim(1-q_{*})^{-a/\kappa}$ and more generally $P(1,h)\sim(1-q_{*})^{(1-\kappa)/\kappa}p_{-}((1-q_{*})^{(1-\kappa)/\kappa}h)$ for negative values of $h$. We obtain therefore

which concludes the analysis of the matching between the two scaling solutions on the UNSAT side of the transition.

V.5 Scaling relations between exponents

The relation $a=1-\kappa/2$, and its consequences that are given in Eq. (80), has been found in Ref. Charbonneau et al. (2014a) within numerical precision by solving the equations for the critical exponents derived in Sec. V.2.2. Eq. (80) is physically very important: in fact, the relation $\gamma=1/(2+\theta)$ has been proven in Wyart (2012) to be a direct consequence of marginal stability (in the perceptron, the same marginal stability argument is discussed in Franz and Parisi (2016)). It would therefore be nice to have a more direct analytical proof of the relation $a=1-\kappa/2$, that does not rely on a numerical calculation.

While in principle it must be possible to obtain such a proof directly from the properties of the equations of Sec. V.2.2, that define all the critical exponents, here we give an independent argument by showing how the matching condition between the two asymptotic solutions allows one to derive analytically this relation. To this aim, we focus on the pressure $p=-[h]$. In Appendix D, we show that $p\propto 1/\chi^{2}$ in the jamming limit $\chi\rightarrow\infty$. Using now $P(1,h)\approx(1-q_{*})^{(1-\kappa)/\kappa}p_{-}((1-q_{*})^{(1-\kappa)/\kappa}h)$ for $h<0$, which was derived in Sec. V.4, together with Eq. (85), we find

Because we know that $[h]\propto 1/\chi^{2}$, we must have

This is compatible with (94) only if $a=1-\kappa/2$, which gives an independent proof of Eq. (80). Finally, recalling that $P(q_{*},h)\sim P(1,h)$, and using Eq. (96) into Eq. (68), we obtain

which holds on the UNSAT side upon approaching jamming.

V.6 Scaling of several interesting observables at the jamming transition

We have now fully characterized the scaling of the basic quantities, $x(q)$, $\lambda(q)$, $f(q,h)$, and $P(q,h)$, in the vicinity of the jamming transition, both on the SAT and on the UNSAT side. On the SAT side, our main results are Eqs. (67) and (68), together with Eqs. (75), (76) and (80). On the UNSAT side, our main results are Eqs. (82), (84), (85) and (97). From these results, the scaling of all the interesting observables can be derived. In this section, we focus on some of the most studied observables in the context of jamming, and we show that all the known results are reproduced by the scaling solution.

We start by considering the energy, the pressure and the number of contacts defined in Eq. (10). From Eqs. (85) and (97), we deduce that in the jamming limit from the UNSAT (jammed) phase,

In particular, the energy and pressure scale as

Note also that from Eq. (84) we have for the isostatic index:

We deduce that at jamming the system is isostatic, and that the excess of contacts in the jammed phase scales like $c-1\propto\chi^{-1}\propto p^{1/2}$. To obtain the complete phenomenology of jamming, one should also prove that the pressure $p\propto\chi^{-2}$ vanishes linearly in the distance from jamming in the $(\alpha,\sigma)$ plane: we did not find a proof of this relation, which at present must be derived from the numerical solution of the equations. Apart from that, the scaling relations $e\propto p^{2}$ and $c-1\propto p^{1/2}$ perfectly agree with numerical observations in jammed sphere packings O’Hern et al. (2003); Liu et al. (2011).

In the SAT phase, at zero temperature, there are by definition no negative gaps, $\rho(h<0)=0$; consequently pressure, energy, and contacts all vanish. However, one can study the limiting values of these quantities when $T\rightarrow 0$. As an example, let us consider the pressure, which is a standard observable in particle systems. In the perceptron, it can be written as the derivative of the free energy with respect to $\sigma$, as it can be checked directly from the definition of the partition function in Eq. (16). Using the replica expression of the free energy in Eq. (30), one then obtains:

Using the equations for $P$ and $f$, it is possible to show that

The scaling form in Eq. (68) is cutoff when $1-q\sim\epsilon^{\kappa}$, and therefore the same scaling holds for $P(1,h)$ and $m(1,h)$ with $1-q\rightarrow\epsilon^{\kappa}$. Then, in the $\epsilon\rightarrow 0$ limit the region $h>0$ does not contribute to the integral because $m(1,h)\rightarrow 0$ while $P(1,h)$ stays finite. One can check that the matching region also gives a subdominant contribution. The leading term is associated to the $h<0$ region, where $m(1,h)\rightarrow-h$ and

We conclude that the reduced pressure $p/T$ diverges proportionally to $\epsilon^{-1}$, and therefore $1-q_{M}\propto(p/T)^{-\kappa}$, a result that has been numerically tested in hard sphere systems Charbonneau et al. (2014a). This provides a physical meaning for the exponent $\kappa$. In a similar way one can study the scaling of the energy and the number of contacts in the SAT phase.

V.6.2 Force and gap distributions

In the UNSAT (jammed) phase, positive gaps correspond to satisfied constraints, while negative gaps can be associated to contact forces according to Eq. (11). The distribution of small positive gaps and small contact forces has been associated with important properties of the packings, including marginal stability Wyart (2012); Lerner et al. (2013). In this section we discuss these distributions. They can be straightforwardly derived from Eq. (85) and Eq. (97) which together imply, for $\chi\gg 1$:

At jamming, when $\chi\rightarrow\infty$, the distribution of gaps concentrates on $h\geq 0$, where it is given by

We conclude by mentioning that one can also study the vibrational spectrum, both in the SAT and UNSAT phases Franz et al. (2015); Altieri et al. (2016), to discuss the presence of soft modes associated to both fullRSB and jamming, and connect the density of states with these critical exponents as it is also observed in sphere packings DeGiuli et al. (2014). Also, it has been shown in Franz and Spigler (2016) that the asymptotic behavior of the function $y(q)$, related to the exponent $\kappa$, controls the scaling of the avalanches at the jamming point, which is therefore different than the scaling in the UNSAT phase.

VI Conclusions

In this work we have formulated the jamming problem as a general constraint satisfaction problem with continuous variables (Sec. II). We have then specialized on the random perceptron, a well known machine learning model, which is a prototype of this class of problems (Sec. III). In the non-convex regime, the model shows a complex zero temperature phase diagram in the plane of the two control parameters $(\alpha,\sigma)$, which has been fully characterized in Sec. IV. In particular, we have shown that for $\sigma<0$ and large enough $|\sigma|$, the phase diagram as a function of $\alpha$ shows the characteristic phenomenology associated with the Random First Order Transition (RFOT) mean field theory of glasses Berthier and Biroli (2011); Cavagna (2009); Parisi and Zamponi (2010); Wolynes and Lubchenko (2012); Kirkpatrick and Thirumalai (2015). Our main result is that the jamming transition, which can be seen as the SAT/UNSAT threshold, is always associated with full replica symmetry breaking in the non-convex regime. In Sec. V, we have thus discussed the scaling behavior of the fullRSB equations around the jamming transition. We have shown that approaching jamming from the SAT phase, one obtains the same critical exponents of the jamming transition of hard spheres in high dimension, thus reproducing the results of Charbonneau et al. (2014a). Furthermore we have extended the study of Charbonneau et al. (2014a) by analyzing the model in the UNSAT phase, where we have obtained the scaling solution of the fullRSB equations, showing that it reproduces the critical behavior of soft spheres, when approaching the transition from the jammed phase O’Hern et al. (2003). We have provided a complete matching between the scaling solutions in the SAT and UNSAT phases, and derived scaling relations between the critical exponents, also showing their physical interpretation as the exponents that control the gap and force distributions at jamming Wyart (2012). Our results are also consistent with the scaling analysis of Goodrich et al. (2016). These results, together with the results on the vibrational spectrum obtained in Franz et al. (2015); Altieri et al. (2016), and the study of avalanches performed in Franz and Spigler (2016), provide a complete study of all the properties of the random perceptron that are relevant for the study of the glass and jamming transition at the mean field level.

This work opens the way to the study of the jamming transition in other ensembles of random constraint satisfaction problems with continuous variables. The outcome of this study is that the non-convex jamming transition lies always in a fullRSB phase and we conjecture that this happens in a large class of CCSP. Although we are not able to prove it we note that this is what happens not only in the model we have analyzed but also in Hard-Spheres in high dimension Charbonneau et al. (2014a). Another natural question that arises is to what extent the scaling behavior that we have found is universal. The fact that the critical exponents at the jamming transition for the random perceptron and hard spheres coincide support the conjecture that in non-convex random CSPs, the jamming point is highly universal. This conjecture has been tested and confirmed in generalized perceptron models with multibody interaction Yoshino (2017). There could be, however, different universality classes: understanding which features of the models determine them is a very important direction for future work. Going beyond the mean field, infinite dimensional level, for example by considering random dilute versions of the perceptron, is another interesting direction for future research.

References

Appendix A Derivation of the replicated free energy

In this Appendix, we give a detailed derivation of the replica equations for the perceptron model. Starting from the partition function in Eq. (16), and neglecting proportionality constants, the replicated partition function can be written as

where $a=1\cdots n$ and $\mu=1\cdots M$; one can check that integrating over $\hat{r}_{a}^{\mu}$ produces delta functions that fix $r_{a}^{\mu}$ as in the original partition function. Next, the Gaussian integral over the quenched disorder $\vec{\xi}^{\mu}$ gives

where we made a change of variables from $\vec{X}_{a}$ to $Q_{ab}$, with Jacobian $\exp(\frac{N}{2}\log\det Q)$. Finally, one can easily show by developing both sides in powers of $Q_{ab}$ that

where to obtain the last line we integrated over $\hat{r}_{a}$ to obtain a $\delta(r_{a}-k_{a})$ and then integrated over $r_{a}$. Introducing $h_{a}=k_{a}-\sigma$ and setting $\alpha=M/N$, we obtain the final result for the replicated partition function

Appendix B Numerical solution of the equations

In this Appendix, we derive the equations corresponding to a discrete $K$RSB ansatz, illustrated in Fig. 5. This is particularly useful for the numerical solution.

For a function $q(x)$ with $K$RSB form, the free energy of Eq. (30) is

and the matrix $M_{ab}$ is also a hierarchical matrix whose components can be written as

which provides a discrete version of Eq. (33). Note that all the $P(m_{i},h)$ are normalized to 1. The saddle point equation for a hierarchical matrix are therefore

To close the equations, we now need to express $q_{i}$ as a function of $q^{-1}_{i}$.

The above equalities are obtained by showing that the derivatives with respect to $x$, as well as the values in $x=0$ or $x=1$ coincide. From these one can derive several useful identities:

From the last Eq. (120) we get $(q^{-1})(0)=-\frac{q(0)}{\lambda(0)^{2}}$. Collecting these two relations we get

This relations allows one to reconstruct $\lambda(x)$ from $(q^{-1})(x)$ and using the first Eq. (120) we can obtain $q(x)$ from $\lambda(x)$. In the discrete case one has $[q]_{i}=\sum_{j=0}^{i-1}m_{j}(q_{j+1}-q_{j})$ and then these equations become

The procedure to solve these equations is therefore (for a fixed grid of $m_{i}$): (i) start from a guess for $q_{i}$, (ii) solve Eqs. (115) and (118) to obtain $f(m_{i},h)$ and $P(m_{i},h)$, (iii) from Eq. (119) compute $q^{-1}_{i}$, and (iv) use Eqs. (123) to compute the new $q_{i}$.

Appendix C Asymptotic behavior in the UNSAT phase

In this Appendix we discuss the asymptotic behavior of $f(q_{M},h)$ in the UNSAT phase where $q_{M}=1-\chi T$. This is given by the zero temperature limit of

For $T\rightarrow 0$ and $h\sim{\cal O}(1)$ we can compute the integral using a saddle point approximation. The saddle point equation is

Appendix D Scaling of the pressure in the UNSAT phase

In this Appendix, we prove that the pressure satisfies the relation $p=-[h]\propto 1/\chi^{2}$, a relation used to prove the scaling relation $a=1-\kappa/2$ in Sec. V.5. To this aim, let us modify the problem by replacing the hard spherical constraint by a Lagrange multiplier $\mu$, and consider the following exact relation

We will prove below that in the zero temperature limit and in the fullRSB phase, $\mu=1/\chi^{2}$, so that the following exact zero temperature relation holds:

Because close to jamming $[h^{2}]\ll[h]$ this reduces to $\lim_{\chi\rightarrow\infty}\chi^{2}\sigma[h]=1$ at jamming, which proves $[h]\propto 1/\chi^{2}$.

We now compute $\mu$ using the replica method. Because $X^{2}$ is not constrained, we have that $Q_{aa}=q_{d}$ and we need to find an additional variational equation for $q_{d}$. Let us write down the free energy as a function of $q_{d}$ and $\mu$. We have

The variational equation for $q(x)$, $f(q,h)$ and $P(q,h)$ do not change except for the initial condition for $f$ which is

and $\lambda(q)$ is given by Eq. (132). Note that the variational equation with respect to $\mu$ fixes the spherical constraint $q_{d}=1$. At this point we can take the variation with respect to $q_{d}$ to get the following equation

It is very easy to show that the last term of the equation above can be rewritten as

where again $m(q_{M},h)=\partial f(q_{M},h)/\partial h$. At this point we can use the saddle point equations (34) and $\eqref{eq:qx2}$ to write

Using that the variational equation over $\mu$ gives $q_{d}=1$, that $q_{M}=1-\chi T$ and that