Gradient-Free Sampling from Generative Models via Stochastic Bounded Extremum Seeking

Alexander Scheinker

stat.MLOct 3, 2026

Gradient-Free Sampling from Generative Models via Stochastic Bounded Extremum Seeking

Authors: Alexander Scheinker

Organizations: Instrumentation and Controls Group Los Alamos National Laboratory Los Alamos, NM 87544, USA

Abstract

We introduce a sampling approach for energy- and score-based generative models that requires no gradient evaluations of the model. Replacing the drift term that would normally contain the score $\nabla_\mathbf{x} \log p_θ(\bf{x})$ with a high-frequency dithered cosine of the model's value, $\sqrt{αω}\,\cos(ωt + k \log p_θ(\bf{x}))$ , produces, in the high-frequency averaging limit, Langevin Markov chain Monte Carlo for energy-based models and the reverse-time SDE of score-based diffusion. We prove that trajectories of the dithered Itô SDE converge to those of the target SDE, driven by the same Brownian motion, uniformly on compact time intervals in probability, by an averaging argument that extends bounded extremum seeking (ES) to Itô processes, with an explicit $O(ω^{-1/2})$ mean-square rate under global bounds. The approach is not confined to smooth targets: it extends to $C^{1,1}$ energies with discontinuous curvature (without ellipticity requirement) and to Sobolev energies whose Hessians exist only off measure zero sets; for Lipschitz energies with gradient kinks the averaged limit remains well posed; the Krylov-Röckner integrability class is the boundary of provability. The approach provides a hard a priori bound on the per-step update rate and applies to explicitly time-varying targets on finite horizons. Gradient-free pixel-space sampling is not competitive with well-tuned backpropagation-based samplers at practical evaluation budgets; the regime where the approach offers an advantage is latent-space sampling when the model is a black box and the target drifts in time. We demonstrate latent-space tracking for time-varying images on CelebA-HQ ( $256{\times}256$ ) from limited 1D projection measurements and latent-space EBM sampling on CIFAR-10.

Figures & tables

Figure 1: Bounded ES dynamics results in various samplers depending on the choice of $(\mathbf{f},\bm{\sigma},k,\alpha,J)$ , where $\bm{\sigma}(\mathbf{x},t)=\operatorname{diag}\!\big(g_{1}(\mathbf{x},t),\dots,g_{d}(\mathbf{x},t)\big)$ , $\mu_{\omega,i}(\mathbf{x},t)=f_{i}(\mathbf{x},t)\;+\;g_{i}(\mathbf{x},t)\,\sqrt{\alpha\omega_{i}}\,\cos\!\big(\omega_{i}t+k\,J(\mathbf{x},t)\big)$ . Averaged dynamics not changed by rough energies, but the provable results are. Above $C^{2}$ means the a.e.-equivalence class of $J$ contains a $C^{2}$ function. $J$ can have measure zero jumps.

Instantiation	$\mathbf{f}$	$\bm{\sigma}$	$k$	$\alpha$	averaged SDE / stationary law
(I) EBM Langevin	$\mathbf{0}$	$\sqrt{2}\,I$	$1$	$1$	$d\bar{\mathbf{x}}=-\nabla\mathcal{E}_{\theta}\,dt+\sqrt{2}\,d\mathbf{w}$ ; $\rho_{\infty}\propto e^{-\mathcal{E}_{\theta}}$
(I λ ) tempered	$\mathbf{0}$	$\sqrt{2}\,I$	$\lambda$	$1$	$d\bar{\mathbf{x}}=-\lambda\nabla\mathcal{E}_{\theta}\,dt+\sqrt{2}\,d\mathbf{w}$ ; $\rho_{\infty}\propto e^{-\lambda\mathcal{E}_{\theta}}$
(II) reverse SDE ( $\tau{=}T{-}t$ )	$-\mathbf{f}(\mathbf{x},T{-}\tau)$	$g(T{-}\tau)I$	$-1$	$2$	$d\bar{\mathbf{x}}=[-\mathbf{f}+g^{2}\nabla\log p_{T-\tau}]d\tau+g\,d\mathbf{w}$
(III) classifier	$\mathbf{0}$	$\sqrt{2}\,I$	$1$	$1$	(I) with $\mathcal{E}_{\theta}(\mathbf{x}\,\|\,y)=-f_{\theta}(\mathbf{x})[y]$

