stat.MLOct 7, 2026

Transformed Samplers with Variance Reduction

Authors: Siran Liu, Michalis Tisias, Petros Dellaportas

Organizations: Department of Statistical Science, University College London, UK. · Google DeepMind, UK. · Department of Statistical Science, University College London, UK, and Department of Statistics, Athens University of Economics and Business, Greece.

Read PDF
Ask AI
Cite

Copy citation

DOINot available for this paper

Abstract

Markov chain Monte Carlo (MCMC) methods are the standard tool for computing expectations under complex probability distributions. Control variates reduce the variance of the resulting estimates, but a good control variate requires solving the Poisson equation of the sampler, which rarely admits a closed-form solution. Exact solutions are available when the sampler's kernel has a known spectral decomposition on a simple reference density. In our work, we extend these solutions to general targets through a learned change of variables. A bijection, such as a normalizing flow, is trained so that the target becomes close to the reference in a latent space, and we show that Markov kernels and their Poisson solutions are transformed by any bijection. Running such samplers in the latent space then yields explicit control variates, and the estimator is consistent under mild tail conditions on the map and target. Importance sampling (IS) from the flow is the limiting case of the same construction and the control variates apply to it as well. Experiments on synthetic targets and real posteriors compare the procedure against state-of-the-art samplers and control variates.

Figures & tables

Figure 1: When to sample instead of reweighting. Each point is one target of Table 1 . The ratio is the precision per second of the best chain-based estimator on the main test function divided by that of SNIS+CV, and the marker shows which chain is best; above the dotted line a chain is more efficient than SNIS+CV. The horizontal axis is the largest normalized importance weight of the flow times NN , over the 100100 replicates of N=106N=10^{6} draws.
TargetSNISTGI-RWMTGI-MALATESSNeuTra
GARCH9.7e7 † (29)8.4e7 † (37)2.2e7 (64)2.4e7 (10)5.0e7 (109)
One-comp.8.9e5 (1.6)6.7e5 (1.7)4.6e5 (12)4.7e5 (1.6)5.8e5 (5.9)
Lotka–Volterra1.0e8 (28)2.3e7 (14)5.6e6 (57)8.9e6 (6.4)2.1e7 (176)
Banana4.1e7 (26)3.4e7 (35)1.1e7 (117)5.9e6 (13)1.7e7 (104)
Dugongs2.0e7 (5.1)5.4e6 (3.5)5.1e4 (0.6)9.2e6 (4.9)1.3e6 (3.3)
Mixture3.9e6 (3.3)1.5e6 (2.3)–1.6e6 (2.4)1.6e4 (1.9)
Table 1: Precision per second on the main test function of each target (Tables S8 and S9 ), with the protocol of Appendix S4 . In parentheses, the variance reduction factor of the control variates, the mean squared error without them over that with them on the same run. Targets are ordered by the largest normalized importance weight of the flow times NN . Entries within 10%10\% of the best in a row are in bold. †\dagger : the Monte Carlo error of the reference value accounts for at least half of the mean squared error; –: RMSE above one posterior standard deviation.
Appendix figures & tables16 assets

Supplementary material from the paper’s appendix.

Appendix

SamplerMeanCovarianceDerivativesCost per step
TGI-RWMβz\beta z(1−β2)I(1-\beta^{2})Inone11
TGI-MALAβz+γ∇log⁡w(z)\beta z+\gamma\nabla\log w(z)(1−β2)I(1-\beta^{2})Igradient22
TGI-pMALAz+γAz∇log⁡πz(z)z+\gamma A_{z}\nabla\log\pi_{z}(z)(1−β2)Az(1-\beta^{2})A_{z}diagonal Hessian2+2d2+2d
Appendix
Table S1: The TGI samplers on the Gaussian reference. The proposal is N(mean,cov)\mathcal{N}(\text{mean},\text{cov}) in the latent space. Cost counts one evaluation of log⁡π\log\pi as one unit and a joint evaluation with its gradient as two.
 KK1246810
