Turing in the Shadows of Nobel and Abel: An Algorithmic Story Behind Two Recent Prizes

The 2021 Nobel Prize in physics was awarded to Giorgio Parisi “for the discovery of the interplay of disorder and fluctuations in physical systems from atomic to planetary scales,” and the 2024 Abel Prize in mathematics was awarded to Michel Talagrand “for his groundbreaking contributions to probability theory and functional analysis, with outstanding applications in mathematical physics and statistics.” What remains largely absent in the popular descriptions of these prizes, however, is the profound contributions the works of both individuals have had to the field of algorithms and computation. The ideas first developed by Parisi and his collaborators relying on remarkably precise physics intuition, and later confirmed by Talagrand and others by no less remarkable mathematical techniques, have revolutionized the way we think algorithmically about optimization problems involving randomness. This is true both in terms of the existence of fast algorithms for some optimization problems, but also in terms of our persistent failures of finding such algorithms for some other optimization problems.

The goal of this article is to highlight these developments and explain how the ideas pioneered by Parisi and Talagrand have led to a remarkably precise characterization of which optimization problems admit fast algorithms, versus those which do not, and furthermore to explain why this characterization holds true.

1. Optimization of Random Structures

The works of Parisi and Talagrand—which were devoted to understanding a mysterious object in statistical physics known as spin glass—offered an entirely novel way of understanding algorithmic successes and failures in tackling optimization problems involving randomness. This progress was propelled by understanding the “physics” properties of the underlying problems, namely the geometry of the solution space. We will illustrate these notions using three examples, all of which fit into a general framework of optimization problems of the form

$$\begin{align} \eta _{\mathrm{OPT}} \triangleq \max _{\sigma \in \Theta } H(\sigma ,\mathcal{R}). {\tag{1}} \end{align}$$

Here $H$ is the function to be maximized, which we call the “Hamiltonian” for reasons to be explained later, $\Theta$ is the space of potential solutions, and $\mathcal{R}$ is a random variable describing the source of randomness. Our three examples all exhibit a persistent computation-to-optimization gap, namely a gap between the optimal value $\eta _{\mathrm{OPT}}$ and the best guaranteed value that can be computed by known efficient algorithms—where “efficient” means “polynomial-time” in the parlance of the theory of computational complexity.

Example I. Largest cliques and independent sets in random graphs

The first of our three examples is the problem of finding a largest clique in a dense Erdös–Rényi graph and the related problem of finding a largest independent set in a sparse Erdös–Rényi graph.

Imagine a club with $n$ members where every pair of members are friends or nonfriends with equal likelihood $1/2$, independently across all ${n\choose 2}$ member pairs. What is the size of the largest group of members who are all friends with each other? This quantity is known as the clique number of the graph whose nodes are the members and whose edges correspond to the pairs of friends. This club of random friend pairings is just a vernacular description of a model of a random graph known as an Erdös–Rényi graph, denoted $\mathbb{G}(n,1/2)$. Paul Erdös introduced this model in 1948 as a very effective tool to study and solve problems in the field of extremal combinatorics. The question of estimating the largest clique size of an Erdös–Rényi graph was solved successfully in the 1970s and has a remarkably precise answer: With high probability (namely with probability approaching $1$ as $n\to \infty$, denoted as “whp”), the size is $2(1+o(1))\log _2 n$. Here $o(1)$ denotes a term which converges to zero as $n\to \infty$. This is an instance of our optimization problem 1 in which $\Theta$ is the set of all subsets of the node set $[n]\coloneq \{1,\ldots ,n\}$ of the graph; $\mathcal{R}$ is the random graph itself; and $H(\sigma ,\mathcal{R})$ is the cardinality of $\sigma$ when $\sigma$ is a clique of $\mathcal{R}$, and $H=0$ otherwise.

The asymptotic value $2(1+o(1))\log _2 n$ has a very simple intuition. Fix an integer value $x$, and let $Z(x)$ denote a random variable corresponding to the number of subsets of $[n]$ which have cardinality $x$ and are cliques. By the linearity of expectation, the expected value of this random variable is ${n\choose x}\left({1\over 2}\right)^{x\choose 2}$. By straightforward computation, this quantity converges to zero when $x$ is asymptotically larger than $2\log _2 n$, which provides a simple proof of the upper bound. The matching lower bound requires more work involving estimation of the variance of $Z(x)$ which is now standard material in many textbooks on random graphs.

What about the problem of efficiently finding the largest clique algorithmically in at most polynomial time, when the graph is given to the algorithm designer as an input? Here and elsewhere we adhere to the gold standard of algorithmic efficiency when requiring that the running time of the algorithm grows at most polynomially in the size of the input. That is, we require that the algorithm belongs to $\mathcal{P}$—the class of algorithm whose running time grows at most polynomially in the size of the problem description, which in our case means it grows at most polynomially in $n$. Note that the largest clique problem can be trivially solved in time $n^{O(\log n)}$ by checking all subsets with cardinality $2(1+o(1))\log _2 n$; this growth rate is super-polynomial, and thus this exhaustive search algorithm is not in $\mathcal{P}$.

This is where our story begins. In 1976, Richard Karp—one of the founding fathers of the modern algorithmic complexity theory—presented a very simple algorithm that computes a clique whose size is guaranteed to be (roughly) half-optimal, namely $(1+o(1))\log _2 n$ Kar76. Given that his algorithm was not particularly “clever,” Karp challenged the community to produce an algorithm with a better performance guarantee, with the stipulation that the algorithm still remains in $\mathcal{P}$. Remarkably, in the almost-fifty years since the problem was introduced, there has been no successful improvement of Karp’s rather naïve algorithm, including notably the algorithm based on Markov Chain Monte Carlo updates, as shown by Jerrum in 1992. This persistent failure further motivated Jerrum to introduce one of the most widely studied problems in random structures and high-dimensional statistics, namely the ”infamous” Hidden Clique Problem. The problem of finding large cliques in $\mathbb{G}(n,1/2)$ thus exhibits a computation-to-optimization gap.

