Continuing on the program that Jonathan and I started with David in PHA 1 and PHA 2, with Ewan and Holden, we put out a PHA 3 preprint on the arXiv (Potential Hessian Ascent III: Sampling the Sherrington-Kirkpatrick Model at \(\beta\) < 1/2, DLSS26). Shortly after PHA 3, Holden, Jonathan and I put out a follow-up, creatively termed (Potential Hessian Ascent IV: Sampling the Sherrington-Kirkpatrick Model at \(\beta\) < 1, LSS26) that builds on the cavity interpolation theory and free probability toolkit developed in PHA 3 to introduce a “local” SDE error analysis which only uses (and proves) “local” regularity properties for the driver and drift terms underlying the matricial functions of certain quantities that drive the SDEs.

Together, these works yield a \(o_n(1)\) TVD sampler for the SK model up to the replica-symmetric threshold \(\beta < 1\), above which sampling from the Gibbs measured is conjectured to be hard. The algorithm combines algorithmic stochastic localization (ASL) with rejection sampling over path-space via Jarzynski’s equality (JE). The analysis uses local regularity properties of the TAP Hessian and the PHA framework, developing new cavity interpolation theory and an extension to the free-probability toolkit introduced in PHA 1, Section-4, which is then combined with delicate SDE error analysis and the ability to efficiently sample from sufficiently “localized” distributions, i.e., those with large enough external field induced by a stochastic localization (SL) process that’s been run long enough.

I have decided to write, with two guest contributions from Holden, a 5-part blog post explaining the conceptual background and development of the algorithm, the main proof skeletons, and the main technical innovations.

• In this first blog post, I will explain the connection between SL and Hessian dynamics (HD), introduce the TAP free energy, discuss the over-determined system of ASL, TAP and PHD SDEs, and the natural development of the ‘‘desiderata’’ that must be proved for the algorithm to succeed (via a running example about bounding the \(W_2\) and KL-divergence between two SDEs). Holden will then introduce Jarzynski’s inequality, and discuss how it can be used to do rejection sampling to suppress the TVD error. I will end with three consequences (and perspectives) of the result – one for theoretical physicists, another for analysts/probabilists, and a final one for theoretical computer scientists.
• In the second blog post, I will spend time developing the cavity interpolation theory that gives exact moment estimates and ‘‘better-than-naive’’’ overlap-concentration along the SL process. This is necessary to prove the desideratum that the algorithmic covariance (Hessian of the TAP free energy) stays close to the true covariance along the SL path. It will turn out to be the case that this theory holds for \(\beta < 1\) even though this was not explicitly written down in PHA 3, and a couple of final details for this were not fleshed out till PHA 4. The key technical contribution here is to run a cavity interpolation for the planted SL model with SL tilt and compute “self-stability” estimates on doing Taylor expansions for bulk-deviation quantities (such as overlaps and magnetizations) to get exact moment estimates.
• In the third blog post, I will spend time explaining the extension to the analysis of the non-commutative interpolation developed in PHA 1, Section-4 to control the diagonal entries of the (squared) algorithmic covariance – this will then imply various regularity estimates and show that the ASL-TAP process closely tracks the PHD process. We will first do this for \(\beta < 1/2\), and then I will then show that, in fact, one can reason about a “regularized” resolvent which agrees with the original resolvent only over a (possibly non-convex) subset of \((-1,1)^n\) and, simultaneously, significantly simplify the free interpolation analysis in PHA 3 by using a result of Bandeira, Boedihardjo and van-Handel, 2023. This allows one to reason about the concentration and free limits of the diagonal sub-algebra of the TAP Hessian (a resolvent) all the way to \(\beta < 1\) where it can be singular in a subset of \((-1,1)^n\).
• In the fourth blog post, Holden will show how the localized distribution concentrates on a wedge of the hypercube (provided it is run for sufficiently long) and how a 2-stage decomposition on this wedge allows one to show the entropy contraction property, which is the final desideratum that needs to be proved – this is the approach adopted in PHA 3. After that, I will briefly overview how a result of Kumar et al, 2026, which allows one to sample from a Gibbs distribution with quadratic potential and sufficiently strong external field on a wedge of \(\lbrace -1,1\rbrace^n\) using an annealing-type procedure, is used in PHA 4 to directly sample from the localized distribution after invoking the fact that SL run for large time will cause the distribution to be heavilty concentrated inside a wedge of sufficient size.
• In the final blog post, I will talk about how to get access to a “safe set” of points \(S_{A}(c) \subset \lbrace -1,1\rbrace^n\) when one fixes a “good” event over the input \(A\) so that, with high probability over the SL paths, the TAP Hessian evaluated at any \(m \in S_{A}(c)\) is \(c\)-strongly convex. This is the place where, perhaps unsurpringly, we rely on a result about AMP state evolution. Essentially, we combine the fact that posterior means of the SL tilted measures are “stable” in time for most SL paths with the \(c\)-strong convexity of Hessians in small balls around AMP iterates, to build “tubes” around good SL paths that we then union to be contained in \(S_{A}(c)\). From this, one straightforwardly obtains that, with high probability over \(A\) and the SL paths, uniformly in time \(t \in [0,T]\), small enough balls around the SL-tilted means \(\lbrace m_t\rbrace_{0 \le t \le T}\) lie in \(S_{A}(c)\).
  

