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

TL;DR

This paper investigates the universality of the SAT-UNSAT (jamming) threshold in non-convex continuous constraint satisfaction problems, revealing critical behavior and universal scaling near the transition, similar to jamming in soft spheres.

Contribution

It demonstrates that non-convex continuous CSPs exhibit a universal critical behavior at the SAT-UNSAT threshold, akin to jamming transitions, and provides detailed phase diagrams and critical exponents.

Findings

Critical behavior at the SAT/UNSAT threshold matches jamming of soft spheres.

The critical exponents are computed for the transition.

A large universality class of such problems is conjectured.

Abstract

Random constraint satisfaction problems (CSP) have been studied extensively using statistical physics techniques. They provide a benchmark to study average case scenarios instead of the worst case one. The interplay between statistical physics of disordered systems and computer science has brought new light into the realm of computational complexity theory, by introducing the notion of clustering of solutions, related to replica symmetry breaking. However, the class of problems in which clustering has been studied often involve discrete degrees of freedom: standard random CSPs are random K-SAT (aka disordered Ising models) or random coloring problems (aka disordered Potts models). In this work we consider instead problems that involve continuous degrees of freedom. The simplest prototype of these problems is the perceptron. Here we discuss in detail the full phase diagram of the model. In the regions of parameter space where the problem is non-convex, leading to multiple disconnected clusters of solutions, the solution is critical at the SAT/UNSAT threshold and lies in the same universality class of the jamming transition of soft spheres. We show how the critical behavior at the satisfiability threshold emerges, and we compute the critical exponents associated to the approach to the transition from both the SAT and UNSAT phase. We conjecture that there is a large universality class of non-convex continuous CSPs whose SAT-UNSAT threshold is described by the same scaling solution.

Figures

Equations

$$h_{\mu }(\vec{X})>0,\mu =1,...,M.$$

$$H[\vec{X}]=\sum _{\mu =1}^{M}v(h_{\mu }(\vec{X}))=\frac{1}{2}\sum _{\mu =1}^M{h}_{\mu }(\vec{X}{)}^2\theta (-h_{\mu }(\vec{X})),v(h)=\frac{1}{2}{h}^2\theta (-h).$$

$$Z=\int \mathcal{D}\vec{X}{e}^{-\beta H[\vec{X}]}=\int \mathcal{D}\vec{X}{e}^{-\frac{\beta }{2}{\sum }_{\mu =1}^M{h}_{\mu }(\vec{X}{)}^2\theta (-h_{\mu }(\vec{X}))},$$

$$\mathrm{f}=-\frac{T}{N}\overline{\log Z},$$

$$\rho (h)=\overline{⟨\hat{\rho }(h)⟩},\hat{\rho }(h)=\frac{1}{M}\sum _{\mu =1}^M\delta (h-h_{\mu }(\vec{X})),$$

$$z\equiv {\int }_{-\infty }^0\mathrm{d}h\rho (h),$$

$$Z=\int \mathcal{D}\vec{X}{e}^{-\beta H[\vec{X}]}=\int \mathcal{D}\vec{X}{e}^{M\int \mathrm{d}h\hat{\rho }(h)[-\beta v(h)]}.$$

$$\frac{\mathrm{d}\mathrm{f}}{\mathrm{d}v(h)}=\frac{1}{N}\frac{\mathrm{d}\overline{\log Z}}{\mathrm{d}[-\beta v(h)]}=\alpha \overline{⟨\hat{\rho }(h)⟩}=\alpha \rho (h).$$

$$[\mathcal{O}(h)]=\overline{⟨\mathcal{O}⟩}\equiv \frac{1}{M}\sum _{\mu =1}^M\overline{⟨\mathcal{O}(h_{\mu })\theta (-h_{\mu })⟩}=\int \mathrm{d}h\rho (h)\mathcal{O}(h)\theta (-h).$$

$$z=[1],p=-[h],e=\alpha [v(h)]=\alpha [h^2]/2.$$