Table 1: Parameter dictionary: choices of

(\mathbf{f},\bm{\sigma},k,\alpha,J)

in Eq. ( 4 )–Eq. ( 5 ) and the resulting averaged dynamics Eq. ( 6 ).

\mathcal{E}_{\theta}

denotes an energy,

p

a density,

\lambda>0

a tempering exponent.

Figure 2: Latent-space tracking from 1D RGB projections of time-varying images.

NMSE	$t{=}1170$	$2340$	$3510$	$4680$	$5850$	$7020$	$8191$
Val, drift (vs. $t{=}0$ )	$.025{\pm}.0006$	$.078{\pm}.0022$	$.146{\pm}.0040$	$.221{\pm}.0061$	$.305{\pm}.0086$	$.400{\pm}.0119$	$.510{\pm}.0167$
Val, bounded ES	$.010{\pm}.0002$	$.019{\pm}.0004$	$.029{\pm}.0005$	$.040{\pm}.0007$	$.052{\pm}.0009$	$.067{\pm}.0013$	$.083{\pm}.0017$

Table 2: Tracking study:

500

held-out pairs

\times

3

seeds, mean

\pm

SEM over the

1{,}500

runs.

Metric	Validation ( $500$ pairs)
PSNR (dB) $\uparrow$	$21.08\pm 3.30$
SSIM $\uparrow$	$0.855\pm 0.095$
FID $\downarrow$	$8.68\pm 0.11$
KID $\times 10^{3}$ $\downarrow$	$5.23\pm 0.47$

Table 3: Tracked vs. target image quality (left), and its evolution over the tracking horizon (right).

Method	Meas./update	Updates	Final NMSE (mean $\pm$ std)	median (IQR)
Bounded ES (ours)	$1$	$8192$	$\mathbf{0.080\pm 0.035}$	$\mathbf{0.072}$ ( $0.055$ – $0.101$ )
SPSA	$2$	$4096$	$0.142\pm 0.059$	$0.131$ ( $0.104$ – $0.181$ )
FD (central)	$512$	$16$	$0.276\pm 0.312$	$0.219$ ( $0.169$ – $0.307$ )
GS-ES (antithetic)	$2m$	$8192/2m$	$1.168\pm 0.627$	$1.028$ ( $0.773$ – $1.348$ )
AutoDiff oracle	—	—	$0.048\pm 0.032$	$0.041$ ( $0.024$ – $0.063$ )

Table 4: Final image-NMSE at

t{=}8191

\pm

is the std over all

219\times 3

runs; median (IQR) across the

219

pairs after seed averaging (the per-pair distributions are right-skewed). “Updates” is the number of parameter updates each method uses,

m

is the tuned number of GS-ES antithetic pairs per update.

Panel of Figure 3 — Figure 3: ES dynamics results for tracking time-varying CIFAR-10 images based on 1D projections.

Tracker	$N_{s}$	$N$	L1 final $\downarrow$	RMSE $\downarrow$	SSIM $\uparrow$	FID (floor) $\downarrow$
Frozen	–	5000	0.568 $\pm$ 0.151	0.342 $\pm$ 0.082	0.046 $\pm$ 0.068	9.95 (9.95)
Bounded ES	200	5000	0.160 $\pm$ 0.070	0.110 $\pm$ 0.044	0.556 $\pm$ 0.130	34.5 (9.95)

Table 5: Final errors

(\mu\pm\sigma)

between tracked and true images, and FID between final tracked and final true images at the sample size

N

of the run, with the same-

N

floor (FID between two independent samples of the true distribution). L1 on