Table of Contents

1. Algorithm design
   • Stochastic localization and Hessian dynamics
   • Algorithmic surrogates via the TAP free energy
   • ASL-TAP and boosting to a KL divergence bound
2. Jarzynski equality
   • Rejection sampling over path space
   • Final desiderata
3. Conclusions
   • Certifying the overlap distribution
   • Weak functional inequalities and resolvent driven SDEs
   • Sampling without functional inequalities
4. Footnotes
   

Algorithm Design

Let us first reason about how to bound the error between two continuous diffusion processes run for large constant time. Think of one process as representing an “ideal” process that, if run for infinite time, will localize on a point in the support of the Gibbs measure, and think of the other as one that attempts to proxy it with a covariance that can be algorithmically computed (but incurs approximation errors).

Denote the ideal process as \(dm_t = Q(m_t)dB_t\) and the algorithmic process as \(d\hat{m}_t = \hat{Q}(\hat{m}_t)dB_t\). Then, by a simple argument using the fact that the Wasserstein-\(2\) (\(W_2\)) distance between two Gaussians is upper bounded by the Frobenius norm of the difference between their covariances, we have

This already hints at why one of the most crucial desiderata in the entire program is going to be an error estimate tracking how well a “surrogate” covariance \(\hat{Q}(\cdot)\) tracks the “ideal” covariance \(Q(\cdot)\) across SL paths \(\lbrace m_t \rbrace_t\).

It turns out, we will need to compute the error between two sets of Ito processes. The ideal processes will be a set of two coupled processes, and the algorithmic case will be similar (though not necessarily coupled). These processes are SDEs that also have drift terms. In addition to the diffusion process before, the ideal process will also include \(dy_t = m_tdt + dB_t\) and the algorithmic estimate will likewise include \(d\hat{y}_t = \hat{f}_t(\hat{m}_t)dt + dB_t\). Now, we set \(y_0 = m_0 = \hat{y}_0 = \hat{m}_0 = 0^n\) (which is the relevant initialization for us) and, assuming that the solutions exist and are well-posed, obtain the following representation for the ideal and algorithmic processes

and

We can now upper-bound the \(W_2\) distance between the ideal and algorithmic processes using the cumulative difference process \(\text{err}_t = \begin{pmatrix} y_t - \hat{y}_t \\ m_t - \hat{f}_t(\hat{m}_t) \end{pmatrix}\). Indeed, a simple application of Ito’s lemma tells us that the rate of the average instantaneous error between the two processes is

whereupon an application of a \(c\)-weighted AM-GM inequality on the first term followed by some triangle inequalities and the fact that \((a+b)^2 \le 2a^2+2b^2\) tells us that

At this point, it is clear that if we are only to run our algorithmic processes for finite time \(T\) and the final four terms in the bound above are of order \(O_{t,\beta}(1)\) at every \(0 \le t \le T\), a Gronwall’s inequality bound will immediately give that \(\mathbb{E}\|\text{err}_T\|^2_2 = O_{T,\beta}(1)\).

Now, if we stop our algorithmic processes at some large constant time \(T\), we must either assert that some deterministic and efficient function \(h(\hat{m}_T,\hat{y}_T)\) outputs a sample with \(o_n(1)\) TVD error, or we can use the pair \((\hat{m}_T,\hat{y}_T)\) as a “warm start” to another efficient algorithm which outputs a sample \(\sigma \in \lbrace -1,1\rbrace^n\) that is \(o_n(1)\) close in TVD error to the target measure. We will choose the latter approach.

Between the bounds on the four final four quantities for the cumulative error process, the required sampler with the warm start mentioned above, and the fact that the algorithmic surrogate \(\hat{Q}(\cdot)\) must be a valid covariance, we already have a list of the desiderata that the algorithm requires, at least to have \(W_2\) error that is \(O_{T,\beta}(1)\)¹:

