Statistical mechanics of the maximum-average submatrix problem

We study the maximum-average submatrix problem, wherein given an N × N matrix J, one needs to find the k × k submatrix with the largest average number of entries. We investigate the problem for random matrices J, whose entries are i.i.d. random variables, by mapping it to a variant of the Sherrington–Kirkpatrick spin-glass model at fixed magnetisation. We analytically characterise the phase diagram of the model as a function of the submatrix average and the size of the submatrix k in the limit N→∞ . We consider submatrices of size k=mN with 0<m<1 . We find a rich phase diagram, including dynamical, static one-step replica symmetry breaking (1-RSB) and full-step RSB. In the limit of m → 0, we find a simpler phase diagram featuring a frozen 1-RSB phase, where the Gibbs measure comprises exponentially many pure states each with zero entropy. We discover an interesting phenomenon, reminiscent of the phenomenology of the binary perceptron: there are efficient algorithms that provably work in the frozen 1-RSB phase.

We consider the maximum-average submatrix (MAS) problem, i.e. the problem of finding the k × k submatrix of an N × N matrix J with the largest average of entries.This is a natural combinatorial optimization problem that has been studied in the mathematical and data science literature [1], mainly in the context of biclustering [2].Theoretical works focused on the case where J is a random matrix (i.i.d.standard Gaussian entries), the size of J is large (N → ∞) and the size of the submatrix k ≪ N [3][4][5].
From a statistical physics point of view, the MAS problem is a natural variant of the well-known Sherrington-Kirkpatrick model [6] with spins σ i ∈ {0, 1} at fixed magnetization.Statistical physics of disordered systems and the related replica method [7] have been used widely to study other combinatorial optimization problems such as graph partitioning [8], matching [9], graph colouring [10], K-satisfiability of Boolean formulas [11], and many others.As far as we are aware, the maximum-average submatrix problem has not been studied from the statistical physics point of view.Filling this gap is the main purpose of the present paper.
Our results reveal the exact phase diagram of the MAS problem when k = mN , N is large and m finite.We unveil that at large values of m as the submatrix average increases the system undergoes a continuous phase transition to a full replica symmetry breaking (RSB) phase [7].At intermediate values of m the phase transition becomes discontinuous, passing through a dynamical onestep RSB and static one-step RSB phases to a full-RSB one [7].At yet lower m, the full-RSB phase then vanishes and the maximum average is given by the one-step RSB solution.In the limit of m → 0, the MAS problem behaves in a way related to the random energy model [12] presenting frozen one-step RSB [13,14].
We also find that in the limit m → 0 the phase diagram presents a region where polynomial algorithms are proven to work [4] yet according to our results the equilibrium behaviour of the problem is given by the frozen one-step RSB phase that is considered algorithmically hard [15].One other such problem is known in the literature -the binary perceptron.For the binary perceptron, an explanation of the discrepancy between equilibrium properties and algorithmic feasibility has been proposed in relation to out-of-equilibrium large-local-entropy regions of the phase space that are not described within the standard replica solution [16].This finding has been used to discuss learning in artificial neural networks [17] and to propose new algorithms [18].The analogy of behaviour between the binary perceptron and the MAS problem is therefore interesting as it may serve to shed more light on the fundamental question of algorithmic hardness.From the point of view of the mathematics of spin glasses, the perceptron problem is difficult to handle due to the effective bipartite structure of the correlations.The MAS problem belongs instead to a class of problems for which the exactness of the replica calculation has been established rigorously in [19].
We now review the mathematical results which we later connect to our analysis.All these results hold in the regime k ≪ N .In [3], the authors proved that the globally optimal submatrix has an average equal to A opt = 2 log N/k (log is the natural logarithm).They also conjectured, and it was later proven by [4], that the Largest Average Submatrix (LAS) algorithm -an efficient iterative row/column optimization scheme -fails to reach the global optimum, as its fixed point -akin to local minima -has with high probability average equal to A LAS = 2 log N/k.The fact that the LAS algorithm fails bears a natural question: is A LAS an algorithmic threshold signalling the onset of a hard phase?In [4], the authors introduced a new algorithm called Incremental Greedy Procedure (IGP) which is able to produce submatrices with average A IGP = 4/3 2 log N/k > A LAS .Additionally, they proved that for averages larger than at least A OGP = 10/(3 √ 3) log N/k > A IGP the problem satisfies the Overlap Gap Property (OGP) [20].This means that i) the LAS threshold A LAS seems to be In the central and right panel we rescale the sub-matrix average as a/ log(1/m)/(m) to highlight the convergence to the limit.We identify five distinct phases.In the RS phase (green) the system is replica symmetric.In the 1-RSB phases replica symmetry is broken to one step and two sub-phases exists, a dynamical 1-RSB with an extensive number of equilibrium pure states (blue) and a static 1-RSB with only finitely many pure state (purple).All the phase boundaries are exact in the thermodynamic limit except the boundary between full-RSB (orange) and UNSAT (red) which would require solving the full-RSB equations.In the full-RSB phase (orange) replica symmetry is completely broken and the set of pure states manifests ultrametricity.In the unsatisfiable phase (UNSAT, red) no submatrix exists with the given values of a and m.The transition from RS and 1-RSB to full-RSB is continuous and caused by an instability of the 1-RSB ansatz (dashed line), while the other transitions are discontinuous.We observe two tricritical points, one at (mc, ac) where the system shows coexistence of RS, 1-RSB and full-RSB phase (white marker), and one at (m * , a * ) where the largest-average submatrices become 1-RSB stable and the full-RSB region ceases to exist (black marker).In the limit m → 0, we observe only the RS, 1-RSB and UNSAT phases.The 1-RSB phase is frozen, meaning that the internal entropy of each pure state goes to zero in the m → 0 limit.an algorithm-specific threshold, and not a more general trace of an intrinsic computational-to-statistical gap, and ii) that the problem will likely exhibit a hard phase preventing algorithms to find submatrices with averages larger than at least A OGP .[21] discusses additional results for k = N γ , 0 < γ < 1.