$$\begin{array}{rl}{F}_i& =-\frac{\mathrm{d}H}{\mathrm{d}{X}_i}=\sum _{\mu =1}^M[-h_{\mu }\theta (-h_{\mu })]\frac{\mathrm{d}{h}_{\mu }}{\mathrm{d}{X}_i}=\sum _{\mu \in \mathcal{C}}{\text{SS}}_{\mu i}{f}_{\mu },\\ f_{\mu }& =-h_{\mu }\theta (-h_{\mu }),\\ {\text{SS}}_{\mu i}& =\frac{\mathrm{d}{h}_{\mu }}{\mathrm{d}{X}_i},\end{array}$$

$${\mathcal{H}}_{ij}=\frac{{\mathrm{d}}^{2}H}{\mathrm{d}{X}_i\mathrm{d}{X}_j}=\sum _{\mu \in \mathcal{C}}\frac{\mathrm{d}{h}_{\mu }}{\mathrm{d}{X}_i}\frac{\mathrm{d}{h}_{\mu }}{\mathrm{d}{X}_j}+\sum _{\mu \in \mathcal{C}}{h}_{\mu }\frac{{\mathrm{d}}^2{h}_{\mu }}{\mathrm{d}{X}_i\mathrm{d}{X}_j}.$$

$$h_{\mu }(\vec{X})=|{\vec{x}}_i-{\vec{x}}_j|-{\sigma }_{ij},{\sigma }_{ij}=\frac{{\sigma }_i+{\sigma }_j}{2}.$$

$$z_p=\frac{2Mz}{n}=\frac{2Nc}{n}=2dc.$$

$$h_{\mu }(\vec{X})=\frac{1}{\sqrt{N}}\vec{X}\cdot {\vec{\xi }}^{\mu }-\sigma$$

$$Z=\int \mathcal{D}\vec{X}{e}^{-\beta H[\vec{X}]}=\int \mathcal{D}\vec{X}\prod _{\mu =1}^M\left[\int \mathrm{d}{r}_{\mu }{e}^{-\beta v(r_{\mu }-\sigma )}\delta \left(r_{\mu }-N^{-1/2}\vec{X}\cdot {\vec{\xi }}^{\mu }\right)\right],$$

$$\mathrm{f}=-\frac{T}{N}\overline{\log Z}=-\frac{T}{N}\underset{n\to 0}{\lim }{\partial }_n\overline{Z^n}.$$

$$Q_{ab}=\frac{1}{N}\overline{⟨\overrightarrow{X_a}\cdot \overrightarrow{X_b}⟩},$$

$$\mathrm{f}=-T\underset{n\to 0}{\lim }s.p.{\partial }_{n}S(Q),$$

$$\begin{array}{r}S(Q)=\frac{1}{2}\log \det Q+\alpha \log \left({e^{\frac{1}{2}{\sum }_{ab}^n{Q}_{ab}\frac{{\partial }^2}{\partial h_a\partial h_b}}\prod _{a=1}^n{e}^{-\beta v(h_a)}|}_{h_a=-\sigma }\right)\end{array}$$

$$0={\left[Q^{-1}\right]}_{ab}+\alpha \left({\beta }^2⟨v^{\prime }(h_a)v^{\prime }(h_b)⟩-\beta {\delta }_{ab}⟨v^{\prime \prime }(h_a)⟩\right)$$

$$⟨O(h_a,h_b,\dots )⟩=\frac{{e^{\frac{1}{2}{\sum }_{ab}^n{Q}_{ab}\frac{{\partial }^2}{\partial h_a\partial h_b}}O(h_a,h_b,\dots ){\prod }_{a=1}^n{e}^{-\beta v(h_a)}|}_{h_a=-\sigma }}{{e^{\frac{1}{2}{\sum }_{ab}{Q}_{ab}\frac{{\partial }^2}{\partial h_a\partial h_b}}{\prod }_{a=1}^n{e}^{-\beta v(h_a)}|}_{h_a=-\sigma }}.$$

$$\underset{n\to 0}{\lim }\frac{1}{n}\log \det Q=\log (1-q_M)+\frac{q_m}{\lambda (0)}+{\int }_0^1\mathrm{d}x\frac{\dot{q}(x)}{\lambda (x)},$$

$$\lambda (x)=1-xq(x)-{\int }_x^1\mathrm{d}yq(y).$$