• \(\hat{Q}(\cdot)\) is a regular and valid covariance, namely \(c(\beta)I_n \preceq \hat{Q}(m) \preceq C(\beta)I_n\) at all \(m \in \lbrace -1,1\rbrace^n\)².
• The covariance errors are small, that is \(\mathbb{E}\|\hat{Q}(m_t) - Q(m_t)\|^2_F \le O_{t,\beta}(1)\).
• The functions \(\hat{f}_t(\cdot)\) are \(C_{\beta,t}\)-Lipschitz with respect to \(\|\cdot\|_2\) inside the solid cube.
• The error PHD-TAP drift error for magnetization is bounded, that is \(\mathbb{E}\|\hat{f}_t(m)-m\|^2_2 \le O_{t,\beta}(1)\) for every \(m \in \lbrace -1,1\rbrace^n\)³.
• After running the (discretized) algorithmic processes for \((\hat{m}_t,\hat{y}_t)\) for time \(T\), there is a sampler that samples from a certain simpler (“stochastically localized”) distribution with \(o_n(1)\) TVD error in polynomial time.

While it is not clear right now, the covariance and PHD-TAP drift error estimates will only permit error of the desired order when \(\beta < 1\) – this will be discussed in the second and third blog post. Additionally, the PSD-ness property of the covariance will hold pointwise inside the cube only for \(\beta < 1/2\). To make the PSD-ness property, and the PHD-TAP drift error estimate (which relies on it in an indirect way), work for \(\beta < 1\), we will need to refine these desiderata to hold only at a certain “safe” set \(S_A(c) \subset \lbrace -1, 1\rbrace^n\) – the PHD-TAP error part of this will be the subject of the third blog post, and obtaining PSD-ness over the “safe” set will be established in the final blog post. Lastly, the Lipschitz error for \(\hat{Q}(\cdot)\) will be a consequence of the Loewner regularity of \(\hat{Q}(\cdot)\) combined with the resolvent structure it has based on the explicit choice we use – see (1.2).

We are now two steps away from getting a \(o_n(1)\)-TVD sampler, if we can show the desiderata outlined above. First, we need to reason about one more source of \(L_2\)-error to go from \(O_{T,\beta}(1)\)-\(W_2\) error to \(O_{T,\beta}(1)\)-KL divergence error via an application of Girsanov’s theorem. At that point, Pinsker’s inequality yields a \(O_{T,\beta}(1)\)-TVD error, but that is still not \(o_n(1)\). The second step, which achieves this, is to use rejection sampling (but over path space) to suppress the error further to \(o_{n}(1)\)-TVD error – this is where Jarzynski’s equality enters the picture. We will now introduce the ideal processes (SL/HD), then define the TAP free energy and use it to derive the algorithmic surrogates. At that point, we will be able to complete the ASL-TAP and PHD \(L_2\)-error bound needed to apply Girsanov’s theorem. We will then move on to defining Jarzynski’s equality, see how it allows us to do rejection sampling, and briefly overview how it suppresses the TVD error further. Doing the last step will incur one more desiderata, at which point we will conclude with our final list of desiderata.

Stochastic localization and Hessian dynamics

Given a measurable space \((\Omega, \mathcal{B}(\Omega))\) a localization process is a stochastic process \(\lbrace\nu_t(\cdot)\rbrace_{t\ge 0}\) over the space of probability measures on \((\Omega, \mathcal{B}(\Omega))\). The process has two distinct properties: 1) for any \(A \in \mathcal{B}(\Omega)\), \(\lim_{t\to\infty} \nu_t(A) \in \lbrace0,1\rbrace\), and 2) \(\nu_t(\cdot)\) is a martingale.

The first property tells us that the stochastic process “localizes” at one particular event \(A\) in the measure space, and the fact that it is a martingale implies that the process is “stationary” upon averaging. The latter property allows us to think of the localization process as a convex decomposition of the measure \(\nu_0(\cdot)\), weighted by the “tilted” measures along the way. For us, there are two important consequences of this:

1. If we can algorithmically simulate the process by “following” the tilts that generate the sequence of tilted measures for a very large (but constant) time \(T\), then we will come close to a sample from the target measure, and
2. We can choose any measure-valued process that localizes and is a martingale based on its convenience for simulating efficiently, and analzying its paths to prove our desiderata.