γ=0.3\gamma=0.3Predicted15.62701.28×1051.28\times 10^{5}1.03×1081.03\times 10^{8}1.30×10111.30\times 10^{11}2.39×10142.39\times 10^{14}
 Empirical15.12681.29×1051.29\times 10^{5}1.17×1081.17\times 10^{8}1.77×10111.77\times 10^{11}4.58×10144.58\times 10^{14}
γ=1\gamma=1Predicted8.351023.35×1043.35\times 10^{4}2.27×1072.27\times 10^{7}2.63×10102.63\times 10^{10}4.66×10134.66\times 10^{13}
 Empirical8.1097.53.24×1043.24\times 10^{4}2.35×1072.35\times 10^{7}3.03×10103.03\times 10^{10}8.13×10138.13\times 10^{13}
Appendix
Table S2: VRF of TGI-RWM on LogNormal(0,0.52)\mathrm{LogNormal}(0,0.5^{2}) with the exact map, F(x)=xF(x)=x and θ=(1,−1)\theta=(1,-1) . Predicted: Proposition 5.1 . Empirical: ratio of mean squared errors over 10001000 stationary chains of 50005000 steps. At γ=1\gamma=1 the chain draws independently from ρ\rho .
kk123456
λk=E[cos⁡kϑ]\lambda_{k}=\mathbb{E}[\cos^{k}\vartheta]00.500.37500.3125
Autocorrelation of Hk(Zt)H_{k}(Z_{t})0.0010.499−0.002-0.0020.368−0.009-0.0090.274
KK1246810
VRF, Proposition 5.13.581173.94×1043.94\times 10^{4}2.71×1072.71\times 10^{7}3.17×10103.17\times 10^{10}5.62×10135.62\times 10^{13}
VRF, simulation3.621093.81×1043.81\times 10^{4}2.91×1072.91\times 10^{7}4.69×10104.69\times 10^{10}1.09×10141.09\times 10^{14}
Appendix
Table S3: Elliptical slice sampling on LogNormal(0,0.52)\mathrm{LogNormal}(0,0.5^{2}) with the exact map, 10001000 stationary chains of 50005000 steps. Top: eigenvalues of the candidate kernel on HkH_{k} and the lag-one autocorrelation of Hk(Zt)H_{k}(Z_{t}) . Bottom: VRF of F(x)=xF(x)=x with the dictionary {H1,…,HK}\{H_{1},\dots,H_{K}\} and θ=(1,−1)\theta=(1,-1) , from Proposition 5.1 and from the ratio of mean squared errors.
ReferenceMapβ\betaAcceptanceESS /n/nVRFESS ×\times VRF /n/nTail gap
Gaussianmoment-matched0.50.780.111.90.2014.4
Gammamoment-matched0.10.920.60169.70.37
Gaussianexact01.000.979.3×1039.3\times 10^{3}9.0×1039.0\times 10^{3}0
Gammaexact0.21.000.601.1×1041.1\times 10^{4}6.9×1036.9\times 10^{3}0
Appendix
Table S4: Gaussian and Gamma references on 12Gamma(2,1)+12Gamma(5,1)\tfrac{1}{2}\mathrm{Gamma}(2,1)+\tfrac{1}{2}\mathrm{Gamma}(5,1) , F(x)=xF(x)=x , K=6K=6 , over 10001000 chains of 50005000 steps started from exact draws of πz\pi_{z} . For each row, β∈[0,0.95]\beta\in[0,0.95] is chosen on a grid to maximize ESS ×\times VRF. The last column is the largest value of log⁡w\log w minus its 99.9%99.9\% quantile on 2×1052\times 10^{5} target draws; for the Gaussian moment-matched map log⁡w\log w is unbounded and this gap grows with the number of draws.
    ESS /n/nVRF