[-1,1]

pixel scale, RMSE/SSIM on

[0,1]

images.

Regularity of $J(\cdot,t)$	Permits	Drift $\nabla J$	Status
$C^{2}$ (Assumption 5 )	classical smoothness	continuous	Thm 1 ; rate $\omega^{-1/2}$ (Thm 3 )
$C^{1,1}=W^{2,\infty}$	bounded, non- $C^{1}$ Hessian	Lipschitz	Prop. 1 (proved)
$W^{2,p}_{\mathrm{uloc}}$ , $p>d$ , $\nabla J$ bounded	unbounded $L^{p}$ Hessian	continuous	Thm 2 (proved)
Lipschitz, $W^{1,\infty}$	gradient kinks	a.e., $\in L^{\infty}$	Prop. 3 ; convergence, no rate
$\nabla J\in L^{q}_{t}L^{p}_{\mathbf{x}}$ , $\frac{d}{p}{+}\frac{2}{q}<1$	unbounded gradient; non-Lipschitz cusps	$\in L^{p}$	limit well posed (Krylov– Röckner); averaging expected
$\nabla J$ distributional, $C^{-\beta}$	jumps in $J$	surface measure	open frontier

Table 6: Provability degrades gracefully as

J

roughens. Statements and proofs in Appendix B .

Appendix figures & tables17 assets

Supplementary material from the paper’s appendix.

Appendix

Figure 4: Latent-space tracking from 1D RGB projections of time-varying images. The top row shows the true images, the second row shows decent tracking, and the bottom row shows catastrophic instability and failure when the ES gain $k\alpha>0$ is extremely large.

$d$	$c=5$	$c=10$	$c=20$	$c=100$
2	$-0.410\pm 0.026$	$-0.489\pm 0.048$	$-0.487\pm 0.057$	$-0.489\pm 0.056$
4	$-0.494\pm 0.035$	$-0.587\pm 0.042$	$-0.594\pm 0.043$	$-0.596\pm 0.038$
8	$-0.574\pm 0.020$	$-0.655\pm 0.014$	$-0.650\pm 0.014$	$-0.655\pm 0.015$
16	$-0.652\pm 0.005$	$-0.705\pm 0.013$	$-0.703\pm 0.018$	$-0.701\pm 0.018$
32	$-0.641\pm 0.006$	$-0.645\pm 0.005$	$-0.649\pm 0.006$	$-0.652\pm 0.006$

Appendix

Table 7: Fitted log–log slope of the coupled sup-error vs.

\omega

;

\sigma=0

(deterministic ES vs. gradient flow), window

\omega\in[10,10^{3}]

N=1.75

. Mean

\pm

std over 5 seeds.

$d$	$c=5$	$c=10$	$c=20$	$c=100$
2	$-0.395\pm 0.037$	$-0.499\pm 0.050$	$-0.504\pm 0.055$	$-0.448\pm 0.058$
4	$-0.465\pm 0.023$	$-0.550\pm 0.050$	$-0.552\pm 0.041$	$-0.571\pm 0.023$
8	$-0.578\pm 0.025$	$-0.653\pm 0.016$	$-0.661\pm 0.021$	$-0.657\pm 0.024$
16	$-0.646\pm 0.014$	$-0.705\pm 0.019$	$-0.709\pm 0.014$	$-0.720\pm 0.017$
32	$-0.635\pm 0.008$	$-0.642\pm 0.007$	$-0.645\pm 0.007$	$-0.646\pm 0.006$

Appendix

Table 8: As Table 7 , for

\sigma=0.1

$d$	$c=5$	$c=10$	$c=20$	$c=100$
2	$-0.387\pm 0.031$	$-0.459\pm 0.021$	$-0.466\pm 0.009$	$-0.452\pm 0.011$
4	$-0.448\pm 0.018$	$-0.500\pm 0.013$	$-0.493\pm 0.010$	$-0.474\pm 0.016$
8	$-0.522\pm 0.011$	$-0.526\pm 0.012$	$-0.514\pm 0.014$	$-0.527\pm 0.009$
16	$-0.558\pm 0.006$	$-0.550\pm 0.007$	$-0.532\pm 0.010$	$-0.539\pm 0.008$
32	$-0.544\pm 0.003$	$-0.530\pm 0.004$	$-0.509\pm 0.009$	$-0.511\pm 0.007$

