Conceptio › Archive › arXiv CS

arXiv CSopen access

High-Magnetization Sampling at Low Temperatures: Ising Models and Bayesian Sparse Linear Regression

· arxiv_cs

arXiv CS · Papers · License: Open Access

Open Source ↗Direct PDF ↓

neural-networks

machine learning, deep learning, neural networks

High-Magnetization Sampling at Low Temperatures: Ising Models and Bayesian Sparse Linear Regression Syamantak Kumar∗

Purnamrita Sarkar†

Kevin Tian‡

Yusong Zhu§

arXiv:2609.08873v1 [cs.DS] 8 Sep 2026

Abstract Sparsity is a powerful structural resource in optimization and statistics. We develop frameworks for leveraging sparsity in sampling problems over the Hamming slice Xkd := {x ∈ {±1}d : |{i : xi = 1}| = k}, in high-dimensional regimes where k ≪ d (i.e., where Xkd is highlymagnetized ). We use our frameworks to design improved samplers for canonical problems in the study of Ising models and Bayesian sparse linear regression. Our first main result considers the Sherrington-Kirkpatrick (SK) model, restricted to fixedmagnetization slices Xkd . We give a polynomial-time sampler for fixed-magnetization SK models at any inverse temperature β > 0, under arbitrary external fields, provided k ≤ cβ d for an appropriate constant cβ . By combining this result with an annealing strategy for estimating normalizing constants, this yields polynomial-time samplers for the SK model at arbitrarily low temperatures, under a sufficiently strong external field strength h. In the large β limit, our framework permits sampling at field strengths within constant factors of the AlmeidaThouless line delineating the replica symmetric and replica symmetry breaking regions [dAT78], improving polynomially over the h(β) required by recent work of [BAR26]. Our second main result concerns the measurement complexity of polynomial-time Bayesian sparse linear regression. Recent work by [KSTZ25] shows how to sample from the canonical Gaussian spike-and-slab posterior model with expected sparsity k, at any signal-to-noise ratio, given n ≳ k 3 log3 d Gaussian measurements. We improve this to n ≳ k 1.5 log2 d + k log3 d, using a common sparsity-aware framework underlying our results on Ising models.

∗

University of Texas at Austin, [email protected] University of Texas at Austin, [email protected] ‡ University of Texas at Austin, [email protected] § University of Texas at Austin, [email protected] †

Contents 1 Introduction 1.1 Our results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Our techniques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Related work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1 3 4 7

2 Preliminaries 2.1 Notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Markov chains . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Statistical models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8 8 10 12

3 Sparse Dobrushin Condition 13 3.1 Basic analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.2 Fixed-magnetization Ising models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 4 Spectral Mixing via Trickle Down 18 4.1 Trickle down framework . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 4.2 Simple sufficient conditions for fast mixing . . . . . . . . . . . . . . . . . . . . . . . . 22 4.3 Linear magnetization in low-temperature SK models . . . . . . . . . . . . . . . . . . 25 5 Annealing 28 5.1 Estimating normalizing constants . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 5.2 Approximate bounded-magnetization sampling . . . . . . . . . . . . . . . . . . . . . 30 5.3 Bounded-magnetization Ising models . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 6 Bayesian Sparse Linear Regression 33 6.1 Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 6.2 Trickle down for support posterior . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 6.3 Main result . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 AI Disclosure

A Sampling Near the Almeida–Thouless Line 51 A.1 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 A.2 Sampling around 1d . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 A.3 Sampling around the mean . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 B Scaling of Infinite ∆-Regular Tree Threshold

Introduction

Harnessing sparsity is a central theme in modern high-dimensional optimization and statistics. From a sample complexity perspective, sparsity is often a blessing: classical results on Gelfand widths [Kas77, GG84] imply that, in the well-studied sparse linear regression (SLR) problem, only n ≈ k log d noisy measurements y = Xθ ⋆ + ξ are information-theoretically sufficient to estimate a k-sparse signal θ ⋆ ∈ Rd up to the noise level ∥ξ∥2 . From an algorithm design perspective, however, imposing sparsity can create highly complex, nonconvex landscapes. Seminal work in compressed sensing overcame this tension by identifying structural conditions, such as the restricted isometry property, under which sparse recovery is possible in polynomial time [CT05, CT06, CRT06, Don06], leading to a broad theory of sparsity-aware optimization [Wai19]. Our goal is to develop an analogous framework for exploiting sparse structure in high-dimensional sampling. In the problems we study, sparsity restricts the number of simultaneously “active” coordinates without fixing their locations. A central motivating problem in this work is sampling from the Sherrington-Kirkpatrick (SK) model (cf. Model 1) under sparsity constraints. Interestingly, understanding sampling algorithms for the SK model under sparsity has other consequences, including improved algorithms for posterior sampling in Bayesian sparse linear regression. SK model. The SK model is a well-studied special case of the Ising model, a measure over vectors of spins from the hypercube x ∈ X d := {±1}d . Such x can be identified with a set S ⊆ [d], the locations of positive spins xi = 1. In Ising models, the underlying measure is 1 ⊤ ⊤ µ(x) ∝ exp β x Jx + h x , x ∈ X d, (1) 2 governed by an interaction matrix J ∈ Rd×d , an external field h ∈ Rd , and an inverse temperature β > 0. Ising models are a canonical testbed across statistical physics, machine learning, and theoretical computer science [WJ08, LPW09, Tal10]. For suitable J, the Gibbs landscape in (1) is known to undergo qualitative phase transitions as β varies [Tal10]. Efficient algorithms exist for sampling under the SK model, where J is drawn from the Gaussian orthogonal ensemble (GOE), for β up to a universal constant c [EKZ22, AJK+ 22, AMS22, AKV24, DLSS26], while conditional hardness has been demonstrated at a higher constant threshold c′ > c [AMS22]. Sparsity from field strength. As intuition for our central high-magnetization SK model, to be introduced next, we describe a well-studied analog: the SK model under a strong positive external field h = θ1d . When θ is large, (1) is biased toward 1d , driving the minority-spin set to become sparse. The high-field strength model thus exhibits a soft form of sparsity. On the algorithmic side, [BAR26] proved polynomial-time mixing of the Glauber dynamics for the SK model at every inverse temperature β > 0, under a high enough field strength h := βθ. They derive their sampler as a consequence of a more general result that proves mixing in Ising models (1), under a condition on the sparse operator norm of J at a constant sparsity scale k = Θ(n). On the geometric side, the Almeida–Thouless (AT) line is the standard benchmark for replica symmetry in the (β, h)-plane. This line, derived by [dAT78] using the replica method (and described in Appendix A), predicts a qualitative shift in the behavior of (1) as h grows compared to β. This line provides a natural field strength scale against which to compare [BAR26] to our results.

The role of sparsity in [BAR26] is implicit: a large field strength h causes independent replicas to have high overlap, so their disagreement set is typically sparse. The relevant interaction depends on a principal submatrix of J, giving a dependence on sparse operator norms. This suggests a complementary question: can high magnetization itself, imposed as a hard constraint rather than induced by an external field, make sampling in the SK model tractable? Sparsity from high magnetization. High-magnetization regimes have long played an important role in statistical physics, e.g., the study of spontaneous magnetization and large deviations [Yan52, Bon14, ADCS14, Ell12]. For measures on X d , high magnetization is a natural analog of sparsity. Under our convention, the positive spins are the active coordinates: an element x of o n (2) Xkd := x ∈ X d : {i : xi = 1} = k has exactly k positive spins and, for k ≤ d2 , magnetization magnitude d − 2k. Taking k ≪ d imposes a sparsity constraint, while maintaining an outcome space of exponential size kd ≈ exp(k log kd ). We study the SK model restricted to both the fixed-magnetization slice Xkd , and its boundedSk d d := magnetization counterpart X≤k i=0 Xi . We now state our first central problem. For every fixed β < ∞, is there a polynomial-time sampler from the fixed-magnetization SK model whenever k ≤ cβ d for a fixed cβ > 0? Theorem 1 answers this question affirmatively, taking the algorithm as the down-up walk. Bayesian SLR. Our second motivating example is a Bayesian variant of SLR. SLR is often phrased as an optimization problem: given noisy measurements (X, y = Xθ ⋆ + ξ), return a k-sparse θb ∈ Rd b (approximately) minimizing the residual norm ∥Xθ−y∥ 2 . When X is RIP, low residual error implies accurate estimation [CRT06], so optimization produces a good point estimate of θ ⋆ . In many applications, however, it is preferable to sample θb from a distribution over plausible signals, e.g., to model uncertainty in support estimation (variable selection) [MB88, GM93]. In the wellestablished Bayesian SLR model, the noise ξ is Gaussian with coordinatewise variance σ 2 , where σ −1 > 0 is a signal-to-noise ratio. For an appropriate prior π over sparse signals θ ⋆ , the goal is then to sample from the posterior induced by the observations: θb ∼ π (· | X, y) , where y = Xθ ⋆ + ξ, θ ⋆ ∼ π, ξ ∼ N 0n , σ 2 In . (3) Samples from the posterior density can then be used in downstream tasks, e.g., constructing credible intervals. In the statistics literature, π is often taken to be the spike-and-slab prior, O k k π= 1− δ0 + N (0, 1), (4) d d i∈[d]

where each signal coordinate θi⋆ is independently set to 0 except with low probability (so the expected sparsity is k). The induced spike-and-slab posterior sampling problem is often called the “theoretical gold standard” for modeling uncertainty in variable selection [JS04, CPS09, IR11, CvdV12, Roc18, PS19]. Unfortunately, this task poses a notorious computational challenge [CSHVdV15], and many heuristics have been developed as approximations [BRG21]. 2

Recently, several works developed provable methods for this sampling problem [YWJ16, MW26], including an algorithm by [KSTZ25] which samples from the posterior density (3), (4) given a sublinear-in-d, n ≳ k 3 log3 d measurements. The counterpart optimization problem is known to be feasible even when n ≈ k log d, motivating our second central problem. What is the measurement threshold n at which spike-and-slab posterior sampling admits polynomial-time algorithms, for any signal-to-noise ratio? We show that the complex SLR posterior density admits enough structure to be captured by our analysis framework for the high-magnetization SK model, and give a polynomial-time posterior sampling algorithm at n ≳ k 1.5 polylog(d) measurements (Theorem 2).

1.1

Our results