Notably, while improving the half-optimality for the clique problem appears algorithmically hard, this hardness is only evidential but not (at least not yet) formal. Namely, we do not have any proof of even conditional algorithmic hardness for this problem say by assuming the validity of the “$\mathcal{P}\ne \mathcal{NP}$” conjecture. This is in stark contrast with the worst-case setup. Namely, consider the problem of finding the largest clique for an arbitrary graph, in particular for graphs which are not necessarily generated according to the $\mathbb{G}(n,1/2)$ distribution. For this version of the problem, it is known that no polynomial time algorithm can produce a clique within any multiplicative factor $n^{1-\epsilon }$ from optimality (in particular, much larger than factor $1/2$) unless $\mathcal{P}\ne \mathcal{NP}$. We see that, in particular, the half-optimality gap emerging in the random case does not align with the classical algorithmic complexity theory based on the $\mathcal{P}\ne \mathcal{NP}$ conjecture.

A very similar story takes place when the underlying model is a sparse random graph defined next. Fix a constant $d$. Consider a random graph denoted by $\mathbb{G}(n,d/n)$ in which every pair of nodes is an edge with probability $d/n$, independently across all pairs. We now switch to the mirror image of the clique problem, known as the independent set problem, defined as follows. A subset of nodes $I\subset [n]$ is called an independent set if it contains no edges. Namely for all $i,j\in I$, the pair $(i,j)$ is not an edge. Frieze proved Fri90 that the size of the largest independent set is $\eta _{d,{\mathrm{OPT}}} \triangleq 2(1+o_d(1)){\log d\over d}n$ whp as $n\to \infty$, where here $o_d(1)$ denotes a constant vanishing in $d$ as $d$ increases.

Consider now the problem of identifying the largest independent set algorithmically in polynomial time. Using the formalism described in 1, $\Theta$ is the set of all subsets of the node set $[n]\coloneq \{1,\ldots ,n\}$ of the graph; $\mathcal{R}$ is the random graph $\mathbb{G}(n,d/n)$ itself; and $H(\sigma ,\mathcal{R})$ is the cardinality of $\sigma$ when $\sigma$ is an independent set in the graph $\mathbb{G}(n,d/n)$, and $H=0$ otherwise. An algorithm similar to Karp’s clique algorithm computes an independent set whose size is roughly half-optimal, namely $(1+o_d(1)){\log d\over d}n$, also whp. No improvement of this algorithm is known and no formal hardness is known either. Thus, this problem is also evidentially hard and exhibits a computation-to-optimization gap.

Example II. Ground states of spin glasses

Our second example is in fact a probabilistic model of spin glasses themselves. Consider a tensor of random variables $J=\left(J_{i_1,\ldots ,i_p}, 1\le i_1<i_2<\cdots < i_p\le n\right)\in \mathbb{R}^{n\otimes p}$, where each entry $J_{i_1,\ldots ,i_p}$ of the tensor has a standard normal distribution and all entries are probabilistically independent. The optimization problem of interest is

$$\begin{align} \eta _{n,p}\triangleq \max _{\sigma \in \{-1,1\}^n} n^{-{p+1\over 2}}\langle J,\sigma ^{\otimes p}\rangle , {\tag{2}} \end{align}$$

where for every $\sigma \in \{-1,1\}^n$,

$$\begin{align} \langle J,\sigma ^{\otimes p}\rangle =\sum _{1\le i_1<i_2<\cdots < i_p\le n}J_{i_1,\ldots ,i_p}\sigma _{i_1}\sigma _{i_2}\cdots \sigma _{i_p}. {\tag{3}} \end{align}$$

The optimal solution $\sigma ^*$ of this problem (which is unique up to a sign flip) is called the ground state of the model using physics jargon, and the optimal value is called the ground state value. The normalization $n^{-{p+1\over 2}}$ is chosen for convenience (see below). In our formalism 1, $\Theta =\{-1,1\}^n$, $\mathcal{R}$ is the tensor $J$, and $H(\sigma ,\mathcal{R})$ is the associated inner product 3. The special case $p=2$ corresponds to the well-known and widely studied Sherrington–Kirkpatrick model.

Notice how the signs of $J_{i_1,i_2}$ impact the optimal choices for $\sigma _i$. When $J_{i_1,i_2}>0$, $\sigma _{i_1}$ and $\sigma _{i_2}$ prefer to have the same sign so that the product $J_{i_1,i_2}\sigma _{i_1}\sigma _{i_2}$ is positive. On the other hand, when $J_{i_1,i_2}<0$, $\sigma _{i_1}$ and $\sigma _{i_2}$ prefer to have opposite signs so that the product $J_{i_1,i_2}\sigma _{i_1}\sigma _{i_2}$ is again positive. What is the best arrangement of signs to $\sigma _i$ which optimizes this tradeoff? This is the core of the optimization problem, which turned out to be quite challenging.

In 1980, Parisi solved optimization problem 2 by introducing a mathematically nonrigorous but remarkably powerful physics-style method called replica symmetry breaking (RSB) Par80. The importance of the RSB-based technique is hard to overstate: It is the method of choice for computing the ground state values and related quantities for complex statistical physics models. Parisi’s solution of problem 2 is as follows. Let $\mathcal{U}$ be the set of nonnegative, nondecreasing, and right-continuous functions $\mu$ on $[0,1]$ satisfying $\int _0^1\mu (t)dt<\infty$. Consider the solution $\Psi _\mu (t,x), (t,x)\in [0,1)\times \mathbb{R}$ of the following PDE with boundary condition $\Psi (1,x)=|x|$:

$$\begin{align*} \partial _t \Psi _\mu =-{p(p-1) \over 2} t^{p-2}\left(\partial ^2_x\Psi _\mu +\mu \left(\partial _x \Psi \right)^2\right). \end{align*}$$