TargetMapStepAcceptanceobservedpredictedobservedpredicted
Circleexactτ=0.5\tau=0.51.000.150.1519001810
Circleexactτ=2\tau=21.000.520.53502518
Circleconcentrations ×0.5\times 0.5τ=2\tau=20.800.31–5.4–
Circleconcentrations ×2\times 2τ=2\tau=20.770.39–2.8–
Sphereexactφ=0.6\varphi=0.61.000.110.12783678
Appendix
Table S5: The MH sampler proposing from QQ on the circle ( 0.7 vM(0.5,4)+0.3 vM(3.5,2)0.7\,\mathrm{vM}(0.5,4)+0.3\,\mathrm{vM}(3.5,2) , F=cos⁡θF=\cos\theta , Fourier degree 66 ) and on the sphere ( vMF\mathrm{vMF} with κ=3\kappa=3 , F=x⋅μF=x\cdot\mu , harmonics of degree at most 44 ), with 10001000 chains of 50005000 steps. Each map is the quantile map of a model π′\pi^{\prime} , so that T#ρ=π′T_{\#}\rho=\pi^{\prime} and w=(π/π′)∘Tw=(\pi/\pi^{\prime})\circ T : on the circle, π′\pi^{\prime} is the target mixture with both concentrations scaled; on the sphere, π′=vMF(κ′)\pi^{\prime}=\mathrm{vMF}(\kappa^{\prime}) . Predicted values are from Proposition 5.1 for the exact maps.
Varianttrain costESS IS/N{}_{\rm IS}/NF\mathcal{F}tail log⁡w\log wVRFmax ∣z∣|z|bias
B5.1×1065.1\times 10^{6}0.9970.02877.634.19 / 2.92208−1.0×10−4-1.0\times 10^{-4}
B + tail7.4×1067.4\times 10^{6}0.9930.0361133.2 / 29.44+5.0×10−5+5.0\times 10^{-5}
A5.5×1065.5\times 10^{6}0.9970.01715.468.95 / 14.596−9.9×10−5-9.9\times 10^{-5}
A + tail (default)7.8×1067.8\times 10^{6}0.9980.01470.69951.5 / 48.25−4.6×10−6-4.6\times 10^{-6}
A + tail, SMALL7.8×1067.8\times 10^{6}0.9960.02450.64840.5 / 39.89−6.5×10−5-6.5\times 10^{-5}
A + tail, 1/10 cost7.8×1057.8\times 10^{5}0.9500.3291.885.66 / 5.768+3.4×10−5+3.4\times 10^{-5}
Appendix
Table S6: Training ablation, banana: every variant trained once for 10410^{4} steps, except the 1/10-cost one, a complete training of 10310^{3} steps with every schedule scaled to its length; flow diagnostics at the last checkpoint (tail log⁡w\log w : 99.9% quantile over the radial band relative to the median over N(0,I)\mathcal{N}(0,I) ) and TGI-RWM chains at JJ ’s γ\gamma (two seeds of 2000 chains, 500 + 5000 steps; main test function). max ∣z∣|z| : the largest ∣μ^r−median μ^∣/(1.4826 MAD)|\hat{\mu}_{r}-\mathrm{median}\,\hat{\mu}|/(1.4826\,\mathrm{MAD}) over the chains’ control-variate estimates μ^r\hat{\mu}_{r} (about 4 for Gaussian spread; larger values mean that a few chains stuck). bias: (μˉ−μ)/σπ(\bar{\mu}-\mu)/\sigma_{\pi} of the pooled estimate; ∗ : more than four standard errors sd(μ^r)/R\mathrm{sd}(\hat{\mu}_{r})/\sqrt{R} .
Varianttrain costESS IS/N{}_{\rm IS}/NF\mathcal{F}tail log⁡w\log wq(x1>0)q(x_{1}>0)VRFmax ∣z∣|z|bias
A7.8×1067.8\times 10^{6}0.9861011.070.9623.5 / 23.67+0.871 ∗
A a without anchor8.2×1068.2\times 10^{6}0.9203881.70.563.59 / 3.118+6.5×10−5+6.5\times 10^{-5}
A a (default)8.2×1068.2\times 10^{6}0.9596.3×1036.3\times 10^{3}1.630.523.59 / 3.8210−5.1×10−4-5.1\times 10^{-4}
A a without tail5.9×1065.9\times 10^{6}0.9325.6×1035.6\times 10^{3}8.370.513.2 / 1.4297+1.6×10−4+1.6\times 10^{-4}
A a , 1/10 cost8.6×1058.6\times 10^{5}0.4891964.620.641.37 / 1.2638+3.8×10−3+3.8\times 10^{-3} ∗
B7.4×1067.4\times 10^{6}0.98585.31.120.969.81 / 8.5216+0.871 ∗
Appendix
Table S7: Training ablation, Gaussian mixture: every variant trained once for 10410^{4} steps, except the 1/10-cost one, a complete training of 10310^{3} steps with every schedule scaled to its length (pilot 8×1048\times 10^{4} ); flow diagnostics at the last checkpoint (tail log⁡w\log w : 99.9% quantile over the radial band relative to the median over N(0,I)\mathcal{N}(0,I) ) and TGI-RWM chains at JJ ’s γ\gamma (two seeds of 2000 chains, 500 + 5000 steps; main test function). max ∣z∣|z| : the largest ∣μ^r−median μ^∣/(1.4826 MAD)|\hat{\mu}_{r}-\mathrm{median}\,\hat{\mu}|/(1.4826\,\mathrm{MAD}) over the chains’ control-variate estimates μ^r\hat{\mu}_{r} (about 4 for Gaussian spread; larger values mean that a few chains stuck). bias: (μˉ−μ)/σπ(\bar{\mu}-\mu)/\sigma_{\pi} of the pooled estimate; ∗ : more than four standard errors sd(μ^r)/R\mathrm{sd}(\hat{\mu}_{r})/\sqrt{R} .
TargetddSchemeTrain. costmax⁡w⋅N\max w{\cdot}NMain FFLLγ\gamma
Banana8A\mathrm{A}7.8×1067.8\times 10^{6}6565x2x_{2}10.90 / 0.20
Funnel10A\mathrm{A}7.8×1067.8\times 10^{6}1.6×1031.6\times 10^{3}vv21.00 / 0.10
Gaussian mixture10Aa\mathrm{A}^{\mathrm{a}}8.2×1068.2\times 10^{6}1.2×1031.2\times 10^{3}x1x_{1}20.95 / 0.20
tt mixture10Aa\mathrm{A}^{\mathrm{a}}8.2×1068.2\times 10^{6}1.3×1051.3\times 10^{5}x1x_{1}21.00 / 0.10
Appendix
Table S8: Settings for the synthetic targets. All flows use the LARGE network, 10410^{4} training steps, batch size 256256 and the tail penalty; the multi-start pilot of the mixtures costs 4×1054\times 10^{5} and is included in the training cost. max⁡w⋅N\max w\cdot N is the largest normalized importance weight of the flow times N=106N=10^{6} , the largest over 100100 replicates. LL is the number of leapfrog steps of NeuTra and γ\gamma the step size of TGI-RWM / TGI-MALA, both chosen by pilot runs.
Figure S1: Analytic transports of the banana, main test function x2x_{2} , from the chain experiments of Appendix S4 . From left to right: the best ESS per unit cost over the step size of TGI-RWM and TGI-MALA, and ESSIS/N\mathrm{ESS}_{IS}/N , against the Fisher divergence of the shear map TcT_{c} ; the ESS, the VRF and their product against the step size γ\gamma , relative to their maxima, for TGI-RWM and TGI-MALA on the accurate member c=0.01c=0.01 (dashed: the step chosen by JJ ; dotted: by the speed measure); and the VRF (solid) and the number of kept Hermite terms (dotted) against the maximal degree KK for each pruning rule, with TGI-MALA on the same member.
Figure S2: The banana with the training of the flow charged. Left: RMSE of the main test function in posterior standard deviations against the total number of target evaluations, with bootstrap 95%95\% bands; each flow-based method shows the lower envelope over training checkpoints and sampling budgets. Right: a grid of the latent space pushed forward by the flow after 100100 , 10001000 and 10410^{4} training steps, over the target density (grey; contours hold 5050 , 9090 and 99%99\% of the mass).
PosteriorddSchemeNetworkTailTrain. costmax⁡w⋅N\max w{\cdot}NMain FFLLγ\gamma
GARCH(1,1)4A\mathrm{A}SMALLyes7.8×1067.8\times 10^{6}33β1\beta_{1}10.90
One-compartment PK4A\mathrm{A}SMALLyes7.8×1067.8\times 10^{6}99log⁡Km\log K_{m}10.99
Dugongs4B\mathrm{B}SMALLyes7.4×1067.4\times 10^{6}3.5×1023.5\times 10^{2}β\beta10.80
Lotka–Volterra8A\mathrm{A}SMALLyes7.8×1067.8\times 10^{6}4242α\alpha11.00
GP Poisson13A\mathrm{A}LARGEno5.5×1065.5\times 10^{6}7.7×1037.7\times 10^{3}α\alpha20.99
Eight schools10A\mathrm{A}LARGEno1.6×1071.6\times 10^{7}6.8×1036.8\times 10^{3}θ1\theta_{1}10.90
Appendix
Table S9: Settings for the real posteriors: dimension, training scheme, network, tail penalty, training cost in target evaluations, importance sampling effective sample size of the flow, main test function, leapfrog steps LL selected for NeuTra and step size γ\gamma of TGI-RWM, both chosen by pilot runs. All flows are trained for 10410^{4} steps with batch size 256256 , except eight schools ( 3×1043\times 10^{4} steps) and the sparse logistic regression ( 3×1053\times 10^{5} steps, batch size 10241024 , a LARGE network of width 256256 , marked ∗ ); for the ϕ4\phi^{4} field the multi-start pilot uses 6464 NUTS chains and 4×1064\times 10^{6} evaluations, included in the training cost.
Figure S3: The 100100 estimates of the main test function of every target, in posterior standard deviations from the reference value (dashed line). Panels are ordered by max⁡w⋅N\max w\cdot N ; circles are estimates beyond 1.51.5 interquartile ranges from the box.
PosteriorFFNUTS+ZVSNIS+CVTGI-RWM+CVTESS+CVNeuTra+ZV
GARCHβ1\beta_{1}2.1e-31.9e-42.0e-43.5e-42.9e-4
One-compartmentlog⁡Km\log K_{m}1.3e-21.3e-31.3e-31.5e-31.3e-3
Lotka–Volterraα\alpha1.6e-31.9e-43.1e-44.7e-42.7e-4
Dugongsβ\beta5.0e-35.0e-48.5e-45.8e-42.0e-3
GP Poissonα\alpha1.5e-21.8e-31.1e-31.0e-37.5e-4
Eight schoolsθ1\theta_{1}1.4e-21.6e-31.4e-31.2e-32.1e-3
Appendix
Table S10: RMSE in posterior standard deviations of the main test function at B=106B=10^{6} target evaluations, without the training of the flow, over 100100 estimates for the flow-based estimators and 5050 for NUTS. NUTS+ZV was more accurate than plain NUTS on every posterior where both were run. †\dagger : extrapolated from the reference run as its Monte Carlo standard error times Cref/B\sqrt{C_{\rm ref}/B} ; ‡\ddagger : NUTS stays in one mode. Entries within 10%10\% of the best in a row are in bold.
  SNISTGI-RWMTGI-RWM δTGI-MALATESSNeuTra