$$\underset{n\to 0}{\lim }\frac{1}{n}\log {e^{\frac{1}{2}{\sum }_{ab}{Q}_{ab}\frac{{\partial }^2}{\partial h_a\partial h_b}}\prod _a{e}^{-\beta v(h_a)}|}_{h_a=-\sigma }={\gamma }_{q_m}\star f(0,h){|}_{h=-\sigma }$$

$${\gamma }_q\star g(h)={\int }_{-\infty }^{\infty }\frac{\mathrm{d}y}{\sqrt{2\pi q}}{e}^{-\frac{y^2}{2q}}g(h-y).$$

$$\dot{f}(x,h)=-\frac{1}{2}\dot{q}(x)\left[f^{\prime \prime }(x,h)+x{f}^{\prime }(x,h{)}^2\right],x_m<x<x_M,$$

$$f(1,h)=\log {\gamma }_{1-q_M}\star e^{-\beta v(h)}.$$

$$\{\begin{array}{l}f(q_M,h)=\log {\gamma }_{1-q_M}\star e^{-\beta v(h)},\\ \dot{f}(q,h)=-\frac{1}{2}\left[f^{\prime \prime }(q,h)+x(q)f^{\prime }(q,h{)}^2\right],q_m<q<q_M.\end{array}$$

$$-\beta \mathrm{f}[x(q)]\equiv \underset{n\to 0}{\lim }{\partial }_{n}S(Q)=\frac{1}{2}\left[\log (1-q_M)+\frac{q_m}{\lambda (q_m)}+{\int }_{q_m}^{q_M}\frac{\mathrm{d}q}{\lambda (q)}\right]+\alpha {\gamma }_{q_m}\star f(q_m,h){|}_{h=-\sigma }$$

Full text

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

II.1 Thermodynamic free energy, and the space of solutions

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

[TABLE]

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

[TABLE]

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

[TABLE]

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,

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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:

[TABLE]

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):

[TABLE]

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

[TABLE]

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]/[1]=p/[1]$, which indeed vanishes at jamming. For contacts, one can thus define scaled forces $f^{s}_{\mu}=[1]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

[TABLE]

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,

[TABLE]

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.

III.1 Definition of the model and free energy

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

[TABLE]

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

[TABLE]

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:

[TABLE]

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”

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

where we have defined

[TABLE]

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[0,1]$ whose general form is plotted in Fig. 1. For a hierarchical matrix, the first term of $S(Q)$ reads Mézard and Parisi (1991)

[TABLE]

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

[TABLE]

Note that $\dot{\lambda}(x)=-x\dot{q}(x)$, and therefore when $q(x)$ is constant also $\lambda(x)$ is constant; furthermore, $\lambda(0)=\lambda(x_{m})=1-\int_{0}^{1}\mathrm{d}xq(x)$, while $\lambda(1)=\lambda(x_{M})=1-q_{M}$. The second term can be written as Duplantier (1981)

[TABLE]

where we have defined the convolution

[TABLE]

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

[TABLE]

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:

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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:

[TABLE]

and

[TABLE]

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)$:

[TABLE]

and

[TABLE]

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}$; 2. 2.

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

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

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

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

Iterate from 2 until convergence.

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

[TABLE]

Note that $\rho(h)$ is correctly normalized so that $\int\mathrm{d}h\rho(h)=1$. Note also that for an observable that is a function of the gaps, we have from Eq. (9):

[TABLE]

In particular the fraction of contacts is given by

[TABLE]

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.

IV.1 The replica symmetric ansatz

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

[TABLE]

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.

IV.1.1 The SAT phase

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

[TABLE]

where

[TABLE]

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

[TABLE]

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:

[TABLE]

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)

[TABLE]

Therefore, the ground state energy is

[TABLE]

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

[TABLE]

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

[TABLE]

For fixed $\sigma$ and $\alpha=\alpha_{\mathrm{J}}(\sigma)+\delta\alpha$, Eq. (47) shows that $\chi^{-1}\sim\delta\alpha$ when $\delta\alpha\rightarrow 0$, and therefore, close to the jamming line, the energy goes as $e\sim\delta\alpha^{2}$. A similar result is obtained at fixed $\alpha$ as a function of $\delta\sigma$: the energy always scales as the square of the distance from the jamming line.

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):