Appendix

Table 9: As Table 7 , for

\sigma=1

(the sampling setting: the averaged dynamics target

p

exactly).

$\sigma=0$
	local slope per segment
$d$	$10{\to}31.6$	$31.6{\to}100$	$100{\to}316$	$316{\to}1000$	top-3 fit
2	$-0.47$	$-0.51$	$-0.50$	$-0.40$	$-0.45$
4	$-0.63$	$-0.64$	$-0.51$	$-0.59$	$-0.55$
8	$-0.53$	$-0.77$	$-0.68$	$-0.54$	$-0.61$
16	$-0.43$	$-0.67$	$-0.89$	$-0.74$	$-0.82$

Appendix

Table 10: Local (segment) slopes of the coupled sup-error between consecutive sweep frequencies

\omega\in\{10,31.6,100,316,1000\}

(segment centers

\approx 17.8,\ 56,\ 178,\ 562

), and the fit over the top three points. Mean over 5 seeds;

N=1.75

	window A ( $\sigma{=}1$ )		window B
$d$	$\Delta_{\min}$	slope	$\Delta_{\min}$	slope, $\sigma{=}0$	slope, $\sigma{=}0.1$	slope, $\sigma{=}1$
2	750	$-0.459\pm 0.021$	10000	$-0.469\pm 0.023$	$-0.440\pm 0.031$	$-0.431\pm 0.005$
4	250	$-0.500\pm 0.013$	3333	$-0.529\pm 0.010$	$-0.497\pm 0.016$	$-0.467\pm 0.019$
8	107	$-0.526\pm 0.012$	1429	$-0.558\pm 0.008$	$-0.557\pm 0.026$	$-0.498\pm 0.009$
16	50	$-0.550\pm 0.007$	667	$-0.615\pm 0.022$	$-0.642\pm 0.004$	$-0.520\pm 0.012$
32	24	$-0.530\pm 0.004$	323	$-0.673\pm 0.013$	$-0.662\pm 0.016$	$-0.523\pm 0.006$

Appendix

Table 11: Whole-window slopes vs. dimension for the two sweep windows, with the top-of-window minimum gap frequency

\Delta_{\min}(\omega_{\mathrm{hi}})=(N{-}1)\,\omega_{\mathrm{hi}}/(d{-}1)

. Window A:

\omega_{\mathrm{hi}}=10^{3}

N=1.75

(

\sigma=1

shown; the

d\leq 32

block coincides with the

c=10

column of Table 9 ). Window B:

\omega_{\mathrm{hi}}=10^{4}

N=2

, all three noise levels. Mean

\pm

std over 5 seeds.

$d$	slope of $t_{\rm conv}$	slope of evals-to-settle	ES/averaged at $\omega{=}10^{3}$
2	$+0.02\,/\,+0.02\,/\,+0.02$	$+1.02\,/\,+1.02\,/\,+1.02$	$1.00\,/\,1.00\,/\,0.99$
4	$+0.02\,/\,+0.02\,/\,+0.07$	$+1.02\,/\,+1.02\,/\,+1.07$	$1.00\,/\,1.00\,/\,1.01$
8	$+0.01\,/\,+0.01\,/\,+0.02$	$+1.01\,/\,+1.01\,/\,+1.02$	$1.00\,/\,1.00\,/\,1.00$
16	$+0.05\,/\,+0.05\,/\,+0.08$	$+1.05\,/\,+1.05\,/\,+1.08$	$1.00\,/\,1.00\,/\,1.00$
32	$+0.04\,/\,+0.04\,/\,+0.04$	$+1.04\,/\,+1.04\,/\,+1.04$	$0.99\,/\,0.99\,/\,1.00$

