Scalable Maximum Entropy Population Synthesis via Persistent Contrastive Divergence

Mirko Degli Esposti

cs.LGMar 28, 2026

Scalable Maximum Entropy Population Synthesis via Persistent Contrastive Divergence

Authors: Mirko Degli Esposti

Organizations: Department of Physics and Astronomy, City Science Lab University of Bologna, Bologna, Italy

Abstract

Maximum entropy (MaxEnt) modelling provides a principled framework for generating synthetic populations from aggregate census data, without access to individual-level microdata. The bottleneck of exact-enumeration approaches is expectation computation by explicit summation over the full tuple space $\cX$ , which becomes infeasible for more than $K \approx 20$ categorical attributes; sampling-based alternatives exist but rely on Metropolis-type schemes that require proposal tuning and rejection steps. We propose GibbsPCDSolver, a stochastic replacement for this computation based on Persistent Contrastive Divergence (PCD): a persistent pool of $N$ synthetic individuals is updated by Gibbs sweeps at each gradient step, providing a stochastic approximation of the model expectations without ever materialising $\cX$ . We validate the approach on controlled benchmarks and on Syn-ISTAT, a $K{=}15$ Italian demographic benchmark with analytically exact marginal targets derived from ISTAT-inspired conditional probability tables. Scaling experiments across $K \in \{12, 20, 30, 40, 50\}$ confirm that GibbsPCDSolver maintains $\MRE \in [0.010, 0.018]$ while $|\cX|$ grows eighteen orders of magnitude, with runtime scaling as $O(K)$ rather than $O(|\cX|)$ . On Syn-ISTAT, GibbsPCDSolver reaches $\MRE{=}0.03$ on training constraints and -- crucially -- produces populations with effective sample size $\Neff = N$ versus $\Neff \approx 0.012\,N$ for generalised raking, an $86.8{\times}$ diversity advantage that is essential for agent-based urban simulations.

Figures & tables

Parameter	Recommendation	Rationale
Pool size $N$	$25\,000$ – $100\,000$	MRE floor $\approx 1/(2\sqrt{N})$ ; sweet spot at $N{=}25$ K
Sweeps $s$	$1$ (sparse, arity $\leq 2$ ); $5$ (mixed arity $\leq 3$ )	Extra sweeps reduce mixing bias (not gradient variance)
Learning rate $\eta$	$0.01$ for $K\geq 15$	$\eta{=}0.05$ oscillates at $K{=}15$ with arity-3 constraints
Window $W$	$50$	Longer window avoids premature stopping during slow descent
Initialisation	$\bm{\lambda}=\mathbf{0}$ , pool uniform	Max-entropy prior; data-based initialisation may introduce CD bias

Table 1: Practical hyperparameter guidelines for GibbsPCDSolver.

$N$	MRE	$\\|\Delta\bm{\lambda}\\|/\\|\bm{\lambda}^{*}\\|$	KL	Time (s)
$500$	0.101	0.344	$5.9{\times}10^{-3}$	0.6
$2{,}000$	0.057	0.341	$1.9{\times}10^{-3}$	1.2
$10{,}000$	0.021	0.349	$4.8{\times}10^{-4}$	6.1
$100{,}000$	0.010	0.347	$5.6{\times}10^{-5}$	86.7

Table 2: Experiment A0 (

K{=}6

|\mathcal{X}|{=}216

m{=}33

s{=}3

\eta{=}0.05

). Exact L-BFGS:

\mathrm{MRE}{=}3{\times}10^{-5}

t{=}0.06

\|\Delta\bm{\lambda}\|

: relative parameter distance

\|\bm{\lambda}_{\mathrm{MCMC}}{-}\bm{\lambda}^{*}\|/\|\bm{\lambda}^{*}\|

Figure 1: Experiment A0 ( $K{=}6$ , $|\mathcal{X}|{=}216$ ). Left : MRE convergence curves for four pool sizes (log scale); dashed line is exact L-BFGS reference. Centre : scatter of $\hat{\bm{\lambda}}$ vs. $\bm{\lambda}^{*}$ at $N{=}10{,}000$ ; the systematic offset from the diagonal reflects a gauge degree of freedom in $\bm{\lambda}$ (see footnote), not solver error. Right : estimated frequencies $\hat{\alpha}_{j}$ vs. targets $\alpha_{j}$ at $N{=}10{,}000$ ; points lie close to the diagonal across the full range $[0.004,\,0.906]$ .