[TABLE]

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)

[TABLE]

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:

[TABLE]

This expression, computed on the value of $q_{M}$ that satisfies the RS saddle point equation, gives a condition on $(\alpha,\sigma)$ that must be satisfied at the dAT instability, hence defining the instability line $\alpha_{\mathrm{dAT}}(\sigma)$ in the phase diagram. Note that computing Eq. (36) on this line gives the “breaking point”, i.e. the value of $x$ at which the instability occurs.

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

[TABLE]

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

[TABLE]

and using Eq. (47) that gives the value of $\chi$, we get that the transition line corresponds to $\sigma=0$, $\forall\alpha>\alpha_{\mathrm{J}}(\sigma=0)=2$. Thus, there are two dAT transition lines where replica symmetry spontaneously breaks, one in the SAT and one in the UNSAT phase; they are reported in Fig. 3. An important conclusion of this study, which is consistent with Ref. Gardner and Derrida (1988), is that in the whole region $\sigma>0$ the RS solution is stable; in particular, in that region the jamming transition happens in the RS phase and the system is hypostatic at jamming.

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 first step is to determine the point at which the RS solution is unstable; this is the dAT line $\alpha_{\mathrm{dAT}}(\sigma)$ determined in Sec. IV.1.3. At the dAT line, the function $q(x)$ is not constant anymore. As mentioned in Sec. IV.1.3, the breaking point $m$ of the RSB solution, i.e. the point $x=m$ where $\dot{q}(x)$ becomes different from zero, can be computed plugging the RS ansatz into Eq. (36). 2. 2.

The next step is to evaluate whether $m<1$ or $m>1$. In fact, because the function $q(x)$ is defined for $x\in[0,1]$, 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. 3. 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 function111See 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

[TABLE]

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.1 The SAT phase

In the SAT phase, the breaking point is

[TABLE]

where $f_{\rm SAT}(q_{M},h)$ is defined in Eq. (41), and $q_{M}$ is the solution of the saddle point equation (43) computed at $\alpha=\alpha_{\mathrm{dAT}}(\sigma)$. A numerical evaluation of $m$ as a function of $\sigma<0$, on the dAT line in the SAT phase, shows that close enough to $\sigma=0$ one has $m<1$, while $m$ increases upon decreasing $\sigma$ towards more negative values. At $\sigma=\sigma_{\mathrm{RFOT}}$ one has $m=1$. Beyond that point, for $\sigma<\sigma_{\mathrm{RFOT}}$, one has $m>1$ and the transition becomes discontinuous, in the RFOT universality class. In order to compute the properties of the phase diagram for $\sigma<\sigma_{\mathrm{RFOT}}$ we need to study the 1RSB solution more carefully, see Sec. IV.3.

Furthermore, the analysis of $\dot{q}(m)$ on the same transition line shows that close enough to $\sigma=0$, for $\sigma>\sigma_{\mathrm{1RSB}}$, one has $\dot{q}(m)>0$ and the dAT instability leads to a fullRSB phase, while for $\sigma\in[\sigma_{\mathrm{RFOT}},\sigma_{\mathrm{1RSB}}]$ one has $\dot{q}(m)<0$ and the instability leads to a continuous transition from a RS phase towards a 1RSB phase. These results are shown in the phase diagram of Fig. 3.

IV.2.2 The UNSAT phase

Defining

[TABLE]

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

[TABLE]

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

[TABLE]

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

In Sec. IV.2 we have shown that for $\sigma<\sigma_{\mathrm{RFOT}}$ the solution for $q(x)$ in the SAT phase becomes of the RFOT type (i.e. a discontinuous 1RSB solution) Castellani and Cavagna (2005). In order to characterize this part of the phase diagram we need to consider explicitly a general 1RSB ansatz. The 1RSB free energy is characterized by three numbers $q_{0}=q_{m},\;q_{1}=q_{M}$ and $m$, the function $q(x)$ being equal to $q_{m}$ for $x\in[0,m]$ and to $q_{M}$ for $x\in[m,1]$ (see Appendix B for a general $K$RSB solution). At zero temperature, in the SAT phase where both $q_{0}<1$ and $q_{1}<1$, the function $-\beta{\rm f}[x(q)]$ has a finite limit, which coincides with the 1RSB entropy, given by