TargetFFplain+CVplain+CVplain+CVplain+CVplain+CVplain+ZV
GARCHβ1\beta_{1} ⋆\star3.3e69.7e7 †2.3e68.4e7 †1.6e66.3e7 †3.4e52.2e72.3e62.4e74.6e55.0e7
 μ\mu3.3e65.3e82.3e62.6e81.6e61.8e83.4e53.6e72.5e63.6e74.4e55.7e8
 α0\alpha_{0}3.1e63.2e7 †2.2e63.2e7 †1.5e62.6e7 †3.5e51.4e72.1e61.7e74.7e53.3e7
 α1\alpha_{1}3.4e64.8e8 †2.3e62.4e8 †1.6e61.5e8 †3.4e52.8e72.3e63.2e74.6e58.1e7
One-comp.log⁡Km\log K_{m} ⋆\star5.7e58.9e54.0e56.7e53.9e56.4e54.0e44.6e52.9e54.7e59.9e45.8e5
Appendix
Table S11: Precision per second on every test function, extending Table 1 ; ⋆\star marks the main test function. For the sparse logistic regression we give log⁡τ\log\tau , w1w_{1} and the median, maximum and minimum over the 2525 coefficients wjw_{j} . Entries within 10%10\% of the best in a row are in bold. †\dagger : the Monte Carlo error of the reference value accounts for at least half of the mean squared error; ‡\ddagger : the errors of this target are limited by a bias of about 0.010.01 posterior standard deviations in log⁡τ\log\tau shared by the flow-based estimators; –: RMSE above one posterior standard deviation.
Methodmean gapfactorsdworst rowwithin 2×2\timesbestfailed
(A) main test function, 1212 functions
SNIS−0.36-0.362.30.65−1.74-1.74 ( tt mixture, x1x_{1} )9/128/120
TGI-RWM−0.28-0.281.90.20−0.65-0.65 (Lotka–Volterra, α\alpha )6/121/120
TGI-MALA−1.98-1.98962.49−6.46-6.46 (Mixture, x1x_{1} )3/122/123
TESS−0.39-0.392.40.31−1.06-1.06 (Lotka–Volterra, α\alpha )5/121/120
NeuTra−0.83-0.836.80.87−2.40-2.40 (Mixture, x1x_{1} )5/122/120
Appendix
Table S12: Stability of the estimators: the gap log⁡10(precision/best of the row)\log_{10}(\text{precision}/\text{best of the row}) of each method with its control variates to the most precise one, over all test functions of Table S11 ( 00 : the best; −1-1 : ten times less precise). Mean gap, the corresponding factor 10−mean10^{-\text{mean}} , standard deviation and worst row; number of rows within a factor 22 of the best and within 10%10\% of it (the bold rule); failures (RMSE above one posterior standard deviation) are counted at RMSE =1=1 . (A) The main test function of each target. (B) All 70 test functions; the mean and standard deviation are over the 1212 per-target means, so that every target has the same weight, the worst row and the counts over all rows.
   TGI-RWMTGI-MALA