There is just one remaining subtlety – while we can choose a localization process that we can analyze and simulate, where do we start it? As we will see, we choose the “linear-tilt” localization scheme which essentially starts at a sample \(x_0\) drawn from the target measure, and consists of noisy observations through time \(t\) that ultimately become more informative and localize as a Dirac measure \(\delta_{x_0}\)⁴. This scheme starts at \(y_0 = 0^n\) and reveals itself at time \(t\) as

where \(x_0 \sim \nu_0\) and \(B_t \sim \mathcal{N}(0,t)\). The measure \(\nu_t(\cdot)\) induced by this scheme is

and it is easy to see that \(\lim_{t \to\infty} \nu_t(\cdot) \to \delta_{x_0}\) almost-surely and that \(\lbrace\nu_t(\cdot)\rbrace_{t\ge 0}\) is a martingale process since they form a Doob martingale. At this point, using the many known equivalent characterizations of stochastic localization (see [Sections 1 & 2, STZ26]) one can also rewrite the linear-tilt process as

where \(m_t = \mathbb{E}_{x\sim\nu_t}\left[x\right]\). A main conceptual insight in our work is that there is yet another equivalent rewrite for a stochastic localization process, and this rewrite tracks the evolution of the magnetizations/averages \(\lbrace m_t\rbrace_{t\ge 0}\) of the “tilted” measures \(\nu_t(\cdot)\) under the linear-tilt localization scheme \(\lbrace y_t\rbrace_{t\ge 0}\). This process is

where \(\mathsf{Cov}(\nu_t)\) is the covariance matrix for the measure \(\nu_t\). It is not difficult to see that this is the SDE that the magnetization process should obey from its definition as the mean of the tilted measures plus an application of Ito’s lemma in conjunction with the fact that the magnetizations must form a martingale. It will turn out to be the case (for convenience more than anything else) that, though the actual algorithm will not run an algorithmic version of this particular version of SL, it will still be very useful to analyze it to prove some of the desiderata⁵. The algorithmic version of this process is

and now a critical part of proving our desiderata and writing down the final algorithm relies on gaining access to an efficiently computable sequence of matrix-valued functions \(\hat{Q} : [-1,1]^n \to M_n(\mathbb{R})_{\text{sa}}\) – this is exactly where the contiguity to a planted model and the TAP free energy will be of assistance.

Algorithmic surrogates via the TAP free energy

Before we introduce the TAP free energy for the SK model, as is relevant in the high-temperature regime, let us quickly remind ourslves of the (random) Gibbs measure associated with it. The target/Gibbs measure for the SK model is given as

where \(A\) is a \(\mathsf{GOE}(n)\) matrix, \(Z\) is the normalizing constant called the partition function, and the density is defined for every \(\sigma \in \Sigma_n := \lbrace -1,1 \rbrace^n\).

The SK model has a particularly nice form for its free energy in the high-temperature regime (\(\beta < 1\)) whose validity was rigorously established in a series of papers – see, for instance [CP19]. The free energy is equivalent to the largest value of the TAP functional evaluated as a supremum over all possible magnetizations. The TAP functional for the SK model at some magnetization \(m \in [-1,1]^n\) is given as

where \(h(\cdot) : [-1,1] \to [0,1]\) is the entropy functional for a Bernoulli random variable (see PHA 3, (2.2)). This structure of the TAP free energy is closely related to taking a Legendre transform of the standard definition of the free energy, and consists of adding the final term as a “correction term” to the Gibbs variational principle over product measures. For our purposes, we don’t just need the TAP free energy of the standard SK model, but we need some analogue of it for any sequence of tilts \(y_t\) generated by the linear-tilt localization scheme we have – this is easily accomplished by adding an external field term. Namely, setting \(\nu_0 = \mu_{\beta A}\) and using the fact that \(\nu_t(m) \propto e^{\langle y_t, m \rangle}\mu_{\beta A}\) gives us the following modification to the TAP free energy, given a fixed tilt \(y \in \mathbb{R}^n\)

We will be interested in viewing this as a Fenchel legendre transform, and given some algorithmic estimate \(\hat{y}_t\) of \(y_t\) (with sufficiently good approximation) we will be able to compute \(\hat{m}_t\) using the strong convexity of this TAP free energy (since that will guarnatee the uniqueness of \(\hat{m}_t\) ), which will in turn allow us to compute the next tilt \(\hat{y}_{t+\delta t}\).

This is fantastic for the algorithmic procedure, and it also suggests a natural covariance \(\hat{Q}(m)\) that can be defined everywhere⁶ in \((-1,1)^n\) – the inverse of the Hessian of the TAP free energy. This is so since \(y_t\) and \(m_t\) are dual to each other, and the Crouzeix identity tells us that the second-derivative with respect to one variable (\(y_t\)) is the functional inverse of the other (\(m_t\)) and this gives