I. THE MODEL
We consider an N r × N c random matrix J composed of i.i.d.Gaussian entries with zero mean and unit variance.A k r × k c submatrix σ of J is defined by two arbitrary subsets of rows and columns I r , I c such that J ij belongs to the submatrix σ if and only if i ∈ I r and j ∈ I c , and the cardinalities satisfy |I r | = k r and |I c | = k c .There are three versions of the MAS problem: • Rectangular MAS: N r , N c , k r and k c are unconstrained.This is the problem relevant for applications [1,2].
• Square MAS of a square matrix: This version is studied in the mathematical literature [3,4,21].
• Principal MAS of a symmetric matrix: T is symmetric and we consider only principal submatrices, i.e. submatrices for which I r = I c .
In the following, we focus on the case of the principal MAS of a symmetric random matrix motivated by its close relation to the SK model.In the SM, we sketch the corresponding solution for the rectangular MAS of a random matrix and realize with surprise that the equations leading to the phase diagram of the square MAS of a square random matrix are exactly the same as the ones for the principal MAS of a symmetric random matrix.This allows us to compare our results directly to the mathematical literature, and it also means that the phase diagram we provide applies to the non-symmetric case.We believe that the physics of the rectangular MAS problem has two more hyperparameters N r /N c and k r /k c leading to a 4-dimensional phase diagram.Its exploration is left for future work.We encode principal submatrices, i.e. their row/column index set I, as Boolean vectors σ = {σ i } N i=1 ∈ {0, 1} N such that, if i ∈ I then σ i = 1 and vice versa.We fix the size of the submatrix to k = mN , which in the Boolean representation translates to the condition i σ i = mN .We call m ∈ (0, 1) magnetization.The average of the entries of a submatrix σ can be then expressed as We define a = A √ N , and we will see that a is of order one in the thermodynamic limit.
We probe the energy landscape of the MAS problem by studying the associated Gibbs measure where β is an inverse temperature that we use to fix the average energy, h is a magnetic field that we use to fix the magnetization m and Z(β, h) is the partition function.
The uncommon plus sign in front of the inverse temperature is due to the fact that the problem is a maximisation problem.Thus, small temperatures correspond to large positive energies in this model.The energy function is defined as which, modulo subleading contributions coming from the diagonal term, is a multiple of the submatrix average a.
As the considered model resembles the classic SK model, note that the mapping of the Boolean spins to ±1 spins, i.e. s = 2σ − 1, leads to an SK model in a random magnetic field, with couplings correlated to the magnetic field (see SM).Such a model has not been considered in the physics literature as far as we are aware.
We will compute all thermodynamic observables through the quenched free entropy, i.e.Φ = lim N →∞ E J log Z(β, h)/N where E J denotes averaging over the distribution of J.We will see that the free entropy can be expressed as a variational problem for the overlap order parameter, which is defined as q = N −1 N i=1 σ a i σ b i ∈ [0, m] for two replicas of the system σ a and σ b .