Define $\mathcal{P}(\mu )=\Psi _\mu (0,0)-{p(p-1)\over 2}\int _0^1 t^{p-2}\mu (t)dt$.

Theorem 1.1.

For every $p$, the ground state value satisfies

$$\begin{align} \lim _n \eta _{p,n}= \inf _{\mu \in \mathcal{U}}\mathcal{P}(\mu )\triangleq \eta _{p, \mathrm{OPT}}, {\tag{4}} \end{align}$$

whp as $n\to \infty$.

The optimal solution $\mu _{p,\mathrm{OPT}}$ of the variational problem above is known to be unique and defines a cumulative distribution function of a probability measure, which has been called since then the Parisi measure. The identity in 4 was derived by Parisi heuristically using the RSB method. It was turned into a theorem by Talagrand in 2006 in his breakthrough paper Tal06. (Thus one could say that the scientific value of the formula in 4 is worthy of two major prizes!) Talagrand’s approach relied on a powerful interpolation bound developed earlier by Guerra Gue03.

What about algorithms? Here the goal is to find a solution $\sigma$ such that $n^{-{p+1\over 2}}\langle J, \sigma ^{\otimes p}\rangle$ is as close to $\eta _{p, \mathrm{OPT}}$ as possible when the tensor $J$ is given as an input. Interestingly, major progress on this question was achieved only very recently, and it deserves a separate discussion found in Section 5. Jumping ahead though, the summary of the state-of-the-art algorithms is as follows. For $p=2$, near optimum solutions can be found by a polynomial time algorithm. On the other hand, when $p>2$, this is evidently not the case. In summary, the largest value $\eta _{p,\mathrm{ALG}}$ achievable by known poly-time algorithms satisfies $\eta _{2,\mathrm{ALG}}=\eta _{2,\mathrm{ OPT}}$, $\eta _{p,\mathrm{ALG}}<\eta _{p,\mathrm{OPT}}, p>2$. Thus similarly to Example I, the model exhibits a computation-to-optimization gap between the optimal values and the values found by the best-known polynomial time algorithms.

Example III. Satisfiability of a random K-SAT formula

Our third example, the final family of examples, consists of randomly generated constraint satisfaction problems, specifically the so-called random K-SAT problem. An instance of a K-SAT model is described by $n$ Boolean variables $x_1,\ldots ,x_n$ (variables with $0/1$ values) and a conjunction $\Phi \triangleq \wedge _{1\le j\le m} C_j$ of $m$ Boolean functions $C_1,\ldots ,C_m$, each of which is a disjunction of precisely $K$ variables from the list $x_1,\ldots ,x_n$ or their negation. Here is an example for $n=10$ and $k=3$: $\Phi =\left( x_3\vee \neg x_7 \vee \neg x_8\right) \wedge \left( \neg x_1\vee \neg x_2 \vee x_7\right) \wedge \left( \neg x_3\vee \neg x_7 \vee \neg x_9\right)$. An assignment of variables $x_1,\ldots ,x_n$ to values $0$ and $1$ is called a satisfying assignment if $\Phi =1$ under this assignment. For example, the assignment $x_2=0,x_3=1,x_9=0$ (and any assignment for the remaining $x_i$s) is a satisfying assignment. The formula $\Phi$ is called satisfiable if there exists at least one satisfying assignment and is called nonsatisfiable otherwise.

A random K-SAT formula is generated according to the following probabilistic construction. For each $j=1,\ldots ,m$, independently across $j$ we select $K$ variables $x_{j,1},\ldots ,x_{j,K}$ uniformly at random from $x_1,\ldots ,x_n$ and attach the negation sign $\neg$ to each of them with probability $1/2$, also independently across $j=1,2,\ldots ,m$ and $k=1,\ldots ,K$. The resulting formula $\boldsymbol{\Phi }_{n,m}$ is a random formula whose probability distribution depends on $n$ and $m$ only. In the context of optimization problem 1, $\Theta =\{0,1\}^n$, $\mathcal{R}$ is the description of the formula $\boldsymbol{\Phi }_{n,m}$, and $H(\sigma ,\mathcal{R})$ is the number of clauses satisfiable by the assignment $\sigma \in \Theta$. In particular, the formula is satisfiable if and only if $\max _{\sigma }H(\sigma ,\mathcal{R})=m$.

Whether or not $\boldsymbol{\Phi }_{n,m}$ is satisfiable or not is now an event, which we write symbolically as $\boldsymbol{\Phi }_{n,m}\in {\mathrm{SAT}}$. The probability of this event $\mathbb{P}\left(\boldsymbol{\Phi }_{n,m}\in {\mathrm{SAT}}\right)$ is our focus. There exists a (conjectural) critical threshold $c_{K,\mathrm{SAT}}$ such that $\lim _n \mathbb{P}\left(\boldsymbol{\Phi }_{n,cn}\in {\mathrm{SAT}}\right)=1$ when $c<c_{K,\mathrm{SAT}}$ and $=0$ when $c>c_{K,\mathrm{SAT}}$. Informally, $c_{K,\mathrm{SAT}}$ represents the largest fraction of clauses to variables for which the solution exists whp asymptotically in $n$. The existence of such a threshold $c_{K,\mathrm{SAT}}$ was recently proven for large enough $K$ DSS22.

The constant $c_{K,\mathrm{SAT}}$ has a very simple asymptotic expansion for large $K$: $c_{K,\mathrm{SAT}}=2^K\log 2(1+o_K(1))$, where $o_K(1)$ hides a function vanishing to zero as $K\to \infty$. The intuition for this asymptotics is rather straightforward. Observe that any fixed assignment of $x_i, 1\le i\le n$, is satisfiable with probability $(1-1/2^K)^m$. Thus the expected number of satisfying solutions is $2^n(1-1/2^K)^m$. The value $2^K\log 2$ is precisely the critical value of $c$ at which this expectation drops from $\exp (\Theta (n))$ (exponentially large) to $\exp (-\Theta (n))$ (exponentially small).