Figure 2: Experiment A1a ( $K{=}8$ , $|\mathcal{X}|{=}6{,}144$ , $m{=}61$ , $s{=}5$ ). Left : MRE convergence curves for four pool sizes (log scale); dashed line is the exact L-BFGS reference ( $\mathrm{MRE}{=}4{\times}10^{-5}$ ). Centre : scatter of $\bm{\lambda}_{\mathrm{MCMC}}$ vs. $\bm{\lambda}^{*}$ at $N{=}50{,}000$ ; the constant offset from the diagonal is the gauge degree of freedom induced by unary constraints and does not affect the distribution $p_{\bm{\lambda}}$ . Right : final MRE vs. $N$ (log–log); the dashed reference line follows $1/\sqrt{N}$ , confirming the theoretical variance floor.

Figure 3: Experiment A1b ( $K{=}10$ , $|\mathcal{X}|{=}59{,}049$ , $m{=}30$ binary constraints, $s{=}5$ ). Left : KL divergence $\mathrm{KL}(p_{\bm{\lambda}^{*}}\|p_{\bm{\lambda}_{\mathrm{MCMC}}})$ (red, left axis) and relative parameter distance $\|\Delta\bm{\lambda}\|/\|\bm{\lambda}^{*}\|$ (blue, right axis) vs. pool size $N$ (log–log); the dotted line shows the $O(1/N)$ reference. Centre : scatter of $\bm{\lambda}_{\mathrm{MCMC}}$ vs. $\bm{\lambda}^{*}$ at $N{=}50{,}000$ ; unlike A0–A1a, no gauge offset is present and points cluster tightly around the diagonal ( $\|\Delta\bm{\lambda}\|/\|\bm{\lambda}^{*}\|{=}0.016$ ). Right : estimated constraint frequencies $\hat{\alpha}_{j}$ vs. exact targets $\alpha^{*}_{j}$ at $N{=}50{,}000$ , for binary constraints (blue).

Exp.	Method	$N$	$s$	MRE	$\\|\Delta\bm{\lambda}\\|/\\|\bm{\lambda}^{*}\\|$	KL
A1a	Exact L-BFGS	—	—	$4\times 10^{-5}$	—	—
(K=8,	Gibbs	1,000	5	0.0991	0.498	0.00898
$\|\mathcal{X}\|$	Gibbs	5,000	5	0.0373	0.505	0.00098
$=$ 6144)	Gibbs	20,000	5	0.0220	0.503	0.00029
	Gibbs	50,000	5	0.0111	0.506	0.00018
	Gibbs	20,000	10	0.0231	0.505	0.00049

Table 3: Algorithmic validation: GibbsPCDSolver vs. exact L-BFGS. MRE: mean relative error on training constraints.

\|\Delta\bm{\lambda}\|

: relative parameter distance (computed only for A1b where

\bm{\lambda}^{*}

is known exactly). KL:

\mathrm{KL}(p_{\bm{\lambda}^{*}}\|p_{\bm{\lambda}_{\mathrm{MCMC}}})

(computed only when

|\mathcal{X}|

is enumerable).

$K$	Method	$\|\mathcal{X}\|$	MRE	$t$ (s)	$N_{\mathrm{eff}}$	$N_{\mathrm{eff}}/N$	$H$ (nats)
12	Gibbs PCD	$8.3{\times}10^{4}$	0.012	607	100,000	100%	9.40
12	Raking	—	$0^{\dagger}$	7	7,138	7.1%	9.19
20	Gibbs PCD	$2.7{\times}10^{7}$	0.010	440	100,000	100%	11.43
20	Raking	—	$0^{\dagger}$	16	3,102	3.1%	9.37
30	Gibbs PCD	$4.7{\times}10^{11}$	0.010	965	100,000	100%	11.51
30	Raking	—	$0^{\dagger}$	16	790	0.8%	8.36

Table 4: A2 scaling experiments: GibbsPCDSolver vs. raking across

K\in\{12,20,30,40,50\}

with ternary constraints only,

N{=}100{,}000

s{=}5

. For

K{=}12

the exact L-BFGS reference is available (

\mathrm{MRE}{=}1.3{\times}10^{-4}

t{=}2.4

s); for