TargetMemberF\mathcal{F}γJ\gamma_{J}γ⋆\gamma^{\star}acc. γJ\gamma_{J} / 11R(1,γJ)\mathrm{R}(1,\gamma_{J})R(γJ,γ⋆)\mathrm{R}(\gamma_{J},\gamma^{\star})γJ\gamma_{J}γ⋆\gamma^{\star}acc. γJ\gamma_{J} / 11R(1,γJ)\mathrm{R}(1,\gamma_{J})R(γJ,γ⋆)\mathrm{R}(\gamma_{J},\gamma^{\star})
Bananac=0.01c=0.017×10−47\times 10^{-4}0.951.000.99/0.991.010.990.570.501.00/0.990.600.99
 c=0.3c=0.30.630.630.910.600.81/0.810.880.870.370.350.92/0.690.011.00
 c=1c=17.07.00.480.450.60/0.550.620.820.140.070.82/0.280.120.40
Funnelϱ=0.01\varrho=0.011.3×10−31.3\times 10^{-3}0.990.850.99/0.991.001.000.400.300.99/0.980.220.90
 ϱ=0.3\varrho=0.31.21.20.990.100.70/0.701.000.440.220.300.94/0.450.080.99
Appendix
Table S13: The step γJ\gamma_{J} selected by JJ against γ=1\gamma=1 and against the best step γ⋆\gamma^{\star} on a grid of 2525 values, on the analytic families: acceptance rate at γJ\gamma_{J} and at 11 , and the ratios R(γ,γ′)\mathrm{R}(\gamma,\gamma^{\prime}) of ESS ×\times VRF of the main test function at γ\gamma to that at γ′\gamma^{\prime} , from 500500 chains of 50005000 steps with common random numbers.