Turning to algorithms, the best-known one succeeds in finding a satisfying assignment only when $c\le c_0 2^K {\log K \over K}\triangleq c_{K, \mathrm{ALG}}$ CO10. Here $c_0$ denotes a universal constant. Namely, the best-known algorithm succeeds only when the clause to variables density is nearly a factor $K/\log K$ (!) smaller than the density for which the solutions are guaranteed to exist whp. This is an even more dramatic gap than the half-optimality gap seen earlier for cliques/independent sets in the $\mathbb{G}(n,1/2)$ and $\mathbb{G}(n,d/n)$ models. Random K-SAT is thus another instantiation of the computation-to-optimization gap, similar to the ones found in Examples I and II.

In the remainder of the article, we will explain the nature of these gaps. In particular, we will explain the values $1/2, \eta _{p, \mathrm{ALG}}$ and $c_{K, \mathrm{SAT}}$ by leveraging powerful ideas first developed toward understanding statistical physics of materials called spin glasses.

2. Physics Matters: Spin Glasses

Statistical physics studies properties of matter using probabilistic models which provide statistical likelihoods one finds physical systems in different states. These likelihoods in their turn are linked with the energy of these states: As the energy decreases, the likelihood increases. One can apply this view point to our optimization problem $\max _\sigma H(\sigma ,\cdot )$, thinking of $H(\sigma ,\cdot )$ as (negative) “energy” (Hamiltonian) of the state $\sigma$. The ground state is then a state with the lowest energy. The link between energies and likelihoods is formalized by introducing the Gibbs measure on the state space $\Theta$, which assigns probability measure $\langle \sigma \rangle =Z^{-1}\exp (\beta H(\sigma ,\cdot ))$ to each configuration $\sigma$. Here $Z$ is the normalizing constant called the partition function. The parameter $\beta$ is called inverse temperature. At low temperature, namely large values of $\beta$, the Gibbs measure becomes concentrated in near ground states: states with nearly the lowest energy. The identity in 4 was first derived for the case of finite $\beta$ and then extended to the case $\beta =\infty$, corresponding to the ground state value. The added value of introducing the Gibbs measure is that it allows one to consider the entire space $\Theta$ of solutions in one shot where the solutions are stratified by their energy levels. This is the so-called “solution space geometry” approach, which turned out to be very fruitful for the analysis of algorithms.

Spin glasses are particularly complex examples of matters which exhibit properties of being liquid and solid at the same time. For a recent treatment of the subject, see CMM$^{+}$23. Mathematically, spin glasses correspond to the fact that the associated Hamiltonian $H$, which is a sum $\sum _j H_j$ of many local Hamiltonians representing interactions of nearby atoms, is “frustrated.” Put differently, it is impossible to find a spin configuration $\sigma$ which simultaneously minimizes the energy of each component $H_j$. This makes it very similar to problems we have considered in our examples.

The idea of repurposing the techniques from the theory of spin glasses for the field of optimization of combinatorial structures goes back to the mid-1980s. In particular, Fu and Anderson in their 1986 paper “Application of statistical mechanics to NP-complete problems in combinatorial optimisation” say, “…we propose to study the general properties of such problems (combinatorial optimization) in the light of …the structure of solution space, the existence and the nature of phase transitions ….”

3. Clustering Phase Transition in Random K-SAT

The first fruitful and mathematically verified link between the solution space geometry and phase transition on the one hand and the algorithmic tractability/hardness on the other hand appeared in concurrent works of Achlioptas and Ricci-Tersenghi ART06 and Mézard, Mora, and Zecchina MMZ05 in the context of the random K-SAT model introduced earlier. Given a random formula $\boldsymbol{\Phi }_{n,m}$, let $\mathcal{SAT}_{n,m}\subset \{0,1\}^n$ denote the (random) set of satisfying assignments. RSB calculations, which were originally introduced to study spin glasses, suggest heuristically that the bulk of this space should be highly clustered at some large enough clause density $c<c_{K,\mathrm{SAT}}$. This is what these papers established in a mathematically rigorous way. In particular, they proved the following theorem.

Theorem 3.1.

For every $K\ge 8$, there exists $c_{K, \mathrm{Clust}}<c_{K,\mathrm{SAT}}$ such that the following holds when $c\in (c_{K, \mathrm{Clust}},c_{K,\mathrm{SAT}})$: whp as $n\to \infty$ there exists a disoint partition $\mathcal{SAT}_{n,m}=\left(\cup _{\ell } \mathcal{SAT}_\ell \right) \cup \mathcal{SAT}_o, \ell \ge 2$, such that the Hamming distance between any two distinct sets $\mathcal{SAT}_{\ell _1}$ and $\mathcal{SAT}_{\ell _2}, \ell _1\ne \ell _2$ is $\Omega (n)$, and $|\mathcal{SAT}_o |/ |\mathcal{SAT}_{n,m}|=\exp (-\Omega (n))$.

Figure 1.

Clustering phase transition at $\alpha _{\mathrm{Clust}} \approx \alpha _{\mathrm{ALG}}$. Below $\alpha _{\mathrm{Clust}}$, a bulk of the solution space (blue) is a giant connected set. Above $\alpha _{\mathrm{Clust}}$, a bulk of the set is partitioned into (blue) clusters, with exponentially small exceptions (gray).

This is illustrated in Figure 1. The authors argued that the presence of such a clustering picture suggests algorithmic hardness of finding satisfying assignments, in particular by means of some local search type algorithms. Furthermore, as subsequent works showed, intriguingly, the onset of this clustering property phase transition coincides with the evidential algorithmic hardness discussed earlier, asymptotically as $K$. More precisely, it was verified that $c_{K, \mathrm{Clust}}=(1+o_K(1))2^K{\log K\over K}=(1+o_K(1))c_{K,\mathrm{ALG}}$! This was the first (though purely hypothetical) link between the solution space geometry and the algorithmic hardness.

