Equilibrium Training of Energy-Based Models with Parallel Trajectory Tempering Nicolas Béreux,1 Aurélien Decelle,2 Cyril Furtlehner,1 and Beatriz Seoane3 1
Université Paris-Saclay, CNRS, INRIA, LISN, 91190 Gif-sur-Yvette, France 2 Escuela Técnica Superior de Ingenieros Industriales, Universidad Politécnica de Madrid, 28006 Madrid, Spain 3 Departamento de Fı́sica Teórica & IPARCOS, Universidad Complutense de Madrid, 28040 Madrid, Spain.
arXiv:2607.27077v1 [cs.LG] 29 Jul 2026
Energy-Based Models (EBMs) provide an interpretable framework for generative modeling of scientific data, but poor Markov Chain Monte Carlo mixing often limits their reliability. We introduce a training algorithm based on Parallel Trajectory Tempering (PTT), which exploits the continuity of the optimization path to maintain equilibrium sampling throughout learning. This enables stable and fast training on highly multimodal and data-scarce scientific datasets. Combined with reservoir sampling and adaptive optimization, PTT has a computational cost comparable to Persistent Contrastive Divergence, making it a practical replacement for standard training methods. It also provides direct estimates of thermalization times, equilibrium samples from trained models, and accurate log-likelihoods at essentially no additional cost. Experiments on Restricted Boltzmann Machines show that PTT consistently outperforms existing EBM training approaches. On discrete tabular data, it also surpasses state-of-the-art deep generative models, yielding higherquality samples and greater robustness to overfitting and limited data. Our results make equilibrium maximum-likelihood training of EBMs practical and computationally efficient.
Advances in imaging, genome sequencing, and highthroughput experiments are generating vast highdimensional datasets across biological scales, from molecules and cells to neural circuits and populations. These data create new opportunities to study living and cognitive systems as many-body systems shaped by collective behavior, disorder, fluctuations, and nonequilibrium dynamics. Realizing this potential requires generative models that combine sufficient expressivity to capture complex dependencies with enough interpretability to reveal the organizing principles of biological systems. Modern generative models have achieved remarkable success in modeling high-dimensional data, but many do not provide explicit access to the probability distributions they learn, limiting their value for scientific discovery, while remaining largely hard to interpret. In contrast, Energy-Based Models (EBMs) represent the data distribution as a Boltzmann distribution, pθ (v) =
exp [−Eθ (v)] , Zθ
Zθ =
X
exp [−Eθ (v)] ,
{v}
(1) where Eθ (v) is a neural-network energy and (Zθ ) the partition function. EBMs connect directly to Boltzmann statistical mechanics [1]: learning infers an effective many-body Hamiltonian whose low-free-energy states encode dominant data structures. Since any energy function can be exactly decomposed into irreducible interactions of increasing order [2], EBMs provide a principled route to higher-order dependencies [3, 4]. For shallow EBMs, the energy landscape can also be tracked during training [5–7], enabling principled pattern extraction and interpretable analyses of experimental data [8], with applications in computational biology [5, 9, 10], neuroscience [11–13], statistical physics [14–16], and quantum physics [17, 18]. Despite these advantages, EBMs have
been largely overshadowed in modern machine learning, not for lack of expressive power, but because they remain difficult to train [19, 20]. Maximum-likelihood learning requires model samples, typically generated by MCMC. As training progresses, EBMs develop increasingly structured energy landscapes and undergo cascades of second-order phase transitions [21–23], whose critical slowing down makes equilibrium sampling and reliable model evaluation increasingly difficult [24]. Training therefore often relies on out-ofequilibrium chains, yielding biased gradients, unstable optimization, and poor generalization, particularly for high-dimensional and multimodal data [19, 25–30]. Contrastive Divergence (CD) [31] made EBM training practical by replacing equilibrium sampling with short chains, while Persistent Contrastive Divergence (PCD) [32] improves on this by maintaining chains across updates. Yet, on structured datasets, persistent chains may eventually fail to track the evolving model distribution and fall out of equilibrium, thereby biasing the likelihood gradient [28, 30]. Similar failures arise in datascarce regimes, where persistent chains become trapped near individual training examples, causing mixing times to diverge and training to break down abruptly, as illustrated for MNIST in the Materials & Methods. Such failure modes are often overlooked in likelihood-based comparisons, despite their strong impact on both model behavior and generative quality: poorly trained models typically exhibit anomalously slow relaxation, increased overfitting, memorization and mode collapse effects, and unreliable evaluation. Controlling equilibration throughout training is therefore particularly important in scientific applications, where interpretability is a central advantage of EBMs. Low-rank pretraining can delay the breakdown of persistent chains [30, 33], while constrained MCMC [28], population annealing [34], reweighting [35] improve ex-
2 where ⟨·⟩M denotes the average with respect to the model distribution [1]. The parameters are then updated according to
Tr ai
ni
ng
tim
e
θ ← θ + η∇θ L,
FIG. 1. Schematic illustration of the Parallel Trajectory Tempering (PTT) algorithm. As training progresses (arrow), successive model checkpoints are added to the ladder. Replicas perform a random walk between neighboring checkpoints via swap moves, ensuring efficient equilibration despite the increasingly structured energy landscape.
ploration at the cost of substantial computation or reduced effective sample size. Parallel Tempering (PT) [36–40] where the system is replicated at a given lader of temperatures {βt }, exchanges replicas across neighboring temperatures and is effective on many datasets [24, 41, 42], but requires many intermediate temperatures and fails in highly structured landscapes, where first-order transitions suppress exchanges and trap replicas [30, 33]. Motivated by a theoretical analysis of RBM learning dynamics, Parallel Trajectory Tempering (PTT) was recently introduced [30]. Rather than tempering in temperature, PTT exchanges replicas between neighboring models along the learning trajectory (Fig. 1), following Hamiltonian-exchange MCMC [43]. In this work, we make PTT practical for EBM training by combining reservoir sampling with adaptive optimization, achieving a computational cost comparable to PCD while maintaining equilibrium sampling even for highly structured or very small datasets. This largely removes the main bottleneck in EBM training—the negative phase—and enables equilibrium training, sampling, and evaluation at near-standard computational cost. I. A.
BACKGROUND
Maximum-likelihood training
EBMs are typically trained by maximizing the loglikelihood (LL) of a training dataset D = {v (µ) }M µ=1 , where M is the number of training examples. The objective function is L = ⟨log pθ (v)⟩D ,
(2)
where ⟨·⟩D denotes the empirical average over the empir PM 1 (µ) ical training distribution pD (v) = M . µ=1 δ v − v Starting from an initial set of parameters, typically chosen at random, the model parameters are updated by gradient ascent on the log-likelihood: ∇θ L = ⟨−∇θ E⟩D − ⟨−∇θ E⟩M ,
(3)
(4)
where η is the learning rate. The data term, ⟨−∇θ E⟩D , is directly estimated from the training set, whereas the model term requires averaging over an exponentially large state space and is therefore approximated by MCMC. This estimate is unbiased only after equilibration, but the thermalization time of Alternating Gibbs Sampling (AGS) typically grows rapidly during training [19, 28], making equilibrium maximum-likelihood learning impractical and motivating heuristic approaches such as CD and PCD. B.
Restricted Boltzmann Machines (RBMs)
In this work, we use RBMs [1, 44] as a paradigmatic EBM. Despite their simple bipartite architecture, RBMs are universal approximators [45] and capture high-order dependencies between visible variables through binary hidden units [46]. The absence of intra-layer couplings makes AGS particularly efficient. The RBM joint energy reads X X X Hθ (v, h) = − Wiµ vi hµ − ai vi − bµ hµ . (5) i,µ
µ
i
With θ = {W, a, b} After marginalizing the binary hidden variables, P one obtains the RBM energy function pθ (v) = Z1 h e−Hθ (v,h) = Z1 e−Eθ (v) , with " !# Nh X X X Eθ (v) = − ai vi − log 1+exp bµ + Wiµ vi . (6) i
µ=1
i
Despite their simplicity, RBMs retain the main sampling challenges of EBMs, including multimodality, slow mixing, and training-induced phase transitions. Yet, as shown below, they achieve strong generative performance on discrete and mixed-type tabular data, outperforming state-of-the-art methods scientific datasets. II.
NEW TRAINING METHOD
Recent work [28] introduced Parallel Trajectory Tempering (PTT), exploiting the smooth evolution of the model distribution during training. Rather than tempering in temperature, PTT uses a ladder of model checkpoints {Et } (also called replica) along the training trajectory, interpolating between initialization and the trained model (see sketch in Fig. 1. Configurations are exchanged between neighboring replicas according to the Metropolis rule pacc (xt ↔ xt−1 ) = min 1, e∆Et (xt )−∆Et (xt−1 ) (7)
3 where ∆Et (x) ≡ Et (x) − Et−1 (x). In practice, each PTT sweep combines one swap proposal between nearby replicas followed by k AGS steps in each model. Because PTT exploits the smooth evolution of the model distribution during training, the replica ladder can be built adaptively. Starting from the random initialization, a checkpoint is frozen whenever the swap acceptance with the previous replica falls below a threshold α, ensuring sufficient overlap between neighbors. If this occurs after a single gradient update, the learning rate is temporarily reduced; a cosine-similarity criterion is used to prevent it from becoming excessively small (see Materials & Methods). A naive implementation would simulate all replicas, with a cost growing linearly with ladder size. We avoid this by introducing a reservoir of Nres equilibrium samples. When a new checkpoint t is added, replica t − 1 is thermalized for more than 20τexp , after which samples are stored every 2τint until the reservoir is filled. Thereafter, only replicas from t onward are evolved, while exchanges with t − 1 use configurations drawn from the reservoir. Since new checkpoints are created only O(10) times during training, the full ladder is rarely simulated; most updates require evolving only the last two replicas. Beyond substantially accelerating sampling, PTT provides, essentially for free, an accurate estimate of the log-likelihood throughout and after training [30], once the replica ensemble has equilibrated. This is achieved by recursively evaluating the partition function, D E Zt = eEt−1 −Et Zt−1 , (8) t−1
where ⟨·⟩t denotes an expectation over the equilibrium distribution of replica t. Since these expectations are computed from the same Markov chains used for training, or to sample the different models after training using PTT, no additional computation cost is required. PTT also provides a direct thermalization diagnostic from replica diffusion across the ladder [47]. At equilibrium, each replica visits the Nm models uniformly, while the checkpoint-index of each replica autocorrelation function, C(t) ∼ e−t/τexp , yields the relaxation time τexp and, through its integral, τint , setting the 2τint spacing between statistically independent samples (see Materials & Methods).
III.
RESULTS
We compare PTT with standard PCD for RBM training and with state-of-the-art generative models, including Bayesian Flow Networks (BFNs) [48], across diverse scientific datasets: (i) the 2D Ising model, (ii) mouse neural spike recordings, (iii) homologous protein families, (iv) medical patient data, and (v) low-data image datasets. These problems combine strong multimodality, limited sample sizes, and strong temporal correlations. In such regimes, PCD chains can remain trapped
(a)
(b)
(c)
(d)
FIG. 2. Learning the 2D Ising model at different temperatures. Comparison of training-data statistics with equilibrium samples from PCD- and PTT-trained RBMs. (a) Magnetization distributions for L = 20: PCD exhibits mode collapse at low temperatures, whereas PTT reproduces both modes. (b) Disconnected susceptibility versus β for PCD(down) and PTT-trained (up) models across lattice sizes, compared with the data (black). (c) Number of intermediate models (temperatures) required by PTT (PT) versus β. (d) Finite-size scaling of the disconnected susceptibility for PTTtrained models using the exact 2D Ising critical exponents. The inset shows the Binder cumulant versus β, showing the expected crossing at the critical point.
in individual modes, preventing accurate estimation of their relative weights. We show that PTT-trained RBMs achieve higher test log-likelihoods, generate more diverse and higher-quality samples, are easier to evaluate and sample from than PCD-trained models, and display markedly greater robustness to overfitting. They also outperform competing generative models across these scientific datasets. For all results, 40% of the data was held out as a test set to assess generalization and detect overfitting, and models were selected at the checkpoint with the highest test log-likelihood. When generating samples from the trained RBMs, thermalization was systematically verified before saving configurations, following the protocol described in Materials & Methods.
A.
Ising 2D
Several studies have reported that standard RBM training fails to faithfully reproduce the low-temperature distribution of the two-dimensional Ising model [49–51], despite the existence of an exact RBM representation with only 2L2 hidden units [46]. As shown below, this limitation arises from inaccurate sampling during conventional CD/PCD training rather than insufficient representational power. To compare PTT and PCD, we train RBMs on the ferromagnetic Ising model on square lattices of size L×L, with L ∈ {8, 12, 16, 20}, for temperatures spanning both sides of the critical point βc = 0.440687. The distribution evolves from unimodal at high temper-
4 atures to bimodal at low temperatures, as shown by the L = 20 magnetization histograms in Fig. 2a. For each (L, β), we generate 106 equilibrium configurations using Swendsen–Wang sampling [52] and train RBMs with Nh = 2L hidden units, sufficient for an exact sparse representation [46], using PTT and PCD. At the maximum test log-likelihood, we generate Mgen = 5000 equilibrium samples using PTT for PTT-trained models and conventional PT for PCD-trained models, for which a suitable PTT ladder could not be reconstructed a posteriori. The PT ladder is built using the same overlap criterion as PTT, and the same thermalization protocol is applied. We assess generative quality by comparing thermodynamic observables between training and generated samples acrossPL and β. From the magnetization per spin, m = L−2 i si , we compute ⟨m2 ⟩, the disconnected susceptibility χdis = L2 ⟨m2 ⟩, and the Binder ratio ⟨m4 ⟩/⟨m2 ⟩2 . Error bars are estimated by bootstrap. Figure 2a compares the magnetization distributions of the data and samples generated by PCD- and PTTtrained RBMs. At low temperatures, PCD samples collapse onto a single magnetization sector (negative), whereas PTT samples correctly reproduce both modes. Figure 2b compares the disconnected susceptibility across lattice sizes: both methods agree in the paramagnetic phase, but PCD deviates sharply near and below the critical temperature, while PTT remains in excellent agreement throughout. In Fig. 2d, the disconnected susceptibility of PTT-trained RBMs follows the expected finitesize scaling of the 2D Ising model [53], while the Binder cumulant exhibits the expected crossing at the critical temperature (inset of Fig. 2d), confirming that PTT reproduces the critical behavior. Finally, we compare the efficiency of PTT with standard PT, which is expected to perform well for the Ising model, where sampling only requires crossing a second-order transition, unlike more complex multimodal datasets involving first-order barriers [30]. Figure 2c shows the number of PTT checkpoints and PT temperature replicas required to reach the same swap acceptance, α = 0.3, for the same PTT-trained RBMs. Below the critical temperature, both ladder sizes initially increase similarly with β and L, but at lower temperatures the PTT ladder saturates, while the number of PT replicas continues to grow with β.
B.
Human Genome Dataset
The Human Genome Dataset (HGD) [54] is a binary dataset derived from the 1000 Genomes Project [55], encoding the presence or absence of mutations relative to a reference genome across 805 genes selected for their high variability between geographical regions. HGD and related variants have been widely used to assess privacy We next consider a multiple sequence alignment (MSA) of homologous sequences from the beta-lactamase
leakage in generative models trained on sensitive genomic data [54, 56–58]. This dataset exhibits strong clustering along its first principal components, Fig. 3a, making training with conventional sampling particularly difficult [28, 30, 59]; its limited sample size further makes it a challenging benchmark for generative models, as shown in previous studies [54, 56]. Early in maximum-likelihood training, the phase space fragments into well-separated clusters, trapping Gibbs sampling and causing PCD to misestimate their relative weights. PTT instead exploits the continuity of the model distribution along the training trajectory to sample intermediate models efficiently. We first compare samples generated by PTT- and PCD-trained RBMs with those from BFN, a state-ofthe-art method for discrete data generation, in Fig. 3a. In each case, we select the model with the highest test log-likelihood and project generated samples onto the first two principal components of the training data (black points), with marginal one-dimensional histograms shown alongside. Both BFN and PCD either exhibit mode collapse or fail to reproduce all modes of the data distribution, and neither accurately captures the pairwise statistics (insets of Fig. 3a). By contrast, the PTT-trained RBM reproduces both the multimodal structure and pairwise statistics with high accuracy. We next use the PRIVET method [58] to quantify underfitting, overfitting, and memorization. PRIVET compares the distribution of nearest-neighbor (NN) distances between generated samples and the training set with the distribution of NN distances within the training set itself, which serves as the reference distribution [58] (More details are given in the Materials & Methods). NN distances systematically smaller than the reference indicate overfitting and memorization, whereas larger distances indicate underfitting. The objective is therefore to reproduce the reference distribution. As shown in Fig. 3b, the BFN underfits the data, while the PCD-trained RBM, owing to its mode collapse, strongly overfits the corresponding cluster. By contrast, the PTT-trained RBM almost perfectly matches the reference distribution at all distances. To further investigate the structure learned by the RBM, we analyze the organization of the free-energy landscape using the hierarchical clustering procedure introduced in Ref. [7]. As shown in Fig. 3c, samples assigned to the same free-energy minimum cluster together, revealing a clear hierarchical organization that first separates continental populations and subsequently resolves finer subpopulation structure. This indicates that the learned energy landscape captures meaningful genetic relationships without using any label information.
C.
Protein Family Sequence Data
domain family PFAM ID:PF13354, whose strongly clustered distribution makes EBM training particularly
5 (a)
575
562
556
550
545
540
972
958
944
904
809
749
683
680
493
349
276
232
194
172
1
978
962
937
896
867
827
810
797
788
756
726
77
160
71
68
62
61
60
53
52
38
26
21
20
17
12
999
992
981
980
977
974
971
968
7
955
947
946
943
933
932
931
920
918
915
0
905
899
898
894
893
888
883
878
875
869
856
850
830
817
795
790
778
773
759
757
745
737
733
730
719
718
704
703
700
685
678
668
666
662
657
642
635
619
611
601
600
715
723
193
713
229
651
274
646
301
643
319
641
338
639
340
636
360
629
370
627
410
622
420
606
478
605
510
589
631
569
659
547
691
516
701
509
507
502
492
489
476
460
445
428
422
415
409
385
359
333
332
312
306
280
255
237
234
221
219
212
180
179
154
141
127
118
107
99
98
93
80
47
44
32
30
10
6
597
325
320
304
283
279
277
271
270
257
254
253
239
236
235
227
213
209
206
203
197
195
190
184
183
178
173
168
157
156
153
147
140
135
96
78
128
124
123
122
112
121
528
525
524
515
501
487
482
464
462
457
455
448
441
438
433
414
405
402
400
FIN
397
YRI
388
South Asian
387
CEU
374
MSL
365
KHV
European
353
JPT
LWK
344
CHS
GWD
East Asian
326
ESN
American
595
Subpopulations
African
579
Continent
108
(c)
103
(b)
137
744
136
761
57
768
46
781
989
792
935
794
934
816
927
826
913
836
895
840
890
843
851
855
841
863
837
873
829
880
825
886
823
911
789
914
784
919
783
942
754
945
747
950
716
985
709
988
694
993
677
994
676 652
91
130
645
175
632
210
618
223
598
251
586
295
576
357
571
403
568
566
512
541
539
531
800
520
818
514
834
504
892
496
909
495
470
444 437
436 435 425 424
416 412
392
376
373 372 369
350
345
342
317 307
298 291
248 243 208
207 200 199
8
186
13
170
14
165
18
158
24
146
40
119
51
115
63
110
74
106
75
104
76
95
88
92
90
70
97
69
101
65
105
50
111
39
113
33
117
28
132
22
133
2
143
982
169
952
177
ACB
948
182
GBR
936
215
929
218
928
233
926
240
925
249
924
267
921
285
907
288
885
292
870
297
868
300
859
305
857
324
849
335
838
346
833
351
808
352
798
361
785
362
782
363
777 774
366
765
378
758
389
755
404
734
406
ASW
IBS
725
407
722
413
706
418
699
430
690
431
686
447
675
450
674
453
663
459
656
466
648
468
647
474
644
475
637
479
607
481
603
484
584
498
583
499
581
505
573
527
572
530
560
534
558
544
548
551
546
552
CLM
TSI
522
561
503
564
500
574
497
578
490
580
471
582
469
587
461
590
456
604
451
613
446
617
443
621
440
624
411
653
391
654
390
661
386
684
364
688
341
693
339
698
334
708
323
710
322
720
313
721
311
739
309
741
MXL
BEB
286
746
282
762
281
764
272
772
263
775
252
779
250
796
247
803
244
805
228
812
202
824
196
839
174
844
166
864
163
866
162
874
149
876
144
877
129
881
114
884
100
891
87
902
82
906
79
916
67
917
64
922
29
930
PEL
GIH
PUR
ITU
CDX
PJL
CHB
STU
25
941
11
949
4
959
997
960
979
964
973
976
887
984
865
998
821
16
814
43
740
54
738
72
724
84
696
86
667
89
557
164
513
188
266
214
155
231
145
258
134
261
957
262
951
296
939
318
903
331
835
337
687
355
673
356
671
375
634
377
626
393
616
396
609
399
602
426
592
429
485
477
463
491
454
494
452
519
401
542
394
563
382
599
379 615
290
623
222
628
217 665
45 670
41 697
31 705
9
731
991
732
961 750
872 770
853
822
742
846
536
847
532 848
442
862
299
860
871
287
882
278
889
58
996
766
36
763
66
711
81
537
94
486
125
473
131
343
439
458
465
508
518
538
543
577
588
591
614
638
712
748
793
820
897
900
908
954
966
19
969
37
995
55
56
148
181
187
191
201
204
211
264
275
303
314
329
383
395
419
220
238
293
315
316
327
367
553
565
682
767
845
854
912
983
987
205
488
567
594
608
692
707
769
771
970
171
321
511
689
717
832
852
34
224
246
268
348
139
273
347
398
434
241
330
549
612
625
226
336
35
151
523
3
559
15
83
150
161
176
198
225
256
259
358
371
727
729
736
743
751
752
753
776
791
807
815
828
842
923
967
23
27
417
42
308
73
85
265
302
408
432
555
585
649
728
760
801
804
811
901
799
813
940
953
185
986
159
49
152
230
142
289
138
294
421
427
449
483
506
526
529
535
554
570
593
596
620
630
633
640
650
658
660
669
695
702
714
735
780
786
787
802
806
819
831
858
861
879
910
938
956
963
965
975
990
5
48
59
102
109
116
120
126
167
189
192
216
242
245
260
269
284
310
328
354
368
380
381
384
423
467
472
480
517
521
533
610
655
664
672
679
681
FIG. 3. Comparison The same behavior is recovered with AIS usingof the sampling quality on HGD. (a) Scatter plots of the training dataset and the samples generated by different models projected onto the first two principal components of the dataset. The inset shows scatter plots of the empirical two-body correlations measured on the training dataset (x-axis) and on the generated samples (y-axis). The black dashed line indicates the identity. (b) Distribution of nearest-neighbor distances between generated samples and the training dataset. (c) Hierarchical clustering of test samples based on their closest free-energy minima across training of the PTT model, following Ref. [7]. Leaves denote test samples and are colored by continental origin (inner ring) and subpopulation (outer ring), revealing clear organization across continents and several subpopulations.
challenging, with potentially very long thermalization times [60]. Pairwise EBMs, often referred to as direct coupling analysis (DCA) models, have long been used to generate protein sequences and to infer residue– residue couplings that are predictive of contacts in threedimensional structures. RBMs generalize this framework by capturing effective multibody interactions [3]. This dataset therefore provides a stringent benchmark of both generative accuracy and interaction inference. We show that PTT-trained RBMs better reproduce the empirical sequence distribution than competing models while improving contact prediction over edDCA [61, 62]. We compare samples generated by PTT-trained RBMs, BFN, and edDCA after projection onto the first two principal components of the data in Fig. 4a. All three models recover the main clusters, although BFN and edDCA reproduce their relative weights less accurately. To quantify these differences, we examine the distributions of the errors in two- and three-body correlations, Fig. 4b. Despite its qualitatively accurate low-dimensional projection, edDCA shows substantially larger correlation errors than BFN, while the PTT-trained RBM outperforms both and approaches the test-set baseline. Consistently, the Earth Mover Distance [63] between generated and training samples, Fig. 4c, is lowest for PTT, although it remains above the test–train reference value. The same
trend appears in the nearest-neighbor distance distribution Fig. 4d: all models remain in the underfitting regime, but PTT lies significantly closer to the data reference. Finally, we extract the effective two-body interactions learned by the RBM from its weights [3] and compare the resulting couplings with the PF13354 contact map. We restrict the comparison to edDCA, for which pairwise couplings are explicit model parameters. Performance is quantified by the Positive Predictive Value (PPV), i.e., the fraction of true contacts among the strongest predicted couplings. Both methods recover true contacts at the top of the ranking, but as more couplings are included, the PPV of edDCA decreases more rapidly, while the PTT-trained RBM retains higher predictive accuracy.
D.
Neural recordings
We next consider Neuropixels recordings from an Allen Institute experiment [64], comprising the simultaneous activity of 2,028 neurons across cortical and subcortical regions in mice performing a visual change-detection task with familiar and novel images. Spike trains are binned into 20 ms intervals and binarized according to whether each neuron fired at least once; preprocessing details are given in Ref. [12]. A key challenge is the strong temporal structure of this dataset, which induces correlations between consecutive samples, substantially reducing the effective number of independent observations, while generating a highly multimodal distribution that makes training particularly difficult. We find that PCD performance depends strongly on training hyperparameters, in particular the number of parallel chains and the minibatch size used to estimate the likelihood gradient. To illustrate this, we train identical RBMs with PCD and PTT using 1000 parallel chains, a standard choice in practice and the value used throughout this work. As shown in Fig. 5a, PCD fails to reproduce the two- and three-point correlations of the data, whereas PTT accurately captures them. Figure 5b compares the distribution P (K) of the number K of simultaneously active neurons, a nonlinear observable sensitive to higher-order statistics. PCD exhibits a clear shift in the high-probability region relative to the train/test data, while PTT closely matches the empirical distribution, with only a small deviation in the low-probability tail. In Fig. 5, we train identical RBMs with PCD and PTT while varying only the number of parallel chains used to estimate the negative phase. Performance is quantified by the mean-squared error between the two- and three-body correlations of generated and training samples. PTT remains stable across all chain counts and consistently achieves lower errors. By contrast, PCD degrades markedly at small numbers of chains and approaches PTT performance only when using about 5000 chains.
6 (d) (a)
(b)
(c)
(e)
(f )
FIG. 4. β-Lactamase (PF13354). (a) Generated and training samples projected onto the first two principal components. (b) Distribution of errors in two- and three-body correlations relative to the training data. (c) Earth Mover Distance between generated and training samples. (d) Nearest-neighbor distance distributions to the training set. (e) Positive Predictive Value (PPV) as a function of the number of top-ranked pairwise couplings. (f ) Evolution of τexp as a function of training time for PCD- and PTT-trained models. (a) (a)
(b)
(b)
(c)
FIG. 5. Neural activity. Impact of the number of persistent chains during training. (a) MSE of two- and three-body correlations for PCD- and PTT-trained RBMs versus the number of parallel chains used to estimate the likelihood gradient during training. (b) Distribution P of the number of simultaneously active neurons, K = i ni , for models trained with 1000 chains. (c) Two- and three-body correlations of generated versus training samples for models trained with 1000 chains. The dashed line denotes identity; Pearson correlations are reported in the legend. All statistics are computed by comparing 15,000 samples from each dataset.
E.
Medical Recommendation Data
As a final benchmark, we consider the Medical Recommendation Dataset [65], a tabular dataset containing 241 patient profiles described by symptoms, causes, diagnosed illnesses, and recommended treatments.The dataset is moderately sparse and combines heterogeneous categorical variables with structured dependencies between symptoms, diagnoses, and treatments. It therefore provides a representative benchmark for tabular gen-
FIG. 6. Sampling quality on Medical data. (a) Scatter plots of generated data against training data on the first two principal components. The inset shows the two-body correlations estimated on the generated data against the train dataset. (b) Train (full line) and test (dashed line) loglikelihood of the PCD- and PTT-trained model during training.
erative modeling. Modeling tabular data remains challenging for modern generative models due to heterogeneous features, imbalanced and multimodal marginals, and complex cross-variable dependencies [66, 67]. We train RBMs with PCD and PTT, together with a BFN model, on this dataset. In order to binarize the dataset, each of the categorical variable was one hot encoded in order to binarize the dataset. As shown in Fig. 6, PTT reaches a substantially higher test loglikelihood than PCD while avoiding the undesirable dynamical behavior observed during PCD training. We also project generated samples onto the first two principal components of the data. While PTT reproduces the main structure of the empirical distribution, BFN fails to generate realistic samples and does not capture the observed modes.
IV.
CONCLUSION
In this work, we introduced a training algorithm based on an optimized sampling scheme inspired by paral-
7
Although EBM training has long been regarded as a major challenge, we showed that our approach remains reliable even for strongly multimodal and clustered dis-
tributions and in low-data regimes. Beyond stable optimization, it enables accurate log-likelihood estimation, allowing controlled studies of convergence and hyperparameter selection, including the number of hidden units and early stopping. It also provides direct diagnostics of sampling quality, since loss of equilibration can be detected from the replica-exchange dynamics. We expect these advances to facilitate the broader use of EBMs in scientific and applied settings, where their compactness, interpretability, and access to accurate likelihood estimates provide distinctive advantages over many alternative generative approaches.
[1] G. Hinton, Nobel lecture: Boltzmann machines, Reviews of Modern Physics 97, 030502 (2025). [2] N. Bulso and Y. Roudi, Restricted Boltzmann machines as models of interacting variables, Neural Computation 33, 2646 (2021). [3] A. Decelle, A. d. J. Navas Gómez, and B. Seoane, Inferring higher-order couplings with neural networks, Physical Review Letters 135, 207301 (2025). [4] A. Decelle, A. d. J. N. Gómez, and B. Seoane, Distributional simplicity bias and effective convexity in energy based models, arXiv preprint arXiv:2605.07844 (2026). [5] J. Tubiana and R. Monasson, Emergence of compositional representations in restricted Boltzmann machines, Physical review letters 118, 138301 (2017). [6] J. Tubiana, S. Cocco, and R. Monasson, Learning protein constitutive motifs from sequence data, Elife 8, e39397 (2019). [7] A. Decelle, B. Seoane, and L. Rosset, Unsupervised hierarchical clustering using the learning dynamics of restricted Boltzmann machines, Phys. Rev. E 108, 014110 (2023). [8] G. di Sarra, B. Bravi, and Y. Roudi, The unbearable lightness of restricted Boltzmann machines: Theoretical insights and biological applications, Europhysics Letters 149, 21002 (2025). [9] B. Bravi, J. Tubiana, S. Cocco, R. Monasson, T. Mora, and A. M. Walczak, Rbm-mhc: a semi-supervised machine-learning method for sample-specific prediction of antigen presentation by hla-i alleles, Cell systems 12, 195 (2021). [10] B. Bravi, Development and use of machine learning algorithms in vaccine target selection, npj Vaccines 9, 15 (2024). [11] T. L. van der Plas, J. Tubiana, G. Le Goc, G. Migault, M. Kunst, H. Baier, V. Bormuth, B. Englitz, and G. Debrégeas, Neural assemblies uncovered by generative modeling explain whole-brain activity statistics and reflect structural connectivity, Elife 12, e83139 (2023). [12] N. Béreux, G. Catania, A. Decelle, F. Mignacco, A. d. J. N. Gómez, and B. Seoane, Uncovering statistical structure in large-scale neural activity with restricted Boltzmann machines, arXiv preprint arXiv:2603.11032 (2026). [13] M. Dommanget-Kott, J. Fernandez-de Cossio-Diaz, G. Faye-Bédrin, G. Debrégeas, and V. Bormuth, Crossindividual translation of spontaneous zebrafish brain activity through a shared latent representation, Pro-
ceedings of the National Academy of Sciences 123, e2529064123 (2026). [14] R. G. Melko, G. Carleo, J. Carrasquilla, and J. I. Cirac, Restricted Boltzmann machines in quantum physics, Nature Physics 15, 887 (2019). [15] A. Barra, G. Genovese, P. Sollich, and D. Tantari, Phase diagram of restricted Boltzmann machines and generalized hopfield networks with arbitrary priors, Physical Review E 97, 022310 (2018). [16] A. Decelle and C. Furtlehner, Restricted Boltzmann machine: Recent advances and mean-field theory, Chinese Physics B 30, 040202 (2021). [17] G. Carleo and M. Troyer, Solving the quantum manybody problem with artificial neural networks, Science 355, 602 (2017). [18] Y. Nomura, A. S. Darmawan, Y. Yamaji, and M. Imada, Restricted Boltzmann machine learning for solving strongly correlated quantum systems, Physical Review B 96, 205152 (2017). [19] A. Decelle, C. Furtlehner, and B. Seoane, Equilibrium and non-equilibrium regimes in the learning of restricted Boltzmann machines, Advances in Neural Information Processing Systems 34, 5345 (2021). [20] R. Liao, S. Kornblith, M. Ren, D. J. Fleet, and G. Hinton, Gaussian-Bernoulli rbms without tears, arXiv preprint arXiv:2210.10318 (2022). [21] A. Decelle, G. Fissore, and C. Furtlehner, Spectral dynamics of learning in restricted Boltzmann machines, Europhysics Letters 119, 60001 (2017). [22] A. Decelle, G. Fissore, and C. Furtlehner, Thermodynamics of restricted Boltzmann machines and related learning dynamics, Journal of Statistical Physics 172, 1576 (2018). [23] D. Bachtis, G. Biroli, A. Decelle, and B. Seoane, Cascade of phase transitions in the training of energy-based models, NeurIPS (2024), arXiv:2405.14689 (2024). [24] O. Krause, A. Fischer, and C. Igel, Algorithms for estimating the partition function of restricted Boltzmann machines, Artificial Intelligence 278, 103195 (2020). [25] E. Nijkamp, M. Hill, S.-C. Zhu, and Y. N. Wu, Learning non-convergent non-persistent short-run mcmc toward energy-based model, Advances in Neural Information Processing Systems 32 (2019). [26] E. Nijkamp, M. Hill, T. Han, S.-C. Zhu, and Y. N. Wu, On the anatomy of MCMC-based maximum likelihood learning of energy-based models, in Proceedings of
lel tempering and designed to exploit the learning trajectories of EBMs, which have recently been shown to undergo second-order phase transitions during training. The method enables high-precision equilibrium training and allows compact RBMs to outperform more advanced generative models across a broad range of scientific datasets, despite their substantially smaller architectures and lower training cost.
8 the AAAI Conference on Artificial Intelligence, Vol. 34 (2020) pp. 5272–5280. [27] E. Agoritsas, G. Catania, A. Decelle, and B. Seoane, Explaining the effects of non-convergent sampling in the training of energy-based models, arXiv preprint arXiv:2301.09428 (2023). [28] N. Béreux, A. Decelle, C. Furtlehner, and B. Seoane, Learning a restricted Boltzmann machine using biased Monte Carlo sampling, SciPost Physics 14, 032 (2023). [29] A. Carbone, A. Decelle, L. Rosset, and B. Seoane, Fast and functional structured data generators rooted in outof-equilibrium physics, IEEE Transactions on Pattern Analysis and Machine Intelligence (2024). [30] N. Béreux, A. Decelle, C. Furtlehner, L. Rosset, and B. Seoane, Fast training and sampling of restricted Boltzmann machines, in 13th International Conference on Learning Representations-ICLR 2025 (2025). [31] G. E. Hinton, Training products of experts by minimizing contrastive divergence, Neural computation 14, 1771 (2002). [32] T. Tieleman, Training restricted Boltzmann machines using approximations to the likelihood gradient, in Proceedings of the 25th international conference on Machine learning (2008) pp. 1064–1071. [33] A. Decelle and C. Furtlehner, Exact training of restricted Boltzmann machines on intrinsically low dimensional data, Physical Review Letters 127, 158303 (2021). [34] O. Krause, A. Fischer, and C. Igel, Populationcontrastive-divergence: Does consistency help with RBM training?, Pattern Recognition Letters 102, 1 (2018). [35] D. Carbone, M. Hua, S. Coste, and E. Vanden-Eijnden, Efficient training of energy-based models using Jarzynski equality, Advances in Neural Information Processing Systems 36 (2024). [36] E. Marinari and G. Parisi, Simulated tempering: a new Monte Carlo scheme, Europhysics letters 19, 451 (1992). [37] A. Lyubartsev, A. Martsinovski, S. Shevkunov, and P. Vorontsov-Velyaminov, New approach to Monte Carlo calculation of the free energy: Method of expanded ensembles, The Journal of chemical physics 96, 1776 (1992). [38] C. J. Geyer and E. A. Thompson, Annealing Markov chain Monte Carlo with applications to ancestral inference, Journal of the American Statistical Association 90, 909 (1995). [39] M. Tesi, E. Janse van Rensburg, E. Orlandini, and S. Whittington, Monte Carlo study of the interacting self-avoiding walk model in three dimensions, Journal of statistical physics 82, 155 (1996). [40] K. Hukushima and K. Nemoto, Exchange Monte Carlo method and application to spin glass simulations, Journal of the Physical Society of Japan 65, 1604 (1996). [41] R. R. Salakhutdinov, Learning in Markov random fields using tempered transitions, Advances in neural information processing systems 22 (2009). [42] G. Desjardins, A. Courville, Y. Bengio, P. Vincent, and O. Delalleau, Tempered markov chain Monte Carlo for training of restricted Boltzmann machines, in Proceedings of the thirteenth international conference on artificial intelligence and statistics (JMLR Workshop and Conference Proceedings, 2010) pp. 145–152. [43] E. Rosta, M. Nowotny, W. Yang, and G. Hummer, Catalytic mechanism of rna backbone cleavage by ribonuclease h from quantum mechanics/molecular mechanics simulations, Journal of the American Chemical Society
133, 8934 (2011). [44] P. Smolensky, In parallel distributed processing: Volume 1 by d. rumelhart and j. mclelland (MIT Press, 1986) Chap. 6: Information Processing in Dynamical Systems: Foundations of Harmony Theory. [45] N. Le Roux and Y. Bengio, Representational power of restricted Boltzmann machines and deep belief networks, Neural computation 20, 1631 (2008). [46] A. Decelle, C. Furtlehner, A. d. J. Navas Gómez, and B. Seoane, Inferring effective couplings with restricted Boltzmann machines, SciPost Physics 16, 095 (2024). [47] R. Alvarez Baños, A. Cruz, L. Fernandez, J. Gil-Narvion, A. Gordillo-Guerrero, M. Guidetti, A. Maiorano, F. Mantovani, E. Marinari, V. Martin-Mayor, et al., Nature of the spin-glass phase at experimental length scales, Journal of Statistical Mechanics: Theory and Experiment 2010, P06026 (2010). [48] A. Graves, R. K. Srivastava, T. Atkinson, and F. Gomez, Bayesian flow networks, arXiv preprint arXiv:2308.07037 (2023). [49] D. Yevick and R. Melko, The accuracy of restricted Boltzmann machine models of ising systems, Computer Physics Communications 258, 107518 (2021). [50] J. Gu and K. Zhang, Thermodynamics of the ising model encoded in restricted Boltzmann machines, Entropy 24, 1701 (2022). [51] M. A. Valle, The capabilities of Boltzmann machines to detect and reconstruct ising system’s configurations from a given temperature, Entropy 25, 1649 (2023). [52] R. H. Swendsen and J.-S. Wang, Nonuniversal critical dynamics in Monte Carlo simulations, Physical review letters 58, 86 (1987). [53] D. J. Amit and V. Martin-Mayor, Field theory, the renormalization group, and critical phenomena: graphs to computers (World Scientific Publishing Company, 2005). [54] B. Yelmen, A. Decelle, L. Ongaro, D. Marnetto, C. Tallec, F. Montinaro, C. Furtlehner, L. Pagani, and F. Jay, Creating artificial human genomes using generative neural networks, PLoS genetics 17, e1009303 (2021). [55] . G. P. Consortium et al., A global reference for human genetic variation, Nature 526, 68 (2015). [56] B. Yelmen, A. Decelle, L. L. Boulos, A. Szatkownik, C. Furtlehner, G. Charpiat, and F. Jay, Deep convolutional and conditional neural networks for large-scale genomic data generation, PLOS Computational Biology 19, e1011584 (2023). [57] B. Yelmen and F. Jay, An overview of deep generative models in functional and evolutionary genomics, Annual Review of Biomedical Data Science 6, 173 (2023). [58] A. Szatkownik, A. Decelle, B. Seoane, N. Béreux, L. Planche, G. Charpiat, B. Yelmen, F. Jay, and C. Furtlehner, PRIVET: Privacy metric based on extreme value theory, arXiv preprint arXiv:2510.24233 (2025). [59] A. Szatkownik, L. Planche, M. Demeulle, T. Chambe, M. C. Ávila-Arcos, E. Huerta-Sanchez, C. Furtlehner, G. Charpiat, F. Jay, and B. Yelmen, Diffusion-based artificial genomes and their usefulness for local ancestry inference, bioRxiv , 2024 (2024). [60] A. P. Muntoni, A. Pagnani, M. Weigt, and F. Zamponi, adabmDCA: adaptive Boltzmann machine learning for biological sequences, BMC bioinformatics 22, 528 (2021). [61] P. Barrat-Charlaix, A. P. Muntoni, K. Shimagaki,
9 M. Weigt, and F. Zamponi, Sparse generative modeling via parameter reduction of Boltzmann machines: application to protein-sequence families, Physical Review E 104, 024407 (2021). [62] L. Rosset, R. Netti, A. P. Muntoni, M. Weigt, and F. Zamponi, adabmdca 2.0—a flexible but easy-to-use package for direct coupling analysis, in Protein Evolution: Methods and Protocols (Springer, 2026) pp. 83–104. [63] Y. Rubner, C. Tomasi, and L. J. Guibas, A metric for distributions with applications to image databases, in Sixth international conference on computer vision (IEEE Cat. No. 98CH36271) (IEEE, 1998) pp. 59–66. [64] N. A. Steinmetz, C. Aydin, A. Lebedeva, M. Okun, M. Pachitariu, M. Bauza, M. Beau, J. Bhagat, C. Böhm, M. Broux, S. Chen, J. Colonell, R. J. Gardner, B. Karsh, F. Kloosterman, D. Kostadinov, C. Mora-Lopez, J. O’Callaghan, J. Park, J. Putzeys, B. Sauerbrei, R. J. J. van Daal, A. Z. Vollan, S. Wang, M. Welkenhuysen, Z. Ye, J. T. Dudman, B. Dutta, A. W. Hantman, K. D. Harris, A. K. Lee, E. I. Moser, J. O’Keefe, A. Renart, K. Svoboda, M. Häusser, S. Haesler, M. Carandini, and T. D. Harris, Neuropixels 2.0: A miniaturized high-density probe for stable, longterm brain recordings, Science 372, eabf4588 (2021), https://www.science.org/doi/pdf/10.1126/science.abf4588. [65] Medical recommendation system, https: //www.kaggle.com/datasets/joymarhew/ medical-reccomadation-dataset/data (2023), accessed 2026-07-24. [66] L. Xu, M. Skoularidou, A. Cuesta-Infante, and K. Veeramachaneni, Modeling tabular data using conditional gan, Advances in neural information processing systems 32 (2019). [67] A. Kotelnikov, D. Baranchuk, I. Rubachev, and A. Babenko, Tabddpm: Modelling tabular data with diffusion models, in International conference on machine learning (PMLR, 2023) pp. 17564–17579. [68] R. A. Banos, A. Cruz, L. Fernandez, J. Gil-Narvion, A. Gordillo-Guerrero, M. Guidetti, A. Maiorano, F. Mantovani, E. Marinari, V. Martin-Mayor, et al., Nature of the spin-glass phase at experimental length scales, Journal of Statistical Mechanics: Theory and Experiment 2010, P06026 (2010). [69] Y. LeCun, L. Bottou, Y. Bengio, P. Haffner, et al., Gradient-based learning applied to document recognition, Proceedings of the IEEE 86, 2278 (1998).
MATERIALS AND METHODS Hyperparameters
All models are trained using two PTT sweep attempts between the reservoir and the last two checkpoints in the ladder. Each PTT sweep consists of attempting swaps of configurations between adjacent models, followed by 10 alternating Gibbs updates, during which the visible and hidden layers are updated successively. The learning rate is initialized at 10−3 and is automatically adjusted throughout training. New checkpoints are added to the ladder whenever the swap acceptance rate between the reservoir and the last checkpoint falls below α = 0.3.
Each time a new checkpoint is added, a PTT simulation involving the entire ladder is run until thermalization is reached (see below for the thermalization criterion), after which a new equilibrium reservoir is generated. The swap acceptance rates across the full PTT ladder are continuously monitored. If any of them falls below α−0.2, training is stopped, as this indicates that the Markov chains have likely drifted out of equilibrium. In this case, the maximum learning rate should be reduced before restarting training.
Adaptive learning rate
The parameters of the EBM are optimized by stochastic gradient ascent. Since PTT provides low-noise gradient estimates once the sampler has equilibrated, the relative orientation of successive gradients carries useful information about the optimization dynamics. We exploit this information to adapt the learning rate during training. Let gt denote the gradient of the log-likelihood at optimization step t. We quantify the alignment between two consecutive gradients through their cosine similarity,
cossim(gt , gt−1 ) =
gt⊤ gt−1 . ∥gt ∥ ∥gt−1 ∥
(9)
When consecutive gradients are well aligned (cossim ≈ 1), the optimization consistently follows the same direction, indicating that a larger learning rate can safely accelerate convergence. Conversely, negative cosine similarity signals that successive updates cancel each other out, indicating instability and the need to reduce the learning rate. A value close to zero suggests orthogonal updates, meaning there is no clear direction to follow, and the learning rate should remain stable. We therefore adjust the learning rate dynamically according to the measured cosine similarity, increasing it when gradients remain aligned, decreasing it when they become anti-aligned and keeping it stable when the updates are decorrelated. Finally, in the case of RBMs, the cosine similarity is computed separately for the visible bias, hidden bias and weight matrix to avoid a noisy gradient in one component disrupting the others. An additional safeguard is provided by the PTT dynamics itself. The swap acceptance between neighboring checkpoints measures the overlap of their equilibrium distributions. If it falls below a prescribed threshold, that we fix always to α = 0.3, the optimization step is rejected, the parameters are restored to the previous checkpoint, and the learning rate is halved before training resumes. This prevents overly large parameter updates from disrupting the PTT ladder while automatically adapting the learning rate to the local optimization landscape.
10 Log-likelihood estimation
Beyond substantially accelerating sampling, PTT provides, at essentially no additional computational cost, an accurate estimate of the log-likelihood throughout and after training [30], provided that the replica ensemble has equilibrated. The key observation is that the partition functions of two consecutive models in the ladder satisfy X X Zt = e−Et (v) = e−Et−1 (v) eEt−1 (v)−Et (v) v
v
= Zt−1 eEt−1 −Et t−1 ,
(10)
where ⟨·⟩t−1 denotes an expectation over the equilibrium distribution of model t − 1. Starting from the first checkpoint, whose partition function is known analytically, the partition function of every subsequent model is obtained recursively as log Zt = log Zt−1 + log eEt−1 −Et t−1 .
(11)
The expectation is estimated using the equilibrium configurations already generated by PTT. Since neighboring checkpoints are introduced only when their overlap is sufficiently large, the energy difference Et − Et−1 remains small and the exponential reweighting factor has low variance. Once log Zt is known, the average log-likelihood of a dataset D = {v µ }M µ=1 follows directly from M M 1 X 1 X Lt = log pt (v µ ) = − Et (v µ ) − log Zt . (12) M µ=1 M µ=1
Furthermore, PTT also provides a simple diagnostic of thermalization by monitoring the diffusion of replicas across the ladder [47] (see Materials & Methods). Thermalization criterion and reservoir creation
To ensure equilibrium sampling, we follow the approach of Ref. [68]. We run nchains = 1000 independent PTT simulations (in this work referred to as parallel chains), each containing Nm replicas evolving simultaneously at different model checkpoints. During the simulation, replicas diffuse along the ladder through swap moves. In equilibrium, every checkpoint is visited with equal probability. Denoting by ni (t) the checkpoint occupied by replica i at Monte Carlo time t, we compute the normalized autocorrelation function C(t) =
which defines the exponential autocorrelation time τexp , corresponding to the slowest relaxation mode of the random walk through the ladder. Whenever a new checkpoint is added, the nchains parallel chains are first evolved for at least 20 τexp to ensure thermalization. After thermalization, statistically independent configurations are generated by recording one configuration every 2τint Monte Carlo steps, where the integrated autocorrelation time is obtained from the self-consistent relation
⟨(ni (t + t0 ) − ⟨n⟩) (ni (t0 ) − ⟨n⟩)⟩ D E , 2 (ni (t0 ) − ⟨n⟩)
(13)
where ⟨·⟩ denotes an average over all replicas, parallel chains, and time origins t0 , and ⟨n⟩ = (Nm − 1)/2. At long times, C(t) ∼ A e−t/τexp ,
(14)
τint =
6τint 1 X + C(t). 2 t=0
(15)
Whenever a new checkpoint is added to the PTT ladder during training, this protocol is used to construct a reservoir of Nres = 10Nchains thermalized and statistically independent configurations. The same thermalization and sampling protocol is used throughout this work to generate all reported samples.
PRIVET assessment of under- and overfitting
To assess under- and overfitting, we use PRIVET [58], a nearest-neighbor method based on extreme-value statistics. For a sample x and a reference set D, we define the nearest-neighbor distance δ(x, D) = min d(x, y). y∈D
We compare the distribution of δ(x, Dtrain ) for generated samples with a reference distribution obtained from distances between two disjoint subsets of the real data. If generated samples are systematically farther from the training data than this reference, the model underfits; if both distributions agree, the model is consistent with good generalization; if generated samples are systematically closer, this indicates overfitting and possible memorization. PRIVET formalizes these deviations using an extreme-value fit to the nearest-neighbor distribution, allowing statistically significant excesses of unusually small distances to be identified.
Problems with PCD in data-scarce regimes
Here we illustrate the failure of PCD-trained EBMs in data-scarce regimes using binarized MNIST [69]. While PCD performs well on the full 50,000-sample training set, reducing the dataset size exposes severe sampling pathologies. We train PCD-RBMs on subsets of size M and save checkpoints throughout training. These checkpoints are assembled a posteriori into a PTT ladder like in Ref. [30] to accelerate sampling, assess thermalization, and estimate the log-likelihood. Figure 7a shows the training and test log-likelihoods estimated from short 1000-sweep
11 PTT runs. Both eventually decrease sharply, including the training likelihood despite it being the optimized objective. The same behavior is recovered with AIS using sufficiently many chains and intermediate temperatures. We also verified this behavior for RBMs with sufficiently few hidden units to allow exact log-likelihood computation by enumeration, confirming that the trainLL collapse persists in low-data regimes. By contrast, PTT-trained RBMs display standard overfitting, with the training likelihood continuing to increase after the test likelihood decreases (Fig. 7b). The breakdown is accompanied by a sharp growth of the exponential autocorrelation time τexp (Fig. 7c). Beyond the red star, equilibration cannot be achieved within 5 × 104 PTT sweeps, whereas relaxation remains controlled for PTT-trained models. Figure 7d further shows that the freezing corresponds to trajectories becoming trapped near individual training examples and generating memorized copies after only a few Gibbs steps. This sudden arrest in the relaxation dynamics is commonly observed in PCD-trained models as we show for instance in Fig. 4–f for proteins. (a)
(b) LL
0
LL
−100
0 −100 −200 −300
M 125 500 2000 −300 103 104 105 Training time (gradient updates)
PTT τexp
1000 −200
100
(d) 101 AGS steps
train
test
10 1
103 AGS steps
(c) 103 Gradient updates
PTT PCD
105 AGS steps
104
Closest sample in dataset
FIG. 7. Failure modes of PCD training in datascarce regimes. (a) Training (solid) and test (dashed) loglikelihoods during PCD training for different dataset sizes M . (b) For M = 100, comparison between PCD (red, learning rate 0.01) and PTT (green, adaptive learning rate with 0.01 maximum learning rate). Red stars mark the onset of the dynamical slowdown, where thermalization requires more than 5 × 104 PTT sweeps, according to our 20τexp criterion. Dark curves use thermalized PTT estimates, whereas shaded curves use short 1000-sweep runs. (c) Corresponding exponential autocorrelation times τexp along training. (d) Five independent alternating Gibbs trajectories of 105 sweeps, initialized from the random configurations shown on the left and using the PCD model marked by the red star in (b). The trajectories rapidly become trapped near individual training examples and remain confined there.