In this section, we overview our main results. To obtain these results, we develop a suite of technical tools for exploiting sparsity in sampling, described at more depth in Section 1.2. High-magnetization SK models. Our first main result concerns sampling in high-magnetization SK models, where J is GOE and the external field h is arbitrary. As a benchmark, [AMS22] showed conditional hardness for sampling in SK models at large inverse temperatures β > c′ for a constant c′ . We show that after fixing the magnetization k ≤ cβ d of the spin vector x ∈ Xkd , the SK model admits polynomial-time sampling at any temperature. Theorem 1 (Informal; see Theorem 4). For any fixed β > 0, h ∈ Rd , and J ∼ GOE(d), there is a constant cβ > 0 such that if k ≤ cβ d, we can sample from the SK model (1) restricted to the Hamming slice Xkd in polynomial time, with high probability over J. The algorithm in Theorem 1 is the canonical down-up walk over the Hamming slice Xkd , and we e 2 + dk 2 max(1, β ∥h∥ )).1 Beyond Theorem 1, Section 4 bound its runtime by a relatively mild O(d ∞ develops a more general framework for proving Poincaré inequalities on the down-up walk for Ising models, using the trickle down (“local-to-global”) theorem of [Opp18, AL20]. For example, we state an analogous result for Gaussian Hopfield models in Corollary 4. Due to Theorem 1, we recover a variant of the main result of [BAR26] by exploiting the relationship between high field strength and high magnetization.2 We do note that [BAR26] directly analyze the Glauber dynamics, perhaps the simplest sampling algorithm over X d , whereas we use an annealing scheme (Section 5) to reduce bounded-magnetization sampling to fixed-magnetization sampling. Our framework yields other interesting consequences beyond the linear sparsity regime in Theorem 1. For example, for k = O(1) independent of d → ∞, our results imply polynomial-time p sampling in high-magnetization SK models for polynomial inverse temperatures β, up to O( d/ log(d)) (cf. Corollary 3). Prior results exploiting notions of sparsity (albeit different from ours) to sample from Ising models at subconstant temperatures only tolerated β ≈ log d [CDKP22, KPPY25]. This comparison is discussed at more length in Section 1.3 and Appendix B. The Almeida-Thouless line. From a quantitative perspective, the AT line is a useful benchmark for comparing our sampling result with [BAR26]. The AT line predicts a phase transition in the 1 2

e to suppress logarithmic factors. In this introduction only, we use O Formally, this follows by combining Corollary 5 and Lemma 24.

geometry of the SK model under an external field h = βh 1d , once the field strength h ≥ hAT (β) is √ large enough as a function of β, where hAT (β) = (1 + o(1))β 2 log β (Lemma 23). Leveraging a simple reduction from high field strength sampling to high-magnetization sampling (Lemma 24), we show in Theorem 6 that √Theorem 1 implies sampling from the SK model whenever h ≥ hHM (β) for a threshold hHM (β) = 2(1 + o(1))hAT (β), i.e., within constant factors of the AT line. We also show in Theorem 7 that, if given access to the signs of the mean vector E[x] under the SK model, a modification of our sampler succeeds at any field √ strength h ≥ (1 + o(1))hAT (β). By comparison, the result of [BAR26] applies whenever h = Ω(β 2 log β), a polynomial factor larger (see discussion in Appendix A). The tighter range of h tolerated by our framework is a result of our basic sampler in Theorem 1 applying for an arbitrary external field h. Bayesian SLR. Finally, we apply our sparse sampling frameworks to spike-and-slab posterior sampling. We obtain a state-of-the-art measurement complexity for a canonical variant of the problem, stated in (3), (4), and studied by [KSTZ25, MW26]. Theorem 2 (Informal; see Theorem 5). Let X ∈ Rn×d have i.i.d. entries distributed as N (0, n1 ). If 1.5 ! 1 1 2 d 3 d , n=Ω k + log log + k + log log δ δ δ δ then, for every σ > 0, there is a polynomial-time algorithm that samples from π(· | X, y), defined in (3) and (4), within total variation distance δ, with probability at least 1 − δ over the model. As in prior work, there are two sources of failure in Theorem 2. Namely, the model (3), (4) may fail to produce a sparse signal θ ⋆ or bounded noise ξ (inhibiting tractability of the problem), and the sampling algorithm itself has an approximation error. Our formal result, Theorem 5, is more general and can handle nonuniform inclusion weights in the prior (see Model 3). Moreover, the runtime of Theorem 2 is relatively practical, e.g., it scales linearly in nd. Our n ≈ k 1.5 log2 d requirement improves quadratically in its dependence on k upon the previous state-of-the-art sampler by [KSTZ25], which uses n = Ω(k 3 log3 d) measurements (see also [MW26], who gave a result in the regime n = Ω(d)). Interestingly, the k 1.5 bottleneck appears inherent to our approach (discussed in the following Section 1.2). This motivates the tantalizing open question of whether spike-and-slab posterior sampling is tractable at n = Ω(k log d) measurements, which would close the gap between frequentist and Bayesian SLR.

1.2

Our techniques

In this section, we overview the main proof ideas behind our sparse Dobrushin (Section 3) and trickle down frameworks (Section 4), our annealing reduction from bounded-magnetization to fixedmagnetization sampling (Section 5), and our application to Bayesian SLR (Section 6). Sparse Dobrushin condition. In Section 3, we give a warm-up path coupling analysis of the down-up walk illustrating why restricting k ≪ d can make low-temperature sampling easier. For a measure supported on k-sized subsets, we compare the conditional laws of the up step from neighboring (k − 1)-sized cores, after excluding the coordinates on which the two cores differ. If the resulting total variation discrepancy is O( k1 ), the walk contracts in Hamming distance and mixes in

O k log kε steps (Lemma 3). This framework depends on a sparse variant of the classical Dobrushin influence matrix [Dob68], so we term it a sparse Dobrushin condition (Definition 1). For fixed-magnetization Ising models, this discrepancy is controlled by β maxi̸=j |Jij |, independently of the external field h. Sparsity sets the required discrepancy bound at ≈ k1 , i.e., a relaxed bound at p higher magnetizations. The maximum entry magnitude under the SK model scales as ≈ log d/d, so thisp warm-up result already shows a variant of Theorem 1 at the higher magnetization level k ≤ cβ d/ log d. The sparse Dobrushin condition is simple and broadly applicable, but it cannot exploit cancellation among signed interactions. This motivates our main technique for extending to the range k = Θ(d), based on spectral expansion and the trickle down theorem. Spectral mixing via trickle down. To obtain sharper parameter ranges in fixed-magnetization Ising models, we leverage the trickle down theorem [Opp18, AL20], a foundational result in the U study of high-dimensional expansion. For a measure π over Uk , and a core R ∈ k−2 , define the link graph of R on U \ R to have edge weights Wij := π(R ∪ {i, j}), and let P be the corresponding random walk matrix. The trickle down theorem (Lemma 6) shows that if we can show λ2 (P) = O( k1 ) for all cores R, then the down-up walk satisfies a O( k1 )-Poincaré inequality. Thus, the global mixing problem reduces to proving uniform spectral expansion of all the induced P. For a fixed-magnetization Ising model, each link has a particularly useful form. Fix a core R, and let r ∈ {±1}d be the spin vector whose positive coordinates are R. For some interaction matrix K, Lemma 8 shows that the link weights satisfy, for an appropriate vector a ∈ RU \R , Wij ∝ ai aj exp(Kij ),

Kij := 4βJij

(i ̸= j).

We view W as a bounded perturbation (parameterized by K) of the rank-one weights aa⊤ : when K is the all-zeroes matrix, λ2 (W) ≤ 0 follows simply by a rank argument.3 Our main trickle down framework gives a tighter characterization of λ2 (P) as a function of K. Specifically, Lemma 9 shows using a second-order Taylor expansion of the exponential that if u⊤ Ku ≤ τ ∥u∥1 ∥u∥2 + τ 2 ∥u∥21 ,

(5)

then λ2 (P) = O(τ 2 ). This sets the required bound on τ at ≈ k −1/2 for the trickle down argument. As a point of comparison, the τ 2 ≈ k1 term in this argument, along with u⊤ Ku ≤ (max |Kij |) ∥u∥21 , already qualitatively recovers the sparse Dobrushin condition (Lemma 10). We next show that the mixed norm condition (5) allows more fine-grained control of K, in terms of its sparse operator norms ∥K∥s,op , i.e., the largest ∥KS×S ∥op among any s-sized sets S. This argument proceeds by applying a shelling decomposition, a classic technique from the sparse recovery ∥u∥2 literature [CRT06] that places the “effective sparsity” of a vector u on the scale of ∥u∥21 . 2 √ 2 A basic application of this strategy (Lemma 11) shows that if ∥K∥s,op ≲ τ s + τ s at all scales s ∈ [d], then (5) holds. Plugging in standard bounds on the sparse operator norms of various matrix ensembles now already gives Theorem 4 up to a logarithmic loss in the tolerated k range (Corollary 3), as well as our strongest conclusion for Gaussian Hopfield Ising models (Corollary 4). Our final application to SK models in the linear magnetization regime k = Θ(d) (Theorem 4) uses more fine-grained estimates of sparse operator norms for GOE matrices. 3

Formally, W = aa⊤ − diag (a)2 after removing self-loops, but removing diag (a)2 can only improve λ2 (P).

Annealing. We take a brief detour to discuss a complementary part of our framework: a technique for lifting fixed-magnetization samplers to the bounded-magnetization setting (measures supported on S ⊆ U : |S| ≤ k). Our approach is based on the fact that bounded-magnetization measures are mixtures of fixed-magnetization measures, with weights proportional to normalizing constants. That is, to approximate a global density π on U≤k to total variation δ, it is enough to estimate X Zi := π(ω) for all 0 ≤ i ≤ k ω∈Ui

to multiplicative error O(δ). We formalize this argument in Lemma 14. Conveniently, there is a rich literature in theoretical computer science reducing between the problems of counting (e.g., normalization constant estimation) and sampling. An existing result, Theorem 6 in [Kol18], is essentially black-box applicable to our setting, giving δ-multiplicative estimates to any Zi by using oracle calls to fixed-magnetization samplers. By building upon the estimator of [Kol18], we state our generic reduction from bounded-magnetization sampling to fixed-magnetization sampling in Lemma 15, and give an example of its use for the SK model in Corollary 5. Spike-and-slab posterior sampling. We finally turn to our second main application: Bayesian SLR with a spike-and-slab prior, as studied by [KSTZ25, MW26]. The primary challenge is to correctly sample the support S := supp(θ ⋆ ) ⊆ [d] from the posterior density πsupp (S | X, y) := π({θ : supp(θ) = S} | X, y). As derived in prior work (cf. Fact 3), this density is proportional to a closed-form expression: πsupp (S) ∝

k d−k

|S|

exp

1 ∥bS ∥2A−1 S 2

√

1 , det AS

(6)

where A, b are induced by the measurements (X, y) and defined in (12). Note that the first term in the above expression can be absorbed into an external field. Our posterior sampler has three components. The first follows [KSTZ25] and uses a sparse recovery preprocessing step (Proposition 3) that identifies coordinates whose inclusion is nearly deterministic from (X, y), providing regularity to the residual posterior density. The second applies our annealing procedure from Section 5 to reduce the problem to a fixed-magnetization variant of (6). The remaining step is to use our trickle down framework to demonstrate mixing of the down-up walk on fixed-magnetization slices of the residual posterior. The induced weight matrices W after pinning a core are more complex than in the SK setting, because of the inverse and determinantal terms in (6). In particular, the resulting perturbation K depends on the pinned core through a Schur complement correction, which yields various dependencies. Due to a union bound over all possible cores, our uniform estimate of τ in (5) scales as n−1/2 + nk (e.g., see (53) in Lemma 20), and setting this√to k −1/2 as required by the trickle down framework gives a n ≳ k 1.5 bottleneck. Removing this k factor beyond the measurement complexity of (frequentist) sparse recovery is an exciting problem, that we leave open as a testbed of “average-case” trickle down theorems avoiding the union bound over worst-case dependencies suffered by our approach.

1.3

Related work

Fixed- and bounded-magnetization sampling. The most conceptually relevant prior algorithmic works are by [CDKP22, KPPY25], both of which study fixed-magnetization problems that exhibit improved phase transitions or critical β as the sparsity (positive spin count) k becomes small. For example, Theorem 2 in [CDKP22] shows that for ferromagnetic Ising models with bounded de1 gree ∆, there is a critical “tree threshold” βc (∆) = O( ∆ ) such that for β > βc (∆), the model undergoes a computational phase transition at a certain magnetization level. Notably, the sparsity tradeoff required by their paper beyond β > βc (∆) decays exponentially in β: we give a more formal derivation in Appendix B, but e.g., for constant ∆, they require k = d exp (−Θ∆ (β)) , which only implies a nontrivial setting (k ̸= 0) up to β ≈ log d. Similarly, Theorem 1 of [KPPY25] shows that the Kawasaki dynamics [Kaw66] (a standard magnetization-conservative dynamics) exhibits a phase transition at β > βc (∆), but their magnetization thresholds inherit the same exponential decay in β, and hence also apply only up to β ≈ log d. On the other hand, our results are derived through more general sparsity-aware analysis frameworks (Definition 1) that continue to improve as k → 1, tolerating β up to poly(d). This improvement is closer in spirit to existing results in the sparse recovery literature, where general structural conditions (e.g., RIP) yield sample complexity requirements that improve monotonically with the sparsity level k. There have been other works that study sampling in high-magnetization regimes, that are motivated by giving analysis frameworks for natural sampling dynamics at fixed magnetization (e.g., the down-up walk or Kawasaki dynamics), rather than obtaining improved thresholds as k → 1. For example, [BBD24] study the local Kawasaki dynamics on random ∆-regular graphs, and prove rapid mixing for β = O(∆−1/2 ); notably, their thresholds do not improve as k → 1. Similarly, a line of work has obtained estimates for the mixing time of the down-up walk that improve as k decreases [ALGV19, CGM19, AKV24], e.g., that it mixes in O(k log k) steps for strongly log-concave measures, but their focus was not the relationship between k and the allowable inverse temperature β. High-field strength spin glasses. The Almeida-Thouless line [dAT78] is the canonical benchmark for the replica symmetry phase transition in the SK model with a homogeneous external field θ1d . From a geometric perspective, a recent work of [Lop26] established replica symmetry throughout the AT region, and [KN26] obtained quantitative overlap concentration in the strict AT region. Stronger forms of quantitative overlap concentration in a more restrictive region (cf. Definition 3) have also appeared in the literature [Tal11, JT17]; in particular, we use a bound from [RW26] in our reduction from high-field strength sampling to high-magnetization sampling (Appendix A). On the algorithmic side, [BAR26] prove polynomial-time mixing of the Glauber dynamics in the SK model at sufficiently high field strength, by exploiting replica overlap and sparse operator norm bounds. We quantitatively improve upon the field strength tolerance of [BAR26] to within constant factors of the AT line, or within 1 + o(1) if granted the signs of the mean vector. Other applications of sparsity in sampling. Recent works by [BAR26, DLSS26] both use sparse operator norm bounds to derive mixing times for sampling from Ising models, related to the analytical core of our work (e.g., Lemma 11). Theorem 1 of [BAR26] shows that if the sparse operator norm is bounded at a sparsity level k = Θ(d), then the Glauber dynamics mixes in polynomial time under a sufficiently large external field. In a similar spirit, Section 7 of [DLSS26] 7

uses sparse operator norm bounds to prove rapid mixing of the SK model on small Hamming balls, but only tolerates β < 21 , as opposed to the arbitrary β > 0 handled by our Theorem 1. More generally, a recurring theme in high-dimensional sampling is that sparsity and fixed-cardinality constraints can influence algorithmic landscapes. Recent work on discrete sampling has obtained sharp thresholds and improved mixing guarantees for fixed-size structures such as independent sets and matchings. In particular, [DP23] identified the computational threshold for approximately counting and sampling independent sets of a given size in bounded-degree graphs, [JMPV23] proved optimal O(k log n) mixing of the down-up walk for independent sets of size k, and [JM24] established polynomial-time mixing of the down-up walk for matchings of size k. Recent progress on Ising models has also extended beyond dense mean-field settings: for example, [LMRW24] proves near-linear time mixing of the Glauber dynamics for sparse random Ising models, including the Viana–Bray spin glass, and also treats certain interaction matrices arising from stochastic block models. Bayesian sparse linear regression. Spike-and-slab priors and their induced posteriors are classical tools for Bayesian variable selection [MB88, GM93]. A large statistical literature studies posterior contraction and uncertainty quantification for sparse priors, including spike-and-slab formulations and closely related continuous relaxations [CSHVdV15, Roc18]. Notably, in the moderate signal-tonoise ratio (SNR) regime, the support posterior measure is known to be highly multimodal, which poses an algorithmic challenge [CSHVdV15]. In comparison, the modes collapse at a high SNR, and the posterior nearly reduces to the prior under a low SNR. On the computational side, many works target posterior modes or tractable approximations, including expectation maximization, Lasso-based procedures, and variational Bayes approximations [RG14, RG18, RS22, MS22]. While these works demonstrate the empirical performance of their algorithms, they do not yield end-to-end algorithmic guarantees for the sampling problem we consider. Among provable posterior samplers, [YWJ16] analyzes a Metropolis-Hastings chain under a truncated sparsity prior, which applies only under high or low signal-to-noise ratio regimes (see discussion after their Eq. (10)), while [MW26] gives a polynomial-time sampler via measure decomposition when the number of measurements n grows at least linearly with the ambient dimension d. A complementary diffusion-based approach for linear inverse problems was recently developed in [BH24] (and analyzed in continuous time), though their framework does not seem to directly apply to spike-and-slab posterior sampling at n ≪ d. The closest prior work is [KSTZ25], which gave the first provable spike-and-slab posterior samplers that apply at arbitrary SNRs while allowing n to remain sublinear in d. Our result quadratically improves upon the k dependence of [KSTZ25].

Preliminaries

In this section we develop preliminaries for the rest of the paper. Section 2.1 provides notation used throughout. Section 2.2 gives basic notation and facts about Markov chains used in our analysis. Section 2.3 introduces the main statistical models we consider for our applications.

2.1

Notation

General notation. We use ≲, ≳, and ≈ in informal exposition only to suppress polylogarithmic factors in problem parameters; formal statements have all dependences explicitly stated.

For n ∈ N we let [n] := {i ∈ N : i ≤ n}. We reserve uppercase and lowercase boldface for matrices and vectors respectively. We use 1d and 0d to denote the all-ones and all-zeroes vectors in Rd , Id to denote the identity in Rd , and 0m×n is the m × n all-zeroes matrix. We denote the entrywise (Hadamard) product of equal-length vectors a, b by a ◦ b. Sd×d and Sd×d ⪰0 respectively denote the symmetric and positive semidefinite d × d matrices, where ⪯ is the Loewner partial order. We use nnz to denote the number of nonzero entries of a vector or matrix, and supp denotes the corresponding index set. ei denotes the ith standard basis vector, and we define the row and column selectors M:i := Mei , Mi: := M⊤ ei . For M ∈ Rm×n and (S, T ) ⊆ [m] × [n], MS×T means the ⊤ appropriate submatrix; our convention is to transpose before indexing, so M⊤ T ×S = (MS×T ) . For 1 ≤ p ≤ ∞, ∥·∥p denotes the vector ℓp norm, and for 1 ≤ p, q ≤ ∞ and a matrix M, the associated operator norm is ∥M∥p→q := maxv∈Rd :∥v∥p ≤1 ∥Mv∥q . For any M ∈ Sd×d ⪰0 we define ∥v∥2M := v⊤ Mv. We use ∥·∥F and ∥·∥op to denote the Frobenius and (2 → 2) operator norms. For M ∈ Sd×d we let λ(M) be its eigenvalues sorted so λ1 (M) ≥ . . . ≥ λd (M); we define σ to similarly return the sorted singular values of its input. We define ω < 2.373 [ADV+ 25] so that multiplying, inverting, and eigendecomposing d × d matrices takes O(dω ) time [Str69, PC99]. We frequently use the following “restricted” or “sparse” quantities to parameterize our results: for k ∈ [d], the top-k norm of a vector v ∈ Rd and matrix M ∈ Sd×d are ∥v∥k,p := ∥M∥k,op :=

sup S⊆[d]:|S|=k

sup

∥vS ∥p for all p ≥ 1,

v⊤ Mv =

v∈Rd

sup S⊆[d]:|S|=k

∥MS×S ∥op .

(7)

∥v∥2 ≤1,nnz(v)≤k

For two subsets S, T with the same size of the same universe U, we use ∆Ham (S, T ) ∈ N ∪ {0} to mean the Hamming distance between S, T , which we define as half the number of differing elements. Probability. For a state space Ω, we let P(Ω) denote all probability measures over Ω. For an event E ⊆ Ω, we let I(E) denote the corresponding 0-1 indicator random variable, and µ(E) denote the probability of the event. We let E[·] and Var[·] denote the expectation and variance. For jointly distributed scalar random variables X, Y , Cov(X, Y ) := E[XY ] − E[X]E[Y ] denotes their covariance; if the random variables are instead vector-valued, Cov(X, Y ) is a matrix of appropriate dimension. We abbreviate Cov(X) := Cov(X, X). We frequently use the following distances between µ, ν ∈ P(Ω), where dω denotes a counting measure if Ω is discrete: Z 1 TV (µ, ν) := |µ(ω) − ν(ω)| dω = sup µ(E) − ν(E), 2 Ω E⊆Ω 2 Z µ(ω) − 1 ν(ω)dω. χ2 (µ∥ν) := Ω ν(ω) We denote the set of couplings of µ ∈ P(Ω), ν ∈ P(Ω′ ) (joint measures on Ω × Ω′ whose marginals agree with µ, ν), by Γ(µ, ν). It is standard that when Ω = Ω′ , an alternative definition of the total variation distance is TV (µ, ν) = inf γ∈Γ(µ,ν) P(ω,ω′ )∼γ [ω ̸= ω ′ ] (Proposition 4.7, [LPW09]). We denote the multivariate normal distribution with specified mean and covariance by N (µ, Σ). N We let Bern(p) be the distribution on {0, 1} with EX∼Bern(p) [X] = p. We use i∈[n] πi to denote 9

a product measure with specified marginals, and we use δω to mean a Dirac measure at ω. When Z ∈ RI≥0 are indexed by a set I, we let Multinomial(Z) denote a draw i ∈ I where Law(i) ∝ Z.

2.2

Markov chains

Let T = {Tω }ω∈Ω be a set of transition distributions for a Markov chain on a state space Ω. For an arbitrary measure µ ∈ P(Ω) we let T µ denote the marginal law of ω ′ where ω ∼ µ and ω ′ ∼ Tω . We say that π is a stationary measure for T if T π = π, i.e., Z Tω′ (ω)π(ω ′ )dω ′ = π(ω), for all ω ∈ Ω. We call a Markov chain T reversible if it has stationary measure π, and π(ω)Tω (ω ′ ) = π(ω ′ )Tω′ (ω), for all (ω, ω ′ ) ∈ Ω × Ω. We often consider sampling from discrete measures over the hypercube with fixed-magnetization or P bounded-magnetization. For shorthand we always let X := {±1}, and Xkd := {x ∈ X d : i∈[d] xi = P d −d + 2k}, where the quantity i∈[d] xi = 1⊤ d x is the magnetization. In other words, Xk is the slice of the hypercube X d corresponding to elements with exactly k copies of 1 andSd − k copies of −1. k d d := For the analogous bounded-magnetization problem, we similarly define X≤k j=0 Xj . It is often helpful to associate elements of Xkd with subsets of [d] with size k (the locations of the 1s). We define set(x) ⊆ [d] for x ∈ X d and vec(S) ∈ X d for S ⊆ [d] in the natural way. We frequently consider the (k-)down-up walk Markov chain. This Markov chain can be defined for any measure π supported on Xkd where k ∈ [d]. We denote its transitions by T DU,π and k is inferred from the definition of π. The transition TxDU,π is defined as follows for x ∈ Xkd , where S := set(x). 1. (Down step.) A uniformly random T ⊆ S with |T | = k − 1 is chosen. 2. (Up step.) A set S ′ ⊇ T with |S ′ | = k is sampled ∝ π(S ′ ), and we step to x′ = vec(S ′ ). It is standard that T DU,π is reversible with stationary measure π (Definitions 6 and 7, and Corollary 11, [KO20]). More formal pseudocode is provided in Algorithm 1. Spectral theory. Let T be a reversible Markov chain, and have stationary measure π ∈ P(Ω). We define the associated Dirichlet form by its action on two functions f, g : Ω → R:4 Z ZZ ET (f, g) := f (ω)g(ω)π(ω)dω − f (ω)g(ω ′ )π(ω)Tω (ω ′ )dωdω ′ , and note that the following identity holds: ZZ 1 ET (f, f ) = (f (ω) − f (ω ′ ))2 π(ω)Tω (ω ′ )dωdω ′ . 2

(8)

For λ ∈ (0, 2], we say that T satisfies a λ-Poincaré inequality if, for all f : Ω → R, Varπ [f ] ≤

1 ET (f, f ). λ

For brevity we will omit discussion of integrability issues throughout, but always restrict the functions under consideration to a class where the integration makes sense.

Bounding the Poincaré constant λ results in an estimate of the mixing time of the Markov chain through the comparison inequality χ2 (µ∥π) ≥ 4TV (µ, π)2 , and the following fact. Lemma 1 (Chapters 12 and 13, [LPW09]). Let T be a reversible Markov chain with stationary measure π, and let lazy(T ) denote the Markov chain that in each step transitions according to T with probability 21 , and otherwise does not transition. Then if T satisfies a λ-Poincaré inequality, if we let πt be the law of an iterate taking t ∈ N steps of lazy(T ) starting from π0 , we have λ t 2 2 χ (πt ∥π) ≤ 1 − χ (π0 ∥π). 2 Path coupling. In Section 3, we develop a new tool for bounding mixing times on fixed-magnetization or bounded-magnetization measures. This tool is based on path coupling, which we introduce here. Lemma 2 (Path coupling). Let T be a Markov chain with stationary measure π ∈ P(Ω). Let m : Ω × Ω → N ∪ {0} be an integer-valued metric. Assume that for any ω, ω ′ ∈ Ω with m(ω, ω ′ ) = k, there exists a path ω = ω0 , ω1 , . . . , ωk = ω ′ such that for every i ∈ [k], m(ωi−1 , ωi ) = 1

and

Tωi−1 (ωi ) > 0.

(9)

Assume furthermore that there exists α ∈ (0, 1) such that for any (ω, ω ′ ) ∈ Ω × Ω with m(ω, ω ′ ) = 1, there exists γ ∈ Γ(Tω , Tω′ ) with E(ψ,ψ′ )∼γ [m(ψ, ψ ′ )] ≤ 1 − α. Then for any π0 ∈ P(Ω) and ϵ ∈ (0, 21 ), 1 diam(Ω) T TV T π0 , π ≤ ϵ, for T ≥ log , where diam(Ω) := sup m(ω, ω ′ ). α ϵ ω,ω ′ ∈Ω×Ω Proof. Fix any pair of initial states ω, ω ′ ∈ Ω and let k := m(ω, ω ′ ). Choose a path ω = ω0 , ω1 , . . . , ωk = ω ′ meeting the conditions (9). For all i ∈ [k] let γi ∈ Γ(Tωi−1 , Tωi ) satisfy E(ψ,ψ′ )∼γi m(ψ, ψ ′ ) ≤ 1 − α. From these couplings we can define a joint measure µ over the product space of ψi ∼ Tωi for all 0 ≤ i ≤ k, such that each marginal on adjacent pairs (ωi−1 , ωi ) agrees with γi . This construction is standard and follows from the “gluing lemma” on couplings: see e.g., Lemma 14.3, [LPW09]. Now draw (ψ0 , . . . , ψk ) ∼ µ. We have   X X Eµ [m(ψ0 , ψk )] ≤ Eµ  m(ψi−1 , ψi ) = Eγi [m(ψi−1 , ψi )] ≤ (1 − α)m(ω, ω ′ ). i∈[k]

i∈[k]

This gives a one-step coupling showing a contraction in m. Thus, after T ≥ α1 log diam(Ω) steps, if we ϵ independently draw (ω0 , ω0′ ) ∼ π0 × π, and then iterate the above construction to produce coupled ′ iterates (ωt , ωt′ ) for all t ∈ [T ] such that ωt ∼ Tωt−1 and ωt′ ∼ Tωt−1 , we have E m(ωT , ωT′ ) ≤ (1 − α)T E m(ω0 , ω0′ ) ≤ (1 − α)T diam(Ω) ≤ ϵ. Because m takes values in N ∪ {0}, and we have m(ω, ω ′ ) ≥ 1 whenever ω ̸= ω ′ , I(ωT ̸= ωT′ ) ≤ m(ωT , ωT′ ). The conclusion follows from the coupling definition of total variation, as ωT ∼ T T π0 , ωT′ ∼ T T π. 11

2.3

Statistical models

We describe the statistical models that induce the main structured distributions we consider. Ising model. An Ising model is specified by an interaction matrix J ∈ Sd×d and optionally, an external field h ∈ Rd . It induces a Gibbs measure π over X d , with an additional parameter β > 0 governing the inverse temperature: 1 ⊤ ⊤ π(x) ∝ exp β x Jx + h x · I(x ∈ X d ). (10) 2 We often refer to fixed-magnetization or bounded-magnetization Ising models, where the magnetization level k ∈ [d] is clear from context. Under these models, our goal is to sample from the Gibbs d respectively. measure π in (10), further conditioned on x ∈ Xkd or x ∈ X≤k The literature on Ising models considers a variety of statistical models for the interaction matrix. We summarize a few standard parameterizations here. Model 1 (Sherrington-Kirkpatrick model, [SK75]). In the Sherrington-Kirkpatrick (SK) model, J is drawn from the Gaussian orthogonal ensemble GOE(d): J ∈ Sd×d has Jij ∼i.i.d. N (0, d1 ) for all (i, j) ∈ [d] × [d] with i < j, and Jii ∼i.i.d. N (0, d2 ) for all i ∈ [d].5 Model 2 (Gaussian Hopfield model, [BvEN99, BG08, HRP+ 20]). In the Gaussian Hopfield model,6 J is drawn from the Wishart ensemble with n degrees of freedom: J ∈ Sd×d has J = n1 G⊤ G, where G ∈ Rn×d has Gij ∼i.i.d. N (0, 1) for all (i, j) ∈ [n] × [d]. Our algorithms’ guarantees will depend on an appropriate norm of J (and sometimes, h). Here we present some standard estimates for the random matrix ensembles in Models 1, 2. Fact 1 (Section 2.5, [Ver18] and Theorem 2.3.5, [AGZ10]). For any δ ∈ (0, 21 ), the following hold under Model 1 for an appropriate constant C > 0, simultaneously with probability ≥ 1 − δ. p 1. max(i,j)∈[d]×[d] |Jij | ≤ C log(d/δ)/d. p 2. ∥J∥op ≤ 2 + C log(1/δ)/d. Fact 2 (Lemma 6.26, [Wai19] and Exercise 4.7.3, [Ver18]). For any δ, ϵ ∈ (0, 12 ), the following hold under Model 2 for an appropriate constant C > 0, simultaneously with probability ≥ 1 − δ. 1. If n ≥ C log( dδ ) · ϵ12 , max(i,j)∈[d]×[d] |Jij − I(i = j)| ≤ ϵ. 1 1 2. If n ≥ C(k log( ed k ) + log( δ )) · ϵ2 for k ∈ [d], ∥J∥k,op ≤ 1 + ϵ.

Remark 1. We focus on Gaussian G in Model 2 for simplicity, although Fact 2 holds for any J = n1 G⊤ G where the entries of G are drawn i.i.d. from a 1-sub-Gaussian distribution (see Section 2.5, [Ver18]). Our analyses only rely on the properties of Model 2 in Fact 2, so they apply to any sub-Gaussian ensemble as well. This captures other common Hopfield model instances, e.g., Rademacher G as often considered in the associative memory literature [Lit74, PF77, Hop82]. 5 Different sources parameterize the diagonal of J differently in the SK model, but measures on X d are invariant P to the diagonal of J, as x2i = 1 for all i ∈ [d], xi ∈ X , so i∈[d] Jii x2i is constant. 6 Hopfield models [Hop82] describe Ising models where J is a Gram matrix, and have been enormously influential in machine learning [AGS85, AGS87]. We follow this naming convention for the case when J is Wishart.

Bayesian sparse linear regression. Our second main application considers Bayesian sparse linear regression, i.e., sampling from the posterior distribution of a sparse linear model. We focus on a canonical parameterization of the problem, induced by the spike-and-slab prior [MB88, GM93] (see also [Chi96, Gew96]) and Gaussian measurement noise, a standard formulation recently studied by the sampling algorithms community [KSTZ25, MW26]. Model 3 (Spike-and-slab posterior sampling). Let q ∈ (0, 1)d and σ > 0 be known, and let k̄ := ∥q∥1 . Let X ∈ Rn×d have entries ∼i.i.d. N (0, n1 ), and suppose that we observe (X, y) where O ((1 − qi )δ0 + qi N (0, 1)) , ξ ∼ N (0n , σ 2 In ), y = Xθ ⋆ + ξ, (11) θ ⋆ ∼ π := i∈[d]

and θ ⋆ and ξ are independent. Our goal is to sample from the posterior π(· | X, y). When designing algorithms for Model 3, there is a reparameterization in terms of πsupp the distribution of supp(θ ⋆ ) | X, y. Indeed, the main algorithmic challenge is sampling from πsupp . Fact 3 (Lemma 8, [KSTZ25]). θ ⋆ ∼ π(· | X, y) in Model 3 can be equivalently generated as follows. 1. First, S ⊆ [d] is sampled from πsupp where ! qi 1 1 2 πsupp (S) ∝ ∥bS ∥A−1 √ , exp S 1 − qi 2 det AS i∈S i 1 1 h + IS , bS ∈ RS := 2 X⊤ y. AS ∈ RS×S := 2 X⊤ X σ σ S: S×S Y

(12)

−1 2. Second, θ ⋆ | S, X, y is sampled from N (A−1 S bS , AS ).

Sparse Dobrushin Condition

In this section, we give our first technical result: a simple sufficient condition for rapid mixing of the down-up walk, patterned off of the classical Dobrushin uniqueness condition [Dob68, Wu06]. The applications of our sparse Dobrushin framework to Ising models (Theorem 3 and Corollary 1) in this section generally obtain weaker parameter tradeoffs than our framework in Section 4 does. Nonetheless, this section serves as a useful proof-of-concept of the temperature improvements achievable in fixed-magnetization settings. We include Theorem 3 both due to its ease of applicability, and because the samplers resulting from it are formally incomparable to those in Section 4, as it trades off a faster mixing time for a stricter requirement on relevant parameters. To ease notation, this section works in a more abstract formulation than in Section 2.2, where we use the down-up walk to sample subsets S ⊆ U, from a distribution π supported on Uk := {S ⊆ U : |S| = k} .

(13)

Here, U is a discrete universe of candidate elements. This straightforwardly captures the setting of the down-up walk in Section 2.2 by equating U ≡ [d] and S ≡ x := vec(S). We provide pseudocode implementing a one-step transition of the down-up walk in Algorithm 1. For brevity, we define down(S) := {T ∈ Uk−1 : T ⊂ S} for all S ∈ Uk , 13

up(T ) := {S ∈ Uk : S ⊃ T } for all T ∈ Uk−1 . Algorithm 1: DU(S, U , π) 1 Input: S ∈ Uk , discrete universe U , π ∈ P(Uk ) DU,π 2 Output: Sample S ′ ∈ Uk from TS 3 T ∼unif. down(S) 4 S ′ ∼ π(· | · ∈ up(T )) 5 return S ′

We state our general framework in Section 3.1 and apply it to Ising models in Section 3.2.

3.1

Basic analysis

We begin by stating our new sparse Dobrushin condition. Definition 1 (Sparse Dobrushin condition). Let π ∈ P(Uk ) with full support where U is a discrete universe with |U| ≥ k, and let α ∈ (0, 1). We say π satisfies an α-sparse Dobrushin condition if TV πU ∥V , πV ∥U ≤ α, for all (U, V ) ∈ Uk−1 × Uk−1 with ∆Ham (U, V ) = 1, where for all (U, V ) ∈ Uk−1 × Uk−1 , we define πU ∥V (i) ∈ P(U \ (U ∪ V )) by π(U ∪ {i}) j∈U \(U ∪V ) π(U ∪ {j})

πU ∥V (i) := P

(14)

The utility of Definition 1 reveals itself through the following bound. Lemma 3. Assume that π ∈ P(Uk ) satisfies an α-sparse Dobrushin condition, and let (S, T ) ∈ Uk × Uk satisfy ∆Ham (S, T ) = 1. Then there exists a coupling γ ∈ Γ(TSDU,π , TTDU,π ) satisfying 1 E(S ′ ,T ′ )∼γ ∆Ham (S ′ , T ′ ) ≤ 1 − + α. k Proof. Let W := S ∩ T and assume S = W ∪ {s} and T = W ∪ {t}. We construct the coupling explicitly. First, for the down step (Line 3) we couple the transitions as follows. 1. Draw u ∼unif. W and r ∼unif. [0, 1] independently. 2. If r < k1 , we set S ↓ = T ↓ = W . If r ≥ k1 , we set S ↓ = S \ {u} and T ↓ = T \ {u}. If S ↓ = T ↓ , we can perfectly couple the subsequent up step (Line 4), yielding ∆Ham (S ′ , T ′ ) = 0. In the other case, S ↓ ̸= T ↓ . Let ρ denote the optimal coupling of (πS ↓ ∥T ↓ , πT ↓ ∥S ↓ ) (see Definition 1) inducing their TV distance. Define the marginal probabilities of selecting the disjoint elements as π(S ↓ ∪ {t}) , ↓ j ∈S / ↓ π(S ∪ {j})

π(T ↓ ∪ {s}) . ↓ j ∈T / ↓ π(T ∪ {j})

qS ↓ = P

qT ↓ = P

Without loss of generality (by symmetry of the statement), assume qS ↓ ≥ qT ↓ . For the up step (Line 4) in the case S ↓ ̸= T ↓ , we couple the transitions as follows. 14

1. With probability qT ↓ , set S ′ = T ′ = S ↓ ∪ {t} = T ↓ ∪ {s}. 2. With probability qS ↓ − qT ↓ , set S ′ = S ↓ ∪ {t}, draw j ∼ πT ↓ ∥S ↓ , and set T ′ = T ↓ ∪ {j}. 3. With probability 1 − qS ↓ , draw (j, j ′ ) ∼ ρ, and set S ′ = S ↓ ∪ {j} and T ′ = T ↓ ∪ {j ′ }. We verify that this is a coupling. The first marginal sets S ′ = S ↓ ∪ {t} with probability qS ↓ , and otherwise samples from the correct conditional distribution over U \ T ∪ {u}. Similarly, the second marginal sets T ′ = T ↓ ∪ {s} correctly, and otherwise samples from the correct conditional distribution over U \ S ∪ {t}. Also, in the third case, the probability ∆Ham (S ′ , T ′ ) ̸= 1 is P j ̸= j ′ ≤ α. Finally, we can compute the total expected Hamming distance: 1 ′ ′ E(S ′ ,T ′ )∼γ ∆Ham (S , T ) ≤ 1 − (qT ↓ · 0 + (qS ↓ − qT ↓ ) · 1 + (1 − qS ↓ ) · ((1 − α) + 2α)) k 1 1 ≤ 1− (1 + α) ≤ 1 − + α. k k

(15)

It is clear that the down-up walk over Uk satisfies the condition in (9) with m = ∆Ham . Thus, applying Lemma 3 within the framework of Lemma 2 gives the following result. 1 -sparse Dobrushin condition, and let ϵ ∈ (0, 12 ). Proposition 1. Let π ∈ P(Uk ) satisfy a α ≤ 2k k Then if T = Ω(k log ϵ ) for a sufficiently large constant, we have for any π0 ∈ P(Uk ), TV (T DU,π )T π0 , π ≤ ϵ.

Of course, the parameter α in Proposition 1 can be taken to be any constant factor smaller than k1 . Proposition 1 also explains our naming choice, as its conclusion becomes stronger (as a function of the sparse Dobrushin condition parameter) as k gets smaller, i.e., π is supported on sparser sets.

3.2

Fixed-magnetization Ising models

In this section, we demonstrate how to apply Proposition 1 to fixed-magnetization Ising models: 1 ⊤ ⊤ π(x) ∝ exp β x Jx + h x · Ix∈X d . (16) k 2 We first require a helper tool to control the sparse Dobrushin condition parameter. Lemma 4. Let π ∈ P(Ω), µ ∈ P(Ω) have π ∝ P and µ ∝ Q for unnormalized densities P, Q. Then sup log ω∈Ω

P (ω) ≤ ∆ =⇒ TV (π, µ) ≤ ∆. Q(ω)

Proof. First observe that R Q(ω ′ )dω ′ P (ω) P (ω) π(ω) ≤ sup log + log RΩ ≤ 2 sup log ≤ 2∆. sup log ′ ′ µ(ω) Q(ω) Q(ω) ω∈Ω ω∈Ω ω∈Ω Ω P (ω )dω Let L(ω) := π(ω) µ(ω) so exp(−2∆) ≤ L(ω) ≤ exp(2∆) for all ω ∈ Ω. Using π = L µ, 1 TV (π, µ) = 2

1 |π(ω) − µ(ω)|dω = 2 Ω

Z µ(ω) |L(ω) − 1| dω = Ω

1 Eµ |L − 1| . 2

Next we use the following convexity bound: if ϕ is convex on [a, b] and random variable X ∈ [a, b], then writing X = ta + (1 − t)b with t = b−X b−a ∈ [0, 1] gives ϕ(X) ≤ tϕ(a) + (1 − t)ϕ(b). Taking expectations yields E[ϕ(X)] ≤

E[X] − a b − E[X] ϕ(a) + ϕ(b). b−a b−a

We apply this bound with X = L, Eµ [L] = 1, ϕ(u) = |u − 1|, and [a, b] = [exp(−2∆), exp(2∆)]. Since a < 1 < b, ϕ(a) = 1 − a and ϕ(b) = b − 1, so 1 (b − 1)(1 − a) TV (π, µ) = Eµ [ϕ(L)] ≤ . 2 b−a Now set r := ab > 1 and write b = ar. Since a ≤ 1 ≤ ar, we have a ∈ [ 1r , 1]. Define fr (a) :=

(ar − 1)(1 − a) (ra − 1)(1 − a) = . ar − a a(r − 1)

A direct derivative computation shows that fr is maximized at a = r−1/2 , and hence √ r−1 1 TV (π, µ) ≤ fr (r−1/2 ) = √ = tanh log r ≤ ∆. 4 r+1

Lemma 4 lets us conclude that for fixed-magnetization Ising models (i.e., (10) restricted to Xkd ), the sparse Dobrushin condition parameter can be controlled by the largest off-diagonal entry of J. Lemma 5. Let π be a fixed-magnetization Ising model (16) with β

max (i,j)∈[d]×[d] i̸=j

Then π satisfies a 8α-sparse Dobrushin condition.

|Jij | ≤ α.

Proof. Write J ← βJ for simplicity in this proof, so we will prove the result when β = 1 without loss d ×Xd of generality. Let (vec(U ), vec(V )) ∈ Xk−1 k−1 have ∆Ham (U, V ) = 1, and denote W := U ∩ V , U = W ∪ {u}, and V = W ∪ {v}. Note that the distributions πU ∥V and πV ∥U defined in (14) are supported on the same set Ω := [d] \ (U ∪ V ), so we are in the setting of Lemma 4. Let i ∈ [d] \ (U ∪ V ), X := U ∪ {i}, Y := V ∪ {i}, and x := vec(X) = 21X − 1d , y := vec(Y ) = 21Y − 1d . Observe that by expanding definitions and cancelling similar terms, ! exp( 21 x⊤ Jx + h⊤ x) 1 ⊤ 1 ⊤ ⊤ ⊤ log = x Jx + h x − y Jy + h y 2 2 exp( 12 y⊤ Jy + h⊤ y) 1 ⊤ ⊤ ⊤ 1 J1 − h 1 = 21⊤ J1 + 2 (h − J1 ) 1 + X X d d d X 2 d 1 ⊤ ⊤ ⊤ ⊤ − 21Y J1Y + 2 (h − J1d ) 1Y + 1d J1d − h 1d 2 = 4 (ei + eW )⊤ J(eu − ev ) + 2 (Juu − Jvv ) + 2(h − J1d )⊤ (eu − ev ) ⊤ = 4e⊤ i J(eu − ev ) + 2(Juu − Jvv ) + 2 (h + J(21W − 1d )) (eu − ev ).

Further, we have 1 ⊤ 1 ⊤ ⊤ ⊤ ⊤ x Jx + h x ∝ exp x Jx + h x − 2Juu − 2(h + J(21W − 1d )) eu , πU ∥V (X) ∝ exp 2 2 1 ⊤ 1 ⊤ ⊤ ⊤ ⊤ y Jy + h y ∝ exp y Jy + h y − 2Jvv − 2(h + J(21W − 1d )) ev . πV ∥U (Y ) ∝ exp 2 2

Above we used that W and u are fixed in the definition of πU ∥V (·), so we can bring terms involving only these indices into the unnormalized density; a similar argument holds for πV ∥U (·). Thus, Lemma 4 applies with (P, Q) set to the right-hand sides above, so its conclusion holds with ∆=

4e⊤ i J(eu − ev ) ≤ 8α.

max (u,v,i)∈[d]×[d]×[d] (u,v,i) distinct

Importantly, the bound in Lemma 5 is independent of h. Intuitively, this follows because the conditional distributions in (14) exclude all elements in U ∪ V , including the non-shared elements, which induce the only difference in the external field. By combining Lemma 5 and Proposition 1, we thus obtain a mixing time bound for fixed-magnetization Ising models. Theorem 3. Let δ ∈ (0, 12 ) and let π be induced by a fixed-magnetization Ising model (16) satisfying β

max (i,j)∈[d]×[d] i̸=j

|Jij | ≤

1 . 16k

Then if T = Ω(k log kδ ) for a sufficiently large constant, we have for any π0 ∈ P(Xkd ), TV (T DU,π )T π0 , π ≤ δ. 17

(17)

We conclude by briefly stating example implications of Theorem 3 and Lemma 5 for sampling from the Gibbs distributions induced by Models 1 and 2. Corollary 1. Let δ ∈ (0, 12 ) and let π be defined as in (16). If π0 ∈ P(Xkd ) and T = Ω(k log kδ ) for an appropriate constant, TV((T DU,π )T π0 , π) ≤ δ, under any of the following conditions. p 1. Under Model 1 with probability ≥ 1 − δ, if β = O( k1 d/ log(d/δ)) for an appropriate constant. 2. Under Model 2 with probability ≥ 1 − δ, if n = Ω(max(1, (βk)2 ) log dδ ) for an appropriate constant. Proof. It suffices to use high-probability bounds on the maximum off-diagonal entry of J (via Facts 1 and 2), combined with Theorem 3 and Lemma 5. We give a brief discussion of the runtime of Algorithm 1 in the Ising model setting. Remark 2 (Runtime of down-up walk for Ising model). For the Ising model, Algorithm 1 can be implemented in time O(d) per iteration, after O(d2 ) time preprocessing. We believe this is folklore, but sketch a proof here. Our implementation maintains a set S ⊆ [d] and applies Algorithm 1 to it. Clearly, Line 3 is implementable in O(k) time. Next, the law of {i} = S ′ \ T in Line 4 is ⊤ β ∝ exp 21T ∪{i} − 1d J 21T ∪{i} − 1d + β h, 21T ∪{i} − 1d 2 ∝ exp (β (2Jii + 4 ⟨Jei , 1T ⟩ − 2 ⟨Jei , 1d ⟩ + 2hi )) . Therefore, it is enough to maintain the quantities, for all i ∈ [d] \ T , X ⟨Jei , 1T ⟩ = Jij , Jii , ⟨J1d , ei ⟩ , hi . j∈T

After O(d2 ) time preprocessing, we store the latter three (constant) quantities, and we can update the first in O(d) time per call to Algorithm 1, since at most 2 coordinates change in T . Finally, given these values the sampling on Line 4 can be performed in time O(d).

Spectral Mixing via Trickle Down

In this section, we develop a second approach to prove mixing bounds on fixed-magnetization measures, based on the trickle down theorem [Opp18, AL20] from the literature on high-dimensional expanders (we recommend Section 5 of the excellent survey [GK23] as an introduction to this topic). This framework, summarized abstractly in Section 4.1, achieves tighter parameter tradeoffs in our applications to Ising models than its sparse Dobrushin counterpart in Section 3. In Section 4.2, we begin by giving a more interpretable sufficient condition for our framework, and two applications, as warmups. Specifically, we show that our trickle down framework qualitatively subsumes the sparse Dobrushin condition, and implies fast mixing for the SK model (Model 1) in the near-proportional magnetization regime, k = Oβ ( logd d ). We also derive an analogous result for the Gaussian Hopfield model (Model 2). Finally, we conclude with our strongest result on the SK model, handling the proportional regime k = Oβ (d), in Section 4.3. 18

4.1

Trickle down framework

We develop our framework for analyzing fast mixing on fixed-magnetization Ising models in three parts. We begin by recalling preliminaries on spectral graph theory, and a statement of the trickle down theorem of [Opp18, AL20]. We then state a sufficient condition (21) for applying the trickle down theorem, when the weights of the distribution in question are governed by a small perturbation of a product graph. Finally, we specialize this framework to Ising models. Trickle down. Let [d] index a finite vertex set. For an edge weight matrix W ∈ Rd×d ≥0 , we define its associated random walk matrix P to be the degree-normalized W, i.e., P := D−1 W, where D := diag (W1d ) .

(18)

This section only considers reversible P, associated with symmetric edge weight matrices W. We next state the trickle down theorem [Opp18, AL20], which bounds the Poincaré constant of the down-up walk induced by a measure in P(Xkd ), in terms of the worst spectral gap among certain restrictions of the measure. It is proven by using the law of total variance, which gives recursive relationships among the spectral gaps of various restricted down-up walks. We defer additional background to [Opp18, AL20], and simply state a sufficient form for our purposes here. For a measure π ∈ P(Xkd ) and a set R ⊆ [d] with |R| = k − 2, define the link graph of R to be the weighted graph on vertices [d] \ R with edge weight matrix W given by Wij = π(R ∪ {i, j}) for all (i, j) ∈ ([d] \ R) × ([d] \ R), i ̸= j. Here we associate the set R ∪ {i, j} with an element of Xkd with those positive coordinates, per our convention. We denote the associated random walk matrix by PR , following (18). Lemma 6 (Theorem 2.5, [Opp18] and Theorem 3.1, [AL20]). For π ∈ P(Xkd ) with full support, if α [d] λ2 (PR ) ≤ for all R ∈ , k−1 k−2 for some α < 1, then T DU,π satisfies a λ-Poincaré inequality, for λ = 1−α k . Spectral gap for rank-one perturbations. To apply Lemma 6, we require tools for bounding [d] the spectral gaps of the link graphs induced by each R ∈ k−2 . The next piece of our framework, Lemma 7, shows such a spectral gap for graphs where the edge weight matrix W is induced by an appropriately-bounded perturbation K of a rank-one matrix aa⊤ . Remark 3. An instructive warmup is when K is the all-zeroes matrix in Lemma 7, in which case W = aa⊤ − diag (a)2

(19)

is a rank-one matrix with its diagonal removed. Because λ2 (aa⊤ ) = 0 and diag (a)2 ∈ Sd×d ⪰0 , the min-max characterization of eigenvalues gives λ2 (W) ≤ 0. The same strategy, applied to the similar matrix D−1/2 WD−1/2 , implies λ2 (P) ≤ 0, following the notation (18). Lemma 7 robustly extends this bound to the case where W in (19) is perturbed by a bounded matrix K.

Lemma 7. Let a ∈ Rd>0 , and let K ∈ Sd×d satisfy Kii = 0 for all i ∈ [d]. Define ( ai aj exp (Kij ) i ̸= j Wij = , for all (i, j) ∈ [d] × [d]. 0 i=j

(20)

Also, let m(K) := max(i,j)∈[d]×[d] |Kij |, let X ∈ Sd×d have Xij = exp(Kij ) − 1 entrywise, and let ! u⊤ (X − Id ) u . (21) ρ(K) := max 0, sup 2 2 u∈Rd :nnz(u)>1 ∥u∥1 − ∥u∥2 Then, following the notation (18), λ2 (P) ≤ exp(m(K))ρ(K). Proof. Let N := D−1/2 WD−1/2 , so that λ2 (P) = λ2 (N). We first write N as the sum of a rank-one matrix and a correction. For A := diag (a), 1

N = D− 2 aa⊤ D− 2 + D− 2 A (X − Id ) AD− 2 . Because the min-max characterization of eigenvalues gives 1 1 1 1 λ2 (N) ≤ λ2 D− 2 aa⊤ D− 2 + λ1 D− 2 A (X − Id ) AD− 2 , it is enough to bound the second term above. We next have 1

λ1 D

− 21

− 12

A (X − Id ) AD

y⊤ D− 2 A(X − Id )AD− 2 y = sup ∥y∥22 y∈Rd :y̸=0d u⊤ (X − Id )u = sup . ⊤ −2 u∈Rd :u̸=0d u A Du

(22)

If the supremum is achieved by a 1-sparse u, then the lemma statement holds, because ρ(K) ≥ 0, and any 1-sparse u has a negative numerator in (22), because X has an all-zeroes diagonal. We now bound the denominator of (22) when nnz(u) > 1:   X u2 Dii X u2 X i  i u⊤ A−2 Du = = aj exp(Kij ) ai a2i i∈[d]

≥ exp (−m(K))

i∈[d]

j∈[d]:j̸=i

X u2 (∥a∥ − ai ) i

i∈[d]

 = exp (−m(K)) 

X u2 ∥a∥

i∈[d]



− ∥u∥22  ≥ exp (−m(K)) ∥u∥21 − ∥u∥22 .

The last line applied the Cauchy-Schwarz inequality. Plugging this into (22) gives the result. As a simple application of Lemma 7, we rederive a standard mixing result on the Curie-Weiss model. Model 4 (Curie-Weiss model, [Wei07, Ell12]). In the Curie-Weiss model, J = d1 1d 1⊤ d. 20

Corollary 2. Let π ∈ P(Xkd ) be induced by a fixed-magnetization Ising model (16), under the Curie-Weiss model (Model 4). Then for any β ∈ R, T DU,π satisfies a k1 -Poincaré inequality. ⊤ Proof. Recall that the Curie-Weiss interaction matrix is J = d1 1d 1⊤ d . Because 1d x is constant over x ∈ Xkd , identifying each x with its set S of positive coordinates, we have π(S) ∝ exp 2βh⊤ 1S . (23) [d] Now for the random walk matrix PR associated with the link graph of R ∈ k−2 , we have that the associated edge weights W follow (20) with K set to the all-zeroes matrix, and

ai = exp (2βhi ) for all i ∈ [d] \ R. Therefore Lemma 7 applies with K = 0 and gives a bound of α = 0 for use with Lemma 6. Corollary 2 rephrases the following proof: the fixed-magnetization Curie-Weiss model is a product distribution restricted to a Hamming slice (23), which Theorem 1.1, [ALGV19] proves a Poincaré inequality for. We include this example to illustrate a trivial case of our framework. Specialization to Ising model. We next derive a generic application of the framework given by Lemmas 6 and 7 to Ising models. Consider a fixed-magnetization Ising model (16), where 1 ⊤ ⊤ π(x) ∝ exp β x Jx + h x · Ix∈X d . k 2 Ising models are invariant to changes in the diagonal of J, so without loss of generality, we explicitly assume in this section that J has zero diagonal, i.e., Jii = 0 for all i ∈ [d]. Lemma 8. For a fixed-magnetization Ising model π (16), and following the notation (21), if ρ(4βJ) exp(m(4βJ)) ≤

1 , 2(k − 1)

1 then T DU,π satisfies a 2k -Poincaré inequality.

Proof. Following the notation of Lemma 6, it is enough to show that for all R ∈ λ2 (PR ) ≤

[d] k−2 ,

1 . 2(k − 1)

We prove this using Lemma 7. Fix some R for the remainder of this proof, and let S := [d] \ R. Let d r ∈ Xk−2 have positive coordinates with indices R, and consider some x = r + 2(ei + ej ) ∈ Xkd , i.e., corresponding to the set R ∪ {i, j}. We have 1 ⊤ 1 ⊤ ⊤ ⊤ β x Jx + h x = β r Jr + h r + 2β (h + Jr)⊤ (ei + ej ) + 4βJij . 2 2 The first term is a constant for all pairs (i, j) ∈ S × S. Therefore, up to a proportionality constant, the edge weights are given by (20), where for all (i, j) ∈ S × S, ai = exp 2β (h + Jr)⊤ ei , Kij = 4βJij . The conclusion now follows from Lemma 7 and the assumption, because excluding R from the coordinates can only decrease both m(4βJ) and ρ(4βJ). 21

4.2

Simple sufficient conditions for fast mixing

Section 4.1 gives a generic strategy for sampling in fixed-magnetization Ising models. By combining Lemmas 6, 7, and 8, our task reduces to bounding ρ(K) and m(K) for K ← 4βJ. The bottleneck is typically to control ρ(K), whose definition (21) is somewhat opaque. In Lemma 9, we give a more interpretable sufficient condition for applying this framework. We show that one specialization of this condition qualitatively recovers the sparse Dobrushin condition, and that another implies improvements in the same applications as considered in Corollary 1. Mixed-norm quadratic form bound. Our first strategy for controlling ρ(K) decomposes Xij = exp(Kij ) − 1 into a linear term in Kij , and a high-order term. The high-order contribution to the quadratic form in X is folded into our assumption (24), which implies a bound on ρ(K). Lemma 9. In the setting of Lemma 8, define K := 4βJ. Then, if u⊤ Ku ≤ τ ∥u∥1 ∥u∥2 + τ 2 ∥u∥21 for all u ∈ Rd ,

(24)

for some τ ∈ [0, 14 ], we have ρ(K) exp(m(K)) ≤ 9τ 2 . Proof. First, note that (24) implies a bound on m(K): taking u ← ei + ej gives √ |Kij | ≤ 2τ + 2τ 2 ≤ 2τ. Therefore, m(K) ≤ 2τ . Further, if we decompose X = K + R, then √ 0 ≤ Rij ≤ sup exp(x) − 1 − x ≤ 2 eτ 2 , for all (i, j) ∈ [d] × [d]. |x|≤2τ

Now by the triangle inequality, and Young’s inequality applied to (24), u⊤ Xu ≤ u⊤ Ku + u⊤ Ru √ 1 1 1 2 ≤ ∥u∥2 + 1 + + 2 e τ 2 ∥u∥21 ≤ ∥u∥22 + 5τ 2 ∥u∥21 . 2 2 2 We thus have

1 u⊤ Xu − ∥u∥22 ≤ − ∥u∥22 + 5τ 2 ∥u∥21 ≤ 5τ 2 ∥u∥21 − ∥u∥22 . 2 Therefore, ρ(K) ≤ 5τ 2 , and the conclusion follows from supτ ∈[0, 1 ] exp (2τ ) ≤ 59 . 4

Recovering a sparse Dobrushin condition. We next observe that the τ 2 ∥u∥21 term alone in (24) already qualitatively recovers the sparse Dobrushin condition of Theorem 3. 1 1 Lemma 10. In the setting of Lemma 8, if m(βJ) ≤ 72k , T DU,π satisfies a 2k -Poincaré inequality. 1 Proof. By combining Lemmas 8 and 9, it suffices to show that (24) holds with τ 2 = 18k . Under the assumption on m(βJ), we have the desired

u⊤ Ku ≤ m(K) ∥u∥21 = m(4βJ) ∥u∥21 ≤

1 ∥u∥21 . 18k

We remark that compared to Theorem 3, Lemma 10 loses a constant factor in the allowable temperature range, and only implies mixing in χ2 divergence, which typically loses a ≈ k factor compared to analogous mixing time bounds in Hamming distance. Applications to Models 1 and 2. Our second application of Lemma 9 controls the τ required in (24) via the largest ∥KS×S ∥op , appropriately normalized by |S|, over all S ⊆ [d]. As intuition for why, if u is the 0-1 indicator vector for some S, ∥KS×S ∥op |u⊤ Ku| 1 p = , ≤ |S| ∥KS×S ∥op · ∥u∥1 ∥u∥2 |S|1.5 |S| ∥KS×S ∥op |u⊤ Ku| 1 . 2 ≤ |S| ∥KS×S ∥op · |S|2 = |S| ∥u∥1 Lemma 11 uses a shelling decomposition to make this intuition rigorous. Interestingly, the shelling decomposition is a standard strategy for passing to continuous notions of sparsity [CRT06], further strengthening connections between our paper’s toolkit and the broader literature. Lemma 11. Let K ∈ Sd×d , and suppose that √ ∥K∥s,op ≤ α s + α2 s for all s ∈ [d].

(25)

Then (24) holds with τ := 32α. Proof. Fix a vector u ̸= 0d , and let

& σ :=

∥u∥21

∥u∥22

which can intuitively be thought of as a numerical analog of the sparsity of u. Sort the coordinates of u (relabel [d] by a permutation) so that |u1 | ≥ |u2 | ≥ . . . ≥ |ud |. Now partition [d] into consecutive blocks B1 = [σ], B2 = [2σ] \ B1 , . . . of size at most σ each. For each ℓ ≥ 2, monotonicity of the coordinates gives ∥uBℓ ∥2 ≤

√

1 σ ∥uBℓ ∥∞ ≤ √ uBℓ−1 1 , σ

so summing, we have X ℓ

1 ∥uBℓ ∥2 ≤ ∥u∥2 + √ ∥u∥1 ≤ 2 ∥u∥2 . σ

Now we decompose u⊤ Ku blockwise, and control each block using α: because each union B

ℓ ∪ B ℓ′

is 2σ-sparse, and each blockwise contribution is only supported on this set, X √ u⊤ Ku ≤ 2(α σ + α2 σ) ∥uBℓ ∥2 uBℓ′ 2 ℓ,ℓ′

√

!2 2

≤ 2(α σ + α σ)

∥uBℓ ∥2

ℓ

Finally, using ∥u∥2 ≤ 2σ −1/2 ∥u∥1 gives the claim. 23

≤ 8(α

√

(26) σ + α σ) ∥u∥22 . 2

√ We now derive an application to the SK model (Model 1). This application only uses the α s term in (25), due to the operator norm behavior of entrywise sub-Gaussian matrices. Per our convention in this section, we use Jii = 0 for all i ∈ [d], which does not affect the Gibbs measure (16). Corollary 3. Let δ ∈ (0, 21 ) and let π be induced by a fixed-magnetization Ising model (16) under the SK model (Model 1). Further, assume that for an appropriate constant, ! s d . β=O k log dδ 1 Then with probability ≥ 1 − δ, T DU,π satisfies a 2k -Poincaré inequality.

Proof. We first claim that for Model 1 with zero diagonal, s ∥J∥s,op log dδ ≤C max √ , d s s∈[d]

(27)

with probability ≥ 1 − δ, for a universal constant C. To see this, fix s ∈ [d]. With probability [d] δ ≥ 1 − ds+1 , Corollary 3.9, [BvH16] shows that for a fixed S ∈ s , s ∥JS×S ∥op ≤ C

s log dδ . d

(28)

Now a union bound over the ≤ ds possible S, and the d possible s ∈ [d], shows (27). Finally, for K = 4βJ, the assumed range on β implies that, following the notation (25), 32α ≤ 1 (18k)−1/2 . Also, Lemma 11 implies that (24) holds with τ = 32α. Combining gives 9τ 2 ≤ 2k in Lemma 9, and then the claim follows from Lemma 8. Corollary 3 directly improves the allowable temperature range in Corollary 1’s SK model special√ ization by a k factor. Unfortunately, it does not permit taking arbitrary β = O(1) unless k is sufficiently sublinear in d, i.e., smaller than logd d . This is an inherent artifact of using the bound √ (27), because the maximum of Ω(d) Gaussians grows with log d, causing an obstruction at |S| = 2. However, it is not inherent to the SK model, and in Section 4.3, we show how to further shave this extraneous logarithmic factor, by more directly controlling ρ(4βJ) in Lemma 8. We conclude the section with a similar improvement upon Corollary 1’s Gaussian Hopfield model specialization. Crucially, to be compatible with the two-regime sub-exponential concentration of the operator norms of Wishart matrices, our proof uses both of the terms in (25). Corollary 4. Let δ ∈ (0, 21 ) and let π be induced by a fixed-magnetization Ising model (16) under the Gaussian Hopfield model (Model 2). Further, assume that for an appropriate constant, d 2 n = Ω max(β, β )k log . δ 1 Then with probability ≥ 1 − δ, T DU,π satisfies a 2k -Poincaré inequality.

Proof. We claim that for Model 2 with zero diagonal, for each fixed S ⊆ [d] with |S| = s, s  d d s log δ s log δ , ∥JS×S ∥op ≤ C  + n n

(29)

δ with probability ≥ 1 − ds+1 , for a universal constant C. At this point, the proof follows identically to that of Corollary 3, using the tighter bound in (25). To see that (29) holds, let J = n1 G⊤ G − D for G ∈ Rn×d where the Gij are i.i.d. Gaussian, and D is a diagonal matrix agreeing with the diagonal of n1 G⊤ G. Then (29) follows because with probability ≥ 1 − δ, s  1 1 s + log δ s + log δ 1 ⊤ , [G G]S×S − IS = O + n n n op (30) ! r s s log δ log δ + , ∥DS×S − IS ∥op = O n n

where the first bound above uses Exercise 4.7.3 of [Ver18], and the second uses Theorem 3.1.1 of [Ver18] with a union bound over all s of the diagonal coordinates. Now (29) follows by substituting δ δ ← ds+1 above and applying the triangle inequality.

4.3

Linear magnetization in low-temperature SK models

We conclude with a tighter analysis of the quantities required by Lemma 7 for the SK model, that removes the extraneous logarithmic factor from Corollary 3. The results in this section hold assuming a high-probability event under Model 1, captured in the following lemma. Lemma 12. Let δ ∈ (0, 12 ) and C > 0 be a sufficiently large universal constant. Then the following events simultaneously hold with probability ≥ 1 − δ over Model 1. 1. For all S ⊆ [d] with |S| = s ∈ [d], s ∥JS×S ∥op ≤ C 

2. We have

s max (i,j)∈[d]×[d]

|Jij | ≤ C

 d s log ed + log s δ . d

log dδ , d

r ∥J1d ∥∞ ≤ C

d log . δ

∥u∥2

3. For all nonzero u ∈ Rd with q := ∥u∥22 , 1

r u⊤ Ju ≤ C ∥u∥21 

s q log(edq) +q d



log dδ  d

4. Let H := J ◦ J − d1 (1d 1⊤ d − Id ). Then, s ∥H∥op ≤ C 

log dδ d



log dδ  d

Proof. We allot a 3δ failure probability for Items 1, 2, and 4, and Item 3 will follow from Item 1. Item 1 follows from the calculation in (27), using the tighter estimate ds = exp(O(s log ed s )). Both parts of Item 2 follow from standard bounds on the maximum of poly(d) i.i.d. Gaussians, where the variance of each entry of J1d is at most 1. Item 3 follows from the same shelling decomposition argument as in Lemma 11. Concretely, perform the same decomposition into blocks B1 , B2 , . . . of size at most σ = ⌈ 1q ⌉. Then the same argument as in (26), combined with the estimate on ∥J∥2σ,op already derived in Item 1, yields s u⊤ Ju ≤ 8C 

s log(edq) + qd

 log dδ  ∥u∥2 . 2 d

The claim then follows by adjusting the constant C, and substituting the definition of q. There is an edge case when 2σ ≥ d, but in this case ∥J∥op satisfies the required bound (see Fact 1). Finally, Item 4 asks to bound the operator norm of a matrix with i.i.d. sub-exponential entries. P 2 − 1 )E . ⊤ for each 1 ≤ i < j ≤ d. Then H = (J Concretely, let Eij := ei e⊤ + e e ij j i j 1≤i<j≤d ij d A straightforward calculation shows that the sub-exponential matrix Bernstein inequality (e.g., Theorem 6.2 in [Tro12] with σ 2 = R = O( d1 )) now applies, which concludes the proof. We now use Items 2, 3, and 4 of Lemma 12 to derive estimates on the parameters ρ(K) and m(K) defined in Lemma 7, when K = 4βJ as derived in Lemma 8. Lemma 13. Assume the success of the events in Lemma 12, and let β ∈ R+ , β̄ := max(β, 1), and K := 4βJ. Then for a sufficiently large universal constant C ′ > 0, assuming s log dδ log dδ 1 ∆ := + ≤ ′ 2 , (31) d d C β̄ log(eβ̄) we have ρ(K) ≤

C ′ β̄ 2 log(eβ̄) , d

m(K) ≤ C ′ β̄ 2 ∆.

Proof. The bound on m(K) follows from the first condition in Item 2 of Lemma 12, and any C ′ ≥ 4C.

To bound ρ(K), we follow the notation of Lemma 7, so Xij = exp(Kij ) − 1 entrywise. We also ∥u∥2

decompose X = K + R as in Lemma 9. For the linear term, letting u ∈ Rd have q := ∥u∥22 , 1

r



log dδ  d C ′′ β̄ 2 log(eβ̄) 2 q , ≤ ∥u∥1 + 4βCq∆ + 4 d

u⊤ Ku = 4β u⊤ Ju ≤ 4βC ∥u∥21 

q log(edq) +q d

(32)

where the second line used the scalar inequality, for an appropriate C ′′ depending on C, p p s 4βC s log(es) ≤ 4β̄C s log(es) ≤ + C ′′ β̄ 2 log(eβ̄), 4 valid for any s, β̄ ≥ 1. For the residual term, we first establish the entrywise bound 0 ≤ exp(Kij ) − Kij − 1 = Rij ≤ K2ij ≤ 16β 2 J2ij , as we have shown m(K) ≤ C ′ β̄ 2 ∆ ≤ 1 already. Thus, letting wi := |ui | for all i ∈ [d], u⊤ Ru ≤ 16β 2 w⊤ (J ◦ J) w 16β 2 ⊤ w − I w 1d 1⊤ (33) d d d 16β 2 ∥u∥21 ≤ 16β 2 C∆ ∥u∥22 + d where the second line used the definition of H from Item 4, and the last line applied Item 4 and our definition of ∆. Finally, by combining (32) and (33), ′′ 2 3 C β̄ log(eβ̄) + 16β 2 2 2 ⊤ 2 u Xu − ∥u∥2 ≤ − + 4βC∆ + 16β C∆ ∥u∥2 + ∥u∥21 4 d C ′ β̄ 2 log(eβ̄) 1 ∥u∥21 , ≤ − ∥u∥22 + 2 d = 16β 2 w⊤ Hw +

′ 2

β̄) for an appropriate C ′ . We now obtain the desired bound on ρ(K) by noting C β̄ log(e ≤ 12 . d

Theorem 4. Let δ ∈ (0, 12 ), and let π be induced by a fixed-magnetization Ising model (16), where d k=O β̄ 2 log(eβ̄) for a sufficiently small constant, defining β̄ := max(1, β). Also, assume that the condition (31) holds. With probability ≥ 1 − δ over the SK model (Model 1), if !!! r 1 d T = Ω k log + k 2 log(d) + β ∥h∥∞ + log δ δ for a sufficiently large constant, we have for any π0 ∈ P(Xkd ), T DU,π TV lazy(T ) π0 , π ≤ δ. 27

Proof. The failure probability comes from Lemma 12, so henceforth condition on its success. Com1 bining Lemmas 7, 8, and 13 implies that T DU,π satisfies a 2k -PI, for the assumed range on k. 2 Lemma 1 now shows the χ divergence of the lazy down-up walk contracts by a factor of Ω( k1 ) in each iteration. The conclusion follows if we can bound the initial χ2 divergence: 1 π0 2 2 − 1 ≤ 2 , where πmin := min π(x). χ (π0 ∥π) = Eπ π πmin x∈Xkd Thus it remains to control πmin . Identify x ∈ Xkd with S ∈

[d] k , so x = 21S − 1d . Then for

⊤ ⊤ F (x) = 21⊤ S J1S − 21S J1d + 2h 1S ,

the Ising measure is ∝ exp(βF ). On the events in Items 1 and 2 of Lemma 12, ! r d 21⊤ 21⊤ , S J1d = O k log S J1S = O (k(1 + ∆)) = O(k), δ where we used that (31) gives ∆ = O(1). Hence, by bounding the range of βF over Xkd , !!! r 1 d πmin ≥ [d] exp −O k + βk ∥h∥∞ + log δ k !!! r d = exp −O k log(d) + βk ∥h∥∞ + log . δ 1 ) and simplifying. The conclusion follows by taking T = Ω(k log πmin

We make two brief remarks on the statement of Theorem 4. First, the condition (31) is a lower bound on the failure probability δ that must apply for Theorem 4 to hold. This condition is relatively mild, and even permits taking δ = exp(−Ω(d)) for appropriate constants. Second, the use of β̄ = max(β, 1) in the statement makes the k upper bound uniform over all β ≥ 0. In particular, if the assumptions hold at a target β ⋆ ≥ 0, they hold simultaneously for every intermediate β ∈ [0, β ⋆ ], which becomes relevant in our applications of annealing in Section 5.

Annealing

Let f : X d → R and β ≥ 0 parameterize a Gibbs measure πβ ∝ exp(βf ). In this section, we develop a general framework for sampling from bounded-magnetization restrictions over X d , π≤k,β (x) ∝ exp (βf (x)) · Ix∈X d , ≤k

(34)

leveraging samplers for fixed-magnetization measures for 0 ≤ s ≤ k, πs,β (x) ∝ exp (βf (x)) · Ix∈Xsd ,

(35)

e.g., those constructed in Sections 3 and 4, as black boxes. Our approach decomposes the task into two stages: we first construct an approximate sampler for the cardinality s, and then, conditioned on this size, invoke the corresponding fixed-size sampler for πs,β to obtain a sample ∼ π≤k,β . 28

In Section 5.1, we employ an annealing-based technique (adapted from [Kol18]) to estimate the normalizing constants of the fixed-magnetization measures πs,β . In Section 5.2, we integrate these estimates into a unified framework for approximately sampling from bounded-magnetization measures π≤k,β . Finally, in Section 5.3, we instantiate our framework for bounded-magnetization variants of the Ising models in Sections 3 and 4, and derive the resulting mixing time bounds.

5.1

Estimating normalizing constants

We first recall an estimation procedure for normalization constants of a discretely-supported measure π ∈ P(Ω), based directly on [Kol18]. For some fixed f : Ω → R, where |Ω| < ∞, let X Z(β) := exp(βf (ω)) (36) ω∈Ω

to be the normalizing constant of the tempered Gibbs distribution πβ ∝ exp(βf ) at inverse temperature β ≥ 0. Observe that Z(0) = |Ω| is known exactly. For a target inverse temperature β ⋆ and b such that error tolerance ϵ ∈ (0, 1), our goal is to produce an estimate Z (1 − ϵ)Z(β ⋆ ) ≤ Zb ≤ (1 + ϵ)Z(β ⋆ ).

(37)

We now state a consequence of the estimation procedure of [Kol18]. Proposition 2. Let f : Ω → [−R, R] for R ≥ 0 and let β ⋆ ≥ 0, (δ, ϵ) ∈ (0, 21 )2 . There is an algorithm EstimateZ(f, β ⋆ , δ, ϵ, A), where A is an algorithm that samples from πβ ∈ P(Ω) ∝ exp (βf ) for any β ∈ [0, β ⋆ ]. The output of EstimateZ satisfies (37) with probability ≥ 1 − δ. Further, it uses 1 1 + β⋆R log (38) N =O 2 ϵ δ calls to A. ⋆

Proof. We first obtain (37) with probability > 12 using O( 1+Rβ ) calls to A, at which point the ϵ2 result follows by taking the median of O(log( 1δ )) copies via a standard Chernoff bound argument. If β ⋆ = 0 or R = 0, the claim is immediate by outputting |Ω|. Otherwise, define H(ω) :=

2R − f (ω) , R

t⋆ := Rβ ⋆ ,

P so that if we define ZH (t) := ω∈Ω exp(−tH(ω)) to be the corresponding normalizing constant for −H at inverse temperature t, then t ZH (t) = exp (−2t) Z , for all t ∈ [0, t⋆ ]. R Thus to produce the required estimate (37), it is enough to estimate Q :=

ZH (0) Z(0) = exp (2Rβ ⋆ ) ⋆ ZH (t ) Z(β ⋆ ) 29

to multiplicative error 1 ± ϵ. Also, observe that q := log Q satisfies q ∈ [t⋆ , 3t⋆ ], because the range of H(ω) is [1, 3]. Theorem 6 in [Kol18] with n = 3 and q = O(1 + Rβ ⋆ ), and its suggested parameters d = 64, m = O(1), and r = O(ϵ−2 ), now gives the claim. We note that Theorem 6 in [Kol18] is stated with an expected query complexity, but Markov’s inequality converts this into a deterministic runtime with a constant failure probability that can be folded into the estimator’s failure.

5.2

Approximate bounded-magnetization sampling

In this section, we develop a framework for approximate sampling from π≤k,β ⋆ , defined in (34), assuming access to samplers for the densities πs,β (35) for all 0 ≤ s ≤ k and 0 ≤ β ≤ β ⋆ . The key observation is that π≤k,β ⋆ admits the decomposition π≤k,β ⋆ =

k X

αi πi,β ⋆ ,

i=0

Zi (β ⋆ ) αi = Pk , ⋆ j=0 Zj (β )

Zi (β) :=

exp (βf (x)) .

(39)

x∈Xid

We first establish Lemma 14, which shows that accurate estimates of the mixture weights α and component distributions πi,β ⋆ suffice to guarantee an accurate approximation of π≤k,β ⋆ . Lemma 14. Let I be an index set and assume that π, π ′ ∈ P(Ω) admit the following decompositions: π = Ei∼ρ [µi ], Then,

π ′ = Ej∼ρ′ [µ′j ],

ρ, ρ′ ∈ P(I),

µi , µ′i ∈ P(Ω) for all i ∈ I.

TV π, π ′ ≤ TV ρ, ρ′ + sup TV µi , µ′i i∈I

Proof. Let πm := Ei∼ρ [µ′i ]. By the triangle inequality of TV distance, TV π, π ′ ≤ TV (π, πm ) + TV πm , π ′ . For the first term, convexity of | · | yields # " X 1X 1 TV (π, πm ) = Ei∼ρ [µi (ω) − µ′i (ω)] ≤ Ei∼ρ µi (ω) − µ′i (ω) ≤ sup TV µi , µ′i . 2 2 i∈I ω∈Ω

ω∈Ω

For the second term, couple i ∼ ρ, i′ ∼ ρ′ to minimize P[i ̸= i′ ], and then couple the draws from µ′i = µ′i′ whenever i = i′ . This produces different samples with probability ≤ TV (ρ, ρ′ ). We now present our sampler for bounded-magnetization measures (34) in Algorithm 2, and establish correctness in Lemma 15. An obstacle to directly applying Proposition 2 is that only approximate samplers π̂i,β are available in place of πi,β ; this is addressed via a coupling argument. d → [−R, R] for R ≥ 0 and let β ⋆ ≥ 0, δ ∈ (0, 1 ). For all 0 ≤ s ≤ k, Lemma 15. Let f : X≤k 2 0 ≤ β ≤ β ⋆ , assume that algorithm As (β, δ ′ ) returns a sample within δ ′ total variation distance of δ πs,β defined in (35). Algorithm 2 uses N calls, each to some As (β, 6N ) where ! k(1 + β ⋆ R) log( kδ ) N =O . (40) δ2

Further, its output x satisfies TV (Law(x), π≤k,β ⋆ ) ≤ δ. 30

δ )). Fix an optimal coupling Proof. For simplicity, in this proof we denote π̂s,β := Law(As (β, 6N between π̂s,β and πs,β for each oracle call to As . By assumption, each call produces an exact sample δ from πs,β with probability at least 1 − 6N . A union bound over the N calls implies that, with δ probability at least 1 − 6 , all oracle samples are exact. δ Next, note that our choice of N satisfies (38) with ϵ ← 6δ and δ ← 12k . Thus, all of the k + 1 ≤ 2k k b estimates {Zs }s=0 computed on Line 4 are correct with probability at least 1 − 6δ . Altogether, by Proposition 2, we have that with probability at least 1 − 3δ that each Zbs satisfies δ Ẑs δ ∈ 1 − ,1 + . Zs 6 6

Applying Lemma 4 and the estimate x ∈ [± 6δ ] =⇒ log(1 + x) ∈ [± 3δ ] yields δ TV (Law(s), Multinomial(Z)) ≤ , 3 where s is the sampled index on Line 6, and Zi := Zi (β ⋆ ) for all 0 ≤ i ≤ k as defined in (39). Combining this with Lemma 14 and the failure probability of 3δ on the final sample, we obtain TV (Law(x), π≤k,β ⋆ ) ≤

2δ 3

on the above high-probability event. Finally, the total failure probability from earlier was 3δ , and contributes additively in total variation. Therefore, TV (Law(x), π≤k,β ⋆ ) ≤ δ.

Algorithm 2: BMGibbsSampler(k, β ⋆ , f, {As }ks=0 , δ) d → [−R, R], 1 Input: Magnetization bound k ∈ [d], target inverse temperature β ⋆ ≥ 0, f : X≤k approximate samplers {As }ks=0 such that for all 0 ≤ β ≤ β ⋆ , δ ′ ∈ (0, 12 ), As (β, δ ′ ) returns a sample x with TV (Law(x), πs,β ) ≤ δ ′ , failure probability δ ∈ (0, 21 ) 2 Output: Sample x such that TV (Law(x), π≤k,β ⋆ ) ≤ δ where π≤k,β ⋆ ∝ exp(β ⋆ f ) · I·∈X d ≤k

3 for i = 0, 1, . . . , k do δ δ , 6δ , Ai (·, 6N )) where fi = f with domain Xid , and N is as defined Zbi ← EstimateZ(fi , β ⋆ , 12k in (40) 5 end b 6 s ∼ Multinomial(Z) 4

δ 7 return x ∼ As (β ⋆ , 3 )

5.3

Bounded-magnetization Ising models

In this section, we show how to apply Lemma 15 to bounded-magnetization Ising models, where 1 f (x) = x⊤ Jx + h⊤ x 2 31

d . We as in (34). To apply our framework, we require an upper bound on |f | over the domain X≤k derive such a bound, which is slightly tightened by the observation that shifting the potential by a constant does not affect the Gibbs measure. The same strategy can be applied in any setting where f has large magnitude but small variation over its domain.

Lemma 16. We have d for all x ∈ X≤k .

|f (x) − f (−1d )| ≤ (2k + 4d) ∥J∥2k,op + 2 ∥h∥k,1 ,

Proof. Let S := set(x) and v := 1S . Because x − (−1d ) = 2v, x + (−1d ) = 2v − 21d , we obtain |f (x) − f (−1d )| ≤ 2 v⊤ Jv + 2 (h − J1d )⊤ v ≤ 2k ∥J∥2k,op + 2 ∥h∥k,1 + 2 1⊤ d Jv , where the last step exploits the fact that v is k-sparse and ∥v∥2 ≤ follows from the bound X 1⊤ Jv ≤ 1⊤ Si Jv ≤ 2d ∥J∥2k,op , d

√

k. Finally, the conclusion

i∈[m]

for an arbitrary partition of [d] into k-sparse sets S1 ∪ S2 ∪ · · · ∪ Sm , where m ≤ kd + 1 ≤ 2d k . By instantiating Lemma 15 with the fixed-magnetization samplers from Sections 3 and 4, we derive samplers for the analogous bounded-magnetization Gibbs measures. For brevity, we only present the bounded-magnetization generalization of Theorem 4 here. d instead of Corollary 5. In the setting of Theorem 4, let π be as in (16) with restriction set X≤k Xkd . There is an algorithm that returns a sample x with TV (Law(x), π) ≤ δ, with probability ≥ 1−δ over the SK model (Model 1), in time ! dkρ log kδ kρ d O · k log + βk 2 ∥h∥∞ + log , where ρ := d + β ∥h∥k,1 . δ2 δ δ

Proof. The algorithm is Algorithm 2 with β ⋆ ← β, where we use

T lazy T DU,πs,β π0

δ as As (β, 6N ) for all 0 ≤ s ≤ k, for an arbitrary π0 ∈ P(Xsd ). Theorem 4 states that we need N d 2 T = O k log + βk ∥h∥∞ + log , δ δ

merging terms for simplicity. Under the success of Item 1 in Lemma 12, the assumed range on k in Theorem 4, and the condition (31), we have β∥J∥2k,op = O(1). Thus, Lemma 16 gives ⋆

β R = O(ρ) =⇒ N = O

kρ log δ2

k δ

in our application of Lemma 15. The conclusion follows from the implementation in Remark 2. 32

Bayesian Sparse Linear Regression

In this section, we apply the local-to-global framework developed earlier to design an approximate sampler for the spike-and-slab posterior in Model 3. As in [KSTZ25], we first use a sparse recovery preprocessing step to remove coordinates whose inclusion is determined from the observations, up to negligible posterior mass. We then sample directly from the preprocessed support posterior. In Section 6.1, we recall the preprocessing reduction of [KSTZ25]. In Section 6.2, we prove rapid mixing of every fixed-size restriction of the exact support posterior and lift these samplers to the bounded-size posterior using tools from Section 5. Finally, Section 6.3 combines the pieces.

6.1

Setup

Throughout this section, we follow the shorthand E := X⊤ X − Id ,

L := log

d δ

to simplify the statement of bounds, where δ is a specified failure probability. We begin with a technical lemma used to decompose the output of the [KSTZ25] reduction. Lemma 17. Let M ∈ Sd×d and let v ∈ Rd satisfy |supp(v)| ≤ r. Then following the notation (7), ∥Mv∥s,2 ≤ ∥M∥r+s,op ∥v∥2 for all s ∈ [d − r]. Proof. Fix T ⊆ [d] with |T | ≤ s, and let S := T ∪ supp(v) so |S| ≤ r + s. Then, ∥[Mv]T ∥2 ≤ ∥MS×S ∥op ∥v∥2 ≤ ∥M∥r+s,op ∥v∥2 .

We next state a variant of the preprocessing strategy used by [KSTZ25]. To obtain our improved sample complexity, we require a somewhat more fine-grained guarantee on its output than the ℓ∞ error bound used in prior works [KSTZ25, CLTZ26]. We state our required property in the form of a decomposition (42), which splits the preprocessing output vector z into two terms: a short vector (with bounded sparse ℓ2 norm), and a flat vector (with bounded ℓ∞ norm). Notably, this strategy is reminiscent of a similar decomposition used algorithmically by [KLL+ 23]. Proposition 3 (Section 3.1, Lemma 8, [KSTZ25] and Theorem 3, [CLTZ26]). In the setting of Model 3, let δ ∈ (0, 1). Then with probability ≥ 1 − δ over the randomness in Model 3, if k := 24(k̄ + log 1δ ), X ∼i.i.d. N (0, n1 ) and n = Ω(k log dδ ) for a large enough constant, there is a subset d U ⊆ [d] and a vector z ∈ Rd that can be computed in time O(nd log δ min{1,σ} ), satisfying |U c | = O(k),

TV (π̂supp , πsupp ) ≤ δ,

where π̂supp is supported on S ∪ U c where S ⊆ U , U c := [d] \ U, with ! Y qi 1 1 c exp ∥zS∪U c ∥2A−1 c √ π̂supp (S ∪ U ) ∝ · I(|S| ≤ k). S∪U 1 − qi 2 c det A S∪U i∈S 33

(41)

Moreover, z admits a decomposition z = z(0) +z(1) , such that for any constant a, there exist constants C0 , C1 > 0 (where C1 depends only on a) with 1 √ C1 kL (0) ≤ √ . ≤ C0 σ + L, z(1) z (42) σ σ n ak,2 ∞ Proof. We explain how to derive this result from [KSTZ25, CLTZ26], as it is not stated in this form. b where θb is the estimator used by Lemma 6 of [KSTZ25] satisfying First, U c is set to supp(θ) √ (43) =O σ L , θb − θ ⋆ ∞

√

with probability ≥ 1 − 8δ . Also, ∥θ ⋆ ∥∞ = O( L) under Model 3 with probability ≥ 1 − 8δ , so √ (44) = O (1 + σ) L . θb ∞

The existence of such an estimator θb that takes inputs (X, y), runs within the stated runtime, and b = O(k) follows from Theorem 3, [CLTZ26] for Gaussian ensembles. We satisfies (43) and |supp(θ)| note that the bound in (43) is obtained by using the tighter error bound for Gaussian observation matrices, discussed at the end of Page 15, [CLTZ26]. The closeness of πsupp and π̂supp then follows from Lemma 6 in [KSTZ25], where the form of π̂supp comes from Fact 3. Next, following Eq. (23) in [KSTZ25], we let 1 z := 2 X⊤ Xθ ⋆ + ξ − Xθb − θb σ 1 ⊤ 1 = 2 X ξ + θ ⋆ − θb − θb + 2 E θ ⋆ − θb . |σ {z } |σ {z } :=z(0)

:=z(1)

b because y = Xθ ⋆ + ξ is given as an input. We Clearly z can be computed given knowledge of θ, now verify the conditions (42). For z(0) , the √ stated bound in (42) follows by combining Lemma 1 of [KSTZ25], which gives ∥X⊤ ξ∥∞ = O(σ L) except with probability 4δ , with (43) and (44). For z(1) , we first condition on nnz(θ ⋆ ) ≤ k, which occurs with probability ≥ 1 − 4δ by Corollary 1, [KSTZ25]. Thus, θ ⋆ − θb is bk-sparse for a constant b. The bound in (42) then follows from b and s ← ak, where Lemma 17 with v ← θ ⋆ − θ, r ! √ kL . θ ⋆ − θb = O(σ kL), ∥E∥(a+b)k,op = O n 2 The first inequality above follows from (43) and our sparsity bound. For the second, following the proof of Corollary 4 up until (30), and adjusting the failure probability to union bound over subsets as in that proof, shows that for all s ∈ [d], with probability ≥ 1 − 4δ , ! r sL sL ∥E∥s,op = O + . (45) n n Finally, the claim follows from a union bound over all five random events in the proof. 34

Remark 4. For notational simplicity, in Section 6.2 we present the argument assuming C := U c = ∅ in Proposition 3. The general case of C is identical up to a universal constant-factor enlargement of the sparsity k. Indeed, write the full-support potential in (41) as X qi 1 1 F (T ) := log + z⊤ A−1 zT − log det AT for all C ⊆ T ⊆ [d]. 1 − qi 2 T T 2 i∈T

Then the posterior on the collapsed support ∅ ⊆ S ⊆ U is ∝ exp(FC (S))I(|S| ≤ k), where FC (S) := F (C ∪ S). All conclusions in Section 6.2 go through unchanged after lifting by C and adjusting constant factors to account for |C| = O(k). For example, Lemma 19 applies verbatim after replacing R ← C ∪ R, and since we only considered |R| ≤ k, the same sparsity bounds hold up to constant factors after this adjustment. Similarly, Lemma 20 considers differences between potentials F (S), which goes through unchanged after replacing F ← FC and S ← C ∪ S. S In the rest of the section, we fix the universe U returned by Proposition 3. We write U≤k := ks=0 Us , analogously to (13). Conditioned on the success of Proposition 3, it is enough to sample from the density π̂supp in (41). We write the (unnormalized) density as exp(F (S)), where we define F (S) =

log

i∈S

qi 1 1 −1 log det AS . + z⊤ S AS zS − 1 − qi 2 2

We now manipulate this expression to be on a more convenient scale. In particular, let γ :=

σ2 , 1 + σ2

τ :=

1 , 1 + σ2

s :=

Under this scaling, and defining s(0) := s(0)

∞

√

LS := γAS = IS + τ ES for all S ⊆ [d].

γτ z,

γτ z(0) and s(1) :=

√ ≤ C0 L,

s(1)

√

(46)

γτ z(1) , (42) implies

C1 kL ≤ √ . n ak,2

(47)

Moreover, again using the notation (46) and combining all linear terms, we obtain F (S) = h⊤ 1S +

1 1 ⊤ −1 qi 1 sS LS sS − log det LS , where hi := log + log γ. 2τ 2 1 − qi 2

(48)

Finally, for β ∈ [0, 1] and 0 ≤ s ≤ k, we define ν≤s,β (S) ∝ exp(βF (S))I(S ∈ U≤s ),

νs,β (S) ∝ exp(βF (S))I(S ∈ Us ).

(49)

With this notation, our target support posterior is ν = ν≤k,1 .

6.2

Trickle down for support posterior

In this section, we apply the trickle down framework of Section 4 to sampling under Model 3. We begin by collecting several additional estimates we require of our draw from the model. Lemma 18. Under Model 3, assume n = Ω(kL) for an appropriate constant. Under the success of the event in Proposition 3, we have the following additional guarantee, for a universal constant p C > 0, and α := L/n. Simultaneously for every R ∈ U≤k , B := U \ R, and u, v ∈ Rd , √ √ ∥E∥2k,op ≤ Cα k, max |Eij | ≤ Cα, ∥ER×B uB ∥2 ≤ Cα k ∥u∥2 + ∥u∥1 , (i,j)∈[d]×[d]

and

u⊤ Ev ≤ C α (∥u∥1 ∥v∥2 + ∥u∥2 ∥v∥1 ) + α2 ∥u∥1 ∥v∥1 .

(50)

Proof. The only event we condition on in this proof is that (45) holds. The first two inequalities then follow by plugging s = 2k and s = 2 into (45), and using our lower bound on n. For inequality (50), assume u, v ̸= 0d , else the claim is immediate. Then let q ∥v∥1 1 2 t := =⇒ t ∥u∥1 + ∥v∥1 = 2 ∥u∥1 ∥v∥1 . ∥u∥1 t Since 1 u Ev = 4 ⊤

1 ⊤ 1 ⊤ 1 1 1 tu + v tu − v E tu + v − E tu − v , t t 4 t t

the claim follows from Lemma 11 applied to the two terms, the triangle inequality, and 1 1 t ∥u∥1 + ∥v∥1 t ∥u∥2 + ∥v∥2 = 2 (∥u∥1 ∥v∥2 + ∥u∥2 ∥v∥1 ) , t t 2 1 t ∥u∥1 + ∥v∥1 = 4 ∥u∥1 ∥v∥1 . t ⊤ For the third inequality, applying (50) and √ the fact that ∥ER×B u∥2 = v EuB for some √ unit vector v supported only on R, such that ∥v∥1 ≤ k, yields the claim upon simplifying with α k = O(1).

Conditioned on the success of the event in Lemma 18, we next show how to apply our framework in Section 4 to sample from the posterior density ν≤k,1 (49). We begin with an exact characterization of the potential gain by including a set of coordinates. Lemma 19. For a fixed R ⊆ U, let B := U \ R and define RR := EB×B − τ EB×R L−1 R ER×B ,

rR := sB − τ EB×R L−1 R sR .

Then, for every T ⊆ B, letting RT := [RR ]T ×T and rT := [rR ]T , F (R ∪ T ) − F (R) = h⊤ 1T +

1 ⊤ 1 rT (IT + τ RT )−1 rT − log det (IT + τ RT ) . 2τ 2

Proof. With the coordinates ordered as R, T , LR τ ER×T LR∪T = . τ ET ×R LT Its Schur complement with respect to the R block is LT − τ 2 ET ×R L−1 R ER×T = IT + τ RT . Consequently, det LR∪T = det(LR ) det (IT + τ RT ). The block inverse formula also gives ⊤ −1 −1 ⊤ −1 s⊤ (IT + τ RT )−1 sT − τ ET ×R L−1 R∪T LR∪T sR∪T = sR LR sR + sT − τ ET ×R LR sR R sR −1 −1 ⊤ = s⊤ rT . R LR sR + rT (IT + τ RT )

Substituting both of the above displays into (48) now proves the claim. 36

(51)

We now give the main technical result in this section, which provides various estimates required to apply the parameter bounds from Lemma 9 to an appropriate interaction matrix. Lemma 20. Assume the events of Proposition 3 and Lemma 18 hold. Let R ∈ U≤k and B := U \ R, and following notation of Lemma 19, let K ∈ RB×B have zero diagonal and satisfy Kij := F (R ∪ {i, j}) − F (R ∪ {i}) − F (R ∪ {j}) + F (R) for all (i, j) ∈ B × B, i ̸= j. p Then the following bounds hold for a universal constant C ′ > 0, letting α := L/n. 1. For all u ∈ RB , following the notation (51), Rij ui uj ≤ C ′ α2 k ∥u∥22 + α ∥u∥1 ∥u∥2 + α2 ∥u∥21 .

X (i,j)∈B×B i̸=j

2. For all (i, j) ∈ B × B, i ̸= j, Kij = −Rij ri rj + ϵij , for |ϵij | ≤ C ′ α2 + α4 k 2

1 + r2i + r2j .

3. For all u ∈ RB , ∥r ◦ u∥1 ≤ C

′

√

kL L ∥u∥1 + √ ∥u∥2 , n

and X

|ui |r2i ≤ C ′

i∈B

∥r ◦ u∥2 ≤ C

′

√

kL L+ √ n

∥u∥2 ,

k 2 L2 L ∥u∥1 + ∥u∥2 . n

Proof. We proceed with the three claims in order. First, observe that by Lemma 18 and |R| ≤ k, √ −1 L−1 ≤ 1 − Cτ α k ≤ 2, R op

(52)

using τ ≤ 1 and taking n large enough. Also, again by applying Lemma 18, ∥R∥max :=

max (i,j)∈B×B

|Rij | ≤

max (i,j)∈B×B

|Eij | + L−1 R op ER×{i} 2 ER×{j} 2

2 √ 1 ≤ Cα + 2 Cα( k + 1) ≤ Cα + 8C 2 α2 k ≤ . 2 The last inequality again used our lower bound on n. Hence,

Rij ui uj ≤ u⊤ Ru + ∥R∥max ∥u∥22

(i,j)∈B×B i̸=j

(53)

≤ u⊤ Ru + 8C 2 α2 k ∥u∥22 + Cα ∥u∥1 ∥u∥2 , and by again applying Lemma 18, we have u⊤ Ru ≤ u⊤ Eu + 2 ∥ER×B u∥22 ≤ 2C α ∥u∥1 ∥u∥2 + α2 ∥u∥21 + 4C 2 α2 k ∥u∥22 + α2 ∥u∥21 . Combining the above two displays proves Item 1. Next, for all i ∈ B, let di := 1 + τ Rii , so that di ∈ [1 ± τ ∥R∥max ] ⊆ [ 21 , 32 ]. Also, let Dij := det I{i,j} + τ R{i,j} = di dj − τ 2 R2ij so that I{i,j} + τ R{i,j}

−1

di τ Rij

τ Rij dj

−1

1 = Dij

dj −τ Rij

−τ Rij . di

Then, applying Lemma 19 with T = ∅, {i}, {j}, {i, j}, the linear terms cancel, so ! 2 2 r 1 1 r 1 1 j i Kij = dj r2i + di r2j − 2τ Rij ri rj − + − log (Dij ) + log (di dj ) 2τ Dij 2τ di dj 2 2 r2 r2 ! τ R2ij dii + djj − 2Rij ri rj τ 2 R2ij 1 = − log 1 − 2Dij 2 di dj ! ! τ R2ij r2i r2j τ 2 R2ij Rij ri rj 1 =− + + − log 1 − . Dij 2Dij di dj 2 di dj

(54)

Using the estimates |di − 1|, |dj − 1|, |Dij − 1| = O(∥R∥max ) = O(α + α2 k) and Taylor expanding each error term (since (53) bounds ∥R∥max by an arbitrary constant), gives for large enough C ′ , 2 2 C′ 1 α + α2 k ri + r2j , −1 ≤ Dij 6 ! 2 2 2 τ Rij ri rj 2 2 C′ + ≤ α + α2 k ri + r2j , 2Dij di dj 6 ! τ 2 R2ij 2 C′ log 1 − ≤ α + α2 k . di dj 6

|Rij ri rj | ·

Combining these bounds within (54), and using (a + b)2 ≤ 2(a2 + b2 ), then yields Item 2. (1)

To conclude, let c := EB×R L−1 and s = s(0) +s(1) , R sR and d := sB −τ c, so that following Lemma 19√ (0) (0) we have rR = sB + d. Recall from (47) that sB has all entries bounded by O( L), such that √ ∥s(0) ◦ u∥1 = O( L) ∥u∥1 ,

√ ∥s(0) ◦ u∥2 = O( L) ∥u∥2 ,

X i∈B

(0)

|ui |(si )2 = O (L) ∥u∥1 .

By choosing C ′ large enough in Item 3, all of the above contributions fit asymptotically within the claimed budgets. It thus suffices to prove Item 3 using d to reweight u rather than r. Next, we bound ∥d∥k,2 . Let T ⊆ B index the largest k coordinates of d by magnitude. By (47) and n = Ω(kL), for every R ∈ U≤k , √ √ = O( kL). + s(1) ∥sR ∥2 ≤ k s(0) ∞

k,2

Then, (52) and the third bound in Lemma 18 imply, for a large enough constant C ′′ , r √ √ C ′′ ∥cT ∥2 ≤ 2Cα k · O( kL) ≤ kL . n q ′′ Combining with (47) and τ ≤ 1 gives ∥dT ∥2 ≤ kL Cn after adjusting C ′′ . Thus, letting the ith largest magnitude amongst d’s coordinates be denoted d(i) , we have shown d2(i) ≤

1 C ′′ k 2 L2 · , for all i ∈ |B|. n min(i, k)

Now, similarly denoting the ith largest magnitude amongst u’s coordinates by u(i) ,   r ′′ 2 2 X u u C k L  (i) (i) √ +√  ∥d ◦ u∥1 ≤ · n i k i∈[|B|] v  u r u X ′′ 2 2 √ √ C k L u C ′ kL1.5 1 ·t ∥u∥22 ≤ C ′ L ∥u∥1 + √ ≤ C ′ L ∥u∥1 + ∥u∥2 , n i n

(55)

(56)

i∈[|B|]

where the second inequality was by Cauchy-Schwarz. Similarly,   2 2 ′′ k 2 L2 X u u C ′′ k 2 L2 C (i)  (i) · + ≤ ∥u∥22 . ∥d ◦ u∥22 ≤ n i k n i∈[|B|]

Finally, recalling n = Ω(kL), and applying the Cauchy-Schwarz inequality as in (56),   ′′ 2 2 X X X u u C ′ k 2 L2 C k L  (i) (i)  |ui |d2i ≤ u(i) d2(i) ≤ · + ≤ ∥u∥2 + C ′ L ∥u∥1 . n i k n i∈B

i∈[|B|]

We are finally ready to give our application of the trickle down framework. Proposition 4. Assume the events of Proposition 3 and Lemma 18 hold, and that n = Ω k 1.5 L2 + kL3 for a sufficiently large constant. Then for all 2 ≤ s ≤ k and β ∈ [0, 1], following the notation (49), 1 -Poincaré inequality. T DU,νs,β satisfies a 2k 39

Proof. Let R ∈ Us−2 , and define K ∈ RB×B following the notation in Lemma 20. We first claim u⊤ Ku ≤

9 1 ∥u∥22 + ∥u∥21 , for all u ∈ RB . 16 8k

(57)

To see this, decompose K as in Item 2 of Lemma 20, and by homogeneity, assume that ∥u∥1 = 1 and ∥u∥2 = q. By combining the estimates in Items 1 and 3, and using that X (i,j)∈B×B i̸=j

k 2 L2 2 Rij ri rj ui uj = O α k L + q2 n √ √ kqL kL +O α L 1+ √ q L+ √ n n k 2 L2 q 2 2 +O α L 1+ n 2 3 k L3 kL2.5 k 2 L3 k 2 L4 kL =O + 2 + + 1.5 + 2 q2 n n n n n 1.5 2 L kL2 L √ + +O q +O n n n 1 1 1 1 1 ≤ q2 + , ≤ q2 + √ q + 4 32k 2 16k 16 k

where the hidden constants above depend only on C ′ . The last line then follows by taking n large enough as stated. Next, for the error term in Item 2, recalling ∥u∥1 = 1, X

ϵij ui uj = O α2 + α4 k 2 ·

(i,j)∈B×B i̸=j

|ui uj | 1 + r2i + r2j

(i,j)∈B×B i̸=j

! 2

4 2

=O α +α k

1+2

|ui |r2i

i∈B

k 2 L2 q 2 4 2 L+ =O α +α k n 2 3 2 4 4 k L k L L k 2 L3 + 2 =O + 3 q +O n2 n n n 1 1 1 2 1 ≤ √ q+ ≤ q + , 32k 16 16k 16 k

(58)

where we used Item 3 in the third line, and again simplified by taking n large enough. Combining the above two displays gives the desired (57). Next, we observe that the link graph induced by R, in the sense of Lemma 6, has edge weight matrix W exactly given in the form required by Lemma 7, with K defined in Lemma 20 and 1 for all i ∈ B. ai := exp β F (R ∪ {i}) − F (R) 2 40

Indeed, for (i, j) ∈ B × B with i ̸= j, ai aj exp (βKij ) = exp (βF (R ∪ {i, j})) as is required by the link graph. We now provide bounds on m(βK) and ρ(βK) as used in Lemma 7. (0)

To bound m(βK), recall the proof of Lemma 20 (i.e., rR = sB + d, (47), and (55)) that we √ from kL showed ∥rR ∥∞ = O( L + √ ). Thus, combining with (53) and Item 2 shows that n m(βK) ≤ m(K) = O

k 2 L2 k 2 L2 1 2 4 2 α+α k L+ +O α +α k L+ ≤ . n n 16 2

To bound ρ(βK), it suffices to restrict to ∥u∥1 = 1, ∥u∥2 = q by homogeneity. We follow the proof of Lemma 9, and write X = βK + Q where Xij = exp(βKij ) − 1 entrywise. By Taylor expansion, 2 2 we have that |Qij | ≤ β 2 K2ij . Again, using Item 2, denoting A := α2 + α4 k 2 = Ln + knL2 ≤ 1, K2ij = O Ar2i r2j + A2 1 + r4i + r4j . Thus, applying the bounds from Item 3, and simplifying, !! 4 L4 q 2 X k + A2 1 + ∥rR ∥2∞ 1+ |ui |r2i u⊤ Qu = O A L2 + n2 i∈B 2 2 4 4 2 k L k 2 L2 q k L q 2 2 +A L+ L+ =O A L + n2 n n 2 2 3 4 4 2 4 4 A k L A k L Ak L 2 q +O + q =O 2 n n n2 A2 k 2 L3 + O AL2 + A2 L2 + n 1 1 1 1 1 ≤ q2 + √ q + ≤ q2 + . 32 32k 16 16k 16 k Above, the first line used a similar simplification as in (58), as well as our bound on ∥rR ∥∞ . The last inequality used our lower bound on n to simplify the various terms. In conclusion, combining with (57), we have shown that 1 1 u⊤ (X − IB ) u ≤ ∥u∥21 − ∥u∥22 =⇒ ρ(βK) ≤ . 4k 4k The rest of the proof follows analogously to Lemma 8, using our bounds on ρ(βK), m(βK). We conclude the section by bounding the range of the potential, for use with Section 5. Lemma 21. Assume the events of Proposition 3 and Lemma 18 hold, and that n = Ω k 1.5 + kL3 for a sufficient constant. Then following the notation (48), max |F (S)| = O k 1 + σ 2 L + k ∥h∥∞ . S∈U≤k

2 Proof. For every S ∈ U≤k , Lemma 18 gives ∥ES×S ∥op ≤ 12 , so the spectrum of L−1 S is in [ 3 , 2]. √ Moreover, (47) and n = Ω(kL) give ∥sS ∥2 = O( kL). Therefore,

1 ⊤ −1 1 sS LS sS ≤ ∥sS ∥22 = O k(1 + σ 2 )L . 2τ τ Further, |h⊤ 1S | ≤ k ∥h∥∞ . The claim follows, as |log det LS | = O(k) using our spectrum bound. For completeness, we record that under the event of Lemma 18, combining the mixing guarantee in Lemma 1, the spectral gap in Proposition 4, and the Lemma 21, shows that it suffices to take 1 2 2 T =Ω k 1 + σ L + ∥h∥∞ log (59) δ′ steps of lazy(T DU,νs,β ), for a sufficient constant, to sample from within δ ′ TV from νs,β (49).

6.3

Main result

In this section, we finally put together the pieces to give our main sampling result for spike-and-slab posterior densities under Model 3. We use the notation FC (S) := F (C ∪ S) defined in Remark 4, C (S) ∝ exp(βF (S))I(S ∈ U ). and for notational simplicity, we also let νs,β s C Algorithm 3: SASPosteriorSampler(X, y, σ, q, δ) 1 1 Input: X ∈ Rn×d , y ∈ Rn , σ > 0, q ∈ (0, 1)d , and δ ∈ (0, 2 ) 2 Output: Sample θ such that TV (Law(θ), π(· | X, y)) ≤ δ with probability ≥ 1 − δ under

Model 3 δ 3 (z, U) ← output of Proposition 3, with error parameter 3 3 4 k ← 24(k̄ + log( δ )) 5 C ← Uc 6 A0 ← algorithm that always outputs ∅ 7 for s ∈ [k] do C 8 As (β, δ ′ ) ← (lazy(T DU,νs,β ))T applied to an arbitrary start, for T as in (59) 9 end δ 10 S ∼ BMGibbsSampler(k, 1, FC , {As }k s=0 , 3 ) e ← S ∪ Uc 11 S −1 −1 12 return θ ∼ N (A e b e, A e ) following Fact 3 S S

Theorem 5. Let δ ∈ (0, 21 ), and following the notation of Model 3, let k = 24(k̄ + log( 3δ )). Suppose 1.5 2 d 3 d n = Ω k log + k log δ δ for a sufficiently large constant, and assume that q ∈ [η, 1 − η]d for η > 0. Then, with probability ≥ 1 − δ over the randomness of Model 3, the output of Algorithm 3 satisfies TV (Law(θ), π(· | X, y)) ≤ δ. 42

1 ), the algorithm runs in time Defining Q := k(1 + σ 2 ) log( dδ ) + k log( η1 + ησ

ndk 3 Q2 log δ2

k kQ log . δ δ

Proof. Throughout the proof, condition on the success of Proposition 3 and the event in Lemma 18, which give the failure probability over Model 3. Under these events, the extended draw Se in Algorithm 3 is within total variation distance δ from the support posterior π(supp(·) | X, y). This is because the density π̂supp in (41) is within TV 3δ of the support posterior by Proposition 3 and Remark 4, and the guarantees of BMGibbsSampler in Lemma 15 imply S is within TV 3δ of an exact sample from π̂supp . Finally, the sample θ | Se is exact (Fact 3) and cannot increase TV. It remains to bound the implementation cost. Observe that Lemma 21 shows that the potential FC is bounded by O(Q) under the assumption q ∈ [η, 1 − η]d , so Lemma 15 uses ! kQ log( kδ ) N =O δ2 δ calls to algorithms As (β, 6N ), each using T steps of the down-up walk where T is defined in (59). We claim that each step can be implemented in O(ndk) time, which gives the runtime claim. U To prove the implementation cost of As for s ∈ [k], fix some R ∈ s−1 that is the result of dropping an element (i.e., after Line 3 of Algorithm 1 has completed). For simplicity, we consider the case when C = ∅. In the general case every core R appearing below is replaced by C ∪ R, whose size remains O(k). To implement Line 4, Lemma 19 gives

FC (R ∪ {j}) − FC (R) = hj +

1 r2j 1 − log dj , 2τ dj 2

where we followed the notation of (54), so rj = sj − τ E{j}×R L−1 R sR ,

dj := 1 + τ Rjj = 1 + τ Ejj − τ 2 E{j}×R L−1 R ER×{j} .

We can compute EB×R = [X⊤ X − Id ]B×R in time O(ndk), and similarly we can compute and invert LR = IR + τ ER in this time. We can also compute all diagonal entries of E in time O(nd). Thus, computing all of the rj takes time O(ndk), and computing each dj takes time O(k 2 ) = O(nk) to compute the relevant quadratic form, so computing all of them takes time O(ndk) as well.

Acknowledgments We thank Thuy-Duong (June) Vuong for her participation at an earlier stage of this project, as well as Sidhanth Mohanty for several helpful conversations. We also thank an anonymous FOCS reviewer for making a suggestion that led to our strategy in Section 4. SK gratefully acknowledges funding support from the Amazon AI PhD Fellowship. KT and YZ thank the NSF AI Institute for Foundations of Machine Learning (IFML) for supporting this project. 43

AI Disclosure A preliminary version of this paper, consisting of Sections 3, 5, and 6 (at a measurement complexity n ≳ k 2 ), was previously submitted to FOCS, with all ideas contributed by the authors. Based on a reviewer’s suggestion, the authors used GPT 5.6 Pro to explore applications of the trickle down theorem to sharpen our results. Specifically, the key perturbation strategy in Lemma 7 was suggested by GPT, which led to a weaker variant of Theorem 1 (Corollary 3). The authors then built on this approach in our final applications in Theorems 1 and 2. We also acknowledge the use of LLMs in understanding the literature on the Almeida-Thouless line, as well as the prior works [CDKP22, KPPY25], which helped us prepare Appendices A and B. The manuscript was written solely by the authors, who take full responsibility for the organization and presentation of all results.

References [ADCS14]

Michael Aizenman, Hugo Duminil-Copin, and Vladas Sidoravicius. Random currents and continuity of ising model’s spontaneous magnetization. Communications in Mathematical Physics, 334(2):719–742, July 2014.

[ADV+ 25]

Josh Alman, Ran Duan, Virginia Vassilevska Williams, Yinzhan Xu, Zixuan Xu, and Renfei Zhou. More asymmetry yields faster matrix multiplication. In Proceedings of the 2025 Annual ACM-SIAM Symposium on Discrete Algorithms, SODA 2025, pages 2005–2039. SIAM, 2025.

[AGS85]

Daniel J Amit, Hanoch Gutfreund, and Haim Sompolinsky. Storing infinite numbers of patterns in a spin-glass model of neural networks. Physical review letters, 55(14):1530, 1985.

[AGS87]

Daniel J Amit, Hanoch Gutfreund, and Haim Sompolinsky. Statistical mechanics of neural networks near saturation. Annals of physics, 173(1):30–67, 1987.

[AGZ10]

Greg W Anderson, Alice Guionnet, and Ofer Zeitouni. An introduction to random matrices. Cambridge university press, 2010.

[AJK+ 22]

Nima Anari, Vishesh Jain, Frederic Koehler, Huy Tuan Pham, and Thuy-Duong Vuong. Entropic independence: optimal mixing of down-up random walks. In STOC ’22: 54th Annual ACM SIGACT Symposium on Theory of Computing, pages 1418– 1430. ACM, 2022.

[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, STOC 2024, pages 1094–1105. ACM, 2024.

[AL20]

Vedat Levi Alev and Lap Chi Lau. Improved analysis of higher order random walks and applications. In Proceedings of the 52nd annual ACM SIGACT symposium on theory of computing, pages 1198–1211, 2020.

[ALGV19]

Nima Anari, Kuikui Liu, Shayan Oveis Gharan, and Cynthia Vinzant. Log-concave polynomials II: high-dimensional walks and an FPRAS for counting bases of a matroid.

In Moses Charikar and Edith Cohen, editors, Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing, STOC 2019, pages 1–12. ACM, 2019. [AMS22]

Ahmed El Alaoui, Andrea Montanari, and Mark Sellke. Sampling from the sherrington-kirkpatrick gibbs measure via algorithmic stochastic localization. In 63rd IEEE Annual Symposium on Foundations of Computer Science, FOCS 2022, pages 323–334. IEEE, 2022.

[BAR26]

Afonso S. Bandeira, Ahmed El Alaoui, and Almut Rödder. Mixing of glauber dynamics on high overlap gibbs measures, 2026.

[BBD24]

Roland Bauerschmidt, Thierry Bodineau, and Benoit Dagallier. Kawasaki dynamics beyond the uniqueness threshold. Probability Theory and Related Fields, 192(1–2):267– 302, 2024.

[BG08]

Adriano Barra and Francesco Guerra. About the ergodic regime in the analogical hopfield neural networks: moments of the partition function. Journal of mathematical physics, 49(12), 2008.

[BH24]

Joan Bruna and Jiequn Han. Provable posterior sampling with denoising oracles via tilted transport. In Advances in Neural Information Processing Systems 38: Annual Conference on Neural Information Processing Systems 2024, NeurIPS 2024, 2024.

[Bon14]

Claudio Bonati. The peierls argument for higher dimensional ising models. European Journal of Physics, 35(3):035002, March 2014.

[BRG21]

Ray Bai, Veronika Rockova, and Edward I. George. Spike-and-slab meets lasso: A review of the spike-and-slab lasso. Handbook of Bayesian Variable Selection, pages 81–108, 2021.

[BvEN99]

Anton Bovier, Aernout CD van Enter, and Beat Niederhauser. Stochastic symmetrybreaking in a gaussian hopfield model. Journal of statistical physics, 95(1):181–213, 1999.

[BvH16]

Afonso S. Bandeira and Ramon van Handel. Sharp nonasymptotic bounds on the norm of random matrices with independent entries. The Annals of Probability, 44(4):2479– 2506, 2016.

[BY22]

Christian Brennecke and Horng-Tzer Yau. The replica symmetric formula for the sk model revisited. Journal of Mathematical Physics, 63(7), 2022.

[CDKP22]

Charlie Carlson, Ewan Davies, Alexandra Kolla, and Will Perkins. Computational thresholds for the fixed-magnetization ising model. In Proceedings of the 54th Annual ACM SIGACT Symposium on Theory of Computing, STOC 2022, page 1459–1472. Association for Computing Machinery, 2022.

[CGM19]

Mary Cryan, Heng Guo, and Giorgos Mousa. Modified log-sobolev inequalities for strongly log-concave distributions. In 60th IEEE Annual Symposium on Foundations of Computer Science, FOCS 2019, pages 1358–1370. IEEE Computer Society, 2019.

[Chi96]

H. Chipman. Bayesian variable selection with related predictors. The Canadian Journal of Statistics, 24:17–36, 1996. 45

[CLTZ26]

Ziyun Chen, Jerry Li, Kevin Tian, and Yusong Zhu. Separating oblivious and adaptive models of variable selection. arXiv preprint arXiv:2602.16568, 2026.

[CPS09]

Carlos M. Carvalho, Nicholas G. Polson, and James G. Scott. Handling sparsity via the horseshoe. In Proceedings of the Twelfth International Conference on Artificial Intelligence and Statistics, AISTATS 2009, volume 5 of JMLR Proceedings, pages 73–80. JMLR.org, 2009.

[CRT06]

Emmanuel J Candès, Justin Romberg, and Terence Tao. Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information. IEEE Transactions on information theory, 52(2):489–509, 2006.

[CSHVdV15] Ismael Castillo, Johannes Schmidt-Hieber, and Aad Van der Vaart. Bayesian linear regression with sparse priors. The Annals of Statistics, 43(5):1986–2018, 2015. [CT05]

Emmanuel J Candes and Terence Tao. Decoding by linear programming. IEEE transactions on information theory, 51(12):4203–4215, 2005.

[CT06]

Emmanuel J Candes and Terence Tao. Near-optimal signal recovery from random projections: Universal encoding strategies? IEEE transactions on information theory, 52(12):5406–5425, 2006.

[CvdV12]

Ismaël Castillo and Aad van der Vaart. Needles and straw in a haystack: Posterior concentration for possibly sparse sequences. The Annals of Statistics, 40(4):2069–2101, August 2012.

[dAT78]

Jairo RL de Almeida and David J Thouless. Stability of the sherrington-kirkpatrick solution of a spin glass model. Journal of Physics A: Mathematical and General, 11(5):983–990, 1978.

[DLSS26]

Ewan Davies, Holden Lee, Juspreet Singh Sandhu, and Jonathan Shi. Potential hessian ascent III: sampling the sherrington-kirkpatrick model at beta < 1/2. CoRR, abs/2605.03718, 2026.

[Dob68]

PL Dobruschin. The description of a random field by means of conditional probabilities and conditions of its regularity. Theory of Probability & Its Applications, 13(2):197– 224, 1968.

[Don06]

David L Donoho. Compressed sensing. IEEE Transactions on information theory, 52(4):1289–1306, 2006.

[DP23]

Ewan Davies and Will Perkins. Approximately counting independent sets of a given size in bounded-degree graphs. SIAM Journal on Computing, 52(2):618–640, 2023.

[EKZ22]

Ronen Eldan, Frederic Koehler, and Ofer Zeitouni. A spectral condition for spectral gap: fast mixing in high-temperature ising models. Probability theory and related fields, 182(3):1035–1051, 2022.

[Ell12]

Richard S Ellis. Entropy, large deviations, and statistical mechanics. Springer Science & Business Media, 2012.

[Gew96]

J. Geweke. Variable selection and model comparison in regression. Bayesian Statistics, 5:609–620, 1996.

[GG84]

Andrei Yur’evich Garnaev and Efim Davydovich Gluskin. The widths of a euclidean ball. In Doklady Akademii Nauk, volume 277, pages 1048–1052. Russian Academy of Sciences, 1984.

[GK23]

Roy Gotlib and Tali Kaufman. Nowhere to go but high: a perspective on highdimensional expanders. In International Congress of Mathematicians, pages 4842– 4871. European Mathematical Society-EMS-Publishing House GmbH, 2023.

[GM93]

Edward I George and Robert E McCulloch. Variable selection via gibbs sampling. Journal of the American Statistical Association, 88(423):881–889, 1993.

[Hop82]

John J Hopfield. Neural networks and physical systems with emergent collective computational abilities. Proceedings of the national academy of sciences, 79(8):2554– 2558, 1982.

[HRP+ 20]

Firas Hamze, Jack Raymond, Christopher A Pattison, Katja Biswas, and Helmut G Katzgraber. Wishart planted ensemble: A tunably rugged pairwise ising model with a first-order phase transition. Physical Review E, 101(5):052102, 2020.

[IR11]

Hemant Ishwaran and J. Sunil Rao. Consistency of spike and slab regression. Statistics & Probability Letters, 81(12):1920–1928, 2011.

[JM24]

Vishesh Jain and Clayton Mizgerd. Rapid Mixing of the Down-Up Walk on Matchings of a Fixed Size. In Approximation, Randomization, and Combinatorial Optimization. Algorithms and Techniques (APPROX/RANDOM 2024), volume 317 of Leibniz International Proceedings in Informatics (LIPIcs), pages 63:1–63:13, 2024.

[JMPV23]

Vishesh Jain, Marcus Michelen, Huy Tuan Pham, and Thuy-Duong Vuong. Optimal mixing of the down-up walk on independent sets of a given size. In 2023 IEEE 64th Annual Symposium on Foundations of Computer Science (FOCS), pages 1665–1681. IEEE, 2023.

[JS04]

Iain M. Johnstone and Bernard W. Silverman. Needles and straw in haystacks: Empirical bayes estimates of possibly sparse sequences. Annals of Statistics, 32(4):1594– 1649, 2004.

[JT17]

Aukosh Jagannath and Ian Tobasco. Some properties of the phase diagram for mixed p-spin glasses. Probability Theory and Related Fields, 167(3):615–672, 2017.

[Kas77]

Boris Sergeevich Kashin. Diameters of some finite-dimensional sets and classes of smooth functions. Izvestiya Rossiiskoi Akademii Nauk. Seriya Matematicheskaya, 41(2):334–351, 1977.

[Kaw66]

Kyozi Kawasaki. Diffusion constants near the critical point for time-dependent ising models. ii. Physical Review, 148(1):375, 1966.

[KLL+ 23]

Jonathan A. Kelner, Jerry Li, Allen Liu, Aaron Sidford, and Kevin Tian. Semirandom sparse recovery in nearly-linear time. In The Thirty Sixth Annual Conference

on Learning Theory, COLT 2023, volume 195 of Proceedings of Machine Learning Research, pages 2352–2398. PMLR, 2023. [KN26]

Seiichiro Kusuoka and Shuta Nakajima. A quantitative replica-symmetric bound of sherrington–kirkpatrick model in the entire de almeida–thouless region. arXiv preprint arXiv:2608.23413, 2026.

[KO20]

Tali Kaufman and Izhar Oppenheim. High order random walks: Beyond spectral gap. Combinatorica, 40(2):245–281, 2020.

[Kol18]

Vladimir Kolmogorov. A faster approximation algorithm for the gibbs partition function. In Conference On Learning Theory, COLT 2018, volume 75 of Proceedings of Machine Learning Research, pages 228–249. PMLR, 2018.

[KPPY25]

Aiya Kuchukova, Marcus Pappik, Will Perkins, and Corrine Yap. Fast and slow mixing of the kawasaki dynamics on bounded-degree graphs. Random Structures & Algorithms, 67(4), 2025.

[KSTZ25]

Symantak Kumar, Purnamrita Sarkar, Kevin Tian, and Yusong Zhu. Spike-and-slab posterior sampling in high dimensions. In The Thirty Eighth Annual Conference on Learning Theory, volume 291 of Proceedings of Machine Learning Research, pages 3407–3462. PMLR, 2025.

[Lit74]

William A Little. The existence of persistent states in the brain. Mathematical biosciences, 19(1-2):101–120, 1974.

[LMRW24]

Kuikui Liu, Sidhanth Mohanty, Amit Rajaraman, and David X. Wu. Fast mixing in sparse random ising models. In 65th IEEE Annual Symposium on Foundations of Computer Science, FOCS 2024, pages 120–128. IEEE, 2024.

[Lop26]

Patrick Lopatto. Replica symmetry up to the de almeida-thouless line in the sherrington-kirkpatrick model. arXiv preprint arXiv:2604.11921, 2026.

[LPW09]

David Asher Levin, Yuval Peres, and Elizabeth Wilmer. Markov Chains and Mixing Times. American Mathematical Society, 2009.

[Lyo89]

Russell Lyons. The ising model and percolation on trees and tree-like graphs. Communications in Mathematical Physics, 125(2):337–353, 1989.

[MB88]

Toby J Mitchell and John J Beauchamp. Bayesian variable selection in linear regression. Journal of the American Statistical Association, 83(404):1023–1032, 1988.

[Mos06]

Elchanan Mossel. Ising model on trees. Lecture notes for STAT 206A, University of California, Berkeley, 2006. Lecture 20.

[MS22]

Sumit Mukherjee and Subhabrata Sen. Variational inference in high-dimensional linear regression. Journal of Machine Learning Research, 23:1–56, 2022.

[MW26]

Andrea Montanari and Yuchen Wu. Provably efficient posterior sampling for sparse linear regression via measure decomposition. Journal of the American Statistical Association, pages 1–19, 2026.

[Opp18]

Izhar Oppenheim. Local spectral expansion approach to high dimensional expanders part i: Descent of spectral gaps. Discrete & Computational Geometry, 59(2):293–330, 2018.

[PC99]

Victor Y. Pan and Zhao Q. Chen. The complexity of the matrix eigenproblem. In Proceedings of the Thirty-First Annual ACM Symposium on Theory of Computing, pages 507–516. ACM, 1999.

[PF77]

Leonid A Pastur and Alexander L Figotin. Exactly soluble model of a spin glass. Soviet Journal of Low Temperature Physics, 3(6):378–383, 1977.

[PS19]

Nicholas G. Polson and Lei Sun. Bayesian ℓ0 -regularized least squares. Applied Stochastic Models in Business and Industry, 35(3):717–731, 2019.

[RG14]

Veronika Ročková and Edward I. George. Emvs: The em approach to bayesian variable selection. Journal of the American Statistical Association, 109(506):828–846, 2014.

[RG18]

Veronika Ročková and Edward I. George. The spike-and-slab lasso. Journal of the American Statistical Association, 113(521):431–444, 2018.

[Roc18]

Veronika Rockova. Bayesian estimation of sparse signals with a continuous spike-andslab prior. Annals of Statistics, 46(1):401–437, 2018.

[RS22]

Kolyan Ray and Botond Szabó. Variational bayes for high-dimensional linear regression with sparse priors. Journal of the American Statistical Association, 117(539):1270–1281, 2022.

[RW26]

Amit Rajaraman and David X Wu. Markov chains approximate message passing. In Proceedings of the 58th Annual ACM Symposium on Theory of Computing, pages 1192–1199, 2026.

[SK75]

David Sherrington and Scott Kirkpatrick. Solvable model of a spin-glass. Physical review letters, 35(26):1792, 1975.

[Str69]

Volker Strassen. Gaussian elimination is not optimal. 13:354–356, 1969.

[Tal10]

Michel Talagrand. Mean field models for spin glasses: Volume I: Basic examples, volume 54. Springer Science & Business Media, 2010.

[Tal11]

Michel Talagrand. Mean Field Models for Spin Glasses: Volume II: Advanced ReplicaSymmetry and Low Temperature, volume 55. Springer Science & Business Media, 01 2011.

[Tro12]

Joel A Tropp. User-friendly tail bounds for sums of random matrices. Foundations of computational mathematics, 12(4):389–434, 2012.

[Ver18]

Roman Vershynin. High-dimensional probability: An introduction with applications in data science, volume 47. Cambridge university press, 2018.

[Wai19]

Martin J Wainwright. High-dimensional statistics: A non-asymptotic viewpoint, volume 48. Cambridge university press, 2019.

Numerische Mathematik,

[Wei07]

Pierre Weiss. L’hypothèse du champ moléculaire et la propriété ferromagnétique. Journal de Physique Théorique et Appliquée, 6(1):661–690, 1907.

[WJ08]

Martin J Wainwright and Michael I Jordan. Graphical models, exponential families, and variational inference. Foundations and Trends® in Machine Learning, 1(1-2):1– 305, 2008.

[Wu06]

Liming Wu. Poincaré and transportation inequalities for gibbs measures under the dobrushin uniqueness condition. The Annals of Probability, 34(5), 2006.

[Yan52]

Chen Ning Yang. The spontaneous magnetization of a two-dimensional ising model. Physical Review, 85(5):808, 1952.

[YWJ16]

Yun Yang, Martin J Wainwright, and Michael I Jordan. On the computational complexity of high-dimensional bayesian variable selection. The Annals of Statistics, 2016.

Sampling Near the Almeida–Thouless Line

In this section, we give an application to sampling near the Almeida-Thouless (AT) line. The AT line is parameterized by β > 0 and a field strength h > 0, and considers the specialized SK model β ⊤ ⊤ πβ,h (x) ∝ exp x Jx + h1d x , x ∈ X d , (60) 2 where J ∼ GOE(d) follows Model 1, and we set the diagonal of J to zero without loss of generality. Thus, h denotes the coefficient of the linear term ∝ 1d . The AT line delineates a region in R2≥0 , representing a pair of parameters (β, h) in (60). This line was originally derived in [dAT78] via the replica method, with the prediction that models induced by h above the line (i.e., large enough as a function of β) are replica symmetric, and that pairs below the line exhibit replica symmetry breaking. This prediction was recently established rigorously for the SK model with a homogeneous external field [Lop26]. However, our application requires a stronger quantitative concentration estimate, so our results hold above the slightly more restrictive weak AT line (Definition 3), leveraging bounds by [Tal11, JT17] as presented by [RW26]. For large h, (60) favors x with more positive spins. Following this convention, we write n o d,− := x ∈ X d : N− (x) ≤ k . N− (x) := |{i ∈ [d] : xi = −1}| , X≤k Our earlier results use the number of positive spins as the sparsity parameter; the two conventions are equivalent under the global spin flip x → −x. Moreover, on every fixed-magnetization slice, h1⊤ d x is constant, so all fixed-size mixing analyses for (60) are independent of h. We now define the regions determined by the AT line. Definition 2 (AT condition, [dAT78]). For β, h > 0, let q = q(β, h) be the unique solution to √ q = E tanh2 (β qZ + h) , Z ∼ N (0, 1). Define

√ αAT (β, h) := β 2 E sech4 (β qZ + h) .

We say that (β, h) satisfies the AT condition if αAT (β, h) ≤ 1. The AT boundary hAT (β) is characterized by αAT (β, hAT (β)) = 1. We also define qAT (β) := q(β, hAT (β)). Definition 3 (Weak AT condition, [BY22]). With q = q(β, h) as in Definition 2, define √ αwAT (β, h) := β 2 E sech2 (β qZ + h) = β 2 (1 − q). We say that (β, h) satisfies the weak AT condition if αwAT (β, h) ≤ 1. The weak AT boundary hwAT (β) is characterized by αwAT (β, hwAT (β)) = 1. We also define qwAT (β) := q(β, hwAT (β)) Since sech4 (u) ≤ sech2 (u) for all u, we have that αAT (β, h) ≤ αwAT (β, h) for all (β, h). Thus, the weak AT region (above the weak AT line) is a subset of the AT region.

A.1

Preliminaries

In this section, we prove two key preliminary results. The first (Lemma 23) computes asymptotics of the boundaries hAT , hwAT as a function of β. The second (Lemma 24) formalizes a reduction from large field strength h to concentration on high-magnetization states. Boundary asymptotics. We first recall a common Laplace estimate that will be used repeatedly. 2

Lemma 22. Let ϕ(t) := √12π exp(− t2 ) be the Gaussian density, and let p ∈ {2, 4}. Suppose that h

there are positive {qβ , hβ }β>0 satisfying, as β → ∞, qβ → 1, bβ := ββ → ∞, and ββ → 0. Then √ E sechp β qβ Z + βbβ = where

Ip √ ϕ β qβ

bβ √ qβ

(1 + o(1)) ,

(61)

sech (u) du = 2,

I2 =

I4 =

4 sech4 (u) du = . 3 R

√ Proof. We first perform a change of variables u = β( qβ z + bβ ), which gives √ Z ϕ b / qβ bβ u u2 β √ p p − 2 E sech β qβ Z + βbβ = sech (u) exp du. √ β qβ βqβ 2β qβ R b

b u

(62)

Since qβ → 1 and ββ → 0, the factor exp( βqβ β − 2βu2 q ) converges pointwise to one. Moreover, for β

sufficiently large β, βqββ → 0, so we have exp

bβ u u2 − 2 βqβ 2β qβ

≤ exp

|u| 2

For p ∈ {2, 4}, sechp (u) exp( |u| 2 ) is integrable. Applying dominated convergence then gives (61). We now derive the (identical) asymptotics of hAT and hwAT . Lemma 23. As β → ∞, hAT (β) =

√

p 1 2 β log β 1 + O , log β

hwAT (β) =

√

p 1 2 β log β 1 + O . log β

Proof. We begin by verifying the conditions of Lemma 22. For ⋆ ∈ {AT, wAT}, write q⋆ := q⋆ (β),

b⋆ :=

h⋆ (β) . β

We also drop the index β from these two sequences for simplicity. First, along the weak AT boundary, 1 − qwAT = β −2 , so qwAT → 1. Similarly, along the AT boundary, 1 √ E sech4 (β qAT Z + hAT (β)) = 2 . β √ Since 1 − qAT = E[sech2 (β qAT Z + hAT (β))], Cauchy-Schwarz gives 1 − qAT ≤ β −1 , so qAT → 1. 52

Next, for either value of ⋆, the boundary equation also implies b⋆ → ∞. Indeed, if b⋆ were bounded along a subsequence, then restricting the integral (62) to u ∈ [−1, 1] would already give a lower bound of order β −1 , because q⋆ → 1, b⋆ is bounded, and sechp (u) = Ω(1) for u ∈ [−1, 1]. This would contradict the boundary equations, which say that (62) evaluates to β −2 . Finally, we claim that bβ⋆ → 0. Otherwise, b⋆ ≥ ϵβ along a subsequence for some ϵ > 0. Now, √ √ split the integral (62) along the events Z ≥ −b⋆ /2 q⋆ or Z ≤ −b⋆ /2 q⋆ . In the former region, 2 the change of variables gives u ≥ ϵβ2 , so the corresponding integral is exp(−Ω(β 2 )). Similarly, the latter region has a probability bounded by exp(−Ω(β 2 )), and sech is pointwise bounded. Thus, the entire integral is exp(−Ω(β 2 )), again contradicting that it equals β −2 by definition. Lemma 22 therefore applies and we obtain that √ 3 qAT bAT ϕ √ = (1 + o(1)) , qAT 4β

bwAT √ qwAT

√

qwAT (1 + o(1)) . 2β

Expanding with the definition of ϕ, and then taking logarithms, then gives for ⋆ ∈ {AT, wAT}, b2⋆ = log β + O(1) =⇒ b2⋆ = 2q⋆ log β + O(1) = 2 log β − 2(1 − q⋆ ) log β + O(1). 2q⋆ Since 1 − qAT ≤ β −1 and 1 − qwAT = β −2 , the desired claims follow: p 1 . b2⋆ = 2 log β + O(1) =⇒ b⋆ = 2 log β 1 + O log β

Magnetization from field strength. To conclude the section, we show that taking h large in (60) implies that πβ,h is concentrated on high-magnetization states. For s ∈ [0, 1], we let H2 (s) := −s log s − (1 − s) log(1 − s) denote the binary entropy, with the convention 0 log 0 = 0. Lemma 24. Fix δ ∈ (0, 12 ). With probability at least 1 − δ over J in Model 1, the following holds simultaneously for every ρ ∈ (0, 12 ) and γ ≥ 0. If s H2 (ρ) 2(1 − ρ) 1 d h≥ +β H2 (ρ) + log + γ, (63) 2ρ ρ d δ then Px∼πβ,h [N− (x) > ρd] ≤ d exp (−2γρd) . Proof. For a fixed S ⊆ [d], define YS := YS ≥ −t

|S| d

i∈S Jij . We claim that simultaneously for all S ⊆ [d], j ∈S /

d 1 , where t(s) := d 2s(1 − s) H2 (s) + log , d δ

(64)

with probability ≥ 1 − δ. To see this, fix some r ∈ [d]. Under Model 1, we have YS ∼ N (0, ds(1 − s)) for all |S| = r and s := dr . Then the standard Gaussian tail bound gives δ P[−YS ≥ t(s)] ≤ exp (−dH2 (s)) . d 53

Since dr ≤ exp(dH2 (s)), we conclude that (64) holds for all |S| = r with probability ≥ 1 − dδ , and then a union bound over all nontrivial layers r ∈ [d] gives the claim. w(S) Next, letting w be the unnormalized weight in (60), a direct calculation gives log w(∅) = −2h|S| − r 2βYS . Thus, again letting s = d for some r ∈ [d], P H2 (s) βt(s) S:|S|=r w(S) ≤ exp 2sd + −h . (65) w(∅) 2s sd

It remains to bound the right-hand side uniformly over s ≥ ρ. First, we can directly check that log(1−s) d H2 (s) < 0, so H2s(s) is decreasing in s. This implies ds ( s ) = s2 s 1 2(1 − s) d t(s) H2 (s) + log = sd s d δ is decreasing in s, because every term in the square root is decreasing. Hence for every s ≥ ρ, the assumed lower bound (63) implies H2 (s) βt(s) H2 (ρ) βt(ρ) + −h≤ + − h ≤ −γ. 2s sd 2ρ ρd Summing (65) over all r > ρd and using that the partition function is ≥ w(∅) proves the claim.

A.2

Sampling around 1d

In this section, we give the basic variant of our result for sampling from (60). This variant shows that √ as β → ∞, when the field strength h exceeds the threshold hAT (β) in Lemma 23 by roughly a 2 factor, we can sample from πβ,h in polynomial time. Theorem 6. Let β ≥ 1 and δ ∈ (0, 12 ), and suppose that (31) holds with δ ← 2δ . If s H2 (ρβ ) 1 − ρβ 1 2d 1 2d h≥ +β 2 H2 (ρβ ) + log + log , 2ρβ ρβ d δ 2ρβ d δ

(66)

c for a universal constant c > 0, then with probability at least 1 − δ over the SK where ρβ := β 2 log(eβ) model (Model 1), there is a polynomial-time algorithm that outputs x satisfying TV (Law(x), πβ,h ) ≤ δ. Furthermore, for δ = poly( d1 ) and fixed β, as d → ∞, the right-hand side of (66) converges to s p H2 (ρβ ) H2 (ρβ ) hHM (β) := + β 2(1 − ρβ ) = 2β log β(1 + o(1)). (67) 2ρβ ρβ

Proof. Throughout this proof, set k := ⌊ρβ d⌋,

γ :=

1 2d log . 2ρβ d δ

Applying Lemma 24 with failure probability 2δ , sparsity lower bound parameter ρβ , and γ as defined above then implies that with probability ≥ 1 − 2δ over J, δ Px∼πβ,h (N− (x) > k) ≤ d exp(−2γρβ d) ≤ . 2 54

Next, to sample from πβ,h conditioned on N− (x) ≤ k, assume k ≥ 1, else it suffices to output 1d . By choosing c sufficiently small, k ≤ ρβ d is in the range required by Corollary 5. Then applying Corollary 5 with accuracy and failure probability 2δ gives the desired sample x, upon flipping 1s and −1s consistently with our convention in this section. A union bound then gives both the total failure probability of δ over the draw J, and the overall accuracy δ to the target πβ,h . Finally, we prove the asymptotic claims. For fixed β and δ = poly( d1 ), the finite-d corrections in (66) vanish as d → ∞, so the limit is s H2 (ρβ ) H2 (ρβ ) hHM (β) = + β 2(1 − ρβ ) . 2ρβ ρβ The conclusion follows because as β → ∞, we have H2 (ρβ ) 1 = log + 1 + O(ρβ ), ρβ ρβ

A.3

log

1 = 2 log β + o(log β). ρβ

Sampling around the mean

√ We next give a stronger variant of Theorem 6 that removes the 2 factor overhead in our lower bound on h, assuming the ability to compute the signs of the mean magnetization vector. Concretely, Theorem 6 is not adapted to the actual realization of J. Instead, we show that knowledge of τi⋆ := sign(mi ), where mi := Eπβ,h [xi ] for all i ∈ [d], for a fixed J, where sign(0) := 1, allows us to recenter the algorithm and improve our h range. Roughly speaking, the idea is to use overlap concentration to redefine our notion of sparsity as disagreement with τi⋆ , rather than disagreement with 1d as used in Section A.2. We next set up some notation for this section. For τ ∈ X d , let Nτ (x) := |{i ∈ [d] : xi ̸= τi }| . ⊗2 Also, for two i.i.d. draws (x(1) , x(2) ) ∼ πβ,h , we define the overlap quantity

R1,2 :=

x(1) , x(2) . d

Observe that independence gives E[R1,2 ] = d1 ∥m∥22 , which is a quantity depending on the realized J. We next state a stronger result from [RW26] which shows that above the weak AT line, R1,2 concentrates around the deterministic quantity q(β, h), independent of J. Proposition 5 (Lemma 4.14, [RW26]). Suppose that (β, h) satisfies the weak AT condition in Definition 3. Then there exists Cβ,h > 0 such that !# " d (R1,2 − q(β, h))2 ≤ 2. (68) EJ Eπ⊗2 exp β,h Cβ,h 55

As a corollary, we upgrade Proposition 5 into a high-probability distance bound to τ ⋆ , replacing the direct sparsity notion from Lemma 24. Our strategy is to first relate the random quantity ∥m∥22 to q(β, h) using Proposition 5, and then to relate ∥m∥22 to the Hamming distance Nτ ⋆ using (70). Lemma 25. Suppose that (β, h) satisfies the weak AT condition. Then there exists Cβ,h > 0 such that for every ρ, δ ∈ (0, 12 ), with probability at least 1 − δ over J, s 1 − q(β, h) 1 Cβ,h log 2δ πβ,h (Nτ ⋆ > ρd) ≤ + . (69) 2ρ 2ρ d Proof. Set qd (J) :=

∥m∥22 ⊗2 [R1,2 ] = qd (J). We first derive d , and recall that Eπβ,h d

1 1 X 1 Eπβ,h [Nτ ⋆ (x)] = (1 − |mi |) ≤ (1 − qd (J)) . d 2d 2

(70)

i=1

Applying Jensen’s inequality conditionally on J to (68) then gives !# " d (qd (J) − q(β, h))2 ≤ 2. EJ exp Cβ,h Hence, for every t > 0, dt2 PJ [qd (J) < q(β, h) − t] ≤ 2 exp − Cβ,h

On the complementary event, (70) and Markov’s inequality under πβ,h give πβ,h (Nτ ⋆ > ρd) ≤ q Taking t =

Cβ,h log 2δ d

1 − q(β, h) + t . 2ρ

proves the claim.

Lemma 25 shows that to apply our bounded-magnetization SK sampler (Corollary 5), we have reduced the problem to making 1 − q(β, h) smaller than the sparsity parameter ρβ ≈ (β 2 log β)−1 . The weak AT condition only gives 1 − q(β, h) ≤ β −2 , so we choose a slightly larger field strength, still asymptotic to the weak AT scale in Lemma 23: p hOC (β) := β 2 log β + 2 log log β + 2 log log log β, qOC (β) := q (β, hOC (β)) . (71) Lemma 26 next shows that the field strength in (71) satisfies the weak AT condition, and thus we can bound the sparsity of its induced model using Lemma 25. Lemma 26. There exist universal constants C, β0 > 0 such that, for every β ≥ β0 , and h ≥ hOC (β), 1 − q(β, h) ≤

C β 2 log β log log β

Consequently, (β, h) satisfies the weak AT condition. 56

(72)

Proof. We first claim that q(β, h) is nondecreasing in h if (β, h) strictly satisfies the weak AT √ condition, so it suffices to prove the result for h = hOC (β). Let T (q, h) := E tanh2 (β qZ + h), so that by definition, q = T (q, h). Then, performing a Gaussian integration by parts gives ∂ √ √ T (q, h) = β 2 E sech2 (β qZ + h) 1 − 3 tanh2 (β qZ + h) ≤ β 2 (1 − q) < 1, ∂q where the last inequality used the weak AT condition. Similarly, ∂ √ √ T (q, h) = 2E tanh (β qZ + h) sech2 (β qZ + h) > 0, ∂h √ because tanh(·) sech2 (·) is odd and positive on R>0 , and β qZ + h is centered around the positive value h > 0. Now, implicit differentiation gives q ′ (h) =

∂ T (q(h), h) ∂ ∂ T (q(h), h)q ′ (h) + T (q(h), h) =⇒ q ′ (h) = ∂h ∂ > 0, ∂q ∂h 1 − ∂q T (q(h), h)

as desired. For the rest of the proof we take h = hOC (β). Write bβ := hOCβ(β) . We first verify Lemma 22’s hypotheses. From (71), the conditions bβ → ∞ and bβ β → 0 are immediate. Moreover, from the fixed-point equation defining qβ := qOC (β),

1 − qβ = E sech

√

β qβ Z + βbβ

≤ sech

βbβ 2

bβ +P Z <− → 0, 2 b

where the inequality split the expectation based on whether Z ≥ − 2β or not, and used that sech2 ≤ 1 pointwise and qβ ≤ 1. Thus, Lemma 22 applies with p = 2, yielding bβ 2 C 1 − qβ = √ ϕ √ (1 + o(1)) ≤ ϕ(bβ ), β qβ qβ β for a universal constant C, where the last inequality used that qβ ≤ 1, ϕ is decreasing on R≥0 , and qβ → 1 as β → ∞. Finally, substituting the definition (71) gives ϕ(bβ ) = √

1 , 2πβ log β log log β

C and combining the above two displays proves (72). Finally, β 2 (1 − qβ ) ≤ log β log log β < 1 as β → ∞, which verifies the weak AT condition with strict inequality. Since we earlier showed q ′ (h) > 0 whenever β 2 (1 − q(h)) < 1, all h > hOC (β) also strictly satisfy the weak AT condition.

We are finally ready to give our main result. Theorem 7. Let β ≥ 1 and δ ∈ (0, 12 ), and suppose that (31) holds with δ ← 2δ . Assume that c h ≥ hOC (β) defined in (71), and define ρβ := β 2 log(eβ) for a universal constant c > 0. Also assume log log β ≥

2C , δ

d≥

4Cβ,h log 4δ , δ 2 ρ2β

(73)

where C is a universal constant and Cβ,h is from Lemma 25, and that we are given τi⋆ := sign Eπβ,h [xi ] for all i ∈ [d]. Then with probability at least 1−δ over the SK model (Model 1), there is a polynomial-time algorithm that outputs x satisfying TV (Law(x), πβ,h ) ≤ δ. Further, √ p hOC (β) = 2β log β (1 + o(1)) . Proof. Let D := diag (τ ⋆ ) throughout. We first claim that the fixed-magnetization sampler in Theorem 4 (and hence, the bounded-magnetization sampler in Corollary 5) is invariant to replacing J ← DJD. To see this, all but one of the estimates in Lemma 12 remain unchanged under this replacement, because Du has the same ℓ2 and ℓ1 norms as u for any vector u, and (DJD)◦(DJD) = J ◦ J. The only difference is that ∥DJD1 √ d ∥∞ could be larger because D1d is no longer independent of J, but the same bound holds up to a d factor. This factors into Theorem 4’s initial χ2 divergence bound, which only affects the claimed runtime by a polynomial factor after taking logarithms. To complete the proof, take k := ⌊ρβ d⌋. Lemmas 25 and 26, and our assumed bounds (73), imply δ Px∼πβ,h [Nτ ⋆ (x) > k] ≤ . 2 The rest of the proof is identical to Theorem 6, under the change of variables x ← Dx. √ In summary, Theorem 7 covers a range of h with a lower bound hOC (β) roughly a 2 factor smaller than Theorem 6’s hHM (β). In particular, up to a 1 + o(1) factor, hOC (β) matches the AT and weak AT thresholds hAT (β) and hwAT (β) derived in Lemma 23. In comparison, a recent work by [BAR26] derives a similar polynomial-time sampling result, but examining their proof (particularly, Corollary 1.2 combined with the √ improvement in Remark 3.3) implies their threshold on the field strength h scales as h = Ω(β 2 log β), i.e., roughly a β factor larger than hAT (β). We remark that Theorem 7 is a conditional result that requires access to τ ⋆ , the signs of the mean magnetization vector; we leave the efficient computation of τ ⋆ as an important open problem. Additionally, our definition of hOC (β) includes potentially unnecessary low-order terms, arising due to a discrepancy between the weak AT region’s definition (which gives 1 − q(β, h) ≤ β −2 ) with the requirements of our sampler in Corollary 5 (which requires a sparsity parameter ρ ≈ (β 2 log β)−1 ). It is also an interesting open problem to improve the thresholds imposed by our approach, with the goal of efficient sampling across the entire AT region.

Scaling of Infinite ∆-Regular Tree Threshold

In this section, we provide a calculation that explains asymptotics induced by the “tree threshold” ∆ βc (∆), which parameterizes the results of [CDKP22, KPPY25]. Concretely, βc (∆) = log( ∆−2 ) is a critical inverse temperature under which the Ising model on the infinite ∆-regular tree undergoes a phase transition. When β < βc (∆), there is a unique fixed point of a certain message-passing recursion on the infinite tree, and when β > βc (∆), two new fixed points appear. We defer an overview to [Lyo89, Mos06]; here we calculate consequences of this threshold asymptotically.

The main algorithmic result of [CDKP22] at low temperatures β > βc (∆) is their Theorem 2(a). + They show that there is a corresponding critical magnetization η∆,β,1 at which the fixed-magnetization + Ising model undergoes a computational phase transition. This η∆,β,1 is the expected spin at the root of the infinite ∆-regular tree, at one of the new fixed points emerging when β > βc (∆) (see Section 4, [CDKP22] for these calculations). Note that [CDKP22] defines the magnetization η ∈ [−1, 1] as (in our notation) −1 + 2k d , so that k = d corresponds to η = 1 (and k = 0 corresponds to η = −1). Rearranging, Theorem 2(a) of [CDKP22] applies at sparsity levels + η < −η∆,β,1 =⇒ k <

d + 1 − η∆,β,1 . 2

+ It remains to understand 1 − η∆,β,1 . From Section 1.1, [CDKP22], letting L be the largest root of

β , L = (∆ − 1) arctanh tanh(L) tanh 2 + denoting ηβ := η∆,β,1 for short,

β ηβ = tanh L + arctanh tanh(L) tanh . 2 1 Let A := arctanh(tanh(L) tanh( β2 )) = ∆−1 L. Then,

β β tanh A tanh A = tanh ((∆ − 1)A) tanh =⇒ tanh = . 2 2 tanh ((∆ − 1)A) Now using the approximation tanh(c) = 1 − Θ(exp(−2c)) for large c, we obtain A = β2 + O∆ (1). Finally, plugging this back into our definition of ηβ , we have 1 − ηβ = 1 − tanh (∆A) = Θ (exp (−2∆A)) = exp (−∆β + O∆ (1)) . Therefore, for the sparsity regime k < d exp (−∆β + O∆ (1))

(74)

to be meaningful (i.e., the inequality above does not hold only when k = 0), the result of [CDKP22] is limited to inverse temperatures of β = O(log d) for constant-degree graphs. The main algorithmic result of the subsequent work [KPPY25] is their Theorem 1.1, which is parameterized at a slightly different critical magnetization ηβ,a , which is at least as large as the ηβ from before. Therefore, Theorem 1.1 of [KPPY25] is also restricted to the regime (74). Finally, we note that the comparison in this section is purely a statement about the allowable temperatures that our sparsity-aware framework tolerates (vs. the prior works [CDKP22, KPPY25]), for our models of interest. For example, directly applying our results to the graph-based Ising models in these prior works only permits rapid mixing in the high-temperature regime β = O( k1 ). However, for other well-studied models, e.g., Models 1 and 2, our results allow taking inverse temperatures as large as β = poly(d) when k is sufficiently small.

Record · ID 668045 · SHA-256 17264fce081f1d86
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.