Very quickly a similar result was established also for the independent set problem in $\mathbb{G}(n,d/n)$ and several other similar problems: For all $c\in (1,2)$, independent sets of size $c(1+o_d(1){\log d\over d}n$ (namely more than half optimum) are clustered, and for all $c\in (0,1)$, they are not clustered.

What about the clustering property in spin glasses? It does take place, and this is in fact what originally prompted Parisi to develop his RSB approach. Clustering in spin glasses is more intricate and comes in the form of a nested hierarchy of subclusters arranged in a particular way called an ultrametric tree. We delay the discussion of this structure until Section 5 when we discuss how this ultrametric structure of clusters successfully guided the creation of the state-of-the-art algorithm for problem 2.

4. Overlap Gap Property: A Provable Barrier to Algorithms

As discussed earlier, the clustering property exhibited by random K-SAT and similar other problems provided unfortunately only a hypothetical explanation of the algorithmic hardness, but not a proof of it. The property, which turned out to be a provable barrier to algorithms, was inspired by the clustering picture but came in the different form of the overlap gap property (OGP), which we now formally introduce in the context of problem 1. For this purpose, we assume that the solution space $\Theta$ is equipped with some metric $d$ (for example, Hamming distance when $\Theta$ is binary cube).

Definition 4.1.

An instance $H(\cdot ,\mathcal{R})$ of optimization problem 1 exhibits an OGP with parameters $\eta _{\mathrm{OGP}}\le \eta _{\mathrm{OPT}}, \nu _1<\nu _2$ if for every $\sigma ,\tau \in \Theta$ satisfying $H(\sigma ,\mathcal{R}),H(\tau ,\mathcal{R})\ge \eta _{\mathrm{OGP}}$ it holds that $d(\sigma ,\tau )\notin (\nu _1,\nu _2)$.

This property has a simple interpretation. The model exhibits the OGP if the distance between every pair of “near”-optimal solutions (solutions with values at least $\eta _{\mathrm{OGP}})$ is either “small” (has a high overlap represented by $d(\sigma ,\tau )\le \nu _1$) or is “large” (has a low overlap represented by $d(\sigma ,\tau )\ge \nu _2$).

We need to extend Definition 4.1 in two ways. First involving ensembles of instances $\mathcal{R}_1,\ldots ,\mathcal{R}_m$, and second involving overlaps of more than two solutions.

Definition 4.2.

Fix a positive integer $M$ and a subset $\boldsymbol{\nu }\subset \Theta ^M$. A family of problem instances $\mathcal{R}_1,\ldots ,\mathcal{R}_T$ exhibits an OGP with parameters $(\eta _{\mathrm{OGP}}, \boldsymbol{\nu })$ if for every $1\le \ell _1\le \ell _2\le \cdots \le \ell _M\le T$ and every $\sigma _1,\ldots ,\sigma _M \in \Theta$ satisfying $H(\sigma _1, \mathcal{R}_{\ell _1}), \ldots , H(\sigma _M, \mathcal{R}_{\ell _M})\ge \eta _{\mathrm{OGP}}$ it holds that $\left(\sigma _1,\ldots ,\sigma _M\right)\notin \boldsymbol{\nu }$.

Definition 4.1 is recovered as a special case when $T=1, M=2$, and

$$\begin{align} \Theta ^2 \supset \boldsymbol{\nu }=\{(\sigma ,\tau ): d(\sigma ,\tau )\in (\nu _1,\nu _2)\}. {\tag{5}} \end{align}$$

We stress that the randomness $\mathcal{R}$ appearing in $H(\sigma _j, \mathcal{R}_{\ell _j})$ is index-specific, representing the fact that $\sigma _j$ is near-optimal (achieves a value of at least $\eta _{\mathrm{OGP}})$ specifically for instance $\mathcal{R}_{\ell _j}$. While the intuition behind Definition 4.1 is fairly direct, the one for Definition 4.2 is rather opaque and requires some elaboration, to be done shortly.

It turns out that the OGP does hold in our models and is a provable obstruction to classes of algorithms. The first of these facts (in the context of Definition 4.1) was already implicitly established in ART06MMZ05 as a method for proving the clustering property for the random K-SAT model. The fact though that the OGP is a provable obstruction to algorithms was established first in GS17 in the context of independent sets in sparse Erdös–Rényi graphs searchable by local algorithms. Loosely speaking, an algorithm is local if for each node $i\in [n]$ the decision to include it into an independent set is conducted entirely based on a constant size graph-theoretic neighborhood around $i$, simultaneously for all $i$, ignoring the rest of the graph. The result in GS17 says the following.

Theorem 4.3.

For every $\epsilon >0$ and sufficiently large $d$, there exist $0< \nu _1<\nu _2<1$ such that whp as $n\to \infty$ the set of independent sets in $\mathbb{G}(n,d/n)$ exhibits an OGP (in the sense of Definition 4.1) with parameters $\left(\eta _{\mathrm{OGP}}\triangleq \left({1\over 2}+{1\over 2\sqrt {2}}+\epsilon \right)\eta _{d,\mathrm{OPT}}; \nu _1n, \nu _2n \right)$. As a result, no local algorithm (appropriately defined) exists for constructing an independent set with size larger than $\eta _{\mathrm{OGP}}$.

Namely, every two independent sets with cardinality which is multiplicative factor ${1\over 2}+{1\over 2\sqrt {2}}$ close to optimality are either at a “small” distance from each other (at most $\nu _1 n$) or at a “large” distance from each other (at least $\nu _2 n$). Theorem 4.3 means that beating a factor ${1\over 2}+{1\over 2\sqrt {2}}$ optimality threshold for this problem will require going beyond the class of local algorithms. Local algorithms do in fact achieve half-optimality (which we recall is the evidential barrier for this problem).

The gap between half-optimality and ${1\over 2}+{1\over 2\sqrt {2}}$ was eliminated in RV17 for the same class of local algorithms and later for a larger class of algorithms (see Theorem 4.4). This was achieved by turning to multi-overlaps of solutions and resorting to Definition 4.2. A method to address broader classes of algorithms requires interpolating between instances, that is having $T\ge 2$ in Definition 4.2, as was introduced in CGPR19.

The broadest currently known class of algorithms ruled out by the OGP is the class of stable (noise insensitive) algorithms. We will not give one formal definition of this class, as it comes in different context-dependent variants. Loosely speaking, the stability property means the following. We think of an algorithm as a map $\mathcal{ALG}: \mathcal{R}\to \Theta$, which maps instances $\mathcal{R}$ of random input (graphs, tensors, etc.) to solutions. We say that the algorithm $\mathcal{ALG}$ is stable if a small perturbation of the input $\mathcal{R}$ results in a small perturbation of the solution $\mathcal{ALG}(\mathcal{R})$. Namely, if $\mathcal{R}\approx \mathcal{R}'$, then $\mathcal{ALG}(\mathcal{R})\approx \mathcal{ALG}(\mathcal{R}')$. A local algorithm is an example of a stable algorithm: deleting or adding one edge of a graph will affect the decisions associated only with nodes which contain this edge in their constant size neighborhood. This is so because the algorithm performs these decisions concurrently (in parallel). Thus deleting or adding one edge changes the independent set produced by an algorithm in at most $O(1)$ many nodes. Other classes of algorithms proven to be successful for these types of problems, including approximate message passing (see Section 5), stochastic localizations, low-degree polynomials, and Boolean circuits with small depth, turn out to be stable as well.