[TABLE]

The variational equations always admit the RS solution $q_{1}=q_{0}$. For any given $\sigma$ this is the unique solution for sufficiently small $\alpha$. For fixed $\sigma<\sigma_{\mathrm{RFOT}}$ and upon increasing $\alpha$, there is a point, called the dynamical transition $\alpha_{\rm dyn}(\sigma)$, where a non-trivial solution with $q_{1}\neq q_{0}$ appears for $m\rightarrow 1$. We do not discuss in details the properties of this transition, for which we refer to Monasson (1995); Krzakala et al. (2007); Castellani and Cavagna (2005). In short, at $\alpha_{\rm dyn}(\sigma)$ the Boltzmann-Gibbs measure on the space of solutions splits into an exponential number of clusters of solutions Krzakala et al. (2007). These clusters can be classified according to their internal entropy that measures their typical size. At the dynamical transition point, the clusters that dominate the Boltzmann-Gibbs measure are the most numerous ones, meaning the ones with higher complexity (or configurational entropy). It is important to note that the free energy remains analytic upon crossing $\alpha_{\rm dyn}(\sigma)$, which is thus not a thermodynamic phase transition Krzakala et al. (2007); Castellani and Cavagna (2005). Increasing $\alpha$, the complexity of the relevant clusters decreases and vanishes at the Kauzmann (or condensation) transition $\alpha=\alpha_{\textrm{K}}(\sigma)$ Monasson (1995); Krzakala et al. (2007); Castellani and Cavagna (2005). At this point, the clusters that dominate the Boltzmann-Gibbs measure are no longer exponentially many, and the system undergoes a thermodynamic phase transition. Note that this two-transition scenario, also called RFOT, is at the basis of the mean field theory of glasses (see Berthier and Biroli (2011); Cavagna (2009); Parisi and Zamponi (2010); Wolynes and Lubchenko (2012); Kirkpatrick and Thirumalai (2015) for reviews).

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

[TABLE]

with $s_{RS}$ given in Eq. (40), and

[TABLE]

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:

[TABLE]

where

[TABLE]

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.

V.1 Asymptotes of $f(q,h)$ and $P(q,h)$

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

[TABLE]

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

[TABLE]

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

[TABLE]

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.

V.2.1 Scaling form of the solution close to jamming

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:

[TABLE]

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:

[TABLE]

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

[TABLE]

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):

[TABLE]

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

[TABLE]

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

[TABLE]

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):

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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:

[TABLE]

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

[TABLE]

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

[TABLE]

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)

[TABLE]

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):

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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):

[TABLE]

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

[TABLE]

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

[TABLE]

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$:

[TABLE]

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

[TABLE]

which satisfies the condition $\int\mathrm{d}h\,\delta P(h)=0$ and shows that there are $\sqrt{1-q}$ corrections only in the region $h<0$.

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):

[TABLE]

If $\chi$ is finite and $q\rightarrow 1$, then in the denominator the factor 1 dominates and $y/\widehat{\lambda}\sim(1-q)^{-1/2}$: in this regime we recover the “regular” zero temperature solution of the UNSAT phase described in Sec. V.3. Conversely, for $\chi\rightarrow\infty$, the coefficient $y_{\chi}\rightarrow\infty$. In fact, from Eq. (88) we obtain $y_{\chi}\propto\chi P(1,0)$ observing that ${\cal M}(t)$ is finite and that $\int_{-\infty}^{0}\mathrm{d}hP(1,h)$ must remain finite because $P(1,h)$ is normalized to 1. In this case, for all $q<1$, at large enough $\chi$ the second term in the denominator of Eq. (91) dominates and one has

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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.

V.6.1 Energy, pressure, number of contacts

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,

[TABLE]

In particular, the energy and pressure scale as

[TABLE]

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

[TABLE]

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:

[TABLE]

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

[TABLE]