K\geq 20

|\mathcal{X}|

is not enumerable. Raking

\mathrm{MRE}{=}0

for

K\leq 30

is exact on training constraints by algebraic construction and is not a meaningful accuracy metric; for

K\geq 40

raking fails to satisfy the constraints (

\mathrm{MRE}>0

Figure 4: A2 scaling experiments ( $N{=}100{,}000$ , arity=3). Left : MRE of GibbsPCDSolver (solid blue) vs. exact MaxEnt of Pachet and Zucker (2026) (dashed blue, their Table 2). The two curves use different problem instances (WuGenerator vs. NPORS); the comparison is qualitative—both achieve MRE in the $1$ – $2\%$ range—and is not a claim that GibbsPCDSolver outperforms exact enumeration. At $K{=}30$ only GibbsPCDSolver is available ( $|\mathcal{X}|\approx 5{\times}10^{11}$ ). Centre : raking $N_{\mathrm{eff}}/N$ (red diamonds, log scale, left axis) and Shannon entropy $H$ (dotted lines, right axis) for both methods. Right : fit time vs. $K$ (log scale).

T1: edu $\times$ emp $\to$ income; T2: area $\times$ car $\to$ transport; T3: sex $\times$ age $\to$ employment.
Method	$\mathrm{MRE}_{\text{train}}$	$\mathrm{MRE}_{\text{T1}}$	$\mathrm{MRE}_{\text{T2}}$	$\mathrm{MRE}_{\text{T3}}$
Gibbs PCD	0.018	0.349	0.169	0.197
Raking	0.000	0.347	0.139	0.214

Table 5: A-ISTAT-2: held-out MRE on ternary constraints. Both solvers trained on 28 unary+binary constraints (

N{=}100{,}000

). Bold : better result.

Method	$N_{\mathrm{eff}}$	$N_{\mathrm{eff}}/N$	$H$ (nats)	Unique profiles
Gibbs PCD	$100{,}000$	$100\%$	$11.47$	$97{,}025$
Raking	$1{,}152$	$1.2\%$	$8.32$	$25{,}786$
Ratio (G/R)	$86.8{\times}$	—	$+3.15$	$3.8{\times}$

Table 6: Population diversity: GibbsPCDSolver vs. raking (

N{=}100{,}000

, Syn-ISTAT,

K{=}15

H_{\max}=\log|\mathcal{X}|=18.93

nats.

Figure 5: Population diversity: GibbsPCDSolver vs. generalised raking ( $N{=}100{,}000$ , Syn-ISTAT, $K{=}15$ , full training on 31 constraints). Left : raking weight distribution (log scale); the dashed line marks the uniform weight $1/N$ of GibbsPCDSolver. After convergence, raking concentrates almost all mass on $N_{\rm eff}{=}1{,}152$ effective individuals ( $1.2\%$ of $N$ ). Centre : Lorenz curve of raking weights (Gini ${=}0.951$ ); the diagonal is perfect equality. The shaded area between the two curves quantifies the weight concentration. Right : normalised diversity metrics for both solvers; ratios above each pair indicate the GibbsPCD advantage. GibbsPCDSolver maintains uniform pool weights by construction ( $N_{\rm eff}{=}N$ ), achieving an $86.8{\times}$ diversity advantage in effective sample size and $3.15$ nats higher Shannon entropy than raking—essential for realistic agent-based dynamics.

Binary constraints
	$s=1$		$s=5$		$s=10$
$N$	MRE	$t$	MRE	$t$	MRE	$t$
$1{,}000$	0.097	3	0.076	11	0.062	22
$5{,}000$	0.040	8	0.052	38	0.039	75
$20{,}000$	0.030	30	0.016	141	0.023	279
$50{,}000$	0.012	78	0.015	364	0.014	728

Table 7: A1c: GibbsPCDSolver sensitivity to pool size

N

and sweep count

s

on the Wu benchmark (

K{=}8

|\mathcal{X}|{=}6{,}144

). Top : binary constraints (

m{=}61

); exact reference

\mathrm{MRE}{=}4{\times}10^{-5}

. Bottom : ternary constraints (

m{=}73

); exact reference

\mathrm{MRE}{=}10^{-4}

. Time in seconds. Bold : Pareto-optimal entry (lowest MRE per unit time) in each block.

Method	$N$	MRE	Iters	Time (s)
Gibbs	$10{,}000$	0.092	323	119
Gibbs	$25{,}000$	0.050	435	418
Gibbs	$50{,}000$	0.032	566	1,057
Gibbs	$100{,}000$	0.032	541	2,074
Gibbs	$200{,}000$	0.021	590	4,818
Raking	$100{,}000$	$0.000^{\dagger}$	60	29

Table 8: A-ISTAT-3: GibbsPCDSolver pool-size sensitivity on Syn-ISTAT (

K{=}15

, 31 constraints,

s{=}5

\eta{=}0.01

). Raking MRE

{}_{\text{train}}{=}0

is exact by construction and not directly comparable; see table 5 .

Figure 6: A-ISTAT-3 — pool-size sensitivity on Syn-ISTAT ( $K{=}15$ , $|\mathcal{X}|{\approx}1.7{\times}10^{8}$ , 31 constraints, $s{=}5$ , $\eta{=}0.01$ ). Left : training MRE vs. pool size $N$ (log–log); the dotted line is the theoretical variance floor $1/(2\sqrt{N})$ ; observed MRE lies above the floor at all $N$ , confirming that mixing bias dominates over gradient variance in this regime. Right : Pareto frontier of MRE vs. fit time; colour encodes $\log_{10}(N)$ . The dashed vertical line marks the fit time of generalised raking ( $29$ s); raking achieves $\mathrm{MRE}_{\text{train}}{=}0$ by algebraic construction ( $\dagger$ ) and is not directly comparable to the stochastic Gibbs estimate — see table 5 for a meaningful accuracy comparison on held-out constraints.

Appendix figures & tables6 assets

Supplementary material from the paper’s appendix.

Appendix

Index	Attribute	$d_{k}$	Domain
1	sex †	2	F, M
2	age †	4	0–24, 25–49, 50–66, 67+
3	marital status	5	NeverMarried, Married, Separated, Divorced, Widowed
4	education †	4	LessThanHS, HighSchool, SomeCollege, Bachelor+
5	employment	3	Employed, Unemployed, NotInLF
6	income	5	None, Low, Medium, UpperMedium, High

Appendix

Table 9: The 15 categorical attributes of Syn-ISTAT. Anchor variables (fixed prior marginals) are marked

\dagger

|\mathcal{X}|=\prod_{k}d_{k}\approx 1.7\times 10^{8}

Attribute	Category	$P$
sex	F	0.510
sex	M	0.490
age	0–24	0.240
	25–49	0.340
	50–66	0.240
	67+	0.180

Appendix

Table 10: Unary marginal distributions for all 15 attributes. Anchor distributions are shown in italics .

#	Table	Conditioning	Conditioned
B1	age_marital	age	marital status
B2	age_employment	age	employment
B3	education_employment	education	employment
B4	employment_income	employment	income
B5	household_children	household size	has children
B6	area_transport	residence area	main transport

Appendix

Table 11: Index of the 13 binary contingency tables.

Education	Employment	None	Low	Medium	UpperMedium	High
LessThanHS	Employed	0.030	0.320	0.400	0.180	0.070
	Unemployed	0.300	0.480	0.180	0.030	0.010
	NotInLF	0.300	0.380	0.220	0.080	0.020
HighSchool	Employed	0.020	0.210	0.430	0.240	0.100
	Unemployed	0.240	0.470	0.220	0.060	0.010
	NotInLF	0.240	0.330	0.260	0.120	0.050

Appendix

Table 18

Area	Car access	Car	PublicTr.	Walking	Bike	Mixed
Urban	NoCar	0.010	0.480	0.250	0.100	0.160
	SharedCar	0.330	0.240	0.100	0.070	0.260
	OwnCar	0.560	0.150	0.080	0.050	0.160
Suburban	NoCar	0.030	0.430	0.220	0.100	0.220
	SharedCar	0.460	0.180	0.070	0.040	0.250
	OwnCar	0.720	0.080	0.040	0.030	0.130

Appendix

Table 19

Sex	Age	Employed	Unemployed	NotInLF
F	0–24	0.330	0.140	0.530
	25–49	0.610	0.100	0.290
	50–66	0.520	0.070	0.410
	67+	0.050	0.010	0.940
M	0–24	0.430	0.100	0.470
	25–49	0.870	0.060	0.070

Appendix

Table 20

Synthetic Population Generation

See all 18 papers

Mar 23, 2026

Maximum Entropy Relaxation of Multi-Way Cardinality Constraints for Synthetic Population Generation

Generating synthetic populations from aggregate statistics is a core component of microsimulation, agent-based modeling, policy analysis, and privacy-preserving data release. Beyond classical census marginals, many applications require matching heterogeneous unary, binary, and ternary constraints derived from surveys, expert knowledge, or automatically extracted descriptions. Constructing populations that satisfy such multi-way constraints simultaneously poses a significant computational challenge. We consider populations where each individual is described by categorical attributes and the target is a collection of global frequency constraints over attribute combinations. Exact formulations scale poorly as the number and arity of constraints increase, especially when the constraints are numerous and overlapping. Grounded in methods from statistical physics, we propose a maximum-entropy relaxation of this problem. Multi-way cardinality constraints are matched in expectation rather than exactly, yielding an exponential-family distribution over complete population assignments and a convex optimization problem over Lagrange multipliers. We evaluate the approach on NPORS-derived scaling benchmarks with 4 to 40 attributes and compare it primarily against generalized raking. The results show that MaxEnt becomes increasingly advantageous as the number of attributes and ternary interactions grows, while raking remains competitive on smaller, lower-arity instances.

François Pachet, Jean-Daniel Zucker

ImagineAllThePeople, Paris, France · IRD, UMMISCO, Sorbonne Universit´e, France

Synthetic Population Generation Synthetic Data Generation Constrained Generative Modeling Computational Social Science

Ask PDF

Feb 19, 2026

MGD: Moment Guided Diffusion for Maximum Entropy Generation

Generating samples from limited information is a fundamental problem across scientific domains. Classical maximum entropy methods provide principled uncertainty quantification from moment constraints but require sampling via MCMC or Langevin dynamics, which typically exhibit exponential slowdown in high dimensions. In contrast, generative models based on diffusion and flow matching efficiently transport noise to data but offer limited theoretical guarantees and can overfit when data is scarce. We introduce Moment Guided Diffusion (MGD), which combines elements of both approaches. Building on the stochastic interpolant framework, MGD samples maximum entropy distributions by solving a stochastic differential equation that guides moments toward prescribed values in finite time, thereby avoiding slow mixing in equilibrium-based methods. We formally obtain, in the large-volatility limit, convergence of MGD to the maximum entropy distribution and derive a tractable estimator of the resulting entropy computed directly from the dynamics. Applications to financial time series, turbulent flows, and cosmological fields using wavelet scattering moments yield estimates of negentropy for high-dimensional multiscale processes.

Etienne Lempereur, Nathanaël Cuvelle--Magar, Florentin Coeurdoux +2

D´epartement d’informatique, ENS, Universit´e PSL, Paris, France · Capital Fund Management, Paris, France · Coll`ege de France, Paris, France +3

Generative Modeling Constrained Generative Modeling Diffusion Sampling

Ask PDF

Oct 2, 2025

Uniform-in-time convergence bounds for Persistent Contrastive Divergence algorithms

We propose a continuous-time formulation of a noisy persistent contrastive divergence (PCD)-like method for maximum likelihood estimation (MLE) of unnormalised densities. Our approach couples parameter updates and sampling of the parametrised density in a multiscale system of stochastic differential equations (SDEs). From this formulation, we derive non-asymptotic bounds for weak test-function errors between the resulting numerical schemes and the MLE point target. The error is decomposed into numerical discretisation, slow-fast averaging, and finite-temperature concentration terms. We also introduce an efficient implementation based on explicit stabilized integrators and establish corresponding long-time error estimates. This leads to a novel method for training energy-based models (EBMs) with quantitative error guarantees.

Paul Felix Valsecchi Oliva, O. Deniz Akyildiz, Andrew Duncan

Department of Mathematics, Imperial College London

Parameter Estimation Stochastic Approximation Energy-Based Models Convergence Analysis

Ask PDF

Abstract

Synthetic Population Generation

Explore similar work

Maximum Entropy Relaxation of Multi-Way Cardinality Constraints for Synthetic Population Generation

MGD: Moment Guided Diffusion for Maximum Entropy Generation

Uniform-in-time convergence bounds for Persistent Contrastive Divergence algorithms