It was established in a series of followup works that the OGP is a barrier to stable algorithms GJW24GMZ22. For our problems outlined as Examples I and III, it provides tight algorithmic bounds as described in the following theorem, stated informally. The sequence of interpolated instances involving these examples is constructed as follows. For our Example I, involving sparse random graph $\mathbb{G}(n,d/n)$, consider any ordering $\mathcal{P}$ of pairs $(i,j), 1\le i<j\le n$, and consider the process of resampling each edge according to the same Bernoulli distribution with values $1,0$ with probabilities $d/n,1-d/n$, respectively, one pair at a time according to the order $\mathcal{P}$. At each step $t$, we obtain a random graph $\mathbb{G}_t$, which is distributionally identical to $\mathbb{G}(n,d/n)$. After $T={n\choose 2}$ steps, we obtain a fresh copy of $\mathbb{G}(n,d/n)$, which is probabilistically independent from the starting copy. This is our interpolated sequence $\mathcal{R}_1,\ldots ,\mathcal{R}_T$ in the context of Example I. It is constructed similarly for the case of random K-SAT (Example III).

Theorem 4.4.

For every large enough $d$, there exist a large enough $M$ and a subset $\boldsymbol{\nu }$ such that whp the set of independent sets in the interpolation sequence $\mathbb{G}_t, 0\le t\le T={n\choose 2}$, of graphs $\mathbb{G}(n,d/n)$ exhibits an OGP with parameters $\eta _{\mathrm{OGP}}=1/2(1+o_d(1))$ and $\boldsymbol{\nu }$. As a result, the largest independent set which can be constructed by stable algorithms is asymptotically $\eta _{\mathrm{OGP}}=1/2(1+o_d(1))$.

Similarly, for every large enough $K$, there exists a large enough $M$ and a subset $\boldsymbol{\nu }$ such that whp the set of satisfying assignments in the interpolation sequence $\boldsymbol{\Phi }_t$ of a random K-SAT formula $\boldsymbol{\Phi }_{n,\alpha n}$ exhibits an OGP with parameters $\eta _{\mathrm{OGP}}=\eta _{\mathrm{OPT}}$ and $\boldsymbol{\nu }$, when $\alpha >\alpha _{K,\mathrm{ALG}}$. As a result, the largest constraint-to-variables density for which stable algorithms can find a solution of a random K-SAT problem is $\alpha _{K,\mathrm{ALG}}$.

In particular, in both cases, the OGP provides a tight characterization of the power of stable algorithms.

For the second part of the theorem, we set $\eta _{\mathrm{OGP}} =\eta _{\mathrm{OPT}}$ since for constraint satisfaction problems there is no notion of approximation to optimality, as we consider exclusively satisfying assignments. Variants of Theorem 4.4 were proven in RV17BH22.