Explore similar work

CardsList
Jun 26, 2025stat.ML

Gaussian Invariant Markov Chain Monte Carlo

We develop sampling methods, which consist of Gaussian invariant versions of random walk Metropolis (RWM), Metropolis adjusted Langevin algorithm (MALA) and second order Hessian or Manifold MALA. Unlike standard RWM and MALA, we show that Gaussian invariant sampling can lead to ergodic estimators with improved statistical efficiency. This is due to a remarkable property of Gaussian invariance that allows us to obtain exact analytical solutions to the Poisson equation for Gaussian targets. These solutions can be used to construct efficient and easy to use control variates for variance reduction of estimators under any intractable target. We demonstrate the new samplers and estimators in several examples, including high dimensional targets in latent Gaussian models where we compare against several advanced methods and obtain state-of-the-art results. We also provide theoretical results regarding geometric ergodicity, and an optimal scaling analysis that shows the dependence of the optimal acceptance rate on the Gaussianity of the target.

Michalis K. Titsias, Angelos Alexopoulos, Siran Liu +1

Google DeepMind, UK · Athens University of Economics and Business, Greece · UCL, UK

Markov Chain Monte CarloVariance Reduction
AskPDF
Jun 25, 2024math.ST

Variance Reduction for Independent Metropolis

