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

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.


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 [1] 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 [2][3][4].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 [5,6], 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 [1,[5][6][7][8][9][10][11][12].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 [11], other exponents, e.g. the ones describing the distributions of forces and the interparticle distances have non trivial, presumably non-rational values [13].
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 [14] it was suggested to study the random perceptron problem [15] 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 [16] 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 [17] it has been shown how to study systematically the free energy landscape of the model using a Thouless-Anderson-Palmer approach [18] to obtain, in particular, the vibrational spectrum in the SAT phase.In [19], 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. 2 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. 3 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. 4 we discuss the zero-temperature phase diagram of the model, and in Sec. 5 we completely characterize the SAT/UNSAT transition (jamming) line and its critical properties.Finally, we present concluding remarks and perspectives for future work.

Continuous constraint satisfaction problems 2.1 Thermodynamic free energy, and the space of solutions
In an abstract form continuous CSPs can be formulated in the following way: find a Ndimensional vector The constraints are specified by some real functions h µ ( X ) : N → , 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 [5] 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 [15], from the study of the partition function where the measure X may include some additional normalization constraint.
One is typically interested in the limit N → ∞, with the number of constraints M scaled appropriately to have a non-trivial limit.In this limit, a sharp SAT-UNSAT phase transition emerges [2][3][4], 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 → ∞ 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 [20], i.e. its fluctuations due to the disorder vanish for N → ∞.In general, f = e−Ts, where e is the thermodynamic energy and s the thermodynamic entropy.In particular, in the SAT phase, when T → 0, e → 0 and s has a finite limit; as a result −βf → s.On the contrary, in the UNSAT phase, when T → 0, e has a finite limit, and Ts → 0, and as a result 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 [4].One can restrict the study to one of these states (also called the "state following" formalism [21][22][23][24]) 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 [24], hence in the following we restrict our study to the equilibrium setting.

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 [7].
A gap variable is just a constraint function h µ ( X ) -the name "gap" originates from the fact that h µ = 0 corresponds to a constraint being on the verge of becoming unsatisfied (or a "contact"), and the value of h µ is thus the distance ("gap") to this configuration.The gap probability distribution ρ(h) is defined as where the brackets denote a thermodynamical average, while the overline denotes an average over quenched disorder, if present.From ρ(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 In the SAT phase, there are no unsatisfied constraints, and ρ(h) = 0 for h < 0, so that z = 0, while in the UNSAT phase ρ(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.M z = N , then the system is said to be isostatic [6].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 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 = H[ X ] /N , the "pressure" p, and the fraction of contacts:

Some consequences of isostaticity: force distribution and soft modes
We now formulate in our abstract continuous CSP setting some well-known consequences of isostaticity [8,10], 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 X -dependence of the gaps h µ and we define, with respect to the Hamiltonian given in Eq. ( 2): where the set of "contacts" = {µ : h µ < 0} is the set of unsatisfied constraints, of size , the matrix has dimension M z × N , and the "contact force" vector f has dimension M z because we consider only the non-zero components f µ .We therefore have F = T f .At zero temperature, we are especially interested in minima of the Hamiltonian, which correspond to vanishing "total forces" F = T f = 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 µ < 0 in the UNSAT phase close to the transition, must then become marginally satisfied, i.e. h µ = 0, right at the transition.Moreover, the average contact force is proportional to the pressure, [1], which indeed vanishes at jamming.For contacts, one can thus define scaled forces f s µ = [1] f µ /p, that remain finite at the jamming transition.These scaled forces still satisfy the condition T f s = 0, where the matrix is calculated at the jamming point.Although by definition both and f are fully determined by the configuration X , one can ask in general how many solutions f to T f = 0 can be found, for fixed .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 T f s = 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 remains finite and determining its unique zero mode requires much lower precision [10].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 [25].Then, one finds that the scaled contact forces also converge to the zero mode of when approaching jamming from the SAT phase.Finally, one can consider the Hessian matrix at jamming.The second term in Eq. ( 13) vanishes because h µ = 0 for contacts, and therefore = T .From this one can deduce that in general 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") [8].
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 , 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. [10, Supplementary Information] for a detailed review of the properties of in the context of spheres, Ref. [16] for a study of the Hessian in the perceptron, and Ref. [7,8] for a discussion of soft modes in sphere packings.

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, hence the total number of degrees of freedom is N = d n.There are M = n(n−1)/2 constraints corresponding to all possible distinct particle pairs µ = 〈i, j〉 (e.g. with i < j), of the form In this case, σ i can be interpreted as the diameter of particle i, and then h 〈i, j〉 is precisely the physical gap between particles i, j.In particular, if h 〈i, j〉 < 0, then the two particles i, j overlap and thus feel a repulsive interaction.This justifies the name "gap" given to h µ , 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[ X ] in Eq. ( 2), with the choice v(h) = h 2 θ (−h)/2, is precisely the Hamiltonian of a system of soft harmonic repulsive spheres [5,6]; the partition function in Eq. ( 3) is the thermodynamical canonical partition function of the model, and 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 µ defined in Eq. ( 12) is precisely the modulus of the contact force associated to the particle pair 〈i j〉, and the matrix encodes the network of contacts in the particle packing [10].In this case, one is then interested in the limit n, |V | → ∞ with fixed density n/|V |.
Note that in this case, in the SAT phase at T → 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 σ i j = σ (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 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 [26].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 [6,27].

Definition of the model and basic equations
Following [14], 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.

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 with the normalization X • X = N .The vectors ξ µ , that following the neural network literature [15] we call "patterns", have components ξ µ i which are independent Gaussian variables of zero average and unit variance.The partition function reads where the measure X contains the constraint X • X = N , i.e. it is the uniform measure on the N -dimensional sphere of radius N .For positive σ the model can be interpreted as a classifier of the random patterns ξ µ , and defines a convex optimization problem.For negative σ the model cannot be interpreted as a classifier.It still is a legitimate CSP, but it is nonconvex.It has been observed in [14] that this is the interesting regime to describe jamming and glassy phases.In both cases one considers the thermodynamic limit N , M → ∞ for values of α = M /N fixed.The disorder average of the free energy can be computed via the replica method: The free energy f is obtained through standard manipulations [20] 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 × n "replica overlap matrix" where X a are replicas of the configurations of the system.Notice that the spherical contraint on the 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 , a = b, is thus where we have defined 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 → 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.[20].

Hierarchical ansatz for the saddle point solution
The general solution of Eq. ( 23), in the limit n → 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 ∈ [0, 1] whose general form is plotted in Fig. 1.For a hierarchical matrix, the first term of S(Q) reads [28] where the dot denotes differentiation with respect to x, and Note that λ(x) = −x q(x), and therefore when q(x) is constant also λ The second term can be written as [29] The expected form of the function q(x) in the fullRSB phase.
where we have defined the convolution The function f (x, h) verifies the Parisi equation (the dot represents a x-differentiation, while the prime a h-differentiation): while for x / ∈ [x m , x M ] one has q(x) = 0 and ḟ (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 ∈ [q m , q M ] where q m and q M are defined in Fig. 1.Defining f (q, h) ≡ f (x(q), h), Eqs.(28,29) 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. ( 30) to compute f (q, h).Then the replicated free energy computed on this particular function is where we have introduced λ(q) = λ(x(q)), given by The next step is to find the variational equations for the function x(q), or equivalently q(x).

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. ( 31) satisfies Eq. (30).A simple way to impose Eq. ( 30) 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. ( 30).Now we can take the variational equations with respect to f (q, h), f (q m , h), and x(q) [20].The resulting equations are: and Note that P(q, h) is normalised to 1 for all q.Whenever x(q) has a continuous part (i.e.ẋ(q) = 0), one can differentiate Eq. ( 35) w.r.t.q; the first and second derivatives lead, respectively, to explicit expressions of λ(q) and x(q) as functions of f (q, h) and P(q, h): and We stress once again that Eqs. ( 36) and ( 37) only hold whenever ẋ(q) = 0.

Iterative solution of the saddle point equations
The thermodynamic value of the free-energy at given values of the parameters (σ, α) can be obtained from the numerical solution of the variational equations.This can be obtained by iteration according to the following procedure: 1. Start with a guess for x(q), q m ≤ q ≤ q M ; 2. Use Eqs.(30) and (34) to obtain an estimate of f (q, h) and P(q, h); 3. Get q m from the ratio of Eqs.(35) and (36), computed in q = q m ; 4. Obtain a new guess for q M from Eqs. (36) and (32) computed in q = q M ; 5. Use Eq. (32) to obtain λ(q) and Eq.(37) to obtain x(q); 6. Iterate from 2 until convergence.
The numerical solution of the partial differential equations ( 30) and ( 34) can be obtained by discretizing the profile x(q) through a stepwise function with K steps (this is known as the KRSB approximation) and then increasing K until convergence.The details of this procedure are described in Appendix B.

Distribution of gaps and contacts
From the general relation in Eq. ( 8), and using Eq. ( 33) we get Note that ρ(h) is correctly normalized so that dhρ(h) = 1.Note also that for an observable that is a function of the gaps, we have from Eq. ( 10): In particular the fraction of contacts is given by

The zero temperature phase diagram
The formulae derived in Sec. 3 hold for any value of the control parameters: the density of constraints α, the parameter σ that enters into the constraints, and the inverse temperature β.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 σ ≥ 0 in [15], and for σ < 0 (but |σ| not too large) in [16].Here, we discuss the complete phase diagram for all values of σ and α, with the result given in Fig. 3, and characterize all the phases and phase transitions that appear.

The replica symmetric ansatz
The solution of replica equations like the ones derived in Sec. 3 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 h) and P(q, h) = P(q M , h) = γ q M (h + σ).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 [4,20].In the SAT phase the value of q M remains finite for T → 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 → 0 and therefore one has q → 1 in that limit.

The SAT phase
Taking the zero temperature limit with constant q M < 1, we get where Note that in this case the function −βf RS (q M ) in Eq. ( 41) has a finite limit for T → 0, which gives the entropy of the system, s(q M ) = lim T →0 [−β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. ( 41) with respect to q M or by considering explicitly the RS ansatz in Eq. ( 35).In the last case we obtain q M (1 One can solve Eq. ( 44) 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 → 1. Taking this limit in Eq. ( 44), using the asymptotic properties of the error function, we get the equation for the critical satisfiability threshold α J (σ) within the replica symmetric ansatz: For σ > 0, our optimization problem is convex, the RS solution is always stable (see Sec. 4.1.3),and Eq.(45) gives the correct result for the SAT-UNSAT transition line [15], as shown in Fig. 3.

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 → 1.At very low temperature, the system performs harmonic vibrations around that minimum, and in that case one can show that ).We thus take the T → 0 limit with 1 − q M = χ T → 0 at the same time.Plugging this scaling in Therefore, the ground state energy is and the saddle point equation for χ is given by where α J (σ) is given in Eq. (45).Eq. ( 48) has a solution for χ only for α > α J (σ), which indeed defines the UNSAT phase.In the SAT phase, q M remains less than one for T → 0, while in the UNSAT phase we have q M = 1 − χ T .To match the two scalings, when α → α J (σ) from the UNSAT phase, we should have that χ → ∞, which indeed follows from Eq. ( 48).Furthermore, inserting the saddle point equation ( 48) in the expression of the energy Eq. ( 47), we obtain For fixed σ and α = α J (σ) + δα, Eq. (48) shows that χ −1 ∼ δα when δα → 0, and therefore, close to the jamming line, the energy goes as e ∼ δα 2 .A similar result is obtained at fixed α as a function of δσ: the energy always scales as the square of the distance from the jamming line.
Plugging the RS ansatz in Eq. ( 40), and using the asymptotic scaling of f (q M , h) in Eq. ( 46), one gets the fraction of contacts in the UNSAT phase (at the replica symmetric level): The isostaticity condition is c = αz = 1, and the isostaticity index c is reported in Fig. 2 along the jamming line.For σ > 0 the system is hypostatic on the jamming line, while it becomes isostatic at σ = 0 [16].

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 [20], or to the discontinuous appearance of a RSB solution, usually called a "Random First Order Transition" (RFOT) [30].
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 σ.
The dAT continuous instability of the RS solution can be discussed by computing the eigenvalues of the Hessian matrix, defined as [20] on the RS saddle point solution q a =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. 3. 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 q(x) is continuously becoming different from zero.Slightly in the unstable phase, around this point x, Eq. ( 36) holds.Upon approaching the instability point from the unstable phase, the function q(x) tends to a constant and Eq. ( 36) reduces to its expression computed on the RS solution: This expression, computed on the value of q M that satisfies the RS saddle point equation, gives a condition on (α, σ) that must be satisfied at the dAT instability, hence defining the instability line α dAT (σ) in the phase diagram.Note that computing Eq. ( 37) on this line gives the "breaking point", i.e. the value of x at which the instability occurs.
In the SAT phase (T → 0 at finite q M ), Eq. (52) reduces to In the UNSAT phase (T → 0 with q M = 1 − χ T ), it instead becomes and using Eq. ( 48) that gives the value of χ, we get that the transition line corresponds to σ = 0, ∀α > α J (σ = 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. [15], is that in the whole region σ > 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.

The nature of the RSB phase
We have seen in Sec.4.1.3that 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: 1.The first step is to determine the point at which the RS solution is unstable; this is the dAT line α dAT (σ) determined in Sec.4.1.3.At the dAT line, the function q(x) is not constant anymore.As mentioned in Sec.4.1.3,the breaking point m of the RSB solution, i.e. the point x = m where q(x) becomes different from zero, can be computed plugging the RS ansatz into Eq.(37).
2. The next step is to evaluate whether m < 1 or m > 1.In fact, because the function q(x) is defined for x ∈ [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 [4,30].We will discuss further this situation in Sec.4.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 ∈ [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 q(m).In fact, a fullRSB solution requires q(m) > 0, because q(x) must be an increasing function1 of x.If this is not the case, then the transition is a continuous transition to a 1RSB solution.To compute q(m), we derive Eq. ( 37), expressed as a function of x, with respect to x.We get By evaluating Eq. ( 55) 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.

The SAT phase
In the SAT phase, the breaking point is where f SAT (q M , h) is defined in Eq. ( 42), and q M is the solution of the saddle point equation (44) computed at α = α dAT (σ).A numerical evaluation of m as a function of σ < 0, on the dAT line in the SAT phase, shows that close enough to σ = 0 one has m < 1, while m increases upon decreasing σ towards more negative values.At σ = σ RFOT one has m = 1.Beyond that point, for σ < σ 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 σ < σ RFOT we need to study the 1RSB solution more carefully, see Sec. 4.3.Furthermore, the analysis of q(m) on the same transition line shows that close enough to σ = 0, for σ > σ 1RSB , one has q(m) > 0 and the dAT instability leads to a fullRSB phase, while for σ ∈ [σ RFOT , σ 1RSB ] one has 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.

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 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 ∝ 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. 5.Moreover, the zero temperature limit of the slope q(m) at the breaking point, in the UNSAT phase, is given by q(m) which is of order 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 [16], see Fig. 3.We will see in Sec. 5 that the scaling with T of both m and q(m) is in perfect agreement with the complete scaling form of q(x) in the UNSAT phase.

The 1RSB free-energy in the SAT phase
In Sec.4.2 we have shown that for σ < σ RFOT the solution for q(x) in the SAT phase becomes of the RFOT type (i.e. a discontinuous 1RSB solution) [30].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 ∈ [0, m] and to q M for x ∈ [m, 1] (see Appendix B for a general KRSB solution).At zero temperature, in the SAT phase where both q 0 < 1 and q 1 < 1, the function −βf[x(q)] has a finite limit, which coincides with the 1RSB entropy, given by The variational equations always admit the RS solution q 1 = q 0 .For any given σ this is the unique solution for sufficiently small α.For fixed σ < σ RFOT and upon increasing α, there is a point, called the dynamical transition α dyn (σ), where a non-trivial solution with q 1 = q 0 appears for m → 1.We do not discuss in details the properties of this transition, for which we refer to [4,30,32].In short, at α dyn (σ) the Boltzmann-Gibbs measure on the space of solutions splits into an exponential number of clusters of solutions [4].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 α dyn (σ), which is thus not a thermodynamic phase transition [4,30].Increasing α, the complexity of the relevant The zero temperature phase diagram of the perceptron.α J (σ) is the jamming transition (or SAT/UNSAT threshold) within the replica symmetric approximation.The replica symmetric solution is stable for all σ > 0. Thus the jamming transition line is exact only for σ > 0. α dAT (σ) are the lines where replica symmetry breaks.In the non-convex region and for σ ∈ [σ 1RSB , 0], α dAT is a transition from a replica symmetric phase to a continuous fullRSB one.For σ ∈ [σ RFOT , σ 1RSB ] the dAT line is a continuous transition from a replica symmetric phase to a stable 1RSB phase.For σ < σ RFOT a RFOT type phenomenology is observed.Keeping σ fixed and increasing α, one has first a dynamical transition for α = α dyn (σ), then a Kauzmann transition for α = α K (σ), and finally a Gardner transition for α = α G (σ).In the UNSAT phase, the instability line α dAT (σ) represents the transition from a replica symmetric phase to a continuous fullRSB one for all α > 2. Therefore, fullRSB occurs in the region delimited by the α dAT (σ) and α G (σ) lines, which contains the whole jamming line for σ < 0. clusters decreases and vanishes at the Kauzmann (or condensation) transition α = α K (σ) [4,30,32].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 [27,[33][34][35][36] for reviews).
To compute α dyn (σ) and α K (σ), we are therefore led to consider the 1RSB entropy for m ≈ 1, where we can write to the leading order s 1RSB (q 1 , q 0 , m) = s RS (q 0 ) + (m − 1)δs RSB (q 0 , q 1 ) , (61) with s RS given in Eq. ( 41), and When m → 1, q 0 can be obtained by setting to zero the derivative of Eq. ( 41), which gives the saddle point equation ( 44).The dynamical transition point α dyn (σ) corresponds to the lowest value of α where the variational equation for q 1 , obtained by setting the derivative of Eq. ( 62) to zero, admits a solution q 1 > q 0 .At the Kauzmann transition point α K (σ), the breaking point m, obtained by setting the derivative with respect to m of Eq. ( 60) to zero, first becomes smaller than one.Upon further increasing α, it can be shown that the 1RSB solution becomes unstable and undergoes a Gardner transition towards a fullRSB phase [37].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.4.1.3for 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 [37] (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 α G (σ) of the 1RSB transition is given by the following condition: where The dynamical, the Kauzmann and Gardner transition lines are plotted in Fig. 3.

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 σ > 0 has been discussed already in [14,15] 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 [7] does not lead to marginal stability.In this section we want to describe the jamming transition for σ < 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 [14].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 [13] 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.
5.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 → ±∞ which are generically valid, in particular in the scaling solutions of our interest.For h → ±∞ the boundary condition for f (q, h) reduces to Using Eq. ( 30) one readily finds For h → ∞, Eq. ( 34) becomes simply Ṗ = P /2, which has a unique solution compatible with the boundary condition For h → −∞ the equation is slightly more complicated; but still, a Gaussian form of the kind P(q, h) ∼ D(q)e −D(q)h 2 satisfies Eq. ( 34) for large negative h, where D(q) verifies Ḋ(q) = −2D(q) 2 + 2D(q) β x(q)

Approaching jamming from the SAT phase
In the SAT phase we can take the limit T = 1/β → 0 by simply replacing e −β v(h) → θ (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 [20].In the jamming limit, the cluster volumes go to zero and therefore the self-overlap of a cluster, q M , goes to one.

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 [13] (see Sec. 3.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 ε is the linear distance from the jamming line (either in α or σ).In addition we assume that q M = 1 − ε κ where κ is an exponent to be determined by the equations.In the jamming limit ε → 0, one has y ∈ [0, 1/ε] → [0, ∞) and q ∈ [q m , q M ] → [q m , 1].We wish to show that a scaling solution exists in this limit, when y is large and q ∼ 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 σ and α, on the jamming line the function p 0 as well as the exponents a and κ are universal (i.e.independent of the precise location on the jamming line).In Sec.5.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 ε, the scaling solution (70) is cutoff when q ∼ q M = 1−ε κ and y ∼ y M = Y ε −1 .Finally, note that while we will be able to prove the existence of the scaling solution, and compute the values of a and κ, we will not be able to prove that ε is proportional to the distance from the jamming line: at present, this follows from the numerical solution of the equations [13].

Proof of the scaling form
The scaling analysis is carried out along the lines of [13,14,38].
Scaling of f (q, h) -The function m(q, h) introduced in Eq. ( 70) satisfies the equation where the differential equation comes from Eq. ( 30) and the asymptotes from Eq. (66).Let us inspect the value of for ε → 0 and q → 1, according to Eq. (70): This expression has a crossover for 1 − q ∼ ε κ , when the two terms in the denominator are of the same order.The scaling form is obtained for 1 − q ε κ , in which case we have Plugging this result and the scaling form (70) in Eq. (5.2.2), we obtain a scaling equation for the function (t): This non-linear equation, with boundary conditions at t → ±∞, admits a unique solution for each value of κ.

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. ( 67) it follows that for h → ∞, P(q, h) remains a finite Gaussian, P(q, h → ∞) ∼ γ q (h + σ).For finite h > 0, the Gaussian will be deformed to a finite function p + (h) as it appears in Eq. ( 70).For h → −∞, on the other hand, we have from Eq. ( 68) that for β → ∞ and in the scaling regime of Eq. ( 73): If κ < 2 (which, we will see, is the case), this equation admits a scaling solution where the term 2D(q) 2 is negligible.One concludes therefore that in the scaling regime, P(q, h → −∞) ∼ D(q)e −D(q)h 2 has the form which has been proven asymptotically for h → −∞ but can be extended to the whole regime where q ∼ 1 and |h| ∼ (1 − q) (κ−1)/κ , as in Eq. ( 70).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 [38].
The "matching" regime with p 0 (t) must be introduced in Eq. (70) 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/ 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. ( 34), to a non-trivial equation for p 0 (t), namely: Recall that in Eq. (5.2.2), (t) depends on κ; it turns out that Eq. (5.2.2) also admits a unique solution for p 0 (t) satisfying the correct asymptotic conditions, but only for a given choice of a = a(κ).For a given κ, Eqs.(5.2.2) and (5.2.2) thus determine (t), p 0 (t) and a.
Determination of the exponent κ -The exponent κ can be fixed using Eq. ( 37) which can be equivalently written as In the scaling regime, Eq. ( 73) gives the left hand side.In the right hand side, we can note that m (q, h) and m (q, h) 2 [1 + m (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 ∼ 1 − q.We obtain Because both (t) and p 0 (t) depend on κ, this is an equation for κ.The numerical solution of the system of Eqs.(5.2.2), (5.2.2) and (82), gives κ = 1.41574 . . .while within the numerical precision one finds the relation The importance of these "scaling" relations will be further discussed in Sec.5.5 and Sec.5.6.

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 → 0 in the UNSAT phase.In Sec.4.2 we have shown that in the UNSAT phase, close to the dAT instability line (σ = 0, α > 2) of the RS solution, the solution q(x) has the following properties: • The Edwards-Anderson order parameter is q M = 1 − χ T • The breaking point at the instability transition line is m = m T • The slope of q(x) at the breaking point on the instability line is q(m) ∼ T Additionally, here we want to show that, like in the Sherrington-Kirkpatrick model [20, 39]: • P(q, h) is smooth and regular at zero temperature • The function β 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 (x) defined in Eq. (57).In the fullRSB region, as in the RS case, we define q M = 1−χ T and we introduce As in Sec.5.2.2, we introduce m(q, h) = λ(q) f (q, h) = λ(q) f (q, h), which satisfies Eq. (5.2.2) with the modification m(q, h → −∞) = −hχ λ(q)/(1 + χ λ(q)) and initial condition The equation for P(q, h) is still Eq.(34).Finally, Eq. ( 36) becomes and Eq. ( 37) becomes identical to Eq. ( 81).The distribution of gaps, given in Eq. ( 38), becomes Because the breaking point m = x(q M ) ∝ T , one has y(q M ) ∼ 1/ T and thus y(q) extends up to infinity in the zero temperature limit.The scaling behavior of y(q) for q → 1 is expected to be [39]: We now check that this scaling is consistent.First of all, this implies that q(x) ∼ 1 − AT 2 /x 2 for some constant A. For x = m T , we get q M = 1 − χ T which matches the behavior on the instability line.Also, one has q(x) ∝ T 2 /x 3 and then q(m) ∝ T , as it is the case along the instability line.The final check can be obtained from Eq. ( 81).We assume that the scaling behavior of m(q, h) for q → 1 is which agrees with the initial condition (86).Plugging this ansatz inside Eq. ( 81) we get which confirms the consistency of Eq. ( 89) and provides an explicit expression of y χ .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 → ∞ we have P(q, h) → γ q (h + σ), which suggests that for h → ∞, P(1, h) is finite and has smooth corrections in q.Next, we can observe that Eq. (68) becomes for q → 1: which therefore indicates that for h → −∞, P(1, h) is finite and has corrections proportional to 1 − q.Plugging P(q, h) ∼ P(1, h) + 1 − q δP(h) in Eq. ( 34), and using Eq.(90), we obtain which satisfies the condition dh δP(h) = 0 and shows that there are 1 − q corrections only in the region h < 0.

The jamming limit from the UNSAT phase
The jamming limit from the UNSAT phase is obtained by considering the T = 0 equations of Sec.5.3 in the limit χ → ∞, 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 − χ T for T → 0: matching the two regimes for T ∼ 0 requires the divergence of χ at the jamming point.
We know that in the UNSAT phase for q → 1 we have, from Eq. ( 89) and (85): If χ is finite and q → 1, then in the denominator the factor 1 dominates and y/ λ ∼ (1 − q) −1/2 : in this regime we recover the "regular" zero temperature solution of the UNSAT phase described in Sec.5.3.Conversely, for χ → ∞, the coefficient y χ → ∞.In fact, from Eq. ( 91) we obtain y χ ∝ χ P(1, 0) observing that (t) is finite and that 0 −∞ dhP(1, h) must remain finite because P(1, h) is normalized to 1.In this case, for all q < 1, at large enough χ the second term in the denominator of Eq. ( 94) dominates and one has , for q ∼ 1 and In this regime, we must have a different scaling solution in which y/ λ ∝ 1/(1 − q): but this is exactly the jamming scaling solution that was already derived in Sec.5.2, which indeed must emerge from the regular zero temperature UNSAT solution upon approaching the jamming point.
To summarize, when χ → ∞ and in the region of q → 1, we expect two different scaling solutions: when 1 − q 1 − q * , we have the "regular" UNSAT scaling of Sec.5.3; for 1 − q 1 − q * we have instead the "jamming" scaling solution of Sec.5.2.The matching point q * is determined by the condition that χ P(1, 0) 1 − q * ∼ 1, and the value of y(q * ) is then simply Note that in the jamming solution, Eq. (70), we have P(q * , 0) ∼ (1 − q * ) −a/κ , and we know from Eq. ( 93) that in the regular solution P(q, h) changes by a very small amount, ∼ 1 − q, when q → 1.We conclude that P(1, 0) ∼ (1 − q * ) −a/κ and more generally 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.

Scaling relations between exponents
The relation a = 1 − κ/2, and its consequences that are given in Eq. ( 83), has been found in Ref. [13] within numerical precision by solving the equations for the critical exponents derived in Sec.5.2.2.Eq. ( 83) is physically very important: in fact, the relation γ = 1/(2+θ ) has been proven in [7] to be a direct consequence of marginal stability (in the perceptron, the same marginal stability argument is discussed in [14]).It would therefore be nice to have a more direct analytical proof of the relation a = 1−κ/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.5.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 ∝ 1/χ 2 in the jamming limit χ → ∞.Using now P(1, h) ≈ (1 − q * ) (1−κ)/κ p − ((1 − q * ) (1−κ)/κ h) for h < 0, which was derived in Sec.5.4, together with Eq. (88), we find Because we know that [h] ∝ 1/χ 2 , we must have This is compatible with (98) only if a = 1−κ/2, which gives an independent proof of Eq. ( 83).

Scaling of several interesting observables at the jamming transition
We have now fully characterized the scaling of the basic quantities, x(q), λ(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.( 69) and ( 70), together with Eqs. ( 78), ( 79) and (83).On the UNSAT side, our main results are Eqs.(85), (87), (88) and (101).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.

Energy, pressure, number of contacts
We start by considering the energy, the pressure and the number of contacts defined in Eq. (11).From Eqs. (88) and (101), 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. (87) 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 ∝ χ −1 ∝ p 1/2 .To obtain the complete phenomenology of jamming, one should also prove that the pressure p ∝ χ −2 vanishes linearly in the distance from jamming in the (α, σ) 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 ∝ p 2 and c − 1 ∝ p 1/2 perfectly agree with numerical observations in jammed sphere packings [5,6].
In the SAT phase, at zero temperature, there are by definition no negative gaps, ρ(h < 0) = 0; consequently pressure, energy, and contacts all vanish.However, one can study the limiting values of these quantities when T → 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 σ, as it can be checked directly from the definition of the partition function in Eq. (17).Using the replica expression of the free energy in Eq. ( 31), one then obtains: Using the equations for P and f , it is possible to show that The proof is obtained by writing the derivative as dh( Ṗ f + P ḟ ), using the equations for P and f and integrating by parts [24].Therefore, in the SAT phase the pressure vanishes proportionally to T , a standard result for hard spheres (which correspond to the T → 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. ( 70), noting that λ(1) = 1 − q M = ε κ , one has The scaling form in Eq. ( 70) is cutoff when 1 − q ∼ ε κ , and therefore the same scaling holds for P(1, h) and m(1, h) with 1 − q → ε κ .Then, in the ε → 0 limit the region h > 0 does not contribute to the integral because m(1, h) → 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) → −h and We conclude that the reduced pressure p/T diverges proportionally to ε −1 , and therefore 1 − q M ∝ (p/T ) −κ , a result that has been numerically tested in hard sphere systems [13].This provides a physical meaning for the exponent κ.In a similar way one can study the scaling of the energy and the number of contacts in the SAT phase.

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. ( 12).The distribution of small positive gaps and small contact forces has been associated with important properties of the packings, including marginal stability [7,8].In this section we discuss these distributions.They can be straightforwardly derived from Eq. (88) and Eq. ( 101) which together imply, for χ 1: At jamming, when χ → ∞, the distribution of gaps concentrates on h ≥ 0, where it is given by ρ where /α according to Eq. ( 104).Isostaticity is expressed by the presence of a delta peak in ρ(h) whose weight is the fraction of contacts.Moreover, this implies that the distribution of small positive (satisfied) gaps, according to Eq. ( 78), behaves as a power-law ρ(h → 0 + ) ∼ h −γ with γ = (2 − κ)/κ = 0.41269 . ... Similarly, the normalized probability distribution of scaled contact forces f s µ = −[1]h µ /p, with p/ [1] given by Eq. ( 103), can be deduced from ρ(h < 0) and reads This implies, according to Eq. ( 79), that the distribution of small forces behaves as a power-law, p( f → 0 + ) ∼ f θ , with θ = 1/γ−2 = 0.42311 . . .These results give a physical interpretation to the critical exponents γ and θ , and connect with the results of [7][8][9][10]12].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 [9,25,38].As a check of the theoretical predictions we have extended the numerical simulations of [14] to large systems: in Fig. 4    We conclude by mentioning that one can also study the vibrational spectrum, both in the SAT and UNSAT phases [16,17], 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 [40].Also, it has been shown in [19] that the asymptotic behavior of the function y(q), related to the exponent κ, controls the scaling of the avalanches at the jamming point, which is therefore different than the scaling in the UNSAT phase.

Conclusions
In this work we have formulated the jamming problem as a general constraint satisfaction problem with continuous variables (Sec.2).We have then specialized on the random perceptron, a well known machine learning model, which is a prototype of this class of problems (Sec.3).In the non-convex regime, the model shows a complex zero temperature phase diagram in the plane of the two control parameters (α, σ), which has been fully characterized in Sec. 4. In particular, we have shown that for σ < 0 and large enough |σ|, the phase diagram as a function of α shows the characteristic phenomenology associated with the Random First Order Transition (RFOT) mean field theory of glasses [27,[33][34][35][36].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. 5, 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 [13].Furthermore we have extended the study of [13] 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 [5].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 [7].Our results are also consistent with the scaling analysis of [11].These results, together with the results on the vibrational spectrum obtained in [16,17], and the study of avalanches performed in [19], 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 [13].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 [41].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.

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. ( 17), and neglecting proportionality constants, the replicated partition function can be written as ∝ e where to obtain the last line we integrated over ra to obtain a δ(r a − k a ) and then integrated over r a .Introducing h a = k a − σ and setting α = M /N , we obtain the final result for the replicated partition function that is reported in Eq. ( 21).

B Numerical solution of the equations
x q 1 q 2 q 3 q 4 q 0 m 4 = 1 In this Appendix, we derive the equations corresponding to a discrete KRSB ansatz, illustrated in Fig. 5.This is particularly useful for the numerical solution.
For a function q(x) with KRSB form, the free energy of Eq. ( 31) is s[q(x)] = lim with 〈q〉 = 1 0 dx q(x), and λ(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. ( 30): The saddle point equation is and the matrix M a b is also a hierarchical matrix whose components can be written as with P(m 0 , h) = γ q 0 (h + σ) We will prove below that in the zero temperature limit and in the fullRSB phase, µ = 1/χ 2 , so that the following exact zero temperature relation holds:
We now compute µ 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 µ.We have s[q d , q(x)] = 1 2 log (q d − q M ) + q m q d − 〈q〉 + 1 0 dx q(x) λ(x) + αγ q m f (0, h)| h=−σ + 1 2 βµ(q d − 1) − α dh P(q M , h) [ f (q M , h) − log γ q d −q M e −β v(h) ] + α dh q M q m dq P(q, h) [ ḟ (q, h) + 1 2 f (q, h) + x(q) f (q, h) 2 ]. λ(q) = q d − q M + q M q dp x(p) . (139) 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 f (q M , h) = log γ q d −q M e −β v(h) (140) and λ(q) is given by Eq. ( 139).Note that the variational equation with respect to µ 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 (141) It is very easy to show that the last term of the equation above can be rewritten as where again m(q M , h) = ∂ f (q M , h)/∂ h.At this point we can use the saddle point equations ( 35) and ( 36) to write q m λ 2 (q m ) + q M q m dp λ 2 (p) = α dhP(q M , h)m 2 (q M , h) , 1 (q d − q M ) 2 = α dhP(q M , h) m (q M , h) 2 . (143) Using that the variational equation over µ gives q d = 1, that q M = 1 − χ T and that we get µ = 1/χ 2 .

Figure 2 :
Figure2: Isostaticity index c on the jamming line, computed on the RS solution for σ ≥ 0, and on the fullRSB solution for σ < 0. The system is hypostatic for σ > 0 and isostatic for σ ≤ 0.

Figure 3 :
Figure3: The zero temperature phase diagram of the perceptron.α J (σ) is the jamming transition (or SAT/UNSAT threshold) within the replica symmetric approximation.The replica symmetric solution is stable for all σ > 0. Thus the jamming transition line is exact only for σ > 0. α dAT (σ) are the lines where replica symmetry breaks.In the non-convex region and for σ ∈ [σ 1RSB , 0], α dAT is a transition from a replica symmetric phase to a continuous fullRSB one.For σ ∈ [σ RFOT , σ 1RSB ] the dAT line is a continuous transition from a replica symmetric phase to a stable 1RSB phase.For σ < σ RFOT a RFOT type phenomenology is observed.Keeping σ fixed and increasing α, one has first a dynamical transition for α = α dyn (σ), then a Kauzmann transition for α = α K (σ), and finally a Gardner transition for α = α G (σ).In the UNSAT phase, the instability line α dAT (σ) represents the transition from a replica symmetric phase to a continuous fullRSB one for all α > 2. Therefore, fullRSB occurs in the region delimited by the α dAT (σ) and α G (σ) lines, which contains the whole jamming line for σ < 0.
we compare the cumulative distribution C f = f 0 d f p( f ) for sizes from N = 50 to N = 1600 with the theoretical prediction.

Figure 4 :
Figure 4: Plot of the cumulative of the force distribution C f = f 0 d f p( f ) for N = 50, . . .1600.The line is proportional to the theoretical prediction C f ∼ f 1+θ .

1 2 1 2 1 2 2 ∂ 1 2 2 ∂
−1/2 X a • ξ µ ) aµ e −β v(r µ a −σ) ,(112)where a = 1 • • • n and µ = 1 • • • M ; one can check that integrating over rµ a produces delta functions that fix rµ a as in the original partition function.Next, the Gaussian integral over the quenched disorder ξ µ gives e −N −1/2 aµ ir µ a X a • ξ µ = e − dr a e a ir a r a − ab ra rb Q ab −β a v(r a −σ) dr a e a ir a r a − ab ra rb Q ab −β a v(r a −σ) a change of variables from X a to Q ab , with Jacobian exp( N 2 log det Q).Finally, one can easily show by developing both sides in powers of Q ab that e − ab ra rb Q ab = e k a ∂ k b e − a k a ir a dr a e a ir a r a − ab ra rb Q ab −β a v(r a −σ) k a ∂ k b e a (r a −k a )ir a −β a v(r a −σ)

Figure 5 :
Figure 5: An example of the parametrization of the matrix Q ab for a 4RSB case.The reader should keep in mind that we use the convention m 0 ≡ m, m k+1 ≡ 1.