Why is the OGP a barrier to stable algorithms? The intuitive reason for this is a rather simple continuity-type argument and best illustrated when $K=2$ in Definition 4.2 and $\boldsymbol{\nu }$ is of the form in 5. Suppose that, in addition to the properties outlined in Definition 4.2, two additional properties hold. First, the sequence of instances $\mathcal{R}_1,\ldots ,\mathcal{R}_T$ is “continuous,” in the sense that $\mathcal{R}_t\approx \mathcal{R}_{t+1}$ for all $t$. Second, suppose that when $\ell _1=1$ and $\ell _K=\ell _2=T$ the case $d(\sigma _1,\sigma _2)\le \nu _1$ does not occur. Namely, while per Definition 4.2, for every pair $\sigma _{t_1},\sigma _{t_2}$ of solutions with values $\ge \eta _{\mathrm{OPT}}$ in the interpolated sequence $\mathcal{R}_t, t=1,\ldots ,T$, both $d(\sigma _{t_1},\sigma _{t_2})\le \nu _1$ (small distance) and $\ge \nu _2$ (large distance) situations are possible; when the two solutions correspond to the beginning (that is, $t_1=1$) and the end (that is, $t_2=T$) of the interpolation sequence, only the latter situation $d(\sigma _{t_1},\sigma _{t_2})\ge \nu _2$ (large distance) is possible. This can indeed be established in many special cases using the same techniques which verify the OGP. Then the argument establishing the OGP as an obstruction to stable algorithms goes informally as follows. Consider implementing a stable algorithm $\mathcal{ALG}$ on a sequence of interpolated instances $\mathcal{R}_t, 1\le t\le T$, resulting in a series of solutions $\mathcal{ALG}(\mathcal{R}_t)$. Since the sequence $\mathcal{R}_t$ is “continuous” (changes slowly as $t$ increases), and the algorithm is stable, then $d\left(\mathcal{ALG}(\mathcal{R}_1),\mathcal{ALG}(\mathcal{R}_t)\right)$ changes “continuously” as well, as $t$ changes from its initial value $d\left(\mathcal{ALG}(\mathcal{R}_1),\mathcal{ALG}(\mathcal{R}_1)\right)=0\le \nu _1$ to its final value $d\left(\mathcal{ALG}(\mathcal{R}_1),\mathcal{ALG}(\mathcal{R}_T)\right)$ which by a property above must be $\ge \nu _2$. For any $t$, by the OGP, we have $d\left(\mathcal{ALG}(\mathcal{R}_1),\mathcal{ALG}(\mathcal{R}_t)\right)$ is either $\le \nu _1$ or $\ge \nu _2$. Therefore, there must exist $\tau$ such that $d\left(\mathcal{ALG}(\mathcal{R}_1),\mathcal{ALG}(\mathcal{R}_\tau )\right)\le \nu _1$ and $d\left(\mathcal{ALG}(\mathcal{R}_1),\mathcal{ALG}(\mathcal{R}_{\tau +1})\right)\ge \nu _2$. By the triangle inequality, this means $d\left(\mathcal{ALG}(\mathcal{R}_\tau ),\mathcal{ALG}(\mathcal{R}_{\tau +1})\right)\ge \nu _2-\nu _1$. But this contradicts stability of $\mathcal{ALG}$. The argument is illustrated in Figure 2.

Figure 2.

OGP is an obstruction to stable algorithms. The algorithmic path (blue curve) has to jump (red section) from the small-distance ($\le \nu _1$) region to the large-distance ($\ge \nu _2$) region in the space of solutions $\Theta$.

Recall now the clustering property discussed earlier. Could one build an algorithmic barrier type argument directly from the clustering property? It appears the answer is no, as evidenced in a different model called the symmetric binary perceptron model GKPX22. In this model, which also involves constraints-to-variables parameter $\alpha$, algorithms find a satisfying solution up to some value $\alpha _{\mathrm{ALG}}>0$ and an OGP above $\alpha _{\mathrm{OGP}} \approx \alpha _{\mathrm{ALG}}>0$ is exhibited. Yet, at the same, the model exhibits clustering at every positive value $\alpha >0$. It thus appears that geometric structures which are barriers to algorithms are more involved than clustering and have to rely on more complex structures such as multioverlaps $\boldsymbol{\nu }$. We will see this in an even more pronounced form in Section 6 in the context of spin glasses.

5. Ultrametricity Guides Algorithms for Spin Glasses

We now turn to the state-of-the-art algorithms for spin glasses, our Example II, specifically algorithms for solving optimization problem 1. A dramatic progress on this problem and a fairly complete understanding of its algorithmic tractability was achieved only recently in the works of Subag Sub21, in the context of a related so-called spherical spin glass model, and of Montanari Mon19, in the context of precisely problem 1 in the special case $p=2$ (namely, the Sherrington–Kirkpatrick model). Montanari has constructed a polynomial time algorithm which achieves any desired proximity $\eta _{p, {\mathrm{OPT}}}-\epsilon$ to optimality whp, assuming an unproven though widely believed conjecture that the cumulative distribution function associated with the Parisi measure $\mu _{p,{\mathrm{OPT}}}$ solving variational problem 4 is strictly increasing within its support. Namely, the support of the measure induced by $\mu _{p,{\mathrm{OPT}}}$ is a connected subset of $[0,1]$. We will connect this property with the OGP shortly.

Soon afterward, the result was extended to general $p$ in EAMS21, and a new algorithmic threshold $\eta _{p, \mathrm{ALG}}\le \eta _{p,\mathrm{OPT}}$ corresponding to the best-known stable algorithms was identified. This is summarized in Theorem 5.1, stated informally.

Theorem 5.1 (Informal).

Let $\mathcal{L}\supset \mathcal{U}$ be the set of functions on $[0,1]$ with total variation appropriately bounded. Let $\eta _{p, \mathrm{ALG}}\triangleq \inf _{\mu \in \mathcal{L}}\mathcal{P}(\mu )\le \eta _{p, \mathrm{OPT}}=\inf _{\mu \in \mathcal{U}}\mathcal{P}(\mu )$. There exists a polynomial time stable algorithm which constructs a spin configuration $\sigma _{\mathrm{ALG}}$ achieving $\approx \eta _{p, \mathrm{ALG}}$ whp as $n\to \infty$.

When $p=2$, it is conjectured that the variational problem $\inf _{\mu \in \mathcal{L}}\mathcal{P}(\mu )$ is solved by a strictly increasing function, and thus modulo this conjecture, Montanari’s algorithm achieves near optimality. It is known however CGPR19 that this is not the case when $p\ge 4$ and $p$ is even. In this case,

$$\begin{align} \eta _{p, \mathrm{ALG}}=\inf _{\mu \in \mathcal{L}}\mathcal{P}(\mu )< \eta _{p, \mathrm{OPT}}=\inf _{\mu \in \mathcal{U}}\mathcal{P}(\mu ). {\tag{6}} \end{align}$$