where \(D(m)\) is the diagonal matrix with entries \(\left\lbrace \frac{1}{1-(m_i)^2}\right\rbrace_{i=1}^n\). A pleasant consequence of the choice of \(\hat{Q}(m)\) is that it is actually independent of the tilt \(y\). Another immediate consequence is that, under the Loewner order sandwich on \(\hat{Q}(\cdot)\) (and consequently \(D(\cdot)\hat{Q}(\cdot)\)) required by the first desideratum for the covariance matrix, one can easily obtain that

where \(R(a) = \left(\beta^2\mathsf{tr}_n[D^{-1}(a)]I_n-\beta A - \frac{2\beta^2}{n}aa^T\right)^{-1}\). Then, conditioning on the event that \(\|A\|_\infty \le 2+\delta_\beta\) and doing some elementary estimates using the fact that \(a,b \in (-1,1)^n\), yields that \(\|R(a)-R(b)\|^2_F \le C(\beta)\|a-b\|^2_2\). Substituting this into the bound above shows that \(\hat{Q}(\cdot)\) is Lipschitz with respect to \(\|\cdot\|_F\) and proves one of the four regularity properties in the desiderata (simply as a consequence of the first desideraturm and surrogate choice of covariance).

I mentioned before that while the algorithmic version of HD, called potential Hessian dynamics (PHD), which we can now see will start at \(\hat{m}_0 = 0^n\) with updates

will not be used directly, it will be useful in the analysis. One can ask, then, how so? One answer is the reason we can’t just boost the \(O_{T,\beta}(1)\) \(W_2\) bound we can assume (given our desiderata, which we will systematically prove) to a \(O(1)\)-KL divergence bound via Girsanov’s theorem – we now focus on why this is, and how to address it.

ASL-TAP and boosting to a KL divergence bound

In a nutshell, the main reason why the \(O_{t,\beta}(1)\) \(W_2\) error can’t be boosted to a KL-divergence error via Girsanov’s theorem is that the algorithmic processes for \(\hat{y}_t\) and \(\hat{m}_t\) are not coupled in any way, but they should satisfy the constraint that \(\nabla\mathcal{F}_{\mathsf{TAP}}(\hat{m}_t,\hat{y}_t) = 0\). This requirement is the algorithmic analogue of the fact that \(y_t\) and \(m_t\) satisfy the relationship \(m_t = \mathbb{E}_{x\sim\nu_t}[x]\). If not, \(\hat{m}_t\) and \(\hat{y}_t\) are simply not coupled in a way the ideal processes are, and this is what stops us from applying Girsanov’s theorem to boost the \(W_2\) error to a KL-divergence error. It turns out that we cannot satisfy all three algorithmic processes at the same time – the algorithmic processes based on the TAP surrogacy are over-determined. We need to satisfy the TAP fixed-point, so we can either choose PHD-TAP

or ASL-TAP

where

can be explicitly derived via an application of Ito’s lemma to the function that gives a precise relationship between \(\hat{y}_t\) and \(\hat{m}_t\), which itself can be straightforwardly derived from the TAP fixed-point equation.

For the algorithm, using the strong convexity of \(\hat{Q}(\cdot)\), we choose to satisfy the TAP fixed-point constraint and run the (discretized) ASL process \(d\hat{y}_t = \hat{m}_tdt + dB_t\). For the analysis, we actually show that, under the desiderata mentioned before, PHD-TAP \((\hat{m}_t, \hat{y}_t)\) stays close to the ideal processes \((y_t,m_t)\) in the sense that

This follows by a triangle inequality, the Lipschitz-ness of \(\hat{m}_t\) (which can be straightofrwardly concluded from the fact that \(\hat{Q}\) has bonuded operator norm by the first desideratum) and the fact that

After this, one can apply Girsanov’s theorem (after an application of the DPI) to show that ASL-TAP stays close to the ideal SL process for \(y_t\).

Consequently, we must account for one additional source of error, and that is the expected squared error between the ASL-TAP and PHD process. The final desiderata to obtain \(O(1)\) KL divergence error will come from this part of the argument.

Jarzynski equality

Stuff and things

Rejection sampling over path space

Stuff and things

Final desiderata

Stuff and things

Conclusions

Stuff and things

Certifying the overlap distribution

Stuff and things

Weak functional inequalities and resolvent driven SDEs

Stuff and things

Sampling without functional inequalities

Stuff and things

FOOTNOTES