II. REPLICA ANALYSIS OF THE FREE ENTROPY
We compute the quenched free entropy Φ using the replica formalism [7], i.e. by using the replica trick E J log Z = lim n→0 (E J Z n − 1)/n, computing E J Z n for integer values of n and performing an analytical continuation to take the n → 0 limit.We perform the analytical continuation under the one-step Replica Symmetry Breaking (1-RSB) ansatz, in which we assume that the Gibbs measure decomposes into a 2-level hierarchy of pure states.This hierarchy is characterised by two overlaps: the average overlap between microstates belonging to the same pure state q 1 , and the one between microstates belonging to different pure states q 0 .The Parisi parameter p acts as a temperature controlling the trade-off between the free entropy of a single pure state, and the entropy of pure states [22].After a derivation that follows steps standard to the replica method [7,23], we obtain the following variational free entropy: where Du and Dv denote integration against a standard Gaussian measure.The variational free entropy depends on the submatrix size/magnetization m, the intra-state overlap q 1 , the inter-state overlap q 0 and the Parisi parameter p.To obtain the equilibrium free entropy we extremize the variational free entropy over m, q 0 and q 1 .The Parisi parameter must be set to one if the resulting complexity (whose definition we provide in the following) is positive, otherwise, the variational free entropy must also be extremized over p.
Under the 1-RSB ansatz we have the following expressions for the observables, to be evaluated at the equilibrium values of the order parameters.The average energy density (and the submatrix average) equals the total entropy (logarithm of the number of microstates) equals and the complexity (logarithm of the number of pure states contributing to the Gibbs measure) equals If the complexity is non-zero, the total entropy decomposes as s total = s internal + Σ, where the internal entropy s internal is the logarithm of the number of microstates contributing to each of the exponentially many pure states.Finally, we need to investigate whether the 1-RSB result is exact in the thermodynamic limit.This is done by analyzing the stability of the 1-RSB ansatz against perturbation of higher-order RSB nature.We perform the so-called Type-I stability analysis (see SM Section I.d), obtaining the stability condition where ℓ(x) = 1/(1 + exp(−x)).When this condition is not satisfied then the 1-RSB results are just an approximation and more steps of RSB need to be taken into account.
We derived the variational free entropy using the replica trick.Note, however, that the proof of the full-RSB free entropy giving the exact solution in the thermodynamic limit from [19] applies to the MAS problem and thus our setting.In order to apply their result to our model we note that [19] constrains the free-entropy to fixed self-overlap q self = N −1 i σ 2 i , while we constrain the model to fixed magnetization m = N −1 i σ i .Due to the choice of Boolean spins, we have that σ 2 i = σ i , so that the two constraints coincide (this is not true in general).

III. THE PHASE DIAGRAM
After solving the above equations, we identify five distinct phases for finite m ∈ (0, 1), and we plot them in Figure 1.
RS phase -For small submatrix-average (corresponding to large temperatures) we observe a replicasymmetric (RS) phase, in which the extremum of the variational free entropy is attained at q * 0 = q * 1 .In this phase, the complexity is zero, while the total entropy is strictly positive.As the submatrix average increases (i.e. the temperature is lowered), the system undergoes a phase transition to an RSB phase.The nature of the transition is different for m ≤ m c and m ≥ m c .We start discussing the former case.
Dynamical 1-RSB phase -For m ≤ m c , we observe a discontinuous transition at a value a dynamic (m) of the sub-matrix average from the RS phase to a dynamical 1-RSB phase, in which the extremum of the variational free entropy satisfies q * 0 ̸ = q * 1 .This transition is identified by a sharp jump of the complexity from zero to a positive value, meaning that the measure shatters into an exponential number of pure states each with non-zero entropy.As m decreases, the appropriately rescaled internal entropy decreases, suggesting that in the m → 0 limit, this phase becomes a frozen 1-RSB phase, see Figure S2 in SM.
Static 1-RSB phase -For m ≤ m c , in the dynamical 1-RSB phase, the complexity continuously decreases as the submatrix average increases.At the value of submatrix average a static (m) > a dynamic (m) at which the complexity vanishes we observe a first-order phase transition, from the dynamical 1-RSB phase to a static 1-RSB phase (q * 0 ̸ = q * 1 , zero complexity, positive entropy).Full-RSB phase -For values m ≥ m c we observe a continuous phase transition from RS to full-RSB phase at a stability (m).The transition happens when the Type-I 1-RSB stability condition (8) fails.We conjecture that this phase is fully replica-symmetry-broken by the similarity to the SK model, in which a similar continuous phase transition to full-RSB occurs.For values of m c ≥ m ≥ m * , we observe a Gardner-like phase transition at a value a stability (m) of the sub-matrix average from the 1-RSB phase to a full-RSB phase.This transition is reminiscent of the one known from the Ising p-spin model [24].
UNSAT phase -As the average of the matrix increases, we encounter a point at which the total entropy vanishes, denoting that the sub-matrix average has reached its maximum value a max (m).After this point, the total entropy becomes negative and we observe an unsatisfiable (UNSAT) phase, where no submatrix with that value of the submatrix-average exists.For m ≥ m * , we provide only an approximate estimate of this transition line (while all other transitions presented are exact up to the precision of the numerical solver).We estimated a max (m) in the 1-RSB solution, even though this ansatz is unstable in this phase, by computing the 1-RSB entropy, finding the temperature at which it vanishes, and computing the corresponding submatrix average.
It is often the case that the 1-RSB prediction for the maximum energy is numerically very close to the full-RSB prediction.Evaluation of the the full-RSB equations, which are proven to give the correct maximum average (analogous to the ground-state energy in the SK model) [19], is left for future work.
Tricritical points -The phase diagram features two tricritical points.The first one, at (m c , a c ) ≈ (0.09 − 0.1, 9.3 − 9.7) marks the coexistence of the RS, 1-RSB and full-RSB phases.It can be pinpointed by finding the intersection of the stability and the static transition lines.As m → m − c , the static and dynamic transitions approach very quickly so that it is very difficult to distinguish them numerically.The second tricritical point is at (m * , a * ) ≈ (0.0001 − 0.002, 100 − 650), marking the crossing between the stability and the UNSAT transition lines.For m ≤ m * the 1-RSB phase is stable up to the maximum average a max (m).This second tricritical point is hard to pinpoint numerically accurately.In the SM, we show analytically that at least for m → 0 the 1-RSB phase is indeed stable up to a max (m).Thus, by continuity, this second tricritical point must exist.