We now explain the algorithmic threshold $\eta _{p, \mathrm{ALG}}$ and the successful algorithm for achieving it. The Parisi measure $\mu _{p,\mathrm{OPT}}$ is conjecturally a measure describing random overlaps of nearly optimal solutions in the following sense. Fix a positive inverse temperature parameter $\beta$, and consider sampling two spin configurations $\sigma _1,\sigma _2$ from the associated Gibbs measure independently. Then $\mu _{p,\mathrm{OPT}}$ is the distribution of the inner product (overlap) $n^{-1}\langle \sigma _1,\sigma _2\rangle$ in the limit $\beta \to \infty$. Furthermore, suppose instead $m$ samples $\sigma _j, j\in [m]$ are generated according to the same Gibbs distribution. Then the distribution of the set of overlaps $n^{-1}\langle \sigma _{j_1},\sigma _{j_2}\rangle$ is believed to be described by distributions over set of $m$ values in $[0,1]$ arranged on a tree. As such, it satisfies a strong form of the metric triangle inequality (with max replacing the addition operator) called ultrametricity, depicted in Figure 3. When the Parisi measure $\mu _{p,\mathrm{OPT}}$ is strictly increasing and there is no gap in its support, the branching points on the tree are located in a continuous way without gaps. It turns out that this coincides with the case $\inf _{\mu \in \mathcal{L}}\mathcal{P}(\mu )=\inf _{\mu \in \mathcal{U}}\mathcal{P}(\mu )$ implying $\eta _{p, \mathrm{ALG}}=\eta _{p, \mathrm{OPT}}$, and this is the case when the algorithm reaches near optimality. The algorithm is constructed as a walk in the interior of $[-1,1]^n$ which starts from the center of the cube and follows the ultrametric tree as a guide. Specifically, at every stage, it makes a small incremental step in the direction (a) maximally improving the objective function and (b) orthogonal to the previous direction. Absence of gaps implies that the algorithm proceeds without interruption until hitting a corner of $[-1,1]^n$ provably reaching near optimality. In a sense, the ultrametric tree “guides” the successive steps of the algorithm.

Figure 3.

Ultrametric tree of solutions in $[-1,1]^n$. End points of the curve correspond to corners of the cube. The Parisi measure has no gap in the support, reflected in the tree branching continously.

When the Parisi measure is not strictly increasing, the authors show that the algorithm reaches the largest energy level achievable by ultrametric trees without gaps. This was shown to be the value of the extended variational problem $\inf _{\mu \in \mathcal{L}}\mathcal{P}(\mu )$ which per 6 is strictly smaller than $\inf _{\mu \in \mathcal{U}}\mathcal{P}(\mu )$.

6. Branching-OGP: Tight Characterization of Algorithmically Solvable Spin Glasses

A final missing piece completing the algorithmic tractability versus hardness picture for spin glasses was delivered by Huang and Sellke HS22. When the Parisi measure $\mu _{p,\mathrm{OPT}}$ exhibits a gap in its support and as a result we have strict inequality 6, this inequality can be turned into an obstruction to algorithms along the lines of Theorems 4.3 and 4.4. Specifically, in the spirit of Definition 4.2, let $\boldsymbol{\nu } _{p,\mathrm{gap}}$ be the set of $M$-tuples $(\sigma _j, j\in [M]) \in \Theta ^M$ such that the set of the associated overlaps $\left(n^{-1}\langle \sigma _i,\sigma _j\rangle , i,j\in [M]\right)$ represents a continuous branching tree (in some precise sense, the kind we see in Figure 3). Then the model exhibits an OGP above the value $\eta _{p,\mathrm{ALG}}$ with $\boldsymbol{\nu }$ serving as the multi-overlap obstruction. More specifically, we have the following theorem.

Theorem 6.1.

Optimization problem 2 exhibits an OGP in the sense of Definition 4.2 with parameters $\left(\eta _{p,\mathrm{OGP}}\triangleq \eta _{p,\mathrm{ALG}},\boldsymbol{\nu }_{p,\mathrm{gap}}\right)$. This property is an obstruction to stable algorithms in finding spin configurations with values larger than $\eta _{p,\mathrm{ALG}}$, whp as $n\to \infty$.

Hence (when combined with Theorem 5.1) $\eta _{p,\mathrm{ALG}}$ represents a tight algorithmic value achievable by stable algorithms.

As a corollary, a $p$-spin model can be solved to optimality by stable algorithms when $p=2$ (modulo the strict monotonicity conjecture discussed above). On the other hand, when $p\ge 4$ is even, the largest value achievable by stable algorithms is strictly suboptimal. It is conjectured that this remains the case for all $p\ge 3$ regardless of the parity and has been verified for the case when $p$ is sufficiently large. The multi-OGP obstruction set $\boldsymbol{\nu }_{p,\mathrm{gap}}$ is called branching-OGP (B-OGP). It subsumes the obstruction sets $\boldsymbol{\nu }$ appearing in the context of Theorems 4.3 and 4.4 as special cases, and it represents the most general known form of OGP-based obstructions.

Figure 4.

Ultrametric tree of solutions when the Parisi measure has a gap in the support leading to gaps in branching locations.

The proof of Theorem 6.1 follows a path similar to the ones used in Theorems 4.3 and 4.4 and is also based on the interpolation scheme between instances, using in a crucial way the approaches developed by Talagrand in Tal06 and earlier by Guerra in Gue03. In particular, one creates an interpolation between many instances $(H(\cdot ,J_j), j\in [M])$ of the Hamiltonian arranged in a particular tree-like continuous way. One then argues that a stable algorithm implemented on this “continuous tree” of instances results in a continuous (due to stability) tree of solutions—the kind represented in Figure 3. Then a separate argument is used to show that every “continuous tree” of spin configurations $(\sigma _j, j\in [M])$ will have at least one solution with energy at most $\eta _{p, \mathrm{ALG}}$, and thus (using an additional concentration argument) stable algorithms cannot beat this threshold. Putting it differently, $M$-tuples of spin configurations achieving simultaneous values larger than $\eta _{p, \mathrm{ALG}}$ have associated ultrametric trees with gaps—as represented by the red segments in Figure 4.

It is reasonable to believe that the scope of algorithms limited by the value $\eta _{p, \mathrm{ALG}}$ is much broader than the class of stable algorithms and quite possibly includes the class of all polynomial time algorithms for these and similar problems. Such a result however is out of reach for now due to its very strong algorithmic complexity implications.