Parallel Sampling from the Ising 𝑝-Spin Model Nima Anari1 , Aniket Das1 , and Alireza Haqi1 1 Stanford University, {anari,aniketd,ahaqi}@stanford.edu
arXiv:2607.12348v1 [cs.DS] 14 Jul 2026
Abstract We study the parallel complexity of sampling from the high-temperature Ising mixed 𝑝-spin Gibbs measure, a canonical instance of a mean-field spin glass on the hypercube {±1}𝑛 . We propose two different algorithms for this problem, corresponding to two different regimes of accuracy. Our first algorithm is a parallel implementation of a Markov chain known as block dynamics, combined with an approximate rejection sampling step that uses an Ising model in a novel way as a proposal distribution to approximate the quadratic interaction terms of the 𝑝-spin Hamiltonian. For any 𝜀 > 0, this algorithm 1 runs in 𝑛 /3 polylog(𝑛/𝜀) parallel time with poly(𝑛, log(1/𝜀)) work, and outputs a sample whose law is 𝜀-close to the 𝑝-spin measure in total variation distance. Our second algorithm uses Picard iterations to parallelize the Algorithmic Stochastic Localization (ASL) process of El Alaoui, Montanari, and Sellke (2025), and for any 𝜀 > 𝜀𝑛 , takes polylog(𝑛/𝜀) parallel time and poly(𝑛/𝜀) work to produce a sample that is 𝜀-close to the 𝑝-spin measure in the normalized 2-Wasserstein metric. Here, 𝜀𝑛 > 0 is a threshold that goes to 0 as 𝑛 → ∞. Our result constitutes a doubly exponential improvement in the 𝜀 dependence of the runtime and an exponential improvement in the 𝜀 dependence of the total work when compared to naïve ASL, whose runtime scales as exp(poly(1/𝜀)).
Contents 1
Introduction 1.1 Our Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Related Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Overview of Our Techniques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2 3 4 5
2
Preliminaries 2.1 Notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Functional Inequalities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Stochastic Localization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Picard Iteration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 TAP Free Energy and the Tilted Mean Estimator . . . . . . . . . . . . . . . . . . . . . . . . . .
9 9 10 11 12 13
3
High Accuracy Parallel Sampler for the 𝑝-Spin Model 3.1 Technical Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Bounds on the Hamiltonian of the Conditional Posterior . . . . . . . . . . . . . . . . . . . . . 3.3 Analysis of Parallel Rejection Sampling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 Proof of Theorem 20 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Proof of Theorem 15 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Proof of Theorem 22 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
14 14 15 20 22 23 23
4
Low Accuracy Parallel Sampler for the 𝑝-Spin Model 4.1 Sign Rounding for Algorithmic Stochastic Localization . . . . . . . . . . . . . . . . . . . . . . 4.2 Parallel Time Analysis of Algorithm 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 TAP-AMP in polylog(𝑛) Parallel Time . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Lipschitz Property of the TAP-AMP Mean Approximation for the SK Model . . . . . . . . . .
25 25 26 29 29
1
A.1 Sign Rounding Analysis for Algorithmic Stochastic Localization
36
A.2 Parallelization of TAP-AMP Mean Approximation Algorithms
38
A.3 Lipschitz Property of the TAP Fixed-Point Computation Algorithm for Mixed 𝑝-Spin Models
41
A.4 exp(poly(1/𝜀)) Steps for Algorithmic Stochastic Localization
46
1
Introduction
Efficient sampling from a probability distribution which is specified only up to an unknown normalizing constant is a central algorithmic problem across computer science, statistical physics, and machine learning. Generally, given a finite state space Ω and a Hamiltonian 𝐻 : Ω → ℝ, the target distribution is a Gibbs measure of the following form: 𝜇(𝑥) =
1 exp(𝐻(𝑥)), 𝑍
𝑍=
Õ
exp(𝐻(𝑥))
(1)
𝑥∈Ω
where the normalizing constant 𝑍, also known as the partition function, is unknown and typically hard to compute. This setting includes classical problems in approximate counting [JVV86] and randomized algorithms such as approximating the volume of convex bodies [DFK91] and approximating matrix permanents [JSV04], Markov-chain Monte Carlo methods in statistical physics [Par81], and posterior sampling in high-dimensional statistical models. Our work studies the problem of sampling from the Gibbs measure of the Ising mixed 𝑝-spin model, which is a Gibbs measure on the hypercube {±1}𝑛 whose Hamiltonian is a random polynomial with Gaussian coefficients. 𝐻(𝑥) =
p 𝑃 Õ 𝛽 𝑝 𝑝! 𝑝−1 𝑝=2 𝑛 2
·
Õ
𝐺𝐽
𝐽⊆[𝑛],|𝐽|=𝑝
Ö
𝑥 𝑘 + ⟨𝑥, ℎ⟩,
𝑘∈𝐽
𝜇(𝑥) =
1 exp(𝐻(𝑥)) 𝑍
(2)
𝑖.𝑖.𝑑
Here, 𝑃 ≥ 2 is an integer, the disorder 𝐺 𝐽 ∼ 𝒩 (0, 1) for every choice of 𝐽, and ℎ is an arbitrary external field independent of 𝐺. The 𝑝-spin Hamiltonian can be equivalently characterized as a Gaussian process Í on {±1}𝑛 with 𝔼𝐺 [𝐻(𝑥)] = ⟨ℎ, 𝑥⟩ and Cov[𝐻(𝑥)𝐻(𝑦)] = 𝑁𝜉(⟨𝑥,𝑦 ⟩/𝑁 ), where 𝜉(𝑡) = 𝑃𝑝=2 𝛽 2𝑝 𝑡 𝑝 denotes the mixture function. The mixed 𝑝-spin model is a canonical instance of a mean-field spin glass, which plays a central role in several areas of probability, mathematical physics, and average-case complexity [Méz+88; Tal10; Pan12; MS24]. Notable special cases include the pure 𝑝-spin model where only one 𝛽 𝑝 is nonzero, and the Sherrington-Kirkpatrick (SK) model where 𝑃 = 2. We consider the problem of approximate sampling from the 𝑝-spin Gibbs measure in the quenched setting, where the disorder 𝐺 is sampled once and then fixed, and the algorithm is tasked with approximately sampling from 𝜇 given random access to the disorder coefficients (which can be stored with poly(𝑛) memory). Given an error parameter 𝜀 > 0, the goal is to design an algorithm that runs in time poly(𝑛, 1/𝜀) for any ˆ 𝜇) ≤ 𝜀 for typical realizations of 𝐺. realization of 𝐺, and produces an output distributed as 𝜇ˆ such that dist(𝜇, Here, dist denotes any appropriate metric between probability distributions. Under this setting, [Adh+22] and [Ana+24] analyze Glauber dynamics [Gla63], a canonical Markov chain for sampling from discrete distributions, which, at every step, resamples a randomly chosen coordinate according to its conditional marginal, while fixing the remaining coordinates. The authors identify a high-temperature regime for the 𝑝-spin model p under which Glauber dynamics efficiently samples from 𝜇. Specifically, they Í𝑃 show that when 𝑝=2 𝛽 𝑝 𝑝 3 ln(𝑝) = 𝑂(1), Glauber dynamics takes 𝑂(𝑛 ln(𝑛/𝜀)) time to output a sample ˆ 𝜇) ≤ 𝜀 with high probability over the randomness of 𝐺. Since the error distributed as 𝜇ˆ satisfying 𝑑TV (𝜇, parameter 𝜀 > 0 can be made as small as desired, we refer to such a guarantee as a High Accuracy Guarantee. In another line of work, [EMS25] propose a different approximate sampler for the 𝑝-spin model with a zero external field (i.e., ℎ = 0 in Eq. (2)) based on an approximate discretization of an Itô diffusion process known 2
as Eldan’s stochastic localization [Eld13]. Their algorithm, named Algorithmic Stochastic Localization (ASL), offers a weaker sampling guarantee, albeit supporting a broader range of temperatures than that covered by the high-temperature condition defined above. Concretely, there exists an error threshold 𝜀𝑛 = 𝑜 𝑛 (1) such that for any 𝜀 > 𝜀𝑛 , ASL runs in time 𝑂(poly(𝑛) exp(poly(1/𝜀))) and with high probability over 𝐺, approximately samples from 𝜇 up to error 𝜀 in the normalized 2-Wasserstein metric. In other words, the output 𝑥ˆ of their algorithm can be coupled with some 𝑥★ ∼ 𝜇 such that 𝑛 −1 𝔼[∥𝑥ˆ − 𝑥★∥] ≤ 𝜀2 . We call such a guarantee a Low Accuracy Guarantee for two key reasons. Firstly, for a fixed problem instance, the accuracy of ASL can only be certified above a threshold 𝜀 > 𝜀𝑛 . Secondly, their guarantees are presented in the normalized 2-Wasserstein metric, which is strictly weaker than the total variation metric for the range of 𝜀 tolerated by their algorithm. A common drawback of both algorithms considered above is that they are inherently sequential. Indeed, Glauber dynamics sequentially updates one coordinate after another in a random order, while Algorithmic Stochastic Localization discretizes the trajectory of a continuous Markov process, which naturally evolves in a sequential fashion. Consequently, these algorithms fail to utilize the full potential of modern hardware, which is generally capable of massive parallelism, both in terms of parallel computation and parallel memory access. This motivates the central problem of our work: Can we design an efficient approximate sampling algorithm for the Ising 𝑝-spin model that can take advantage of parallel hardware? Specifically, we consider the problem of sampling from the 𝑝-spin Gibbs measure in the PRAM model of computation, where the algorithm is allowed to query the disorder 𝐺 in multiple rounds, with each round comprising polynomially many synchronous queries that are executed in parallel. Given an error tolerance 𝜀 > 0, our goal is to design an algorithm that approximately samples from 𝜇 up to error 𝜀 in some suitable metric for typical realizations of 𝐺, while also optimizing for two criteria: Parallel Runtime or Round Complexity, which denotes the number of rounds of parallel queries made to 𝐺, and Work, which denotes the total number of queries made to 𝐺. Ideally, for every realization of 𝐺, we want the parallel runtime to scale at most sublinearly in 𝑛 and 𝜀−1 , and the work to scale at most polynomially in 𝑛 and 𝜀−1 . Compared to the complexity of sequential sampling, the theory of parallel sampling for mean-field spin glasses is far less developed. To our knowledge, the only prior parallel algorithm for sampling from the 𝑝-spin Gibbs measure is the work of [Lee23], which considers the 𝑝-spin model in the high-temperature regime and √ achieves a high-accuracy guarantee with a parallel runtime of 𝑂( 𝑛 · polylog(𝑛/𝜀)) and 𝑂(poly(𝑛/𝜀)) work. Beyond this, for the special case of the high-temperature SK model (i.e., 𝑃 = 2 and 𝛽 2 < 1/4), [Che+25] show that the Restricted Gaussian Dynamics algorithm of [CE22] is a high-accuracy approximate parallel sampler with parallel runtime 𝑂(polylog(𝑛/𝜀)) and 𝑂(poly(𝑛/𝜀)) work.
1.1
Our Contributions
In this work, we develop two different algorithms for parallel sampling from the Ising mixed 𝑝-spin Gibbs measure. The first is a high-accuracy algorithm which samples from the 𝑝-spin model up to 𝜀 total variation error in 𝑂(𝑛 1/3 polylog(𝑛/𝜀)) parallel time and 𝑂(poly(𝑛, log(1/𝜀))) work. To our knowledge, this constitutes an improvement both in terms of parallel runtime √ and work compared to the previous best known result of [Lee23], which exhibits a parallel runtime of 𝑂( 𝑛polylog(𝑛/𝜀)) with 𝑂(poly(𝑛/𝜀)) work. Theorem 1 (High-Accuracy Parallel 𝑝-spin Sampler, Informal Version of Theorem 15). Let 𝜇 denote the Ising 𝑝-spin p Gibbs measure as defined in Eq. (2). Assume the model is in the high-temperature regime, i.e., Í𝑃 𝛽 𝑝 3 ln(𝑝) ≤ 𝐶 for some 𝐶 = Θ(1). Then, there exists a parallel sampling algorithm (Algorithm 1) satisfying 𝑝 𝑝=2 the following for any 𝜀 > 0.
e 1/3 ln(1/𝜀)) with 𝑂(poly(𝑛, log(1/𝜀))) work. • The algorithm runs in parallel time 𝑂(𝑛 • With probability at least 1 − exp(−𝑐𝑛) over the randomness of the disorder 𝐺, the algorithm outputs 𝑥ˆ ∼ 𝜇ˆ such ˆ 𝜇) ≤ 𝜀. that 𝑑TV (𝜇, Our second algorithm obtains a low-accuracy approximate sampling guarantee in the normalized 2Wasserstein metric, and enjoys a parallel runtime of 𝑂(polylog(𝑛/𝜀)) with 𝑂(poly(𝑛/𝜀)) work. This significantly
3
improves upon the 𝑂(poly(𝑛) · exp(poly(1/𝜀))) runtime of [EMS25] and constitutes the first parallel algorithm for sampling from the 𝑝-spin model with polylog(𝑛/𝜀) parallel runtime in any metric. Theorem 2 (Low-Accuracy Parallel 𝑝-spin Sampler, Informal Version of Theorem 27). Let 𝜇 denote the Ising 𝑝-spin Gibbs p measure as defined in Eq. (2) with ℎ = 0. Assume the model is in the high-temperature regime, i.e., Í𝑃 3 𝑝=2 𝛽 𝑝 𝑝 ln(𝑝) ≤ 𝐶 for some 𝐶 = Θ(1). Then, there exists an 𝜀𝑛 > 0 such that 𝜀𝑛 → 0 as 𝑛 → ∞, and a parallel sampling algorithm (Algorithm 4) which satisfies the following for any 𝜀 > 𝜀𝑛 : • The algorithm runs in parallel time 𝑂(polylog(𝑛/𝜀)) with 𝑂(poly(𝑛/𝜀)) work. • With probability at least 1 − 𝑜(1) over the randomness of the disorder 𝐺, the output 𝑥ˆ of the algorithm can be coupled to some 𝑥★ ∼ 𝜇 such that 𝑛 −1 𝔼[∥𝑥ˆ − 𝑥★∥2 ] ≤ 𝜀2 . For the special case of the SK model, our low-accuracy algorithm satisfies the above guarantee for 𝛽2 ≤ 0.375. To the best of our knowledge, this constitutes the only parallel sampler for the SK model running in 𝑂(polylog(𝑛/𝜀)) time with 𝑂(poly(𝑛/𝜀)) work in this regime of inverse temperatures.
1.2
Related Work
Sampling from mean-field spin glasses. For the Sherrington–Kirkpatrick model (𝑃 = 2), classical mixing e 2 ) mixing time of Glauber dynamics only up to the criteria such as Dobrushin’s condition [Dob68] imply 𝑂(𝑛 1/2 − vanishing window 𝛽 2 = 𝑂(𝑛 ). This falls short of the dimension-free 𝛽2 regime predicted in the statistical physics literature and leads to algorithmically unsatisfactory guarantees. The first major breakthrough was e 2 ) time the work of [EKZ22] which used stochastic localization to prove that Glauber dynamics mixes in 𝑂(𝑛 e up to 𝛽 2 < 1/4. The mixing time was subsequently sharpened to 𝑂(𝑛) by [Ana+22] and [CE22], and the range e 2 ) mixing of 𝛽 2 was improved to 𝛽 2 < 0.295 by [AKV24]. For the Ising 𝑝-spin model, [Adh+22] proved an 𝑂(𝑛 p Í𝑃 time for Glauber dynamics in the high-temperature regime 𝑝=2 𝛽 𝑝 𝑝 3 ln(𝑝) = 𝑂(1). The mixing time was
e later improved to 𝑂(𝑛) by [Ana+24]. In a different direction, [EMS22] (see also [Cel24]) introduced Algorithmic Stochastic Localization (ASL) for the SK model with zero external field. For any 𝜀 > 𝜀𝑛 for a fixed 𝜀𝑛 = 𝑜 𝑛 (1), their algorithm runs in time 𝑂(poly(𝑛) exp(poly(1/𝜀))) and obtains a low-accuracy guarantee: sampling up to 𝜀-error in the normalized 2-Wasserstein distance. Their algorithm works up to 𝛽2 < 1, which is conjectured to be the sharp threshold up to which the SK model admits a polynomial-time approximate sampler. Their framework was later extended to the Ising 𝑝-spin model in [EMS25], obtaining a similar low-accuracy guarantee but covering a broader range of temperatures than that of [Adh+22] and [Ana+24]. There also exists a parallel line of work on the spherical 𝑝-spin model, where the Gibbs measure 𝜇 is supported on the sphere {𝑥 ∈ ℝ𝑛 | ∥𝑥∥2 = 𝑛}, a setting that is often more analytically tractable than the hypercube. In the high-temperature regime, [GJ19] proved that a continuous-time diffusion process known e time. Beyond the as Langevin dynamics samples from the spherical 𝑝-spin model in total variation in 𝑂(1) high-temperature regime considered in [GJ19], the recent breakthrough of [HMP24] designed an improved version of ASL for the spherical 𝑝-spin model with zero external field that samples up to 𝑜(1) error in total variation in 𝑂(poly(𝑛)) time beyond the high-temperature regime considered by [GJ19]. Building upon their techniques, [Hua+25] proved that continuous-time annealed Langevin dynamics samples up to total variation error exp(−Ω(𝑛 −1/5 )) in 𝑂(poly(𝑛)) time beyond the high-temperature regime. Parallel sampling. Parallel sampling has been a longstanding topic of interest in theoretical computer science, and more recently, in machine learning, where diffusion-based generative models are often deployed on massively parallel hardware [Lai+25; Shi+23; Hu+25]. In theoretical computer science, the problem dates back to at least [MVV87], who developed an RNC algorithm (i.e., running in 𝑂(polylog(𝑛)) parallel time on a PRAM with 𝑂(poly(𝑛)) work) for finding a perfect matching, and asked whether one can also sample a uniform random perfect matching in parallel; this remains a major open problem to date. Recent works study parallel sampling assuming access to global counting oracles. One such model isthe weighted counting oracle, where one has access to the log Laplace transform 𝐿(𝑤) = ln 𝔼𝜇 exp(⟨𝑤, 𝑥⟩) . A
4
weaker model is the unweighted counting oracle, which returns conditional marginals ℙ𝑋∼𝜇 [𝑋𝑖 = 𝑎 | 𝑋𝑆 = 𝜔𝑆 ] for 𝑆 ⊆ [𝑛], 𝑖 ∉ 𝑆, 𝑎 ∈ {±1}, and 𝜔𝑆 ∈ {±1}𝑆 . Such models are of interest for several combinatorial distributions such as spanning trees, Eulerian tours, planar perfect matchings, and determinantal point processes, where such oracles can be implemented on a PRAM in 𝑂(polylog(𝑛)) time and 𝑂(poly(𝑛)) work [Csa75]. Assuming access to weighted counting oracles, [Ana+23; ACV24] obtained RNC samplers for distributions satisfying ∥∇2 𝐿(𝑤)∥op = 𝑂(1). For more general distributions that do not satisfy such a regularity condition, [Hu+25] e 2/3 ) parallel time with polynomial work. The parallel runtime was developed an algorithm running in 𝑂(𝑛 √ e 𝑛) by [Ana+25]. In the weaker unweighted counting model, [AGR24] developed an later improved to 𝑂( √ 𝑂(𝑛 2/3 ) time parallel sampler with 𝑂(poly(𝑛)) work. The parallel runtime was later improved to 𝑂( 𝑛) by [Ana+25]. Contrary to the above, our work falls under the more challenging local oracle setting, where one has to design a parallel sampler for a Gibbs measure given access to its Hamiltonian. This constitutes a far weaker form of access compared to the global counting oracles described above. In fact, the problem of computing a counting oracle of a Gibbs measure given its Hamiltonian is NP-hard in the worst case [GŠV20], and highly nontrivial even for specialized Hamiltonians like the SK model. In this setting, [FHY21; LY22] develop RNC samplers for Ising models satisfying a Dobrushin-type condition. For mean-field spin glass models, this only covers a vanishing window of temperatures, e.g., 𝛽 2 = 𝑂(𝑛 −1/2 ) for the SK model. More recently, [Che+25] show that the Restricted Gaussian Dynamics algorithm of [CE22] is an RNC sampler for the SK model. For general Ising 𝑝-spin models, much less is known. To our knowledge, the only parallel sampler in this setting √ e 𝑛) parallel time with polynomial work. prior to our work is [Lee23], which runs in 𝑂(
1.3
Overview of Our Techniques
Algorithm 1: 𝑠-Glauber Dynamics with Parallel Rejection Sampling Sample 𝑋 (0) from 𝜇0 ∝ 𝑒 ⟨ℎ,𝑥⟩ for 𝑡 = 0, . . . , 𝑇 − 1 do Sample 𝑆 ∼ Uniform [𝑛] 𝑠 (𝑡+1)
𝑋𝑆 𝑋
(𝑡+1) 𝑆
(𝑡)
← ParallelRejectionSampler(𝑆, 𝑋 , 𝜀/10𝑇 ) 𝑆
←𝑋
(see Algorithm 2)
(𝑡) 𝑆
return 𝑋 (𝑇) Algorithm 2: Approximate Parallel Rejection Sampling with an Ising Proposal (ParallelRejectionSampler) Input: 𝑆, pinning 𝜏 on 𝑆, 𝜀step Set 𝑦 = 𝑋𝑆 and decompose the conditional Hamiltonian 𝐻𝑆,𝜏 (𝑦) = 𝐻(𝑦, 𝜏) on 𝑆 as follows 𝐻(𝑦) = ⟨𝑏 𝑆,𝜏 , 𝑦⟩ +
1 ⊤ 𝑦 𝐴𝑆,𝜏 𝑦 + 𝑅 𝑆,𝜏 (𝑦). 2
Set the rejection parameter 𝑐¯ = Θ(1), 𝐿 = Θ( 𝑐¯ ln(1/𝜀step )), 𝜀RGD = Θ(𝜀step/𝐿), 𝜂 = Θ(1), and 𝛿 = Θ(𝜀step ) ℎ̂ = CenteringExternalField(𝐴𝑆,𝜏 , 𝑏 𝑆,𝜏 , 𝑅 𝑆,𝜏 , 𝜂, 𝛿) (see Algorithm 1) for ℓ = 1 to 𝐿 do In Parallel 𝑈ℓ , 𝑉ℓ ∼ RGD 𝐴𝑆,𝜏 , 𝑏 𝑆,𝜏 + ℎ̂, 𝜀RGD iid
𝑊ℓ ∼ Uniform[0, 1] Accept 𝑈ℓ if 𝑊ℓ ≤ min 1, 𝑐¯−1 exp(𝜓(𝑈ℓ ) − 𝜓(𝑉ℓ )) where 𝜓(𝑦) = 𝑅 𝑆,𝜏 (𝑦) − ℎ̂, 𝑦 return the first accepted sample; declare FAILURE otherwise High Accuracy Sampler Based on the prior work of [Lee23], our high-accuracy sampler for the 𝑝-spin model, stated in Algorithm 1, is a computationally efficient implementation of a Markov chain called 𝑠-Glauber
5
dynamics (also known as 𝑠-block dynamics), whose kernel 𝑃𝑠 (𝑋 (𝑡) , ·) is defined as follows: 1. Sample 𝑆 uniformly at random from all subsets of [𝑛] of size 𝑠.
(𝑡+1)
2. Sample 𝑋 (𝑡+1) from the posterior 𝜇 𝑋 (𝑡+1) 𝑋𝑆
(𝑡)
= 𝑋𝑆
The above Markov chain is reversible with respect to 𝜇 and, for 𝑠 = 1, it corresponds to Glauber dynamics, a canonical sampling algorithm for discrete spaces. When 𝜇 is the 𝑝-spin Gibbs measure at high temperature, e 𝑛/𝑠 ) [Ana+24; Lee23]. Although this mixing time 𝑠-Glauber dynamics is known to exhibit a mixing time of 𝑂( improves upon increasing 𝑠, the computational complexity of implementing each step of 𝑠-Glauber dynamics worsens significantly as 𝑠 grows. Indeed, a naïve implementation of the posterior sampling step (i.e., Step 2 above) takes 𝑂(2𝑠 ) time, which is computationally prohibitive, and completely negates the potential benefits of any mixing time improvements for 𝑠 > 1. To circumvent this, our Algorithm 1 replaces the posterior sampling step in 𝑠-Glauber dynamics with an approximate rejection sampler (Algorithm 2), whose proposal distribution is a carefully chosen Ising model. The idea of using rejection sampling in the posterior sampling step also appears in the work of [Lee23], which approximates the posterior sampling step with a product proposal. While this enables computational efficiency (as product measures can be sampled in 𝑂(1) parallel time), their rejection sampling procedure √ e 1/2 ). In fails to correctly approximate the posterior when 𝑠 = 𝜔( 𝑛), and obtains a parallel runtime of 𝑂(𝑛 contrast, the choice of our Ising proposal, which serves as the key technical innovation underlying our parallel speedup, can be made to approximate the posterior to arbitrary accuracy as long as 𝑠 = 𝑂(𝑛 2/3 ), e 1/3 ). leading to an improved parallel runtime of 𝑂(𝑛 To motivate the construction of our Ising proposal, we observe that the conditional posterior distribution 𝜇𝑆,𝜏 (𝑦) ≔ 𝜇(𝑥 𝑆 = 𝑦 |𝑥 𝑆¯ = 𝜏) ∝ exp(𝐻𝑆,𝜏 (𝑦)) is also a 𝑝-spin Gibbs measure whose Hamiltonian, 𝐻𝑆,𝜏 (𝑦) = 𝐻(𝑥) for 𝑥 𝑆 = 𝑦 and 𝑥 𝑆¯ = 𝜏, is obtained by restricting 𝐻 to the coordinates in 𝑆 and pinning the coordinates in 𝑆¯ according to 𝜏. 𝐻𝑆,𝜏 is then a polynomial of the underlying disorder 𝐺, and admits the following decomposition: 𝐻𝑆,𝜏 (𝑦) = 𝑏 𝑆,𝜏 , 𝑦 +
1 𝑦, 𝐴𝑆,𝜏 𝑦 + 𝑅 𝑆,𝜏 (𝑦) 2
(3)
where 𝑏 𝑆,𝜏 , 𝐴𝑆,𝜏 , and 𝑅 𝑆,𝜏 depend on 𝐺. Motivated by this decomposition, our rejection sampler uses a proposal distribution 𝜋𝑆,𝜏 such that ln 𝜋𝑆,𝜏 serves as a quadratic approximation to 𝐻𝑆,𝜏 , i.e., 𝜋𝑆,𝜏 (𝑥) ∝ exp 21 𝑦, 𝐴𝑆,𝜏 𝑦 + 𝑏 𝑆,𝜏 + ℎ̂, 𝑦
1 . The performance and accuracy of our algorithm are governed by two
key problems: 1. How easily can we sample from the proposal? 2. How well does our rejection sampler approximate the posterior? Below, we discuss how we address each of these questions in our proof of Theorem 1. Analyzing the parallel complexity of our algorithm is a delicate problem due to the choice of our proposal distribution. Indeed, sampling from an Ising model is in itself a nontrivial problem, which can even be NP-hard in the worst case [SS12; GKK24]. To resolve this, we use the concentration properties of the 𝑝-spin Hamiltonian to prove that sampling from 𝜋𝑆,𝜏 is computationally tractable. In particular, the Hamiltonian of 𝜋𝑆,𝜏 satisfies a spectral condition which ensures that a Markov chain known as Restricted Gaussian Dynamics e mixing time [CE22]. (or RGD, see Algorithm 3) approximately samples from 𝜋𝑆,𝜏 in total variation, with 𝑂(1) e Moreover, each step of RGD can be implemented in 𝑂(1) parallel time with 𝑂(poly(𝑛)) work [Che+25]. Thus, to ensure parallel efficiency, our rejection sampler uses approximate samples from the proposal distribution, e parallel calls to RGD, thereby implying an 𝑂(1) e parallel runtime for Algorithm 2, which in drawn via 𝑂(1) e 𝑛/𝑠 ) parallel runtime for Algorithm 1. turn leads to an 𝑂( Finally, the optimal choice of 𝑠 is determined by the approximation error of our rejection sampler. There are two main sources of this error, namely, the approximate sampling error of RGD (which can be controlled to arbitrary precision without sacrificing algorithmic efficiency), and the error due to mismatch between 1 As explained in Section 3.3, the nonzero external field ℎ̂ centers our rejection sampling proposal, leading to improved accuracy.
6
Algorithm 3: Restricted Gaussian Dynamics (RGD) [CE22; Che+25] Input: symmetric interaction matrix 𝐽 ∈ ℝ𝑚×𝑚 , external field ℎ ∈ ℝ𝑚 , tolerance 𝜀 Output: Sample 𝑋 satisfying 𝑑TV Law(𝑋), 𝜈𝐽,ℎ ≤ 𝜖 where 𝜈𝐽,ℎ (𝑥) ∝ exp( 21 ⟨𝑥, 𝐽𝑥⟩ + ⟨ℎ, 𝑥⟩) Initialize 𝑋 (0) ∈ {±1}𝑚 arbitrarily Set 𝑇 = Θ(log(𝑚/𝜖)) for 𝑡 = 1, . . . , 𝑇 do Sample 𝑦 ∼ 𝒩 (𝑋 (𝑡) , 𝐽 −1 ) Sample 𝑋 (𝑡) from the product measure 𝑝(𝑥) ∝ exp ⟨𝑥, ℎ + 𝐽 𝑦⟩ return 𝑋 (𝑇) the proposal 𝜋𝑆,𝜏 and the conditional posterior 𝜇𝑆,𝜏 , which is governed by the log-density 𝜓(𝑦) = ln d𝜋𝑆,𝜏 (𝑦). 𝑆,𝜏 We show that this error can be controlled by proving exponential tail bounds on 𝜓 under 𝜋𝑆,𝜏 , which we establish via the concentration properties of the 𝑝-spin Hamiltonian and the isoperimetric properties of the proposal, and conclude that Algorithm 2 uniformly approximates the posterior sampling step of 𝑠-Glauber e 1/3 ) for Algorithm 1. dynamics for 𝑠 ≤ 𝑂(𝑛 2/3 ), leading to an overall runtime of 𝑂(𝑛 d𝜇
Low Accuracy Sampler Adapting the framework of [EMS22; EMS25], our low-accuracy parallel sampler, presented in Algorithm 4, is based on Eldan’s Stochastic Localization (SL) [Eld13], which we briefly introduce. For any probability measure 𝜇 supported on {±1}𝑛 , the associated SL process (𝑦𝑡 )𝑡≥0 is defined by the following stochastic differential equation on ℝ𝑛 .
𝔼𝑥∼𝜇 𝑥 · exp( 𝑦, 𝑥 ) 𝑚(𝑦) = 𝔼𝑥∼𝜇 exp( 𝑦, 𝑥 )
d𝑦𝑡 = 𝑚(𝑦𝑡 ) d𝑡 + d𝐵𝑡 , 𝑦0 = 0,
(4)
Here, (𝐵𝑡 )𝑡≥0 denotes the standard Brownian motion on ℝ𝑛 , and 𝑚(𝑦) corresponds to the mean of the exponential tilt of the measure 𝜇, which we call the tilted mean. A key property of SL is that 𝑡 −1 𝑦𝑡 is marginally distributed as 𝑥★ + 𝜁 𝑡 where 𝑥★ ∼ 𝜇 and 𝜁 𝑡 ∼ 𝒩 (0, 𝑡 −1 𝐼). Intuitively, 𝑡 −1 𝑦𝑡 is approximately distributed as 𝜇 for large enough 𝑡, and thus, a suitable discretization of Eq. (4), such as the Euler discretization below, can be used to approximately sample from 𝜇.
e 𝑦(𝑘+1)ℎ = e 𝑦 𝑘 ℎ + ℎ𝑚(e 𝑦 𝑘 ℎ ) + (𝐵(𝑘+1)ℎ − 𝐵 𝑘 ℎ ), 0 ≤ 𝑘 < 𝐾,
𝑦0 = 0
(5)
Being the Euler discretization of a continuous Markov process, the discrete dynamical system specified by Eq. (5) naturally evolves in a sequential fashion, and thus, naively simulating 𝐾 steps requires Θ(𝐾) calls to an oracle for the tilted mean 𝑚, even on a PRAM. To circumvent this, we employ a parallelization technique based on the concept of Picard iteration, a classical technique for proving the existence and uniqueness of ODE solutions [Per13], which was later used algorithmically by [ACV24] to develop parallel algorithms for sampling from continuous log-concave densities in ℝ𝑛 . The key idea behind Picard iteration is to view the trajectory defined by Eq. (5) as the solution to the following fixed point problem defined by an operator 𝒯 : ℝ𝑛×𝐾 → ℝ𝑛×𝐾 over discrete trajectories: ′
𝒯 (𝑧) ≔ 𝑧 ;
𝑧 ′𝑘 ℎ = ℎ
𝑘−1 Õ
𝑚(𝑧 𝑖 ℎ ) + 𝐵 𝑘 ℎ ∀ 𝑘 ∈ [𝐾],
𝑧 0′ = 0
(6)
𝑖=0
By unrolling the recurrence, one can observe that the trajectory defined by Eq. (5) is a fixed point of the operator 𝒯 , i.e., 𝒯 (e 𝑦) = e 𝑦 . Consequently, the discretized SL process can be approximated via the following fixed point iteration over discrete trajectories, which we call Picard iteration. (𝑟) e 𝑦𝑘 ℎ = ℎ
𝑘−1 Õ 𝑖=0
(𝑟−1)
𝑚(e 𝑦𝑖 ℎ
) + 𝐵 𝑘 ℎ , ∀ 𝑘 ∈ [𝐾],
7
(𝑟) e 𝑦0 = 0
(7)
Algorithm 4: Picard Algorithmic Stochastic Localization Input : Parallel depth 𝑅, steps 𝐾, step size ℎ, parameters (𝜂, 𝐾AMP , 𝐾NGD ) 𝐵0 ← 0 for 𝑘 ← 0 to 𝐾 − 1 do In Parallel 𝐵(𝑘+1)ℎ ← 𝐵 𝑘 ℎ + 𝜀𝑘 , where 𝜀𝑘 ∼ 𝒩 (0, ℎ𝐼) (0)
𝑦ˆ 𝑘 ℎ ← 0 for 𝑟 ← 1 to 𝑅 do (𝑟) 𝑦ˆ0 ← 0 for 𝑘 ← 0 to 𝐾 − 1 do In Parallel (𝑟−1) (𝑟−1) ˆ 𝑦ˆ 𝑘 ℎ 𝑚 ← TAP-AMP 𝑦ˆ 𝑘 ℎ , 𝜂, 𝑞★(𝑘 ℎ), 𝐾AMP , 𝐾NGD (see Algorithm 5) and Eq. (39) for 𝑘 ← 1 to 𝐾 do In Parallel Í (𝑟) (𝑟−1) ˆ 𝑦ˆ 𝑗 ℎ 𝑦ˆ 𝑘 ℎ ← ℎ 𝑘−1 + 𝐵𝑘 ℎ 𝑗=0 𝑚 (𝑅)
Output : 𝑥ˆ ← sign 𝑦ˆ 𝐾 ℎ
The key algorithmic benefit of this reformulation lies in its inherent parallelizability, albeit at the cost of introducing some mild approximation error. Indeed, Eq. (7) turns the problem of sequentially evaluating the recurrence in Eq. (5) into a highly parallelizable fixed point computation over the whole trajectory: given access to a parallel oracle for the tilted mean 𝑚, each Picard iteration can be implemented efficiently in (𝑟−1) parallel across all time points. In particular, at depth 𝑟, (e 𝑦 𝑘 ℎ ) 𝑘∈[𝐾] can be computed simultaneously in one parallel oracle call, after which the cumulative drift terms in Eq. (7) can be computed via parallel prefix sums in 𝑂(log(𝑘)) time and 𝑂(𝑘) work. However, to guarantee the effectiveness of this approach as a parallel sampler, one needs to quantify how well the Picard iterate e 𝑦 (𝑅) approximates its limiting trajectory e 𝑦 . To this end, one can show (see Lemma 12) that Lipschitzness of the tilted mean (which can be guaranteed in the high-temperature regime) suffices to ensure that the Picard iterates converge to the discretized SL process e ℎ) iterations suffice to approximate the discrete SL trajectory in exponentially fast. In particular, 𝑅 = Θ(𝐾 Eq. (5) to arbitrary accuracy. The key challenge in algorithmically implementing the Picard iteration in Eq. (7) arises from the approximate computation of the tilted mean 𝑚, which is highly nontrivial and is known to be NP-hard even for worst-case Hamiltonians [GŠV20]. To circumvent this, we make use of the TAP-AMP algorithm of [EMS23] (Algorithm 5), ˆ of the tilted mean for points along the continuous SL trajectory. This which constructs an estimator 𝑚 results in the following inexact fixed point iteration, which we call Picard Algorithmic Stochastic Localization (Algorithm 4). (𝑟) 𝑦ˆ 𝑘 ℎ = ℎ
𝑘−1 Õ
(𝑟−1)
ˆ 𝑦ˆ 𝑖 ℎ 𝑚(
(𝑟)
) + 𝐵 𝑘 ℎ , ∀ 𝑘 ∈ [𝐾], 𝑦ˆ0 = 0
(8)
𝑖=0
While standard perturbation-based analyses of inexact fixed point iterations rely on uniformly controlling the deviation between the ideal and the inexact iterates, such a strategy is not applicable to Eq. (8) due to the ˆ In fact, the current best known analysis of the absence of uniform guarantees on the estimation error of 𝑚. ˆ 𝑡 ) − 𝑚(𝑦𝑡 )∥2 ≤ 𝑛𝜀𝑛2 for points tilted mean estimator can certify only a very weak error bound of the form ∥𝑚(𝑦 along the continuous SL trajectory. Here, 𝜀𝑛 is a fixed error threshold which converges to 0 as 𝑛 → ∞. This subtle distinction constitutes a major obstacle in our analysis since the Picard trajectories at low depth (i.e., low values of 𝑟) are expected to deviate significantly from the continuous SL trajectory. Consequently, the (𝑟) (𝑟) ˆ 𝑦ˆ 𝑘 ℎ ) − 𝑚( 𝑦ˆ 𝑘 ℎ )∥. currently available guarantees are insufficient for controlling ∥𝑚( To circumvent this obstacle, we depart from the standard perturbation-based analysis of inexact fixed point iteration, and instead analyze Eq. (8) directly. In particular, we prove that in the high-temperature regime, ˆ the tilted mean estimator 𝑚(𝑦) itself is uniformly 𝑂(1)-Lipschitz. Hence, the inexact Picard iteration in Eq. (8)
8
Algorithm 5: AMP Algorithm for Tilted Mean Estimation (TAP-AMP) Input : Iterations 𝐾 AMP , 𝐾NGD , step size 𝜂, overlap parameter 𝑞 ˆ −1 ← 0, 𝑧 0 ← 0 𝑚 for 𝑘 ← 0 to 𝐾AMP − 1 do ˆ 𝑘 ← tanh(𝑧 𝑘 ) 𝑚 Í 𝑞ˆ 𝑘 ← 𝑛1 𝑛𝑖=1 tanh2 (𝑧 𝑖,𝑘 ) 𝑏 𝑘 ← (1 − 𝑞ˆ 𝑘 ) 𝜉′′ ( 𝑞ˆ 𝑘 ) ˆ 𝑘 ) + 𝑦 − 𝑏𝑘 𝑚 ˆ 𝑘−1 𝑧 𝑘+1 ← 𝛽∇𝐻(𝑚
[EMS25, Algorithm 1]
ˆ +,0 ← 𝑚 ˆ 𝐾AMP 𝑢 0 ← 𝑧 𝐾AMP , 𝑚 for 𝑘 ← 0 to 𝐾NGD − 1 do ˆ +,𝑘 ; 𝑦, 𝑞) 𝑢 𝑘+1 ← 𝑢 𝑘 − 𝜂∇ℱ̂TAP (𝑚 (see Eq. (43)) ˆ +,𝑘+1 ← tanh(𝑢 𝑘+1 ) 𝑚 ˆ +,𝐾NGD Output : 𝑚
e ℎ) ˆ 𝑘 ℎ ) + (𝐵(𝑘+1)ℎ − 𝐵 𝑘 ℎ ) in 𝑅 = 𝑂(𝐾 converges exponentially fast to the discrete process 𝑦ˆ(𝑘+1)ℎ = 𝑦ˆ 𝑘 ℎ + ℎ 𝑚(𝑦 ˆ and iterations. We then prove a sign rounding stability guarantee for the limiting discrete trajectory 𝑦, (𝑅) ˆ 𝑦 show that for any 𝜀 ≥ 𝜀𝑛 , the distribution of sign( 𝐾 ℎ /𝐾 ℎ ) approximates the 𝑝-spin measure 𝜇 up to 𝜀 error in normalized 2-Wasserstein distance whenever 𝐾 ℎ = 𝑂(log(𝑛/𝜀)) and ℎ = 𝑂(poly(1/𝜀)). Since the 𝐾 iterations are executed in parallel in each round, Algorithm 4 exhibits a parallel runtime of 𝑂(polylog(𝑛/𝜀)) with 𝑂(poly(𝑛/𝜀)) work. Compared to the 𝑂(poly(𝑛) exp(poly(1/𝜀))) runtime of Algorithmic Stochastic Localization, this constitutes a doubly exponential improvement in runtime and an exponential improvement in work in terms of 𝜀.
2
Preliminaries
2.1
Notation
For differentiable functions 𝑔 : ℝ𝑑 → ℝ, we use ∇𝑔 and ∇2 𝑔 to denote their Euclidean gradient and Hessian, respectively. For functions 𝑓 : {±1}𝑚 → ℝ, we define the discrete partial derivative, and the associated discrete gradient and discrete Hessian as follows: 𝜕𝑖 𝑓 (𝑥) =
1 · [ 𝑓 (𝑥 ¬𝑖 ; 𝑥 𝑖 = +1) − 𝑓 (𝑥 ¬𝑖 ; 𝑥 𝑖 = −1)], 2
(∇ 𝑓 (𝑥))𝑖 = 𝜕𝑖 𝑓 (𝑥),
For any 𝑥 ∈ ℝ𝑚 and multi-set 𝐽 ⊆ [𝑛], we define 𝑥 𝐽 = space Ω, we denote their total variation distance as 𝑑TV (𝑃, 𝑄) =
Î
(∇2 𝑓 (𝑥))𝑖𝑗 = 𝜕𝑖 𝜕 𝑗 𝑓 (𝑥)
(9)
𝑘∈𝐽 𝑥 𝑘 . For any two distributions 𝑃 and 𝑄 on a state
1Õ |𝑃(𝑥) − 𝑄(𝑥)| 2
(10)
𝑥∈Ω
When clear from the context, for random variables 𝑥, 𝑦 we use 𝑑TV (𝑥, 𝑦) to denote 𝑑TV (Law(𝑥), Law(𝑦)). Definition 3 (Normalized 2-Wasserstein distance). For probability measures 𝜇, 𝜈 on ℝ𝑑 with a finite second moment, which shall always be equipped with the Euclidean metric unless stated otherwise, we define the 2-Wasserstein distance 𝑊2 (𝜇, 𝜈) and its normalized variant 𝑊2,𝑑 (𝜇, 𝜈) as follows: 𝑊2 (𝜇, 𝜈) =
inf 𝒞 ∈Π(𝜇,𝜈)
𝔼(𝑋 ,𝑌)∼𝒞 ∥𝑋 − 𝑌∥
2
1/2
,
1 𝑊2,𝑑 (𝜇, 𝜈) = √ · 𝑊2 (𝜇, 𝜈) 𝑑
where Π(𝜇, 𝜈) denotes the set of all couplings of the measures 𝜇 and 𝜈.
9
(11)
Throughout, we use 𝐻 and 𝜇 to denote the Ising 𝑝-spin Hamiltonian and its associated Gibbs measure as Í defined in Eq. (2). We use 𝜉(𝑡) = 𝑃𝑝=2 𝛽 2𝑝 𝑡 𝑝 to denote its mixture function. We also define the coefficients ℭ(𝛽) and 𝔇(𝛽) as follows: ℭ(𝛽) B
𝑃 Õ
𝛽𝑝
q
𝑝 3 ln(𝑝) ;
𝔇(𝛽) B
𝑝=2
2.2
𝑃 Õ
q
𝛽 𝑝 2𝑝 𝑝 3 ln(𝑝)
(12)
𝑝=2
Functional Inequalities
Definition 4 (Poincaré and Log-Sobolev Inequalities). A measure 𝜋 on {±1}𝑚 is said to satisfy a Poincaré inequality with constant 𝜆PI if for any 𝑓 : {±1}𝑚 → ℝ, Var𝜋 [ 𝑓 ] ≤ 𝜆PI 𝔼𝜋 ∥∇ 𝑓 ∥2 .
(13)
𝜋 is said to satisfy a Log-Sobolev inequality with constant 𝜆LSI if for any 𝑓 : {±1}𝑚 → ℝ, Ent𝜋 [ 𝑓 2 ] ≤ 𝜆LSI 𝔼𝜋 ∥∇ 𝑓 ∥2 .
(14)
Furthermore, LSI implies PI with 𝜆PI ≤ 𝜆LSI . Definition 5 (Approximate Tensorization of Entropy). A measure 𝜋 on {±1}𝑚 is said to satisfy 𝐶-approximate tensorization of entropy (or 𝐶-ATE) if for any 𝑓 : {±1}𝑚 → ℝ, Ent𝜋 [ 𝑓 2 ] ≤ 𝐶
𝑚 Õ
𝔼𝜋 Ent𝜇(𝑥 𝑖 |𝑥−𝑖 ) [ 𝑓 2 ] .
(15)
𝑖=1
It is easy to show that any probability measure on {±1} satisfies an LSI with a constant of 1/2 [BB19]. Hence, any measure 𝜋 satisfying 𝐶-approximate tensorization of entropy also satisfies an LSI with 𝜆LSI = 𝐶/2. The following Lipschitz concentration guarantee for measures satisfying an LSI is standard, and follows directly from the Herbst argument [Van14, Ch. 3]. Theorem 6 (Lipschitz Concentration under LSI). Let 𝜋 be a measure on {±1}𝑚 satisfying an LSI. Then, for any 𝑓 : {±1}𝑚 → ℝ, which satisfies max ∥∇ 𝑓 (𝑥)∥ ≤ 𝐺, the following holds: 𝑥∈{±1}𝑚
𝜋 𝑓 − 𝔼𝜋 [ 𝑓 ] ≥ 𝑡 ≤ 2 exp −
𝑡2
(16)
2𝜆LSI 𝐺 2
The following Bernstein-type concentration bound for measures satisfying an ATE follows from Propositions 2.16 and 2.18 of [SS19]. Theorem 7 (Bernstein-Type Concentration for ATE Measures). Let 𝜋 be a measure on {±1}𝑚 satisfying 𝐶approximate tensorization of entropy. Then, for any 𝑓 : {±1}𝑚 → ℝ,
ª 𝑡2 𝑡 © 𝑐1 𝜋 𝑓 − 𝔼𝜋 [ 𝑓 ] ≥ 𝑡 ≤ 2 exp− min , ®. 𝐶 max ∥∇2 𝑓 (𝑥)∥𝐹 𝔼𝜋 ∥∇ 𝑓 ∥2 𝑥∈{±1} 𝑚 « ¬
(17)
The following theorem establishing ATE for high-temperature Ising models is implicit in [Ana+24] and also appears as Theorem 4.1 in [Lee23]. Theorem 8 (ATE for High-Temperature Ising Models). Let 𝜈𝐽,ℎ (𝑥) ∝ exp 12 ⟨𝑥, 𝐽𝑥⟩ + ⟨ℎ, 𝑥⟩ be an Ising model on {±1}𝑚 with ∥𝐽∥op < 1. Then, 𝜈𝐽,ℎ satisfies 𝐶-ATE with 𝐶 = (1 − ∥𝐽∥op )−1 . Consequently, 𝜈𝐽,ℎ satisfies an LSI (and hence, PI), with 𝜆LSI = 12 · (1 − ∥𝐽∥op )−1 .
10
2.3
Stochastic Localization
Let 𝜇 denote any arbitrary measure on ℝ𝑛 with a finite second moment. Stochastic localization, introduced by [Eld13], denotes the following diffusion process: d𝑦𝑡 = 𝑚(𝑦𝑡 , 𝑡) d𝑡 + d𝐵𝑡 ,
𝑦0 = 0 i √ 𝑚(𝑦, 𝑡) = 𝔼𝑥∼𝜇,𝜁∼𝒩 (0,1) 𝑥 𝑡𝑥 + 𝑡𝜁 = 𝑦
h
(18)
Stochastic localization exhibits several interesting properties that have found innumerable applications in areas such as convex geometry and the analysis of Markov chains [LV24; EKZ22; CE22]. Among these, the following property, which appears in [EM22], will be of particular importance to us. Theorem 9 (Equivalent Characterization of SL [EM22]). Let (𝑦𝑡 )𝑡≥0 denote the stochastic localization process associated with the measure 𝜇, as defined in Eq. (18). Then, there exists an 𝑥★ ∼ 𝜇 and a standard Brownian motion (𝑊𝑡 )𝑡≥0 such that for any 𝑡 ≥ 0, 𝑦𝑡 = 𝑡𝑥★ + 𝑊𝑡 . Consequently, the marginal law of 𝑡 −1 𝑦𝑡 is equal to that of 𝑥★ + 𝜁 𝑡 where 𝜁 𝑡 ∼ 𝒩 (0, 𝑡 −1 𝐼). The following lemmas, which follow from Theorem 9, are vital for analyzing the rounding scheme used in Algorithm 4. Lemma 10 (Rounding SL for the hypercube). Let 𝜇 be a measure supported on {±1}𝑛 , and let (𝑦𝑡 )𝑡≥0 denote the 2 associated SL process. For any 𝑇 > 0, define 𝑥ˆ𝑇 = sign(𝑇 −1 𝑦𝑇 ). Then, 𝑊2,𝑛 (Law( 𝑥ˆ𝑇 ), 𝜇) ≤ 8 exp(−𝑇/2). Proof. By Theorem 9, there exist 𝑥 ∗ ∼ 𝜇 and 𝑧 ∼ 𝒩 (0, 𝑇 −1 𝐼) such that 𝑇 −1 𝑦𝑇 = 𝑥 ∗ + 𝑧. Then,
1 𝔼 ∥𝑥ˆ𝑇 − 𝑥 ∗ ∥2 𝑛 𝑛 2i 1Õ h ∗ 𝔼 𝑥 𝑖 − sign(𝑥 ∗𝑖 + 𝑧 𝑖 ) = 𝑛
2 𝑊2,𝑛 Law( 𝑥ˆ𝑇 ), 𝜇 ≤
=
≤
1 𝑛 1 𝑛
𝑖=1 𝑛 Õ 𝑖=1 𝑛 Õ
h
4ℙ 𝑥 ∗𝑖 ≠ sign(𝑥 ∗𝑖 + 𝑧 𝑖 )
4ℙ |𝑧 𝑖 | ≥ 1
(19) (20)
i (21)
(22)
𝑖=1
≤ 8𝑒 −𝑇/2 .
(23)
Lemma 11 (Stability of SL Rounding). Let 𝜇 be a distribution on {±1}𝑛 , and let (𝑦𝑡 )𝑡≥0 denote the associated SL process. For any 𝑇 > 0, let 𝑧ˆ 𝑇 be a random variable satisfying 𝑊22 (ˆ𝑧𝑇 , 𝑇 −1 𝑦𝑇 ) ≤ 𝑛Δ2 . Then, 2 𝑊2,𝑛 sign(ˆ𝑧𝑇 ), sign(𝑇 −1 𝑦𝑇 ) ≤ 8𝑒 −𝑇/8 + 16Δ2 .
(24)
∗ −1 −1 ∗ −1 Proof. By Theorem 9, there exist 𝑥−1 ∼ 𝜇2and 𝜁 ∼2 𝒩 (0, 𝑇 𝐼) such that 𝑇 𝑦𝑇 = 𝑥 + 𝜁. Couple 𝑧ˆ 𝑇 and 𝑇 𝑦𝑇 𝑊2 -optimally such that 𝔼 ∥ˆ𝑧𝑇 − 𝑇 𝑦𝑇 ∥ ≤ 𝑛Δ . Then,
𝑊22 sign(ˆ𝑧𝑇 ), sign(𝑇 −1 𝑦𝑇 ) ≤ 𝔼 ∥ sign(ˆ𝑧𝑇 ) − sign(𝑇 −1 𝑦𝑇 )∥2
=4
𝑛 Õ ℙ sign(ˆ𝑧𝑇,𝑖 ) ≠ sign(𝑇 −1 𝑦𝑇,𝑖 ) . 𝑖=1
11
(25)
Observe that if |𝑇 −1 𝑦𝑇,𝑖 | > 1/2 and |ˆ𝑧𝑇,𝑖 − 𝑇 −1 𝑦𝑇,𝑖 | < 1/2, then sign(ˆ𝑧𝑇,𝑖 ) = sign(𝑇 −1 𝑦𝑇,𝑖 ). It follows that,
ℙ sign(ˆ𝑧𝑇,𝑖 ) ≠ sign(𝑇 −1 𝑦𝑇,𝑖 ) ≤ ℙ |𝑇 −1 𝑦𝑇,𝑖 | ≤ 1/2 + ℙ |ˆ𝑧𝑇,𝑖 − 𝑇 −1 𝑦𝑇,𝑖 | ≥ 1/2 ≤ ℙ |𝑥 ∗𝑖 + 𝜁 𝑖 | ≤ 1/2 + 4𝔼 (ˆ𝑧𝑇,𝑖 − 𝑇 −1 𝑦𝑇,𝑖 )2 ≤ ℙ |𝜁 𝑖 | ≥ 1/2 + 4𝔼 (ˆ𝑧𝑇,𝑖 − 𝑇 −1 𝑦𝑇,𝑖 )2 ≤ 2𝑒 −𝑇/8 + 4𝔼 (ˆ𝑧𝑇,𝑖 − 𝑇 −1 𝑦𝑇,𝑖 )2 .
(26) (27) (28) (29)
Substituting the above into Eq. (25),
16 𝔼 ∥ˆ𝑧𝑇 − 𝑇 −1 𝑦𝑇 ∥2 𝑛 ≤ 8𝑒 −𝑇/8 + 16Δ2 .
2 𝑊2,𝑛 sign(ˆ𝑧𝑇 ), sign(𝑇 −1 𝑦𝑇 ) ≤ 8𝑒 −𝑇/8 +
2.4
(30) (31)
Picard Iteration
In its canonical form, Picard iteration is a well-known analytic technique for proving the existence and uniqueness of solutions to well-posed ordinary differential equations [Per13]. Recently, it has been used as a technique for parallelizing Langevin-type algorithms for continuous sampling [ACV24] as well as denoising diffusion models in machine learning [Shi+23]. In this section, we present a concise introduction to the discrete-time analog of Picard iterations, which is the version used throughout this work. Consider a discrete-time dynamical system (𝑥 𝑘 ℎ )0≤𝑘≤𝐾 on ℝ𝑛 defined by the following iteration: 𝑥(𝑘+1)ℎ = 𝑥 𝑘 ℎ + ℎ 𝑓 (𝑘 ℎ, 𝑥 𝑘 ℎ ) + 𝑔(𝑘 ℎ),
𝑘 ∈ {0, . . . , 𝐾 − 1}
(32)
where ℎ > 0 denotes the step-size and 𝑇 = 𝐾 ℎ denotes the time horizon, and the functions 𝑓 and 𝑔 can potentially be random. We observe that solving the difference equation above is an inherently sequential task. Concretely, given an initialization 𝑥0 and access to an oracle for computing 𝑓 and 𝑔, computing the trajectory (𝑥 𝑘 ℎ )0≤𝑘≤𝐾 requires Θ(𝐾) oracle calls. The key idea in Picard iteration is to view the entire trajectory (𝑥 𝑘 ℎ )0≤𝑘≤𝐾 as the fixed point of a trajectory-level iteration in ℝ𝑛×𝐾 such that, given access to a parallel oracle for computing 𝑓 and 𝑔, each round of this trajectory-level iteration can be computed in Θ(1) oracle calls. To begin, we unroll the recurrence in Eq. (32) and observe that the solution (𝑥 𝑘 ℎ )0≤𝑘≤𝐾 satisfies the following: 𝑥 𝑘 ℎ = 𝑥0 +
𝑘−1 Õ
ℎ 𝑓 (𝑗 ℎ, 𝑥 𝑗 ℎ ) + 𝑔(𝑗 ℎ) ,
0≤𝑘≤𝐾
(33)
𝑗=0 (𝑟)
The Picard iteration associated with Eq. (32) is the following sequence of trajectories (𝑥 𝑘 ℎ )0≤𝑘≤𝐾 , 𝑟 ∈ ℕ . (𝑟+1)
𝑥𝑘ℎ
= 𝑥0 +
𝑘−1 Õ
(𝑟)
ℎ 𝑓 (𝑗 ℎ, 𝑥 𝑗 ℎ ) + 𝑔(𝑗 ℎ)
(34)
𝑗=0
We make two key observations regarding the above trajectory sequence: (𝑟)
• The solution 𝑥 𝑘 ℎ to Eq. (32) is a fixed point of Eq. (34). Indeed, setting 𝑥 𝑘 ℎ = 𝑥 𝑘 ℎ for every 𝑘 and using (𝑟+1)
Eq. (33), we conclude that 𝑥 𝑘 ℎ
= 𝑥𝑘ℎ .
• Each round of the iteration in Eq. (34) can be computed using Θ(1) calls to a parallel oracle for computing 𝑓 and 𝑔. Indeed, the values of 𝑔(𝑗 ℎ), 0 ≤ 𝑗 < 𝐾 can be precomputed using one parallel oracle call, and (𝑟) in each round, the values of 𝑓 (𝑗 ℎ, 𝑥 𝑗 ℎ ), 0 ≤ 𝑗 < 𝐾 can be computed via one call to a parallel oracle to 𝑓 .
12
The parallel efficiency of Picard iteration then depends on how fast the sequence of trajectories in Picard iteration converges to its fixed point trajectory. Below, we prove that Picard iteration exhibits an exponential convergence rate whenever 𝑓 (𝑡, 𝑥) is uniformly 𝐿-Lipschitz in 𝑥. Lemma 12 (Exponential convergence of Picard iteration under Lipschitzness). Let (𝑥 𝑘 ℎ )0≤𝑘≤𝐾 denote the solution (𝑟) to Eq. (32), where 𝑇 = 𝐾 ℎ, and let (𝑥 𝑘 ℎ )0≤𝑘≤𝐾, 𝑟≥0 denote the associated Picard iteration as defined in Eq. (34). Suppose (0)
𝑓 (𝑡, ·) is uniformly 𝐿-Lipschitz in 𝑥 and define 𝑀 = max0≤𝑘≤𝐾 ∥𝑥 𝑘 ℎ − 𝑥 𝑘 ℎ ∥. Then, for any 𝑟 ≥ 0 and 0 ≤ 𝑘 ≤ 𝐾, (𝑘 ℎ𝐿)𝑟 · 𝑀. 𝑟!
(𝑟)
∥𝑥 𝑘 ℎ − 𝑥 𝑘 ℎ ∥ ≤
(35)
In particular, 𝑅 = max{𝑒 2 𝐿𝑇, ln(𝑀/𝜀)} Picard iterations suffice to ensure the following: (𝑅)
max ∥𝑥 𝑘 ℎ − 𝑥 𝑘 ℎ ∥ ≤ 𝜀.
(36)
𝑘∈{0,...,𝐾}
(𝑘 ℎ𝐿)𝑟
(𝑟)
Proof. Let Δ 𝑘,𝑟 = ∥𝑥 𝑘 ℎ − 𝑥 𝑘 ℎ ∥. We shall prove that Δ 𝑘,𝑟 ≤ 𝑟! · 𝑀 by induction on 𝑟. Clearly, this is true for 𝑟 = 0. Now, suppose this holds for 𝑗 = 1, . . . , 𝑟. Then, by the uniform 𝐿-Lipschitzness of 𝑓 (𝑡, ·), (𝑟+1)
Δ 𝑘,𝑟+1 = ∥𝑥 𝑘 ℎ ≤ℎ
𝑘−1 Õ
− 𝑥𝑘ℎ∥ (𝑟)
∥ 𝑓 (𝑗 ℎ, 𝑥 𝑗 ℎ ) − 𝑓 (𝑗 ℎ, 𝑥 𝑗 ℎ )∥
𝑗=0
≤ ℎ𝐿
𝑘−1 Õ
Δ 𝑗,𝑟
𝑗=0 𝑘−1
(ℎ𝐿)𝑟+1 Õ 𝑟 𝑗 ≤𝑀 𝑟! 𝑗=0
≤𝑀 ≤ Hence, by induction, Δ 𝑘,𝑟 ≤
(ℎ𝐿)𝑟+1
∫ 𝑘
𝑟!
0
𝑡 𝑟 𝑑𝑡
(𝑘 ℎ𝐿)𝑟+1 𝑀. (𝑟 + 1)!
(37)
(𝑘 ℎ𝐿)𝑟 𝑀. Now, for 𝑅 ≥ max{𝑒 2 𝐿𝑇, ln(𝑀/𝜀)}, 𝑟!
max Δ 𝑘,𝑅 ≤ 𝑀 · 0≤𝑘≤𝐾
𝑒𝐿𝑇 𝑅 (𝐾 ℎ𝐿)𝑅 ≤𝑀· ≤ 𝑀𝑒 −𝑅 ≤ 𝜀. 𝑅! 𝑅
(38)
Generally, given some accuracy parameter 𝜀, 𝑇 = Θ̃(1) and ℎ = 𝑂(poly(1/𝜀)). Hence, the naive sequential e e 1/𝜀)) oracle calls. In contrast, when 𝑓 is 𝑂(1)-Lipschitz, algorithm for solving Eq. (32) requires 𝐾 = 𝑂(poly( approximately solving Eq. (32) to 𝑂(poly(𝜀)) error (which generally suffices for downstream applications) e 1/𝜀)) parallel oracle calls, leading to a significant parallel speedup at the cost of a benign requires 𝑅 = 𝑂(ln( (and easily tunable) approximation error.
2.5
TAP Free Energy and the Tilted Mean Estimator
The construction of the tilted mean estimator (Algorithm 5) of El Alaoui, Montanari, and Sellke [EMS25] is motivated by a deep structural property of the 𝑝-spin model: The mean of the tilted Gibbs measure at high temperature is an approximate stationary point of the Thouless–Anderson–Palmer (TAP) free energy [TAP77; Méz+88; Tal10], which we briefly introduce below.
13
Recall that the mixture function of the 𝑝-spin model is defined as 𝜉(𝑡) = asymptotic overlap 𝑞★(𝑡) is defined as follows:
′
𝑞 𝑘+1 (𝑡) = 𝔼𝑊∼𝑁(0,1) tanh 𝑡 + 𝜉 (𝑞 𝑘 ) + 𝑊
p
𝜉′ (𝑞 𝑘 ) + 𝑡
2
,
Í𝑃
𝑞 𝑘 (𝑡) = 0,
2 𝑝 𝑝=2 𝛽 𝑝 𝑡 .
For any 𝑡 ≥ 0, the
𝑞★(𝑡) = lim 𝑞 𝑘 (𝑡)
(39)
𝑘→∞
Since 𝑞 𝑘 (𝑡) is independent of 𝐺, computable to high accuracy via a one-dimensional Gaussian integral, and approaches the limit 𝑞★(𝑡) exponentially fast, we assume, for simplicity, that 𝑞★(𝑡) is known exactly for 𝑡 = 𝑘 ℎ, 𝑘 ∈ [𝐾]. The TAP free energy associated with the 𝑝-spin Hamiltonian in Eq. (2) is defined as follows: ℱTAP (𝑚, 𝑦) = −𝐻(𝑚) − 𝑦, 𝑚 −
𝑛 Õ
ℎ(𝑚 𝑖 ) − ONS(𝑄(𝑚))
(40)
𝑖=1
𝑛 [𝜉(1) − 𝜉(𝑄) − (1 − 𝑄)𝜉′ (𝑄)] 2 1+𝑚 1+𝑚 1−𝑚 1−𝑚 −1 2 𝑄(𝑚) = 𝑛 ∥𝑚∥ , ℎ(𝑚) = − ln ln − 2 2 2 2
ONS(𝑞) =
(41) (42)
For the sake of tractability, El Alaoui, Montanari, and Sellke [EMS25] also define an approximate version of the TAP free energy which replaces the ONS(𝑄(𝑚)) term by a quadratic approximation. ℱ̂TAP (𝑚 ; 𝑦, 𝑞) = −𝐻(𝑚) − 𝑦, 𝑚 −
𝑛 Õ
ℎ(𝑚 𝑖 ) − ONS(𝑞) − ONS′ (𝑞)(𝑄(𝑚) − 𝑞) +
𝑖=1
𝑛Γ (𝑄(𝑚) − 𝑞)2 8
(43)
where 𝑞 is set to be the asymptotic overlap 𝑞★(𝑡) for 𝑡 = 𝑘 ℎ, and Γ = Θ(1) is a tunable parameter depending only on 𝜉. [EMS25] show that the mean of the tilted Gibbs measure continues to be an approximate stationary point of ℱ̂TAP . The tilted mean estimator in Algorithm 5 then works in two phases. The first phase runs ˆ 𝐾AMP for ℱ̂TAP around which it is a series of AMP updates to compute an approximate stationary point 𝑚 locally strongly convex. The second stage exploits this local strong convexity to refine the AMP estimate via natural gradient descent. Lemma 13 (tanh is 1-Lipschitz). For all 𝑎, 𝑏 ∈ ℝ𝑛 , ∥tanh(𝑎) − tanh(𝑏)∥ ≤ ∥𝑎 − 𝑏∥. Lemma 14 (sech2 is globally 1-Lipschitz on ℝ). For all 𝑠, 𝑡 ∈ ℝ, 4 sech2 (𝑠) − sech2 (𝑡) ≤ |𝑠 − 𝑡| √ . 3 3
3
High Accuracy Parallel Sampler for the 𝑝-Spin Model
In this section, we analyze Algorithm 1 and prove the following Theorem. Theorem 15 (Analysis of Algorithm 1). For any 𝜀 > 0, the output 𝑋 (𝑇) of Algorithm 1 run with 𝑠 = Θ(𝑛 2/3 log(𝑛/𝜀)−2/3 ) and 𝑇 = Θ(𝑛 1/3 log(𝑛/𝜀)5/3 ) satisfies 𝑑TV (Law(𝑋 (𝑇) ), 𝜇) ≤ 𝜀 with probability at least 1 − 𝑒 −𝑐𝑛 . Moreover, Algorithm 1 has a parallel runtime of 𝑂(𝑛 1/3 polylog(𝑛/𝜀)) with 𝑂(poly(𝑛, log(𝑛/𝜀))) work.
3.1
Technical Preliminaries
To analyze the mixing time of Algorithm 1, we use the following result on the rapid mixing of 𝑠-Glauber dynamics for the 𝑝-spin model at high temperature. This result is implicit in [Ana+24, Theorem 6] and also appears in [Lee23].
14
Theorem 16 (Rapid Mixing of 𝑠-Glauber Dynamics). Let 𝑃𝑠 denote the kernel of 𝑠-Glauber dynamics for the 𝑝-spin Gibbs measure. Then, there exists an absolute constant 𝐶 > 0 such that if ℭ(𝛽) ≤ 𝐶, then the following holds with probability at least 1 − exp(−Ω(𝑛)) for any initial law 𝜈 𝐷KL 𝜈𝑃𝑠𝑡 || 𝜇 ≤ exp(− 𝛼𝑠𝑡 𝑛 )𝐷KL 𝜈 || 𝜇
(44)
where 𝛼 is a constant depending only on 𝔇(𝛽). We also use the following result, which proves that Restricted Gaussian Dynamics samples from an Ising model on {±1}𝑚 up to 𝜀 accuracy in total variation in polylog(𝑚/𝜀) time and poly(𝑚, log(1/𝜀)) work. This result appears in [Che+25], but with a suboptimal work dependence of poly(𝑚/𝜀). Theorem 17 (RNC Sampling of Ising Models via Restricted Gaussian Dynamics). Consider the Ising model 𝜈𝐽,ℎ (𝑥) ∝ exp( 12 ⟨𝑥, 𝐽𝑥⟩ + ⟨ℎ, 𝑥⟩) where 𝐽 ∈ ℝ𝑚×𝑚 satisfies ∥𝐽∥op ≤ 1 − 𝑐 for some numerical constant 𝑐. Then, Algo rithm 3 outputs a sample 𝑋 (𝑇) satisfying 𝑑TV 𝑋 (𝑇) , 𝜈𝐽,ℎ ≤ 𝜖 in 𝑂(log3 (𝑚/𝜀)) parallel time using 𝑂(poly(𝑚, log(𝑚/𝜀))) work. Proof. By adding a constant multiple of the identity if necessary (which does not change the distribution 𝜈𝐽,ℎ ), we can assume that 2𝑐 𝐼 ⪯ 𝐽 ⪯ (1 − 2𝑐 )𝐼. By Proposition 16, Proposition 27, and Theorem 49 of [CE22], taking 𝑇 = Θ(log(𝑚/𝜀)) suffices to ensure 𝑑TV 𝑋 (𝑇) , 𝜈𝐽,ℎ ≤ 𝜖. By [Csa75], computing 𝐽 −1 and sampling from the Gaussian density 𝒩 (𝑋 (𝑡) , 𝐽 −1 ) can be performed on a PRAM in 𝑂(log2 (𝑚)) parallel time and 𝑂(𝑚 𝜔 ) work, where 2 ≤ 𝜔 ≤ 3 is the matrix multiplication exponent. Sampling from the product measure 𝑝(𝑥) ∝ exp( 𝑥, ℎ + 𝐽 𝑦 ) can be done in 𝑂(1) parallel time on a PRAM with 𝑂(𝑚) work. Hence, the overall parallel runtime is 𝑂(log3 (𝑚/𝜀)) with work 𝑂(poly(𝑚, log(𝑚/𝜀))).
3.2
Bounds on the Hamiltonian of the Conditional Posterior
𝑛 Consider any arbitrary 𝑆 ∈ [𝑛] 𝑠 . For any 𝑥 ∈ {±1} , let 𝑥 = (𝑦, 𝜏) where 𝑦 = 𝑥 𝑆 denotes the block spins in 𝑆 and 𝜏 = 𝑥 𝑆¯ denotes the pinning outside 𝑆. Then the conditional distribution 𝜇𝑆,𝜏 (𝑦) = 𝜇(𝑦 |𝑥 𝑆¯ = 𝜏) of the spins in 𝑆 satisfies 𝜇𝑆,𝜏 (𝑦) ∝ exp(𝐻𝑆,𝜏 (𝑦)) where:
𝐻𝑆,𝜏 (𝑦) = 𝑏 𝑆,𝜏 , 𝑦 + 21 𝑦, 𝐴𝑆,𝜏 𝑦 + 𝑅 𝑆,𝜏 (𝑦)
(45)
Here, 𝐴𝑆,𝜏 is a symmetric matrix with zero diagonals with: (𝐴𝑆,𝜏 )𝑖𝑗 =
p 𝑃 Õ 𝛽 𝑝 𝑝!
Õ
𝑝−1 ¯ 𝑝=2 𝑛 2 𝐾⊆𝑆,|𝐾|=𝑝−2
𝐺 𝐾∪{𝑖,𝑗} 𝜏𝐾 ,
𝑖, 𝑗 ∈ 𝑆, 𝑖 ≠ 𝑗
(46)
The remainder term is given by 𝑅 𝑆,𝜏 (𝑦) =
p 𝑃 Õ 𝛽 𝑝 𝑝!
Õ
𝑛 (𝑝−1)/2 𝐽⊆[𝑛],|𝐽|=𝑝, 𝑝=3
¯
𝐺 𝐽 𝑦 𝐽∩𝑆 𝜏 𝐽∩𝑆
(47)
|𝐽∩𝑆|≥3
As discussed earlier, each step of 𝑠-Glauber dynamics requires sampling from 𝜇𝑆,𝜏 , which can be computationally expensive. We address this via a parallelized implementation of rejection sampling in Algorithm 2 where 𝜇𝑆,𝜏 (𝑦) is approximated by an Ising model of the form 𝜋𝑆,𝜏 (𝑦) ∝ exp( 12 𝑦, 𝐴𝑆,𝜏 𝑦 + 𝑏 𝑆,𝜏 + ℎ̂, 𝑦 ), where ℎ̂ is an external field that satisfies ℎ̂ ≈ 𝔼 𝑦∼𝜋𝑆,𝜏 [∇𝑅 𝑆,𝜏 (𝑦)]. To ensure the parallel efficiency of Algorithm 2, we generate approximate samples from the Ising proposal 𝜋𝑆,𝜏 using Restricted Gaussian Dynamics (Algorithm 3). To this end, the following theorem proves that ∥𝐴𝑆,𝜏 ∥op ≤ 1/4 holds with high probability uniformly over all choices of 𝑆 and 𝜏 whenever 𝑠 = 𝑂(𝑛 2/3 ). Hence, e parallel time. Algorithm 3 can approximately sample from the proposal in 𝑂(1)
15
Theorem 18 (Uniform Operator Norm Bound for Ising Proposal). There exists a numerical constant 𝐶 > 0 such that for ℭ(𝛽) ≤ 𝐶 and 𝑠 = 𝑜(𝑛), the following holds with probability at least 1 − exp(−𝑐𝑛): ∥𝐴𝑆,𝜏 ∥op ≤
1 , 4
∀𝑆∈
[𝑛] ¯ , 𝜏 ∈ {±1}𝑆 𝑠
(48)
Proof. From Eq. (46),
(𝐴𝑆,𝜏 )𝑖𝑗 =
¯
p 𝑃 Õ 𝛽 𝑝 𝑝!
Õ
𝑛 (𝑝−1)/2 𝑝=2
𝐺 𝐾∪{𝑖,𝑗} 𝜏𝐾 ,
𝑖, 𝑗 ∈ 𝑆, 𝑖 ≠ 𝑗.
(49)
𝐾⊆𝑆¯ |𝐾|=𝑝−2
Since 𝜏 ∈ {±1}𝑆 and 𝑔𝐾∪{𝑖,𝑗} ∼ 𝒩 (0, 1) for each 𝐾 and 𝑝, (𝐴𝑆,𝜏 )𝑖𝑗 are independent centered Gaussians with iid
Var (𝐴𝑆,𝜏 )𝑖𝑗 =
𝑃
𝑃 𝛽 2 𝑝! Õ 𝑝 𝑛−𝑠 𝑝=2
𝑛 𝑝−1
ℭ(𝛽)2 1 Õ 2 𝑝! 𝛽𝑝 ≤ . ≤ 𝑛 𝑛 𝑝−2 (𝑝 − 2)!
(50)
𝑝=2
Let 𝒩 be a 1/4-net of the unit sphere in ℝ𝑆 . By a standard covering argument, ∥𝐴𝑆,𝜏 ∥ ≤ 2 max𝑢,𝑣∈𝒩 𝑢 ⊤ 𝐴𝑆,𝜏 𝑣. Since (𝐴𝑆,𝜏 )𝑖𝑗 are independent centered Gaussians for 𝑖 ≠ 𝑗, and (𝐴𝑆,𝜏 )𝑖𝑖 = 0, 𝑢 ⊤ 𝐴𝑆,𝜏 𝑣 is a centered Gaussian for any 𝑢, 𝑣 ∈ 𝒩 with, Var 𝑢 ⊤ 𝐴𝑆,𝜏 𝑣 =
Õ
Var (𝐴𝑆,𝜏 )𝑖𝑗 (𝑢𝑖 𝑣 𝑗 )2 ≤
𝑖≠𝑗
ℭ(𝛽)2 . 𝑛
(51)
Since |𝒩 | ≤ 9𝑠 , by a union bound
𝑛𝑡 2 ℙ ∥𝐴𝑆,𝜏 ∥op ≥ 𝑡 ≤ 92𝑠 exp − ℭ(𝛽) 2 .
(52)
Taking a union bound over all (𝑆, 𝜏) and choosing ℭ(𝛽) = 𝑂(1) sufficiently small,
𝑛 𝑛−𝑠 2𝑠 𝑐𝑛 ℙ ∃ (𝑆, 𝜏) s.t. ∥𝐴𝑆,𝜏 ∥op ≥ 14 ≤ 2 9 exp − ℭ(𝛽) 2 𝑠 𝑐𝑛
≤ exp 𝑛 ln 2 + 𝑐1 𝑠 ln 𝑛 − ℭ(𝛽)2 ≤ exp(−𝑐2 𝑛).
(53) (54)
Naturally, the success probability of our rejection sampling procedure is governed by how well the Ising proposal 𝜋𝑆,𝜏 approximates the conditional distribution 𝜇𝑆,𝜏 , or, in other words, how large the influence of the remainder term 𝑅 𝑆,𝜏 is. We quantify this via a uniform bound on the discrete Hessian of 𝑅 𝑆,𝜏 . Theorem 19 (Uniform Bound on Hessian of Remainder Term). There exists a numerical constant 𝐶 > 0 such that √ for ℭ(𝛽) ≤ 𝐶, the following holds with probability at least 1 − exp(−𝑐𝑛) whenever Ω( 𝑛) ≤ 𝑠 ≤ 𝑜(𝑛) 𝐶1 ℭ(𝛽)2 𝑠 3 ∥∇ 𝑅 𝑆,𝜏 (𝑦)∥2𝐹 ≤ , 2
[𝑛] ∀𝑆∈ , (𝑦, 𝜏) ∈ {±1}𝑛 𝑠
2
𝑛
(55)
Moreover, there exists a 𝐶𝑃 = Θ(1) such that the following holds with probability at least 1 − exp(−𝑐𝑛) ∥∇𝑅 𝑆,𝜏 (𝑦)∥∞ ≤ 𝑛
𝐶𝑃
,
[𝑛] ∀𝑆∈ , (𝑦, 𝜏) ∈ {±1}𝑛 𝑠
16
(56)
Proof. We first derive a uniform Frobenius norm bound on the Hessian. From Eq. (47),
p 𝑃 Õ 𝛽 𝑝 𝑝!
∇ 𝑅 𝑆,𝜏 𝑖𝑗 =
2
Õ
𝑝−1 𝑝=3 𝑛 2 𝐽⊆[𝑛], |𝐽|=𝑝 |𝐽∩𝑆|≥3 {𝑖,𝑗}⊆𝐽
𝐺 𝐽 𝜏 𝐽∩𝑆 𝑦 𝐽∩𝑆\{𝑖,𝑗} .
(57)
Let 𝒥 = {𝐽 ⊆ [𝑛] : 3 ≤ |𝐽| ≤ 𝑃} and let 𝑔 = (𝐺 𝐽 )𝐽∈𝒥 . Note that 𝑔 is a standard Gaussian in ℝ𝒥 . Writing 𝑥 = (𝑦, 𝜏), define the matrix 𝐵 = (𝐵𝐼,𝐽 )𝐼∈(𝑆), 𝐽∈𝒥 by 2
p
𝛽 |𝐽| |𝐽|!
𝐵𝐼,𝐽 =
|𝐽|−1 𝑛 2
1(𝐼 ⊆ 𝐽, |𝐽 ∩ 𝑆| ≥ 3) 𝑥 𝐽 .
(58)
Note that for 𝐼 = {𝑖, 𝑗}, 𝑖, 𝑗 ∈ 𝑆, ∇2 𝑅 𝑆,𝜏 (𝑦) 𝑖𝑗 =
p Õ 𝛽|𝐽| |𝐽|!
1(𝐼 ⊆ 𝐽, |𝐽 ∩ 𝑆| ≥ 3) 𝑥 𝐽\𝐼 𝑔𝐽 = |𝐽|−1
𝐽∈𝒥
𝑛
2
1 1 Õ 𝐵𝐼,𝐽 𝑔𝐽 = 𝐼 (𝐵𝑔)𝐼 . 𝑥 𝐼 𝐽∈𝒥 𝑥
(59)
Since 𝑥 ∈ {±1}𝑛 , ∥∇2 𝑅 𝑆,𝜏 (𝑦)∥2𝐹 = ∥𝐵𝑔∥2 . Let 𝑀 = 𝐵𝐵⊤ . Then, by the Hanson-Wright inequality (Vershynin [Ver18], Theorem 6.2.1),
n Tr(𝑀)2 Tr(𝑀) o , ℙ ∥∇2 𝑅 𝑆,𝜏 (𝑦)∥2𝐹 ≥ 2Tr(𝑀) ≤ exp −𝑐 min . 2 ∥𝑀∥𝐹
∥𝑀∥op
(60)
We now compute Tr(𝑀), ∥𝑀∥op and ∥𝑀∥𝐹 . Since 𝑥 ∈ {±1}𝑛 , we obtain the following from Eq. (58) for any 𝐼1 , 𝐼2 ∈ 𝑆2 ,
𝑀𝐼1 ,𝐼2 =
Õ
𝐵𝐼1 ,𝐽 𝐵𝐼2 ,𝐽 =
𝑃 𝛽 2 𝑝! Õ 𝑝 𝑝=3
𝐽∈𝒥
𝑛 𝑝−1
𝑁𝑝 (𝐼1 , 𝐼2 ),
(61)
o [𝑛] 𝑁𝑝 (𝐼1 , 𝐼2 ) = # 𝐽 ∈ : 𝐼1 ∪ 𝐼2 ⊆ 𝐽, |𝐽 ∩ 𝑆| ≥ 3 . 𝑝 n
(62)
Clearly, 𝑁𝑝 (𝐼1 , 𝐼2 ) depends only on |𝐼1 ∩ 𝐼2 | ∈ {0, 1, 2}. Hence, we consider three cases. Case 1: 𝐼1 = 𝐼2 = 𝐼. Write 𝐽 = 𝐼 ∪ 𝐿 where 𝐿 ∩ 𝐼 = ∅. Since |𝐽 ∩ 𝑆| ≥ 3, we have |𝐿 ∩ 𝑆| ≥ 1, i.e. |𝐿 ∩ (𝑆 \ 𝐼)| ≥ 1. Hence, 𝑁𝑝 (𝐼1 , 𝐼2 ) = #{𝐿 ⊆ [𝑛] \ 𝐼 : |𝐿| = 𝑝 − 2} − #{𝐿 ⊆ [𝑛] \ 𝐼 : 𝐿 ∩ (𝑆 \ 𝐼) = ∅, |𝐿| = 𝑝 − 2}
=
𝑛−2 𝑛−𝑠 − . 𝑝−2 𝑝−2
(63)
Case 2: |𝐼1 ∩ 𝐼2 | = 1. Since 𝐼1 ∪ 𝐼2 ⊆ 𝑆 and |𝐼1 ∪ 𝐼2 | = 3, any 𝑝-set 𝐽 ⊇ 𝐼1 ∪ 𝐼2 satisfies |𝐽 ∩ 𝑆| ≥ 3. Hence, [𝑛] 𝑛−3 𝑁𝑝 (𝐼1 , 𝐼2 ) = #{𝐽 ∈ : 𝐼1 ∪ 𝐼2 ⊆ 𝐽} = . 𝑝 𝑝−3
Case 3: 𝐼1 ∩ 𝐼2 = ∅.
17
(64)
Then 𝐼1 ∪ 𝐼2 ⊆ 𝑆 and |𝐼1 ∪ 𝐼2 | = 4. If 𝑝 = 3, no such 𝐽 can exist. If 𝑝 ≥ 4, any 𝑝-set 𝐽 containing 𝐼1 ∪ 𝐼2 satisfies |𝐽 ∩ 𝑆| ≥ 3. Hence,
( 𝑁𝑝 (𝐼1 , 𝐼2 ) =
𝑛−4 𝑝−4 ,
𝑝 ≥ 4,
0,
otherwise.
(65)
To combine the three cases, we define Γ0 =
𝑃 𝛽 2 𝑝! Õ 𝑝 𝑛−4 𝑝=4
Γ1 =
𝑃 𝛽 2 𝑝! Õ 𝑝 𝑛−3 𝑝=3
Γ2 =
𝑛 𝑝−1 𝑝 − 4
𝑛 𝑝−1 𝑝 − 3
,
(66)
,
(67)
" 𝑃 𝛽 2 𝑝! Õ 𝑝 𝑛−2 𝑝=3
𝑛−𝑠 − 𝑝−2 𝑝−2
𝑛 𝑝−1
#
.
(68)
Then, by Eq. (61), Eq. (63), Eq. (64) and Eq. (65), we obtain the following:
Γ2 , 𝐼 1 = 𝐼 2 , 𝑀𝐼1 ,𝐼2 = Γ1 , |𝐼1 ∩ 𝐼2 | = 1, Γ0 , 𝐼1 ∩ 𝐼2 = ∅.
(69)
𝑠 Γ2 . 2
(70)
It follows that, Tr(𝑀) =
Õ
𝑀𝐼,𝐼 =
𝐼∈ 𝑆2
()
From Eq. (69), we observe that 𝑀 is symmetric with non-negative entries and constant row sums. Hence, by the Perron-Frobenius theorem, ∥𝑀∥op = max 𝐼1
Õ
𝑀𝐼1 ,𝐼2 .
(71)
𝐼2
For any 𝐼1 ∈ 𝑆2 , there are exactly 2(𝑠 − 2) possible 𝐼2 such that |𝐼1 ∩ 𝐼2 | = 1 and 𝐼1 ∩ 𝐼2 = ∅. Hence,
𝑠−2 possible 𝐼2 such that 2
𝑠−2 ∥𝑀∥op = Γ2 + 2(𝑠 − 2)Γ1 + Γ0 . 2
(72)
By a similar argument, ∥𝑀∥2𝐹 =
Õ
𝑀𝐼21 ,𝐼2
(73)
𝐼1 ,𝐼2
"
#
𝑠 𝑠−2 2 = Γ22 + 2(𝑠 − 2)Γ21 + Γ0 . 2 2
(74)
To apply Eq. (60), we need to lower bound Tr(𝑀)/∥𝑀∥op and Tr(𝑀)2 /∥𝑀∥2𝐹 . To this end, define 𝐴3 =
𝑃 Õ
𝛽 2𝑝 𝑝(𝑝 − 1)(𝑝 − 2) ≤ ℭ(𝛽)2 .
𝑝=2
18
(75)
Then, Γ1 =
1 Õ 2 𝑝! 𝐴3 ≤ 2 𝛽 = 2, 𝑝−3 𝑛 𝑝=3 𝑝 (𝑝 − 3)! 𝑛
𝑛 𝑝−1
𝑝=3
𝑃 𝛽 2 𝑝! Õ 𝑝 𝑛−4
Γ0 =
𝑃
𝑃 𝛽 2 𝑝! Õ 𝑝 𝑛−3
𝑝=4
𝑛 𝑝−1 𝑝 − 4
(76)
𝑃
≤
1 Õ 2 𝑝! 𝑃𝐴3 𝛽𝑝 ≤ 3 . 3 𝑛 𝑝=4 (𝑝 − 4)! 𝑛
(77)
To control Γ2 , recall the identity 𝑄 𝑝 :=
𝑠−3
Õ 𝑛−𝑠+𝑖 𝑛−2 𝑛−𝑠 − = . 𝑝−2 𝑝−2 𝑝−3
(78)
𝑖=0
Then, 𝑛−𝑠 𝑛−3 (𝑠 − 2) ≤ 𝑄 𝑝 ≤ (𝑠 − 2) . 𝑝−3 𝑝−3
(79)
It follows that, Γ2 =
𝑃 𝛽 2 𝑝! Õ 𝑝
𝑃 𝛽 2 𝑝! Õ 𝑝 𝑛−3
𝑛 𝑝=3
𝑛 𝑝−1 𝑝 − 3 𝑝=3
𝑄 ≤ (𝑠 − 2) 𝑝−1 𝑝
≤
𝑠𝐴3 . 𝑛2
(80)
Similarly, we can lower bound Γ2 as follows: Γ2 ≥ (𝑠 − 2)
𝑝=3
For 𝑛 ≥ Θ(𝑃), since 𝑠 = 𝑜(𝑛),
𝑝−4 2 𝑃 𝑠 + 𝑗 𝑠 − 2 Õ 𝛽 𝑝 𝑝! Ö ≥ 1 − . 𝑝−3 𝑛 𝑛 2 𝑝=3 (𝑝 − 3)! 𝑗=0
𝑃 𝛽 2 𝑝! Õ 𝑝 𝑛−𝑠 𝑛 𝑝−1
(81)
𝑠+𝑗 1 𝑛 ≤ /2 for any 𝑗 ≤ 𝑃 − 4. Hence, from Eq. (80) and Eq. (81), we conclude:
𝑐𝑃
𝑠𝐴3 𝑠𝐴3 ≤ Γ2 ≤ 2 . 𝑛2 𝑛
(82)
From Eq. (70), Eq. (72), Eq. (73) Eq. (76), Eq. (77) and Eq. (82), we obtain: Tr(𝑀) ≥ 𝑐 𝑃
𝑠 3 𝐴3 , 𝑛2
∥𝑀∥op ≤ 𝐶𝑃
𝑠𝐴3 , 𝑛2
∥𝑀∥2𝐹 ≤ 𝐶𝑃
𝑠 4 𝐴23 𝑛4
.
(83)
It follows that,
n Tr(𝑀) Tr(𝑀)2 o min
∥𝑀∥op
,
∥𝑀∥2𝐹
≥ 𝑐𝑃 𝑠 2 .
(84)
Moreover, by Eq. (70) and Eq. (82) 𝑠 3 ℭ(𝛽)2 𝑠 𝑠 3 𝐴3 Γ2 ≤ ≤ . 2 𝑛2 𝑛2
Tr(𝑀) =
(85)
Substituting Eq. (84) and Eq. (85) into Eq. (60), we obtain:
h 𝑠 3 ℭ(𝛽)2 i ℙ ∥∇2 𝑅 𝑆,𝜏 (𝑦)∥2𝐹 ≥ ≤ exp(−𝑐 𝑃 𝑠 2 ). 2 𝑛
19
(86)
√ Finally, taking a union bound over all 𝑆 ⊆ [𝑛], |𝑆| = 𝑠, and (𝑦, 𝜏) ∈ {±1}𝑛 and using 𝑠 ≥ Ω( 𝑛), we obtain:
" ℙ max
max
) (𝑦,𝜏)∈{±1}
𝑆∈ [𝑛] 𝑠
(
𝑠 3 ℭ(𝛽)2 ∥∇ 𝑅 𝑆,𝜏 (𝑦)∥2𝐹 ≥ 2 𝑛
#
𝑛 𝑛 2 exp(−𝑐 𝑃 𝑠 2 ) 𝑠
≤
2
𝑛
(87)
≤ exp(−𝑐 𝑃 𝑠 2 + 𝑐 1 𝑛 + 𝑐 2 𝑠 log 𝑛)
(88)
≤ exp(−𝑐 𝑃 𝑠 ),
(89)
2
To derive the uniform bound on the gradient, we observe that. 𝑃 Õ 𝑛
≤ 𝐶𝑃 𝑛 𝑃 .
(90)
h i 2 ℙ max |𝐺 𝐽 | > 𝑛 ≤ 𝐶𝑃 𝑛 𝑃 𝑒 −𝑛 /2 ≤ 𝑒 −𝑐𝑛 .
(91)
{𝐽 ⊆ [𝑛] : 2 ≤ |𝐽| ≤ 𝑃} =
𝑝
𝑝=2
Then, by a union bound and Gaussian concentration, 𝐽⊆[𝑛] 2≤|𝐽|≤𝑃
Then, by Eq. (47), 𝜕𝑖 𝑅 𝑆,𝜏 (𝑦) =
p 𝑃 Õ 𝛽 𝑝 𝑝! 𝑝−1 𝑝=3 𝑛 2
Õ
𝐺 𝐽 𝜏 𝐽∩𝑆 𝑦 (𝐽∩𝑆)\{𝑖} .
(92)
𝐽⊆[𝑛], |𝐽|=𝑝 |𝐽∩𝑆|≥3, 𝐽∋𝑖
Thus, on the event max𝐽 |𝑔𝐽 | ≤ 𝑛, |𝜕𝑖 𝑅 𝑆,𝜏 (𝑦)| ≤
p 𝑃 Õ 𝛽 𝑝 𝑝! 𝑝−1 𝑝=3 𝑛 2
Õ
|𝑔𝐽 | ≤ 𝑛
𝐽⊆[𝑛] |𝐽|=𝑝, 𝑖∈𝐽
p 𝑃 Õ 𝛽 𝑝 𝑝! 𝑛 − 1 𝑝−1 𝑝=3 𝑛 2
𝑝−1
≤ 𝑛 𝐶𝑃 .
(93)
Since the above bound holds uniformly in (𝑆, 𝜏, 𝑦, 𝑖), we conclude that there exists a 𝐶𝑃 = Θ(1) such that sup𝑆,𝜏,𝑦 ∥∇𝑅 𝑆,𝜏 (𝑦)∥∞ ≤ 𝑛 𝐶𝑃 .
3.3
Analysis of Parallel Rejection Sampling
In this section, we shall establish the accuracy and parallel efficiency of Algorithm 2. Henceforth, we condition on the event ℰ that the bounds in Theorem 18 and Theorem 19 hold uniformly over all 𝑆 and 𝜏. 𝑆 ˜ 𝑆,𝜏 denote the law of Theorem 20 (Analysis of Parallel Rejection Sampler). For any 𝑆 ∈ [𝑛] 𝑠 and 𝜏 ∈ {±1} , let 𝜇 the sample output by Algorithm 2. Conditioned on the event ℰ, the following holds whenever 𝑠 = Θ(𝑛 2/3 ):
[𝑛] ∀𝑆 ∈ , 𝜏 ∈ {±1}𝑆 . 𝑠
𝑑TV 𝜇˜ 𝑆,𝜏 , 𝜇𝑆,𝜏 ≤ 𝜀step
(94)
Moreover, Algorithm 2 has parallel runtime 𝑂(polylog(𝑛/𝜀step )) with 𝑂(poly(𝑛, log(1/𝜀step ))) work. The proof of Theorem 20, which is presented in Section 3.3.1, involves several technical components, the first of which is an accuracy guarantee for the approximate rejection sampling step in Algorithm 2. This result appears as Lemma 4.4 in Lee [Lee23] and Lemma 2 in [FYC23]. 𝑑𝑃 Lemma 21 (Accuracy of Rejection Sampling). Let 𝑃 and 𝑄 be probability measures satisfying 𝑑𝑄 (𝑦) ∝ 𝑒 𝜓(𝑦) .
Let 𝑌, 𝑍 ∼ 𝑄, and define 𝑅 = exp(𝜓(𝑌) − 𝜓(𝑍)). Consider the following rejection sampling procedure: Draw 𝑈 ∼ Uniform[0, 1] and accept 𝑌 if 𝑈 ≤ min{1, 𝑅/𝑐 } for some 𝑐 ≥ 1. Let 𝑃˜ denote the law of the accepted sample. Then, iid
𝑑TV 𝑃, 𝑃˜ ≤
1 · 𝔼[max{𝑅 − 𝑐, 0}], 𝔼[𝑅]
1 ℙ 𝑌 is accepted = · 𝔼[min{𝑅, 𝑐}] ≥ 2𝑐1 𝑐
20
(95)
Algorithm 1 Centering External Field Input: 𝐴𝑆,𝜏 , 𝑏 𝑆,𝜏 , 𝑅 𝑆,𝜏 , accuracy 𝜂 > 0, failure probability 𝛿 ∈ (0, 1) Set iterations 𝑇 = Θ log(𝑛/𝜂) and sample size 𝐾 = Θ 𝑠𝜂−2 log(𝑠𝑇/𝛿) Initialize 𝑧 0 ← 0 for 𝑡 = 0, . . . , 𝑇 − 1 do for 𝑖 = 1, . . . , 𝐾 do In parallel (𝑡) 𝑈 𝑖 ∼ RGD(𝐴𝑆,𝜏 , 𝑏 𝑆,𝜏 + 𝑧 𝑡 , 𝛿/100𝐾𝑇 ) 𝑧 𝑡+1 = 𝐾1
Í𝐾
(𝑡) 𝑖=1 ∇𝑅 𝑆,𝜏 (𝑈 𝑖 )
return ℎ̂ 𝑆,𝜏 = 𝑧𝑇 In addition, suppose 𝜓(𝑌) − 𝜓(𝑍) exhibits subexponential tails of the form ℙ[𝜓(𝑌) − 𝜓(𝑍) ≥ 𝑡] ≤ 𝑎 exp(−𝑏𝑡) for 𝑏 > 1. Then, the following holds: 𝑑TV 𝑃, 𝑃˜ ≤
𝑎 · 𝑐 −(𝑏−1) 𝑏−1
(96)
In Algorithm 2, (𝑋 , 𝑍) corresponds to (𝑈ℓ , 𝑉ℓ ), both of which are approximately distributed as 𝜋𝑆,𝜏 (modulo 𝑑𝜇 the negligible sampling error of Algorithm 3, which we ignore for now). Since 𝑑𝜋𝑆,𝜏 (𝑦) ∝ 𝑒 𝜓(𝑦) , where 𝑆,𝜏 𝜓(𝑦) = 𝑅 𝑆,𝜏 (𝑦) − ⟨ ℎ̂, 𝑦⟩, Lemma 21 suggests that we can bound the TV error of Algorithm 2 by controlling the
e − 𝜓(𝑉)|, e where 𝑈 e, 𝑉 e ∼ 𝜋𝑆,𝜏 , sampling error of Algorithm 3 and proving a subexponential tail bound for |𝜓(𝑈) or equivalently, for |𝜓(𝑌) − 𝔼[𝜓(𝑌)]| under 𝑌 ∼ 𝜋𝑆,𝜏 . Now, by Theorem 8 and Theorem 18, 𝜋𝑆,𝜏 satisfies 4/3-ATE with high probability. Thus, Theorem 7 gives us the following concentration bound: iid
ℙ𝑌∼𝜋𝑆,𝜏 |𝜓(𝑌) − 𝔼𝜋𝑆,𝜏 [𝜓(𝑌)]| ≤ 2 exp(−𝑐 min{
𝑡2 𝑡 , }) 2 𝔼𝜋𝑆,𝜏 [∥∇𝜓∥ ] max ∥∇2 𝜓(𝑦)∥𝐹
(97)
𝑦∈{±1}𝑆
Since ∥∇2 𝜓(𝑦)∥2𝐹 = ∥∇2 𝑅 𝑆,𝜏 (𝑦)∥𝐹 ≲ ℭ(𝛽) · 𝑛𝑠 2 holds whp uniformly for any 𝑆, 𝜏 and 𝑦 by Theorem 19. Hence, setting 𝑠 = Θ(𝑛 2/3 ) and ℭ(𝛽) sufficiently small, max ∥∇2 𝜓(𝑦)∥𝐹 ≤ 𝛿 for any 𝛿 = 𝑂(1). To control 𝔼𝜋𝑆,𝜏 [∥∇𝜓∥2 ], 3
𝑦∈{±1}𝑆
we note that:
𝔼𝜋𝑆,𝜏 [∥∇𝜓∥2 ] = ∥ ℎ̂ − 𝔼𝜋𝑆,𝜏 [∇𝑅 𝑆,𝜏 ] ∥2 + 𝔼𝜋𝑆,𝜏 ∥ ∇𝑅 𝑆,𝜏 − 𝔼[∇𝑅 𝑆,𝜏 ] ∥2
(98)
Using the Poincare Inequality for 𝜋𝑆,𝜏 and Theorem 19, the second term in Eq. (98) can be made 𝑂(𝛿2 ) for any 𝛿 = 𝑂(1) by setting 𝑠 = Θ(𝑛 2/3 ) and ℭ(𝛽) sufficiently small. The centering field ℎ̂ must be chosen so as to make the first term ∥ ℎ̂ − 𝔼𝜋𝑆,𝜏 [∇𝑅 𝑆,𝜏 ]∥2 as small as possible. Since 𝜋𝑆,𝜏 = 𝜈𝐴𝑆,𝜏 , 𝑏𝑆,𝜏 + ℎ̂ where 𝜈𝐽,𝑣 (𝑦) ∝ exp( 12 ⟨𝑦, 𝐽 𝑦⟩ + ⟨𝑣, 𝑦⟩), the optimal choice of ℎ̂ is given by the solution to the following fixed point problem: ℎ = 𝔼𝜈𝐴𝑆,𝜏 , ℎ+𝑏𝑆,𝜏 [∇𝑅 𝑆,𝜏 ].
(99)
As we shall demonstrate, Eq. (99) admits a unique solution which can be efficiently approximated by an inexact fixed point iteration, namely Algorithm 1, which replaces the expectation in Eq. (99) by an empirical average over approximate samples drawn via RGD. In particular, we prove the following guarantee in Section 3.5 Theorem 22 (Analysis of Algorithm 1). Conditioned on the event ℰ, the output ℎ̂ of Algorithm 1 satisfies ∥ ℎ̂ − 𝔼𝜋𝑆,𝜏 [∇𝑅 𝑆,𝜏 ]∥ ≤ 𝜂 with probability at least 1 − 𝛿. Moreover, Algorithm 1 runs in 𝑂(polylog(𝑛, 𝜂−1 , 𝛿−1 )) parallel time with 𝑂 poly( 𝑛𝜂 , ln(1/𝛿)) work.
21
3.3.1
Proof of Theorem 20
Proof. Recall that 𝜋𝑆,𝜏 (𝑦) ∝ exp 12 ⟨𝑦, 𝐴𝑆,𝜏 𝑦⟩ + ⟨𝑏 𝑆,𝜏 + ℎ̂, 𝑦⟩ where ℎ̂ is the output of Algorithm 1. Let
eℓ , 𝑉 eℓ )ℓ ∈[𝐿] ∼ 𝜋𝑆,𝜏 . Then, by Theorem 17, setting 𝜀RGD = Θ(𝜀step/𝐿) suffices to ensure the following: (𝑈 iid
eℓ , 𝑉 eℓ ) ∀ℓ ≤ 𝐿 ≥ 1 − 𝜀step . ℙ (𝑈ℓ , 𝑉ℓ ) = (𝑈 100
(100)
Henceforth, condition on the above event, and treat (𝑈ℓ , 𝑉ℓ ) as if they were exact i.i.d. samples from 𝜋𝑆,𝜏 . Now, define 𝜓(𝑦) = 𝑅 𝑆,𝜏 (𝑦) − ⟨ ℎ̂, 𝑦⟩. As discussed in Section 3.3, since we condition on the event ℰ that the bounds in Theorem 18 and Theorem 19 hold, the following holds for some universal constant Δ = Θ(1) by setting 𝑠 = (𝑛Δ)2/3 ln(1/𝜀step )−2/3 and ℭ(𝛽) = 𝑂(1) sufficiently small. max ∥∇2 𝜓(𝑦)∥2𝐹 = max ∥∇2 𝑅 𝑆,𝜏 (𝑦)∥2𝐹 ≤ 𝐶1
𝑦∈{±1}𝑠
𝑦∈{±1}𝑠
ℭ(𝛽)2 𝑠 3 Δ2 ≤ 𝑛2 ln(1/𝜀step )2
(101)
By Theorem 22, setting 𝜂 ≤ Δ/8 and 𝛿 = 𝜀step/100,
h ℙ ℎ̂ − 𝔼𝜋𝑆,𝜏 [∇𝑅 𝑆,𝜏 ] ≤
i 𝜀step Δ ≥1− . 100 8 ln(1/𝜀step )
(102)
Henceforth, we condition on the above event as well. Now,
h i 2 2 𝔼𝜋𝑆,𝜏 ∥∇𝜓∥2 = ℎ̂ − 𝔼𝜋𝑆,𝜏 [∇𝑅 𝑆,𝜏 ] + 𝔼 𝑦∼𝜋𝑆,𝜏 ∇𝑅 𝑆,𝜏 (𝑦) − 𝔼𝜋𝑆,𝜏 [∇𝑅 𝑆,𝜏 ] Õ Δ2 ≤
64 ln(1/𝜀step )2
(103)
Var𝜋𝑆,𝜏 [𝜕𝑖 𝑅 𝑆,𝜏 ].
+
(104)
𝑖∈𝑆
By Theorem 18 and Theorem 8, 𝜋𝑆,𝜏 satisfies ATE (and hence, PI) with constant 4/3. Hence,
Õ
Var𝜋𝑆,𝜏 [𝜕𝑖 𝑅 𝑆,𝜏 ] ≤ 34
𝑖∈𝑆
Õ 𝑖∈𝑆
2
Hence, 𝔼𝜋𝑆,𝜏 ∥∇𝜓∥2 ≤ ln(12Δ /𝜀
2 step )
2 𝔼𝜋𝑆,𝜏 ∥∇𝜕𝑖 𝑅 𝑆,𝜏 ∥2 = 43 𝔼𝜋𝑆,𝜏 ∥∇2 𝑅 𝑆,𝜏 ∥2𝐹 ≤ 3 ln(4Δ 1/𝜀
step )
2
(105)
. Since 𝜋𝑆,𝜏 satisfies ATE with constant 4/3, by Theorem 7,
h i 𝑡 2 ln(1/𝜀step )2 𝑡 ln(1/𝜀step ) ℙ𝜋𝑆,𝜏 𝜓(𝑌) − 𝔼𝜋𝑆,𝜏 [𝜓] ≥ 𝑡 ≤ 2 exp −𝐶 min . , 2 Δ
Δ
(106)
Taking a union bound and choosing Δ = Θ(1) small enough, we obtain:
h i −2𝑡 ln(1/𝜀step ) e e . ℙ𝑈e ,𝑉∼𝜋 e 𝑆,𝜏 𝜓(𝑈) − 𝜓(𝑉) ≥ 𝑡 ≤ 4𝑒
(107)
Now consider the parallel rejection sampling step in Algorithm 2. By Lemma 21, setting 𝐿 = Θ(𝑐¯ ln(1/𝜀step )), we conclude the following:
𝐿 𝜀step . ℙ[at least one of the 𝐿 parallel trials is accepted] ≥ 1 − 1 − 21𝑐¯ ≥ 1 − 100
(108)
Henceforth, we condition on the event above. Now, let e 𝜇𝑆,𝜏 denote the conditional law of the first accepted sample among the 𝐿 parallel rejection sampling trials. Then, by Lemma 21, setting 𝑐¯ = Θ(1) yields the following: 𝜀step
𝑑TV (e 𝜇𝑆,𝜏 , 𝜇𝑆,𝜏 ) ≤ 4𝑐¯− ln( /𝜀step ) ≤ 100 1
22
(109)
To conclude, let 𝜇ˆ 𝑆,𝜏 denote the law of the output of Algorithm 2. From Eq. (109), Eq. (108), Eq. (102) and Eq. (100), we conclude the following after taking the appropriate union bounds. 𝑑TV 𝜇ˆ 𝑆,𝜏 , 𝜇𝑆,𝜏 ≤ 𝜀step .
(110)
It remains to analyze the parallel runtime of Algorithm 2. Note that Algorithm 2 involves one call to Algorithm 1, which has parallel runtime 𝑂(polylog(𝑛/𝜀step )) with 𝑂(poly(𝑛, ln(1/𝜀step ))) work as per Theorem 22; and 𝐿 = Θ(ln(1/𝜀step )) parallel calls to Algorithm 3, each of which requires polylog(𝑛/𝜀step ) parallel time with poly(𝑛, ln(1/𝜀step )) work as per Theorem 17. Thus, Algorithm 2 exhibits a parallel runtime 𝑂(polylog(𝑛/𝜀step )) using poly(𝑛, ln(1/𝜀step )) work.
3.4
Proof of Theorem 15
Proof. Let 𝜇𝑡 = Law(𝑋 (𝑡) ). Since 𝜇0 (𝑥) ∝ exp(⟨ℎ, 𝑥⟩), by Lemma 5.5 of [Lee23],
ℙ 𝐷KL 𝜇0 || 𝜇 ≤ 𝐶𝑛 ≥ 1 − exp(−𝑐𝑛).
(111)
Let 𝒢 = ℰ ∩ 𝐷KL 𝜇0 || 𝜇 ≤ 𝐶𝑛 . Then, ℙ[𝒢] ≥ 1 − 𝑒 −𝑐𝑛 , and henceforth, we condition on 𝒢.
e𝑠 denote the Markov kernel for 𝑠-Glauber dynamics and the Markov chain in Algorithm 1 Let 𝑃𝑠 and 𝑃 respectively. By Theorem 20, for 𝑠 = Θ(( log(𝑛𝑛/𝜀) )2/3 ), e𝑠 (𝑥, ·) ≤ 𝜀step = 𝜀 . max 𝑑TV 𝑃𝑠 (𝑥, ·), 𝑃 10𝑇
(112)
𝜀 𝑑TV 𝜇𝑇 , 𝜇0 𝑃𝑠𝑇 ≤ 𝑇𝜀step ≤ 10 .
(113)
𝑥∈{±1}𝑛
Hence,
By Theorem 16, the following holds with probability 1 − 𝑒 −𝑐𝑛 : 𝐷KL 𝜇0 𝑃𝑠𝑇 || 𝜇
≤ exp − 𝑐𝑠𝑇 𝑛
𝑐𝑇 𝐶𝑛 ≤ exp − 1/3 𝑛 log(𝑛/𝜀)2/3
𝐶𝑛.
(114)
Setting 𝑇 = Θ(𝑛 1/3 log(𝑛/𝜀)5/3 ), applying Pinsker’s inequality and taking appropriate union bounds, we conclude that the following holds with probability at least 1 − exp(−𝑐𝑛): 𝜀 𝑑TV (𝜇𝑇 , 𝜇) ≤ 𝑑TV 𝜇𝑇 , 𝜇0 𝑃𝑠𝑇 + 𝑑TV 𝜇0 𝑃𝑠𝑇 , 𝜇 ≤ 10 + 2𝜀 ≤ 𝜀.
(115)
To compute the parallel runtime, note that Algorithm 1 makes 𝑇 calls to Algorithm 2 each of which has parallel runtime 𝑂(polylog(𝑛/𝜀)) with 𝑂(poly(𝑛/𝜀)) processors as per Theorem 20. Hence, Algorithm 1 exhibits a parallel runtime of 𝑂(𝑛 1/3 polylog(𝑛/𝜀)) with 𝑂(poly(𝑛/𝜀)) processors.
3.5
Proof of Theorem 22
For any choice of (𝑆, 𝜏), define the measure 𝛾𝑧 on {±1}𝑠 as 𝛾𝑧 (𝑦) ∝ exp 21 ⟨𝑦, 𝐴𝑆,𝜏 𝑦⟩ + ⟨𝑏 𝑆,𝜏 + 𝑧, 𝑦⟩ . Since we condition on the event ℰ, ∥𝐴𝑆,𝜏 ∥op ≤ 14 . Thus, by Theorem 8, 𝛾𝑧 satisfies LSI with constant 4/3 for any 𝑧 ∈ ℝ𝑠 . 1 Moreover, setting 𝑠 = 𝑂(𝑛 2/3 ) and ℭ(𝛽) = 𝑂(1) sufficiently small, one can set sup 𝑦∈{±1}𝑠 ∥∇2 𝑅 𝑆,𝜏 (𝑦)∥𝐹 ≤ 100 . 𝑠 𝑠 Now, define 𝐹 : ℝ → ℝ as follows:
Í 𝐹(𝑧) = 𝔼𝛾𝑧 ∇𝑅 𝑆,𝜏 (𝑦) =
1 𝑦∈{±1}𝑠 ∇𝑅 𝑆,𝜏 (𝑦) exp 2 ⟨𝑦, 𝐴𝑆,𝜏 𝑦⟩ + ⟨𝑏 𝑆,𝜏 + 𝑧, 𝑦⟩ Í 1 𝑦∈{±1}𝑠 exp 2 ⟨𝑦, 𝐴𝑆,𝜏 𝑦⟩ + ⟨𝑏 𝑆,𝜏 + 𝑧, 𝑦⟩
.
(116)
Taking derivatives with respect to 𝑧 and rearranging terms, ∇𝐹(𝑧) = 𝔼𝛾𝑧 ∇𝑅 𝑆,𝜏 (𝑦)𝑦 ⊤ − 𝔼𝛾𝑧 ∇𝑅 𝑆,𝜏 (𝑦) 𝔼𝛾𝑧 𝑦 ⊤ .
23
(117)
Consider any unit vectors 𝑢, 𝑣 ∈ ℝ𝑠 . Then, by Cauchy–Schwarz and the Poincaré inequality for 𝛾𝑧 , 𝑢 ⊤ ∇𝐹(𝑧)𝑣 = Cov𝛾𝑧 ⟨𝑢, ∇𝑅 𝑆,𝜏 ⟩, ⟨𝑣, 𝑦⟩ ≤
q
Var𝛾𝑧 ⟨𝑢, ∇𝑅 𝑆,𝜏 ⟩ Var𝛾𝑧 ⟨𝑣, 𝑦⟩ ≤
q
1 4 𝔼𝛾𝑧 ⟨𝑢, ∇2 𝑅 𝑆,𝜏 (𝑦)𝑢⟩ < . 3 4 (118)
Since ∥∇𝐹(𝑧)∥op < 14 , by the Banach fixed point theorem, there exists a unique ℎ ∗ satisfying the following fixed point iteration: ℎ ∗ = 𝐹(ℎ ∗ ) = 𝔼𝛾ℎ∗ ∇𝑅 𝑆,𝜏 (𝑦) .
Recall that for 𝑡 ∈ {0, . . . , 𝑇 − 1}, 𝑧 𝑡+1 = 𝐾1
(119)
Í𝐾
(𝑡) (𝑡) 𝛿 , where 𝑈 𝑖 = RGD 𝐴𝑆,𝜏 , 𝑏 𝑆,𝜏 + 𝑧 𝑡 , 100𝐾𝑇 𝑖=1 ∇𝑅 𝑆,𝜏 𝑈 𝑖
. Thus,
(𝑡) iid by Theorem 17, there exist 𝑌𝑖 ∼ 𝛾𝑧𝑡 such that the following holds with probability at least 1 − 𝛿: (𝑡)
(𝑡)
𝑈 𝑖 = 𝑌𝑖
∀ 𝑡 < 𝑇, 𝑖 ∈ [𝐾].
(120)
1 Henceforth, we condition on this event, so that 𝑧 𝑡+1 = 𝐾1 𝐾𝑖=1 ∇𝑅 𝑆,𝜏 𝑌𝑖 . Since ∥∇2 𝑅 𝑆,𝜏 ∥𝐹 ≤ 100 and 𝛾𝑧𝑡 satisfies LSI with constant 4/3, the following holds by Theorem C and a union bound over 𝑗 ∈ 𝑆: (𝑡)
Í
ℙ
𝐾 h 1Õ
𝐾
(𝑡)
∇𝑅 𝑆,𝜏 𝑌𝑖
− 𝔼𝛾𝑧𝑡 ∇𝑅 𝑆,𝜏
i
≥ 𝑢 ≤ 2𝑠𝑒 −𝑐𝐾𝑢 /𝑠 .
2
(121)
𝑖=1
Setting 𝑢 = 16√𝑠 , 𝐾 = Θ 𝑠𝜂−2 ln(𝑠𝑇/𝛿) , and taking a union bound over 𝑡 < 𝑇, we conclude that the following holds with probability at least 1 − 𝛿/2: 3𝜂
∥𝑧 𝑡+1 − 𝐹(𝑧 𝑡 )∥ ≤
3𝜂 16
∀ 𝑡 < 𝑇.
(122)
Let Δ𝑡 = ∥𝑧 𝑡 − ℎ ∗ ∥. Since 𝐹(ℎ ∗ ) = ℎ ∗ and 𝐹 is 1/4-Lipschitz, Δ𝑡+1 = ∥𝑧 𝑡+1 − ℎ ∗ ∥ = ∥𝑧 𝑡+1 − 𝐹(𝑧 𝑡 ) + 𝐹(𝑧 𝑡 ) − 𝐹(ℎ ∗ )∥ ≤
3𝜂 1 Δ𝑡 + . 4 16
(123)
Unrolling the recurrence, Δ𝑇 ≤ 4−𝑇 Δ0 +
𝜂 . 4
(124)
Since 𝑧 0 = 0 and ∥∇𝑅 𝑆,𝜏 ∥∞ ≤ 𝑛 Θ(1) when conditioned on the event ℰ, ∥Δ0 ∥ = ∥ℎ ∗ ∥ = ∥𝐹(ℎ ∗ )∥ = 𝔼𝛾ℎ∗ ∇𝑅 𝑆,𝜏
≤ 𝑛 Θ(1) .
(125)
𝜂
𝜂
= ∥ ℎ̂ − 𝐹( ℎ̂)∥ ≤ ∥ ℎ̂ − ℎ ∗ ∥ + ∥𝐹( ℎ̂) − 𝐹(ℎ ∗ )∥ ≤ 𝜂.
Setting 𝑇 = Θ ln(𝑛/𝜂) suffices to ensure Δ𝑇 ≤ 2 .. Finally, since ℎ̂ = 𝑧𝑇 , we have ∥ ℎ̂ − ℎ ∗ ∥ ≤ 2 which then implies the following:
ℎ̂ − 𝔼𝜋𝑆,𝜏 ∇𝑅 𝑆,𝜏
= ℎ̂ − 𝔼𝛾ℎ̂ ∇𝑅 𝑆,𝜏
(126)
To analyze the parallel runtime and work, note that Algorithm 1 runs for 𝑇 iterations and each iteration 𝛿 performs 𝐾 parallel calls to Algorithm 3 with error tolerance 100𝐾𝑇 . Hence, by Theorem 17, Algorithm 1 has a 𝑛 −1 −1 1 parallel runtime of 𝑂(polylog(𝑛, 𝜂 , 𝛿 )) with 𝑂 poly( 𝜂 , ln( /𝛿)) work.
24
Low Accuracy Parallel Sampler for the 𝑝-Spin Model
4
In this section, we prove that Picard Algorithmic Stochastic Localization (Algorithm 4) samples from the mixed 𝑝-spin model with Wasserstein guarantees. We first illustrate our results for the SK model, a special case of the 𝑝-spin model, because the analysis is simpler and the calculations are closely related. Moreover, for the SK model, we obtain a parallelization result for a broader range of temperatures than the best previous parallel algorithm of Chen, Liu, Yin, and Zhang [Che+25]. The TAP-AMP mean-estimation algorithm for the SK model (Algorithm 6) is slightly different from that for the mixed 𝑝-spin model, but it uses the same AMP iterations followed by NGD iterations. We prove that the TAP-AMP algorithms proposed in prior work for computing mean approximations in the SK and mixed 𝑝-spin models, Algorithm 6 and Algorithm 5, are Lipschitz with respect to the external field and e parallel computable. This enables us to apply Picard iteration to the steps of algorithmic stochastic 𝑂(1) localization Eq. (8), yielding a parallel speedup. For the rest of the section, we first present an improved parameter dependence for algorithmic stochastic localization in Section 4.1. Then, using the improved parameters, we prove the parallelization result for algorithmic stochastic localization in Section 4.2, relying on the following two results:
e parallel time (Section 4.3). • The TAP-AMP mean approximations are executable in 𝑂(1) • The TAP-AMP mean approximations are Lipschitz (Section 4.4).
4.1
Sign Rounding for Algorithmic Stochastic Localization
El Alaoui, Montanari, and Sellke [EMS22; EMS25] provide Wasserstein guarantees for algorithmic stochastic localization using TAP fixed points as mean approximations. Their analysis controls the propagation of the mean-estimation error along the discretized stochastic localization trajectory. They obtain an 𝜀 upper bound in normalized Wasserstein distance between the output of the sampling algorithm and the target distribution by using exp(poly(1/𝜀)) discrete steps to run the stochastic localization up to a time horizon poly(1/𝜀). Proposition 23. For 𝜀𝑛 < 𝜀, El Alaoui, Montanari, and Sellke [EMS22; EMS25] provide algorithmic stochastic localization for a time horizon poly(1/𝜀) with exp(poly(1/𝜀)) discretization steps. Proof. We defer the proof to the appendix, see Section A.4. In this work, in addition to the Picard parallelization and Lipschitz analysis, we improve the dependence of the discretization step size and the time horizon in stochastic localization on the parameter 𝜀 by analyzing sign rounding for algorithmic stochastic localization. Theorem 24. For any 𝑛, there exists a value 𝜀𝑛 such that, for any 𝜀 > 𝜀𝑛 and any high-temperature mixed 𝑝-spin model with Hamiltonian 𝐻𝑛 and no external field, there is an algorithm based on algorithmic stochastic localization using TAP-AMP mean approximations (Algorithm 6 and Algorithm 5) whose output is within 𝜀 of the target distribution in normalized Wasserstein distance with probability 1 − 𝑜(1) over the disorder. Moreover, the algorithmic stochastic localization procedure uses a discretization step of size 𝛿 and time horizon 𝑡 satisfying 𝛿 = poly(𝜀), 𝑡 = 𝑂(log(1/𝜀)). The output of the algorithm is generated by sign rounding, namely sign( 𝑦ˆ𝑡 ), where 𝑦ˆ𝑡 is the output of the discrete algorithmic stochastic localization process at time 𝑡. Furthermore, if 𝑦𝑡 denotes the true stochastic localization process, then 𝑊2,𝑛 (𝑦𝑡 , 𝑦ˆ𝑡 ) ≤ 𝑂(𝜀).
25
Proof. We sketch the proof here; see Section A.1 for the full proof. Let 𝑦𝑡 denote the true stochastic localization process at time 𝑡. Then 𝑦𝑡 satisfies the distributional identity 𝑦𝑡 /𝑡 ∼ 𝑥 ∗ + 𝐵𝑡 /𝑡, where 𝑥 ∗ ∼ 𝜇 and 𝐵𝑡 /𝑡 ∼ 𝒩 (0, 𝐼/𝑡) due to the characterization by Theorem 9. Consider the sign-rounding procedure applied to 𝑦𝑡 . Due to Lemma 10, we have 𝑊2,𝑛 (sign(𝑦𝑡 ), 𝜇) ≤ 𝑂(exp(−𝑡/2)). Therefore, choosing 𝑡 = 𝑂(log(1/𝜀)) makes the rounding error of the true stochastic localization trajectory at most poly(𝜀). Using the error-propagation analysis of [EMS22; EMS25], one can show that the approximate trajectory 𝑦ˆ𝑡 , obtained using TAP fixed points as mean approximations, satisfies 𝑊2,𝑛 (𝑦𝑡 , 𝑦ˆ𝑡 ) ≤ 𝜀, with discretization step size polynomial in 𝜀 and time horizon 𝑡 = 𝑂(log(1/𝜀)). Finally, applying the stability of sign rounding, as in Lemma 11, transfers this trajectory-level approximation to the desired normalized Wasserstein guarantee for the output distribution. Remark 25. In [EMS22; EMS25], the authors prove that there exists a sequence 𝜀𝑛 → 0 such that Algorithmic Stochastic Localization outputs a distribution whose normalized Wasserstein distance from the target measure is at most 𝜀𝑛 . Their guarantee and analysis are asymptotic, and the precise non-asymptotic dependence of 𝜀𝑛 on 𝑛 is left implicit. They show that the error accumulated by the TAP-AMP mean approximation up to a certain time horizon can be controlled. Combining this propagated-error bound with Lemma 11 yields a normalized Wasserstein guarantee. For this work, let 𝜀𝑛 be the smallest value such that
𝑊2,𝑛 𝑦ˆlog(1/𝜀𝑛 ) , 𝑦log(1/𝜀𝑛 ) ≤ 𝜀𝑛 , where 𝑦ˆ𝑡 is generated by the implementation of Algorithmic Stochastic Localization using TAP-AMP mean estimates, with discretization step size 𝛿 = poly(𝜀𝑛 ). Also, take 𝑇𝑛 = 𝑂(log( 𝜀1𝑛 )) as the time horizon where the error between the TAP fixed point and the true mean can be upper bounded by poly(𝜀𝑛 ).
4.2
Parallel Time Analysis of Algorithm 4
In this section, we analyze the parallelization result for algorithmic stochastic localization using Picard iteration. The Picard iterations produce a trajectory that is close, in Wasserstein distance, to the true solution of the discrete stochastic localization process using the TAP-AMP mean approximation. This discrete process is itself close to the true continuous-time stochastic localization process by Theorem 24. Therefore, by the triangle inequality for Wasserstein distance, the input to the rounding step is close to the output of the true stochastic localization process. However, in our proposed algorithm (Algorithm 4), the final output is obtained by applying the sign function. It is not immediate that small Wasserstein distance before rounding implies small Wasserstein distance after rounding. Indeed, if two coordinates are both close to zero but have opposite signs, applying the sign function can amplify their discrepancy. To address this issue, we use the stability result for sign rounding in Lemma 11. This lemma shows that once the Picard trajectory is sufficiently close to the true stochastic localization trajectory, the corresponding sign-rounded outputs remain close in normalized Wasserstein distance.
26
We present our results for the SK model and the mixed 𝑝-spin model separately. For the SK model, the Gibbs measure is defined by 𝛽 𝜇𝐴 (𝑥) ∝ exp ⟨𝑥, 𝐴𝑥⟩ , 2 where 𝐴 ∼ GOE(𝑛). Here GOE(𝑛) denotes the distribution over symmetric matrices whose entries are Gaussian with variance of order 1/𝑛. Although the proof structure is nearly the same, we separate the two cases because, for the SK model, we obtain a parallel sampling algorithm in a temperature regime not covered by previous parallel algorithms. The analysis for the SK model yields the threshold 𝛽 0 ≈ 0.3753 for Lipschitzness and parallelization, improving on the 𝛽 < 1/4 threshold at which the previous polylog(𝑛)-depth algorithm for sampling from the SK model applies [Che+25]. However, unlike the algorithm of Chen, Liu, Yin, and Zhang [Che+25], our low-accuracy algorithm cannot achieve arbitrarily small error. Theorem 26 (SK Model Parallelization). For any 𝜀 ≥ 𝜀𝑛 and inverse temperature 𝛽 < 𝛽0 ≈ 0.3753, there exist parameters (Theorem 24) 𝜂, 𝐾AMP , 𝐾NGD = polylog(𝑛/𝜀), such that the following holds. Let 𝑅 = 𝑂(
𝐾 = poly(1/𝜀),
𝛿 = poly(𝜀),
𝑡 = 𝐾𝛿 = 𝑂(log(1/𝜀))
log2 (1/𝜀) 𝛽0 −𝛽 ), and execute Algorithm 4 on a Hamiltonian of the form
𝐻𝑛 (𝑥) = 𝑥 ⊤ 𝐴𝑥 ˆ with parameters (𝜂, 𝐾AMP , 𝐾NGD , 𝐾, 𝛿). using TAP-AMP (Algorithm 6) as the approximate mean function 𝑚 Then Algorithm 4 outputs a random point xalg ∈ {−1, +1}𝑛 with law 𝜇A such that, with probability 1 − 𝑜(1) over A ∼ GOE(𝑛), alg 𝑊2,𝑛 (𝜇A , 𝜇A ) ≤ 𝑂(𝜀). alg
The parallel runtime of the algorithm is poly(log(𝑛/𝜀)), and its total work is poly(𝑛/𝜀). Proof. First, we argue that with the suitable choice of parameters, the output of Algorithm 4 is close in Wasserstein distance to the output of the sequential algorithmic stochastic localization process Eq. (5). Then we prove that with the parameter choices, the algorithm runs in poly(log(𝑛/𝜀)) parallel time. Since the Brownian motions are sampled and fixed throughout the Picard iteration, Algorithm 4 can be viewed as a differential equation of the form in Eq. (32), where the drift function is the mean-computation ˆ function 𝑚. ˆ given by Algorithm 6, is 𝑂(1/(𝛽 0 − 𝛽))-Lipschitz By Corollary 33, the approximate mean TAP-AMP, i.e. 𝑚, with respect to the tilt, with high probability over the disorder. Hence, we can apply the Picard convergence theorem Lemma 12 to obtain fast convergence of the Picard iteration. (0)
The Picard iteration is initialized with 𝑦ˆ 𝑖𝛿 = 0 for 0 ≤ 𝑖 ≤ 𝐾. Let (1)
∗ 𝑀 = max 𝑦ˆ 𝑖𝛿 − 𝑦ˆ 𝑖𝛿 , 0≤𝑖≤𝐾
∗ where 𝑦ˆ 𝑖𝛿 denotes the true sequential solution of the discrete stochastic localization equation at time 𝑖𝛿. After one Picard iteration, 𝑖−1 Õ √ (1) ∗ ˆ𝑦 𝑖𝛿 − 𝑦ˆ 𝑖𝛿 ˆ 𝑛 , 0) − 𝑚(𝐻 ˆ 𝑛 , 𝑦ˆℓ∗𝛿 ) ≤ 𝛿𝑖 𝑛. = 𝛿 𝑚(𝐻 ℓ =0
Since 𝑖𝛿 ≤ 𝐾𝛿 = 𝑡 = 𝑂(log(1/𝜀)), this gives the uniform bound √ 𝑀 ≤ 𝑂 𝑛 log(1/𝜀) .
27
By the exponential convergence of the Picard iterations Lemma 12, we have
√
(𝑅)
∗ 𝑦ˆ 𝐾𝛿 − 𝑦ˆ 𝐾𝛿 ≤
𝑛
𝐾𝛿 𝛽0 −𝛽
𝑅
𝑅!
.
Since 𝑡 = 𝐾𝛿 = 𝑂(log(1/𝜀)), taking 𝐾𝛿 log(1/𝜀) log2 (1/𝜀) 𝑅≥Ω =Ω 𝛽0 − 𝛽 𝛽0 − 𝛽
!
implies (𝑅)
∗ 𝑦ˆ 𝐾𝛿 − 𝑦ˆ 𝐾𝛿 ≤
√
𝑛 poly(𝜀).
Consequently,
(𝑅)
∗ 𝑊2,𝑛 𝑦ˆ 𝐾𝛿 /𝑡, 𝑦ˆ 𝐾𝛿 /𝑡 ≤ poly(𝜀).
We now use Lemma 11 to bound the normalized ℓ2 Wasserstein distance between the rounded output of the algorithm and the target distribution. Let 𝑦𝑡 denote the true continuous-time stochastic localization process at time 𝑡 = 𝐾𝛿, and let 𝜇𝐴 denote the target distribution. By the triangle inequality,
(𝑅)
(𝑅)
∗ ∗ 𝑊2,𝑛 𝑦ˆ 𝐾𝛿 /𝑡, 𝑦𝑡 /𝑡 ≤ 𝑊2,𝑛 𝑦ˆ 𝐾𝛿 /𝑡, 𝑦ˆ 𝐾𝛿 /𝑡 + 𝑊2,𝑛 𝑦ˆ 𝐾𝛿 /𝑡, 𝑦𝑡 /𝑡 ≤ 𝜀 + poly(𝜀).
∗ Here, 𝑊2,𝑛 𝑦ˆ 𝐾𝛿 /𝑡, 𝑦𝑡 /𝑡 ≤ 𝜀 holds with probability 1 − 𝑜(1) according to Theorem 24. Applying Lemma 11 with (𝑅) Δ = 𝑊2,𝑛 𝑦ˆ 𝐾𝛿 /𝑡, 𝑦𝑡 /𝑡 ≤ 𝜀 + poly(𝜀),
we obtain (𝑅)
𝑊2,𝑛 sign( 𝑦ˆ 𝐾𝛿 ), 𝜇𝑨 ≤ 4(𝜀 + poly(𝜀)) + 𝑒 −Ω(𝑡) = 4𝜀 + poly(𝜀) = 𝑂(𝜀).
Thus, the normalized Wasserstein distance between the output of the parallel algorithm and the target distribution is bounded by 𝑂(𝜀). The stated parallel runtime and total work follow from the parameter choices above and the runtime of the parallel Picard implementation, given the poly(log(𝑛/𝜀)) parallel runtime for the TAP-AMP mean-approximation algorithm (Proposition 28). Theorem 27 (𝑝-spin Model Parallelization). For any 𝜀 ≥ 𝜀𝑛 and temperature coefficients {𝛽 𝑝 }𝑝≥2 satisfying ℭ(𝛽) ≔
𝑃 Õ
q
𝛽 𝑝 𝑝 3 ln(𝑝) < 𝛾0
𝔇(𝛽) =
𝑝=2
𝑃 Õ
q
𝛽 𝑝 2𝑝 𝑝 3 ln(𝑝) < ∞
𝑝=2
there exist parameters (Theorem 24) 𝜂, 𝐾AMP , 𝐾NGD , 𝐾 = poly(1/𝜀),
𝛿 = poly(𝜀),
𝑡 = 𝐾𝛿 = 𝑂(log(1/𝜀)),
such that the following holds. Let 𝑅 = 𝑂ℭ(𝛽) log2 (1/𝜀) be the number of Picard iterations. Given query access to the 𝑝-spin Hamiltonian as defined in Eq. (2), Algorithm 4, run with parameters (𝜂, 𝐾AMP , 𝐾NGD , 𝐾, 𝛿), outputs a random alg point xalg ∈ {−1, +1}𝑛 with law 𝜇𝐻𝑛 such that, with probability 1 − 𝑜(1) over the Gaussian coefficients {𝑔𝐽 }𝐽⊆[𝑛] ,
𝑊2,𝑛 (𝜇𝐻𝑛 , 𝜇𝐻𝑛 ) ≤ 𝑂(𝜀). alg
Here, 𝜀 is the error incurred by the sequential algorithm in Theorem 24, while poly(𝜀) is the additional error incurred by replacing the sequential updates with Picard iteration. The parallel runtime of the algorithm is poly(log(𝑛/𝜀)) and its total work is poly(𝑛, 1/𝜀).
28
Proof. The proof follows the same argument as the SK case Theorem 26; we give only the brief sketch. Fix the Brownian increments used by Algorithm 4. With these increments fixed, the algorithm can be viewed ˆ Since the temperature as a discrete differential equation whose drift is given by the approximate mean map 𝑚. condition holds, Corollary 43 implies that the TAP-AMP mean approximation is Lipschitz with respect to the external field with high probability over the disorder. Then Lemma 12, with the same initialization as in the SK case, √ implies that the Picard trajectory contracts to the discrete stochastic localization trajectory up to error 𝑛 poly(𝜀). Hence, after normalization, (𝑅)
∗ 𝑊2,𝑛 ( 𝑦ˆ 𝐾𝛿 /𝑡, 𝑦ˆ 𝐾𝛿 /𝑡) ≤ poly(𝜀).
The sequential stochastic localization guarantee from Theorem 24 contributes the 𝑂(𝜀) error. Combining this with the Picard error, the triangle inequality gives (𝑅)
(𝑅)
𝑊2,𝑛 ( 𝑦ˆ 𝐾𝛿 /𝑡, 𝑦𝑡 /𝑡) ≤ 𝑊2,𝑛 ( 𝑦ˆ𝑡 /𝑡, 𝑦ˆ𝑡∗ /𝑡) + 𝑊2,𝑛 ( 𝑦ˆ𝑡∗ /𝑡, 𝑦𝑡 /𝑡) ≤ 𝑂(𝜀) + poly(𝜀). Finally, applying the stability rounding lemma (Lemma 11), we have 𝑊2,𝑛 (𝜇𝐻𝑛 , 𝜇𝐻𝑛 ) ≤ 𝑂(𝜀) + poly(𝜀) = 𝑂(𝜀). alg
The stated parallel runtime and total work follow from the parameter choices above and the runtime of the parallel Picard implementation, given the poly(log(𝑛/𝜀)) parallel runtime for the TAP-AMP meanapproximation algorithm (Proposition 28).
4.3
TAP-AMP in polylog(𝑛) Parallel Time
The TAP-AMP mean-approximations can be computed in polylog(𝑛/𝜀) parallel depth. The TAP-AMP algorithms Algorithm 6 and Algorithm 5 consist of AMP iterations followed by NGD iterations. Each individual AMP or NGD iteration is computable in parallel with depth 𝑂(log 𝑛); thus, it suffices to prove e bounds on the numbers of AMP and NGD iterations. Since this proof is rather involved and does not 𝑂(1) affect the rest of the argument, we give only a brief sketch here and defer the details to the appendix; see Section A.2. Proposition 28. The mean-approximation algorithms Algorithm 6 and Algorithm 5 are computable with parallel depth polylog(𝑛/𝜀), given query access to the mixed 𝑝-spin Hamiltonian and the parameter choices specified in Theorem 24 for attaining normalized Wasserstein error 𝜀. Proof. We give a proof sketch here and defer the full proof to Section A.2. The mean-approximation algorithms consist of two phases: AMP iterations followed by NGD iterations. Each AMP iteration is computable in e parallel depth. Thus, the parallel depth of the AMP phase is controlled by the parameter 𝐾AMP , for 𝑂(1) which we give an explicit bound in the full proof. For the NGD phase, given gradient-query access to the Hamiltonian, the parallel depth is controlled by the number of NGD iterations. In the full proof, we show that this number is at most polylog(𝑛).
4.4
Lipschitz Property of the TAP-AMP Mean Approximation for the SK Model
In this section, we prove that, for the SK model at high temperature 𝛽 < 𝛽0 ≈ 0.3753, the TAP-AMP algorithm
(Algorithm 6) for mean approximation is 𝑂 𝛽01−𝛽 -Lipschitz with respect to the tilt parameter, or external field, with high probability. We provide an analogous result for the mixed 𝑝-spin model in Section A.3.
29
Algorithm 6: Mean of the Tilted Ising Measure [EMS22] Input: Data 𝐴 ∈ ℝ𝑛×𝑛 , 𝑦 ∈ ℝ𝑛 , parameters 𝛽, 𝜂 > 0, 𝑞 ∈ (0, 1), iteration numbers 𝐾 AMP , 𝐾NGD 𝑚−1 = 𝑧0 = 0 for 𝑘 = 0, . . . , 𝐾AMP − 1 do 𝛽2 Í
𝑚 𝑘 = tanh(𝑧 𝑘 ), 𝑏 𝑘 = 𝑛 𝑛𝑖=1 (1 − tanh2 ((𝑧 𝑘 )𝑖 )) 𝑧 𝑘+1 = 𝛽𝐴𝑚 𝑘 + 𝑦 − 𝑏 𝑘 𝑚 𝑘−1 𝑢 0 = 𝑧 𝐾AMP for 𝑘 = 0, . . . , 𝐾NGD − 1 do 𝑢 𝑘+1 = 𝑢 𝑘 − 𝜂 · ∇ℱ̂TAP (𝑚 𝑘+ ; 𝑦, 𝑞) + 𝑚 𝑘+1 = tanh(𝑢 𝑘+1 ) return 𝑚𝐾+NGD Algorithm 6 consists of two phases: it first uses AMP iterations attempting to get close to the fixed point of the TAP iteration, and then it uses an NGD algorithm to close the gap to the fixed point as much as needed. Here, we will take a closer look at the two phases of Algorithm 6. Phase 1: AMP phase.
In the first phase, initialize 𝑚−1 = 0,
𝑧 0 = 0,
and for 𝑘 = 0, 1, . . . , 𝐾AMP − 1 define 𝑚 𝑘 = tanh(𝑧 𝑘 ),
(127)
𝛽2 Õ sech2 (𝑧 𝑘 )𝑖 , 𝑛
(128)
𝑧 𝑘+1 = 𝛽𝐴𝑚 𝑘 + 𝑦 − 𝑏 𝑘 𝑚 𝑘−1 .
(129)
𝑏𝑘 =
𝑖∈[𝑛]
Phase 2: NGD phase on TAP free energy.
Initialize 𝑢0 := 𝑧 𝐾AMP . For 𝑡 = 0, 1, . . . , 𝐾 NGD − 1 define
𝑚𝑡+ = tanh(𝑢𝑡 ),
(130)
𝑢𝑡+1 = 𝑢𝑡 − 𝜂 ∇𝑚 ℱ̂TAP (𝑚𝑡+ ; 𝑦, 𝑞).
(131)
The algorithm outputs ˆ 𝑚(𝐴, 𝑦) := 𝑚𝐾+NGD = tanh(𝑢𝐾NGD ). ˆ (𝐾AMP ,𝐾NGD ) (𝐴, 𝑦) be Definition 29. For a given external field 𝑦 ∈ ℝ𝑛 and interaction matrix 𝐴 ∈ ℝ𝑛×𝑛 , let 𝑚 the output of the algorithm after running the AMP iteration for 𝐾 AMP steps and then running NGD for 𝐾 NGD ˆ steps. We refer to 𝑚(𝑦) as the tilted-mean computation for the external field 𝑦 whenever 𝐴, 𝐾 AMP , and 𝐾 NGD are clear from the context. Our goal in this section is to find the temperature 𝛽(𝐴) = 𝛽(∥𝐴∥op ) such that, for any two tilts 𝑦 and e 𝑦 , the TAP fixed point of the Ising model with interaction matrix 𝐴 and temperature 𝛽 satisfies ˆ (𝐾AMP ,𝐾NGD ) (𝐴, 𝑦) − 𝑚 ˆ (𝐾AMP ,𝐾NGD ) (𝐴, e ∥𝑚 𝑦 )∥ ≤
𝛾(𝐴) ∥𝑦 − e 𝑦 ∥. 𝛽(𝐴) − 𝛽
Note that if the matrix 𝐴 is sampled from the GOE(𝑛) distribution, then its maximum eigenvalue concentrates around 2; thus, if we have inverse temperature 𝛽 satisfying 𝛾(2, 𝛽) < 1, then we have the Lipschitz property ˆ with high probability. for the TAP-AMP mean approximation 𝑚 To prove the Lipschitz property of the mean-approximation algorithm (Algorithm 6), we will use the Lipschitz property of the following two functions:
30
• tanh is 1-Lipschitz (Lemma 13). • sech2 is
4 √ -Lipschitz (Lemma 14). 3 3
We will prove that the AMP and NGD iterations are both Lipschitz. Throughout the rest of this section, we fix two external fields 𝑦 and e 𝑦 and use the following notation, where the unadorned variables correˆ (𝐾AMP ,𝐾NGD ) (𝐴, 𝑦) and the tilded variables correspond to the computation of spond to the computation of 𝑚 ˆ (𝐾AMP ,𝐾NGD ) (𝐴, e 𝑚 𝑦 ): • Δ𝑦 := 𝑦 − e 𝑦 (tilt difference) • Δ𝑧 𝑘 := 𝑧 𝑘 − e 𝑧 𝑘 (AMP iterations)
e 𝑘 (AMP iterations) • Δ𝑚 𝑘 := 𝑚 𝑘 − 𝑚 • Δ𝑢 𝑘 := 𝑢 𝑘 − e 𝑢 𝑘 (NGD iterations)
e−1 = 0 for the initialization. Also, Δ𝑧 0 = Δ𝑧−1 = 0 and 𝑚−1 = 𝑚 Theorem 30 (Lipschitzness of AMP phase). Assume 𝛾(𝐴, 𝛽) := 𝛽 ∥𝐴∥op + (1 + ∥Δ𝑧 𝑘 ∥ ≤
1 ∥Δ𝑦∥, 1−𝛾
∥Δ𝑚 𝑘 ∥ ≤
4 √ )𝛽 2 < 1. Then for all 𝑘 ≥ 0, 3 3
1 ∥Δ𝑦∥. 1−𝛾
Proof. Let 𝑉𝑘 := max0≤𝑗≤𝑘 ∥Δ𝑧 𝑗 ∥. We will prove the following two inequalities for each iteration: • ∥Δ𝑧 𝑘+1 ∥ ≤ (𝛽 ∥𝐴∥op +
4𝛽2 √ )∥Δ𝑧 𝑘 ∥ + 𝛽 2 ∥Δ𝑧 𝑘−1 ∥ + ∥Δ𝑦∥ 3 3
• ∥Δ𝑚 𝑘 ∥ ≤ ∥Δ𝑧 𝑘 ∥ The second item follows immediately from the 1-Lipschitz property in Lemma 13. Thus, we focus on proving the first item. By taking the difference of Eq. (129) for the two tilts 𝑦 and e 𝑦 , we have the following identity:
e 𝑘−1 ∥ ∥Δ𝑧 𝑘+1 ∥ = ∥𝛽𝐴Δ𝑚 𝑘 + Δ𝑦 − 𝑏 𝑘 Δ𝑚 𝑘−1 − (𝑏 𝑘 − e 𝑏 𝑘 )𝑚 e 𝑘−1 ∥ ≤𝛽 ∥𝐴∥op ∥Δ𝑚 𝑘 ∥ + ∥Δ𝑦∥ + |𝑏 𝑘 | · ∥Δ𝑚 𝑘−1 ∥ + |𝑏 𝑘 − e 𝑏 𝑘 | · ∥𝑚
(132)
Now, we upper bound each term in Eq. (132). Since tanh is 1-Lipschitz, we have ∥Δ𝑚 𝑘 ∥ ≤ ∥Δ𝑧 𝑘 ∥ ≤ 𝑉𝑘 and e 𝑘−1 ∥. We upper bound the ∥Δ𝑚 𝑘−1 ∥ ≤ ∥Δ𝑧 𝑘−1 ∥ ≤ 𝑉𝑘 . Also, 0 ≤ 𝑏 𝑘 ≤ 𝛽 2 . The remaining term is |𝑏 𝑘 − e 𝑏 𝑘 | · ∥𝑚 absolute value of the scalar 𝑏 𝑘 − e 𝑏 𝑘 as follows: 𝑛
|𝑏 𝑘 − e 𝑏𝑘 | =
𝛽2 Õ sech2 (𝑧 𝑘 )𝑖 − sech2 (e 𝑧 𝑘 )𝑖 𝑛 𝑖=1
≤
≤
𝑛 𝛽2 Õ
𝑛
𝑖=1 𝑛 2 Õ
𝛽 𝑛
sech2 (𝑧 𝑘 )𝑖 − sech2 (e 𝑧 𝑘 )𝑖
4 𝑧 𝑘 )𝑖 | √ |(𝑧 𝑘 )𝑖 − (e 𝑖=1 3 3
𝛽2 4 ≤ √ · √ ∥Δ𝑧 𝑘 ∥, 𝑛 3 3
(Triangle inequality)
(sech2 is
4 √ -Lipschitz) 3 3
(Cauchy–Schwarz)
√ e 𝑘−1 ∥ ≤ 𝑛, this implies |𝑏 𝑘 − e e 𝑘−1 ∥ ≤ 𝛽2 √4 ∥Δ𝑧 𝑘 ∥ ≤ 𝛽2 √4 𝑉𝑘 . Using the inequalities for the Since ∥𝑚 𝑏 𝑘 | · ∥𝑚 3 3 3 3 terms in Eq. (132), we have 4 ∥Δ𝑧 𝑘+1 ∥ ≤ 𝑉𝑘 𝛽 ∥𝐴∥op + ∥Δ𝑦∥ + 𝛽2𝑉𝑘 + 𝛽 2 √ 𝑉𝑘 3 3 = 𝛾(𝐴, 𝛽)𝑉𝑘 + ∥Δ𝑦∥
31
Hence, 𝑉𝑘+1 ≤ max{𝑉𝑘 , 𝛾𝑉𝑘 + ∥Δ𝑦∥}. As long as 𝛾(𝐴, 𝛽) < 1, the uniform bound ∥Δ𝑦∥/(1 − 𝛾) follows by a simple induction.
Lipschitzness of the full AMP+NGD algorithm for the SK model 1 We have proved so far that after an arbitrary number of iterations, ∥Δ𝑚 𝑘 ∥ is upper bounded by 1−𝛾(𝐴,𝛽) ∥Δ𝑦∥. Now, we prove that applying NGD after AMP preserves the same bound for the difference of the mean estimates. For convenience, and to avoid overloading the notation 𝑚 𝑘 from the AMP phase, we analyze the NGD step in the variables 𝑢𝑖 defined by 𝑚 = tanh(𝑢) (equivalently 𝑢 = atanh(𝑚) component-wise).
Rewriting the gradient of TAP in 𝑢-space. The TAP equation is slightly different for the SK model than for the mixed 𝑝-spin model. For a prescribed overlap parameter 𝑞, the SK TAP free energy satisfies ∇ℱTAP (𝑚; 𝑦, 𝑞) = −𝛽𝐴𝑚 − 𝑦 + arctanh(𝑚) + 𝛽 2 (1 − 𝑞)𝑚.
(133)
Plugging 𝑚 = tanh(𝑢) into (133) gives ∇𝑚 ℱ̂TAP (tanh(𝑢); 𝑦, 𝑞) = −𝛽𝐴 tanh(𝑢) − 𝑦 + 𝑢 + 𝛽2 (1 − 𝑞) tanh(𝑢)
= − 𝛽𝐴 − 𝛽 2 (1 − 𝑞)𝐼 tanh(𝑢) − 𝑦 + 𝑢. Define 𝐵 := 𝛽𝐴 − 𝛽2 (1 − 𝑞)𝐼. Then the NGD update (131) becomes 𝑢𝑡+1 = (1 − 𝜂)𝑢𝑡 + 𝜂 𝐵 tanh(𝑢𝑡 ) + 𝜂 𝑦.
(NGD update)
Theorem 31 (TAP Lipschitz property). Assume 𝛾(𝐴, 𝛽) = 𝛽 ∥𝐴∥op + (1 + √4 )𝛽 2 < 1 and 𝑞 ∈ (0, 1). Let 0 < 𝜂 ≤ 1. 3 3
Set 𝑢0 = 𝑧 𝐾AMP and e 𝑢0 = e 𝑧 𝐾AMP . Then for every 𝑡 ≥ 0, ∥𝑢𝑡 (𝑦) − 𝑢𝑡 (e 𝑦 )∥ ≤
1 ∥𝑦 − e 𝑦 ∥. 1 − 𝛾(𝐴, 𝛽)
Consequently, the final output satisfies ˆ ˆ e ∥𝑚(𝐴, 𝑦) − 𝑚(𝐴, 𝑦 )∥ ≤
1 ∥𝑦 − e 𝑦 ∥. 1 − 𝛾(𝐴, 𝛽)
Proof. Recall that Δ𝑢𝑡 := 𝑢𝑡 (𝑦) − 𝑢𝑡 (e 𝑦 ). By Theorem 30 and initialization 𝑢0 = 𝑧 𝐾AMP , ∥Δ𝑢0 ∥ = ∥Δ𝑧 𝐾AMP ∥ ≤
1 ∥Δ𝑦∥. 1−𝛾
Subtract (NGD update) for 𝑦 and e 𝑦: ∥Δ𝑢𝑡+1 ∥ = ∥(1 − 𝜂)Δ𝑢𝑡 + 𝜂 𝐵 tanh(𝑢𝑡 ) − tanh(e 𝑢𝑡 ) + 𝜂 Δ𝑦∥
≤ (1 − 𝜂)∥Δ𝑢𝑡 ∥ + 𝜂 ∥𝐵∥op ∥tanh(𝑢𝑡 ) − tanh(e 𝑢𝑡 )∥ + 𝜂∥Δ𝑦∥ ≤ (1 − 𝜂)∥Δ𝑢𝑡 ∥ + 𝜂 ∥𝐵∥op ∥Δ𝑢𝑡 ∥ + 𝜂∥Δ𝑦∥
(triangle inequality) (tanh is Lipschitz)
∥Δ𝑦∥
Now, it is enough to prove that if ∥Δ𝑢𝑡 ∥ ≤ 1−𝛾(𝐴,𝛽) , then ∥Δ𝑢𝑡+1 ∥ satisfies the same upper bound. Since 𝐵 = 𝛽𝐴 − 𝛽 2 (1 − 𝑞)𝐼, we have ∥𝐵∥op ≤ 𝛽 ∥𝐴∥op + 𝛽2 (1 − 𝑞) ≤ 𝛽 ∥𝐴∥op + 𝛽2 since 𝑞 ∈ (0, 1). Thus, we need to prove that
1 − 𝜂 + 𝜂(𝛽 ∥𝐴∥op + 𝛽 2 )
1 1 +𝜂 ≤ ⇐⇒ 𝜂(𝛽 ∥𝐴∥op + 𝛽 2 − 𝛾(𝐴, 𝛽)) ≤ 0. 1 − 𝛾(𝐴, 𝛽) 1 − 𝛾(𝐴, 𝛽) 32
Therefore, the inequality holds because 𝛾(𝐴, 𝛽) = 𝛽 ∥𝐴∥op + (1 + √4 )𝛽2 ≥ 𝛽 ∥𝐴∥op + 𝛽 2 . This concludes the 3 3 Lipschitz property of Algorithm 6 given 𝛾(𝐴, 𝛽) < 1. Remark 32. The key stability requirement for AMP and NGD iterations is 4 𝛾(𝐴, 𝛽) = 𝛽 ∥𝐴∥op + (1 + √ )𝛽 2 < 1. 3 3 In random matrix settings (e.g. 𝐴 ∼ GOE(𝑛)), one typically has ∥𝐴∥op = 2 + 𝑜(1) with high probability, so any 1 sufficiently small 𝛽 ensuring 𝛾(2, 𝛽) < 1 implies that the TAP-AMP mean approximation is 1−𝛾(2,𝛽)−𝑐 -Lipschitz with high probability for every sufficiently small constant 𝑐. Corollary 33. Fix 𝜀 > 0. If the SK inverse temperature 𝛽 is at most 𝛽0 − 𝜀, where 𝛽 0 ≈ 0.3753 is the positive root of the equation 4 𝛾(2, 𝛽) = 2𝛽 + (1 + √ )𝛽2 = 1, 3 3 then the tilted-mean computation in Algorithm 6 is 𝑂( 1𝜀 )-Lipschitz for the SK model with high probability. Proof. We know that ∥𝐴∥op is concentrated around 2, and the probability that it exceeds 2 + 𝜀/𝛽 0 is 𝑒 −Ω𝜀 (𝑛) [AGZ10]. Recall that 𝛽 0 is the positive solution of the equation 2𝑥 + (1 + √4 )𝑥 2 = 1. Thus, we only need to show that the parameter 𝛾 = 𝛽 ∥𝐴∥op + (1 +
3 3 4 √ )𝛽2 ≤ 1 − Ω(𝜀) with high probability. 3 3
4 4 𝛽 ∥𝐴∥op + (1 + √ )𝛽2 ≤ 𝛽(2 + 𝜀/𝛽 0 ) + (1 + √ )𝛽2 3 3 3 3 4 ≤ 𝜀 + 2𝛽 + (1 + √ )𝛽2 3 3 4 ≤ −𝜀 + 2𝛽 0 + (1 + √ )𝛽 20 = 1 − 𝜀 3 3 Since P{∥𝐴∥op > 2 + 𝜀/𝛽 0 } ≤ 𝑒 −Ω𝜀 (𝑛) , the tilted-mean computation is 𝑂( 1𝜀 )-Lipschitz with high probability if 𝛽 ≤ 𝛽0 − 𝜀.
References [ACV24]
Nima Anari, Sinho Chewi, and Thuy-Duong Vuong. “Fast parallel sampling under isoperimetry”. In: The Thirty Seventh Annual Conference on Learning Theory. PMLR. 2024, pp. 161–185. [Adh+22] Arka Adhikari, Christian Brennecke, Changji Xu, and Horng-Tzer Yau. “Spectral Gap Estimates for Mixed 𝑝-Spin Models at High Temperature”. In: arXiv preprint arXiv:2208.07844 (2022). [AGJ20] Gerard Ben Arous, Reza Gheissari, and Aukosh Jagannath. “Algorithmic thresholds for tensor PCA”. In: The Annals of Probability 48.4 (2020), pp. 2052–2087. [AGR24] Nima Anari, Ruiquan Gao, and Aviad Rubinstein. “Parallel sampling via counting”. In: Proceedings of the 56th Annual ACM Symposium on Theory of Computing. 2024, pp. 537–548. [AGZ10] Greg W. Anderson, Alice Guionnet, and Ofer Zeitouni. An Introduction to Random Matrices. Vol. 118. Cambridge Studies in Advanced Mathematics. Cambridge University Press, 2010. [AKV24] Nima Anari, Frederic Koehler, and Thuy-Duong Vuong. “Trickle-down in localization schemes and applications”. In: Proceedings of the 56th Annual ACM Symposium on Theory of Computing. 2024, pp. 1094–1105. [Ana+22] Nima Anari, Vishesh Jain, Frederic Koehler, Huy Tuan Pham, and Thuy-Duong Vuong. “Entropic independence: optimal mixing of down-up random walks”. In: Proceedings of the 54th Annual ACM SIGACT Symposium on Theory of Computing. 2022, pp. 1418–1430.
33
[Ana+23]
[Ana+24]
[Ana+25] [BB19] [CE22]
[Cel24] [Che+25] [Csa75] [DFK91]
[Dob68] [EKZ22]
[Eld13] [EM22] [EMS22]
[EMS23] [EMS25] [FHY21]
[FYC23]
[GJ19] [GKK24] [Gla63] [GŠV20]
Nima Anari, Yizhi Huang, Tianyu Liu, Thuy-Duong Vuong, Brian Xu, and Katherine Yu. “Parallel discrete sampling via continuous walks”. In: Proceedings of the 55th Annual ACM Symposium on Theory of Computing. 2023, pp. 103–116. Nima Anari, Vishesh Jain, Frederic Koehler, Huy Tuan Pham, and Thuy-Duong Vuong. “Universality of spectral independence with applications to fast mixing in spin glasses”. In: Proceedings of the 2024 Annual ACM-SIAM Symposium on Discrete Algorithms (SODA). SIAM. 2024, pp. 5029–5056. Nima Anari, Carlo Baronio, CJ Chen, Alireza Haqi, Frederic Koehler, Anqi Li, and Thuy-Duong Vuong. “Parallel Sampling via Autospeculation”. In: arXiv preprint arXiv:2511.07869 (2025). Roland Bauerschmidt and Thierry Bodineau. “A very simple proof of the LSI for high temperature spin systems”. In: Journal of Functional Analysis 276.8 (2019), pp. 2582–2588. Yuansi Chen and Ronen Eldan. “Localization schemes: A framework for proving mixing bounds for Markov chains”. In: 2022 IEEE 63rd Annual Symposium on Foundations of Computer Science (FOCS). IEEE. 2022, pp. 110–122. Michael Celentano. “Sudakov–Fernique post-AMP, and a new proof of the local convexity of the TAP free energy”. In: The Annals of Probability 52.3 (2024), pp. 923–954. Xiaoyu Chen, Hongyang Liu, Yitong Yin, and Xinyuan Zhang. “Efficient Parallel Ising Samplers via Localization Schemes”. In: arXiv preprint arXiv:2505.05185 (2025). Laszlo Csanky. “Fast parallel matrix inversion algorithms”. In: 16th Annual Symposium on Foundations of Computer Science (sfcs 1975). IEEE. 1975, pp. 11–12. Martin Dyer, Alan Frieze, and Ravi Kannan. “A random polynomial-time algorithm for approximating the volume of convex bodies”. In: Journal of the ACM (JACM) 38.1 (1991), pp. 1– 17. PL Dobruschin. “The description of a random field by means of conditional probabilities and conditions of its regularity”. In: Theory of Probability & Its Applications 13.2 (1968), pp. 197–224. Ronen Eldan, Frederic Koehler, and Ofer Zeitouni. “A spectral condition for spectral gap: fast mixing in high-temperature Ising models”. In: Probability theory and related fields 182.3 (2022), pp. 1035–1051. Ronen Eldan. “Thin shell implies spectral gap up to polylog via a stochastic localization scheme”. In: Geometric and Functional Analysis 23.2 (2013), pp. 532–569. Ahmed El Alaoui and Andrea Montanari. “An information-theoretic view of stochastic localization”. In: IEEE Transactions on Information Theory 68.11 (2022), pp. 7423–7426. Ahmed El Alaoui, Andrea Montanari, and Mark Sellke. “Sampling from the SherringtonKirkpatrick Gibbs measure via algorithmic stochastic localization”. In: 2022 IEEE 63rd Annual Symposium on Foundations of Computer Science (FOCS). IEEE. 2022, pp. 323–334. Ahmed El Alaoui, Andrea Montanari, and Mark Sellke. “Algorithmic stochastic localization for sampling from the p-spin model”. In: arXiv preprint arXiv:2307.07130 (2023). Ahmed El Alaoui, Andrea Montanari, and Mark Sellke. “Sampling from mean-field Gibbs measures via diffusion processes”. In: Probability and Mathematical Physics 6.3 (2025), pp. 961–1022. Weiming Feng, Thomas P Hayes, and Yitong Yin. “Distributed metropolis sampler with optimal parallelism”. In: Proceedings of the 2021 ACM-SIAM Symposium on Discrete Algorithms (SODA). SIAM. 2021, pp. 2121–2140. Jiaojiao Fan, Bo Yuan, and Yongxin Chen. “Improved dimension dependence of a proximal algorithm for sampling”. In: The Thirty Sixth Annual Conference on Learning Theory. PMLR. 2023, pp. 1473–1521. Reza Gheissari and Aukosh Jagannath. “On the spectral gap of spherical spin glass dynamics”. In: (2019). Andreas Galanis, Alkis Kalavasis, and Anthimos Vardis Kandiros. “On sampling from Ising models with spectral constraints”. In: arXiv preprint arXiv:2407.07645 (2024). Roy J Glauber. “Time-dependent statistics of the Ising model”. In: Journal of mathematical physics 4.2 (1963), pp. 294–307. Andreas Galanis, Daniel Štefankovič, and Eric Vigoda. “The complexity of approximating averages on bounded-degree graphs”. In: 2020 IEEE 61st Annual Symposium on Foundations of Computer Science (FOCS). IEEE. 2020, pp. 1345–1355.
34
[HMP24]
Brice Huang, Andrea Montanari, and Huy Tuan Pham. “Sampling from spherical spin glasses in total variation via algorithmic stochastic localization”. In: arXiv preprint arXiv:2404.15651 (2024). [Hu+25] Hengyuan Hu, Aniket Das, Dorsa Sadigh, and Nima Anari. “Diffusion Models are Secretly Exchangeable: Parallelizing DDPMs via Autospeculation”. In: arXiv preprint arXiv:2505.03983 (2025). [Hua+25] Brice Huang, Sidhanth Mohanty, Amit Rajaraman, and David X Wu. “Weak Poincaré inequalities, simulated annealing, and sampling from spherical spin glasses”. In: Proceedings of the 57th Annual ACM Symposium on Theory of Computing. 2025, pp. 915–923. Mark Jerrum, Alistair Sinclair, and Eric Vigoda. “A polynomial-time approximation algorithm [JSV04] for the permanent of a matrix with nonnegative entries”. In: Journal of the ACM (JACM) 51.4 (2004), pp. 671–697. [JVV86] Mark R Jerrum, Leslie G Valiant, and Vijay V Vazirani. “Random generation of combinatorial structures from a uniform distribution”. In: Theoretical computer science 43 (1986), pp. 169–188. [Lai+25] Chieh-Hsin Lai, Yang Song, Dongjun Kim, Yuki Mitsufuji, and Stefano Ermon. “The principles of diffusion models”. In: arXiv preprint arXiv:2510.21890 (2025). [Lee23] Holden Lee. “Parallelising glauber dynamics”. In: arXiv preprint arXiv:2307.07131 (2023). [LV24] Yin Tat Lee and Santosh S Vempala. “Eldan’s stochastic localization and the KLS conjecture: Isoperimetry, concentration and mixing”. In: Annals of Mathematics 199.3 (2024), pp. 1043–1092. Hongyang Liu and Yitong Yin. “Simple parallel algorithms for single-site dynamics”. In: [LY22] Proceedings of the 54th Annual ACM SIGACT Symposium on Theory of Computing. 2022, pp. 1431– 1444. [Méz+88] Marc Mézard, Giorgio Parisi, Miguel Angel Virasoro, and David J Thouless. Spin glass theory and beyond. 1988. [MS24] Andrea Montanari and Subhabrata Sen. “A friendly tutorial on mean-field spin glass techniques for non-physicists”. In: Foundations and Trends in Machine Learning 17.1 (2024), pp. 1–173. [MVV87] Ketan Mulmuley, Umesh V Vazirani, and Vijay V Vazirani. “Matching is as easy as matrix inversion”. In: Proceedings of the nineteenth annual ACM symposium on Theory of computing. 1987, pp. 345–354. [Pan12] Dmitry Panchenko. “The Sherrington-Kirkpatrick model: an overview”. In: Journal of Statistical Physics 149.2 (2012), pp. 362–383. [Par81] Giorgio Parisi. “Correlation functions and computer simulations”. In: Nuclear Physics B 180.3 (1981), pp. 378–384. [Per13] Lawrence Perko. Differential equations and dynamical systems. Vol. 7. Springer Science & Business Media, 2013. [Shi+23] Andy Shih, Suneel Belkhale, Stefano Ermon, Dorsa Sadigh, and Nima Anari. “Parallel sampling of diffusion models”. In: Advances in Neural Information Processing Systems 36 (2023), pp. 4263–4276. [SS12] Allan Sly and Nike Sun. “The computational hardness of counting in two-spin models on d-regular graphs”. In: 2012 IEEE 53rd Annual Symposium on Foundations of Computer Science. IEEE. 2012, pp. 361–369. Holger Sambale and Arthur Sinulis. “Modified log-Sobolev inequalities and two-level concentra[SS19] tion”. In: arXiv preprint arXiv:1905.06137 (2019). [Tal10] Michel Talagrand. Mean field models for spin glasses: Volume I: Basic examples. Vol. 54. Springer Science & Business Media, 2010. [TAP77] David J Thouless, Philip W Anderson, and Robert G Palmer. “Solution of’solvable model of a spin glass’”. In: Philosophical Magazine 35.3 (1977), pp. 593–601. [Van14] Ramon Van Handel. Probability in high dimension. Tech. rep. 2014. [Ver18] Roman Vershynin. High-dimensional probability: An introduction with applications in data science. Vol. 47. Cambridge university press, 2018.
35
A.1
Sign Rounding Analysis for Algorithmic Stochastic Localization
In this appendix, we show the improved parameter dependence for Algorithmic Stochastic Localization (ASL) using the notation of Algorithm 4 and the low-accuracy discussion in the introduction. We use 𝐾
for the number of localization steps,
𝛿
for the discretization step size,
so that the localization time horizon is 𝐾𝛿. The accuracy parameter 𝜀 below is the target error in normalized 2-Wasserstein distance 𝑊2,𝑛 from Definition 3. The threshold 𝜀𝑛 = 𝑜 𝑛 (1) is the accuracy threshold of the tilted-mean estimator used in the ASL analysis.
b𝑡 (𝑦) be the TAP-AMP estimate of the tilted mean at localization time 𝑡; the dependence on 𝑡 enters Let 𝑚 through the overlap parameter 𝑞★(𝑡) in Algorithm 5. The sequential ASL Euler trajectory is b b 𝑘𝛿 (b 𝑦(𝑘+1)𝛿 = b 𝑦 𝑘𝛿 + 𝛿 𝑚 𝑦 𝑘𝛿 ) + 𝐵(𝑘+1)𝛿 − 𝐵 𝑘𝛿 ,
b 𝑦0 = 0,
√
or equivalently, with 𝑤 𝑘+1 = (𝐵(𝑘+1)𝛿 − 𝐵 𝑘𝛿 )/ 𝛿 ∼ 𝒩 (0, 𝐼𝑛 ), √
b b𝐺,𝑘𝛿 (b 𝑦(𝑘+1)𝛿 = b 𝑦 𝑘𝛿 + 𝛿 𝑚 𝑦 𝑘𝛿 ) + 𝛿 𝑤 𝑘+1 . Lemma 34 (Lemma 4.10 of [EMS22]). Let 𝛽 < 𝛽0 , let 𝑐 = 𝑂(1), and fix a time horizon 𝑇SL < 𝑡 𝑛 and 𝜀 > 𝜀𝑛 . Let 𝐴 ∼ GOE(𝑛), and let 𝑦 ∈ ℝ𝑛 . Let 𝑞★(𝛽, 𝑡) be the AMP state-evolution overlap. Then for 𝜌 ∈ poly(𝜀) there exist 𝐾 AMP = 𝐾 AMP (𝛽, 𝑇SL , 𝜀) ≤ polylog(𝑛/𝜀) and 𝐾 NGD ≤ polylog(𝑛/𝜀) with the following property. With probability 1 − 𝑜 𝑛 (1) over (𝐴, 𝑦(·)), simultaneously for all 𝑡 ∈ (0, 𝑇SL ] the stationary point (fixed point of TAP) satisfies: √ ∥𝑚(𝐴, 𝑦(𝑡)) − 𝑚★(𝐴, 𝑦(𝑡); 𝑡)∥ ≤ 𝜌 𝑡𝑛, Proof. The proof is essentially the same as [EMS22]. For the choices of the parameters 𝐾 NGD and 𝐾 AMP , see Proposition 28. Define the normalized trajectory and mean-estimation errors 1 𝑦 Δ 𝑘 := √ b 𝑦 𝑘𝛿 − 𝑦 𝑘𝛿 , 𝑛
1 b𝐺,𝑘𝛿 (b Δ𝑚 𝑚 𝑦 𝑘𝛿 ) − 𝑚𝐺 (𝑦 𝑘𝛿 ) . 𝑘 := √ 𝑛
Lemma 35 (ASL error propagation, [EMS22]). Fix a time horizon 𝑡 < 𝑡 𝑛 (see remark Remark 25). In the high-temperature regime covered by the ASL analysis of [EMS22; EMS25], there exist constants 𝐶 = 𝑂 𝛽 (1), depending only on the temperature parameters, and a deterministic sequence 𝜉𝑛 → 0 such that the following holds with probability 1 − 𝑜 𝑛 (1) over the disorder. For every 𝑘 ≥ 0 satisfying 𝑘𝛿 ≤ 𝑡 and every 𝛿 ∈ (0, 1), √ √ 𝑦 Δ 𝑘 ≤ 𝐶𝑒 𝐶 𝑘𝛿 𝑘𝛿 𝜌 𝑘𝛿 + 𝛿 + 𝜉𝑛 , and
√ √ √ 𝐶 𝑘𝛿 Δ𝑚 ≤ 𝐶𝑒 𝑘𝛿 𝜌 𝑘𝛿 + 𝛿 + 𝐶𝜌 𝑘𝛿 + 𝜉𝑛 . 𝑘
Here 𝜌 = poly(𝜀) denotes the accuracy parameter in the tilted-mean estimation guarantee. Proof. We imitate the proof of Lemma 4.14 of [EMS22] with new parameters. We couple the continuous stochastic-localization trajectory and the discretized ASL trajectory using the same Brownian increments. Throughout the proof, let 𝜉𝑛 denote a deterministic non-negative sequence tending to zero. Set 𝑡 𝑘 = 𝑘𝛿. The proof is by induction on 𝑘. The case 𝑘 = 0 is immediate since b 𝑦0 = 𝑦0 = 0. Assume the claimed bounds hold up to time 𝑡 𝑘 . Subtracting the two coupled recursions gives
b b𝐺,𝑡 𝑘 (b 𝑦𝑡 𝑘+1 − 𝑦𝑡 𝑘+1 = b 𝑦𝑡 𝑘 − 𝑦𝑡 𝑘 + 𝛿 𝑚 𝑦𝑡 𝑘 ) − 𝑚𝐺 (𝑦𝑡 𝑘 ) + 36
∫ 𝑡 𝑘+1 𝑡𝑘
𝑚𝐺 (𝑦𝑡 𝑘 ) − 𝑚𝐺 (𝑦 𝑠 ) 𝑑𝑠.
Hence 1 𝑦 𝑦 Δ 𝑘+1 ≤ Δ 𝑘 + 𝛿Δ𝑚 𝑘 + √ 𝑛
∫ 𝑡 𝑘+1 𝑡𝑘
∥𝑚𝐺 (𝑦 𝑠 ) − 𝑚𝐺 (𝑦𝑡 𝑘 )∥ 𝑑𝑠.
By the continuity estimate for the tilted mean along the stochastic-localization trajectory [EMS22, Lemma 4.9], 1 √ 𝑛
∫ 𝑡 𝑘+1 𝑡𝑘
∥𝑚𝐺 (𝑦 𝑠 ) − 𝑚𝐺 (𝑦𝑡 𝑘 )∥ 𝑑𝑠 ≤ 𝐶𝛿 3/2 + 𝜉𝑛
with probability 1 − 𝑜 𝑛 (1). Therefore 𝑦
𝑦
3/2 Δ 𝑘+1 ≤ Δ 𝑘 + 𝛿Δ𝑚 + 𝜉𝑛 . 𝑘 + 𝐶𝛿
It remains to provide an analogous inequality for the mean-estimation bound. By the trajectory bound just proved, taking 𝜌, 𝛿 = poly(𝜀) as in Lemma 35, and following the analysis of [EMS22, Lemma 4.14] for mean estimation, we have √ 𝑦 Δ𝑚 𝑘+1 ≤ 𝐶Δ 𝑘+1 + 2𝜌 𝑡 𝑘+1 + 𝜉𝑛 . Here, 𝐶 is the Lipschitz constant of the AMP iteration, which is 𝑂(1) as long as 𝛽 < 𝛽0 is fixed. Now, we need 𝑦 to prove that Δ𝑚 and Δ𝐾 are 𝑂(𝜀) given 𝛿 and 𝜌, where 𝑡 𝐾 = 𝐾𝛿 = 𝑂(log(1/𝜀)). We upper bound Δ𝑚 and 𝐾 𝐾 𝑦 𝑚 Δ𝐾 by writing the following continuous differential-inequality upper bound for the variables 𝑀𝑡 = Δ𝑡 and 𝑦 𝑌𝑡 = Δ𝑡 : √ 𝑀𝑡 ≤ 𝐶𝑌𝑡 + 2𝜌 𝑡 √ 𝑑𝑌𝑡 ≤ 𝑀𝑡 + 𝐶1 𝛿
(134) (135)
Since 𝜌 and 𝛿 are polynomially small in 𝜀 and 𝐶1 is an absolute constant depending only on 𝛽 (see [EMS22, Lemma 4.9 and Lemma 4.5]), 𝑑𝑌𝑡 ≤ 𝐶𝑌𝑡 + 𝛼 where 𝛼 is polynomially small in 𝜀. Given 𝑌0 = 0, using the exact solution as an upper bound, we conclude that 𝑌𝐾𝛿 ≤ 𝐶𝛼 𝑒 𝐶𝐾𝛿 − 𝐶𝛼 = 𝑂(𝛼𝑒 𝐶𝐾𝛿 ) = poly(𝜀).
Theorem 36 (Discretization error at logarithmic localization time). Let the accuracy parameter satisfy 𝜀 > 𝜀𝑛 , and let 𝑡 𝜀 = 𝑂(log( 1𝜀 )) < 𝑡 𝑛 be the chosen localization time. Assume, for notational simplicity, that 𝐾𝛿 = 𝑡 𝜀 . Let 𝑦ˆ 𝐾𝛿 be the output of algorithmic stochastic localization using the parameter choices of Proposition 28 for TAP-AMP and the 𝜌 and 𝛿 of Lemma 35. Then, with probability 1 − 𝑜(1), 𝑊2,𝑛 ( 𝑦ˆ𝑡 𝜀 /𝑡 𝜀 , 𝑦𝑡 𝜀 /𝑡 𝜀 ) ≤ 𝑂(𝜀) Proof. Apply Lemma 35 at 𝑘 = 𝐾, so that 𝐾𝛿 = 𝑡 𝜀 . With 𝜌 = 𝐶𝑡 𝜀 Δ𝑚 𝑡𝜀 𝐾 ≤ 𝐶𝑒
√ √
√
𝛿, we have
√
√ √ 𝛿 𝑡 𝜀 + 𝛿 + 𝐶 𝛿 𝑡 𝜀 + 𝜉𝑛 .
Since 𝑡 𝜀 ≥ 1, the trajectory error is bounded by 1 𝑦 Δ𝐾 /𝑡 𝜀 ≤ poly(𝜀)/log( ) ≤ poly(𝜀) 𝜀 Using sign rounding (Lemma 10) and the stability argument Lemma 11, we obtain 𝑊2,𝑛 (sign( 𝑦ˆ𝑡 𝜀 ), 𝜇) ≤ 𝑂(𝜀) where 𝜇 denotes the true distribution. 37
A.2
Parallelization of TAP-AMP Mean Approximation Algorithms
e parallel-depth claim for the tilted-mean oracle used in Theorem 24. We keep the Here, we prove the 𝑂(1) notation of the main body: 𝐾 and 𝛿 denote, respectively, the number of stochastic-localization steps and the discretization step size, while 𝐾 AMP and 𝐾NGD denote the internal iteration counts of the TAP-AMP tilted mean estimator. The localization horizon is 𝑡 𝜀 = 𝐾𝛿, and by Theorem 24 we may take 𝑇𝜀 = 𝑂 log(1/𝜀) ,
𝛿 = poly(𝜀).
The overlap parameter in the mixed 𝑝-spin algorithm is denoted by 𝑞★(𝑡), as in Eq. (39). In the SK specialization we write 𝑞★SK (𝛽, 𝑡) for the corresponding fixed point in [EMS22]; this avoids overloading the mixed 𝑝-spin notation. Proof of Proposition 28. Both mean-estimation routines in the main body, Algorithm 6 for the SK model and Algorithm 5 for the mixed 𝑝-spin model, have the same structure: an AMP phase followed by a natural-gradient-descent phase on the approximate TAP free energy. We count parallel time in the PRAM model with parallel Hamiltonian-derivative oracle access, as in the oracle model described in the introduction. Equivalently, for the explicit fixed-degree mixed 𝑝-spin polynomial representation, every gradient or Hessianvector computation is a fixed-degree tensor contraction and has polylog(𝑛) parallel depth and polynomial work by parallel reductions. Consider one AMP iteration. The coordinate-wise nonlinearities tanh(·) and sech2 (·) are computed independently over the coordinates. The empirical overlap, Onsager coefficient, and inner products appearing in the update are computed by parallel reductions in 𝑂(log 𝑛) depth. The only model-dependent operation is the evaluation of ∇𝐻(𝑚 𝑘 ), or 𝐴𝑚 𝑘 in the SK specialization, which is one parallel Hamiltonian-derivative query. Hence each AMP iteration has polylog(𝑛) parallel depth and polynomial work. The same argument applies to one NGD iteration: it evaluates ∇ℱ̂TAP (𝑚; 𝑦, 𝑞) and then performs only coordinate-wise operations and parallel reductions. It remains to bound 𝐾 AMP and 𝐾 NGD . For the mixed 𝑝-spin model, the high-temperature TAP-AMP analysis of [EMS25] gives fixed constants, depending only on the mixture 𝜉 and on the requested tilted-mean accuracy, for which the AMP phase enters the locally strongly convex basin of ℱ̂TAP (·; 𝑦, 𝑞★(𝑡)) uniformly along the ASL trajectory. Taking the auxiliary tilted-mean accuracy to be 𝜌 = poly(𝜀) as in Lemma 35 gives 𝐾 AMP = polylog(𝑛/𝜀) for the parameter regime used by Theorem 24. In the SK case, Lemma 37 below gives the sharper explicit bound 𝐾 AMP = 𝑂 𝛽 log(1/𝛼) + log(1 + 𝑡 𝜀 ) , where 𝛼 is the local AMP accuracy parameter. With 𝛼 = poly(𝜀/𝑛) and 𝑡 𝜀 = 𝑂(log(1/𝜀)), this gives a polylog(𝑛/𝜀) bound. A similar upper bound can be proved for the mixed 𝑝-spin case, as long as the temperature is sufficiently high, by essentially imitating the proof for the SK case. For the NGD phase, the TAP analysis supplies a local strong-convexity and smoothness window around the AMP initialization. Namely, in the relevant neighborhood of the TAP fixed point, the objective in the natural parameter 𝑢 = arctanh(𝑚) is 𝜇-strongly convex and 𝑀-smooth, with 0 < 𝜇 ≤ 𝑀 < ∞ depending only on the high-temperature parameters. Taking 𝜂 ≤ 1/𝑀, the standard contraction estimate gives ∥𝑢𝑠 − 𝑢★∥2 ≤ (1 − 𝜂𝜇)𝑠 ∥𝑢0 − 𝑢★∥2 .
√ where 𝑢★ is the minimizer. Thus, reducing the NGD error to the auxiliary accuracy 𝜌 𝑡 𝜀 𝑛 requires 𝐾 NGD = 𝑂 𝛽 log(1/𝜌) + log 𝑛 + log(1 + 𝑇𝜀 ) = polylog(𝑛/𝜀),
again for 𝜌 = poly(𝜀/𝑛). Since every AMP and NGD iteration has polylog(𝑛) parallel depth and the number of such iterations is polylog(𝑛/𝜀), each tilted-mean call in Algorithm 4 is computable in parallel depth polylog(𝑛/𝜀), as claimed.
38
Lemma 37 (Explicit SK AMP iteration bound). Fix 0 < 𝛽 < 1/2, a horizon 𝑇 ≥ 1, and a local AMP accuracy parameter 𝛼 ∈ (0, 1). Let 𝑊 ∼ 𝒩 (0, 1) and define
2 √ mmse(𝛾) := 1 − 𝔼[tanh 𝛾 + 𝛾 𝑊 ],
𝛾 ≥ 0.
Let 𝛾0 (𝛽, 𝑡) = 0 and
𝛾𝑘+1 (𝛽, 𝑡) = 𝛽 2 1 − mmse 𝛾𝑘 (𝛽, 𝑡) + 𝑡 ,
𝑡 ∈ (0, 𝑇].
Let 𝛾★(𝛽, 𝑡) denote the nonnegative fixed point of this recursion, and define 𝑞 SK (𝛽, 𝑡) := 𝑘
𝛾𝑘+1 (𝛽, 𝑡) , 𝛽2
𝛾★(𝛽, 𝑡) . 𝛽2
𝑞★SK (𝛽, 𝑡) :=
Set 𝑎 𝛽 :=
𝛽2 , 1 − 𝛽2
1 − 2𝛽 , 4
𝜆𝛽 :=
and, for 𝑥 ∈ ℝ, write ⌈𝑥⌉+ := max{1, ⌈𝑥⌉}. Define 𝐾 𝑞 :=
log(16𝑎 𝛽 /𝛼)
−2 log 𝛽
, +
√ √ log 16 𝑎 𝛽 /(𝜆𝛽 𝛼) 𝐾 𝑔,1 := , − log 𝛽 + and
√ √ log 24𝑎 𝛽 𝑡 𝜀 /(𝜆𝛽 𝛼) 𝐾 𝑔,2 := . −2 log 𝛽 +
Let SK 𝐾 AMP := max{𝐾 𝑞 , 𝐾 𝑔,1 , 𝐾 𝑔,2 }.
For the SK AMP phase, let (𝑧 𝑘 , 𝑚 𝑘 ) 𝑘≥0 denote the iterates of Algorithm 6 with external field 𝑦(𝑡). Thus 𝑚−1 = 0, and
𝑧 0 = 0,
𝑚 𝑘 = tanh(𝑧 𝑘 ),
𝑛
𝑏 𝑘 :=
𝛽2 Õ sech2 (𝑧 𝑘 )𝑖 , 𝑛
𝑧 𝑘+1 = 𝛽𝐴𝑚 𝑘 + 𝑦(𝑡) − 𝑏 𝑘 𝑚 𝑘−1 ,
𝑖=1
SK where 𝐴 is the SK interaction matrix. Then 𝐾 AMP = 𝐾AMP is sufficient for the large-𝐾 AMP requirements used in [EMS22, Lemmas 4.10 and 4.11]:
|𝑞 SK (𝛽, 𝑡) − 𝑞 SK (𝛽, 𝑡)| 𝑘 𝑘 2
1
sup
sup
𝑡
𝑡∈(0,𝑇] 𝑘1 ,𝑘 2 ≥𝐾 AMP
and sup 𝑡∈(0,𝑇] 𝑞∈[𝑞 SK
≤
𝛼 , 16
√ 𝜆𝛽 𝛼 ∥∇ℱ̂TAP (𝑚𝐾AMP ; 𝑦(𝑡), 𝑞)∥2 sup ≤ . p-lim √ 4 𝑡𝑛 (𝛽,𝑡),𝑞★SK (𝛽,𝑡)] 𝑛→∞
𝐾 AMP
Consequently, SK 𝐾 AMP = 𝑂 𝛽 log(1/𝛼) + log(1 + 𝑇) .
39
Proof. Let 𝑓𝑡 (𝛾) := 𝛽2 1 − mmse(𝛾 + 𝑡) .
The recursion is 𝛾𝑘+1 (𝛽, 𝑡) = 𝑓𝑡 (𝛾𝑘 (𝛽, 𝑡)). Since mmse is decreasing on ℝ≥0 , 𝑓𝑡 is increasing, and hence 0 = 𝛾0 (𝛽, 𝑡) ≤ 𝛾1 (𝛽, 𝑡) ≤ 𝛾2 (𝛽, 𝑡) ≤ · · · ≤ 𝛾★(𝛽, 𝑡). By [EMS22, Lemma 4.5], for 𝛽 < 1 and 𝑡 > 0, 0 ≤ 𝛾★(𝛽, 𝑡) − 𝛾𝑘 (𝛽, 𝑡) ≤ 𝛽 2𝑘 𝛾★(𝛽, 𝑡), and
𝛾1 (𝛽, 𝑡) ≥ 1 − 𝛽2 . 𝛾★(𝛽, 𝑡)
Moreover, 𝛾1 (𝛽, 𝑡) = 𝛽 2 (1 − mmse(𝑡)) ≤ 𝛽 2 𝑡, using 1 − mmse(𝑡) ≤ 𝑡. Therefore 𝛾★(𝛽, 𝑡) ≤
𝛾1 (𝛽, 𝑡) ≤ 𝑎 𝛽 𝑡. 1 − 𝛽2
Thus 0 ≤ 𝑞★SK (𝛽, 𝑡) − 𝑞 SK (𝛽, 𝑡) = 𝑘
𝛾★(𝛽, 𝑡) − 𝛾𝑘+1 (𝛽, 𝑡) 𝛽2
≤ 𝛽2𝑘 𝛾★(𝛽, 𝑡) ≤ 𝑎 𝛽 𝛽 2𝑘 𝑡. (𝛽, 𝑡) is nondecreasing in 𝑘, for any 𝑘1 , 𝑘2 ≥ 𝐾, Since 𝑞 SK 𝑘 |𝑞 SK (𝛽, 𝑡) − 𝑞 SK (𝛽, 𝑡)| ≤ 𝑎 𝛽 𝛽 2𝐾 𝑡. 𝑘1 𝑘2 The first displayed requirement follows as soon as 𝑎 𝛽 𝛽 2𝐾 ≤ 𝛼/16, which is guaranteed by 𝐾 ≥ 𝐾 𝑞 . It remains to control the TAP-gradient term. Let Δ 𝑘 (𝑡) := 𝛾𝑘+1 (𝛽, 𝑡) − 𝛾𝑘 (𝛽, 𝑡). The monotonicity of 𝛾𝑘 and the preceding bound imply 0 ≤ Δ 𝑘 (𝑡) ≤ 𝛾★(𝛽, 𝑡) − 𝛾𝑘 (𝛽, 𝑡) ≤ 𝑎 𝛽 𝛽2𝑘 𝑡. State evolution for the SK AMP iteration gives ∥𝑧 𝑘+1 − 𝑧 𝑘 ∥22 p-lim 𝑛→∞
𝑛
= Δ 𝑘 (𝑡)2 + Δ 𝑘 (𝑡).
Consequently, uniformly over 𝑡 ∈ (0, 𝑇], ∥𝑧 𝑘+1 − 𝑧 𝑘 ∥2 p-lim ≤ √ 𝑡𝑛 𝑛→∞
r
Δ 𝑘 (𝑡)2 + Δ 𝑘 (𝑡) 𝑡 √ p ≤ 𝑎 𝛽 𝛽 𝑘 + 𝑎 𝛽 𝑇 𝛽2𝑘 .
Since tanh is 1-Lipschitz, ∥𝑚 𝑘 − 𝑚 𝑘−1 ∥2 ≤ ∥𝑧 𝑘 − 𝑧 𝑘−1 ∥2 . The TAP-gradient identity in [EMS22, Lemma 4.11] gives, for 𝑞 ∈ [𝑞 SK (𝛽, 𝑡), 𝑞★SK (𝛽, 𝑡)], 𝑘 ∥∇ℱ̂TAP (𝑚 𝑘 ; 𝑦(𝑡), 𝑞)∥2 ∥𝑧 𝑘+1 − 𝑧 𝑘 ∥2 ∥𝑚 𝑘−1 − 𝑚 𝑘 ∥2 ≤ + 𝛽2 √ √ √ 𝑡𝑛 𝑡𝑛 𝑡𝑛 SK SK 𝑞★ (𝛽, 𝑡) − 𝑞 𝑘 (𝛽, 𝑡) + 𝛽2 + 𝑜 𝑛,ℙ (1). √ 𝑡 40
Taking the 𝑝-limit and using the bounds above yields ∥∇ℱ̂TAP (𝑚 𝑘 ; 𝑦(𝑡), 𝑞)∥2 √ 𝑡𝑛 𝑡∈(0,𝑇] 𝑞∈[𝑞 SK (𝛽,𝑡),𝑞★SK (𝛽,𝑡)] 𝑛→∞ 𝑘 p p √ √ √ 𝑎 𝛽 𝛽 𝑘 + 𝑎 𝛽 𝑇 𝛽2𝑘 + 𝛽 2 𝑎 𝛽 𝛽 𝑘−1 + 𝑎 𝛽 𝑇 𝛽2𝑘−2 + 𝛽 2 𝑎 𝛽 𝑇 𝛽2𝑘 . ≤ sup
sup
p-lim
Because 0 < 𝛽 < 1, the right-hand side is at most √ p 2 𝑎 𝛽 𝛽 𝑘 + 3𝑎 𝛽 𝑇 𝛽2𝑘 . It is therefore enough to require 𝐾
2 𝑎𝛽 𝛽 ≤
p
√ 𝜆𝛽 𝛼
√ √ 𝜆𝛽 𝛼 2𝐾 3𝑎 𝛽 𝑇 𝛽 ≤ , 8
,
8
which are guaranteed by 𝐾 ≥ 𝐾 𝑔,1 and 𝐾 ≥ 𝐾 𝑔,2 , respectively. This proves the two large-𝐾 requirements. SK The asymptotic expression for 𝐾 AMP follows because 𝛽 ∈ (0, 1/2) is fixed, so − log 𝛽 is bounded away from zero.
A.3
Lipschitz Property of the TAP Fixed-Point Computation Algorithm for Mixed 𝑝-Spin Models
In this section, we prove that, for the mixed 𝑝-spin model at sufficiently high temperature, the TAP fixed-point computation algorithm is Lipschitz with respect to the tilt parameter 𝑦 with high probability. The proof follows the same two-phase structure as in the SK case: first an AMP phase gets close to the TAP fixed point, and then a natural-gradient descent step in the dual variables closes the remaining gap. We recall the standard notation for the mixed 𝑝-spin glass. First, the Hamiltonian takes the following form: 𝐻(𝑥) =
p 𝑃 Õ 𝛽 𝑝 𝑝! 𝑝−1 𝑝=2 𝑛 2
·
Õ
𝐺𝐽
𝐽⊆[𝑛],|𝐽|=𝑝
Ö
𝑥𝑘 ,
𝜇(𝑥) =
𝑘∈𝐽
1 exp(𝐻(𝑥)) 𝑍
(136)
where we assume zero external field. Also recall the associated mixture function 𝜉(𝑡) =
𝑃 Õ
𝛽 2𝑝 𝑡 𝑝
⟨𝑥, 𝑥 ′ ⟩ 𝔼[𝐻𝑛 (𝑥)𝐻𝑛 (𝑥 )] = 𝑛𝜉 . 𝑛
′
equivalently,
𝑝=2
Throughout this section, we condition on the high-probability event that the Hamiltonian has a dimensionfree Hessian bound. Namely, for a deterministic constant 𝐿𝐻 depending only on 𝜉, we condition on the following event: sup ∥∇2 𝐻𝑛 (𝑚)∥op ≤ 𝐿𝐻 .
(Hessian Bound)
𝑚∈[−1,1]𝑛
Using Gaussian concentration (see [AGJ20, Theorem 3.3]) one can prove that the above event holds with probability at least 1 − 𝑒 −𝑐𝑛 . Definition 38 (Lipschitz threshold for the mixed 𝑝-spin model). Let 4 𝐿0 := √ , 3 3
𝜙(𝑡) := (1 − 𝑡)𝜉′′ (𝑡),
41
𝐿 𝜙 := sup |𝜙 ′ (𝑡)|. 𝑡∈[0,1]
Fix 𝑞 ∈ (0, 1) and the quadratic TAP parameter 𝑇 = 𝑇(𝜉). Define 𝛾AMP (𝐿𝐻 , 𝜉) := 𝐿𝐻 + 𝜉′′ (1) + 𝐿0 𝐿 𝜙 , 3𝑇 + (1 − 𝑞)𝜉′′ (𝑞), 𝛾NGD (𝐿𝐻 , 𝜉, 𝑞, 𝑇) := 𝐿𝐻 + 2
(AMP-Lipchitz) (NGD-Lipchitz)
and 𝛾mix (𝐿𝐻 , 𝜉, 𝑞, 𝑇) := max 𝛾AMP (𝐿𝐻 , 𝜉), 𝛾NGD (𝐿𝐻 , 𝜉, 𝑞, 𝑇) .
The TAP fixed-point computation for the 𝑝-spin model is very similar to that for the SK model, but it is more complicated because it uses derivatives of the mixture function. Here, we review the TAP fixed-point computation used in [EMS23]. Phase 1: AMP phase.
We have the following initialization: 𝑚−1 = 0,
𝑧 0 = 0.
For 𝑘 = 0, 1, . . . , 𝐾AMP − 1, define 𝑚 𝑘 = tanh(𝑧 𝑘 ), 𝑛
𝑞𝑘 =
1 1Õ tanh2 (𝑧 𝑘 )𝑖 = ∥𝑚 𝑘 ∥2 , 𝑛 𝑛 𝑖=1
𝑏 𝑘 = (1 − 𝑞 𝑘 )𝜉′′ (𝑞 𝑘 ), and 𝑧 𝑘+1 = ∇𝐻𝑛 (𝑚 𝑘 ) + 𝑦 − 𝑏 𝑘 𝑚 𝑘−1 . Phase 2: NGD phase on the approximate TAP free energy.
(138)
Initialize
𝑢0 := 𝑧 𝐾AMP . For 𝑡 = 0, 1, . . . , 𝐾NGD − 1, define 𝑚𝑡+ = tanh(𝑢𝑡 ), and 𝑢𝑡+1 = 𝑢𝑡 − 𝜂∇𝑚 ℱ̂TAP (𝑚𝑡+ ; 𝑦, 𝑞).
(139)
The algorithm outputs
b (𝐺, 𝑦) := 𝑚𝐾+NGD = tanh(𝑢𝐾NGD ). 𝑚 Definition 39 (Tilted-mean computation for mixed 𝑝-spin). For a given disorder realization 𝐺, external field 𝑦 ∈ ℝ𝑛 , and iteration counts 𝐾AMP and 𝐾 NGD , let
b(𝐾AMP ,𝐾NGD ) (𝐺, 𝑦) 𝑚 denote the output of the AMP+NGD algorithm above. When 𝐺, 𝐾 AMP , and 𝐾 NGD are clear from context, we b (𝑦). simply write 𝑚 Our goal is to prove that, for any two tilts 𝑦, e 𝑦 ∈ ℝ𝑛 , the output satisfies
b (𝐺, 𝑦) − 𝑚 b (𝐺, e ∥𝑚 𝑦 )∥ ≤
1 ∥𝑦 − e 𝑦 ∥, 1 − 𝛾mix
provided 𝛾mix < 1. The proof uses the following elementary Lipschitz facts: • tanh is 1-Lipschitz; √ • tanh2 is 𝐿0 = 4/(3 3)-Lipschitz; 42
• 𝜙(𝑡) = (1 − 𝑡)𝜉′′ (𝑡) is 𝐿 𝜙 -Lipschitz on [0, 1]. Throughout the proof, fix two external fields 𝑦 and e 𝑦 . Unadorned variables correspond to the computation with tilt 𝑦, and tilded variables correspond to the computation with tilt e 𝑦 . We use the notation Δ𝑦 := 𝑦 − e 𝑦,
Δ𝑧 𝑘 := 𝑧 𝑘 − e 𝑧𝑘 ,
e𝑘 , Δ𝑚 𝑘 := 𝑚 𝑘 − 𝑚
Δ𝑞 𝑘 := 𝑞 𝑘 − e 𝑞𝑘 ,
Δ𝑏 𝑘 := 𝑏 𝑘 − e 𝑏𝑘 ,
Δ𝑢𝑡 := 𝑢𝑡 − e 𝑢𝑡 .
e−1 = 0 by initialization. Also, Δ𝑧0 = 0, Δ𝑧−1 = 0, and 𝑚−1 = 𝑚 Theorem 40 (Lipschitzness of AMP phase for mixed 𝑝-spin). Assume the Hessian bound (137) holds and 𝛾AMP (𝐿𝐻 , 𝜉) = 𝐿𝐻 + 𝜉′′ (1) + 𝐿0 𝐿 𝜙 < 1. Then for all 𝑘 ≥ 0, ∥Δ𝑧 𝑘 ∥ ≤
1 ∥Δ𝑦∥, 1 − 𝛾AMP
∥Δ𝑚 𝑘 ∥ ≤
1 ∥Δ𝑦∥. 1 − 𝛾AMP
Proof. Let 𝑉𝑘 := max ∥Δ𝑧 𝑗 ∥. 0≤𝑗≤𝑘
We first prove the one-step inequality ∥Δ𝑧 𝑘+1 ∥ ≤ 𝐿𝐻 + 𝐿0 𝐿 𝜙 ∥Δ𝑧 𝑘 ∥ + 𝜉′′ (1)∥Δ𝑧 𝑘−1 ∥ + ∥Δ𝑦∥.
Taking the difference of (138) for 𝑦 and e 𝑦 gives
e 𝑘 )] + Δ𝑦 − 𝑏 𝑘 Δ𝑚 𝑘−1 − Δ𝑏 𝑘 𝑚 e 𝑘−1 . Δ𝑧 𝑘+1 = [∇𝐻𝑛 (𝑚 𝑘 ) − ∇𝐻𝑛 (𝑚 Since (137) holds, ∇𝐻𝑛 is 𝐿𝐻 -Lipschitz on [−1, 1]𝑛 . Since tanh is 1-Lipschitz, ∥Δ𝑚 𝑘 ∥ ≤ ∥Δ𝑧 𝑘 ∥,
∥Δ𝑚 𝑘−1 ∥ ≤ ∥Δ𝑧 𝑘−1 ∥.
Moreover, 0 ≤ 𝑏 𝑘 = (1 − 𝑞 𝑘 )𝜉′′ (𝑞 𝑘 ) ≤ 𝜉′′ (1). It remains to control the Onsager coefficient difference. Since 𝑏 𝑘 = 𝜙(𝑞 𝑘 ), |Δ𝑏 𝑘 | ≤ 𝐿 𝜙 |Δ𝑞 𝑘 |. Using the 𝐿0 -Lipschitzness of tanh2 , 𝑛
|Δ𝑞 𝑘 | =
1 Õ tanh2 (𝑧 𝑘 )𝑖 − tanh2 (e 𝑧 𝑘 )𝑖 𝑛 𝑖=1
𝐿0 𝐿0 ≤ ∥Δ𝑧 𝑘 ∥1 ≤ √ ∥Δ𝑧 𝑘 ∥. 𝑛 𝑛
e 𝑘−1 ∥ ≤ Since ∥𝑚
√
𝑛, we get
e 𝑘−1 ∥ ≤ 𝐿0 𝐿 𝜙 ∥Δ𝑧 𝑘 ∥. |Δ𝑏 𝑘 | ∥𝑚
Combining the above estimates with the triangle inequality proves (140). Therefore, ∥Δ𝑧 𝑘+1 ∥ ≤ 𝛾AMP 𝑉𝑘 + ∥Δ𝑦∥, and hence 𝑉𝑘+1 ≤ max{𝑉𝑘 , 𝛾AMP𝑉𝑘 + ∥Δ𝑦∥}. Since 𝛾AMP < 1, induction gives 𝑉𝑘 ≤
1 ∥Δ𝑦∥ 1 − 𝛾AMP
for every 𝑘. The same bound for Δ𝑚 𝑘 follows from the 1-Lipschitzness of tanh. 43
(140)
Lipschitzness of the full AMP+NGD algorithm for mixed 𝑝-spin We next prove that the NGD phase preserves the AMP Lipschitz bound, up to replacing 𝛾AMP by the larger stability parameter 𝛾mix . Rewrite the gradient of TAP in 𝑢-space. We recall several definitions from the preliminaries regarding the TAP equation.
ℱTAP (𝑚, 𝑦) = −𝐻(𝑚) − 𝑦, 𝑚 −
𝑛 Õ
ℎ(𝑚 𝑖 ) − ONS(𝑄(𝑚))
(141)
𝑖=1
𝑛 [𝜉(1) − 𝜉(𝑄) − (1 − 𝑄)𝜉′ (𝑄)] 2 𝑛 ONS′ (𝑞) = − (1 − 𝑞)𝜉′′ (𝑞) 2 1+𝑚 1+𝑚 1−𝑚 1−𝑚 −1 2 𝑄(𝑚) = 𝑛 ∥𝑚∥ , ℎ(𝑚) = − ln − ln 2 2 2 2 ONS(𝑞) =
(142) (143) (144)
A direct calculation gives ∇𝑚 𝐹TAP (𝑚; 𝑦, 𝑞) = −∇𝐻𝑛 (𝑚) − 𝑦 + arctanh(𝑚) + (1 − 𝑞)𝜉′′ (𝑞)𝑚 𝑇 + 𝑄(𝑚) − 𝑞 𝑚. 2 Plugging in 𝑚 = tanh(𝑢) gives ∇𝑚 𝐹TAP (tanh(𝑢); 𝑦, 𝑞) = 𝑢 − 𝑦 − 𝒢(𝑢), where 𝒢(𝑢) := ∇𝐻𝑛 (tanh(𝑢)) − (1 − 𝑞)𝜉′′ (𝑞) tanh(𝑢) −
𝑇 𝐽(tanh(𝑢)), 2
and 𝐽(𝑚) := 𝑄(𝑚) − 𝑞 𝑚.
Thus the NGD update (139) becomes 𝑢𝑡+1 = (1 − 𝜂)𝑢𝑡 + 𝜂𝑦 + 𝜂𝒢(𝑢𝑡 ).
(145)
Lemma 41 (Lipschitzness of the quadratic TAP correction). For 𝐽(𝑚) = 𝑄(𝑚) − 𝑞 𝑚 with 𝑞 ∈ (0, 1),
sup ∥∇𝐽(𝑚)∥op ≤ 3.
𝑚∈[−1,1]𝑛
Proof. Since ∇𝐽(𝑚) = 𝑄(𝑚) − 𝑞 𝐼 +
and 𝑄(𝑚) ∈ [0, 1], 𝑞 ∈ (0, 1), and ∥𝑚∥ ≤
√
2𝑚𝑚 ⊤ , 𝑛
𝑛, we have
∥∇𝐽(𝑚)∥op ≤ 1 +
44
2∥𝑚∥2 ≤ 3. 𝑛
Theorem 42 (TAP Lipschitz property for mixed 𝑝-spin). Assume the Hessian bound (137) holds and 𝛾mix (𝐿𝐻 , 𝜉, 𝑞, 𝑇) < 1. Let 0 < 𝜂 ≤ 1. Set 𝑢0 = 𝑧 𝐾AMP and e 𝑢0 = e 𝑧 𝐾AMP . Then for every 𝑡 ≥ 0, ∥𝑢𝑡 (𝑦) − 𝑢𝑡 (e 𝑦 )∥ ≤
1 ∥𝑦 − e 𝑦 ∥. 1 − 𝛾mix
Consequently, the final output satisfies
b(𝐾AMP ,𝐾NGD ) (𝐺, 𝑦) − 𝑚 b(𝐾AMP ,𝐾NGD ) (𝐺, e ∥𝑚 𝑦 )∥ ≤
1 ∥𝑦 − e 𝑦 ∥. 1 − 𝛾mix
If the NGD phase is run to convergence and the limit is the TAP fixed point 𝑚TAP (𝐺, 𝑦), then the same Lipschitz bound holds for the TAP fixed point map 𝑦 ↦→ 𝑚TAP (𝐺, 𝑦). Proof. By Theorem 40 and the initialization of NGD, ∥Δ𝑢0 ∥ = ∥Δ𝑧 𝐾AMP ∥ ≤
1 1 ∥Δ𝑦∥ ≤ ∥Δ𝑦∥. 1 − 𝛾AMP 1 − 𝛾mix
We now show that the NGD step preserves this bound. First, 𝒢 is 𝛾NGD -Lipschitz. Indeed, by the Hessian bound, the 1-Lipschitzness of tanh, and Lemma 41, ∥𝒢(𝑢) − 𝒢(𝑣)∥ ≤ 𝐿𝐻 ∥𝑢 − 𝑣∥ + (1 − 𝑞)𝜉′′ (𝑞)∥𝑢 − 𝑣∥ +
3𝑇 ∥𝑢 − 𝑣∥ 2
= 𝛾NGD ∥𝑢 − 𝑣∥. Subtracting (145) for 𝑦 and e 𝑦 gives ∥Δ𝑢𝑡+1 ∥ ≤ (1 − 𝜂)∥Δ𝑢𝑡 ∥ + 𝜂∥𝒢(𝑢𝑡 ) − 𝒢(e 𝑢𝑡 )∥ + 𝜂∥Δ𝑦∥ ≤ (1 − 𝜂 + 𝜂𝛾NGD )∥Δ𝑢𝑡 ∥ + 𝜂∥Δ𝑦∥.
Assume inductively that ∥Δ𝑢𝑡 ∥ ≤
1 ∥Δ𝑦∥. 1 − 𝛾mix
Since 𝛾NGD ≤ 𝛾mix , we obtain ∥Δ𝑢𝑡+1 ∥ ≤ 1 − 𝜂 + 𝜂𝛾NGD ≤
1 ∥Δ𝑦∥ + 𝜂∥Δ𝑦∥ 1 − 𝛾mix
1 ∥Δ𝑦∥. 1 − 𝛾mix
Thus the bound holds for all 𝑡. Finally, since tanh is 1-Lipschitz,
b (𝐺, 𝑦) − 𝑚 b (𝐺, e ∥𝑚 𝑦 )∥ ≤ ∥𝑢𝐾NGD (𝑦) − 𝑢𝐾NGD (e 𝑦 )∥ ≤
1 ∥𝑦 − e 𝑦 ∥. 1 − 𝛾mix
If 𝑢𝐾NGD (𝑦) and 𝑢𝐾NGD (e 𝑦 ) converge as 𝐾 NGD → ∞, the same inequality passes to the limit, giving the TAP fixed-point statement. Corollary 43 (High-probability Lipschitzness for fixed mixed 𝑝-spin). Fix a finite mixture 𝜉, 𝑞 ∈ (0, 1), and 𝜀 > 0. Assume that the temperature coefficients {𝛽 𝑝 }𝑃𝑝=2 satisfy ℭ({𝛽 𝑝 }) B
𝑃 Õ
q
𝛽 𝑝 𝑝 3 log 𝑝 < 𝛾0 ,
𝔇({𝛽 𝑝 }) B
𝑝=2
𝑃 Õ 𝑝=2
45
q
𝛽 𝑝 2𝑝 𝑝 3 log 𝑝 < ∞.
Then, with high probability over the disorder realization 𝐺, the TAP-AMP mean computation in Algorithm 5 is Lipschitz. Specifically, for all 𝑦, e 𝑦 ∈ ℝ𝑛 ,
b (𝐺, 𝑦) − 𝑚 b (𝐺, e ∥𝑚 𝑦 )∥ ≤ 𝑂 𝜉,𝑞,𝑇 (1) ∥𝑦 − e 𝑦 ∥. Equivalently, the contraction parameter satisfies 𝛾mix (𝐿𝐻 , 𝜉, 𝑞, 𝑇) B max 𝛾AMP (𝐿𝐻 , 𝜉), 𝛾NGD (𝐿𝐻 , 𝜉, 𝑞, 𝑇) < 1.
Proof. Under the high-temperature assumption, the constants appearing in the TAP-AMP Lipschitz analysis, namely 𝐿 𝜙 , 𝐿𝐻 , 𝑇, and the relevant derivatives of 𝜉, are bounded by a constant 𝑐 = 𝑐(𝜉, 𝑞, 𝑇) which can be made sufficiently small by taking the temperature coefficients small enough. In particular, using the bounds from the TAP-AMP analysis, 𝛾AMP (𝐿𝐻 , 𝜉) = 𝐿𝐻 + 𝜉′′ (1) + 𝐿0 𝐿 𝜙
(146)
4 ≤ 𝑐 2+ √ . 3 3
(147)
Similarly, 𝛾NGD (𝐿𝐻 , 𝜉, 𝑞, 𝑇) = 𝐿𝐻 + ≤
3𝑇 + (1 − 𝑞)𝜉′′ (𝑞) 2
7 𝑐. 2
(148) (149)
Therefore, choosing the high-temperature constant 𝛾0 sufficiently small ensures that both 𝛾AMP (𝐿𝐻 , 𝜉) < 1
and
𝛾NGD (𝐿𝐻 , 𝜉, 𝑞, 𝑇) < 1.
Hence 𝛾mix (𝐿𝐻 , 𝜉, 𝑞, 𝑇) < 1. By the Lipschitzness criterion for the TAP-AMP mean computation, this implies that, with high probability over 𝐺, b (𝐺, 𝑦) − 𝑚 b (𝐺, e ∥𝑚 𝑦 )∥ ≤ 𝑂 𝜉,𝑞,𝑇 (1) ∥𝑦 − e 𝑦∥ for all 𝑦, e 𝑦 ∈ ℝ𝑛 .
A.4
exp(poly(1/𝜀)) Steps for Algorithmic Stochastic Localization
Proof of Proposition 23. Let 𝑡 𝜀 = 𝐾𝛿 be the localization time horizon, using notation consistent with El Alaoui, Montanari, and Sellke [EMS22]. Their Lemma √ 4.14 gives, for the Euler ASL trajectory and with the mean-estimation accuracy parameter chosen as 𝜌 = 𝛿, √ √ √ 𝐵𝐾 ≤ 𝐶𝑒 𝐶𝑡 𝜀 𝑡 𝜀 𝜌 𝑡 𝜀 + 𝛿 + 𝐶𝜌 𝑡 𝜀 + 𝑜 𝑛 (1). Here, the parameter 𝐶 is of order 6𝐾AMP = poly(1/𝜀) (Proposition 28). Hence, for 𝑡 𝜀 ≥ 1, √ 3/2 𝐵𝐾 ≤ 𝐶𝑒 𝐶𝑡 𝜀 𝑡 𝜀 𝛿 + 𝑜 𝑛 (1). −1/2
Combining this with the localization error 𝑡 𝜀 , the EMS proof obtains the pre-rounding bound ([EMS22, Lemma 4.15]) √ −1/2 3/2 b (𝐴, b 𝔼𝑊2,𝑛 𝜇𝐴 , ℒ(𝑚 𝑦𝐾 )) ≤ 𝑡 𝜀 + 𝐶𝑒 𝐶𝑡 𝜀 𝑡 𝜀 𝛿 + 𝑜 𝑛 (1). To make the final rounded output have 𝑊2,𝑛 -error at most 𝜀, their rounding lemma requires the pre-rounding error to be 𝑂(𝜀2 ). Thus, they take their parameters to be √ −3/2 𝑡 𝜀 = poly(1/𝜀), 𝛿 ≤ 𝐶𝜀2 𝑒 −𝐶𝑡 𝜀 𝑡 𝜀 . 46
Equivalently, 𝛿 ≤ 𝐶𝜀4 𝑒 −2𝐶𝑡 𝜀 𝑡 𝜀−3 . Therefore the number of discretization steps satisfies 𝐾=
𝑡𝜀 ≥ 𝐶 −1 𝜀−4 𝑒 2𝐶𝑡 𝜀 𝑡 𝜀4 = exp poly(1/𝜀) . 𝛿
47