Appendix

Table 12: Convergence-speed summary. Fitted log–log slopes of

t_{\rm conv}

and of the evaluations-to-settle vs.

\omega

over

\omega\in[10,10^{3}]

, and the ES-to-averaged settling-time ratio at

\omega=10^{3}

. Entries

\sigma=0\,/\,0.1\,/\,1

; 5 seeds; all

75

configurations settled

5/5

$d$	$\omega{=}10$	$\omega{=}31.6$	$\omega{=}100$	$\omega{=}316$	$\omega{=}10^{3}$
2	$2.38{\pm}0.02$ ; $132$	$2.40{\pm}0.01$ ; $423$	$2.54{\pm}0.01$ ; $1412$	$2.58{\pm}0.02$ ; $4551$	$2.60{\pm}0.02$ ; $14478$
4	$2.33{\pm}0.01$ ; $130$	$2.33{\pm}0.00$ ; $411$	$2.48{\pm}0.00$ ; $1384$	$2.55{\pm}0.01$ ; $4485$	$2.56{\pm}0.01$ ; $14261$
8	$2.58{\pm}0.01$ ; $144$	$2.26{\pm}0.02$ ; $399$	$2.43{\pm}0.00$ ; $1355$	$2.51{\pm}0.00$ ; $4425$	$2.54{\pm}0.00$ ; $14143$
16	$2.09{\pm}0.05$ ; $117$	$2.07{\pm}0.07$ ; $364$	$2.36{\pm}0.00$ ; $1312$	$2.48{\pm}0.00$ ; $4368$	$2.52{\pm}0.00$ ; $14049$
32	$2.11{\pm}0.01$ ; $118$	$2.11{\pm}0.03$ ; $372$	$2.39{\pm}0.01$ ; $1334$	$2.42{\pm}0.00$ ; $4258$	$2.50{\pm}0.00$ ; $13943$

Appendix

Table 13: Full settling measurements,

\sigma=0

(

c=20

; 5 seeds; all cells settled

5/5

, transient

5/5

). Per cell:

t_{\rm conv}

and evaluations-to-settle.

$d$	$\omega{=}10$	$\omega{=}31.6$	$\omega{=}100$	$\omega{=}316$	$\omega{=}10^{3}$
2	$2.15{\pm}0.28$ ; $120$	$2.13{\pm}0.17$ ; $375$	$2.22{\pm}0.29$ ; $1237$	$2.05{\pm}0.08$ ; $3604$	$2.43{\pm}0.19$ ; $13513$
4	$1.68{\pm}0.31$ ; $94$	$2.30{\pm}0.17$ ; $405$	$2.41{\pm}0.19$ ; $1342$	$2.32{\pm}0.13$ ; $4095$	$2.47{\pm}0.19$ ; $13734$
8	$0.46{\pm}0.92$ ; $26$	$2.32{\pm}0.13$ ; $409$	$2.29{\pm}0.12$ ; $1277$	$2.46{\pm}0.14$ ; $4326$	$2.46{\pm}0.07$ ; $13709$
16	$2.15{\pm}0.04$ ; $120$	$1.31{\pm}0.11$ ; $232$	$2.39{\pm}0.07$ ; $1333$	$2.49{\pm}0.12$ ; $4380$	$2.48{\pm}0.06$ ; $13835$
32	$2.13{\pm}0.05$ ; $118$	$2.21{\pm}0.15$ ; $389$	$2.17{\pm}0.05$ ; $1209$	$2.49{\pm}0.12$ ; $4378$	$2.51{\pm}0.08$ ; $14001$

Appendix

Table 14: As Table 13 , for

\sigma=1

(the sampling setting). All cells settled

5/5

; transient

5/5

except

(d{=}8,\ \omega{=}10)

with

2/5

— the below-threshold cell excluded from the slope fits of Table 12 .

	$d=2$		$d=10$
