Asymptotically exact variational flows via involutive MCMC kernels
Authors: Zuheng Xu
Trevor Campbell
Venue: NeurIPS 2025 poster ()
Keywords: Bayesian
variational inference
MCMC
convergence
Primary area: probabilistic_methods
TL;DR: A general recipe for constructing tuning-free, asymptotically exact variational flows from general involutive MCMC kernels.

Abstract: Most expressive variational families---such as normalizing flows---lack practical convergence guarantees, 
as their theoretical assurances typically hold only at the intractable global optimum.
In this work, we present a general recipe for constructing tuning-free, asymptotically exact variational ﬂows on *arbitrary* state
spaces from involutive MCMC kernels.
The core methodological component is a novel representation of general involutive MCMC kernels as invertible, measure-
preserving iterated random function systems, which act as the ﬂow maps of our variational ﬂows. This leads to three new variational families with provable total variation convergence. Our framework resolves key practical limitations of existing variational families with similar guarantees (e.g., MixFlows), while requiring substantially weaker theoretical assumptions. Finally, we demonstrate the competitive performance of our ﬂows across tasks including posterior approximation, Monte Carlo estimates, and normalization constant estimation, outperforming or matching No-U-Turn sampler (NUTS) and black-box normalizing ﬂows.

Asymptotically exact variational flows via involutive MCMC kernels
Zuheng Xu
Trevor Campbell
Department of Statistics
University of British Columbia trevor] @stat .bc.
[zuheng xu
Abstract
Most expressive variational families-suue as normalizing flows Jack-lack practical convergence guarantees, as their theoretical assurances typically hold only at the intractable global optimum. In this work, we present general recipe for con- structing tuning-free, asymptotically exact variational flows on arbitrary state spaces from involutive MCMC kerness. The core methodological component is novel representation of general involutive MCMC kernels as invertible, measure- preserving iterated random function systems, which act as the flow maps of our variational flows. This leads to three new variational families with provable total variation convergence. Our framework resolves key practical limitations of exist- ing variational families with similar guarantees (e.g.. MixFlows), while requiring substantially weaker theoretical assumptions. Finally, demonstrate the com- petitive performance of our flows across tasks including posterior approximation, Monte Carlo estimates, and normalization constant estimation, outperforming matching No-U-Turn sampler (NUTS) and black-box normalizing flows.
Introduction
Variational inference (VI) [1-3] is a general methodology for approximate probabilistic inference, where the goal is approximate target distribution (e.g., Bayesian posterior) within specified variational family. This variational family is typically chosen to be a parametric family that enables tractable inference-alloving for i.i.d. sampling and density evaluation [2-7]. This tractability offers key benefits: it enables the evaluation and optimization of approximation quality via unbiased estimates the evidence lower bound (ELBO) [3], which corresponds the Kullback-Leible (KL) divergence [8] to the target distribution up to constant. Moreover, it facilitates downstream tasks such importance sampling [9, 10] and normalization constant estimation.
The quality of a variational approximation is fundamentally determined by the expressiveness of its variational family. Significant progress has been made in constructing flexible families, includ- ing boosted mixtures [11-16] and normalizing flows [4, 6, 17-20J. These families often exhibit universal approximation guarantees [16. 21, 22]: the number of mixture components flow layers grows, the family can approximate any distribution arbitrarily well under mild assumptions. However, major limitation cemains--thorrtiie approximation, which is rarely obtained in practice due to non-convex optimization. In contrast, Markov chain Monte Carlo (MCMC) [23, 24; 25, Ch. 11, 12] is asymptotically exact, meaning is guaranteed to produce arbitrarily accurate results given sufficient computation for any valid choice of tuning parameters (though some values may yield higher efficiency than others).
guarantees pertain only the optimal variational
This has motivated extensive work aimed at bridging VI and MCMC. In particular, many recent approaches fall under the categories of differentiable annealed importance sampling (DAIS) or dif- ferentiable sequential Monte Carlo (DSMC) [26-36]. These methods can be interpreted as gradient- based tuning of AIS/SMC exploration backward kernels improve approximation quality. How-
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

ever, their theoretical guarantees in the limit of flow length [29, 31, 37] often rely on idealized assumptions suchh as perfect transitions or diminishing stepsizes- that rarely hold in practice or only apply under optimal tuning. Moreover, AIS/SMC methods are known to be highly sensitive to tuning, and DAIS/DSMC methods inherit substantial tuning cost to ensure robust performance [38].
Xu et al. [39] introduced MixFlow, an asymptotically exact variational family that does not require optimal tuning. A MixFlow is constructed by averaging pushforwards of a reference distribution under repeated application of an invertible map. When this map is both ergodic and measure- preserving (e.m.p) with respect to the target distribution म, MixFlows converge to in total vari- ation as the number of steps increases, while retaining the tractability of standard variational in- ference. However, its practical applicability is limited by the challenge of designing an invertible T-e.m.p. map for general continuous targets (several solutions exist for discrete spaces [40, 41]). The main obstacles are: (1) continuous e.m.p. maps often involve simulation of ODEs, which re- quires discretized numerical methods that destroy the e.m.p. property: (2) exactness often requires discrete Serropolis-Hastinn (MH) corrections that are not invertible; and (3) proving ergodicity of such maps is very challenging. For example, Xu et al. [39] proposed a map based on the uncor- rected Hamiltonian Monte Carlo (HMC), which is neither exactly measure- praeerrepess nor provably ergodic. Other existing Hamiltonian- baaed methods [42] also suffer from discretization error and are non-ergodic [43]. Attempts via deterministic Gibbs samplers based on measure- pressrring ODEs [44] or CDF/inverse-CDD the exact transformations. MH corrections used to restore exactness [40, 45] result in non-invertible transformations due to the accept-reject mechanism; recall that invertibility is required by varia- tional flows to enable tractable density evaluation. To date, there is no framework for constructing variational families whose practical implementation achieves MCMC-like asymptotic exactness.
transformations [40] are also limited by the intractability of computing
In this work, we address the challenges mentioned above and propose a new framework for devel- oping practical, asymptotically exact variational flows. Rather than relying on e.m.p dynamics as in MixFlow [39], our framework leverages iterated random functions (IRF)* [47]-a a type of random dynamical system. The main contributions of this work are as follows:
1. We develop a method for deriving exact measure-preserving transformations from general involutive MCMC kernels [48, 49], while preserving invertibility of the transformation.
2. We introduce a more general framework for constructing asymptotically exact flows, lead- ing to three novel variational families beyond the original MixFlow for general state spaces.
3. We establish total variation convergence guarantees for these new families under signifi- cantly weaker assumptions than those required in MixFlow theory [39], notably relaxing the ergodicity conditions of the flow maps.
2 Background
Throughout, let be target distribution a measurable space (X, equipped with a o-finite base measure m. All distributions are assumed to have densities with respect to the base measure on their corresponding spaces, and use the same symbol denote both a distribution and its density. Given a transformation f and a distribution p. we write fp((z) )) for the function f evaluated at p(z), and fp(x) for the density of the pushforward distribution fp evaluated at ...
2.1 Homogeneous MixFlows
A mixed variational flow (MixFlow) [39, 41] is built from a deterministic, -rerodoic (Definition D.1) and measure- measure-preserviee preserving (e.m.p.) diffeomorphism f.. Given such a map and reference distribu- tion go on X' that enables i.i.d. sampling and density evaluation, the MixFlow density is given by
1 Va E X, -r(x)= =
Σ (m((x) -
= 7 Σ go (5'xx)
J(x) |det Vs(x)), (1)
J
where f and f'ro denote mapping x or pushing go through t (t > 0) iterations of f. We use the convention that qo = go (MixFlow of length 0 is just the reference distribution qo). Eq. (1) is tractable if f-1 and the Jacobian J can be evaluated. To generate X ~ TT, we first draw Xo ~
I IRFs are also referred to as iterated function systems (IFS) in some literature, e.g. [46]. 2 is diffeomorphism if is continuously differentiable and has continuously differentiable inverse.

---

Qu and a flow length K ~ 2 and then map Xo through K iterations of f. i... X =f Since qr built from a homooenoous e.m.p dynamical system, we label it a homogeneous MixFlow, to distinguish it from our proposed random dynamical system flows (see Section 4). The asymptotic exactness of homogeneous MixFlows comes from the fact that lim TV T) = O regardless of the tuning of the flow map f [39, Theorem 4.2].
In practice, the map/ is typically designed mimic familiar MCMC kernels [39, 41], so that its tra- jectories have similar statistical behavior to the corresponding Markov chain. Despite this, general constructions of exact e.m.p. MixFlow maps for continuous target distributions remain unavailable. As discussed in the introduction, achieving both exact measure preservation and ergodicity is highly non-trivial in practice. Consequently, practitioners often rely on approximate maps, leading gup between theoretical guarantees and practical implementations. These approximations can introduce numerical instability and degrade performance as increases [39, 50). In Section 4.1, we show how to design homogeneous MixFlows that are exact in practice. Additionally, we present refined char- acterization of the density Tr by leveraging the measure- prrssevvng property of f, which simplifies implementation, improves robustness, and provides a more intuitive convergence analysis.
2.2 Involutive MCMC
An involutive MCMC kernel [48, 49, 51] is a Metroplis-type Markov kernel with a deterministic proposal defined by an involution g. i.e., self-inverse function satisfying g = و g- 'his framework encompasses a broad class of MCMC algorithms, with many popular algorithms appearing as special cases [49, 52-55] (see Appendix A.1). The detailed transition procedure of involutive MCMC is described in Algorithm 1 of Appendix A.2. Consider an auxiliary variable e defined on a space V. with conditional density p(u x) given æ X with respect to base measure my on V, and the augmented target density (2x, v) == π p(r))))) Let m == m x my be the joint base measure on X' V. each transition from state æ proceeds in three steps:
X' V. For an involution g:X
1. Sample an auxiliary variable ~ p(dv 1x);
2. Propose a new state (z', = g(x, v);
3. Accept x' with probability min (1, ale''dd J (x) where J (z, v) = dom
dm (z, u)?
An involutive Markov kernel K defined this way is reversible with respect both the augmented target (II, v) and its marginal (2x) [51, Theorem 2]. Proposition 2.1. The involutive MCMC kernel K(x', o/1 v) (defined in Algorithm 1) satisfies that K('', v/|o,, ())( v) = K(x, ula', , (') ,'), ((a(())((x) = (xx)=')<xx')
where K is the marginalized kernel defined as: [x][z) K(x', of (z)))(()))))de)
2.3 Iterated random functions
An iterated random function (IRF) system [47] on A' consists of sequence of random maps: iid
Vt E N, X, (X..)(.)), E (0.)teN p4,
where (ʃo : X + X E 0} is a set of parametrized functions, with each drawn randomly from a distribution y the parameter space O The above IRF induces a Markov kernel given by:
Vx E X, BB P(x, = / 1 [=(fe(=)))((d).
(3)
This yields simple characterization of the action of the Markov process P on distribution q: Pa(y) == fx P(z. v)g(dx) = =E [Jo(v))), 0 ~ n4 foq: pushforward of q under fo-
(4)
Throughout this work, we focus on IRFs where the family (so : 0 E Q} satisfies Assumption 2.2. Assumption 2.2. For p-a.s. all 0 O. fo(.) is bijective and r-measure-prest preessrr-preserving (T-m.p.). Fur- thermore, is the unique invariant distribution of the Markov kernel P induced by the IRF.
Forr differentiable continuous state spaces (e.g.. RR), Jy (x, v) = /det Vg(z, v) is its Jacobian deter- minant. We adopt the theserc- formulation of Tierney [51] handle arbitrary state spaces.
3

---

Assumption 2.2 implies that the sequence of iterates X, produced by the IRF behave like a m invariant, irreducible Markov chain. Therefore, long-run averages of IRF iterates converge to expec- tations under и, following the standard law of large numbers (LLN) for MCMC [56], also known as the random Birkhoff ergodic theorem the IRF literature [57; 58, Cor. 2.2.]. Theorem 2.3 syn- thesizes these results under Assumption 2.2, providing a unified statement for convenient use in our framework; proof can be found in Appendix D.I.
Theorem 2.3. Suppose that IRF So satisfies Assumption 2.2. Then, given @ [T(T),, we have that
1. for u-..e. x E X and p-almost all (01)LEN, as T + 00:
1
1
Σ
T-1
Φ (fe. o...o fer (x)) [[(()), X X~T;
(5)
T
2. for p-almost all (fe) as 00:
Σ à (So, o... la (x)) ( B[o(X), X~T.
(6)
T c
Moreover, the same results hold for the inverse IRF (5o e)
3 Invertible measung- IRF from involutive MCMC
In this section, we provide concrete, general construction of invertible and exactly measure- preserving IRFs based on involutive MCMC kernels. The key idea, originally developed for the MH sampler [40, 45], is to further augment the space with two additional variables u, E [0, 1J uha [0, 1]. The variable uh pairs with the auxiliary variable v of dimension d, and ua encodes the ran- domness in the accept/reject decision. Let the augmented target मे and space S be defined as:
* (s) = (<))((v) ax) 10.1j (w,) [0.1 (ua).
(2r, U, the Un, Ua) E X x Vx [0, [0.
The two uniform auxiliary variables u, and ua will be refreshed with two random parameters (O., Pa) ~ / Unif(O, 1jd x 1n. Without loss of generality, we describe the IRF con- struction assuming dimensinnsi target T(x). The IRF fo(s) = fo(a,, v, un ua) is defined by the following steps (Algorithm 2):
1. Uniform auxiliary refreshment: the + (up + 0p) mod 1, ua t (ua + +aa) mod 1
2. Update (v, un) pair via CDF/inverse-CDF of p(- 4: u/ +- (p(l))()), T Fotl (uo)
F(r',0') [2( 7)) 3. Propose and compute the MH-ratio: (x', u(,, rt न(.,8)
4. Accept or reject: If uu > reject and stay at the pre- involution state s' (x. U. uho, ua). Otherwise, set u/, < Su and accept the post-involution state s' (a', v' u/ v'a)-
The correspondence with involutive MCMC (Algorithm 1) is: Step 2 simulates v ~ p(dvlz) via in- verse CDF sampling, Step 3 mirrors the involution and MH ratio computation, and Step 4 performs the accept/reject step while explicitly tracking the randomness th involved in the decision. Further- more, as mentioned in Section 4.1, one can use this map fixing f., Og to some pre-specified constant values (rather than sampling from μı). And finally, using the same Jacobian computation as in [45, Eq. (25)), one can show that 40 (0z, 0a) O, the IRF fo (Algorithm 2 of Appendix B) is π --measure- pressevvn......ing...
homogeneous MixFlow by simply
Proposition 3.1. The map given by Algorithm 2 satisfies Assumption 2.2 for т. its induced Markov kernel P is irreducible.
One must be able to compute f (-) if fo is to be used as flow layer in a MixFlow. Steps 1- 3 are straightforward to invert. The main challenge lies in inverting the accept/reject Step we need to recover the accept/reject decision based solely on the output state s'. Depending differ- ent decisions, s' could either be the pre-involution state (a, U, uhe, ua) or the post-involution state (a', v', u/g, v')). Since the transformation v/a ua/r only present in the acceptance branch, infer- ring the branch incorrectly would lead to the failure of recovering Ua (hence the entire state s).
'Typically, p(v lies in simple family; for instance, in HMC with diagonal mass matrix, is a diagonal Gaussian, whose CDF and inverse-CDD CDF can be computed stably. For multidimensional v, a Gibbs-style update on the conditionals of p(-(2) can be used.

---

meaiibitty
invertiiity
Scsett
Figure 1: Inversion error of IRFs (based on HMC, uncorrected HMC, MALA, and RWMH) over increasing flow length T. Verticle axis shows the 2-norm error of reconstructing
(x, up Ut)
(s = (ar v) for the uncorrected HMC IRF) sampled from go by the composing the forward simulation for ....... fo, (s) and its inverse. The lines indicate the mean, and error regions indicate the standard deviation over 32 independent initializations from go.
We address this challenge (pseudocode in Algorithm 3 of Appendix B) by exploiting the self-inverse
property of the involution 9. First note that g(r, V) = (x',0') and g(æ', v²) (ar, 7). Suppose that s' = (I# w#, u ux {g(a* 1 ), (I#H, wt }} is exactly the unordered pair {(2, ū), (a', v')). Then from the property of the Jacobian of g (i.e., (2 () = Jy (x', v/) and vice versa), we observe
T(X' J (2, (ス,, 8)
[") Jo ( ,2,)')
(1x 6)
Hence, recomputing the MH-ratio as in Step yields
(x* wtt) I (g(x# x#))
Jn (g(** w+)) E (r,r
where r corresponds to the true MH-ratio as computed in the forward pass. The key observation to infer the accept/reject decision then follows: If u 7 < 1, then the acceptance branch was taken, so Ua = uu F;; otherwise the move was rejected as Ua cannot be larger than 1.
Fig. empirically verifies that one can successfully invert the proposed IRF map for four MCMC- based IRFs HMC[5, 60], uncorrected HMC [61]. MALA[62), and RWMH [63] on four syn- thetic targets defined in Appendix E. The same hyperparameters are used for every example: each (uncorrected) HMC transition consists of 50 leapfrog steps with step size 0.02; MALA uses step size 0.25; RWMH uses step size 0.3. We evaluate the 2-norm error of reconstructing sampled from mean-field Gaussian variational approximation go by the composing the forward simulation far ....o ible up to T R 200 iterations, while the RWMH IRF remain invertible up to T x 1000 iterations. Notably, the corrected HMC IRF is consistently more stable than its uncorrected counterpart used in past MixFlows work; the additional MH step discards trajectories with large numerical error that would otherwise cause the dynamics to diverge. Although Algorithm and Algorithm 3 are exact inverses in theory, floating-poit point round-off accumulates with T and exact reconstruction can fail [39, 50). In practice, however, the resulting statistical error in downstream variational inference is often negligible, thanks to the shadowing property of chaotic dynamical systems [50].
(x, v, uy Ua)
o fo, (s) and its inverse. Both HMC variants and MALA remain reliably invert-
Variational flows based on IRFs
In this section, we present a methodology that transforms any IRF system satisfying Assumption 2.2 into asymptotically exact variational family. Alongside the exact homogeneous MixFlows derived from IRFs and their refined analysis, we introduce three additional families constructed from the same IRF but combined differently and show that all converge to the target in total variation. Proofs are deferred to Appendix D. For simplicity, we present the methodology and theory using IRFs defined directly the original space X' rather than the augmented space S. We also assume access to a reference distribution go supporting i.i.d. sampling and tractable density evaluation.
4.1 Improved homogeneous MixFlows
As reviewed in Section 2.1, the homogeneous MixFlow @r is defined as 'T = + Sesn-''a f'go with the convention go = -o- Given an IRF sa satisfying Assumption 2.2, one can construct a homogeneous flow map f by fixing the parameter 0 to a constant value 0* (e.g.. r/16), rather than sampling
5

---

from the distribution μ. This provides generic way of building exact n-m.p. flow maps. A key property not noted in prior MixFlow work [39, 41] simplified expression for the density of gr on arbitrary state spaces, enabled by measure-theoren a measure-pees paeeurrep map. Specifically, for any x-m.p. bijection f. foo(a (a) = (a) 20 ( (f- as introduced in Appendix C. This yields a simplified form for the density of @r (in contrast to Eq. (1)):
formulation of the pushforward density under
1 T 7 E (ox)) = =xxx== 7 Σ (F(()), E
T
Cr (r(x)
In practice, this expression is particulary useful for evaluating the flow density; practitioners can evaluate (r))) without tracking the Jacobians of f explicitly, which simplifies implementation and avoids numerical instability from accumulating Jacobians over long trajectories.
Moreover, the explicit expression Eq. (7) offers an intuitive understanding of why qr converges. While the original convergence result in [39, Theorem 4.2] relied on general operator theory for e.m.p. systems [64], the density-based perspective is more transparent. ʃ is m-..m.pp the Birkhoff ergodic theorem [65; 64, p. 212] implies that +. Σi=1 @
& (x)) - 1. Consequently. for т-aee X. (x) - (x) as T + This enables substantially simplified proof of the convergence of homogeneous Mixflow. The proof of Theorem 4.1 can be found in Appendix D.2.
Theorem 4.1. Suppose that is a TT..mm.p diffeomorphism, and qo < π. Then, as T 00,
qr(x) + [(x), T--.e. X, and (V πT 0.
is worth noting that Assumption 2.2 does not guarantee the ergodicity of specific for leaving a gap between theory and the practical implementation of homogeneous MixFlows. In the remainder of this section, we introduce three new MixFlow families designed address this limitation.
4.2
IRF MixFlows
An IRF MixFlow is mixture of pushforwards of reference go through an IRF sequence:
+
Σ Jor o...o fo Yo. with the convention that go -0+
where .I..... Or is a cached i.i.d. sequence drawn from μ.. When constructing the flow, we first sample and freeze the random stream fl...., fr. yielding an inhomogeneous sequence of param- eterized bijections. Then to draw X ~ qī, we treat qi as mixture of T distributions:
K (1,
Xo ~000 40
X fex o...o fo, (Xo)
Note crucially that each sample X is generated using the same frozen sequence 01,..., Gr. For density evaluation, we compute the inverse IRF før.... f--1 Because each fo is x-m.p., by Proposition C.3, the density takes similar form as in a homogeneous MixFlow (Eq. (7)):
T
(((z) = (z)
go (fo,
o... oo.
VI
t=1
However, note that this density requires simulating the backward process of the inverse IRF ([47)) (x(x)= : So,
So (x) for & [T],
which cannot be computed sequentially. As a result, IRF MixFlows incur a quadratic density evalu- ation cost O(T?). Fortunately, this backward process can be computed in parallel fashion, as the
computation of each Xx( [T'] is independent. We recommend deploying IRF MixFlows on modern parallel hardware (e.g., GPUs) for efficient density evaluation.
IRF MixFlows share the total variation convergence guarantee (Theorem 4.2) of homogeneous MixFlows. The proof (Appendix D.3.1) is similar to the original MixFlow argument [39, Theo- rem 4.2], interpreting the IRF (Eq, (2)) as a time-homogeneous, e.m.p. dynamical system over the joint space ON x X. However, we emphasize that Assumption 2.2 is significantly weaker than the ergodicity assumption of Theorem 4.1. See Section 4.5 for detailed discussion.
Theorem 4.2. Let P denote the joint distribution over the i.i.d. sequence (0.)teN μ If Assump- tion 2.2 holds and do < n, then
TV (qi,
T) =,0 as - oc.
(8)
An implication of this result in continuous state space is that for any x-m x-m.p. diffcomorphism f. the Jacobian determinant must satisfy /det V (x))
z((-*=) established in Proposition C.3.
.())

---

<(M) 4.3 Backward IRF MixFlows
To address the ( density cost of IRF MixFlows, we propose a simple modification: constructing the flow from the backward process. Specifically, we define the backward IRF MixFlow as:
I
E fo ...... fa,90, with the same convention that do = qu- 1=1
T
This construction retains (TT) complexity of sampling X 4r via: K ~ Unif 2, 21. ...,,
40 while reducing the density computation cost to O(T). The density of रr is given by:
Xo
x= So, ..... Sox (Xo),
T
-
Σ E=1
go
(Jo.
.....
E
T
This mirrors the density formula of homogeneous MixFlows (Eq. (7)), enabling the use of the ran- dom ergodic theorem (Theorem 2.3) to establish the same pointwise and total variation convergence. Theorem 4.3. If Assumption 2.2 holds and qu & then π-..e. E X' and p-almost all (6t )eee 0y (z) (x) and TV a) + T 20.
4.4 Ensemble IRF MixFlows
All MixFlow variants discussed above-including homogeneous MixFlows-and are based on ergodic averaging along the flow. This inherently limits their convergence rate to O1/TT), the first com- ponent always retains 1/T mixing weight. In contrast, MCMC methods often exhibit geometric convergence in their marginal distributions under suitable conditions [56; 66, Ch. 15]. Motivated by this, we propose the ensemble IRF MixFlows, which instead uses an ensemble average of the endpoint of multiple IRF trajectories in attempt to match T T'sstep MCMC marginal distribution:
M
M
(M) 9T
Σ m=1
fomm o ..... form9u
M
M
m=1
where each 0f_h)_.
corresponds to an independent IRF realization. As the case of the previous MixFlows, the M streams of randomness p\m) cached frozen) when sampling and computing densities. The resulting density of the ensemble IRF MixFlow is given by:
HM (x) = (<)))
Σ 90
(z)
M
pfm))
whose computation costs (T (or ((T) when parallelized across the M streams). Drawing
X fr (M) takes OCT + M) operations:
K 2n,
X qu
X f o.... fop390
Intuitively, the flow length controls the bias the IRF system, while the ensemble size M controls the variance of the Monte Carlo average. This tradeoff is formalized in the following result. Theorem 4.4. Suppose that Assumption 2.2 holds, and that Va
X, 2o (x) I B co. Then,
Ea
TV
GT
1
E
Varour
Goo (J o...o So (x)) x
VM
a
+ B. E [TV(R?
( 3x, X~T.
where R is the Markov kernel induced by the inverse IRF fo and 8x is the Dirac measure at X.
In the setting where TV 8x π) -Opp") for some p 1), the convergence rate of 4r can be heuristically characterized TV (-(M)
O 0 VM capturing the tradeoff between the bias (via 7) and variance (via M). Given a fixed computational budget, choosing the balance between flow length T' and ensemble size M is critical. In the extreme case of M = convergence will fail entirely- any r-mmesurr-preeeervnn preserving f satisfies TV(fq, т) TV x) [67, Theorem 1]. On the other hand, small leads to high bias due to insufficient mixing. This tradeoff closely relates to recent studies on parallel MCMC algorithms [68, 69].
07 (z(M),=)
7

---

4.5 Discussion
Relaxing ergodicity. A major advantage of IRF-d MixFlows over homogeneous MixFlows is that IRF-based MixFlows require only that the kernel P admits unique invariant distribution (Assumption 2.2), a significantly weaker condition than the ergodicity assumed by homogeneous MixFlows. In fact, whenever the set O* := {0 So rr-eroddic) has positive jı-measure, Assump- tion 2.2 automatically holds [46, Corollary 3.3]. Uniqueness of the invariant distribution is also easily verified by checking that P is irreducible [56, 66). The IRFs we construct in Section 3 corre- spond to involutive MCMC kernels that are known to be irreducible, whereas establishing ergodicity in MixFlows is typically so difficult that it is assumed without proof [39, 42, 70].
Which flow to choose? All four flows are asymptotically exact, yet their density formulae reveal different bias-variance and cost-aceuracy trade-offs. In every case the density ratio takes the form density (x) = N Σ 40 (I)(,)), where Tn is a composition of inverse IRF/ergodic maps, and N can be the flow length or ensemble size. Hence practical convergence of each flow is dictated by how quickly N >m=1 da (x)) converges constant. Empirically (see Appendix E.1.I) we find that IRF MixFlows often reach given accuracy at shorter flow lengths than homogencous or backward IRF MixFlows, but a full theoretical comparison study is deferred future work.
5 Experiments
This section presents an empirical evaluation of the four proposed flows- three IRF variants and homogeneous MixFlows (collectively referred to as "IRF flows" since homogeneous MixFlows can be viewed as a special case). We compare them against two normalizing flows, RealNVP [19]
and Neural Spline Flow (NSF) [71], and against the No-U-Turn Sampler (NUTS) [72]. Variational
methods are assessed by their (i) ELBO and (ii) accuracy of the importance sampling estimate of the normalization constant log for the unnormalized density 7:
(Xn), 4T
(X.)
idd ~ QT:
where or (T,, Tr, 4 Ti T (7
2 j: 7
n=1
and (ii) importance sampling effective sample size (ESS) [73-75]. Sampling methods are evalu- ated via their Monte Carlo estimation error. In all cases, all flows start from the same reference Gaussian trained for 10K Adam steps with batch size 10 and learning
distribution qn: a mean- rate 10-3, All IRF flows are evaluated with 64 i.i.d. draws, while normalizing flows use 1024. Full experimental details appear in Appendix E.
5.1 Synthetic examples
Our synthetic experiments consist of four 2-dimensional targets used by Xu et al. [39]: the Ba- nana [76]. Neal's funnel [77], a cross-shaped Gaussian mixture, and warped Gaussian distribution. Fig. shows a comparison of the original Hamiltonian-MixFUn Mwx--Mi---iiiiiii on an uncorrected HMC kernel-with our corrected version including the MH step. For each target run both flows with identical hyper-paramede parameters (50 leapfrog steps per transition, several step-sizes) and estimate the total- variation (TV) distance to the ground truth using 512 i.i.d. samples. Across all targets and step-sizes, the corrected HMC-based MixFlow consistently achieves lower TV error and remains robust as the step-size grows. In contrast, the uncorrected variant often deteriorates with longer flows because the inexact map error accumulates (e.g.. the green dashed curve in the third panel). At larger step sizes the uncorrected flow frequently diverges, producing NaNs (marked by crosses), whereas the corrected flow remains stable-echoing the inversion stability results in Fig.
We next compare the four IRF flows with RealNVP and NSF. Two IRF variants examined: HMC- based (50 leapfrog steps per transition; T = 200) and RWMH-based (T' = 4000). Each normalizing flow consists 6 flow layers, and is trained via 50,000 Adam steps with batch size 32; we tune the learning rates in the grid {10-4, 10-3 10 10-2). and report the results of the setting with smallest median TV distance over runs. Additional implementation details can be found in Appendix E.
Figs. 3a and 3b display the ELBO and log Z estimates (via importance sampling) for the Banana target; the remaining synthetic cases show the same pattern (Fig. 7 in Appendix 1.1.3). As synthetic targets are normalized, a perfect variational approximation has both metrics near 0. The IRF flows meet this mark consistently across runs, whereas RealNVP and NSF exhibit high variability and often
8

---

Tow Deraetoo
TO Norteatoo
langth
Fax
Fox couleegth
Fox Tocilengh
(a) Banana
(b) Funnel
(c) Cross
(d) Warped Gaussian
Figure 2: Votal-variation error for homogeneous MixFlow built on corrected (solid) versus uncor- rected (dashed) HMC kernels, plotted against flow length T for several step sizes. Each curve the mean over 32 independent runs; shaded bands (+11 SD) show run-to-run variability. A cross marks any setting where at least one run returned a NaN (instability), at which point the trace is terminated.
20
normalization consTait
10
1000
5S5 po0 gampl
800
15
FINSH 057
HINHH MEL
-30
200
perfict OUTple
Racknudd KINoTCA How
Servtion SEAIAAA
(a) ELBO
(b) log estimates.
(c) Per-sample ESS.
Figure 3: Variational approximation quality of IRF Flows versus ReaINVP and NSF. Box plots for IRF flows are based on 32 independent runs, and 10 runs for the normalizing flows. The black dashed line in (c) indicates the optimal ESS of perfect i.d. samples.
produce extreme ELBO log Z values. We restrict the vertical range of the ELBO plot for better visualization; full-range plots are in Fig. 7b. We also note that training instability is common for the normalizing flows: on the Funnel example, 10 of 15 RealNVP runs and all NSF runs diverged.
Fig. 3c further examine the par-sample importance sampling ESS (see Fig. 7d on similar results for other examples), which reflects the x² divergence from the variational distribution to the target [78]. The ESS is orders of magnitude higher for IRF flows than for the normalizing flows. Additionally, provide comparisons among the three ergodic averaging MixFlow variants in Appendix E.1.1, and ensemble-ssizellenghh trade-offs for ensemble IRF MixFlows are explored in Appendix E.
5.2 Real-data experiments
The real-data experiments include the Sregesiot- (TReg; 4-dimensional), and the Sparse linear regression (SparseReg: 83-dimensional) from [39], and latent Brownian motion model (Brownian, :n;; 32-dimensional) and the Log-Gaud -au-siann Cox process model (LGCP; 1600-dimensional) from the Inference Gym library [79]. Each normalizing flow is trained via 50,000 Adam steps of batch size 32; we grid-search both the learning rates {10- 10 10-23 and flow layers {6, 10), and report the configuration with the highest median ELBO over
An additional mean-field Gaussian baseline is optimized for the same number of steps and batch size with learning rate 10-3
All IRF variants use RWMH kernel, with the step size tuned to achieve a 0.8 acceptance rate using bisection search between 0.001 and 10. In each search step, we estimate acceptance rate with 5,000 RWMH-IRF iterations. We set T 5000 for the backward IRF and homogeneous MixFlow and en- semble IRF MixFlow, and set T 4000 for the IRF MixFlow. Normalizing flow results are omitted for LGCP, which did not finish training within 48 hours the same computation cluster. Ground truth values are estimated using AIS with a dense temperature grid; see the details in Appendix E.2.
As in the synthetic experiments, our exact flows match or modestly improve upon-tth best-tuned RealNVP and NSF in both ELBO (Fig. 4a) and log accuracy (Fig. 4b), and outperform the mean- field baseline by wide margin. The per-sample importance-sampling ESS shows the same advan- tage (Fig. 8b). Crucially, normalizing flow training orders magnitude more expensive (Fig. 4d), whereas the exact flows achieve comparable accuracy a fraction of the computational cost.

---

1441
-20
-2600
(a) ELBO
2450
3460
Day
(b) log Z estimates
(c) Maximum absolute error of coordinate wise posterior standard standard-deviation deviation estimates relative to NUTS
THe S0cee3
THE SAGREEO
(d) Computation time (in seconds) for each method. To ensure consistent environment, all timing results obtained by rerunning the methods the local machine (hardware details provided in Appendix MFVI, NSF, and RealNVP were each run once, as their execution times are deterministic given the flow architecture and optimization settings. For IRF flows and NUTS, timing statistics are based on 10 independent runs.
Figure 4: Results on real-data benchmarks (columns, from left right); TReg(d 4). Brownian(d = 32), SparseReg (d = 83), and LGCP (d = 1600).
We further compare coordinate-wise posterior mean estimates (Fig. 8c in Appendix E.2) and stan- dard deviation estimates (Fig. 4c) against NUTS, reporting the maximum absolute error across di- mensions relative the estimated ground truth. NUTS initialized with independent draws from qu and run for 10,000 iterations including 5000 warm-up iterations. IRF flows outperform NUTS on two models and are slightly worse on the other two-ye they do so at generally faster compu- tation time (Fig. 4d). Note that the goal of this work is not to outperform MCMC, but rather to construct a variational family that provides asymptotic exactness and similar sampling performance; IRF MixFlows meet this standard.
Conclusion
We introduced general framework for building asymptotically exact variational families from gen- eral involutive MCMC kernels. By constructing invertible, these kernels, we overcome the main practical limitation of MixFlow [39] and enable the construc- tion of a broad class of exact flows. We also provided a streamlined theoretical analysis for flows based on weasure-preserving transformations and demonstrated their empirical advantages in den- sity approximation and importance sampling. promising direction is pair our framework with recent automatic-tuning MCMC [52-55], developing truly tuning-free exact flows in practice.
maps directly from
10

---

Acknowledgments and Disclosure of Funding
The authors sincerely thank Peter Orbanz for pointing us to the IRF literature, which provided the initial inspiration for this work, and Alexandre Bouchard-Côde for suggesting applying our methods to normalizing constant estimation. T. Campbell and Z. Xu acknowledge support from the NSERC Discovery Grant RGPIN-2025-04200 ing platform at the University of British Columbia, and the compute cluster provided by the Digital Research Alliance of Canada.
We are also grateful for access to the ARC Sockeye comput-
References
[1] Michael Jordan, Zoubin Ghahramani, Tommi Jaakkola, and Lawrence Saul. An introduction variational methods for graphical models. Machine Learning, 37:183 233, 1999.
[2] Martin Wainwright and Michael Jordan. Graphical models, exponential families, and variational inference. Foundations and Trends in Machine Learning, 1(1-2):1-305, 2008.
[3] David Blei, Alp Kucukelbir, and Jon McAuliffe. Variational inference: review for statisticians, Journal of the American Statistical Association, 112(518):559 877, 2017.
[4] Danilo Rezende and Shakir Mohamed. Variational inference with normalizing flows. In International Conference on Machine Learning, 2015.
Rajesh Ranganath, Dustin Tran, and David Blei. Hicrarchical variational models. In International Con- ference on Machine Learning, 2016.
[6] George Papamakarios, Eric Nalisnick, Danilo Jimenez Rezende, Shakir Mohamed, and Balaji Lakshmi- narayanan, Normalizing flows for probabilistic modeling and inference. Journal of Machine Learning Research, 22:1-64, 2021.
[7] Charles Margossian and Lawrence Saul. Variational inference in location-scale families: Exact recovery of the mean and correlation matrix. In Advances in Neural Information Processing Systems, 2024.
[8] Solomon Kullback and Richard Leibler. On information and sufficiency. The Annals of Mathematical Statistics, 22(1):79-86 1951.
[9] Herman Kahn and Andy Marshall. Methods reducing sample size Monte Carlo computations. Journal of the Operations Research Society of America, 1(5):263-278, 1953.
[10] Herman Kahn. Use of diffcrent Monte Carlo sampling techniques. Technical report, Rand Corporation, 1955.
[11] Fangjian Guo, Xiangyu Wang, Kai Fan, Tamara Broderick, and David Dunson. Boosting variational inference. In Advances Neural Information Processing Systems, 2016.
[12] Andrew Miller, Nicholas Foti, and Ryan Adams. Variational hoosting: iteratively refining posterior ap- proximations. In International Conference on Machine Learning. 2017.
[13] Xiangyu Wang. Boosting variational inference: theory and examples. Master's thesis, Duke University, 2016.
[14] Francesco Locatello, Rajiv Khanna, Joydeep Ghosh, and Gunnar Rätsch. Boosting variational inference: an optimization perspective. In International Conference on Artificial Intelligence and Statistics, 2018.
[15] Francesco Locatello, Gideon Dresdner, Rajiv Khanna, Isabel Valera, and Gunnar Rätsch. Boosting black box variational inference. In Advances in Neural Information Processing Systems, 2018.
[16] Trevor Campbell and Xinglong Li. Universal boosting variational inference. In Advances in Neural Information Processing Systems. 2019.
[17] Ivan Kobyzev, Simon Prince, and Marcus Brubaker. Normalizing flows: an introduction and review of current methods. IEEE Transactions on Pattern Analysis and Machine Intelligence, 43(11):3964-33 3979, 2021.
[18] Abhinav Agrawal, Daniel Sheldon, and Justin Domke, Advances in black-box VI: Normalizing flows, importance weighting, and optimization. In Advances in Neural Information Processing Systems, 2020.

---

[19] Laurent Dinh, Jaseha Sohl-Dickstein, and Samy Bengio. Density estimation using Real NVP. In Interna- tional Conference on Learning Represeniations, 2017.
[20] Rianne van den Berg, Leonard Hasenclever, Jakub Tomezak, and Max Welling. Sylvester normalizing flows for variational inference. Conference on Uncertainty in Artificial Intelligence, 2018.
[21] Frederic Kochler, Viraj Mchta, and Andrej Risteski. Representational aspects of depth and conditioning in normalizing flows, In International Conference on Machine Learning, 2021.
[22] Zhifeng Kong and Kamalika Chaudhuri. The expressive power of class of normalizing flow models. In International Conference on Artificial Intelligence and Statistics, 2020.
[23] Christian Robert and George Casella. Monte Carlo Statistical Methods. Springer, 2nd edition, 2004.
[24] Christian Robert and George Casella. A short history of Markov chain Monte Carlo: subjective recollec- tions from incomplete data. Statistical Science, 26(1):102-115, 2011.
[25] Andrew Gelman, John Carlin, Hal Stern, David Dunson. Aki Vehtari, and Donald Rubin. Bayesian data analysis. CRC Press, 3rd edition, 2013.
[26] Tim Salimans, Diederik Kingma, and Max Welling, Markov chain Monte Carlo and variational inference: bridging the gap. In International Conference on Machine Learning, 2015.
[27] Christopher Wolf, Maximilian Karl, and Patrick van der Smagt. Variational inference with Hamiltonian Monte Carlo. arXiv: (609.08203, 2016.
[28] Tomas Geffner and Justin Domke. MCMC variational inference via uncorrected Hamiltonian anncaling. In Advances in Neural Information Processing Systems, 2021.
[29] Guodong Zhang, Kyle Hsu, Jianing Li, Chelsea Finn, and Roger Grosse. Differentiable annealed impor- tance sampling and the perils of gradient noise. In Advances in Neurul Information Processing Systems, 2021.
[30] Achille Thin, Nikita Kotelevskii, Alain Durmus, Maxim Panov, Eric Moulines, and Arnaud Doucet. Monte Carlo variational auto nuto-encoders. International Conference on Machine Learning. 2021.
[31] Martin Jankowiak and Du Phan. Surrogate likelihoods for variational annealed importance sampling. 1rXiv:21 12194, 2021.
[32] Tomas Geffner and Justin Domke. Langevin diffusion variational inference. In International Conference on Artificial Intelligence and Statistics, 2023.
[33] Michael Arbel, Alexander Matthews, and Arnaud Doucet. Annealed flow transport Monte Carlo. In International Conference on Machine Learning, 2021.
[34] Alex Matthews, Michael Arbel, Danilo Jimenez Rezende, and Arnaud Doucet. Continual repeated an- nealed flow transport Monte Carlo. In International Conference on Machine Learning, 2022.
[35] Christian Naesseth, Scott Linderman, Rajesh Ranganath, and David Blei. Variational sequential Monte Carlo. In International Conference on Artificial Intelligence and Statistics, 2018.
[36] Chris J Maddison, John Lawson, George Tucker, Nicolas Heess, Mohammad Norouzi, Andriy Mnih, Arnaud Doucet, and Yee Teh. Filtering variational objectives. Advances in Neural Information Processing Systems, 2017.
[37] Johannes Zenn and Robert Bamler. Differentiable annealed importance sampling minimizes the sym- metrized Kullback-Leibler divergence between initial and target distribution. In Forty-first International Conference on Machine Learning (ICML), 2024.
[38] Kyurae Kim, Zuheng Xu, Jacob Gardner, and Trevor Camphell. Tuning sequential Monte Carlo samplers via greedy incremental divergence minimization. In International Conference on Machine Learning, 2025.
[39] Zuheng Xu. Naitong Chen, and Trevor Campbell. MixFlows: principled variational inference via mixed flows. In International Conference on Machine Learning, 2022.
[40] Radford Neal. How to view an MCMC simulation as permutation, with applications to parallel simulation and improved importance sampling. arXiv: 2012.
12

---

[41] Gian Carlo Diluvi, Benjamin Bloem-Reddy and Trevor Campbell. Mixed variational flows for discrete variables. In International Conference Artificial Intelligence and Statistics, 2024.
[42] Greg Steeg and Aram Galstyan. Hamiltonian dynamics with oo-NNewtonian momentum for rapid sampling. Advances Neural Information Processing Systems, 2021.
[43] Jakob Robnik, Bruno De Luca. Eva Silverstein, and Uroš Seljak. Microcanonical Hamiltonian Monte Carlo. Journal of Machine Learning Research, 24(311):1-34, 2023.
[44] Kirill Neklyudov, Roberto Bondesan, and Max Welling. Deterministic Gibbs sampling via ordinary dif- ferential equations. vrXiv2210.10188 2021.
[45] lain Murray and Lloyd Elliott. Driving Markov chain Monte Carlo with a dependent random stream. arXiv: /204.3187 2012.
[46] Takehiko Morita. Deterministic version lemmas in ergodic theory of random dynamical systems. Hi- roshima mathematical journal, 18:15-29, 1988.
[47] Persi Diaconis and David Freedman. Iterated random functions. SIAM review, 41:45-76, 1999.
[48] Luke Tierney. Markov chains for exploring posterior distributions. The Annals of Statistics, 27(4):1701- 1728, 1994.
[49] Kirill Neklyudov, Max Welling, Evgenii Egorov, and Dmitry Vetrov. Involutive MCMC: a unifying framework. In International Conference on Machine Learning, 2020.
[50] Zuheng Xu and Trevor Campbell, Embracing the chaos: analysis and diagnosis of numerical instability in variational flows. In Advances in Neural Information Processing Systems, 2023.
[51] Luke Tiemey. A note Messopos kemnels for general state spaces. Annals of applied proba- bility, 1998.
[52] Tiange Liu, Nikola Surjanovic, Miguel Biron-Lattes, Alexandre Bouchard-Côté, and Trevor Campbell. AutoStep: locally adaptive involutive MCMC. arXiv:2410.18922 2024.
[53] Miguel Biron-Lattes, Nikola Surjanovic, Saifuddin Syed, Trevor Campbell, and Alexandre Bouchard- Côté. autoMALA: Locally adaptive neetroollis-djjssee Langevin algorithm. In International Conference on Artificial Intelligence and Statistics, 2024.
[54] Nawaf Rab-R, Bob Carpenter, and Milo Marsden. GIST: Gibbs self-tuning for locally adaptive Hamiltonian Monte Carlo. (tti:22004.125333 2024.
[55] Nawaf Rou-Rabce, Bob Carpenter, Tore Selland Kleppe, and Milo Marsden. Incorporating local step-size adaptivity into the No-U-Turn sampler using Gibbs self tuning. arXiv:2408.0822 08259, 2024.
[56] Gareth Roberts and Jeffrey Rosenthal General state space Markov chains and MCMC algorithms. Prob- ability Surveys, 1:20-71, 2004.
[57] Shizuo Kakutani. Random ergodic theorems and Markoff processes with a stable distribution. Pro- ceedings of the Second Berkeley Symposium on Mathematical Statistics and Probability, 1950.
[58] Yuri Kifer. Ergodic theory of random transformations, volume 10. Springer Science & Business Media, 2012.
[59] Radford Neal. MCMC using Hamiltonian dynamics. In Steve Brooks, Andrew Gelman, Galin Jones, and Xiao-Li Meng, editors, Handbook of Markov chain Monte Carlo, chapter 5, CRC Press, 2011.
[60] Simon Duane, Anthony Kennedy, Brian Pendleton, and Duncan Roweth. Hybrid Monte Carlo. Physics letters B, 195(2)211-222, 1987.
[61] Zuheng Xu and Trevor Campbell. The computational asymptotics of variational inference and the Laplace approximation. Statistics and Computing, 32(4):1-37, 2022.
[62] Peter Rossky, Jimmie D Doll, and Harold L Friedman. Brownian dynamics as smart Monte Carlo simulation. The Journal of Chemical Physics, 69(10)44628446 1978.
[63] Samuel Livingstone, Geometric ergodicity of the random walk Metropolis with position-dependee dependent pro- posal covariance. Mathematics, 9(4), 2021.
13

---

[64] Tanja Eisner, Bálint Farkas, Markus Haase. and Rainer Nagel. Operator Theoretic Aspects of Ergodic Theory. Graduate Texts in Mathematics. Springer. 2015.
[65] George Birkhoff. Proof of the ergodic theorem. Proceedings of the National Academy of Sciences, 17 (12):66-6660 660, 1931.
[66] Randal Douc, Eric Moulines, Pierre Priouret, and Philippe Soulier. Markov chains. Springer, 2018.
[67] Yu Qiao and Nobuaki Minematsu. A study on invariance of f-divergence and its application to speech recognition. IEEE Transactions on Signal Processing, 58(7)3388--380, 2010.
[68] Charles C Margossian, Matthew D Hoffman, Pavel Sountsov, Lionel Riou-Durand, Aki Vehtari, and Andrew Gelman. Nested ft: Assessing the convergence Markov chain Monte Carlo when running many short chains. Bayesian Analysis, 1(1):1-28, 2024.
[69] Pavel Sountsov, Colin Carroll, and Matthew Hoffman. Running Markov chain Monte Carlo on modern hardware and software. arXiv:2411.0426, 2024.
[70] Paul Tupper: Ergodicity and the numerical simulation of Hamiltonian systems. SIAM Journal on Applied Dynamical Systems, 4(3):563-587, 2005.
[71] Conor Durkan, Artur Bekasov, lain Murray, and George Papamakarios. Neural spline flows. Advances in Neural Information Processing Systems, 2019.
[72] Matthew Hoffman and Andrew Gelman. The No-U-Tum Sampler: adaptively setting path lengths in Hamiltonian Monte Carlo. Journal of Machine Learning Research, 15(1):1593-1622, 2014.
[73] Augustine Kong. A note on importance sampling using standardized weights. University of Chicago, Dept. of Statistics, Tech. Rep, 348:14, 1992.
[74] Augustine Kong, Jun Liu, and Wing Hung Wong. Sequential imputations and Bayesian missing data problems. Journal of the American statistical association, 89(425):27-225 1994.
[75] Jun Liu. Metropolized independent sampling with comparisons to rejection sampling and importance sampling. Statistics and Computing, 6: 6:113-119, 1996.
[76] Heikki Haario, Ecro Saksman, and Johanna Tamminen. An adaptive Metropolis algorithm. Bernoulli, pages 223-242, 2001.
[77] Radford Neal. Slice sampling. The Annals of Statistics, 31(3):705-785 2003.
[78] Sergios Agapiou, Omiros Papaspiliippoloos Daniel Sanz-Alonso, and Andrew M Stuart. Importance sampling: Intrinsic dimension and computational cost. Staristical Science, pages 405 431, 2017.
[79] Pavel Sountsov. Alexey Radul, and contributors. Inference gym. 2020. URL https://pypi.orz project/inference _gmm.
[80] Wilfred Hastings, Monte Carlo sampling methods using Markov chains and their applications. Biometrika, 57:97-106 1970.
[81] Lars Onsager. Crystal statistics. two-dimensional model with an order-disorder transition. Physical review, 65(3-4), 1944.
[82] Matthew Stephens. Bayesian methods for mixtures of normal distributions. PhD thesis, University of Oxford, 1997.
[83] Edward George and Robert McCulloch, Variable selection Gibbs sampling. Journal of the American Statistical Association, 88(423):881-889 1993.
[84] Kai Xu. Hong Ge, Will Tebbutt, Mohamed Tarek, Martin Trapp. Zoubin Ghahramani. vancedHMC.jl: robust, modular and efficient implementation of advanced HMC algorithms. In Sympo- sium Advances in Appreximate Bayesian Inference. 2020.
[85] Nicolas Chopin, Francesca Crucinio, and Anna Korba. A connection between tempering and entropic mirror descent. In International Conference on Machine Learning, 2024.
[86] Saifuddin Syed, Alexandre Bouchard-Côté, Kevin Chern, and Arnaud Doucet. Optimised annealed se- quential Monte Carlo samplers. arXiv:2408.1 12057, 2024.
14

---

NeurIPS Paper Checklist
1. Claims
Question: the main claims made abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: The main claim of this work a general framework of constructing exact variational families, which is justified with precise theoretical results and empirical demonstrations.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contribu- tions made in the paper and important assumptions and limitations. A No NA question will not be perceived well by the reviewers.
this
The claims made should match theoretical and experimental results, and reflect how much results can he expected to generalize to other settings.
It is fine to include aspirational goals as motivation long is clear that these goals are not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification:
Limitation of the methodology is discussed right after we introduce the method.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are discussed in the paper.
The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point out strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well well-specification. asymptotic approximations only holding locally). The authors should reflect how these as- sumptions might be violated in practice and what the implications would he.
The authors should reflect on a few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
scope the claims made, e.g., the approach was only tested
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low or images taken in low lighting. Or speech-te tochh-tt-ite text system might not be used reliably provide closed captions for online lectures because it fails handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how
they scale with dataset size.
applicable, the authors should discuss possible limitations of their approach address prob-
lems of privacy and fairness.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't It acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparcncy play important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed penalize honesty concerning limitations.
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and com- plete (and correct) proof?
Answer: [Yes]
Justification: We make precise statement of our theoretical results, and defer all the proofs to the Appendix.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross-rfference..
15

---

All assumptions should be clearly stated or referenced in the statement of any theorems.
The proof's can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors encouraged to provide short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Does the paper fully disclose all the information needed to reproduce the main experimen- tal results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: We provide pscudo-code for our algorithm and links to the github that reproduce all experiments in the Appendix.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer this question will not be perceived well by the reviewers: Making the paper reproducible important, regardless whether the code and data are provided or not.
the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For ex- ample. if the contribution is a novel architecture, describing the architecture fully might suffice, or the contribution is specific model and empirical evaluation, it may be necessary to either make possible for others replicate the model with the same dataset, provide access to the model. general. releasing code and data often one good way to accomplish this, but repro- ducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of large language modell, releasing of a model checkpoint, other means that appropriate the research performed.
While NeurIPS does require releasing code, the conference does require all submissions provide some reasonable avenue for reproducibility, which may depend on the nature of contribution. For example
(a) the contribution is primarily new algorithm, the paper should make clear how to reproduce that algorithm.
(b) the contribution is primarily new model architecture, the paper should describe architecture clearly and fully.
(c) If the contribution is way to access this model reproducing results or way reproduce model (c.g., with an open-source dataset or instructions for how to construct the dataset).
model (e.g.. large language model), then there should either
(d) We recognize that reproducibility may be tricky in welcome to describe the particular way they provide for reproducibility. In the case of closed-source models, may be that access to the model is limited some way (e.g.. to registered users), but it should be possible for other rescarchers to have some path to reproducing verifying the results.
cases, which case authors
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described supplemental material?
Answer: [Yes]
Justification: Code and data are provided as Zip file in the supplementary materials; link to the github repo containing all the code and data provided in Appendix.
Guidelincs:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurTPS code and data submission guidelines (httpe::/nips.c// cct.sipps.cppbbbic////cc guides/Code https///oooSuuiiiisooooooooo
for more details.
While we encourage the release of code and data, we understand that this might not be possible, so No is an accepaabl answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
16

---

The instructions should contain the exact command and environment needed to run to repro- duce the results. See the NeurIPS code and data submission guidelines (https://nig) cc/
onPolicy)
more details.
The authors should provide instructions on data access and preparation, including how access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only a subset experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code is permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (c.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: We describe high-level experimental setting in the experiment section (Section 5.), and defer additional details in the Appendix.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to a level of detail that necessary appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, or as supplemental material.
7. Experiment statistical significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate infor- mation about the statistical significance of the experiments?
Answer: [Yes]
Justification: Error bars of each figure is described the caption of the figure or explained in the corresponding text discussing the figure.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" the results are accompanied by error bars, confidence inter- vals, or statistical significance tests, at least for the experiments that support the main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, overall run with given experimental conditions).
The method for calculating the error bars should he explained (closed form formula, call to library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
It should be clear whether the error bar is the standard deviation or the standard error of the mean.
It is OK report I-sigma error bars, but one should state it. The authors should preferably report 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is verified.
For asymmetric distributions, the authors should be careful symmetric error bars would yield results that are out of range (e.g. negative error rates). error bars reported in tables plots, The authors should explain in the text how they calculated and reference the corresponding ligures or tables the text.
show in tables figures
8. Experiments compute resources
Question: For each experiment, does the paper provide sufficient information the computer sources (type of compute workers, memory, time execution) needed to reproduce the experiments?
Answer: [Yes]
Justification: We provide informations of used computational resources Appendix, submitted as PDF in the supplementary materials.
Guidelines:
17

---

The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimen- tal runs as well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (c.g., preliminary or failed experiments that didn't make it into the paper).
9. Code of ethics
Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https: //t////eurri/////eeeeeeeee apeerrip-aa cc/pablicc/tthtcccdl..ee.......e
Question:
Answer:
[Yes]
Justification:
Guidelines:
The answer means that authors have not reviewed the NeurlPS Code of Ethics.
the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
authors should make sure to preserve anonymity (e.g., if there special consideration due to laws or regulations in their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal im- pacts of the work performed?
Answer: [NA]
Justification: This paper focus on theoretical and methodological development in the field of compu- tational Statistics. The societal consequences need not to be dicussed for this work.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No. they should explain why their work has no societal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (c.g., dis- information, generating fake profiles, surveillance), fairness considerations (e.g.. deployment technologies that could make decisions that unfairly impact specific groups), privacy consid- erations, and security considerations.
The conference expects that many papers will be foundational research and not tied to partic- ular applications, let alone deployments. However, if there is a direct path to any negative applications. the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strate- gies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, im- proving the efficiency and accessibility MI.).
11. Safeguards
Does the paper describe safeguards that have been put in place responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: This paper focus on theoretical and methodological development in the field of compu- tational Statisties, which poses no risk for misuse.
Guidelines:
The answer NA means that the paper poses such risks.
18

---

Released models that have high risk for misuse or dual-use should be released with necessary of the model, for example by requiring that users adhere
safeguards allow for controlled usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging. and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Question: Are creators or original owners assels (e.g., code, data, models), used the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected"'
Answer: [Yes]
Justification: We refer all the used models/code in the paper.
Guidelines:
The answer means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (e.g.. CC-BY 4.0) should be included for each asset.
For scraped data from particular source (e.g., website), the copyright and terms of service of that source should he provided.
If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswi thcode datasttts has curated licenses for some datasets. Their licensing guide can help determine the license a dataset.
For existing datasets that re-packaged. both the original license and the license of the derived asset (if it changed) should be provided.
If this information not available online, the authors encouraged reach out the asset's creators.
13. New assets
Question: Are new assets introduced in paper well documented and the documentation provided alongside the assets?
Answer: [Yes]
Justification: Code datasets for reproducing our results is documented submitted a zip as supplementary materials.
Guidelines:
The answer NA means that the paper does not release new assets.
Rescarchers should communicate the details of the dataset/code/mmode as part their submis- sions via structured templates. This includes details about training, license, limitations, etc. used.
The paper should discuss whether how consent was obtained from people whose asset
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and research with human subjects
For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots. if applicable, as well as details about compensation (if any)?
Question:
Answer: [NA]
Justification:
This work does not involve crowdsourcing research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should included main paper.
According to the NeurIPS Code of Ethics, workers involved data collection, curation, other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals equivalent for research with human subjects
19

---

Does the paper describe potential risks incurred by study participants, whether such risks disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equiv- alent approval/reiiew based requirements of your country or institution) were obtained?
Question:
Answer: [NA]
Justification:
This work does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Depending on the country which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and loca- tions, and expect authors adhere the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applica- blc), such the institution conducting the review.
16. Declaration of LLM usage
Question: Does the paper describe usage of LLMs it an important, original, or non -oanstaddard component of the core methods this research? Note that the L.M used only for writing, editing, or formatting purposes and does not impact the core methodology, scientific rigorousness, or originality of the research, declaration is not required.
Answer: [NA]
Justification: The methodological development does not involve LLMs; the only LLM usage is for checking grammer mistakes and typos of the text.
Guidelines:
The answer NA means that the method development in this research does not involve LLMs as any important, original, or non-standard components.
Please refer L.M policy http///aaap...cc//oo/eteeee/e/// should or should not described.
for what
20

---

Contents
1
Introduction
Background
2.1 Homogeneous MixFlows
Involutive MCMC
Iterated random functions
3
Invertible measure-preserving IRF from involutive MCMC
4 Variational flows based on IRFs
4.1 Improved homogencous MixFlows
4.2
IRF MixFlows
6
4.3 Backward IRF MixFlows
4.4 Ensemble IRF MixFlows
7
4.5 Discussion
8
5
Experiments
5.1 Synthetic examples
Real-data experiments
Conclusion
10
Additional content about involutive MCMC
23
A.1 Examples of involutive MCMC
23
A.2 Pseudocode of involutive MCMC
23
B Pseudocode for IRF and inverse IRF based on involutive MCMC
24
Measure-theoretic formulation of pushforward density
25
Proofs
26
D.I Proof of Theorem
26
D.2 Convergence the homogeneous MixFlow
26
D.3 Convergence the IRF MixFlow
27
D.3.1 Convergence the product space
27
D.3.2 From the joint convergence to Theorem 4.2
28
D.4
Convergence the backward IRF MixFlow
................................
29
D.5
Convergence of the ensemble IRF MixFlow
29
D.6 Proof of Proposition 3.1
30
E Additional experimental details
31
E.1
Synthetic experiments
31
E.1.1 Relative performance of homogeneous, IRF, and backward IRF MixFlows
32
E.1.2 Ensemble IRF MixFlows: scaling up M T
33
21

---

E.1.3
Additional results for synthetic examples
34
E.2 Additional results for real-data experiments
35
22

---

A Additional content about involutive MCMC
A.1 Examples of involutive MCMC
Here, illustrate how the generic Metropolis-Hastee Hastinggs (MIT) algorithm [48, 51J, random-walk Metropolis- Hastings (RWMI)) [63, 80), and Hamiltonian Monte Carlo (HMC) [59, 60], into this framework by specity- ing the corresponding auxiliary distribution p(-1z) and the involution map
Example A.1 (MH sampler, Section B.3. of [49]) The Metropolis Hastings sampler with proposal distribution p(dx' (z) can cast an involutive MCMC method by defining the auxiliary distribution as p(dulz) and using the swap involution g (xz, u) + (v, r).
Example A.2 (RWMH sampler; Section 2. of [52]). RWMH with step size is obtained by setting
g(z,
+ ev; eu,,
p(dvla:) (O, I)-
Example A.3 (HMC; [60]). the involutive formulation of HMC, auxiliary variable corresponds to momentum variable, and p(u[x) is the momentum distribution, typically a Gaussian distribution independent of x. The involution map consists of applying ls steps of the leapfrog integrator, followed by a momentun sign flip:
9 -
L
-1
where (x, > (x', ) denotes a single leapfrog step (of step size e) given by
11/2 + v 2 v log)())
z' + EV2/2
v' + 12 v log क (x').
A.2 Pseudocode of involutive MCMC
Algorithm Involutive MCMC kernel K K(x', v' d'|xr
Require: current state z, target n, auxiliary distribution p(dv|z), involution g
1: V ~ p(do(z)
>> sample auxiliary variable generate proposal via the involution
2: (a',')) t g(
3: a min
1, *(=' .0)
Jg [g( m))
> compute the acceptance probability
)(.,))
Accept or reject
4: u Unif 1]
5: if a then
6:
x' r'ta
reject
7: end if
8: return x',
23

---

B Pseudocode for IRF and inverse IRF based on involutive MCMC
Algorithm 2 IRF based on involutive MCMC fo(s)
Require: joint state = (2, D My Ua) random parameters -(e, 8a) update uniform auxiliary variables
1: By <- (up + 0.) mod 1
2: Oa Uua + 0m) mod 1
involutive MCMC with target ((x), auxiliary distribution p(dula)) involution
3: u/ - (t)
4: V + (0(1))(u
5: (z', t g(
6: r - x(2,0) (x, 3)
m(arlar))
>> Compute MH ratio
7: if Uh > r then
8:
return 2,, U, u Ua
reject and return pre-involution state
9: end if
w/g + Mh
>> uu S r implies that u/ E [0. 1]
11: return x', vl M U/a
accept and return after-involution state
Algorithm 3 Inverse IRF based involutive MCMC fū (s')
Require:
joint state s' = (x', v', u',, u(),, random parameters 0 (fu, 0a) recover pre- and post-involution pair g0"", v²)
1: (xr, 7) +-
this will either be in line 6 of Algorithm if accepted, -- otherwise
2: F π(2,8) (2,
[(x',=')
check accept or reject
3: Uu f u/, -
update ua (line 10 Algorithm 2) if the forward pass was an accept
4: if Uu > then
forward pass was a reject (see line 6-7 of Algorithm 2)
-rre-involution state
5: (I,T) t I',
6: Ua +- v/a
7: end if
inverse line 3-4 of Algorithm
8: + F-1 (un)
9: tye < Fp((la)(T)
inverse update of the uniform auxiliary variables (line 1-2 of Algorithm 2)
10: + (up + 1-1 02) mod
Ua + (ua + 0.) mod!
12: return v, U, Ug Ua
24

---

C Measure-theoretic formulation of pushforward density
A fundamental formula when studying variational inference is the the change of variable formula, which char- acterizes the density of transformed distribution. For diffeomorphism X'' + X' on continuous space, the density of X J(Y) Y qo, is given by
Vx E 2, gx(=) fqolz) 90 (x)) J (x)-
J(x) |det Vf()))
However, the assumptions of differentiability and a continuous state space can be restrictive, as many inference problems involve discrete or hybrid spaces (e.g.. Ising model [81], Bayesian Gaussian mixture model [82]. and spike-and-slab model [83]). To handle general state spaces, we adopt measure- measure-theorelic formulation of the pushforward density, stated Proposition C.I, for generic bijection f. This result is well known (see, e.g., 51 for its use in the general involutive MCMC framework), but we include a proof here for completeness.
Proposition C.I. Suppose that f X + X is bijective. For a distribution & n. for all E 2'
d(fg) (x) dm
dg
(5'xx df dx
x)
dx (x).
Proof of Proposition C.I. First, note that if y < n. then fq < J This implies that d(fg)
di (x) = d(fv)
df (z), Vz X.
(z)
dr
dı
remains to show that
d(fg) afx
dr
dg
o
g1 It suffices to show that VA B,
[ dy o - df = / d(S4)
dfr fg(A).
dm
dfm
Note that for all A E B, we have that
J. dq
dq
do (f (f' (()))))( = (()) du
(z)r(dx) = g(5 fc(A),
which completes the proof.
It is worth noting that for Euclidean space X' equipped with the Lebesgue measure and diffcomorphism dfm
J. dm (xr) is precisely the Jacobian determinant |det V (x)].
meassur-ppeeservngg is further x-measure preserving, then dfa density.
yielding a simplified expression the pushforward
dx
Corollary C.2. Suppose that f is bijective and n-measure-preserving For a distribution < n, for all
E X:
a(fg)
dq (S''2). dæ
dæ (x)
Aside from the generality Corollary C.2 over the diffcomorphic case, provides an elegant formula of the pushforward density under a mecsure-preserne map. We invoke Corollary C.2 frequently when developing and analyzing MixFlows.
Beyond extending the diffcomorphic case, Corollary C.2 offers an elegant expression for the pushforward density under measure- preserving measure- map. We frequently invoke this result when developing and analyzing MixFlows. Finally, present specialization of Corollary C.2 for diffeomorphic which provides conve- nient characterization n-meesure-pressevation.
Proposition C.3. Let f f X + X be a diffeomorphism, π be probability distribution on X. with density (denoted by xxx) with respect to duminating measure X. Then,
1. f is trmmassrrepreserving if and only if f1 is m-measure-preserrvng
2. f is tr-maasure-preserving if and only for X-a.e. X, Jy(x) /dt det (z))
s(x) HUF ((r))
Proof of Proposition C.3. By definition, f m-preserving if only if fF(xx) T(x) f 1 1 (xx) Examin- ing the density of the pushforward f via the cf formula, we have
x(x)
Ve E X, (m(() = + (z)) =()))
m(x) < J,(z) (
+x)-
The second claim follows from the fact that = Nz.
25

---

D Proofs
D.1 Proof of Theorem 2.3
As introduced in the main text, the IRF fo induces Markov kernel given by:
Vo E X, VB E B, B)a / (=(f((=)))((d))
This yields simple characterization of action of the Markov process on distribution q: (P4)(u) = / ((r. m)q(dr)= = E [Soq(H)] ~ μs, foq: pushforward of under fo.
We can further characterize the Markov kernel R).:)) induced the inverse IRF So
V VAEB, B, R(x, A) = / 1a(fo (x)))(d))
which is precisely the reversal of P(--)):
P x B) क R( A) / 1 B))(d))) B)p(d)).
(10)
where т @ P(A B) = f4 P(z B)T See Kakutani [57. Eq. (4.5) for the detailed derivation. Notice that P is reversible wrt æ. i.e. π @ P @ R. both the IRF fa and its inverse fø induce the same Markov process P. In other words, = R. From Eq. (10), we can see that sufficient and necessary condition so that
P Q is that
/ ((fo(A) B)u(d)) = (A) (A))Bd)
Proof of Theorem 2.3. From Eq. (4), we see that must admit x as stationary distribution. Douc et al. [66, Theorem 5.2.6] further states that if is the unique invariant probability measure of P, then the Markov process P' ergodic. Therefore, the LLM of ergodic Markov process [66, Theorem 5.29] guarantees Eq. (5), and the random ergodic theorem [58, Cor. 2.2.] ensures Eq. (6).
Then as discussed above, Kakutani [57, Theorem 3.J show that Assumption 2.2 holds for Sa its induced Markov process P if and only if Assumption 2.2 holds for the inverse IRF fø and its induced R. Therefore. the same convergence holds for the inverse IRF
D.2 Convergence of the homogeneous MixFlow
Definition D.1 (Ergodic map 164, pp. 73, 105)), f + X ergodic for ffr all measurable sets X, f(A) = A implies that (AA) E 13.
The most notable implication of a m-e.m.p f is that the long-run average repeated applications of converges to the expectation under ग,, a result known as the Birkhoff ergadic thenrem [65; 64, 212]. The full statement is given in Theorem D.2.
Theorem D.2 (Ergodic Theorem [65; 64, p. 212]). Suppose f X - X is meassare-preserving and ergodic for n. and à L' (T)) Then
1 T
f T- E X.
lim T T0000
Lemma D.3 (Scheffé's Lemma). Let dnn be sequence of integrable funntioon on a measure space (X, B, that convergences w.a.s. to 6. Then
/ ( - ((x))=((x)) +0, n n
if and only if
/ (o)())) ff (o(x)) / n n- oo.
Proof of Theorem 4.1. Note that the Jacobian of the n-e.m.p f is (ss(()//(( (x)) by Proposition C.3, allow- ing the density of qr to be expressed as:
1 T
T
(+())=
Σ lr (x)(
x(x) T
Σ
@
(x)).
Vr
T
t=1
26

---

The pointwise density convergence is the direct consequence of Eq. (11). Specifically, provided 40 < n. we have 40/m (T), the Birkhoff ergodic theorem [65; 64, p. 212] (see Theorem D.2) ensures:
(x) π -- X,
as T +
(11)
total variation convergence is then by the direct application of the Scheffié's lemma Lemma D.3. Notice that
TV m)
/ T (z)
m(dx)
20 (z) x(ds).
To apply Lemma D.3, we set (x) = E.=
d (1o and set o(x) = 1. Because do & n, all on
are T-integrable. Then, for all E N. obtain that
/ (0)())/=(d)))) = Φn (1))
(fo (d))))))))))))))))))))))
π
/ (0( (as for=
(0(=))=(<)),
yielding the second convergence in Lemma D.3.
D.3 Convergence of the IRF MixFlow
As hinted in the main text, the proof of Theorem 4.2 involves interpreting the IRF hom-ehomogeneuus, e.m.p. dynamical system on the joint space @ X. Specifically, we define a map & (Eq. (12)) whose iterates evolve both the state X, and the parameter sequence (0r) Overall, the proof proceeds in two steps. First, we show that the joint law of (0t, Xx) converges in total variation to P @ æ. Second, we deduce marginal convergenee for X,. Scction 1 establishes the joint result, while Section D.3.2 explains why it suffices to prove Theorem 4.2.
D.3.1
Convergence the product space
The key technique for proving the joint convergence is to interpret the iterative process Eq. (2) autonomous, ergodic, measure- preserving dynamical system the joint space e* Given this framework, joint convergence follows immediately, as substantiated by Xu et al. [39, Theorem 4.2) (which is based on the mean ergodic theorem).
For brevity, we define = @º Fr = F 8K and P be the joint distrihution of (0e)ee with independent marginal distribution Define the shift operator a Ω + by
ow (wo, .I....) ++ (wi,,w,
And let (fn)neN be the coordinate process on (2, Fr. P), i.c., for all w = (wo, Wi,...) E S
On (n(()) = Wp.
By definition, we have 0,+1 = 0, o a, and (fen)nEN with (0x)neN μ can be formally understood as (Jen(w)) weN ~ satisfying that Jo,, = Sopoon (w) = Sopcomwy. For the rest of this work, we abuse the notation by writing (o,(w) as far_ all
Now consider the product probability space (52 N.Ja B,P x), where x denotes the joint tribution with independent marginals P and v on Ω and N' respectively. We define the transformation b: Ω x X' + Ω x by
((w, x) = (aw, (x)),, V(u,a) EΩ X.
(12)
Note that Eq. (12) equivalently describes the iterative process Eq. (2) with i..dd (0x)en.. For the rest of the proof, we will focus on the autonomous dynamical system (Ω x X' Fx & B, x n, 4).
Theorem D.4. Under the same assumption of Theorem 4.2, have
1
TV
(Px PXx
+0 as N 00.
(13)
1=1
27

---

Proof of Theorem D.4. We first show that & preserves P PXX, namely, ((Px)) = Pxx. Forall E L L (Pxx).
[ f(w,z)((P x ()
((F (E)
Jax
S (o x a
(14)
1 (t)
Since is meessre-ppreserving for P due to the i.idd. assumption, and + fx) is re-mesuur-pressevving by hypothesis, we obtain that
( x)P)))(()
E(w, ((w.
[ J E(w, ((w))))f=
((w,2))(d))P(d
<(w, (ddo,de)
nxx
(P x (E)
This concludes that (1 x1x X, 2 8, 4) is measure- preserving dynamical system.
We further show that (1 X, Fx @ B, π, 4) is an ergodic dynamical system. Morita [46, Theorem 4.1] shows that it is equivalent to show the ergodicity of the shift dynamical Dystem-(A**, BaN -nn:ruueded by the Markov process associated to Eq. (2). Here Pr is the unique probability measure on (2* _B85) so that the coordinate process (X1, X2....) is Markov chain with kernel P (Eq. (3)) and initial distribution (2,. ). Douc et [66, Theorem 5.2.6] further guarantees that if is the unique invariant probability measure of P, then 1888 T),,T) is both pecasre-ppreservve and ergodic. Hence, the second assertion of Assumption 2.2 guarantees the ergodicity of (12 x X,, Fy & B, क, 4)
æ. and the shift operator on X i... =(.o,
Finally, we apply Theorem 4.2 in al. [39] to finish the proof. Given that measure- measure-preserving and ergodic for P x π. remains to show that < implies that x < P x n. For all B and F E Fix.
0 x x) (F, B) ((F) (B)
P(F) 0 BB) 0.
Since (P x (F, B) = P(F) x q B), PFF = 0, then P(F) x q( 0, and if π = 0, then (B) = 0 by hypothesis and P(F) go B) 0 as well. Therefore, Xu et al. [39, Theorem 4.2] yields the desired result.
D.3.2 From the joint convergence to Theorem 4.2
Finally, we justify why Eq. (13) is sufficient for Fq. (8).
Proof of Theorem 4.2. first derive the explicit expression of 4(P x qo) and examine conditional proba- bility measure. Following the same derivation as Eq. (14), for all & (P x go).
((P go)(E)
[ / Elw, ((c n
(15)
where the second equality is by the fact that is measure-preserving for F. Eq. (15) demonstrates that 4(Px go) can be disintegrated into the marginal distribution P(dw) on Ω and the conditional distribution (daz)) yielding that
X x.I((:) ~ for, fo, for >
where Xo ~ go- Hence, disintegration of N (T x qo) on the slice (0,, 02....) )) is
Σ fan for go.
N
28

---

Then we show that the total variation convergence of the joint distribution (Theorem D.4) implies the total variation convergence of the conditionals (Theorem 4.2). For all N E N,
dP go) d(x
TV
N Σ
(' Px
m(dr))p(dd))
N
Notice that for all N. the Radon-Nikodym derivative and & is P wheasure-preservingg And explicitly, since IP x qu and IP x have same marginal distributions on Ω we have
d#" (xx)) always exists given that IP Qu < P d(2x*)
d (LP x qo) d(P x T)
fa, 2... 10,90
Henee,
N
L
N
fon o .... 90(x)
Σ
(')
4" 4x
=(d=)
T(x)
TV
So o .... So 40:
(0n)Dck ~ IP
Since TV (V(s)) (5-) is always non nog-netaeeve, the left-hand side converges to 0 as N + co yields that the following convergence holds in probability P:
TV
Σ
fon
q0, + 0, as N
N
This completes the proof.
D.4
Convergence of the backward IRF MixFlow
Proof of Theorem The pointwise density convergence the direct consequence of Eq. (9) via Theorem 2.3. The total variation convergence is then established using identical strategy as the proof Theorem 4.1 via Scheffé's ' lemma Lemma D.3.
D.5 Convergence of the ensemble IRF MixFlow
Proof of Theorem 4.4. By the definition of the total variation,
QT M)
TV
4r
(M
(x)
(dw))
go
1o o
of (x) x(dx).
By the triangle incquality,
M
qo
pm))
ao (x)
R²
20 x(dx) + /1 2u (t) (dw) -1 =(dx).
(x)
M
m==1
(16)
We derive upper hounds for two terms on the right-hand side separately.
For the first term, taking the expectation with respect to the randomness of 0 ~ μ., and interchange the order of integrations,
M
E
9o
afmm (ofm) (x)) - R²
(90)
(x) T(dx)
M
mm
Σ
M
q0
mfm) (x)
RD
90
(X) a (1of120))
M
p(m)
Notice that Var E X.
(xr)
RT where the randomness comes from the inde-
((m)
m=1
pendent realization of Os, where the induced the Markov process of fø Therefore, applying Jensen's inequality yields
Vareur
940 (5or
f0r (x)) [x]
VM
L
29

---

For the second term of Eq. (16), since is globally bounded by constant B <
we have that
/ To (v) (d()) -1 x((x)
= // (v) dn (d - x
(v))(dy) T(dx)
- TV(R" dz, ()d(
= E (x N~n.
This completes the proof.
D.6 Proof of Proposition 3.1
Proof of Proposition 3.1. We first verify that map defined Algorithm is fi-meauure- F-mmesune-preserving, invoking the second part of Proposition 3.1. The algorithm has four steps (see Section 3): we compute the Jacobian of each step. Steps 3-4 involve discrete accept/reject decision, branch the transformation is diffeomorphism, making the Jacobian well defined.
treat the two branches separately-witte within
Step describes constant shifts applied to uniform random variables, which preserves (o. and Jnif (dua) with Jacobian 1.
2. Step is the CDF/nnerse-CDD transformation of p(-(z). As long the CDF F(+|z) is well-defined, this step describes diffeomorphism x [0, 1]. The corresponding Jacobian given by:
p xr)
p(v
3. We analyze step 3 and together. In the rejection branch, additional transformation is applied, the Jacobian In the acceptance branch, step involves the involution mapping, with Jacobian
Og(n, and step 4 rescale by the MH-ratto r, yielding combined Jacobian with step
[(,,6)
Henec, in the rejection branch, the combined jacobian of step I evaluated on s' = (xx, 3. u', Ma) is
p x)) (2, v. uhy Ma)
p (2, v, uo u.)
acceptance branch, combined jacobian step 1-4 evaluated on s' (5', thes ul))
p(ū(z) म "(r',') म F(x', ('a',"', uhe ula) p(u(z) ((x 0)
x <(, US Mey ua)
Both satisfy the criterion of Proposition 3.1; 1; the map is therefore F-measure- prresure-reeservigg
Finally, we show uniqueness invariant distribution. By Douc et al. [66, Corollary 9.2.16), an irreducible kernel has at most one invariant distribution. Because each S' preserves म. the induced Markov kernel P must admit I invariant distribution. If P is irreducible, then is its unique invariant distribution.
30

---

E Additional experimental details
For all homogeneous MixFlows variants, the uniform- shift parameters were fixed to 0, = and 0. x/7. For NUTS benchmarks, normalizing flow architectures were implemented as follows. In RealNVP, the affine coupling layers consist of two separate multilayer perceptrons (MLPs)-O0 one for scaling and one for shifting- each with three fully con- nected layers and LeakyReLU activations. For Neural Spline Flows (NSF), we set the spline bandwidth to B = 30, and used K 11 knots. For synthetic examples, the hidden dimension in each MLP was set to 32 for RealNVP and 64 for NSF. For real-data examples, the hidden dimension set min(d, 64)., where d the dimensionality of target posterior distribution.
the Julia package AdvanCedHMC., ji [84] with default settings throughout. The
Experiments are conducted on following platforms: a local machine equipped with an AMD Ryzen 5900X CPU and 64 GB of RAM, the ARC Sockeye computing platform at the University of British Columbia, and the neh-performance compute cluster provided by the Digital Research Alliance of Canada. Code for reproducing the main experimental results is available https: htpp:://gghhul com/ zuhengxu/Mixxlor. j1 git.
E.1 Synthetic experiments
The four target distributions used in this experiment are as follows:
1. the banana distribution [76]:
y1
N
0
100
9)
y1
== 0.1;
y2
822 + byi 1006]
2. Ncals' funnel [77]:
21 N (0,0²),
22
x2 N (0, exp (2))
o² 36;
cross-shaped distribution: in particular, a Gaussian mixture the form
0
o..15²
0
2
x 7
W
0
1
+
N
0
0.15²
N
2
0.15²
0
+
1
W
([-2
0
0.15
0
and warped Gaussian distribution
y1
0
[llll2 cos (atan2 (y2, v1)
1 [lvlla)
=
y2
0.123
[ivllasin (atan2 v1)
_ [|v||==)
where atan2(y, x)) the angle, in radians, between the positive x axis and the ray to the point (x. y).
31

---

E.1.1 Relative performance of homogeneous, IRF, and backward IRF MixFlows
10-
10-
35-
14-
14-
3.-
1006
2300 3800
1004
2500
1000
2:00
3B0
2540
Total-variation error for homogeneous, IRF, and backward IRF MixFlows built on RWMH kernels, plotted against flow length T for the most performant step sizes among {0.05, 0.2, 1.0}. Each curve is the mean over 32 independent runs; shaded bands (#1 SD) show run-to-run variability. Crows RANNN pord homo
Crovs wwwer Lie
Crous_ Fww.H bood io
Cross RWWH mmome
Meon
Meen
Meen
Iteration
Iteration
Iteration
42
Mean
Mean
Mean
Iteration
Iteration
Iteration
Iteration
Mean
Mean
Iteration
Iteration
Iteration
Iteration
Running mean estimates over 3000 iterates from different IRF and MCMC dynamics based on RWMH, eval- uated on the Cross distribution across 32 independent runs. Each line represents the trajectory of a single run. From top to bottom, the rows show the running mean of the test functions (21, + 3ı, (21, 22) + 22, and (211, x22) + dh (x1112) From left to right, the columns correspond to the dynamic of inverse IRF fz the backward process of the inverse IRF, hime-homogeneuus dynamics for, and the standard RWMH MCMC.
Per-sample MCMC effective sample size (ESS) estimates on the test function computed from trajectories generated by various IRF and MCMC dynamics based on HMC, MALA, and RWMH kernels. The trajectory lengths are set to 300 for HMC-based dynamics, 2000 for MALA, and 4000 for RWMH. Each ESS value is computed from single trajectory, and the boxplots summarize the ESS estimates over 32 independent runs per method. The per-sample ESS for i.i.d. samples will be 1.
Figure 5: Results showing difference between homogeneous, IRF, and backward IRF MixFlows
Fig. 5a compares the total variation (TV) errors of homogeneous, IRF, and backward IRF MixFlows con- structed from RWMH kemelss. Overall, homogeneous and backward IRF MixFlows perform similarly, though the latter exhibits slightly improved accuracy longer flow lengths. IRF MixFlow consistently outperforms both, achieving faster TV convergence and lower variability across runs. As discussed in Section 4.5, this im- provement stems from differences in the convergence hehavior of the series π >> 20 (Ta(z)), where Tk represents the sequence of transformations used in the density computation each MixFlow variant.
Fig. 5h further illustrates this effect by showing running mean estimates over 3000 iterations for the Cross distribution. From top to bottom, each row shows the mean of the test functions (I1, + 21, (I), X2) + 22. and (1)x))) - 20 (r). From left right, the columns correspond to the inverse IRF fo (backward IRF MixFlow), the backward process of the inverse IRF (IRF MixFlow), the time- homogeneous flow fa- (homogeneous MixFlow). and standard RWMH MCMC. The backward process exhibits significantly faster convergence in all cases, consistent with the superior TV performance of IRF MixFlows under equal flow lengths. This advantage arises from reduced autocorrelaton in the hackward iterates.
32

---

Fig. Sc reports the per-sample MCMC effective sample size (ESS) for the test function estimated from trajectories generated using various IRF and MCMC dynamics based HMC, MALA, and RWMHL This metric captures degree autocorrelation
2u (2)) across iterations. Backward process dynamics consistently yield ESS values orders of magnitude higher than other methods often approaching the ideal of independent sampling, with relative ESS close to in some cases.
E.1.2 Ensemble IRF MixFlows: scaling up M or T
lmnghh
Fox langth
langt
Fix ensemble size M 30 increase flow length.
Siep V422
stepsiee
SIED
Number Couatiie
Number Cowembbe
Number
Mamber Coaer bin
Fix flow length 100 increase ensemble size.
Figure 6: TV eтror of ensemble IRF MixFlows based on HMC over increasing ensemble size M and flow length T. Each curve is the mean over 32 independent runs; shaded bands (#11 SD) show run-to-run variability.
33

---

E.1.3 Additional results for synthetic examples
0.0
FLA0
LLNU
-03
-30
HNC
-20
HMC ae
-30
Scrticn wndd
(a)
ELBO: (from left right) funnel, cross, warped Gaussian.
23r17"
(b) ELBO in full range: (from left to right) Banana, funnel, cross, warped Gaussian.
normalization constoin
log normaliation constant
log normalication constant
-5-
0.75
2.5
0.50
0.25
HWli logo
Trm
FELEst sperie
0.00
She Egg
(c) log Z estimates: (from left to right) funnel, cross, warped Gaussian.
kg
ESS per somple
ESS per sample
16"1
ESS
10
FINCH Beffect camg
perfect comppl
(d) Per-sample importance sampling ESS: (from left to right) funnel, cross, warped Gaussian.
Figure 7: Variational approximation quality of IRF Flows versus RealNVP and NSF. Box plots for IRF flows are based on 32 independent runs, and 10 runs for the normalizing flows.
34

---

E.2 Additional results for real-data experiments
To approximate the ground truth, ran AIS procedure with 1096 particles with adaptive schedule selection. The initial temperature schedule generated via mirror descent [85] with small step size of 0.005: the schedule was then refined for five rounds using adaptive scheme Syed et al. [86]. yielding more than 1000 annealing steps for each data set. All reference values are taken as the median estimates across 10 independent runs of the above procedure.
134.5
-2450
ELBO in full range.
Per-sample importance sampling ESS.
FENOOT caor gou orrer
840
randoo cocoe
Maximal absolute error the coordinate mean estimation against NUTS
Figure 8: Results real-data benchmarks (columns, from left to right); TReg(d 4).
Brownian(d = 32), SparseReg (d = 83), and LGCP (d = 1600)
35


Reviews:
[{"review_id": "8kmshy4TJN", "reviewer": "Reviewer_UFbd", "summary": "The paper extends the recently introduced class of Mixflows algorithms - to construct variational approximations via ergodic measure preserving transformations - with non-deterministic maps. This allows for some results in terms of total variation convergence to the target distribution. They compare against normalizing flows and MCMC on estimating log normalization constants and certain posterior moments respectively.", "questions": "Additional comments/questions/suggestions for the authors:\n\n- MixFlows and their extension here seem (to me, at least) clearly closely connected to adaptive importance sampling - which has a relatively well-developed literature - where exact density evaluations are provided and final estimators are often able to use all generated samples during the learning process. It would be good to comment on any connections that the authors see or do not see. Adaptive IS literature dates back to at least Pugh (1966), with other common references being Oh et al (1992), Bugallo et al (2017), Portier et al (2018). Of particular relevance here is that there exist adaptive IS algorithms where the proposal learning is driven by MCMC (Martino et al (2017)). Perhaps the authors could expand as it is sometimes possible to obtain central limit theorems for AIS estimates (see Portier et al (2018)), in contrast to estimates with MixFlows it seems. \n- As it seems that all generated samples can be used in the final estimation here, it would be interesting to explicitly write down the Z (or log Z) estimate with importance sampling \n- An empirical comparison with tempered SMC with HMC/NUTS kernels, seen by some as a very good reference for Z estimation, may strengthen further the relevance of the proposed approach ( default implementations exploiting GPUs are available in modern Python libraries, such as BlackJax), and commenting the advantages of the proposed approach (for example, not having to choose a sequence of distributions as in SMC). It may also be worth noting that SMC has been viewed also an adaptive importance sampling algorithm with intermediate target distributions. \n- Aside of / instead of Theorem 2.2, I would have expected a theorem for the convergence of the normalizing constant estimates Z (or log Z thereof). After all $q_T$ for any $T$ does provide exact density evaluation, therefore any valid Z estimate would be naturally based on importance sampling, and this is what is done in the experiments as I understand. Further, the variational family would provide consistent estimates of expectations of integrable test functions even for finite $T$ when IS is used. \n- In practice, is it easier to tune the MCMC kernels in your method than it is to tune the hyperaparameters of the normalizing flows training ? As in a real world scenario, that (difficult to quantify) cost may be also taken into account. \n- Using the ELBO metric seems somewhat redundant (but perhaps the authors can explain) with estimation of log Z with importance sampling - as clearly the ELBO is intractable and has to be MC estimated anyway, so one may as well just refer to the estimation of log Z (which is of course a biased but consistent estimator the true log Z, with an MSE relating also to KL between target and proposal). \n- Is Theorem 3.1. a convergence in distribution ? Please state it clearly.  \n\n\nReferences: \n- Pugh EL. A gradient technique of adaptive Monte Carlo. SIAM Review. 1966 Jul;8(3):346-55.\n- Oh, M.S. and Berger, J.O., 1992. Adaptive importance sampling in Monte Carlo integration. Journal of statistical computation and simulation, 41(3-4), pp.143-168.\n- Bugallo, M.F., Elvira, V., Martino, L., Luengo, D., Miguez, J. and Djuric, P.M., 2017. Adaptive importance sampling: The past, the present, and the future. IEEE Signal Processing Magazine, 34(4), pp.60-79. \n- Portier, F. and Delyon, B., 2018. Asymptotic optimality of adaptive importance sampling. Advances in neural information processing systems, 31.\n- Martino, L., Elvira, V., Luengo, D. and Corander, J., 2017. Layered adaptive importance sampling. Statistics and Computing, 27, pp.599-623.", "limitations": "The authors have a fairly honest presentation and discuss limitations.", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 3, "significance": 4, "originality": 4, "strengths_and_weaknesses": "Strengths: the paper provides a nice continuation of the Mixflows research line with a nontrivial and interesting extension of this class of algorithms. The paper is largely clearly written, and the mathematical notation is generally precise. Experiments touch realistic posterior distributions and carefully take into account total computation cost (for instance, of training the neural networks in normalizing flows). \n\nWeaknesses: Re-emphasizing this is a good paper, there can always be some reservations. Some parts remain partially unclear, although this may be a limitation of my understanding of iterated random function theory and/or Markov chain theory. In 3.2, is the simulation of the backward process required ? Reference [34] is cited, but it seems overkill to have to check it. Section 4 is also quite dense and hard to follow. Further, concrete examples of four MCMC-based IRFs arrive only in page 7. I recommend concrete instances to come earlier in the presentation of the IRF framework, and describe in a concrete way how they can be (approximately) inverted. \nFurther, the comment starting line 230: \"Notably, the corrected HMC IRF is consistently more stable than its uncorrected counterpart used in past MixFlows work, because the additional MH accept reject step discards trajectories with large numerical error that would otherwise cause the dynamics to diverge.\" seems a bit directly at odds with the comment immediately after in line 235. Perhaps the authors can explain better here. Finally, the assumption in Theorem 3.3 of bounded (reverse) importance ratio with respect to the initial distribution seems strong, although I may be wrong. Could you comment in general about the importance of the choice of initial distribution ? (For instance, in SMC it matters quite a bit).", "was_revised": true, "final_justification": "See my full comment below (which I wrote before realising I'd have to write here as well).", "rebuttal": "Thank you for your insightful comments and for acknowledging our work! Please see our point-to-point response below.\n\n> **Weaknesses**: Re-emphasizing this is a good paper, there can always be some reservations. Some parts remain partially unclear, although this may be a limitation of my understanding of iterated random function theory and/or Markov chain theory.\n\nThank you again for your acknowledgement and pointing out texts that deserve further clarification! We would revise the manuscript in camera-ready for a clearer explanation.\n\n> In 3.2, is the simulation of the backward process required ? Reference [34] is cited, but it seems overkill to have to check it. \n\nWe will state this more clearly in revision. For the backward IRF MixFlow, you do not need to simulate the full backward process for either independent sample generation or density evaluation. Because $q_t$, drawing a sample simply involves selecting one mixture component uniformly at random, and applying that component’s forward map, avoiding simulating the whole backward process.\n\n> Section 4 is also quite dense and hard to follow. Further, concrete examples of four MCMC-based IRFs arrive only in page 7. I recommend concrete instances to come earlier in the presentation of the IRF framework, and describe in a concrete way how they can be (approximately) inverted. \n\nThank you for the suggestion. In the revised manuscript, we will reorganize Section 4 to introduce MCMC‑based IRFs up front and include concrete examples—such as RWMH and HMC—to illustrate how these involutive MCMC algorithms admit an IRF representation and can be (approximately) inverted.\n\n> the comment starting line 230: \"Notably, the corrected HMC IRF is consistently more stable than its uncorrected counterpart used in past MixFlows work, because the additional MH accept reject step discards trajectories with large numerical error that would otherwise cause the dynamics to diverge.\" seems a bit directly at odds with the comment immediately after in line 235. \n\nThank you for highlighting this point! We will edit this paragraph in revision to resolve the confusion. To clarify, divergent chains are more unstable and exhibit drastically higher numerical inversion error. Although shadowing property does reduce the impact of such numerical error significantly, applying MH correction to filter out the unstable trajectories keeps the computation stable from the outset, so that all downstream inferences remain reliable without relying on shadowing.\n\n> the assumption in Theorem 3.3 of bounded (reverse) importance ratio with respect to the initial distribution seems strong, although I may be wrong. Could you comment in general about the importance of the choice of initial distribution ? (For instance, in SMC it matters quite a bit).\n\nThis is a good question! We would like to clarify that the bounded density-ratio condition is not required to prove the asymptotic exactness of the ensemble MixFlow. For Theorem 3.3, however, we opted for a statement that intuitively highlights how the flow length $ T $ and the ensemble size $ M $ presents the bias-variance trade-off, rather than the tightest possible bound or the weakest assumption.\n\nBroadly speaking, the role of the initial distribution in MixFlow parallels its role in SMC or general MCMC: the closer $ q_0 $ is to $ \\pi $, the fewer IRF iterations (T) are needed to reach a given approximation quality. Although our current theory establishes asymptotic exactness, it does not yet quantify the impact of $ q_0 $’s quality. We plan to conduct a more refined theoretical analysis that explicitly addresses the significance of the initial distribution. Thank you again for your insightful feedback!\n\n\n> **Questions**: MixFlows and their extension here seem (to me, at least) clearly closely connected to adaptive importance sampling...\n\nThank you for highlighting the rich literature on adaptive multiple importance sampling—we weren’t previously aware of those specific references, and there are indeed many intriguing connections to our work.\n \nAt a high level, both MixFlows and adaptive IS methods build a mixture approximation to the target distribution using proposals informed by MCMC samples. To our understanding, in adaptive IS (e.g., the Markov adaptive IS of Martino et al., 2017), one typically constructs a kernel‐smoothed density around a growing population of MCMC iterates, tuning kernel bandwidths via intermediate importance weights. By contrast, MixFlows directly parameterize a sequence of flow transformations using the Markov chain transitions themselves. This distinction leads to different theoretical guarantees: MixFlows enjoy pointwise density convergence and total‐variation convergence under minimal assumptions (thanks to their asymptotically exact nature), whereas adaptive IS approaches generally require additional conditions—such as an infinite particle limit and vanishing kernel scales—to attain similar convergence results.\n\nThat said, adaptive IS does offer valuable ideas, particularly the use of intermediate weights and particles to adapt proposal hyperparameters. Incorporating a similarly robust adaptation mechanism for MixFlow kernel or step‐size parameters is an exciting direction for future work, and we will discuss this connection and potential extensions in the camera‐ready version.\n\n> As it seems that all generated samples can be used in the final estimation here, it would be interesting to explicitly write down the Z (or log Z) estimate with importance sampling\n\nGood point! We realize that we didn’t explicitly write out our Z estimates and we will make sure to clarify this in revision. For the current version, we just take iid samples from the MixFlow and do the following naive IS estimator based on iid draws: \n$$ \nZ \\approx \\frac{1}{N}\\sum_{n=1}^N  \\frac{\\pi}{q_T}(X_n), \\quad (X_n)_{n=1}^N \\stackrel{iid}{\\sim} q_T\n$$\n\nHowever, we are aware that more sophisticated estimates can be obtained by taking advantage of the pushforward structure. For instance, for the deterministic MixFlow, we can use the trajectory averaged IS estimator (as noted in Algorithm 3 in [26]) that leverages all intermediate samples. We can also do the PISA estimator as investigated in [Andrieu et al, 2025] to use all intermediate samples. In this work, we try to restrict the scope to the construction of the variational family to avoid overcomplicating the narrative. A more comprehensive investigation of how to do estimates with the family is left for future work.\n\n> An empirical comparison with tempered SMC with HMC/NUTS kernels, seen by some as a very good reference for Z estimation, may strengthen further the relevance of the proposed approach...\n\nThank you for the suggestion. We feel it may not fully clarify the relationship between our approach and SMC/annealed IS. On one hand, ensemble MixFlow does resemble annealed importance sampling—both propagate particles via parallel MCMC chains and assign weights, with performance driven by how well particles explore the target.  On the other hand, our method uses a homogeneous Markov kernel throughout, whereas SMC/annealed IS relies on inhomogeneous kernels guided by an annealing path.  This key distinction makes it difficult to disentangle the effect of the annealing path, and thus can make direct performance comparisons less informative. Therefore, understanding the connections to SMC/annealed IS requires a more careful investigation. \n\nThat said, we are actively working on a follow‑up work that incorporates annealing and establishes a direct connection to annealed IS methods. Stay tuned!\n\n> Aside of / instead of Theorem 2.2, I would have expected a theorem for the convergence of the normalizing constant estimates Z (or log Z thereof). After all  for any  does provide exact density evaluation, therefore any valid Z estimate would be naturally based on importance sampling, ...\n\nYes, you are totally correct that all standard guarantees of importance sampling would apply to MixFlows. As aforementioned, we are working on a followup work for a more comprehensive theory analysis and only keep theoretical results justifying the asymptotic exactness of MixFlows in this manuscript.\n\nWe’d like to clarify that the purpose of Theorem 2.2 is to introduce the standard ergodic theorem of the IRF on any L1(pi) functions. This serves as the intuition on why we have pointwise density convergence of IRF MixFlows. For example, in eq(5) setting $ \\phi $ (in the statement of Thm 2.2) to be $ \\frac{q_0}{\\pi} $ shows that $ \\lim_{T \\to \\infty} \\overleftarrow{q_T}(x)= \\pi(x)$. \n\n> Using the ELBO metric seems somewhat redundant (but perhaps the authors can explain) with estimation of log Z with importance sampling...\n\nThis is a fair point. We keep ELBO results in the main text mostly because it's the standard metric that’s used in the flow and VI community. Depending on the space, we might consider moving some ELBO results to Apdx in camera-ready.\n\n> Is Theorem 3.1. a convergence in distribution ? Please state it clearly.\n\nIt is a  convergence in probability (with respect to the probability measure $ \\mathbb{P} $). We will clarify this in camera-ready."}, {"review_id": "XKQyYH75QO", "reviewer": "Reviewer_xH4x", "summary": "The paper proposes a variational family for performing variational inference in Bayesian models, designed to have the theoretical convergence guarantees of MCMC. The key idea is to construct the variational family by generalising mix variational flows, which use Markov kernels to approximate the target distribution. The authors focus on involutive kernels and, to improve expressiveness and theoretical guarantees, the paper composes multiple involutive kernels using iterated random functions (IRFs), making random mix variational flows. These random mix flows relax the ergodic condition for convergence of the flow. Moreover, the method uses the backward process of the IRFs for training, arguing that this reduces computational cost. They further enhance their approach by proposing an ensemble of mix variational flows based on IRFs, and they provide theoretical results suggesting this ensemble can accelerate convergence. They introduce invertible versions of the kernel compositions and evaluate their performance on both synthetic and real datasets.", "questions": "Can you explain exactly how is it that the inverse of the backward IRF MixFlow has lower computational cost than the inverse of the forward IRF MixFlow?\n\nHow do the data benchmarks and computational cost of experiments compare to MALA or RWMH instead of NUTS? Maybe the computational cost figures would be clearer without the baseline. Maybe also present wall clock time?\n\nWhat is the computing cluster used for experiments? is it a cluster of CPUs or GPUs? Using GPUs would parallelise many of the gradient and density evaluations of the computing cost. It would be interesting to know the wall clock time of training and inference of your method in comparison with normalizing flows when using parallel computing.\n\nWhat are the computational costs of sample generation? especially comparing random mix flows with normalizing flows, or MCMC.", "limitations": "There is no discussion on the limitations of the method. More about the computational costs, which are clearly expensive for an ensemble method but at what gain in expressiveness/convergence? and the limitations on architecture design and how this affects the outcome.", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Quality:\nThe paper is technically sound and builds on a foundation of existing theory, the proposal of using involutive kernels and their composition through IRFs is well motivated for the task of building variational families with weaker conditions for convergence. The practical algorithm presented for invertible measure-preserving IRFs makes the method applicable beyond theory. However, the work feels somewhat incomplete in places, particularly around the justification for using the backward IRF process, which is claimed to reduce computational cost but lacks convincing explanation or empirical evidence. The experiments do not fully explore the implications or variations of the proposed method, especially regarding architectural choices and computational costs. The authors acknowledge some limitations but could be more thorough in critically evaluating different design choices.\n\nClarity:\nThe paper is generally well organised and clearly written. However, the rationale and benefits of using the backward process of IRFs instead of the forward require better explanations. The limited discussion on implementation details, such as architecture choices and computational trade-offs for ensemble methods, reduces clarity for readers looking to apply the approach or reproduce results.\n\nSignificance:\nBy linking MCMC kernels, IRFs, and variational inference, the paper opens a promising direction for constructing flexible yet theoretically grounded approximations. The potential impact lies in enabling practitioners and researchers to build variational approximations that converge in total variation to the true posterior, addressing a known limitation of many variational approaches. The open questions regarding computational efficiency, backward vs. forward IRF usage, and ensemble trade-offs leave open how impactful the method will be in practice.\n\nOriginality:\nWhile the paper largely combines and extends several existing ideas, it does so in a novel and meaningful way. The framing of variational families as compositions of involutive kernels via IRFs, particularly using the backward process, is an original contribution that advances understanding of how to construct variational approximations. The ensemble approach to accelerate convergence also adds novelty.", "was_revised": true, "final_justification": "While I believe the experimental section could be strengthened, the core methodological contribution of the paper is both original and well executed. The proposed framework is meaningful in constructing asymptotically exact variational approximations, and for that reason, I believe the paper should be considered for acceptance.", "rebuttal": "Thank you for your review and insightful comments! Please see our point-to-point response below.\n\n> Summary: ... Moreover, the method uses the backward process of the IRFs for training, arguing that this reduces computational cost...\n\nWe hope to clarify that we didn’t use any of the IRF for **training**. The IRF backward process is used to define the backward IRF MixFlow family. We remark that all variants of MixFlows are asymptotically exact regardless of tuning.\n\n> Question: Can you explain exactly how is it that the inverse of the backward IRF MixFlow has lower computational cost than the inverse of the forward IRF MixFlow?\n\nTo make this precise, the density computation of the backward IRF MixFlow (Eq(5)) has lower computational cost than the density computation of the forward IRF MixFlow (Eq(9)). \nAs explained in line 124-126, the latter requires simulating the backward process of the inverse IRF $ \\overleftarrow{X_t}(x):=f_{\\theta_1}^{-1} \\circ \\cdots \\circ f_{\\theta_t}^{-1}(x) \\quad \\text { for } t \\in[T]$. Notice that we cannot compute this sequence sequentially—$ \\overleftarrow{X_t}(x) \\neq f_{\\theta_t} \\circ \\overleftarrow{X_{t-1}}(x) $. Each $ \\overleftarrow{X_t} $ has to be computed independently,  taking $O(T^2)$ operations to compute $ \\\\{\\overleftarrow{X_t}(x): t\\in [T]\\\\} $. \n\nIn contrast, the density computation of the backward IRF MixFlow (eq(9)) requires simulating $  \\\\{f_{\\theta_t}^{-1} \\circ \\cdots \\circ f_{\\theta_1}^{-1}(x): t \\in [T] \\\\}$. Notice that the indexing is different from the aforementioned backward process. This process can be computed sequentially by applying $f_{\\theta_t}$ in order, taking $ O(T) $ operations.\n\nWe think this explanation suffices to justify the difference of the computational cost of the two IRF MixFlows, and are happy to make further clarification in camera-ready. \n\n> Quality: The experiments do not fully explore the implications or variations of the proposed method, especially regarding architectural choices and computational costs. The authors acknowledge some limitations but could be more thorough in critically evaluating different design choices.\n\nWe appreciate your suggestion and fully agree that a rigorous understanding of the relative strengths and weaknesses among the four MixFlow variants is important. As noted in the main text (lines 177–183) and in Appendix E1.1.1–E1.1.2, we have already included empirical comparisons and discussions around the different MixFlow members and design choices. That said, we are happy to expand on these analyses in the camera-ready version if space permits.\n\nWe would also like to clarify that a comprehensive investigation into the design choices and robust implementation strategies for MixFlows—such as more sophisticated tuning strategies for the MCMC kernel, and principled selection of M and T for the ensemble Mixflows—entails deeper methodological work. These are active research areas in Bayesian computation (e.g., [38-41, 47,48]) and are beyond the scope of the current submission. Our main goal in this work is to establish the foundational MixFlow framework and generic recipes of designing its flow maps so that we can have room for further refinement in future work.\n\n> Clarity: the rationale and benefits of using the backward process of IRFs instead of the forward require better explanations. The limited discussion on implementation details, such as architecture choices and computational trade-offs for ensemble methods, reduces clarity for readers looking to apply the approach or reproduce results.\n\nPlease see our repsonse to the previous comment.\n\n> How do the data benchmarks and computational cost of experiments compare to MALA or RWMH instead of NUTS?\n\nHMC (NUTS) is known to outperform MALA in most settings (see, e.g., Chen and Gatmiry, 2023), so we do not expect MALA to offer more competitive performance on our benchmarks. Likewise, while RWMH may incur lower computational cost, it generally produces less accurate posterior estimates than NUTS.\nMore importantly—as we note on lines 300–301—our goal is not to establish superiority over MCMC methods. We use NUTS as a benchmark to demonstrate that our asymptotically exact flows achieve comparable accuracy in estimating posterior summaries, while retaining the tractability of variational inference methods. We believe this choice of experiments is sufficient to support our claims.\n\n> Maybe the computational cost figures would be clearer without the baseline. Maybe also present wall clock time?\n\nThis is a good suggestion! We now have a figure comparing the wall clock computation time of all the methods, which shows a similar phenomenon as the current Figure 4(d). For example, for the Brownian example, the wall clock time (in seconds) for each methods are shown as follows:\n| Deter/Bwd IRF MixFlow | Ensemble MixFlow    |NSF     |  RealNVP |  MFVI   |  NUTS   |\n| :-------- | :------- | :------ | :-------- | :------- |:------- |\n| ~ 2  |  ~14  | ~2600     |  ~ 784|    ~15 |  ~27   |\n\nIn camera-ready, we will update Fig 4(d) with the wall clock time.\n\n\n> What is the computing cluster used for experiments? is it a cluster of CPUs or GPUs? Using GPUs would parallelise many of the gradient and density evaluations of the computing cost. It would be interesting to know the wall clock time of training and inference of your method in comparison with normalizing flows when using parallel computing.\n\nWe conducted all experiments using CPU clusters, consistent with the wall-clock times we linked above. Detailed specifications of the computing clusters will be provided in the camera-ready version. Regarding the computational efficiency comparison between MixFlows and Normalizing Flows (NFs), we would like to clarify that our work does not claim superior computational efficiency for MixFlows. Indeed, NFs can significantly benefit from GPU-based parallelization. However, we emphasize that MixFlows—particularly ensemble MixFlows and IRF MixFlows as highlighted in our paper—can also leverage similar parallelizations.\n\nCrucially, our primary argument is centered around the fundamental difference between these two families of flows: MixFlows are asymptotically exact, whereas NFs are not. As a result, NFs typically require careful architecture selection and extensive hyperparameter tuning to avoid unreliable results or training failures. For instance, as demonstrated in Section 5.2 (Figure 4, third column), the training of RealNVP consistently diverged, returning NaNs across all configurations.  Similar challenges with NFs were also noted in the original MixFlow study (Section 6.3 of [26]) and in the RealNVP paper itself (Section 3.7 of [19]). \n\nIn contrast, tuning of MixFlows is much easier as they are parameterized via MCMC kernels whose hyperparameters (such as step sizes) are more interpretable and mostly low-dimensional. We can simply tune these parameters based on MH acceptance rate (what we did in our experiments). For all our real data experiments, the step size tuning of RWMH based MixFlows is completed within 20 seconds to 1min, while NF training typically takes several hours to more than a day on the same hardware.\n\n\n> What are the computational costs of sample generation? especially comparing random mix flows with normalizing flows, or MCMC.\n\nIf using the same involutive MCMC kernel, sample generation of the IRF MixFlows (simulating the forward process of IRF) is almost identical as running the MCMC iterations, but with the additional capabilities of weighing the samples based on the importance sampling weights with respect to the target.\n\nAfter training, the sample generation of NFs is typically much faster than MixFlows. But as we mentioned above, NFs come with significant initial tuning time and quality of samples relies on careful tuning. For example, in the Brownian example, by the time we have trained RealNVP, we have already run over 10million MCMC iterations (RWMH kernel)."}, {"review_id": "5QOjBdsMRf", "reviewer": "Reviewer_UHcw", "summary": "The authors introduce a IRF-based method for asymptotically exact variational inference that exempts the ergodicity assumptions from prior approaches. Towards this objective, they develop an involutive (self-invertible) MCMC kernel that parameterizes a measure-preserving IRF, which is shown to be invertible (Section 4) and asymptotically exact (Proposition 4.1). An ensemble algorithm with potentially better convergence rates is also proposed (Section 3.3). \n\nIn summary, the paper is well organized and nicely written, though somewhat dense. The contributions are clear and stated in a mathematically precise manner. However, the empirical analysis could be more comprehensive and rigorously executed.", "questions": "See weaknesses. \n\nAlso, traditionally, variational methods are characterized by finding the best tractable approximation to an intractable target distribution via gradient-based search. IRF-based approaches, however, take a different path. Could the authors discuss the connection between the proposed methodology and traditional VI?", "limitations": "Limitations are properly addressed and future directions, e.g., designing automatically tuned and asymptotically exact variational methods, are concretely described.", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "## Strengths \n\n1. The paper is theoretically sound and rigorously built upon the well-founded ergodic theory.\n\n2. The work’s contributions and the limitations of its predecessors are clearly laid out and revisited throughout the paper. \n\n3. Empirical results suggest that the proposed approach converges faster than traditional flow-based methods - while being asymptotically exact. \n\n4. Computer code has been released for improved reproducibility. \n\n5. Proofs for all statements and a simpler proof for the asymptotic exactness of MixFlows have been included in the appendix.\n\n6. The proposed algorithm is, to the best of my knowledge, the first asymptotically exact variational method that doesn’t explicitly require ergodicity. \n\n## Weaknesses\n\nMy concerns are mostly concentrated on the paper’s limited empirical analysis and on isolated claims. I will be happy to increase my score if the issues below are properly addressed.\n\n1. Section 4 states that, “without loss of generality”, the IRF construction is described for a one-dimensional target distribution. It is unclear for me how such a construction can be straightforwardly generalized to higher dimensions. Could the authors elaborate on this extension?\n\n2. Sections 5.1 and 5.2 report the results for the learning rate minimizing the TV distance and maximizing the ELBO, respectively. Although I acknowledge that this promotes a fair assessment of the proposed method, one is often interested in how the algorithm’s sensibility to hyperparametric changes. Authors should consider reporting the average and standard deviation of the selected metrics across runs. \n\n3. Contrarily to NUTS, RealNVP, and NSF, the IRF variants seem to require a warm-up phase to appropriately define a step size for the RWMH kernel. It is unclear, however, whether this initial step is incorporated into the computational cost analysis in Figure 4(d). Also, I could not understand whether this hyperparameter is the same for every experiment or was independently chosen for each target distribution. Authors should clarify the cost of IRF tuning. \n\n4. Synthetic experiments are based on 2-dimensional distributions that can be easily plotted. Could the authors show the resulting samples for each of the considered methods?", "was_revised": true, "final_justification": "I appreciate the authors' clear responses. I have increased my score to 5.", "rebuttal": "Thank you for your insightful comments and for acknowledging our work! Please see our point-to-point response below.\n\n> Strength: ... \n\nThank you again for your acknowledgement!\n\n> Questions:  Could the authors discuss the connection between the proposed methodology and traditional VI?\n\n\nWe appreciate this important question, as it clarifies the positioning of our work. The MixFlow framework still lies in the traditional VI paradigm—i.e., approximating a target distribution within a chosen variational family. However, MixFlows uniquely define a variational family that is asymptotically exact, meaning **every member within the family is guaranteed to converge to the target distribution as flow length increases**, independent of its hyperparameter tuning. Therefore, unlike typical black-box variational families such as normalizing flows, the convergence of MixFlows does not critically depend on meticulous hyperparameter tuning—though different hyperparameter choices can still influence the approximation quality when using a finite number of flows.\nMoreover, since MixFlows are parameterized through MCMC kernels, their associated hyperparameters (e.g., step size) are typically low-dimensional and interpretable. This facilitates intuitive and effective tuning, such as targeting specific MH acceptance rates, which we have found to be consistently reliable across various examples. Nonetheless, developing more sophisticated and robust tuning procedures remains an interesting direction for future investigation.\n\n> Section 4 states that, “without loss of generality”, the IRF construction is described for a one-dimensional target distribution. It is unclear for me how such a construction can be straightforwardly generalized to higher dimensions. Could the authors elaborate on this extension?\n\nThank you for highlighting this important point! We will clarify this extension and provide concrete examples (e.g., MALA, HMC, RWMH) for general dimensions in the final version of the paper (likely in the Appendix due to page limitations).\n\n\nTo briefly illustrate, we elaborate here on how the four-step IRF construction (lines 194–198) naturally generalizes to higher dimensions, using RWMH as a concrete example. Importantly, steps 3 and 4 remain unchanged when moving to higher dimensions; only steps 1 and 2 require minor adjustments:\n- Step 1 (line 194): For a general d-dimensional target, both $u_v , \\theta_v \\in \\mathbb{R}^d$, the step $ u_v \\gets (u_v + \\theta_v) mod 1$ is simply applied coordinate-wise.\n- Step 2 (line 195): In the specific case of RWMH (described in Example B.3, Appendix B.2), the distribution $ \\rho(\\mathrm{d} v |x) $ is a standard normal distribution $ N(0, I) $. Thus, we can implement step 2 in a coordinate-wise fashion as well. Specifically, for \n$$ i \\in 1, \\dots, d,  \\quad [u_v’]i \\gets F_i(v_i), \\quad [\\tilde v]i \\gets F_i^{-1}([u_{v}]i) $$\nWhere $ F_i $ and  $ F^{-1}_i $ represent, respectively, the CDF and inverse CDF of marginal (1-dimensional) standard normal distribution for the i-th coordinate.\n\n> Contrarily to NUTS, RealNVP, and NSF, the IRF variants seem to require a warm-up phase to appropriately define a step size for the RWMH kernel. It is unclear, however, whether this initial step is incorporated into the computational cost analysis in Figure 4(d). Also, I could not understand whether this hyperparameter is the same for every experiment or was independently chosen for each target distribution. Authors should clarify the cost of IRF tuning.\n\nThank you for pointing this out! We've revised the manuscript to clarify this issue explicitly. Figure 4(d) indeed reports the **total computational cost**, including both the tuning phase (tuning step size for the RWMH kernel targeting at a 0.8 MH acceptance rate) and the subsequent inference. Additionally, we now have a figure comparing the total wall-clock time for all methods, showing a similar phenomenon as Figure 4(d). We will include the wall clock time comparison in the camera-ready.\nTo further clarify, we did not reuse a single hyperparameter value across all runs. Instead, each independent run includes its own step-size tuning procedure.\n\n> Sections 5.1 and 5.2 report the results for the learning rate minimizing the TV distance and maximizing the ELBO, respectively. Although I acknowledge that this promotes a fair assessment of the proposed method, one is often interested in how the algorithm’s sensibility to hyperparametric changes. Authors should consider reporting the average and standard deviation of the selected metrics across runs.\n\nAs mentioned in the previous response, in all real‐data experiments (Section 5.2) we evaluated every MixFlow run—32 independent trials—using the same tuning procedure, and we compared these to the normalizing‐flow (NF) results obtained under their best‐tuned settings (thereby giving NFs a practical advantage in our comparisons).\nFor the synthetic benchmarks (Section 3.1), we plotted MixFlow performance across a range of kernel step sizes for the total variation metric (see Figure 2 in the main text and Figure 6 in Appendix C.1.2). Each curve represents the mean over 32 independent runs, and the shaded bands ($ \\pm 1 $ SD) depict run‐to‐run variability. We believe these results sufficiently demonstrate MixFlow’s robustness to hyperparameter variations.\n\n> Synthetic experiments are based on 2-dimensional distributions that can be easily plotted. Could the authors show the resulting samples for each of the considered methods?\n\nYes, we have obtained scatters of 1000 iid samples obtained from the IRF mixflows based on HMC kernels. The sample scatter matches the synthetic distribution almost perfectly. We will include this figure in the revised manuscript (perhaps included in the Apdx).\n\nAs for a reference, the original deterministic MixFlow work (based on uncorrected HMC kernels) shows almost perfect sample scatters (see Figure 1 in [26])."}, {"review_id": "lWdndkbzLa", "reviewer": "Reviewer_zTkn", "summary": "This paper aims at addressing convergence guarantees in variational inference (VI) models. \nThe quality of a variational approximation is fundamentally determined by the expressiveness of the variational model. Nonetheless, an expressive model is not necessarily guaranteed to converge due to the typical non-convexity of the optimization problem. \nConversely, MCMC methods are asymptotically exact but computationally intensive and not always differentiable or easily tunable. The paper tries to bridge this gap with this work. The authors introduce involutive MCMC kernels as invertible, measure-preserving iterated random functions (IRFs) acting as the flow maps of the (exact) variational flows. \n\nThe authors introduce some theorems to motivate their theorethical findings and validate their approach with extensive numerical experiments.", "questions": "- line 36: what do authors mean with \"optimal tuning\"?\n- line 52-53: how about stochastic normalizing flows or recent methods combining MCMC sampling with training of variational families?\n- eq(1) what is the difference between $x$ and $X$?\n- I am having a hard time understanding the rationale behind eq (1): is this essentially representing an ensemble of flows, each of which comes with a different number of transformations/diffeomorphisms? How are these diffeomorphisms constructed in the first place?\n- step 3 in Sec. 2.2: is the second term in the min a second order partial derivative or what does this notation represent?\n- Sec 2.3: eq (2): Isn't this the same as a NF with multiple couplings? Could the authors clearly comment on differences and analogies between traditional normalizing flows and IRFs-based flows?\n- can the author comment on what they mean by $\\mu-a.s.$ and $\\pi-a.e.$?\n- In equation at bottom of page 3: this seems to be just the same as the average of normalizing Flow with increasing number of coupling layers, where $f_\\theta_i$ represents a coupling layer with weights $\\theta_i$. Can the authors comment on this intuition?\n- Eq. 6 is reference on line 124 and later but never appears. \n- I am confused by first equation in Sec. 3.2: Shouldn't $f_\\theta$ act on $x$? How come it acts on the density $q_\\theta$ in that equation? (Please reference all equations with numbers so that it is easier to identify them).\n- First equation below line 151: I find the notation hard to parse and potentially not entirely consistent with the notation before. Shouldn't there be two sums in there? one running over t and one running over m?\n- Lines 157/158: Do the authors there mean to average over different M flows all with fixed depth T? Does T have to be the same for every flow?\n- Sec. 3.4 **relaxing ergodicity**: up to this point in the paper it is still very unclear to me where the MCMC involutive kernel contributes there. I believe it would be crucial for the authors to carefully revise the paper to make this concept emerging much clearer since the earlier sections of the manuscript. \n- line 179: with using flow density and not $q_{\\theta/T}$ here?\n- The confusion stemming from section 3.4 still applies to me in Sec. 4. Can the authors try to re-elaborate in different words on the crucial differences between their approach and standard NFs method? I believe a related work section and also a more focused section on differences/analogies between their deterministic flows and NFs would be crucial to enhance clarity. \n- Line 228-229: I am unsure if I get this statement. What makes an approach *reliably invertible*? To my understanding, only RWMH IRF seems to be reliably invertible, while all the other seems to quickly explode in errors after 100 steps. \n- Line 243-244: Why not start from a standard normal distribution?\n- Fig 4 is extremely hard to parse. Moreover, in some plots, it looks as if some approaches do not even appear in the plots. I wonder if this is a plotting issue or if the training was unsuccessful. If that is the case, I believe a more thorough discussion of these plots is needed. Perhaps it might be a good idea to carefully select a subset of these plots for the main text and thoroughly focus on those.", "limitations": "As far as I can tell, there's no section dedicated to the **Limitations** of the method, neither in the paper nor in the appendix.", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 1, "significance": 2, "originality": 2, "strengths_and_weaknesses": "# **Strengths**\n- The paper addresses an important and timely topic in the field of variational inference.\n- The experimental evaluation is comprehensive, covering a diverse range of benchmarks, including both synthetic and real-world datasets.\n\n# **Weaknesses**\n- The paper is difficult to follow at times and I believe is requires some efforts in enhancing clarity and accessibility.\n- The notation is occasionally inconsistent, and some notations or acronyms are not properly defined.\n- Certain concepts are assumed as known without sufficient explanation or background, which may hinder accessibility to a broader audience.\n- The writing lacks clarity in presenting the main contribution. In particular, the significance of constructing *exact* variational families via involutive MCMC kernels does not come through clearly in the exposition.\n- The visual presentation of experimental results is suboptimal. Figures are too small, as are the axis labels and legends. In addition, some plots lack sufficient detail in their captions (e.g., see Fig. 4).\n- The conclusions drawn from the experimental results appear broader than the initial scope outlined in the abstract. This mismatch may confuse readers—for example, Fig. 2 is more confusing than illuminating, as it does not clearly support the central claim that the proposed methods outperform standard normalizing flows.\n- Additional concerns are discussed further in the responses to the review questions.", "was_revised": true, "final_justification": "The authors' replies to my review and those from other reviewers have clarified many concern I had. \nNonetheless, I still believe the paper is hardly accessible to a broad audience and would greatly benefit from enhancing clarity in many parts of the manuscript (the authors themselves acknowledged that they will make parts of the manuscript more precise in the revised version). For this reason, I still believe the paper in its current form is not suited to the broad audience of NeurIPS. \n\nOn the other end, I was now able to fully appreciate and acknowledge the theoretical contribution of this work, hence the reason for raising my score. \n\nTo summarize, I will not fight to reject this paper should the other reviewers be strongly in favor for its acceptance, but at the same time I won't fight for it to be accepted either, in case the common agreement from the other reviewers is to reject the paper.", "rebuttal": "Thank you for reviewing our work! Please see our point-to-point response below.\n> **Weakness**: \nThe paper is difficult to follow at times... //\nThe notation is occasionally inconsistent... //\nCertain concepts are assumed as known without sufficient explanation or background...\n\nWe would appreciate more specificity from the reviewer regarding the issues with writing,  notations and undefined acronyms. We are happy to revise the manuscript to improve the clarity.\n\n> ... the significance of constructing exact variational families via involutive MCMC kernels does not come through clearly in the exposition\n\nThe main contributions of this work were stated as the 3 bullet points in the end of the introduction. And we point the review to the 3rd paragraph of the introduction on the significance of enabling the construction of MixFlows via general involutive MCMC kernels. As we explained explicitly in line 40-41, it is not clear how to design general invertible measure-preserving transformations for general continuous target distributions. This work addresses the limitation in full generality by constructing such transformations via involutive MCMC kernels.\n\n> The visual presentation of experimental results is suboptimal...\n\nThank you for your feedback. We have improved the visualizations of the figures by adjusting its size and better curation of the reported results, and have included more informative captions for all figures in the manuscript. \n\n> The conclusions drawn from the experimental results appear broader than the initial scope outlined in the abstract. This mismatch may confuse readers—for example, Fig. 2 is more confusing than illuminating, as it does not clearly support the central claim that the proposed methods outperform standard normalizing flows.\n\nWe respectfully disagree that our conclusions exceed the scope outlined in the abstract. First, Figure 2 does not involve any comparison with normalizing flows; as noted on lines 250–251, it simply contrasts the original Hamiltonian‑MixFlow—constructed with an uncorrected HMC kernel—with our corrected variant that includes the Metropolis–Hastings step. Second, our empirical results including comparisons against NFs (Figures 3 and 4) indeed support the abstract’s claim (lines 12–14) that our methods achieve competitive performance relative to standard normalizing flows.\n\nMore crucially, it is not our central claim that our proposed methods outperform standard normalizing flows. The key contribution of this work, as stated in the abstract (line 3-4) is the general construction of an asymptotically exact variational family via involutive MCMC kernels.\n\n> **Questions**: what do authors mean with \"optimal tuning\"?\n\nThis is an admittedly imprecise expression. We will revise this in camera-ready. We meant to state that the convergence of MixFlow does not rely on particular choice of its hyperparameters, reiterating its asymptotic exactness.\n\n> how about stochastic normalizing flows or recent methods combining MCMC sampling with training of variational families?\n\nThis is a good question and we will include the discussion in camera-ready!\n\nThese differentiabl AIS methods (DAIS; e.g., Wolf et al. 2016; Geffner & Domke 2021; Zhang et al. 2021; Thin et al. 2021b; Jankowiak & Phan 2021; Doucet et al. 2022; Geffner & Domke 2023) can be viewed as gradient-based tuning of AIS/SMC kernels or backward kernels. Unlike MixFlows, DAIS doesn’t produce an explicit density approximation, and its theory \n(Jankowiak & Phan 2021; Zhang et al. 2021; Zenn & Bamler 2024) typically operates under idealized assumptions—perfect transitions or optimal backward kernels—that do hold in practice or are only justifiable given an optimal tuning. These guarantees do not have the same generality of MixFlows’ asymptotic exactness.\n\n> eq(1) what is the difference between x and X?\n\nAs defined by $ x \\in \\mathcal{X} $, x here denotes a point of the state space $ \\mathcal{X} $. $ X \\sim \\bar q_T $ indicates that the random variable $X$ is distributed according to $ \\bar q_T  $.\n\n> can the author comment on what they mean by $ \\mu-a.s. $ and $ \\pi-a.e. $ ?\n\nHere, “ \\mu‑a.s.” stands for “\\mu‑almost surely” and “\\pi‑a.e.” for “\\pi‑almost everywhere.” Both are standard shorthand in probability theory, meaning the stated property holds except on a set of measure zero (with respect to $ \\mu $ or $ \\pi $, respectively).\n>  ...the rationale behind eq (1):...\n\nSection 2.1 offers a concise primer on deterministic MixFlows.  Eq (1) formalizes the fact that a deterministic MixFlow is a mixture over flows obtained by applying a single ergodic measure‑preserving (e.m.p.) map $ f $ repeatedly to the base distribution $ q_0 $. In other words, each mixture component corresponds to $ f^t $ for some iteration count $ t $.\n\nFor the concrete construction of the e.m.p. $ f $, as we explained in line 78–80, is not available in the literature for general continuous distribution and is indeed a major contribution of this work.\n\n> step 3 in Sec. 2.2: ...\n\n$ \\left|\\frac{\\partial g(x, v)}{\\partial(x, v)}\\right| $ denotes the determinant of the Jacobian matrix of the involution map $ g $. We will denote this explicitly in camera-ready.\n\n> Sec 2.3: eq (2): Isn't this the same as a NF with multiple couplings? Could the authors clearly comment on differences and analogies between traditional normalizing flows and IRFs-based flows?\n\nWe respectfully acknowledge that we’re unclear on what the reviewer means by “NFs with multiple couplings.” Equation (2) in Section 2.3 simply defines an IRF system, not a specific flow architecture.\n\nBroadly, the key differences between traditional normalizing flows (NFs) and our IRF‐based MixFlows are:\n(1). NFs apply a single invertible neural-network pushforward to a reference distribution; MixFlows form a mixture of multiple pushforwards.\n(2). NFs use target-agnostic flow maps, whereas MixFlows embed target information directly into each map, ensuring exact preservation.\n\nThese design choices result in the crucial difference between these two families of flows: MixFlows are asymptotically exact, meaning that **every member within the family is guaranteed to converge to the target distribution as flow length increases**, independent of its hyperparameter tuning. Traditional NFs do not offer this convergence guarantee.\n\n \n\n\n> shouldn't $f_\\theta$ act on $x$? How come it acts on the density...\n\nIn this work, given a transformation $ f $ and a distribution $ p $, we write $ f p$ (f acts on the distribution $ p $) as the pushforward distribution of $ p $ through $ f $. We defined this notation in line 72. \n\nThe equation $ \\overrightarrow{q_T}:=\\frac{1}{T} \\sum_{t=0}^{T-1} f_{\\theta_t} \\circ \\cdots \\circ f_{\\theta_1} q_0 $ denotes a mixture of $ T $ pushforward distributions $ \\\\{f_{\\theta_t} \\circ \\cdots \\circ f_{\\theta_1} q_0: t = 0, \\dots, T-1 \\\\} $. Same for the first displayed math of Section 3.2. The flow mapping $ f_{\\theta_i} $ will be the invertible measure-preserving IRF modified from involutive MCMC kernels (as we detailed in section 4).\n\n> In equation at bottom of page 3: this seems to be just the same as the average of normalizing Flow with increasing number of coupling layers, where f_\\theta_i represents a coupling layer with weights .\n\n$ f_{\\theta_i} $ is not a coupling layer. \nAs stated in line 101-103, $f_\\theta$ is a set of parameterized bijection that preserves $\\pi$. Thus, the parameterization of $f_\\theta$ has to be informed by each particular $\\pi$, and generic coupling flow layers wouldn’t suffice.\n\n> Eq. 6 is reference on line 124 and later but never appears.\n\nEq (6) appears in the Apdx C.1 (line 827 of the PDF in the submitted Supplementary material). We will clarify this in the main text for explicitness.\n\n> First equation below line 151: I find the notation hard to parse and potentially not entirely consistent with the notation before. Shouldn't there be two sums in there? one running over t and one running over m?\n\nThis is not an inconsistent notation nor a mistake. There is no sum running over $ t $. As mentioned in line 150-151, the ensemble IRF MixFlow is constructed by mixing over the final pushforwards of M IRF trajectories (sum over m), while all other IRF MixFlow variants mix along a single IRF trajectory (sum over t) . \n\n> Lines 157/158: Do the authors there mean to average over different M flows all with fixed depth T? Does T have to be the same for every flow?\n\nYes, we use a fixed depth T for all M flows. Though one can use different length T for each flow if there are good reasons.\n\n> with using flow density and not $q_\\theta/T$  here?\n\nThis expression aims to explain shared patterns of the 4 variants of MixFlows. To avoid the collision of the notation, each MixFlow is denoted differently. And we use “flow density ” to indicate the density expression for all the four flows to avoid redundancy in words. We will edit this in camera-ready to improve the preciseness.\n\n\n> To my understanding, only RWMH IRF seems to be reliably invertible, while all the other seems to quickly explode in errors after 100 steps.\n\nWe stated that HMC and MALA remain reliably invertible up to $ T \\approx 200 $ iterations. Indeed, Figure 1 shows that by $ T = 200 $ the inversion error for both methods stays between $ 10^{-10} - 10^{-2} $. We believe that this aligns with our claim.\n\n> Why not start from a standard normal distribution?\n\nWe adopted the same experimental setting from the original MixFlow work [26]. In general, it is recommended to choose a reference distribution $ q_0 $ to be close to the target distribution to avoid requiring excessive large flow length.  \n\n> Fig 4 is extremely hard to parse...\n\nYes, the excluded methods correspond to failure of the training. Again, we appreciate the reviews’ suggestion on improving the visualizations of the reported results. As responded above, we will work on this in revision."}]

Author rebuttal:

Meta-review: The authors give a general method for asymptotically exact variational inference using families of distributions obtained from normalizing flows, constructed from involutive MCMC kernels through representing them as iterated random function systems. They demonstrate improvement over previous methods. 

Their method is general and flexible, resolving practical difficulties with previous methods.

Reviewers had a lot of questions, which the authors answered satisfactorily. All authors increased their ratings and recommend acceptance.
I remind the authors to incorporate the suggestions into the final manuscript, especially clarifications to be accessible to a broad audience.
Decision: Accept (poster)