IV. THE SMALL MAGNETIZATION LIMIT
We now study the phase diagram in the m → 0 limit, corresponding to the 1 ≪ k ≪ N regime.The limit must be taken carefully in order to preserve the extensivity of the energy function in the thermodynamic limit.Indeed, we have that for fixed m and N where # denotes the logarithm of the number of microstates at fixed m and N .As m → 0, the energy must be rescaled by c(m) = log(1/m)/m, and the entropy and complexity must be rescaled by m log(1/m).This can be achieved by considering the m → 0 limit of (4) at fixed b = β/c(m).We perform analytically the limit in the RS and 1-RSB solutions, leading to the following phase diagram.
RS phase -For sub-matrix average a < a dynamic = √ 2 c(m), we observe a stable RS phase with zero com-plexity and positive total entropy.In this phase, q 0 = q 1 = m 2 .
Frozen 1-RSB phase -For sub-matrix average a dynamic < a < a static = 2 c(m), we observe a stable 1-RSB phase with q 1 = m, q 0 = m 2 and complexity Σ = 1 − b 2 /4.This is a frozen phase, meaning that each of the exponentially-many pure states contributing to the measure has zero internal entropy.
UNSAT phase -For sub-matrix average a > a static , the total entropy is negative, signalling the onset of the UNSAT phase.
The threshold a dynamic and a static coincide with the thresholds proved in [3] for, respectively, the submatrixaverage of the local maxima A LAS = a dynamic / √ N and the maximum submatrix-average achievable A opt = a static / √ N .Thus, as a byproduct of our analysis, we obtain an equilibrium interpretation of A LAS as a freezing transition.
The m → 0 limit of the MAS phase diagram resembles closely that of the Random Energy Model (REM) [12].More precisely, we find that the static threshold, as well as the values of the entropy and complexity, do coincide (in the REM the static transition threshold equals a static,REM = 2, the complexity equals Σ = 1 − b 2 /4 and the internal entropy is zero for all b > 0).This connection is related to the fact that the MAS energy E(σ) is a Gaussian random variable with covariance ⟨E(σ)E(σ ′ )⟩ ∝ q(σ, σ ′ ) ≤ m, which vanishes in the m → 0 limit.The notable difference between the REM and the MAS problem is given by the finite value of the dynamic threshold in the MAS problem, while in the REM the system is frozen at all temperatures.
In [4], the authors introduced an algorithm called IGP, and they proved that it can find submatrices with average A which, following our analysis, are inside the frozen 1-RSB region.This is at odds with the common belief that solutions in frozen states are algorithmically hard to find [15].Another problem in which a similar situation happens is the binary perceptron, where the algorithmic feasibility was explained in relation to out-of-equilibrium dense regions [16,17].We leave for future work a deeper understanding of the out-of-equilibrium properties of the MAS problem, and more generally the relation between them and algorithmic tractability.

