Paper deep dive
Parallel Noising in Neural Markov Logic Networks
Peter Jung, Giuseppe Marra, Ondrej Kuzelka
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 92%
Last extracted: 7/23/2026, 3:27:46 AM
Summary
The paper introduces NMLN*, an enhanced version of Neural Markov Logic Networks (NMLNs) that addresses limitations in expressiveness and inference speed. It replaces fragment-based potentials with global Graph Neural Network (GNN) energies and introduces 'parallel noising,' a replica-exchange MCMC algorithm using a ladder of data-space corruptions (Bernoulli noise) instead of temperature ladders. This approach improves mixing and sample quality for graph generation tasks, particularly molecular structures, outperforming previous NMLN variants and matching specialized recurrent models.
Entities (7)
Relation Signals (6)
NMLN* → uses → Graph Neural Networks
confidence 95% · NMLN* combines two orthogonal changes: global graph neural network (GNN) energies
NMLN* → uses → Parallel Noising
confidence 95% · NMLN* combines... a corruption-defined replica-exchange sampler that we call parallel noising
Parallel Noising → isvariantof → Replica-Exchange MCMC
confidence 92% · Parallel noising is an exact replica-exchange algorithm for an arbitrary ladder of target distributions
Parallel Noising → uses → Bernoulli Corruption
confidence 90% · Parallel noising... replaces temperature ladders with a ladder of data-space corruptions... Bernoulli corruption levels
NMLN* → improvesperformanceon → Molecular Graph Generation
confidence 88% · NMLN* substantially improves mixing and sample quality on molecular graph generation benchmarks
NMLN* → outperforms → Neural Markov Logic Networks
confidence 85% · NMLN* substantially improves mixing and sample quality... compared to prior NMLNs
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Neural Markov Logic Networks (NMLNs) are a flexible neurosymbolic relational model. Previous work has shown that, although NMLNs achieve strong performance as generative models for small relational structures, they underperform diffusion-based generative graph models on larger structures. In this paper, we strengthen NMLNs along two main dimensions: (i) we increase the expressive capacity of their potential functions using graph neural networks, and (ii) we develop a new training and inference algorithm inspired by parallel-tempering Markov chain Monte Carlo methods, which we name parallel noising. Together, these enhancements enable NMLNs to attain strong performance in graph generation relative to general diffusion-based generative graph models. Furthermore, they allow NMLNs to match the performance of specialized text-based recurrent models when generating small molecular structures.
Tags
Links
- Source: https://arxiv.org/abs/2607.19126v1
- Canonical: https://arxiv.org/abs/2607.19126v1
Trouble viewing inline? Open PDF directly →
Full Text
84,278 characters extracted from source content.
Expand or collapse full text
Parallel Noising in Neural Markov Logic Networks Peter Jung1, Giuseppe Marra2, Ondřej Kuželka1 1Czech Technical University, Prague, Czech Republic 2Department of Computer Science, KU Leuven, Leuven, Belgium Abstract Neural Markov Logic Networks (NMLNs) are a flexible neurosymbolic relational model. Previous work has shown that, although NMLNs achieve strong performance as generative models for small relational structures, they underperform diffusion-based generative graph models on larger structures. In this paper, we strengthen NMLNs along two main dimensions: (i) we increase the expressive capacity of their potential functions using graph neural networks, and (i) we develop a new training and inference algorithm inspired by parallel-tempering Markov chain Monte Carlo methods, which we name parallel noising. Together, these enhancements enable NMLNs to attain strong performance in graph generation relative to general diffusion-based generative graph models. Furthermore, they allow NMLNs to match the performance of specialized text-based recurrent models when generating small molecular structures. 1 Introduction Neural Markov Logic Networks (NMLNs) are a flexible neurosymbolic relational model that combines the maximum-entropy semantics of Markov logic (Richardson and Domingos, 2006) with learnable (neural) relational potentials (Marra and Kuželka, 2021; Jung et al., 2024). A persistent obstacle for using NMLNs in practice is inference: likelihood training and downstream tasks require sampling from high-dimensional, multi-modal distributions over possible worlds, yet classical single-chain Gibbs sampling mixes poorly and is brittle in the presence of near-deterministic constraints (Marra and Kuželka, 2021). A natural first attempt to improve mixing is parallel tempering (also known as replica exchange (Swendsen and Wang, 1986; Geyer, 1991; Hukushima and Nemoto, 1996)), which runs multiple chains targeting related “flattened” distributions and swaps their states so that “hot” chains can help the “cold” chain cross energy barriers. In practice, however, as also confirmed by our initial experiments, temperature ladders can be brittle for NMLNs: the effective energy scale depends on model parameters and can drift substantially during training, so a ladder that yields good swap acceptance early on may quickly become ineffective. We therefore propose parallel noising: a replica-exchange MCMC scheme for NMLNs that replaces temperature ladders with a ladder of data-space corruptions. Concretely, we fix Bernoulli corruption levels 0≤νN<⋯<ν1≤10≤ _N<·s< _1≤ 1 and define a ladder of intermediate targets indexed by these noise rates. During inference, each level performs within-level MCMC updates and we regularly attempt swaps of adjacent levels using the standard replica-exchange Metropolis ratio, evaluated under the (unnormalized) NMLN log-potentials at the two noise levels rather than under temperature-rescaled energies. The lowest-noise chain targets the desired NMLN, while higher-noise chains are easier to explore and help the sampler cross energy barriers. Parallel noising offers two advantages that we make precise: (i) it is an exact replica-exchange algorithm for an arbitrary ladder of target distributions, so the lowest-noise marginal is the desired target at stationarity (Proposition 1); and (i) because the ladder is defined by a fixed corruption operator rather than an energy rescaling, it remains well-behaved under time-varying targets (e.g., during training), where temperature schedules are difficult to tune (Appendix A.3 (toy example)). Contributions. • Parallel noising for NMLNs: we introduce a corruption-defined replica-exchange sampler that stays well-behaved under time-varying targets during learning, and combine it with global GNN energies to form NMLN*. • Why it works: we characterize swap acceptance via distributional overlap and derive simple, distribution-free bounds for Bernoulli noising schedules. • Empirical impact: NMLN* substantially improves mixing and sample quality on molecular graph generation benchmarks compared to prior NMLNs and strong generative baselines (Sec. 8). Organization. Section 3 reviews replica-exchange MCMC and the NMLN learning objective. Section 4 defines NMLNs and the Bernoulli corruption operator used to build a noise ladder. Section 5 summarizes the main bottlenecks of existing NMLNs—limited expressiveness of fragment-based potentials and slow-mixing Gibbs inference—and motivates our two remedies. Section 6 introduces NMLN* with global GNN energies and the parallel noising sampler. Section 7 develops theoretical properties and connects parallel noising to parallel tempering. Section 8 presents experiments. 2 Related Work Replica-exchange MCMC and tempering. Replica exchange (parallel tempering) couples chains targeting a ladder of related distributions and proposes swap moves so that exploratory (“hot”) chains help the target (“cold”) chain traverse energy barriers (Swendsen and Wang, 1986; Geyer, 1991; Hukushima and Nemoto, 1996). Most practical variants define the ladder by temperature scaling, and performance depends on maintaining sufficient overlap between adjacent temperatures. For discrete energy-based models, this can be brittle because the effective energy scale can drift during learning, motivating adaptive and non-reversible tempering variants (Desjardins et al., 2010; Miasojedow et al., 2013; Syed et al., 2022). Corruption/noising for discrete structures. Corruption operators are widely used to smooth discrete distributions and stabilize learning, including in denoising/diffusion models for graphs and molecules (Vignac et al., 2023). In NMLNs, injecting noise during training was introduced as a pragmatic way to avoid near-deterministic constraints that can break Gibbs sampling (Marra and Kuželka, 2021). Our contribution is to use corruption to define a replica-exchange ladder: PN remains an exact sampler for the chosen intermediate targets, and these targets remain meaningful even when the underlying model evolves during training. Neurosymbolic probabilistic logic and inference. NeSy systems include constraint- and fuzzy-logic approaches (Diligenti et al., 2017; Van Krieken et al., 2022; Badreddine et al., 2022; Xu et al., 2018) and probabilistic-logic approaches (Manhaeve et al., 2018; Marra et al., 2020; Li et al., 2023; Ahmed et al., 2022; Winters et al., 2022; van Krieken et al., 2025). Within probabilistic NeSy, undirected energy-based models—including Markov logic networks (Richardson and Domingos, 2006) and probabilistic soft logic models (Pryor et al., 2023)—are particularly relevant here; relational neural machines provide another closely related energy-based formulation (Marra et al., 2020). NMLNs (Marra and Kuželka, 2021; Jung et al., 2024) sit in this family but replace hand-specified first-order rules with neural potentials that can learn soft constraints from data. Exact inference in probabilistic NeSy systems such as DeepProbLog uses weighted model counting (Manhaeve et al., 2018), which is #P\#P-hard in general (Abboud et al., 2020). To scale inference and learning, practical systems use several distinct computational strategies: A*-like search over high-probability proofs (Manhaeve et al., 2021); provenance-semiring-based differentiable reasoning (Li et al., 2023); learned neural approximations to symbolic inference (van Krieken et al., 2023); neural approximate model counting (Abboud et al., 2020); sampling-based learning objectives (Verreet et al., 2024); specialized gradient estimators for categorical random variables (De Smet et al., 2023); hybrid approximate Bayesian inference, automated reasoning, and gradient estimation for sequential NeSy models (De Smet et al., 2025); and general MCMC methods (Robert and Casella, 1999), including Gibbs sampling in NMLNs (Marra and Kuželka, 2021). Our focus is on making MCMC-style inference effective even when the learned energy encodes highly deterministic structure. Generative NeSy models. A subset of NeSy models extends these frameworks to generative tasks, e.g., via constrained deep generative models, including GANs (Di Liello et al., 2020; Stoian et al., 2024), VAEs (Misino et al., 2022), sequential neurosymbolic Markov models (De Smet et al., 2025), or diffusion models (Huang et al., 2024; Scassola et al., 2023). These methods span heterogeneous output domains, including tabular records, symbolic music, sequential latent states, and perceptual data. By contrast, NMLNs directly define joint distributions over complete symbolic worlds and can be used as generative models for graphs and relational databases. van Krieken et al. (2025) use diffusion to model dependencies among symbolic concepts in perceptual NeSy pipelines, whereas NMLNs directly model complete symbolic worlds. 3 Preliminaries We briefly review replica-exchange MCMC (parallel tempering) notation and highlight the acceptance mechanism that underlies PN. The key idea is to run a ladder of related MCMC chains, each targeting a related distribution, and allow state swaps so that exploratory chains help the target-level chain cross energy barriers. A more detailed review (including an illustrative example) is in Appendix A.1. 3.1 Replica Exchange and Parallel Tempering Replica exchange runs N Markov chains targeting a ladder of distributions πii=1N\ _i\_i=1^N and regularly proposes swaps between adjacent levels. Writing πi(ω)∝exp(ρi(ω)) _i(ω) ( _i(ω)), the swap of states (ωi,ωi+1)( _i, _i+1) is accepted with probability αi,i+1 _i,i+1 =min1,exp(ρi(ωi+1)−ρi(ωi) = \1, ( _i( _i+1)- _i( _i) (1) +ρi+1(ωi)−ρi+1(ωi+1)). + _i+1( _i)- _i+1( _i+1) ) \. Overlaps between adjacent levels govern swap rates: when πi _i and πi+1 _i+1 concentrate on very different regions, swaps are rarely accepted and the ladder decouples. Practical schedules therefore aim to keep adjacent targets sufficiently close (e.g., by geometric temperature spacing in PT or by gradual corruption increments in PN). 4 Neural Markov Logic Networks 4.1 Model definition and learning In Neural Markov Logic Networks, we consider a function-free first-order logic language ℒL, which is built from a set of constants ℒC_L and predicates ℛℒ=⋃iℛiR_L= _iR_i, where ℛiR_i contains the predicates of arity i. For c1,c2,…,cm∈ℒc_1,c_2,…,c_m _L and R∈ℛmR _m, we call R(c1,c2,…,cm)R(c_1,c_2,…,c_m) a ground atom. We define possible world ω to be the pair (,)(C,A), where ⊆ℒC _L, A is a subset of the set of all ground atoms that can be built from the constants in C and any relation in ℛℒR_L. We define Ωℒ _L to be the set of all possible worlds over ℒL. Intuitively, a given possible world defines a set of truetrue facts one can state using the constants (entities) and the relations of the language ℒL. Let Φ(ω;):Ωℒ→ℝ (ω;w): _L be a parametric relational potential implemented by a relational neural network that maps a possible world ω to a scalar score. Given a set of potential functions Φ1 _1, …, Φm _m, a neural Markov logic network (NMLN) (Marra and Kuželka, 2021) is the parametric exponential-family distribution over possible worlds from a given Ωℒ _L: P(ω)=1Zexp(∑iβiΦi(ω;i)),P(ω)= 1Z ( _i _i _i(ω;w_i) ), where βi _i and iw_i are parameters and Z=∑ω∈Ωℒexp(∑iβiΦi(ω;i))Z= _ω∈ _L ( _i _i _i(ω;w_i) ) is the normalization constant (partition function). Given a training set of possible worlds Y, NMLNs can be learned by maximizing the following log-likelihood by some gradient-based method (Marra and Kuželka, 2021): maxi,βi∑ω^∈∑i=1mβiΦi(ω^;i)−logZ. _w_i, _i \ _ ω _i=1^m _i _i( ω;w_i)- Z \. (2) The gradient of the log-likelihood function of NMLNs, given the data ω ω (here ω ω is an example from the training data, represented as a possible world) takes the form: ∂logP(ω^)∂wi,j ∂ P( ω)∂ w_i,j =βi(∂Φi(ω^;i)∂wi,j−ω∼P[∂Φi(ω;i)∂wi,j]) = _i ( ∂ _i( ω;w_i)∂ w_i,j-E_ω P\! [ ∂ _i(ω;w_i)∂ w_i,j ] ) (3) ∂logP(ω^)∂βi ∂ P( ω)∂ _i =Φi(ω^;i)−ω∼P[Φi(ω;i)] = _i( ω;w_i)-E_ω P\! [ _i(ω;w_i) ] (4) Hence, what we need to be able to compute in order to perform the gradient-based learning are the four types of quantities: Φi(ω^;i) _i( ω;w_i), ∂Φi(ω^;i)∂wi,j ∂ _i( ω;w_i)∂ w_i,j, ω∼P[∂Φi(ω;i)∂wi,j]E_ω P\! [ ∂ _i(ω;w_i)∂ w_i,j ], and ω∼P[Φi(ω;i)]E_ω P[ _i(ω;w_i)]. The first one can be computed by evaluating the potential function Φi _i on the training example ω ω, the second can also be obtained jointly for all j using the standard backpropagation algorithm. The remaining two are more difficult to compute because they involve expectation over samples from the distribution given by the current weights. We can use sampling to approximate these expectations. In particular, Gibbs Sampling is used in the original version of NMLNs, where it is used not only for computing the gradients but also for predictions—both for marginal and conditional inference (i.e., computing marginal and conditional probabilities) and for sampling in generative settings. A significant drawback of Gibbs sampling is that it often requires a large number of steps before converging to the target distribution (Casella and Berger, 2024). 4.2 Bernoulli corruption and corruption ladders A practical difficulty for MCMC in relational models is the prevalence of near-deterministic structure: hard or almost-hard constraints can create extremely sharp modes separated by regions of negligible probability mass. In such regimes, single-chain Gibbs sampling mixes poorly, because most local moves are rejected or trapped in a narrow basin. We therefore introduce a Bernoulli corruption operator on possible worlds. Fix a domain size and identify a possible world ω with the binary vector of its ground atoms in a fixed vocabulary; let d denote the number of such atoms. For a noise parameter ν∈[0,1]ν∈[0,1], define the bit-flip channel KνK_ν by independently flipping each atom with probability ν. That is, for worlds ω,ω′∈Ωω,ω ∈ , Kν(ω→ω′)=∏j=1d[(1−ν)ωj′=ωj+ν 1ωj′≠ωj].K_ν(ω→ω )\;=\; _j=1^d [(1-ν)1\ω _j= _j\+ν\,1\ω _j≠ _j\ ]. (5) For any target distribution π on Ω (e.g., the NMLN defined by learned potentials), the ν-noised distribution is the pushforward πν(ω′)=∑ω∈Ωπ(ω)Kν(ω→ω′). _ν(ω )\;=\; _ω∈ π(ω)\,K_ν(ω→ω ). Corruption for molecular graphs. The theoretical development treats ω as a binary vector of ground atoms. In the molecule experiments, each molecule is represented as a finite relational structure over a fixed vocabulary of predicates encoding (i) node labels (atom types) and (i) edge labels (bond types). Concretely, we use unary predicates t(v) Atom_t(v) for atom type t at node v and binary predicates b(u,v) Bond_b(u,v) for bond type b on an (undirected) edge u,v\u,v\. Categorical attributes are encoded by one-hot groups of ground atoms. Our corruption operator acts independently across these groups but preserves the one-hot constraints within each group: with probability 1−ν1-ν we keep the current category, and with probability ν we resample a category uniformly from the remaining options.111Equivalently, this is a simple symmetric categorical channel; the bit-flip channel is recovered for truly binary predicates. For undirected bonds we corrupt only pairs with u<vu<v and mirror the result to enforce symmetry. We do not hard-enforce chemical validity (valence constraints, aromaticity, etc.) during MCMC; instead, validity is an evaluation property (Section 8). In practice, the learned energy assigns low probability to invalid structures once such constraints are captured by the potentials. Intuitively, increasing ν smooths the distribution and reduces the severity of narrow deterministic basins. In Section 6, we use these noised distributions to construct a ladder of intermediate targets for replica exchange. 5 Limitations of Existing NMLNs Section 4 defined NMLNs and reviewed the Bernoulli corruption operator that has been used in prior work to smooth otherwise sharp energies. While NMLNs are conceptually flexible, existing instantiations and learning pipelines still run into two bottlenecks that have limited their performance on larger structured domains: (i) local potential families whose effective receptive field is controlled by a fragment width parameter, and (i) sampling-based inference that relies on Gibbs-style local updates and becomes brittle as the learned energy approaches deterministic constraints. These limitations motivate two changes: GNN-based relational potentials and the replica-exchange sampler which we describe in Section 6. 5.1 Limited Expressive Power In the original NMLN formulation (Marra and Kuželka, 2021) (and its recent extensions (Jung et al., 2024)), potentials are built by aggregating scores of small induced fragments. This also means that capturing long-range relational regularities (e.g., connectivity or path-based constraints in graphs) typically requires increasing the fragment width, quickly leading to O(nk)O(n^k) fragments and correspondingly expensive computation. This intuition is formalized by Jung et al. (2024): they show (Theorem 1 in (Jung et al., 2024)) that NMLNs with potential functions of width k can represent the same distributions as classical Markov Logic Networks (Richardson and Domingos, 2006) whose formulas use at most k variables, but without quantifiers or constants. The restriction to quantifier-free, bounded-variable logic is limiting in practice; for example, even simple global graph properties such as “no isolated vertices” cannot be captured with constant width (see (Jung et al., 2024) for additional examples). The work of Jung et al. (2024) also proposes more expressive variants and relates them to Quantified Markov Logic Networks (Gutiérrez-Basulto et al., 2018), which allow ∃/∀∃/∀ quantifiers in prenex normal form with at most k variables. However, as soon as the desired dependency spans longer paths, prenex restrictions again force k to grow, and the model size and computation scale in the same prohibitive way as increasing fragment width. These results motivate potential families with hierarchical aggregation and an expanding receptive field at fixed parameter size. In short, fragment-based potential families face a hard tradeoff between expressivity and scalability: capturing long-range relational structure typically forces the fragment width to grow, leading to an O(nk)O(n^k) blow-up in the number of induced fragments and in the cost of evaluating the energy. 5.2 Inference (Sampling) The second limitation, already highlighted in the introduction, is inference. Exact marginal inference and exact sampling in NMLNs are generally intractable under standard complexity assumptions.222This follows from the complexity of first-order model counting (Beame et al., 2015), since WFOMC can be solved given access to an oracle for marginal inference in NMLNs. This follows the same reasoning as similar arguments done before for Markov Logic Networks in the literature, exploiting the fact that classical MLNs can be represented as NMLNs with suitably chosen potential functions. Accordingly, existing NMLN implementations rely on (variants of) Gibbs sampling in practice.333Marra and Kuželka (2021) also introduce variants of blocked Gibbs sampling to speed up inference, but their implementation still relies on Gibbs-style local updates. Unfortunately, Gibbs mixes poorly for the kinds of rugged, multi-modal energies that arise once the model starts to encode near-deterministic constraints. In classical MLNs one can sometimes handle hard constraints using MC-SAT (Poon and Domingos, 2006), but MC-SAT requires an explicit logical encoding of the deterministic clauses. In NMLNs, constraints are encoded implicitly by neural potentials, so this route is unavailable. The original NMLN paper therefore injected noise during training: at the beginning of each epoch, each ground atom of a training world is flipped with probability πn _n (Marra and Kuželka, 2021). While this can make Gibbs sampling feasible, it couples modeling and inference: the model is optimized for a corrupted data distribution, and truly deterministic structure can only be approximated. As a result, the amount of corruption becomes a brittle hyperparameter and sample quality degrades when the task requires crisp constraints. 5.3 Transient corruption as an auxiliary ladder This motivates using the corruption operator KνK_ν (Section 4.2) only as a transient exploration device: we keep the (nearly) clean NMLN at the lowest noise level νN _N as the learning target—typically νN=0 _N=0, but it can also be a small positive noise—and use higher corruption levels to define a sequence of auxiliary intermediate distributions that are easier to mix. In parallel, we alleviate the expressivity bottleneck with global GNN energies (Section 6.1). Building on these two ingredients, parallel noising runs a replica-exchange sampler across noise levels and swaps states, allowing highly corrupted chains to traverse modes while preserving the lowest-noise chain as an exact sampler of the desired target distribution (Section 6.2). In our experiments we set the target level to a small but nonzero noise, νN=10−3 _N=10^-3, so the lowest-noise chain samples from a lightly corrupted target. Setting νN=0 _N=0 recovers the fully clean NMLN when desired. More generally, one can choose νN _N to be any sufficiently small positive noise if the application calls for a slightly smoothed target. Higher-noise levels are auxiliary and exist solely to accelerate exploration of the target chain. 6 NMLN*: Global Potentials and Parallel Noising We propose NMLN*, an NMLN variant that addresses the two bottlenecks from Section 5: (i) limited expressive power of fragment-based potentials, and (i) poor mixing of Gibbs-style inference on sharp, multi-modal energies. NMLN* combines two orthogonal changes: global graph neural network (GNN) energies (Section 6.1) and a corruption-defined replica-exchange sampler that we call parallel noising (Sections 6.2–6.5). 6.1 Global potentials via GNNs Prior NMLN variants commonly rely on local fragment-based potentials (e.g., DeepSet-style aggregation), which can miss global structural signals needed for larger relational structures (and larger molecules). We therefore instantiate Φ(ω;) (ω;w) as a global GNN energy: we convert a world ω into a labeled graph whose node/edge labels correspond to unary/binary predicates (atom/bond types in molecules), apply message passing with relation-specific parameters (e.g., an R-GCN-like layer), pool node representations, and map the pooled vector to a scalar energy via an MLP (Schlichtkrull et al., 2018). This change is independent of PN: any potential family can be used at each noise level, but sharper (more expressive) targets typically make mixing harder, increasing the benefit of replica exchange. Additional architectural details are given in Appendix B. 6.2 Parallel noising: replica exchange over corruption levels Parallel noising follows the standard replica-exchange template (Swendsen and Wang, 1986; Geyer, 1991; Hukushima and Nemoto, 1996) (Section 3.1), but replaces temperature scaling with a ladder indexed by corruption levels. Noise ladder and level-wise targets. Fix noise levels 0≤νN<νN−1<⋯<ν1≤10≤ _N< _N-1<·s< _1≤ 1, where νN _N is the smallest noise (the target level). For each level i∈1,…,Ni∈\1,…,N\ we define a target distribution πi _i with an unnormalized log-density ρi:Ω→ℝ _i: : πi(ω)=exp(ρi(ω))Zi,Zi=∑ω′∈Ωexp(ρi(ω′)). _i(ω)\;=\; ( _i(ω))Z_i, Z_i= _ω ∈ ( _i(ω )). In the learning setup of Section 6.3, each πi _i is an NMLN with its own parameters θi _i trained on data corrupted at rate νi _i; we then write ρi(ω)=ρθi(ω) _i(ω)= _ _i(ω). The sampling algorithm itself only requires that we can evaluate each ρi(ω) _i(ω) up to an additive constant. Within-level updates and swaps. Parallel noising targets the product distribution on ΩN ^N, Π(ω1,…,ωN)=∏i=1Nπi(ωi), ( _1,…, _N)\;=\; _i=1^N _i( _i), and alternates: (i) within-level updates using any Markov kernels KiK_i that leave πi _i invariant (e.g., blocked Gibbs / Metropolis-within-Gibbs on ωi _i), and (i) swap proposals between adjacent levels. A swap proposal exchanges (ωi,ωi+1)↦(ωi+1,ωi)( _i, _i+1) ( _i+1, _i) and is accepted with the Metropolis probability αi,i+1(ωi,ωi+1) _i,i+1( _i, _i+1) =min1,πi(ωi+1)πi(ωi)⋅πi+1(ωi)πi+1(ωi+1) = \1,\; _i( _i+1) _i( _i)· _i+1( _i) _i+1( _i+1) \ (6) =min1,exp(Δi,i+1), = \1,\; ( _i,i+1) \, where Δi,i+1=ρi(ωi+1)−ρi(ωi)+ρi+1(ωi)−ρi+1(ωi+1). _i,i+1= _i( _i+1)- _i( _i)+ _i+1( _i)- _i+1( _i+1). This is the standard replica-exchange ratio (Remark 1), but it does not require defining or tuning temperatures. Algorithm 1 summarizes one sweep. Algorithm 1 One sweep of Gibbs sampling with Parallel Noising (PN) 1: Input: states (ω1,…,ωN)( _1,…, _N), potentials (ρ1,…,ρN)( _1,…, _N), within-level kernels (K1,…,KN)(K_1,…,K_N) 2: for i=1i=1 to N do 3: ωi∼Ki(ωi,⋅) _i K_i( _i,·) within-level update (e.g., blocked Gibbs) 4: end for 5: for phase∈even,oddphase∈\even,odd\ do 6: start←1start← 1 if phase is even else 22 7: for i=starti=start to N−1N-1 step 22 do 8: propose swap (ωi,ωi+1)←(ωi+1,ωi)( _i, _i+1)←( _i+1, _i) 9: accept with probability αi,i+1(ωi,ωi+1) _i,i+1( _i, _i+1) in (6) 10: end for 11: end for 12: Output: updated states (ω1,…,ωN)( _1,…, _N) 6.3 Training objective and estimator A central motivation for PN is to stabilize likelihood training by improving the quality of model samples used in the negative phase. The noise levels define auxiliary targets that help mixing; we fit one NMLN per level (no parameter sharing across levels) and ultimately care about the lowest-noise model. Per-level likelihood objectives. Let pdatap_data denote the empirical distribution over training worlds in Ω . For each level i, define the corrupted-data distribution qi:=pdataKνiq_i:=p_dataK_ _i (Section 4.2). We fit an NMLN at each level i with parameters θi _i and potential ρθi _ _i, inducing πθi(ω)∝exp(ρθi(ω)). _ _i(ω)\; \; ( _ _i(ω) ). Training maximizes the sum of per-level log-likelihoods maxθ1,…,θN∑i=1Nω∼qi[ρθi(ω)]−logZθi, _ _1,…, _N\; _i=1^N\;E_ω q_i [ _ _i(ω) ]\;-\; Z_ _i, (7) where Zθi=∑ω′∈Ωexp(ρθi(ω′))Z_ _i= _ω ∈ ( _ _i(ω )). The gradient for level i has the usual positive/negative-phase form ∇θiℓi=ω∼qi[∇θiρθi(ω)]−ω∼πθi[∇θiρθi(ω)], _ _i _i=E_ω q_i [ _ _i _ _i(ω) ]-E_ω _ _i [ _ _i _ _i(ω) ], where the negative-phase expectation is approximated by MCMC. Persistent PN chains for the negative phase. We maintain R persistent PN replicas (Section 6.2), each containing states (ω1(r),…,ωN(r))(ω^(r)_1,…,ω^(r)_N). At each SGD step we advance each replica by a small, fixed number T of PN sweeps (Algorithm 1), and we approximate the model expectation at level i using the current snapshot ωi(r)r=1R\ω^(r)_i\_r=1^R. We do not store or average over full chain histories; this yields a persistent stochastic gradient estimator analogous in spirit to persistent contrastive divergence (Tieleman, 2008), which maintains Markov chains across parameter updates, rather than ordinary contrastive divergence (Hinton, 2002), which reinitializes chains from data. Algorithm 2 shows one SGD step. Algorithm 2 One SGD step for training per-level NMLNs with PN 1: Input: minibatch ℬB, noise levels ν1:N _1:N, parameters θii=1N\ _i\_i=1^N, persistent states ωi(r)i,r\ω^(r)_i\_i,r, PN steps T 2: for r=1r=1 to R do 3: for t=1t=1 to T do 4: (ω1(r),…,ωN(r))←PN-Sweep((ω1(r),…,ωN(r)),ρθii=1N)(ω^(r)_1,…,ω^(r)_N)← PN-Sweep((ω^(r)_1,…,ω^(r)_N),\ _ _i\_i=1^N) Alg. 1 5: end for 6: end for 7: for i=1i=1 to N do 8: Sample corrupted minibatch ℬ~i B_i by corrupting each ω∈ℬω with KνiK_ _i 9: ℒ^i(θi)←1|ℬ~i|∑ω~∈ℬ~iρθi(ω~)−1R∑r=1Rρθi(ωi(r)) L_i( _i)← 1| B_i| _ ω∈ B_i _ _i( ω)\;-\; 1R _r=1^R _ _i(ω^(r)_i) 10: Update θi←θi+η∇θiℒ^i(θi) _i← _i+η\, _ _i L_i( _i) 11: end for 6.4 Why corruption levels are attractive in changing models During likelihood training, the target distribution changes as parameters are updated. In temperature-based parallel tempering, swap acceptance can deteriorate abruptly when the effective energy scale drifts, so temperature ladders may require frequent retuning. In PN, the ladder is defined by the fixed corruption operator KνK_ν and fixed noise levels νi\ _i\; Section 7 provides overlap and tracking guarantees that depend on the corruption increments rather than on a temperature scale. 6.5 Choosing the noising schedule in practice PN requires choosing the noise levels ν1>⋯>νN _1>·s> _N. As with all replica-exchange methods, performance hinges on overlap between adjacent targets: if πi _i and πi+1 _i+1 concentrate on disjoint regions, swaps are rarely accepted and the ladder decouples. Section 7 makes this connection precise by relating expected swap acceptance to an overlap quantity (Proposition 2) and by providing a distribution-free sanity bound for Bernoulli noising (Proposition 3). In practice, these results suggest using sufficiently fine spacing in ν so that adjacent levels retain non-trivial overlap, and validating this by monitoring swap acceptance rates. 7 Theoretical Properties of Parallel Noising This section summarizes the key guarantees that guide practical PN design; proofs are in Appendix C and Appendix D, except for Proposition 2, which is a standard replica-exchange identity that we cite. Correctness (exactness). PN is standard replica exchange on the product target Π(ω1:N)=∏i=1Nπi(ωi) ( _1:N)= _i=1^N _i( _i): within-level kernels preserve each πi _i, and Metropolis swaps preserve Π . Consequently, in stationarity the marginal of the lowest-noise chain is exactly the desired target πN _N (Proposition 1). Proposition 1 (Correctness of PN). Assume that for each level i, the within-level kernel KiK_i leaves πi _i invariant. Let K=⨂i=1NKiK= _i=1^NK_i be the product update on ΩN ^N and let Si,i+1S_i,i+1 be the Metropolis swap kernel between adjacent levels i and i+1i+1 with acceptance probability as in Eq. (6). Then: 1. K leaves the product distribution Π(ω1,…,ωN)=∏i=1Nπi(ωi) ( _1,…, _N)= _i=1^N _i( _i) invariant; 2. for each i, Si,i+1S_i,i+1 leaves Π invariant; 3. any composition of K and swap kernels (e.g., Algorithm 1) leaves Π invariant; and 4. in stationarity, the marginal distribution of the lowest-noise chain ωN _N is exactly πN _N. Proof. See Appendix D. Relation to temperature-based PT. Parallel tempering corresponds to replica exchange over a ladder of temperature-scaled targets, typically of the form πi(ω)∝exp(βiρ(ω)) _i(ω) ( _iρ(ω)) for a base energy ρ (Swendsen and Wang, 1986; Geyer, 1991; Hukushima and Nemoto, 1996). Under this temperature-scaled choice, our swap acceptance in Eq. (6) reduces to the standard PT ratio. In contrast, PN constructs intermediate targets via Bernoulli corruption in data space; this ladder is generally not equivalent to temperature scaling, so we do not claim that PT can be obtained from Bernoulli-corrupted ladders. Rather, PT and PN are two different instantiations of the same replica-exchange template with different intermediate distributions. Acceptance as distributional overlap. Expected swap acceptance between adjacent levels equals an overlap quantity: it is high exactly when adjacent targets have substantial overlap (Proposition 2). This motivates selecting νi\ _i\ so that adjacent acceptance rates remain non-trivial during training (we monitor this directly). Proposition 2 (Expected acceptance equals overlap). Let μ=πi⊗πi+1μ= _i _i+1 and μswap=πi+1⊗πiμ^swap= _i+1 _i. If (X,Y)∼μ(X,Y) μ, then [αi,i+1(X,Y)]=∫mindμ,dμswap=1−‖μ−μswap‖TV,E [ _i,i+1(X,Y) ]= \dμ,dμ^swap\=1-\|μ-μ^swap\|_TV, i.e., the mean swap acceptance equals an overlap coefficient between adjacent targets (Kofke, 2002). A simple worst-case sanity bound for Bernoulli noising. For Bernoulli corruption on d ground atoms, the total variation distance between adjacent noise levels admits a distribution-free bound proportional to d|ν−ν′|d|ν-ν | (Appendix Proposition 3), yielding a corresponding acceptance lower bound (Appendix Corollary 2). While typically loose, it clarifies that schedules with large jumps in ν can destroy overlap on high-dimensional worlds. Time-varying targets during learning. Because temperature ladders implicitly depend on an evolving energy scale, PT swap rates can degrade as parameters drift. With PN the ladder is defined by a fixed corruption operator. We show in the appendix, on a toy example, that PN may continue to perform well even in situations where PT already degrades. 8 Experiments We evaluate NMLN* (GNN potentials + PN) on four molecular benchmarks (ChEMBL, QM9, ZINC250k, MOSES). We report results only on ChEMBL and defer the full suite, ablations, and implementation details to Appendix B. Protocol and metric. Following Jung et al. (2024), we generate 2M samples per method and report recall curves: among the t most frequently generated unique molecules, how many occur in the test set (frequency as a proxy for model probability). (a) ChEMBL 10 (b) ChEMBL 15 (c) ChEMBL 20 Figure 1: Recall curves on ChEMBL: NMLN* vs. DiGress (discrete and continuous). Baselines. We compare against the strongest available NMLN baseline (DeepSet-NMLN (Jung et al., 2024)) and a general diffusion model for graphs (DiGress; discrete and continuous (Vignac et al., 2023)). Additional molecule-specific generators (MoleculeRNN variants; PaccMann) are reported in Appendix B. Figure 2: ChEMBL-10 recall curves. NMLN* (GNN + PN) improves recall over prior NMLN variants; full results and ablations are in the appendix. Q1: Can NMLN* outperform state-of-the-art diffusion models? We compare NMLN* to DiGress (Vignac et al., 2023) (continuous and discrete). Figure 1 shows that NMLN* substantially improves recall across ChEMBL sizes. Q2: How much do the improvements matter? To disentangle stronger potentials implemented as GNNs from better inference due to parallel noising, we run compute-aware ablations. Baseline (old) NMLNs. We consider NMLN s.1, the DeepSet-NMLN of Jung et al. (2024) (8 parallel chains), and NMLN s.2, which increases this to 40 chains (still practical to train). NMLN* uses N=5N=5 noise levels with R=100R=100 replicas (500 chains total) and reports samples from the coldest level. Scaling DeepSet-NMLN further (e.g., 100 chains) was prohibitively slow due to dense local potentials. Effect of the potential family (sampler fixed). Keeping plain multi-chain Gibbs sampling, NMLN s.3 (GNN + Gibbs) replaces DeepSet potentials with GNN potentials at the default chain count. Effect of the sampler (potential fixed). Fixing the GNN potential, we vary sampling strategy and compute: • NMLN s.4 (GNN + Gibbs, 100 chains) matches the number of coldest-level replicas in NMLN*. • NMLN s.5 (GNN + Gibbs, 500 chains) matches the total chain count (N×RN× R) of NMLN*. • NMLN* (GNN + PN) runs the full parallel-noising ladder (N=5N=5, R=100R=100) with swap moves. Results. Figure 2 shows that both improvements matter, and that their combination is clearly the best-performing configuration. First, replacing DeepSet potentials with GNN potentials already yields a clear improvement under the same plain Gibbs sampler (compare s.3 against the older NMLN baselines s.1 and s.2). Second, improving the sampler with parallel noising (PN) provides additional gains beyond what can be explained by simply running more Gibbs chains. In particular, s.4 and s.5 increase the number of Gibbs chains to match the coldest-level replicas and the total chain count of NMLN*, respectively, yet NMLN* (GNN + PN) still performs better. Actually, simply increasing the number of parallel chains without a replica exchange scheme seems to hurt performance. Therefore, the improvement is not a pure compute/chain-count effect: PN yields better inference than multi-chain Gibbs at matched compute. 9 Conclusions We introduced NMLN*, a strengthened Neural Markov Logic Network that tackles two practical bottlenecks: limited expressivity of fragment-based potentials and brittle, slow-mixing Gibbs inference under near-deterministic learned constraints. On the modeling side, we instantiate NMLN potentials with global GNN energies to capture long-range structure. On the inference side, we propose parallel noising, an exact replica-exchange sampler that couples chains across a ladder of fixed corruption levels rather than temperatures, which stays well-behaved even as the target distribution changes during training. Empirically, NMLN* improves recall and sample quality on molecular generation benchmarks. References R. Abboud, I. Ceylan, and T. Lukasiewicz (2020) Learning to reason: leveraging neural networks for approximate DNF counting. In Proceedings of the AAAI Conference on Artificial Intelligence, Vol. 34, p. 3097–3104. External Links: Document Cited by: §2. K. Ahmed, S. Teso, K. Chang, G. Van den Broeck, and A. Vergari (2022) Semantic probabilistic layers for neuro-symbolic learning. Advances in Neural Information Processing Systems 35, p. 29944–29959. Cited by: §2. S. Badreddine, A. d. Garcez, L. Serafini, and M. Spranger (2022) Logic tensor networks. Artificial Intelligence 303, p. 103649. Cited by: §2. P. Beame, G. Van den Broeck, E. Gribkoff, and D. Suciu (2015) Symmetric weighted first-order model counting. In Proceedings of the 34th ACM SIGMOD-SIGACT-SIGAI Symposium on Principles of Database Systems, p. 313–328. Cited by: footnote 2. J. Born, M. Manica, A. Oskooei, J. Cadow, G. Markert, and M. Rodríguez Martínez (2021) PaccMannRL: de novo generation of hit-like anticancer molecules from transcriptomic data via reinforcement learning. iScience 24 (4), p. 102269. External Links: Document, Link Cited by: 4th item. G. Casella and R. L. Berger (2024) Statistical inference. 2 edition, Chapman and Hall/CRC. External Links: ISBN 9781032593036 Cited by: §4.1. L. De Smet, E. Sansone, and P. Zuidberg Dos Martires (2023) Differentiable sampling of categorical distributions using the CatLog-derivative trick. Advances in Neural Information Processing Systems 36, p. 30416–30428. Cited by: §2. L. De Smet, G. Venturato, L. De Raedt, and G. Marra (2025) Relational neurosymbolic markov models. In Proceedings of the AAAI Conference on Artificial Intelligence, Vol. 39, p. 16181–16189. Cited by: §2, §2. G. Desjardins, A. Courville, Y. Bengio, P. Vincent, and O. Delalleau (2010) Tempered markov chain monte carlo for training of restricted boltzmann machines. In Proceedings of the thirteenth international conference on artificial intelligence and statistics, p. 145–152. Cited by: §A.1, §2. L. Di Liello, P. Ardino, J. Gobbi, P. Morettin, S. Teso, and A. Passerini (2020) Efficient generation of structured objects with constrained adversarial networks. Advances in neural information processing systems 33, p. 14663–14674. Cited by: §2. M. Diligenti, M. Gori, and C. Saccà (2017) Semantic-based regularization for learning and inference. Artificial Intelligence 244, p. 143–165. Cited by: §2. C. J. Geyer (1991) Markov chain monte carlo maximum likelihood. In Proceedings of the 23rd Symposium on the Interface: Computing Science and Statistics, p. 156–163. Cited by: §A.1, §1, §2, §6.2, §7. V. Gutiérrez-Basulto, J. C. Jung, and O. Kuzelka (2018) Quantified markov logic networks. In Principles of Knowledge Representation and Reasoning: Proceedings of the Sixteenth International Conference, KR 2018, M. Thielscher, F. Toni, and F. Wolter (Eds.), p. 602–612. External Links: Link Cited by: §5.1. G. E. Hinton (2002) Training products of experts by minimizing contrastive divergence. Neural Computation 14 (8), p. 1771–1800. External Links: Document Cited by: §6.3. Y. Huang, A. Ghatare, Y. Liu, Z. Hu, Q. Zhang, C. Shama Sastry, S. Gururani, S. Oore, and Y. Yue (2024) Symbolic music generation with non-differentiable rule guided diffusion. In Proceedings of the 41st International Conference on Machine Learning, Proceedings of Machine Learning Research, Vol. 235, p. 19772–19797. External Links: Link Cited by: §2. K. Hukushima and K. Nemoto (1996) Exchange monte carlo method and application to spin glass simulations. Journal of the Physical Society of Japan 65 (6), p. 1604–1608. Cited by: §A.1, §1, §2, §6.2, §7. J. J. Irwin and B. K. Shoichet (2005) ZINC—a free database of commercially available compounds for virtual screening. Journal of Chemical Information and Modeling 45 (1), p. 177–182. Note: PMID: 15667143 External Links: Document, Link Cited by: 3rd item. P. Jung, G. Marra, and O. Kuželka (2024) Quantified neural markov logic networks. International Journal of Approximate Reasoning 171, p. 109172. Cited by: §B.3, §1, §2, §5.1, §5.1, §5.1, §8, §8, §8. D. A. Kofke (2002) On the acceptance probability of replica-exchange monte carlo trials. The Journal of Chemical Physics 117 (15), p. 6911–6914. External Links: Document Cited by: Appendix D, Proposition 2. D. A. Kofke (2004) Erratum: “on the acceptance probability of replica-exchange monte carlo trials” [J. chem. phys. 117, 6911 (2002)]. The Journal of Chemical Physics 120 (22), p. 10852. External Links: Document Cited by: Appendix D. M. Krenn, F. Häse, A. K. Nigam, P. Friederich, and A. Aspuru-Guzik (2020) Self-referencing embedded strings (SELFIES): a 100% robust molecular string representation. Machine Learning: Science and Technology 1 (4), p. 045024. External Links: Document, Link Cited by: 3rd item, §B.1. Z. Li, J. Huang, and M. Naik (2023) Scallop: a language for neurosymbolic programming. Proceedings of the ACM on Programming Languages 7 (PLDI), p. 1463–1487. Cited by: §2, §2. R. Manhaeve, S. Dumančić, A. Kimmig, T. Demeester, and L. De Raedt (2018) DeepProbLog: neural probabilistic logic programming. Advances in Neural Information Processing Systems 31. Cited by: §2, §2. R. Manhaeve, G. Marra, and L. De Raedt (2021) Approximate inference for neural probabilistic logic programming. In 18th International Conference on Principles of Knowledge Representation and Reasoning (KR), p. 475–486. Cited by: §2. G. Marra, M. Diligenti, F. Giannini, M. Gori, and M. Maggini (2020) Relational neural machines. In ECAI 2020, p. 1340–1347. Cited by: §2. G. Marra and O. Kuželka (2021) Neural markov logic networks. In Proceedings of the Thirty-Seventh Conference on Uncertainty in Artificial Intelligence, Proceedings of Machine Learning Research, Vol. 161, p. 908–917. External Links: Link Cited by: §1, §2, §2, §2, §4.1, §4.1, §5.1, §5.2, footnote 3. B. Miasojedow, E. Moulines, and M. Vihola (2013) An adaptive parallel tempering algorithm. Journal of Computational and Graphical Statistics 22 (3), p. 649–664. External Links: Document Cited by: §2. E. Misino, G. Marra, and E. Sansone (2022) VAEL: bridging variational autoencoders and probabilistic logic programming. Advances in Neural Information Processing Systems 35, p. 4667–4679. Cited by: §2. PaccMann contributors (2021) paccmann_rl: pipeline for the PaccMannRL study. Note: https://github.com/PaccMann/paccmann\_rlGitHub repository Cited by: 4th item. D. Polykovskiy, A. Zhebrak, B. Sanchez-Lengeling, S. Golovanov, O. Tatanov, S. Belyaev, R. Kurbanov, A. Artamonov, V. Aladinskiy, M. Veselov, A. Kadurin, S. Johansson, H. Chen, S. Nikolenko, A. Aspuru-Guzik, and A. Zhavoronkov (2020) Molecular Sets (MOSES): A Benchmarking Platform for Molecular Generation Models. Frontiers in Pharmacology 11, p. 565644. External Links: Document, Link Cited by: 4th item. H. Poon and P. M. Domingos (2006) Sound and efficient inference with probabilistic and deterministic dependencies. In Proceedings, The Twenty-First National Conference on Artificial Intelligence and the Eighteenth Innovative Applications of Artificial Intelligence Conference, p. 458–463. Cited by: §5.2. C. Pryor, C. Dickens, E. Augustine, A. Albalak, W. Y. Wang, and L. Getoor (2023) NeuPSL: neural probabilistic soft logic. In Proceedings of the Thirty-Second International Joint Conference on Artificial Intelligence, p. 4145–4153. External Links: Document Cited by: §2. R. Ramakrishnan, P. O. Dral, M. Rupp, and O. A. von Lilienfeld (2014) Quantum chemistry structures and properties of 134 kilo molecules. Scientific Data 1, p. 140022. External Links: Document, Link Cited by: 2nd item. M. Richardson and P. Domingos (2006) Markov logic networks. Machine learning 62 (1), p. 107–136. Cited by: §1, §2, §5.1. C. P. Robert and G. Casella (1999) Monte carlo statistical methods. Springer Texts in Statistics, Springer. External Links: Document, Link Cited by: §2. L. Ruddigkeit, R. van Deursen, L. C. Blum, and J. Reymond (2012) Enumeration of 166 billion organic small molecules in the chemical universe database GDB-17. Journal of Chemical Information and Modeling 52 (11), p. 2864–2875. Note: PMID: 23088335 External Links: Document, Link Cited by: 2nd item. M. S. M. Sajjadi, O. Bachem, M. Lucic, O. Bousquet, and S. Gelly (2018) Assessing generative models via precision and recall. Advances in Neural Information Processing Systems 31. Cited by: §B.3. D. Scassola, S. Saccani, G. Carbone, and L. Bortolussi (2023) Zero-shot conditioning of score-based diffusion models by neuro-symbolic constraints. arXiv preprint arXiv:2308.16534. Cited by: §2. M. Schlichtkrull, T. N. Kipf, P. Bloem, R. van den Berg, I. Titov, and M. Welling (2018) Modeling relational data with graph convolutional networks. In The Semantic Web, Lecture Notes in Computer Science, Vol. 10843, p. 593–607. External Links: Document Cited by: §6.1. W. Shi (2025) Molecule-RNN: a recurrent neural network that generates drug-like molecules for drug discovery. Note: https://github.com/shiwentao00/Molecule-RNNGitHub repository; accessed 2025-04-15 Cited by: 1st item. M. C. Stoian, S. Dyrmishi, M. Cordy, T. Lukasiewicz, and E. Giunchiglia (2024) How realistic is your synthetic data? constraining deep generative models for tabular data. In The Twelfth International Conference on Learning Representations, External Links: Link Cited by: §2. R. H. Swendsen and J. Wang (1986) Replica monte carlo simulation of spin-glasses. Physical Review Letters 57 (21), p. 2607–2609. Cited by: §A.1, §1, §2, §6.2, §7. S. Syed, A. Bouchard-Côté, G. Deligiannidis, and A. Doucet (2022) Non-reversible parallel tempering: a scalable highly parallel mcmc scheme. Journal of the Royal Statistical Society Series B: Statistical Methodology 84 (2), p. 321–350. Cited by: §A.1, §2. T. Tieleman (2008) Training restricted boltzmann machines using approximations to the likelihood gradient. In Proceedings of the 25th International Conference on Machine Learning, p. 1064–1071. External Links: Document, Link Cited by: §6.3. E. Van Krieken, E. Acar, and F. Van Harmelen (2022) Analyzing differentiable fuzzy logic operators. Artificial Intelligence 302, p. 103602. Cited by: §2. E. van Krieken, P. Minervini, E. Ponti, and A. Vergari (2025) Neurosymbolic diffusion models. In Advances in Neural Information Processing Systems, Vol. 38. External Links: 2505.13138, Link Cited by: §2, §2. E. van Krieken, T. Thanapalasingam, J. Tomczak, F. Van Harmelen, and A. Ten Teije (2023) A-NeSI: a scalable approximate method for probabilistic neurosymbolic inference. Advances in Neural Information Processing Systems 36, p. 24586–24609. Cited by: §2. V. Verreet, L. De Smet, L. De Raedt, and E. Sansone (2024) EXPLAIN, AGREE, LEARN: scaling learning for neural probabilistic logic. In ECAI 2024: 27th European Conference on Artificial Intelligence, Frontiers in Artificial Intelligence and Applications, Vol. 392, p. 1349–1356. External Links: 2408.08133, Link Cited by: §2. C. Vignac and P. Frossard (2023) DiGress: discrete denoising diffusion for graph generation. Note: https://github.com/cvignac/DiGressGitHub repository; accessed 2025-04-17 Cited by: 1st item. C. Vignac, I. Krawczuk, A. Siraudin, B. Wang, V. Cevher, and P. Frossard (2023) DiGress: discrete denoising diffusion for graph generation. In The Eleventh International Conference on Learning Representations, External Links: Link Cited by: 1st item, §B.1, §B.5, §2, §8, §8. T. Winters, G. Marra, R. Manhaeve, and L. De Raedt (2022) DeepStochLog: neural stochastic logic programming. In Proceedings of the AAAI Conference on Artificial Intelligence, Vol. 36, p. 10090–10100. Cited by: §2. J. Xu, Z. Zhang, T. Friedman, Y. Liang, and G. Van den Broeck (2018) A semantic loss function for deep learning with symbolic knowledge. In Proceedings of the 35th International Conference on Machine Learning, Proceedings of Machine Learning Research, Vol. 80, p. 5502–5511. External Links: Link Cited by: §2. B. Zdrazil, E. Felix, F. Hunter, E. J. Manners, J. Blackshaw, S. Corbett, M. de Veij, H. Ioannidis, D. Mendez Lopez, J. F. Mosquera, M. P. Magarinos, N. Bosc, R. Arcila, T. Kizilören, A. Gaulton, A. P. Bento, M. F. Adasme, P. Monecke, G. A. Landrum, and A. R. Leach (2024) The ChEMBL database in 2023: a drug discovery platform spanning multiple bioactivity data types and time periods. Nucleic Acids Research 52 (D1), p. D1180–D1192. External Links: ISSN 0305-1048, Document, Link Cited by: 1st item. Appendix Appendix A Parallel Noising versus Parallel Tempering A.1 Parallel Tempering MCMC Parallel tempering MCMC [Swendsen and Wang, 1986, Geyer, 1991, Hukushima and Nemoto, 1996, Desjardins et al., 2010], which is also known as replica exchange MCMC, is a classical method for sampling from complex distributions where plain Gibbs sampling would be inefficient. It improves mixing by simulating several Markov chains in parallel, which may exchange samples (explained below). Each chain is associated with an inverse temperature parameter βi∈(0,1] _i∈(0,1], where the chain with βk=1 _k=1 targets the true distribution we want to sample from—lower βi _i’s correspond to flatter versions of this distribution. Formally, each chain samples from the tempered distribution Pi(x)=1Ziexp(−βiE(x)),P_i(x)= 1Z_i (- _iE(x)), where E(x)E(x) is the energy function of the target distribution and ZiZ_i is the normalization constant, known as partition function. The idea is that chains at lower inverse temperatures (i.e., higher temperatures) are better at exploring the space due to the flattened energy landscape, while the chain at βk=1 _k=1 provides samples from the true target. To allow information sharing between chains, the algorithm performs swap attempts between adjacent chains x(i)x^(i) and x(i+1)x^(i+1). The proposed swap (x(i),x(i+1))↦(x(i+1),x(i))(x^(i),x^(i+1)) (x^(i+1),x^(i)) is accepted with probability α=min1,Pi(x(i+1))⋅Pi+1(x(i))Pi(x(i))⋅Pi+1(x(i+1))=min1,exp(−βiE(x(i+1))+βiE(x(i))−βi+1E(x(i))+βi+1E(x(i+1)))α= \1,\; P_i(x^(i+1))· P_i+1(x^(i))P_i(x^(i))· P_i+1(x^(i+1)) \\\ = \1,\; (- _iE(x^(i+1))+ _iE(x^(i)) .\\ .- _i+1E(x^(i))+ _i+1E(x^(i+1)) ) \ which ensures that each chain samples from the correct corresponding distribution. It turns out that it is also beneficial to perform the swapping in “odd” and “even” phases, where in the “odd” phases, the swaps are attempted for chains with odd number i and their neighbor i+1i+1, and analogically for even phases [Syed et al., 2022]. The motivation behind this technique is that chains at higher temperatures are more likely to traverse low-probability regions and can thus escape local modes. The swap mechanism then allows this exploratory information to be propagated back to the lower-temperature chains. Example 1. Figure 3 shows, as an illustration, results of sampling from a given bimodal density with and without parallel tempering—as can be seen, parallel tempering helps the sampler to cross the region of low density separating the two modes. Figure 3: Sampling from a bi-modal distribution: No PT vs. PT Remark 1. It will be convenient for the exposition further in the paper to make the substitution ρi(x)=−βiE(x) _i(x)=- _iE(x). Here, we just absorbed the inverse temperature parameter and the energy function into the potential function ρi _i. The swap probability in parallel tempering then becomes α=min1,exp(ρi(x(i+1))−ρi(x(i))+ρi+1(x(i))−ρi+1(x(i+1)))α= \1,\; ( _i(x^(i+1))- _i(x^(i))+ _i+1(x^(i))- _i+1(x^(i+1)) )\. A.2 Intuition: Parallel Noising vs. Parallel Tempering Parallel noising (PN) and parallel tempering (PT) are both replica-exchange methods: they run a ladder of auxiliary chains and use swap moves to transfer exploratory states back toward the target distribution. The crucial difference is how the auxiliary distributions are constructed. PT builds its ladder by rescaling the energy with different temperatures. High-temperature chains flatten the energy landscape and can therefore move more easily between modes. This works well when the temperature ladder is well matched to the energy scale of the target: adjacent chains must have enough distributional overlap for swap moves to be accepted. However, this requirement can become problematic when the target is very sharp, highly constrained, or changing during training. In such cases, a fixed temperature ladder may either be too coarse, leading to low swap acceptance, or too mild, giving hot chains that are not exploratory enough. PN instead builds its ladder by corrupting the state space. Rather than changing the energy scale, it defines auxiliary chains corresponding to increasingly noised versions of the structured object. The higher-noise chains explore relaxed versions of the same combinatorial problem, while the lowest-noise chain remains the target of interest. This can be more robust when the difficulty comes from sharp symbolic or combinatorial constraints, because the ladder is tied to the structure of the data rather than to a temperature scale that may be hard to tune. The next two subsections illustrate this distinction in controlled settings. The first diagnostic uses a synthetic Boolean distribution whose energy scale changes over time, showing how PT can become sensitive to the temperature schedule. The second diagnostic uses a pixel-based n-queens generation task with CNN potentials, where the goal is to sample configurations with few or no queen attacks. A.3 Time-varying Boolean tracking example This toy example isolates the issue of changing energy scales during learning. It shows how a fixed temperature ladder can become poorly matched to a target distributions whose low energy barriers grow over time due to learning from deterministic distributions. Target. Let x∈0,1dx∈\0,1\^d and let k=‖x‖1k=\|x\|_1 be its Hamming weight. At time t we define πt(x) _t(x) ∝exp−Et(x), \-E_t(x)\, Et(x) E_t(x) ≡Et(k)=βtmink,d−k−ht(2k−d). ≡ E_t(k)= _t \k,d-k\-h_t(2k-d). The first term creates two wells around 0d0^d and 1d1^d; increasing βt _t raises the barrier between them. The second term tilts the distribution toward one well, and we sweep hth_t from negative to positive so that the preferred well changes during the run. We track the magnetization m(x)=1d∑j=1d(2xj−1)=2k−d∈[−1,1],m(x)= 1d _j=1^d(2x_j-1)= 2k-dd∈[-1,1], whose exact expectation can be computed by summing over Hamming weights with multiplicities (dk) dk. PT and PN configurations. For PT, the replicas target πt,Tℓ(x)∝exp−Et(x)/Tℓ _t,T_ (x) \-E_t(x)/T_ \ using a geometrically spaced temperature ladder Tℓ=Tmaxℓ/(L−1)T_ =T_ /(L-1), with Tmax=4T_ =4 and L=10L=10. For PN, the replicas target corrupted distributions πt,ν=πtKν _t,ν= _tK_ν, where KνK_ν is the independent bit-flip channel. We use the quadratic noise schedule νℓ=12(ℓ/(Lν−1))2 _ = 12( /(L_ν-1))^2, which is denser near the clean target. Both samplers use local bit-flip moves and odd–even adjacent swaps. Result. Figure 4 shows the cold-replica estimate of πt[m]E_ _t[m], averaged across inner-chain samples and random seeds. As the barrier increases, PT with a fixed temperature range increasingly lags behind the changing target. PN tracks the change more closely because high-noise replicas remain close to uniform even when the underlying energy becomes sharp. This supports the use of corruption levels as a robust auxiliary ladder when training gradually learns near-deterministic constraints. Figure 4: Tracking a time-varying Boolean two-well target. Ground truth πt[m]E_ _t[m] is computed exactly by summing over Hamming weights. PT uses a fixed temperature ladder, while PN uses a fixed noising ladder. Curves for PT and PN average the cold-replica magnetization over inner-chain samples and multiple random seeds. A.4 Pixel-based n-queens generation We also evaluate PN and PT on a pixel-based n-queens generation task. A board is represented as a binary image, with a queen indicated by an active pixel. The energy is parameterized by a convolutional neural network, so the task tests whether the sampler can explore a learned, image-level potential whose high-probability states satisfy a global combinatorial constraint. Metric. For each generated board, we count the number of attacking queen pairs, i.e., pairs of active pixels that share a row, column, or diagonal. We report the average number of attacks over generated samples; lower is better, and zero attacks corresponds to a solved board when the board also contains the required number of queens. This metric directly measures how much hard constraint violation remains in the generated samples. Comparison. We use the same CNN potential and the same within-level local updates for both methods. PT uses a temperature ladder beginning at the clean target (T=1T=1) and increasing through 1.0,1.08,1.1664,1.259712,1.360489\1.0,1.08,1.1664,1.259712,1.360489\. PN uses the corruption ladder ν∈0,0.02,0.05,0.1,0.2ν∈\0,0.02,0.05,0.1,0.2\. In both cases, we report the clean/cold chain, and the curves correspond to independent runs. Result. Figure 5 shows that both methods rapidly reduce the number of attacks from the random initialization, but their long-run behavior differs. PT plateaus at a higher number of remaining conflicts, whereas PN continues to reduce the average number of attacks across the run. The final samples in Figure 6 give a qualitative view of the same effect: high-noise PN levels remain exploratory, while the clean PN chain produces less conflicted boards than the corresponding cold PT chain. This diagnostic supports the main claim that noising ladders can be better suited than fixed temperature ladders for learned energies with crisp combinatorial structure. Figure 5: Pixel-based n-queens diagnostic. Average number of attacking queen pairs as a function of generated samples. Lower is better. PN reaches consistently lower attack counts than PT after the initial burn-in, indicating better exploration of low-conflict boards under the learned CNN potential. (a) PT levels. (b) PN levels. Figure 6: Representative final n-queens samples at different temperature/noise levels. The cold/clean level is shown at the top of each panel. Appendix B Additional experimental results and ablations B.1 Compared Models Besides the latest version of NMLN, we compare NMLN* with the following two diffusion-based algorithms from [Vignac et al., 2023]: • DiGress Discrete [Vignac et al., 2023], implemented in [Vignac and Frossard, 2023], utilizes a discrete diffusion process that progressively edits graphs with noise, through the process of adding or removing edges and changing the labels. • DiGress Continuous which is the continuous version using Gaussian noise. We also compare with molecule-specific generators that use string encodings such as SMILES or SELFIES [Krenn et al., 2020]. SELFIES guarantees molecular validity for strings over its alphabet, whereas SMILES does not. These methods exploit chemistry-specific representations and therefore do not apply to arbitrary relational domains, but they provide informative domain-specialized baselines: • Molecule-RNN Char [Shi, 2025] is a recurrent neural network designed to generate novel molecules from the distribution of a molecular training dataset. This method is based on the SMILES representation. • Molecule-RNN Regex A more structured approach using regular expressions, where elements such as square-bracketed atoms (e.g., ”[O-]”) are treated as individual tokens. • Molecule-RNN Selfies [Krenn et al., 2020] uses the Selfies representation of molecules—this model therefore generates only valid molecules. • PaccMann [Born et al., 2021], implemented by the PaccMannRL pipeline [PaccMann contributors, 2021], is a VAE-based generative model guided by reinforcement learning, where a learned reward function predicts molecule–target interactions, enabling the targeted generation of novel compounds for given protein or transcriptomic inputs. B.2 Datasets We report experimental results on four datasets restricted to molecules of sizes 9, 10, 15, and 20. These datasets are: • ChEMBL [Zdrazil et al., 2024], which is a manually curated database of bioactive molecules with drug-like properties, • QM9 [Ruddigkeit et al., 2012, Ramakrishnan et al., 2014] contains stable small organic molecules made up of carbon, hydrogen, oxygen, nitrogen, and fluorine atoms, • ZINC250k [Irwin and Shoichet, 2005], which is a free database of commercially available compounds for virtual screening, and • MOSES [Polykovskiy et al., 2020], which is a curated and cleaned subset of the ZINC database, prepared specifically for benchmarking molecular generative models. The number of molecules in the subsets of these datasets which we use are reported in Table 1. Dataset Size 9 Size 10 Size 15 Size 20 ChEMBL 1 684 2 589 14 658 44 222 QM9 109 813 - - - MOSES - - 579 195 650 ZINC250k 28 90 1 766 8 590 Table 1: Number of molecules in the datasets with a given number of heavy atoms. B.3 Methodology We mainly follow the experimental methodology from [Jung et al., 2024] where it was advocated to use a variant of ROC curves, called coverage or recall curves (we use the latter term in this paper), and motivated by Sajjadi et al. [2018], obtained as follows: We let the model generate N samples and we collect all the unique ones. For each unique sample, we store its frequency, i.e., how many times they were generated. Here, frequency is understood as a proxy for the probability of the molecule given by the learned model. As a generative performance indicator, we compute how many of the t most frequently generated molecules are in the test set. We then plot this measure w.r.t. frequency thresholds t. This metric is related to a ROC curve, but in the generative setting. To ensure a fair comparison, we let all methods generate a fixed number of 2M samples per model and dataset. While some methods may require longer generation times, this can be considered negligible in the context of real-world molecular synthesis where the top-ranking molecules suggested by the models may need to eventually be synthesized and tested. B.4 Detailed Experimental Settings Datasets are partitioned into training and test sets using an 80/20 split. All models were trained using their default hyperparameters as provided in the respective open-source repositories. The only modifications involved adjusting the generation pipeline to produce 2,000,000 samples, saving outputs in the SMILES format, and, in some cases, disabling third-party post-processing filters that would otherwise exclude invalid samples. The GNN potentials in NMLN* consist of 10 R-GCN layers with a hidden dimensionality of 128 and ReLU activations. The update function is implemented as a standard linear layer with a matching hidden size and activation. The output of the GNN is formed by aggregating all intermediate layer outputs and applying a final MLP with three linear layers, each using ReLU activations and the same hidden size. B.5 Experimental Questions and Results In this section, we state several research questions and address them experimentally. Q1: Does NMLN* outperform the latest NMLN? Figure 7: Recall curves for different variants of NMLN on the dataset ChEMBL 10. We compared NMLN* with two baselines derived from the previous model: NMLN s.1, which uses the DeepSet-NMLN settings, and NMLN s.2, which increases the number of parallel chains from 8 in the DeepSet-NMLN to 40, which is still trainable in a reasonable time. In NMLN*, we run R=100R=100 replicas over N=5N=5 noise levels (500 chains total); we report samples from the coldest level (setting it to 100 for the DeepSet-NMLN already lead to unacceptably high runtime, which is a consequence of the dense local potential functions used in the original NMLNs). As shown in Figure 7, the NMLN* models outperform both NMLN competitors already for molecules with 10 atoms by a significant margin, demonstrating the effectiveness of the proposed modifications. Q2: How much do the improvements matter? To separate the two contributions—stronger potentials and better inference—we consider a set of compute-aware ablations, all summarized in Figure 7. Effect of the potential family (sampler fixed). We compare the original DeepSet potentials (NMLN s.1/s.2 from Q1) against GNN potentials while keeping the sampler as plain multi-chain Gibbs: • NMLN s.3 (GNN + Gibbs) replaces the DeepSet potential with the GNN potential but keeps the default chain count. Effect of the sampler (potential fixed). We then keep the GNN potential fixed and vary the sampling strategy and compute budget: • NMLN s.4 (GNN + Gibbs, 100 chains) increases the number of Gibbs chains to match the number of replicas at the coldest level of NMLN*. • NMLN s.5 (GNN + Gibbs, 500 chains) matches the total number of maintained chains (N×RN× R) in NMLN*. • NMLN* (GNN + PN) uses the full parallel-noising ladder with N=5N=5 levels and R=100R=100 replicas (500 chains total), including swap moves. Compute note. Replica exchange incurs overhead (swap evaluations) beyond within-level updates. To avoid overstating improvements, we report results both at matched chain count (s.4) and at matched total chain count (s.5), and we use the same generation budget (2M produced samples) across all methods. To keep within-level updates tractable with a global GNN energy, we batch many candidate local edits on the GPU and reuse intermediate message-passing activations whenever possible; for the molecule sizes considered here (up to 20 heavy atoms), a full forward recomputation per proposed edit is also feasible and remains a small fraction of overall training time. As shown in Figure 7, each individual modification contributes to performance improvements, but the full NMLN* configuration achieves the highest performance overall. Q3: Can NMLN* outperform state-of-the-art diffusion models? To answer this question, we compared NMLN* with DiGress [Vignac et al., 2023], both in the continuous and in the discrete versions. DiGress is another general graph-generation method and therefore a direct competitor to NMLNs. In Appendix Figure 8, 9, 10, and 11, we compare the recall curves. On most of the datasets, NMLN* significantly outperforms diffusion models. On size 20, the recall is small for both models as the number of possible molecules is huge, but NMLN* still achieves a higher recall. We also observed that NMLN* was able to generate a significantly larger number of valid samples on the larger molecules. (a) ChEMBL 9 (b) ChEMBL 10 (c) ChEMBL 15 (d) ChEMBL 20 Figure 8: Recall curves for different sizes on the dataset ChEMBL. (a) MOSES15 (b) MOSES20 Figure 9: Recall curves for different sizes on the dataset MOSES. Figure 10: Recall curves for the QM9. (a) ZINC250k 9 (b) ZINC250k 10 (c) ZINC250k 15 (d) ZINC250k 20 Figure 11: Recall curves for different sizes on the dataset ZINC250k. Q4: How does NMLN* compare to specialized models for molecules? We further compared NMLN* against molecule-specific generative methods, which leverage domain-specific representations such as SMILES or SELFIES. These representations inherently encode chemical rules, enabling such models to achieve higher validity and more efficient exploration of molecular space. In Appendix Figure 12, 13, 14, and 15 we show that NMLN* outperforms molecule-specific baselines on smaller molecule sizes (9 and 10 heavy atoms), which shows their effectiveness even without relying on handcrafted molecular priors. However, at sizes 15 and 20, where both NMLN* and DiGresses begin to exhibit performance degradation, molecule-specific methods continue to cover a substantial portion of the test set. This suggests that the inductive biases encoded in molecular representations become increasingly beneficial as molecular complexity grows. (a) ChEMBL 9 (b) ChEMBL 10 (c) ChEMBL 15 (d) ChEMBL 20 Figure 12: Recall curves for different sizes on the dataset ChEMBL. (a) MOSES15 (b) MOSES20 Figure 13: Recall curves for different sizes on the dataset MOSES. Figure 14: Recall curves for the QM9. (a) ZINC250k 9 (b) ZINC250k 10 (c) ZINC250k 15 (d) ZINC250k 20 Figure 15: Recall curves for different sizes on the dataset ZINC250k. A Peek Into the Chains of NMLNs* Validity Across Noise Levels Intuitively, chains in parallel noising with less noise should produce samples that are valid molecules more often than chains with more noise. We measured the fraction of valid molecules produced at each of the five NMLN* noise levels in parallel noising (Fig. 17). Validity was near zero at the highest noise level (0.1) and increased monotonically as noise decreased, with the lowest level (0.001) yielding almost exclusively valid samples. Recall Across Noise Levels Figure 16 illustrates the recall curves of the five levels of single NMLN* from a single experiment, corresponding to progressively lower noise levels: 0.1, 0.01, 0.005, 0.0025, and 0.001. As expected, the first level (noise = 0.1) demonstrates limited recall due to its high stochasticity and exploratory behavior. In contrast, the final level (noise = 0.001) achieves the highest recall. Figure 16: Recall curve of NMLN* on different levels on the dataset ChEMBL. Figure 17: Validity of NMLN* on different levels on the dataset ChEMBL 10. Appendix C Full Theory for Parallel Noising This section develops properties of parallel noising (PN) that are useful both conceptually and for designing practical ladders. Throughout, Ω is a finite state space of possible worlds (for a fixed domain size) and each level i targets πi(ω)∝exp(ρi(ω)) _i(ω) ( _i(ω)) as in Section 6. C.1 Swap acceptance and distributional overlap Let μ=πi⊗πi+1μ= _i _i+1 and μswap=πi+1⊗πiμ^swap= _i+1 _i be distributions on Ω2 ^2. Denote total variation distance by ∥⋅∥TV\|·\|_TV. Corollary 1 (A simple lower bound). [αi,i+1(X,Y)]≥1−2‖πi−πi+1‖TV.E [ _i,i+1(X,Y) ]≥ 1-2\| _i- _i+1\|_TV. Proof. Write δ=πi−πi+1δ= _i- _i+1. Then μ−μswap μ-μ^swap =πi⊗πi+1−πi+1⊗πi = _i _i+1- _i+1 _i =δ⊗πi+1−πi+1⊗δ. =δ _i+1- _i+1 δ. Using ‖A‖TV=12‖A‖1\|A\|_TV= 12\|A\|_1, the triangle inequality, and ‖δ⊗π‖1=‖δ‖1‖π‖1=‖δ‖1\|δ π\|_1=\|δ\|_1\|π\|_1=\|δ\|_1, ‖μ−μswap‖TV \|μ-μ^swap\|_TV ≤12(‖δ‖1+‖δ‖1) ≤ 12(\|δ\|_1+\|δ\|_1) =‖δ‖1 =\|δ\|_1 =2‖πi−πi+1‖TV. =2\| _i- _i+1\|_TV. Apply Proposition 2. ∎ C.2 A distribution-free bound for Bernoulli noising We now specialize to the Bernoulli corruption operator KνK_ν from Section 4.2. For any distribution p on 0,1d\0,1\^d, let pν=pKνp_ν=pK_ν denote the ν-noised distribution. Proposition 3 (Total variation between adjacent noise levels). For any p and any ν,ν′∈[0,1]ν,ν ∈[0,1], ‖pν−pν′‖TV≤d|ν−ν′|.\|p_ν-p_ν \|_TV≤ d\,|ν-ν |. Interpretation. Proposition 3 is a distribution-free (worst-case) sanity check: it depends on the ambient dimension d of the ground-atom encoding and can therefore be extremely loose for realistic relational or molecular representations (where d may be in the thousands or more). In particular, the implied acceptance lower bound in Corollary 2 becomes vacuous once d|νi−νi+1|d| _i- _i+1| is not small. We therefore do not use this bound to tune noise schedules. Instead, in practice we select and validate ladders using the overlap/acceptance characterization from Section C.1: we monitor adjacent swap acceptance rates during training and choose νi\ _i\ such that adjacent pairs maintain non-trivial acceptance. We also report recall/quality as a function of compute (number of energy evaluations) when comparing different ladders or numbers of replicas (Section 8). Proof. Construct a coupling (Y,Y′)(Y,Y ) with Y∼pνY p_ν and Y′∼pν′Y p_ν as follows: sample X∼pX p, sample i.i.d. U1,…,Ud∼Unif(0,1)U_1,…,U_d (0,1), set Bj=Uj<νB_j=1\U_j<ν\ and Bj′=Uj<ν′B _j=1\U_j<ν \, and define Yj=Xj⊕BjY_j=X_j B_j and Yj′=Xj⊕Bj′Y _j=X_j B _j. Then ℙ(Y≠Y′)≤∑j=1dℙ(Bj≠Bj′)=d|ν−ν′|P(Y≠ Y )≤ _j=1^dP(B_j≠ B _j)=d|ν-ν |. By the coupling characterization of total variation, ‖pν−pν′‖TV≤ℙ(Y≠Y′)≤d|ν−ν′|\|p_ν-p_ν \|_TV (Y≠ Y )≤ d|ν-ν |. ∎ Corollary 2 (Lower bound on expected swap acceptance). Assume (idealized) that two adjacent PN levels target πi=pνi _i=p_ _i and πi+1=pνi+1 _i+1=p_ _i+1 for a common underlying distribution p. Then in stationarity, [αi,i+1(X,Y)]≥1−2d|νi−νi+1|.E [ _i,i+1(X,Y) ]≥ 1-2d| _i- _i+1|. This worst-case bound is informative only when 2d|νi−νi+1|≪12d| _i- _i+1| 1; otherwise the right-hand side becomes vacuous (and may be negative). Proof. Combine Corollary 1 with Proposition 3. ∎ C.3 Invariance under potential shifts Proposition 4 (Shift-invariance of swap decisions). Replacing ρi _i by ρ~i=ρi+ci ρ_i= _i+c_i for any constants ci∈ℝc_i leaves the swap acceptance probability (6) unchanged. Proof. In Δi,i+1 _i,i+1, each ρi _i appears once with a plus and once with a minus sign, so additive constants cancel. ∎ Appendix D Proofs Proof for Proposition 1. Proof. Let Ω be the (discrete) single-replica state space and ΩN ^N the joint state space. For each noise level j∈1,…,Nj∈\1,…,N\, let πj:Ω→[0,1] _j: →[0,1] be the target pmf and Kj:Ω×Ω→[0,1]K_j: × →[0,1] a Markov kernel satisfying invariance ∑ω~∈Ωπj(ω~)Kj(ω~,ω)=πj(ω) _ ω∈ _j( ω)\,K_j( ω,ω)= _j(ω) for all ω∈Ωω∈ . All distributions and kernels are over a discrete state space, so we work with probability mass functions and sums. Let Π(ω1,…,ωN)=∏j=1Nπj(ωj) ( _1,…, _N)\;=\; _j=1^N _j( _j) and K(ω~,ω)=∏j=1NKj(ω~j,ωj).K( ω,ω)\;=\; _j=1^NK_j( ω_j, _j). (1) Product updates. For any ω=(ω1,…,ωN)ω=( _1,…, _N), (ΠK)(ω) ( K)(ω) =∑ω~∈ΩNΠ(ω~)K(ω~,ω) = _ ω∈ ^N ( ω)\,K( ω,ω) =∑ω~1,…,ω~N∏j=1Nπj(ω~j)Kj(ω~j,ωj) = _ ω_1,…, ω_N\; _j=1^N _j( ω_j)\,K_j( ω_j, _j) =∏j=1N(∑ω~j∈Ωπj(ω~j)Kj(ω~j,ωj)) = _j=1^N ( _ ω_j∈ _j( ω_j)\,K_j( ω_j, _j) ) =∏j=1Nπj(ωj)=Π(ω), = _j=1^N _j( _j)= (ω), where we used the assumed invariance of each coordinate kernel, ∑ω~jπj(ω~j)Kj(ω~j,ωj)=πj(ωj) _ ω_j _j( ω_j)K_j( ω_j, _j)= _j( _j), and applied repeated distributivity, which allowed us to change the order of summation and multiplication. Thus K leaves Π invariant. (2) Swap updates. Fix i and define the swap map T on ΩN ^N by exchanging coordinates i and i+1i+1: T(ω1,…,ωi,ωi+1,…,ωN)=(ω1,…,ωi+1,ωi,…,ωN)T( _1,…, _i, _i+1,…, _N)=( _1,…, _i+1, _i,…, _N). For any state x∈ΩNx∈ ^N, Π(T(x))Π(x)=πi(xi+1)πi+1(xi)πi(xi)πi+1(xi+1), (T(x)) (x)= _i(x_i+1)\, _i+1(x_i) _i(x_i)\, _i+1(x_i+1), so the Metropolis acceptance probability α(x)=min1,Π(T(x))/Π(x)α(x)= \1, (T(x))/ (x)\ coincides with (6). The swap proposal is deterministic and involutive (T(T(x))=xT(T(x))=x), hence symmetric, and we have the pointwise identity Π(x)α(x)=minΠ(x),Π(T(x))=Π(T(x))α(T(x)). (x)\,α(x)= \ (x), (T(x))\= (T(x))\,α(T(x)). Since the swap kernel Si,i+1S_i,i+1 only transitions from x to T(x)T(x) (with probability α(x)α(x)) or stays at x, this equality is exactly the detailed balance condition for the pair (x,T(x))(x,T(x)), implying ΠSi,i+1=Π S_i,i+1= . (3) Compositions. If ΠP=Π P= and ΠQ=Π Q= , then Π(PQ)=(ΠP)Q=ΠQ=Π (PQ)=( P)Q= Q= . Therefore any finite composition of product updates K and swap updates Si,i+1S_i,i+1 leaves Π invariant, and so does the overall transition kernel used in Algorithm 1. (4) Marginal correctness of the lowest-noise chain. By construction, the ωN _N-marginal of Π is πN _N. Hence, when the joint chain is at stationarity, ωN∼πN _N _N. ∎ Remark on Proposition 2. This expected-acceptance/overlap identity is standard in replica exchange; see Kofke [2002] and its erratum [Kofke, 2004], which give an exact expression for the mean exchange acceptance probability in terms of the overlap of adjacent energy distributions (Eq. (7) of the original article). Proof for Corollary 1. Proof. Write δ=πi−πi+1δ= _i- _i+1. Then μ−μswap μ-μ^swap =πi⊗πi+1−πi+1⊗πi = _i _i+1- _i+1 _i =δ⊗πi+1−πi+1⊗δ. =δ _i+1- _i+1 δ. Using ‖A‖TV=12‖A‖1\|A\|_TV= 12\|A\|_1, the triangle inequality, and ‖δ⊗π‖1=‖δ‖1‖π‖1=‖δ‖1\|δ π\|_1=\|δ\|_1\|π\|_1=\|δ\|_1, ‖μ−μswap‖TV \|μ-μ^swap\|_TV ≤12(‖δ‖1+‖δ‖1) ≤ 12(\|δ\|_1+\|δ\|_1) =‖δ‖1 =\|δ\|_1 =2‖πi−πi+1‖TV. =2\| _i- _i+1\|_TV. Apply Proposition 2. ∎ Proof for Proposition 3. Proof. Construct a coupling (Y,Y′)(Y,Y ) with Y∼pνY p_ν and Y′∼pν′Y p_ν as follows: sample X∼pX p, sample i.i.d. U1,…,Ud∼Unif(0,1)U_1,…,U_d (0,1), set Bj=Uj<νB_j=1\U_j<ν\ and Bj′=Uj<ν′B _j=1\U_j<ν \, and define Yj=Xj⊕BjY_j=X_j B_j and Yj′=Xj⊕Bj′Y _j=X_j B _j. Then ℙ(Y≠Y′)≤∑j=1dℙ(Bj≠Bj′)=d|ν−ν′|P(Y≠ Y )≤ _j=1^dP(B_j≠ B _j)=d|ν-ν |. By the coupling characterization of total variation, ‖pν−pν′‖TV≤ℙ(Y≠Y′)≤d|ν−ν′|\|p_ν-p_ν \|_TV (Y≠ Y )≤ d|ν-ν |. ∎ Proof for Corollary 2. Proof. Combine Corollary 1 with Proposition 3. ∎ Proof for Proposition 4. Proof. In Δi,i+1 _i,i+1, each ρi _i appears once with a plus and once with a minus sign, so additive constants cancel. ∎