Assume that we would like to estimate the expected value of a function FF with respect to an intractable density ππ, which is specified up to some unknown normalising constant. We prove that if ππ is close enough under KL divergence to another density qq, an independent Metropolis sampler estimator that obtains samples from ππ with proposal density qq, enriched with a variance reduction computational strategy based on control variates, achieves smaller asymptotic variance than i.i.d. sampling from ππ. The control variates construction requires no extra computational effort but assumes that the expected value of FF under qq is analytically available. We illustrate this result by calculating the marginal likelihood in a linear regression model with prior-likelihood conflict and a non-conjugate prior. Furthermore, we propose an adaptive independent Metropolis algorithm that adapts the proposal density such that its KL divergence with the target is being reduced. We demonstrate its applicability in a Bayesian logistic and Gaussian process regression problems and we rigorously justify our asymptotic arguments under easily verifiable and essentially minimal conditions.

Siran Liu, Petros Dellaportas, Michalis K. Titsias

Department of Statistical Science, University College London, UK. · Department of Statistics, Athens University of Economics and Business, Greece. · Google DeepMind, UK.

Markov Chain Monte CarloVariance Reduction
AskPDF
Dec 23, 2025cs.LG

Control Variate Score Matching for Diffusion Models

Sampling from unnormalized probability densities is a pervasive challenge across the computational and physical sciences. Diffusion models provide a powerful generative framework for this task, but their success relies on accurately estimating the score of the perturbed target distribution. Current approaches face a dichotomy between two standard estimation methods: the Denoising Score Identity (DSI) requires data samples and exhibits high variance at low noise levels, whereas the Target Score Identity (TSI) relies on the energy function and suffers from diverging variance at high noise levels. In this work, we reconcile both approaches by introducing the Control Variate Score Identity (CVSI), an unbiased estimator with an analytically optimal, state- and time-dependent control coefficient that theoretically minimizes variance over the entire diffusion process. CVSI serves as a robust plug-in estimator that significantly enhances performance and efficiency in data-free sampler learning and training-free diffusion sampling. These gains scale to complex, high-dimensional energy-based models.

Khaled Kahouli, Romuald Elie, Klaus-Robert Müller +3

1Google DeepMind · 2BIFOLD – Berlin Institute for the Foundations of Learning and Data · 3Machine Learning Group, Technische Universität Berlin +2

Score MatchingDiffusion Model SamplingScore-Based Generative ModelingVariance Reduction
AskPDF