V. CONCLUSIONS
In this paper, we studied the maximum average submatrix problem using tools from the statistical physics of disordered systems, and in particular a mapping onto a variant of the SK model.
We unveiled the phase diagram in the large submatrix regime k = mN , discovering a rich phenomenology including glassy phases, and phases where exponentiallymany pure states contribute to the equilibrium behavior of the system.
By considering the m → 0 limit, we characterized the phase diagram in the small submatrix regime k ≪ N , shedding some light on previous results [3,4] and highlighting a connection to the Random Energy model.We note that there exist efficient algorithms that work in the frozen 1-RSB phase, usually associated with hardalgorithmic phases, similar to what happens in the binary perceptron due to non-equilibrium phenomena.
Our findings leave many questions to be answered, such as i) the study of the out-of-equilibrium properties of the problem and their relation with algorithmic hardness and ii) how for k ≪ N the vanishingly small correlations between the energies combine to shift the dynamical temperature from infinity (REM) to finite (MAS).
We conclude by remarking that our techniques generalise straightforwardly to the case in which the entries of J are non-Gaussian as long as they are i.i.d. with finite first and second moment, and to the rectangular MAS problem, in which both J and the submatrices may be rectangular possibly with different aspect ratios.

Statistical mechanics of the maximum-average submatrix problem Supplemental Materials
The code used to produce the plots can be found at https://github.com/SPOC-group/Maximum-Average-Submatrix.git.

I. MAPPING FROM BOOLEAN TO BINARY SPINS
The energy of the system is defined as where σ is a configuration of N boolean spins σ ∈ {0, 1} N .The straight-forward mapping from boolean to binary spins, i.e. defining s i = 2σ i − 1, leads to the equivalent energy function Thus, written as a function of binary spin configurations s, the energy has the same SK-like interaction term, but develops an additional random magnetic field interaction (the second term).Notice that the random field is not independent from the pair-wise interactions J, leading to a model which is quantitatively different from the usually studied SK model.

II. UNIVERSALITY OF THE PHASE TRANSITIONS
In Figure 1 we presented the phase diagram for the MAS problem for finite m and small m, and observed a variety of phase transitions.In the following table we provide references to a selection of other mean-field models in which the same type of phase transitions arise.We order then in a way that the upper lines are analogous to large m and lower lines to small m in the MAS problem.The phenomenology at each transition of the MAS problem is qualitatively similar to the phenomenology of the corresponding transition in the models listed below.
We notice that a tricritical point similar to the one we found at m = m c ("Tricritical 2" in Figure 1) is present in the perceptron with negative margin.We could not find an example of a model in which a tricritical point akin to the one we find at m = m * ("Tricritical 1" in Figure 1) arises.