The proof is obtained by writing the derivative as $\int\mathrm{d}h(\dot{P}f^{\prime}+P\dot{f}^{\prime})$, using the equations for $P$ and $f$ and integrating by parts Rainone and Urbani (2016). Therefore, in the SAT phase the pressure vanishes proportionally to $T$, a standard result for hard spheres (which correspond to the $T\rightarrow 0$ limit of soft spheres in the SAT phase). However approaching the jamming line the ratio $p/T$ (also called “reduced pressure” or “compressibility factor” in the hard spheres literature) diverges. Inserting the scaling form given by Eq. (68), noting that $\lambda(1)=1-q_{M}=\epsilon^{\kappa}$, one has

[TABLE]

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

[TABLE]

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$:

[TABLE]

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

[TABLE]

where $[1]=\int_{-\infty}^{0}\mathrm{d}h\rho(h)=\int_{-\infty}^{0}\mathrm{d}tp_{-}(t)=1/\alpha$ according to Eq. (100). Isostaticity is expressed by the presence of a delta peak in $\rho(h)$ whose weight is the fraction of contacts. Moreover, this implies that the distribution of small positive (satisfied) gaps, according to Eq. (75), behaves as a power-law $\rho(h\rightarrow 0^{+})\sim h^{-\gamma}$ with $\gamma=(2-\kappa)/\kappa=0.41269\ldots$. Similarly, the normalized probability distribution of scaled contact forces $f^{s}_{\mu}=-[1]h_{\mu}/p$, with $p/[1]$ given by Eq. (99), can be deduced from $\rho(h<0)$ and reads

[TABLE]

Note that $\int_{0}^{\infty}\mathrm{d}fp(f)=\int_{-\infty}^{0}\mathrm{d}fp(f)\,f=1$. This implies, according to Eq. (76), that the distribution of small forces behaves as a power-law, $p(f\rightarrow 0^{+})\sim f^{\theta}$, with $\theta=1/\gamma-2=0.42311\ldots$. These results give a physical interpretation to the critical exponents $\gamma$ and $\theta$, and connect with the results of Wyart (2012); Lerner et al. (2013); Charbonneau et al. (2012, 2015); Kallus (2016). Similar results can be obtained upon approaching the jamming transition from the SAT phase, with the only difference that in that case one has to define the forces by separating the positive gaps into those who become contacts at jamming and those who stay positive Brito and Wyart (2006); Charbonneau et al. (2014b); Charbonneau et al. (2012). As a check of the theoretical predictions we have extended the numerical simulations of Franz and Parisi (2016) to large systems: in Fig. 4 we compare the cumulative distribution $C_{f}=\int_{0}^{f}\mathrm{d}\tilde{f}p(\tilde{f})$ for sizes from $N=50$ to $N=1600$ with the theoretical prediction.

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.

Acknowledgements.

This work was supported by grants from the Simons Foundation (No. 454941, S. F.; No. 454949, G. P.; No. 454955, F. Z.). Ce travail a bénéficié d’une aide Investissements d’Avenir du LabEx PALM (ANR-10-LABX-0039-PALM) (P. Urbani).

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

[TABLE]

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

[TABLE]

We obtain then

[TABLE]

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

[TABLE]

and therefore

[TABLE]

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

[TABLE]

that is reported in Eq. (20).

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

[TABLE]

with $\left\langle q\right\rangle=\int_{0}^{1}\mathrm{d}x\,q(x)$, and $\lambda(x)$ is also a piecewise constant function whose values can be derived from $q(x)$. The function $f(m_{i},h)$ is the discrete version of $f(q,h)$ and it satisfies the discrete version of Eq. (29):

[TABLE]

The saddle point equation is

[TABLE]

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

[TABLE]

with

[TABLE]

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

[TABLE]

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

We start from the exact relations (here $q_{d}=1$ is the diagonal element and $[q](x)=xq(x)-\int_{0}^{x}\mathrm{d}yq(y)$):

[TABLE]

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:

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

Its solutions are given by

[TABLE]

Thus, defining $t=z-z^{*}(h)$ obtain

[TABLE]

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

[TABLE]

from which it follows that

[TABLE]

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:

[TABLE]

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

[TABLE]

where now

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

we get $\mu=1/\chi^{2}$.