Method	NFE $2{\times}10^{3}$	NFE $2{\times}10^{4}$	NFE $2{\times}10^{3}$	NFE $2{\times}10^{4}$
oracle ULA (1 grad)	$0.191\pm 0.038$	$0.215\pm 0.069$	$0.114\pm 0.011$	$0.110\pm 0.018$
oracle MALA (1 grad)	$0.213\pm 0.056$	$0.239\pm 0.053$	$0.113\pm 0.025$	$0.102\pm 0.008$
ZO-ULA (2 ev)	$0.295\pm 0.052$	$0.218\pm 0.041$	$0.123\pm 0.010$	$0.135\pm 0.030$
FD-ULA ( $d{+}1$ ev)	$0.308\pm 0.054$	$0.259\pm 0.095$	$0.307\pm 0.007$	$0.118\pm 0.012$
ES, ours (1 ev)	$0.211\pm 0.026$	$0.232\pm 0.042$	$0.126\pm 0.007$	$0.118\pm 0.022$

Appendix

Table 15: Sliced-

W_{2}

to ground truth on final states of

1024

chains (mean

\pm

std, seeds

1

–

5

). Metric floors:

0.245

(

d{=}2

0.110

(

d{=}10

); entries within

15\%

of the floor are at metric resolution. The ES pilot selected

(T,\omega)=(8,90),(4,1795),(4,180),(16,449)

for the four columns, respectively. All methods recovered

4/4

modes in every run, and unbiased RBF-MMD 2 values were zero within noise throughout.

Method (evals/update)	$d{=}4$	$d{=}16$	$d{=}32$	$d{=}64$	$d{=}128$	$d{=}256$
rand. search, cached ( $1$ )	$100$ (526)	$99.1{\pm}0.3$ (884)	$39.3{\pm}2.6$ (—)	$0$	$0$	$0$
rand. search, re-eval ( $2$ )	$100$ (530)	$100$ (908)	$100$ (1346)	$0$	$0$	$0$
FD-GD ( $2d$ )	$100$ (589)	$100$ (1248)	$0.2{\pm}0.1$ (—) †	$0$	$0$	$0$
ZO-GD, $m{=}10$ ( $20$ )	$100$ (1060)	$100$ (1460)	$100$ (1660)	$79.0{\pm}0.8$ (1920)	$100$ (2440)	$100$ (4784)
SPSA, $\eta_{\mathrm{g}}{=}1/d$ ( $2$ )	$100$ (526)	$100$ (811)	$100$ (967)	$100$ (1246)	$100$ (2111)	$100$ (4303)
bounded ES ( $1$ )	$100$ (608)	$100$ (1036)	$99.8{\pm}0.1$ (1322)	$95.9{\pm}0.7$ (1797)	$99.4{\pm}0.4$ (3084)	$47.7{\pm}1.0$ (—)

Appendix

Table 16: Fixed-target benchmark at

\sigma_{J}=0.03

, EMA readout Eq. ( 86 ): success rate in

\%

(mean

\pm

sd over

5

replicates of

500

chains; sd

0.0

omitted) with the median hitting time

E_{50}

in evaluations in parentheses (seed-to-seed sd of

E_{50}

\leq 18

evaluations in every entry shown); “—” where fewer than half the chains succeed. Budget

B=2000

(

d\leq 64

4000

(

d{=}128

8000

(

d{=}256

). All EMA hitting times include the common filter warm-up measured by the clean-gradient row (which converges in

\leq 40

updates).

\dagger

: warm-up artifact of the filtered readout only — under the last-iterate readout FD-GD succeeds at

d{=}32

(

99.9\%

, final gap

0.267

; Table 18 ) but its

\lfloor B/2d\rfloor=31

updates arrive too late in the budget for the

\tau_{\mathrm{f}}=300

filter to settle; at

d\geq 64

FD-GD fails under both readouts.

Method (evals/update)	$d{=}4$	$d{=}16$	$d{=}32$	$d{=}64$	$d{=}128$	$d{=}256$
rand. search, cached ( $1$ )	$99.8{\pm}0.2$ (540)	$30.5{\pm}1.7$ (—)	$0$	$0$	$0$	$0$
rand. search, re-eval ( $2$ )	$100$ (530)	$100$ (948)	$56.4{\pm}3.5$ (1878)	$0$	$0$	$0$
FD-GD ( $2d$ )	$100$ (586)	$100$ (1248)	$0.1{\pm}0.1$ (—) †	$0$	$0$	$0$
ZO-GD, $m{=}10$ ( $20$ )	$100$ (1060)	$100$ (1464)	$100$ (1672)	$74.0{\pm}0.9$ (1940)	$100$ (2480)	$100$ (5228)
SPSA, $\eta_{\mathrm{g}}{=}1/d$ ( $2$ )	$100$ (532)	$100$ (818)	$100$ (980)	$100$ (1269)	$100$ (2173)	$100$ (4485)
bounded ES ( $1$ )	$100$ (620)	$100$ (1054)	$99.8{\pm}0.2$ (1349)	$90.7{\pm}1.2$ (1850)	$97.4{\pm}1.2$ (3328)	$2.5{\pm}1.0$ (—)

Appendix

Table 17: As Table 16 , at

\sigma_{J}=0.1

(seed-to-seed sd of

E_{50}

\leq 16

evaluations except re-evaluating random search at

d{=}32

\pm 48

Method	$\sigma_{J}{=}0.03$	$\sigma_{J}{=}0.1$	$\sigma_{J}{=}0.3$	$\sigma_{J}{=}1.0$
rand. search, cached	$0.575\to 0.572$	$1.91\to 1.90$	$5.66\to 5.64$	$19.0\to 19.0$
rand. search, re-eval	$0.532\to 0.177$	$1.78\to 0.572$	$5.37\to 1.73$	$18.0\to 5.82$
FD-GD	$0.267\to 0.756^{\dagger}$	$0.282\to 0.766^{\dagger}$	$0.413\to 0.856^{\dagger}$	$1.91\to 1.89$
ZO-GD, $m{=}10$	$0.017\to 0.106^{\dagger}$	$0.037\to 0.117^{\dagger}$	$0.037\to 0.117^{\dagger}$	$0.039\to 0.123^{\dagger}$
SPSA	$0.023\to 0.005$	$0.260\to 0.045$	$0.260\to 0.045$	$0.264\to 0.047$
bounded ES	$0.113\to 0.015$	$0.224\to 0.070$	$1.25\to 0.609$	$25.8\to 21.6$

Appendix

Table 18: Identical-trajectory readout study at

d=32

B=2000

: final median gap

\Delta\mathcal{E}

under the last-iterate (raw) readout

\to

under the EMA readout Eq. ( 86 ), for the same runs (the two scorers observe one trajectory per method per replicate). Mean over

5

replicates of

500

chains; seed-to-seed sd

\leq 0.01

for entries below

1

and

\leq 3\%

of the mean for entries above

1

. The raw

\to

EMA ratio measures each method’s recoverable error — zero-mean scatter about the optimum — as opposed to dynamical failure, which no readout repairs.

\dagger

: EMA exceeds raw due to residual filter warm-up weight

e^{-(B-e_{\mathrm{conv}})/\tau_{\mathrm{f}}}

on pre-convergence states (FD-GD and ZO-GD converge late in their update-sparse budgets), not an error floor.

(a) dimension scan, $\sigma_{J}=0.1$
$d$	clean floor	$F_{\mathrm{osc}}$ pred.	noisy floor	excess (obs.)	excess (pred.)
4	$0.013$	$0.011$	$0.048$	$0.035$	$0.044$
16	$0.042$	$0.045$	$0.124$	$0.082$	$0.088$
32	$0.102$	$0.090$	$0.223$	$0.121$	$0.124$
64	$0.240$	$0.179$	$0.416$	$0.176$	$0.176$
128	$0.456$	$0.358$	$0.712$	$0.256$	$0.249$

Appendix

Table 19: Validation of the floor model Eq. ( 91 ) under the raw (last-iterate) readout. Panel (a): dimension scan at

\sigma_{J}=0.1

(single matched-stream replicate,

2000

chains: the clean and noisy runs share initialization; the gain schedule Eq. ( 89 ) is in force, so the prediction is

0.55\,k(d)\,\sigma_{J}^{2}\,d=2.2\,\sigma_{J}^{2}\sqrt{d}

). Panel (b): noise scan at

d=32

(

5

replicates of

500

chains, mean

\pm

sd; excess relative to the clean floor

0.102

). The

\sigma_{J}{=}1

row lies outside the small-jitter validity of the law (

k\sigma_{J}\approx 0.71

) and is included to locate the breakdown.

Panel of Figure 5 — Figure 5: ES dynamics results for tracking time-varying CIFAR-10 images based on 1D projections.

Figure 6: Additional CELEBA-HQ tracking examples.

Figure 7: Effect on $\rho^{\lambda}$ of $\lambda$ -tempering ES sampling of 2D multimodal target (white dots: visited locations). As $\lambda$ decreases, local peaks and valleys are escaped, and the sampler becomes mass-covering; large $\lambda$ concentrates on modes.

Explore similar work

Energy Generative Modeling: A Lyapunov-based Energy Matching Perspective
May 7, 2026Yixuan Wang, Wenqian Xue, Warren E. DixonEnergy-Based Models
Department of Mechanical and Aerospace Engineering University of Florida
Training-Free Generative Sampling via Moment-Matched Score Smoothing
May 14, 2026Zhenyu Yao, Daniel PaulinLangevin Dynamics Diffusion Sampling
College of Computing and Data Science Nanyang Technological University
Slowly Annealed Langevin Dynamics: Theory and Applications to Training-Free Guided Generation
May 8, 2026Atsushi Nitanda, Dake Bu, Yueming Lyu +1Diffusion Model Guidance Langevin Dynamics
Agency for Science, Technology and Research (A⋆STAR) · Nanyang Technological University · City University of Hong Kong

May 7, 2026J/K move · Enter open · S save

Energy Generative Modeling: A Lyapunov-based Energy Matching Perspective

Yixuan Wang, Wenqian Xue, Warren E. Dixon

Department of Mechanical and Aerospace Engineering University of Florida

Energy-Based Models

Open paper AskPDF arXiv

Generative models based on static scalar energy functions represent an emerging paradigm in which a single time independent potential drives sample generation through its gradient field, eliminating the need for time conditioning entirely. We unify the training and sampling phases of this paradigm, conventionally treated as separate procedures, within a single framework: density transport on the Wasserstein space, cast as a nonlinear control problem in which the Kullback Leibler (KL) divergence serves as a Lyapunov function. Training and sampling are then two instances of this same master dynamics, differing only in initial condition. Within this autonomous framework we develop two analytic results. First, since the Lyapunov certificate is asymptotic, we derive a finite step stopping criterion for Langevin sampling and prove that no Lyapunov certificate exists for the deterministic gradient flow on the same energy landscape. Second, the reformulation brings the toolkit of nonlinear control theory to bear on static scalar energy generative modeling, that is, we show that additive composition of trained scalar energies retains an explicit Gibbs invariant measure and inherits the closed-loop Lyapunov certificate. Beyond these immediate results, this reformulation bridges static scalar energy generative models with the full toolkit of nonlinear control theory, opening the door to barrier functions for constrained generation and contraction metrics for accelerated sampling. Experiments on synthetic distributions validate the theoretical predictions.

Abstract

Explore similar work

Energy Generative Modeling: A Lyapunov-based Energy Matching Perspective

Training-Free Generative Sampling via Moment-Matched Score Smoothing

Slowly Annealed Langevin Dynamics: Theory and Applications to Training-Free Guided Generation