III. THE 1-RSB SOLUTION
a. Variational free entropy -To derive the variational free entropy (4), one can follow the derivation for the SK model presented in [23].The only difference in our case is that σ 2 = σ as σ ∈ {0, 1}, contrary to the usual case where logistic( We highlighted that the extremization condition inn p is equivalent to imposing zero complexity.Their derivation requires the usage of integration by parts repeatedly, in the form b. Observables -To derive the expressions for the energy e (5) and total entropy s (6) in the 1-RSB solution, use the grand-canonical thermodynamic relations e = ∂ β Φ−hm and Φ(β) = s+βe+βhm.Applying these relations to the variational free entropy (4) produces the equations presented in the text, to be evaluated then at the equilibrium values of the order parameters.
c. Complexity -To derive the expression for the complexity Σ ( 7) in the 1-RSB solution, we follow [22].We start by introducing a deformed partition function where p is just a weighting parameter -at p = 1 we are computing the usual partition function Z -and the sum over α runs over all pure states contributing to the Gibbs measure, each with free entropy ϕ α .One can interpret p as being a number of replicas of our system that are constrained to be in the same pure state, and Z(p) would then be the correct partition function for such replicated system.Thus, we have E log Z(p)/N ∼ pΦ 1−RSB (p), where the factor p matches the number of replicas used on the two sides of the equation.Now, the extremization over ϕ on the left-hand side can be performed explicitly, and the extremizer ϕ * (p) satisfies Plugging p = 1, corresponding to the equilibrium partition function, we get d. Stability -The most general stability condition for the k-RSB solution can be derived by computing the Hessian at the k-RSB equilibrium of the n-replicas free-entropy.A detailed derivation for the RS solution (q 0 = q 1 ) of the SK model can be found in [23,Chapter 3].By adapting it to our model, it is easy but tedious to see that the the RS stability condition where H RS (u) = √ qz + h + β 2 (m − q).For the 1-RSB stability, we check a simpler condition, called Type-I stability.It is the linear stability of the 2-RSB extremization conditions (seen as fixed point iterations for the overlaps) around the 1-RSB fixed point under the perturbation q 2−RSB 2 = q 1−RSB 1 + ϵ and q 2−RSB 0,1 = q 1−RSB 0,1 .We conjecture that this threshold is not just a bound, but the correct one, in analogy with what happens in the SK model.In practice, we consider the equation for q 2 in the 2-RSB solution, i.e. where q 1 z, and we evaluate it at the 1-RSB equilibrium value of the order parameters, i.e. q 2−RSB 0,1 . Finally, we expand at first order in ϵ.We report the details of the computation in Section VI.The linear stability threshold is the point at which the O(ϵ) of the r.h.s.equals 1.This corresponds to the condition (here p 2−RSB 2 = 1 gives the strictest condition), giving (8).

IV. NUMERICAL SOLUTION FOR THE EQUILIBRIUM ORDER PARAMETERS
The numerical solution for the equilibrium order parameters is performed by solving the saddle-point equations where logistic(x) = (1 + e −x ) −1 and H(u, v) = h + β 2 (m − q 1 ) + √ q 0 u + √ q 1 − q 0 v.There are several aspects that must be explained.
a. What to do with the Parisi parameter p -The complexity prescribes that at equilibrium p = 1 if the associated complexity is positive, and otherwise p = p * such that Σ(p * ) = 0. We solve the first three equations of (S12) for m, q 0 , q 1 at p = 1 first and compute the associated complexity.If it is negative, we discard the solution and solve also for Σ(p) = 0.
b. How to solve (S12) at fixed p -To solve the first three equations of (S12) for m, q 0 , q 1 at fixed p, we could turn them in a fixed-point iteration scheme, as it is done for the SK model.An added element of complexity is that we actually want to fix m, and solve for h.This cannot be done a posteriori as the m(h) function at the stable fixed point of the equation is not single-valued (while the inverse m(h) is).Thus, we iterate the equations for q 0 and q 1 as fixed-point equations, and after each iteration we solve the first equation at fixed m for h -by bisection for example -using the current value of the other order parameters.Whenever no solution for h is found, we perturb it with random Gaussian noise of small variance, and reiterate the procedure.Convergence is declared when the relative change in Euclidean norm of the order parameters is below a set tolerance, in our case 10 −6 .
c. How to solve (S12) imposing Σ(p) = 0 -We use the same procedure as that described above, but every 3 iterations we also solve Σ(p) = 0 for p, using the current value of the other order parameters.
d. How to stably compute the integrals in (S12) -To compute the integrals in Du and Dv in the saddle-point equations we use Gauss-Hermite integration, which approximates Gaussian expectations as for a specific set of weights w i and base-points z i .We use n = 71, and pre-compute the weights and base-points in advance.
Notice that this integrals can lead to overflows and underflows due to exponential factors.For example, consider the simple RS case q 0 = q 1 .In this case, all integrals in the v variable drop, as there is no dependence on v anymore, leading to a cancellation of the factors (1 + exp(βH)) p .If this cancellation has to happen numerically, it can lead to the aforementioned numerical issues.To solve the equations in a numerically stable way, we work in log-space.We write the integrals as

Du
where stable logsumexp functions are routinely implemented in many programming languages.The idea of the logsumexp trick is to write where x * = max i {x i }, so that all exponentials have negative argument.This allows to compute the integrals quickly and stably.e.How to enter the large inverse temperature 1-RSB phase and the dynamical 1-RSB phase In both the large inverse temperature 1-RSB phase and the dynamical 1-RSB phase, the solution of the 1-RSB equations depends highly on the initialization.In the first case, poor initialization leads to non-convergence, while the the second case it leads to convergence to the RS fixed point q 0 = q 1 , in which the complexity erroneously vanishes.To avoid these problems, we manually find good initializations in the full-RSB region or static 1-RSB respectively, at an inverse temperature just above the relevant phase boundary.We then change slightly β (increasing or decreasing it as needed) and feed the previous solution of the equations as initialization at the new inverse temperature.This allows to go at larger values of β in the first case, and inside the dynamical 1-RSB region in the second case.
f. Extrapolation of the observables to compute the stability and the SAT/UNSAT thresholds In order to compute the stability and SAT/UNSAT thresholds we need to extrapolate the values of the energy, entropy and the stability condition to large values of β, larger than reachable with our numerical methods.All these quantities behave as const + O(β) corrections in the physically relevant region (deep in the unstable phase this is not necessarily true anymore for the stability condition).Thus, we fit the large inverse temperature tails of said observables to the const + O(β) asymptotics, and estimate the stability and SAT/UNSAT thresholds on this extrapolated data as the points at which, respectively, the stability condition is not satisfied and the entropy becomes negative.
g. Estimation of the dynamic and static transitions We estimate the dynamic and static transition by considering the complexity (at p = 1), and computing i) where it first develops a discontinuity jumping from zero to positive as a function of β, and ii) where it firsts continuously evolves from positive to negative.To properly threshold for these events, it is important to rescale the complexity with its natural scaling m log(1/m) (see the main-text discussion for a justification of this scaling).
h. Phase diagram as a function of the inverse temperature For completeness, in Figure S1 we present the phase diagram as a function of m and β.

V. THE SMALL MAGNETIZATION LIMIT FOR THE SYMMETRIC MODEL
In this section we derive the m → 0 behaviour of the symmetric MAS problem.We do this by considering a joint of m → 0 and β → ∞, with a scaling relation tuned to preserve the extensivity of the Hamiltonian in the limit.In particular, we have that for fixed m and N where # denotes the logarithm of the number of microstates at fixed m and N .As m → 0, the energy must be rescaled by c(m) = log(1/m)/m, and the entropy and complexity must be rescaled by m log(1/m).This can be achieved by considering the m → 0 limit of (4) at fixed b = β/c(m).Below, we perform analytically the limit in the RS and 1-RSB solutions, leading to the following phase diagram.Notice that our m → 0 results describe the regime 1 ≪ k ≪ N .Indeed, for k = O(1) one must be careful, as the diagonal terms of the submatrix average i J ii σ i (which we discard as subleading in the thermodynamic limit in our analysis) become important.

A. The scaling limit of the RS solution
We start by considering the RS equations in the scaling limit Then We notice here that if q = cm 1+α , with c > 0 and α > 0 due to 0 ≤ q ≤ m, the second term goes to zero and can be safely neglected.The dependence on u drops, and the saddle point equations imply self-consistently The first equation reads from which The energy density equals (we consider the symmetric model here, but for the non-symmetric model the computations are analogous) Given the a posteriori knowledge that the small temperature phase is a 1-RSB phase, we can also access the values of the observables in the frozen 1-RSB phase by • computing b stability as the valued of b at which the entropy density becomes negative; • computing the observables at b stability , as in the 1-RSB phase they do not depend on the temperature anymore.Thus, we have (recall q ≪ m) and the entropy density Notice that, modulo a factor log 2 due a different normalisation of the entropy, the RS m → 0 limit reproduces exactly the results for the RS solution of the Random Energy Model [12].where and Σ(p) is the complexity.The 1-RSB free entropy potential is given by We consider the scaling limit where at leading order.We also assume that µ > 0, in accordance to the RS solution.
We have that and we see that the only term that is going to zero is the last term.Thus, the dependence on u can be dropped, confirming immediately through the second SP equation that q 0 = m 2 .The solution of the equation is detailed in the next section V C.Here we summarise the results.If c 2 g → 0 -happening whenever c goes to 0 polynomially -then the dependence on v drops, and we get q 1 = m 2 similarly as it happened for q 0 .This retrieves the RS solution computed in the previous Section.
If c = O m (1), then we have g = − log m to all non-vanishing orders in m, and under the condition b ≥ max( √ 2/p, 2/p).The equilibrium solution is given by the condition p = 1 whenever Σ(b, p = 1) > 0, and by the value of p such that Σ(b, p) = 0 otherwise.This gives p = 1 for √ 2 < b < 2, and p = 2/b for b > 2.
• In the region of positive complexity p = 1, we have that exponentially many thermodynamic states with the same free entropy contribute to the thermodynamics of the system.We have that so that confirming that this phase is a frozen 1-RSB phase, and In black, we plot the analytical predictions for the m → 0 limit.We see that both observables converge to the m → 0 limit extremely slowly, in accordance with our analytical analysis that shows that the next-to-leading order in the m ≪ 1 expansion is only logarithmically small, see SM.Moreover, we see that the internal entropy in the dynamical 1-RSB region (where it is different from the total entropy due to the non-vanishing complexity) decreases as m goes to zero, foreshadowing the frozen 1-RSB phase that arises in the limit m → 0.
so that confirming that for b > 2 we have that a = a max = a stability .
We find another 1-RSB solution for 0 < p < 1 that is not stable under 2-RSB perturbations, and that could not be reproduced by us by solving numerically the 1-RSB saddle-point equations, so we discard it as unphysical.
C. Solution of the 1-RSB SP equations in the scaling limit We start by noticing that for c 2 → 0 such that c 2 g → 0 the dependence on v drops, and we get q 1 = m 2 similarly as it happened for q 0 .This retrieves the RS solution.If instead the v dependence does not drop, then we can expect q 1 = c 2 m with c ̸ = 0. Thus, our task is now to solve the remaining three SP equations for µ, c, p at fixed b and for c ̸ = 0, µ > 0.
Notice that, calling we have that the 1-RSB equations involve integrals which are in the setting of Proposition 1, with p > 0 and ℓ = 0, 1, 2.
The condition B 1 ̸ = 0 is equivalent to c ̸ = 0. Using Proposition 1, we then have that and and the equations to be solved are and recall We consider the equations in the three cases 0 < p < 1, 1 ≤ p < 2 and p ≥ 2 separately: 1. Case 0 < p < 1.We have the subcases: leading to the following equations (notice that if −A 2 > 0 than G(−A/B) = 1): 1a) Here all integrals are equal at leading order, so there exists no solution to the SP equations which is compatible with the scaling ansatz.1b) Same as 1a).1c) The equations give c 2 = 1 and e f (p) G A B + pB = m.We have The second equation implies at leading logarithmic order For this to be satisfied, we need to impose for some γ > 0 with ϵ > 0 small, x ∈ (p 1 , 1) and m, q 1 , q 0 , p 1 are the solutions to the 1-RSB SP equations.
We start by expanding The replicated partition function reads (remember that σ 2 = σ as σ = 0, 1 and similar for τ )

FIG. 1 .
FIG.1.The phase diagram of the MAS problem as a function of the submatrix-average a and the submatrix size m = k/N for linear scale in m (left), logarithmic scale in m (center) and m → 0 (right).In the central and right panel we rescale the sub-matrix average as a/ log(1/m)/(m) to highlight the convergence to the limit.We identify five distinct phases.In the RS phase (green) the system is replica symmetric.In the 1-RSB phases replica symmetry is broken to one step and two sub-phases exists, a dynamical 1-RSB with an extensive number of equilibrium pure states (blue) and a static 1-RSB with only finitely many pure state (purple).All the phase boundaries are exact in the thermodynamic limit except the boundary between full-RSB (orange) and UNSAT (red) which would require solving the full-RSB equations.In the full-RSB phase (orange) replica symmetry is completely broken and the set of pure states manifests ultrametricity.In the unsatisfiable phase (UNSAT, red) no submatrix exists with the given values of a and m.The transition from RS and 1-RSB to full-RSB is continuous and caused by an instability of the 1-RSB ansatz (dashed line), while the other transitions are discontinuous.We observe two tricritical points, one at (mc, ac) where the system shows coexistence of RS, 1-RSB and full-RSB phase (white marker), and one at (m * , a * ) where the largest-average submatrices become 1-RSB stable and the full-RSB region ceases to exist (black marker).In the limit m → 0, we observe only the RS, 1-RSB and UNSAT phases.The 1-RSB phase is frozen, meaning that the internal entropy of each pure state goes to zero in the m → 0 limit.

FIG. S1 .
FIG. S1.The phase diagram of the MAS problem as a function of the inverse temperature β and the submatrix size m = k/N for linear scale in m (left) nad logarithmic scale in m (right).On the right we rescaled the inverse temperature to highlight the convergence to the m → 0 limit.

)•
FIG. S2.Behavior of the internal entropy (left panel, solid lines), total entropy (left panel, dashed lines) and the submatrix average (right panel) as a function of the rescaled inverse temperature β/ log(1/m)/m for various values of m.In black, we plot the analytical predictions for the m → 0 limit.We see that both observables converge to the m → 0 limit extremely slowly, in accordance with our analytical analysis that shows that the next-to-leading order in the m ≪ 1 expansion is only logarithmically small, see SM.Moreover, we see that the internal entropy in the dynamical 1-RSB region (where it is different from the total entropy due to the non-vanishing complexity) decreases as m goes to zero, foreshadowing the frozen 1-RSB phase that arises in the limit m → 0.

dm r a dm c a a<b dq r ab dq c ab × exp N r βh r a m r a + N c βh c a m c a + N r N c β 2 2 a m r a m c a + N r N c β 2 a<b
the order parameters using delta functions and their exponential representationE J Z n = a