Embracing the chaos: analysis and diagnosis of numerical instability in variational flows
Authors: Zuheng Xu
Trevor Campbell
Venue: NeurIPS 2023 poster (2023)
Keywords: variational flow
numerical instability
shadowing property

Abstract: In this paper, we investigate the impact of numerical instability on the reliability of sampling, density evaluation, and evidence lower bound (ELBO) estimation in variational flows. We first empirically demonstrate that common flows can exhibit a catastrophic accumulation of error: the numerical flow map deviates significantly from the exact map---which affects sampling---and the numerical inverse flow map does not accurately recover the initial input---which affects density and ELBO computations. Surprisingly though, we find that results produced by flows are often accurate enough for applications despite the presence of serious numerical instability. In this work, we treat variational flows as chaotic dynamical systems, and leverage shadowing theory to elucidate this behavior via theoretical guarantees on the error of sampling, density evaluation, and ELBO estimation. Finally, we develop and empirically test a diagnostic procedure that can be used to validate results produced by numerically unstable flows in practice.

Embracing the chaos: analysis and diagnosis of numerical instability in variational flows
Zuheng Xu
Trevor Campbell
Department of Statistics University of British Columbia [zuheng. xu trevor] @stat wbc.
Abstract
In this paper, we investigate the impact of numerical instability on the reliability sampling, density evaluation, and evidence lower bound (ELBO) estimation variational flows. We first empirically demonstrate that common flows can exhibit catastrophic accumulation of error: the numerical flow map deviates significantly from the exact map map-which affects sampling and the numerical inverse flow map does not accurately recover the initial aaput-whiee affects density and ELBÒ computations. Surprisingly though, we find that results produced by flows are often accurate enough for applications despite the presence of serious numerical instability. In this work, we treat variational flows dynamical systems, and leverage shadowing theory to elucidate this behavior via theoretical guarantees on the error of sampling, density evaluation, and ELBO estimation. Finally, we develop and empirically test a diagnostic procedure that can be used to validate results produced by numerically unstable flows in practice.
Introduction
Variational families of probability distributions play prominent role in generative modelling and probabilistic inference. A standard construction of a variational family involves passing a simple reference distribution-such as standard normal normal-through sequence of parametrized invertible transformations, i.e., flow [1- Modern flow-based variational families have the representational capacity to adapt to highly complex target distributions while still providing computationally tractable i..d. draws and density evaluations, which in turn enables convenient training by minimizing the Kullback-Leibler (KI.) divergence
[ via stochastic optimization methods 11111
In practice, creating flexible variational flow often requires many flow layers, as the expressive power of variational flow generally increases with the number of transformations (5] 12 13]. Composing
many transformations in this manner, however, can create numerical instability (5. 14 i.c., the
tendency for numerical errors from imprecise floating-pons -pint computations to be quickly magnified along the flow. This accumulation of error can result in low-quality sample generation [4], numerical non-invertibility [ [1),, and unstable training [16]. But in this work, we demonstrate that this is not always the case with a counterintuitive empirical example: a flow that exhibits severe numerical instability- with error that grows exponentially in the length of the flow but which surprisingly still returns accurate density values, i.i.d. draws, and evidence lower bound (ELBO) estimates.
Motivated by this example, the goal of this work is twofold: (1) to provide a theoretical understanding of how numerical instability in flows relates to error in downstream results; and (2) to provide diagnostic procedure such that users of normalizing flows can check whether their particular results are reliable in practice. We develop a theoretical framework that investigates the influence numerical instability using shadowing theory from dynamical systems [17-20]. Intuitively, shadowing theory asserts that although a numerical flow trajectory may diverge quickly from its exact counterpart,
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

there often exists another exact trajectory nearby that remains close to the numerical trajectory (the shadowing property). We provide theoretical enror bounds for three core tasks given a flow exhibiting the shadowing property- -ample generation, density evaluation, and ELBO estimation that show that the error grows much more slowly than the error of the trajectory itself. Our results pertain to both normalizing flows 2 3 and mixed variational flows [ We develop diagnostic procedure for practical verification of the shadowing property of a given flow. Finally, we validate our theory and diagnostic procedure MixFlow [5] on both synthetic and real data examples, and Hamiltonian variational flow [6 on synthetic examples.
Related work. Chang et al. [21] connected the numerical invertibility of deep residual networks (ResNets) and ODE stability analysis by interpreting ResNets as discretized differential equations, but did not provide quantitative analysis of how stability affects the error of sample generation or density evaluation. This approach also does not apply to general flow transformations such as coupling layers. Behrmann et al. [14] analyzed the numerical inversion error of generic normalizing flows via bi-Lipschitz continuity for each flow layer, resulting in error bounds that grow exponentially with flow length. These bounds do not reflect empirical results in which the error in downstream statistical procedures (sampling, density evaluation, and ELBO estimation) tends to grow much more slowly. Beyond the variational literature, Bahsoun et al. 22 investigated statistical properties of numerical trajectories by modeling numerical round-off errors as small random perturbations, thus viewing the numerical trajectory as a sample path of Markov chain. However, their analysis was constrained by the requirement that the exact trajectory follows a time-homogeneous, strongly mixing measure- presurving dynamical system, condition that does generally hold for variational flows. Tupper [23] established a relationship between shadowing and the weak distance between inexact and exact orbits in dynamical systems. Unlike our work, their results are limited quantifying the quality of pushforward samples and do not extend to the evaluation of densities and ELBO estimation. In contrast, our error analysis provides a more comprehensive understanding of the numerical error of all three downstream statistical procedures.
2 Background: variational flows
variational family is parametrized set of probability distributions {ąx E A} on some space X' that each enable tractable i.i.d. sampling and density evaluation. The ability to obtain draws and compute density values is crucial for fitting the family data via maximum likelihood in generative modelling. as well as maximizing the ELBO in variational inference.
/ qx (z) log
p(x) QA (z)
ELBO (4x Ip)
dx,
where p(x) is an unnormalized density corresponding to target probability distribution x. In this work focus on X' C Rd endowed with its Borel o-aleebra and Fuclidean norm || |], and parameter set A C RD. We assume all distributions have densities with respect to a common base measure on X, and will use the same symbol to denote distribution and its density. Finally, we will suppress the parameter subscript A, since we consider fixed member of a variational family throughout.
Normalizing flows. One approach to building variational distribution is push reference distribution, qu. through a measurable bijection on X. To ensure the function is flexible yet still enables tractable computation, it is often constructed via the composition of simpler measurable, invertible "layers" Fi...., FN, referred together as a (normalizing) flow 1 2 4. To generate a draw Y q. draw X qo and evaluate Y FN Fx-1 .... F (X). When each F, is diffeomorphism the density of q
(x) 11Th N Jn (F. .... F (x))
4o F.-T
Vr X
((z)
Jn(2) |det (zF,(z)). (1)
An unbiased ELBO estimate can be obtained using draw from go via
X ~ B(X) of go (X)
log
(2)
n=
AA differentiable map with differentiable inverse; we assume all flow maps in this work are diffeomorphisms.

---

MixFlow Fwd Orbit (N 30)
MixFlow numerical error
MixFlow numerical error
30
exact numerical
10
10
20
Error
Error
10
10 -10
10
-10
Rod
Basi
20
30
20
60 100
-20
-10
10
20
#transformations
#transformations
Banana: num. and exact orbits
Banana: flow error
(c) Lin. Reg.: flow error
MixFlow Fwd Orbit (N = 50)
MixFlow numerical error
MixFlow numerical error
exact numerical
10
Two
10
Error
10
10-
-2
-4
40
100
20
100
2
2
4
# #transformations
#transformations
Cross: num. and exact orbits
(e) Cross: flow error
(f) Log. Reg.: flow error
Figure I: MixFlow forward (fwd) and backward (bwd) orbit errors on the Banana and Cross distribu- tions, and
real data examples. The first column visualizes the exact and numerical orbits with the same starting point on the synthetic targets. For a better visualization, we display every 2nd and 4th states of presented orbits (instead of the complete orbits) in Figs. and Id respectively. Figs.
lb
lc
le
and
If show the median and upper/llwer quartile forward error | p and
backward error
|B x B* comparing & transformations of the forward exact/approximate maps
F, F or backward exact/aprroximate maps B, B. Statistics are plotted over 100 initialization draws from the reference distribution Qo- For the exact maps we use 2048-bit BigFloat representation, and for the numerica approximations we use 64-bit Float representation. The "exactness" of BigFloat representation is justified in Fig. 15 in Appendix D.3
Common choices of Fn are RealNVP 24 neural spline 25 Sylvester (26 Hamiltonian-basee 7,7 and planar and radial 2 flows; see Papamakarios et al. for a comprehensive overview.
MixFlows. Instead of using only the final pushforward distribution as the variational distribution Xu et al. [5] proposed averaging over all the pushforwards along flow trajectory with identical flow layers, i.c., Fr. F for some fixed F. When is ergodic and measure- peasure-ppreerving for some target distribution π, Xu et al. [5] show that averaged flows guarantee total variation convergence to æ in the limit of increasing flow length. To generate draw ~ 4. we first draw X qo and flow length K ~ Unif 1.. N), and then evaluate Y = FK (X) = F o F o ..... F(X) via K iterations of the map The density formula is given by the mixture of intermediate flow densities,
......,..
Va E X,
(xx)
Σ
N
No
(x-
J(x) = |det V. F(x))
(3)
N N+1
J(F-1))
n=0
An unbiased ELBO estimate can be obtained using a draw from go via
1
X ~ 4o E(X)
log q (F
(4)
N
Instability in variational flows and layer-wise error analysis
In practice, normalizing flows often involve tens of flow layers, with each layer coupled with deep neural network blocks 2 24 or discretized ODE simulations (6] 7 past work on MixFlows requires thousands of flow transformations
[5] These flows are typically chaotic and sensitive to small perturbations, such as floating-poii point representation errors. As result, the final output of a numerically computed flow can significantly deviate from its exact counterpart, which may cause downstream To understand the accumulation of error,
issues during training or evaluation of the flow 12 14 27

---

MixFlow sample scatters
MixFlow sample scatters
30
20
2
10
-2
-10
-4
-20
-10
10
20
-2
2
Banana sample seatter
(b) Cross sample scatter
(c) Banana log-densities
(d) Cross log-densities
Figure 2: Figs. 2a and 2b respectively show sample scatters produced by the naïve application of MixFlow formulae targeting at the banana and cross distributions, without accounting for numerical crror. Figs, 2c and 2d display comparisons of log-densities on both synthetic examples. The true log- target the left, the exact MixFlow evaluation (computed via 2048-bit BigFloat representation) in the middle, and the numerical MixFlow evaluation is on the right.
Behrmann al. 14 assume that each flow layer Fn and its inverse are Lipschitz and have bounded error, sup, Fn (z)
(x) S 8, and analyze the error each layer individually. This layer-wise analysis tends to yield error bounds that grow exponentially in the flow length N (see Appendix A.T):
N-1
| o ... F((2) FN o. o
<
Σ II Lip(Fs).
(5)
1= j=n+2
Given a constant flow map F, = with Lipschitz constant the bound is )) which suggests that the етrог accumulates exponentially in the length of the flow. In Fig- [] we test this hypothesis empirically using MixFlows on two synthetic targets (the banana and cross distribution) and two real data examples (Bayesian linear regression and logistic problems) taken from past work on MixFlows [ See Appendix D for the details of this test. Figs. Tb] Tc Te and Tf confirms that the exponential growth in the error bound Eq. (5) is reasonable; the error does indeed grow exponentially quickly in practice. And Fig. 16 (b) in Appendix D.3) further demonstrate that after fewer than 100 transformations both flows have error on the same order of magnitude as the scale of the target distribution. Naïvely, this implies that sampling, density evaluation, and ELBO estimation may be corrupted badly by numerical error.
But counterintuitively, Fig. 2 we find that simply ignoring the issue of numerical error and using the exact formulae yield reasonable-lookin of course only qualitative in nature; we provide corresponding quantitative results later in Section 6 But Fig. apppars violate the principle of "garbage in, garbage out:" the buildup of significant numerical error in the sample trajectories themselves does not seem to have strong effect on the quality of downstream sampling and density evaluation. The remainder of this paper is dedicated to resolving this countcrintuitive behavior using shadowing theory (28).
density evaluations and sample scatters. These results
Global error analysis of variational flows via shadowing
4.1 The shadowing property
We analyze variational flows as finite discrete-ting ticrettettmme dynamical systems (see, e.g., (28). In this work. we consider dynamical system be sequence of diffcomorphisms on (X, |l ID. We define the forward dynamics (F,) N= - to be the flow layer sequence, and the backward dynamics (B,) to be comprised of the inverted flow layers B, = F (n- in reverse order. An orbit of dynamical system starting e X' is the sequence of states produced by the sequence of maps when initialized

---

$4
14
$2
$3
(2(
$
50
33
00
00
21
22
34 F4(13)
3
Figure 3: A visualization of a pseudo-orbit and shadowing orbit. Solid arrows and filled dots indicate exact orbits, while dashed arrows and open dots indicate pseudo-orbits (e.g., via numerical compu- tations). Red indicates the exact orbit (zu..... I4) that one intends to compute. Blue indicates the numerically computed d-pseudo-oriit (fg. ................................ &4). Grey indicates the corresponding c-shadowing orbit : 84). At the top right of the figure, || Fa($a)ll < 8 illustrates the 6 numerical enror each step of the pseudo-orbit. The grey dashed circles demonstrate the e-shadowing window.
at Therefore, the forward and backward orbits initialized at X' are
Forward Orbit:
x = 3o
f x1
Exy ON
N-N Bx=F
Pa=FN- B Fel
2 :-1 20
Backward Orbit.
Given numerical implementations Fn. so Fn and Bn Bn with tolerance 0, i.e.
Va X, [F (n( S8, [[ Bm, (x)] <
(6)
define the forward and backward pseudo-dynanine be (Fr) and (B,, W respectively. along with their forward and backward pseudo-orbits initialized at x E Ex:
AIsFi
PugPA 22
FNNFFN AN
Forward Pseudo-Orbit
f1
B --N N
F B
R-1 in = x Backward Pseudo-Orbit.
For notational brevity, we subscripts on elements of X throughout to denote forward/backward orbit states. For example, given a random element Z E X drawn from some distribution, Z₁. is the ttth state generated by the forward dynamics (F,, N initialized Z, and Z- is the 1.th state generated by the backward dynamics (B,), initialized at Z. Hat accents denote the same for the pseudo-dynamics: for example. 2n is the tth state generated by (F.) when initialized at Z.
The forward/backward orbits satisfy J (*() and (1++) (xx=)) exactly at each step. On the other hand, the forward/backward pseudo-orbits incur a small amount of error,
| (+x))l 5 and || (*
which can be magnified quickly along the orbit. However, there often exists another exact orbit starting at some other point E that remains in a close neighbourhood of the numerical orbit. This property, illustrated in Fig- 3 is referred to as the shadowing property [19 31.
Definition 4.1 (Shadowing property (8)). The forward dynamics (Fw) has the (c, 8)-shdowinn property for all
X and all (F,) satisfying Eq. (6), there exists E X such that
Vk 1,1,..,,
|| Exll <.
An analogous definition holds for the backward dynamics (B.) there is shadowing orbit a that is nearby the pseudo-orbit f=---add for the joint forward and backward dynamics, where there is shadowing orbit nearby both the backward and forward pseudo-orbits simultaneously. The key idea in this paper is that, intuitively, if the numerically computed pseudo-orbit is close to
5

---

some exact orbit (the shadowing orbit), statistical computations based on the pseudo-orbie e.g.. sampling, density evaluation, and ELBO cioimation-sholle be close those obtained via that exact orbit. We will defer the examination of when (and to what extent) shadowing holds to Section 5 in this section, we will use it as an assumption when analyzing statistical computations with numerically implemented variational flows.
4.2 Error analysis of normalizing flows via shadowing
Sampling, Our first goal is to relate the marginal distributions of X and XN for X i... to quantity the error sampling due to numerical approximation. Assume the normalizing flow has the (e, -eshaddwwnin property, and let &o be the marginal distribution of the shadowing orbit start point. We suspect that &o & qu in some sense; for example, distance Drr.(S (So. go) is at most due to shadowing. And Eu is indeed an implicit function of qu: it a fixed point of a twice differentiable function involving the whole orbit starting at qu [17 Page. 176]. But the distribution Ço is generally hard to describe more completely, and thus it is common to impose additional assumptions. For example, past work shows that the Sevy-Prokhorov metric DLP .NNNNN is bounded by e, under the assumption that Eo = qo [23]. We provide a more general result (Proposition 4.2) without distributional assumptions Dyp(: (:))). :), as its analysis is simpler and both metrize weak distance.
know that the bounded Lipschitz
Eo. We control Br. (.-)) rather than
Proposition 4.2. Suppose the forward dynamics has the (c, 8) -hhadoowgg property, and X~ Then
Ef(Xx) - BF((NN) = & + 2UDTV (So, 4o).
sup SUALI(())
In particular, with é = U = 1, we obtain that D XN) + (Evv qu).
Recall the layerwise error bound from Eq. 5 which suggests that the difference in orbit and pseudo-orbit grows exponentially in N-and compare to Proposition 4.2] which asserts that the error controlled by the shadowing window size c. This window size may depend on N, but we find in Section it is usually not much larger than 8 in practice, which itself is typically near the precision of the relevant numerical representation. We will show how to estimate e later in Section
Density evaluation. The exact density g(z) follows Eq. L.. while the approximation ((x) is the same except that we use the backward pseudo-dynamics. For E X, a differentiable function g X + R+, define the local Lipschitz constant for the logarithm of around x as:
Lgg (w)
sup
V (u)g(
11
Theorem shows that, given the shadowing property, the numerical error is controlled by the shadowing window size e and the sum of the Lipschitz constants along the pseudo-orbit. The constant Lg.e occurs because we are essentially evaluating ((s) rather than g(z), where E X' is the backward shadowing orbit initialization. The remaining constants occur because of the approximation of the shadowing orbit with the nearby, camerically-comppeee pseudo-orbit in the density formula.
Theorem 4.3. Suppose the backward dynamics has the (c, 6)-shdowiing property. Then
Vr E X, log ((x) log ((x) Se. (maac(e) + ( + Σ LJ To (8-))) )
The exact ELBO estimation function is given in Eq. @2); the numerical ELBO estimation function B(x)) is the same except that we use the forward pseulo-dynamics. The quantity E(X), X ~ qo is an unbiased estimate of the exact ELBO, while É(X) biased by numerical error; Theorem 4.4 quantifies this error. Note that for simplicity we assume that the initial state distributions qo. Eu described earlier are identical. It is possible to obtain a bound including Drv(2o {o) term, but this would require the assumption that log p(x)//(//// (q(x) is uniformly bounded X.
ELBO estimation.
Theorem 4.4. Suppose the forward dynamics has the (c, 60shhwwowing property; and Eo = go- Then
N
ELBO (gllp) E[β(X)I S (xe(XN)) + Lon (X) Σ
Ls (Xa)
for X ~ go.
n=1

---

4.3 Error analysis of MixFlows via shadowing
The error analysis for MixFlows for finite length E N parallels that of normalizing flows, except that F, F. B,, = B, and J, = J for = ......N. However, when F and B are ergodic and preserrii for the target æ, we provide asymptotically simpler bounds the large N limit that do not depend on the difficult-O--nallzze details of the Lipschitz constants along pseudo-orbit. These results show that the error of sampling tends be constant flow length, while the error of density evaluation and ELBO estimation grows at most linearly. We say the forward dynamics has the infinite (e, d) -shadowinng property if (Fr) has the (e, -shadowwi property [32 Analogous definitions hold for both the backward and joint forward/backward dynamics.
Sampling. Similar to Proposition 4.2 for normalizing flows, Proposition 4.5 bounds the error in exact draws and approximate draws Y from the MixFlow. The result demonstrates that error does directly depend the flow length N, but rather on the shadowing window In addition, in the setting where the flow map F is æ-ergodic and peeasure-preserving we can (asymptotically) remove the total variation term between the initial shadowing distribution ço and the reference distribution go-
Note that in the asymptotic bound, the distributions of Y and Y are functions of N.
Proposition 4.5. Suppose the forward dynamics has the (c, () -hhadowing property, and X q0- Let Y XK. γ Xx for K Unif N.. Then
Ef(y) - Ef(r) $ sle+ 2UDTY (Eo. 90)
sup AASISULLPp(SS
In particular, é = 1, we obtain that DRI. (Y, S 2Dvv(Eo, qo). If additionnlly the forward dynamics has the infinite (e, 8) -sshdwwwigg property, -eerodic and meeasre-preservvnng and Eo. qu < then
lim sup D BL N N
(1,)) SE.
direct corollary of the second result in Proposition 4.5 is that for W ~ π., limn-+ox BEL (W, S e,
which results from the fact that DBr (W, (W, Y) - 0 15 Theorem 4.1]. In other words, the bounded Lipschitz distance of the approximated MixFlow and the target æ is asymptotically controlled by the shadowing window c. This improves upon earlier theory governing approximate MixFlows, for which the best guarantee available had total variation error growing as O(N) (5] Theorem 4.3].
Density evaluation. The exact density follows Eq. (3); the numerical approximation is the same except that we use the backward pseudo-dynamics. We obtain a similar finite-N result as in Theo- rem 4.3 Further, we show that given infinite shadowing-wherr the shadowing window size is independent of flow length N-the numerical approximation error asymptotically grows at most linearly in N, in proportion to e.
Theorem 4.6. Suppose the (e, Össhadowwingg property holds for the backward dynamics. Then
N
VI X, [log((x) - log g(x)l s. (Ln,((x)+ + max La (ā_ +
Σ LJ. (i)c((n)) n=
OSAEN
If additionally the backward dynamics has the infinite (e, 8--shhdoownng property, F is --ergodic and measure- prresuri...ee.vingg. Lg,e(x) o(N) as N sc. and Eu < n, then for glmalmost every x E X,
lim sup N N
log ((x) log (x=) S [Lgg, Lx)(()),, X~m.
The exact ELBO formula for the MixFlow is given in Eq. Note that here we do not simply substitute the forward/backward pseudo-orbits as needed; the naive approximation of the terms q( (F)x) would involve n applications of f² followed by N applications of B, which do not exactly invert one another. Instead, we analyze the method proposed in [5), which involves simulating a single forward pseudo-orbit .,,..., &n and backward pseudo-orbit ........ &-N starting from E X, and then caching these and using them needed in the exact formula.
ELBO estimation.
7

---

Theorem 4.7. Suppose the joint forward and backward dynamics has the (c, (C.))saahdowing property, and Eo q0- Then for X q0-
N Σ n=0
ELBO (gllp) (E [(x) E.
max oo (Xn) + Σ (xn( (xne(Xn)+
N
-NSnSN
R=-N
If additionally the joint forward and backward dynamics has the infinite (e, 5) 8)-shddowing property, F is x-ergodic and meeasre-ppresente and for some S m1 < oo and 1/m1 + 1/m2 1 Lpp.s Leguess LIe E [= (T) and dgo dı E (T), then
lim sup N
ELBO (gllp)
[(X)
S 2E 2E (e + (x,(()),, X~m.
N
5 Computation of the shadowing window size
We have so far assumed the (c, --sshddonnng property throughout. To make use of our results in practice, is crucial understand the shadowing window size e. (.. conditions for finite dynamical system to have the shadowing property, and characterizes the size of the shadowing window c. Note that throughout this section we focus on the forward dynamies; the backward and joint dynamics can be treated identically. Let || denote the spectral norm of matrix.
sufficient
Theorem 5.1 (Finite shadowing theorem). Suppose the dynamics (Fa) are C2 diffeomorphisms on X, and the sneudo-dynamice (Fa)m= satisfy Eq. 16). For a given E X, define the operator A: XN++1 + X by
(Au)s up+1 FFk+1 for u = (2o, U1,.... uN) E XN, k
-
Let
{ Suu saXs 11 Fn+ (fn, + w)ll 1,,,,...,NN -1] where X Amin (AA7)-1/2
M max
1f2MX²8 < then the pseudo-orbit starting at x is shadowed with 2Ad.
Proof. This result follows directly by following the proof of [18, Theorem 11.3] with nonconstant flow maps Fn. and then using the right-inverse norm from 33 Corollary 4.2].
In order to apply Theorem 5.1 we need to (1) estimate e.g., by examining numerical error of one step of the map in practice; (2) compute Amin ( (AA"); and (3) estimate M, e.g., by bounding the third derivative of F,, While estimating M and ō are problem-specifie, one can employ standard procedures for computing Amin (AA"). The matrix AAA has -oock-triiaagonll form.
D +1
-DT D2D +1
-D2
-DJ
AAT
-DN-1 DN-IDE +1 -D
DT DN +1
where D2 = FF E Rddd Vk E [N]. Notice that AAT is a symmetric positive definite sparse matrix with bandwidth and so has O(Nd²) entries. The inherent structured sparsity can be leveraged to design efficient eigenvalue computation methods, e.g., the inverse power method or tridiagonalization via Lanczos iterations (35) followed by divide-and-conques algorithms However, in our experiments, directly calling the e ggiin function provided in Julia suffices; as illustrated in Figs. [4] and [9 the shadowing window computation requires only a few seconds for low dimensional synthetic examples with hundreds flow layers, and several minutes for real data examples. Hence, we didn't pursue specialized methods, leaving that for future work. also noteworthy that practical computations of Dr can introduce floating-poin errors, influencing the accuracy of X. To address this, one might consider adjusting the shadowing window size. We explain how to manage this in Appendix B.I and why such numerical discrepancies minimally impact results.
34
36

---

sampling
sampling
sampling errer
sampling
000000 1250000
10"
Rei or
Rell er
10
10
10"
10 8
& 10
10
10
10
200 a00
dtarstormations
atransformations
#transformations
#transtormations
(a) Banana
(b) Cross
(c) Lin. Reg.
(d) Log. Reg.
Figure 4: MixFlow relative sample average computation error on three test functions: D= E (sin()) + 11; and Σ (s))moid The lines indicate the median, and er- ror regions indicate 25th to t5th percentile from 20 runs.
MixFlow og de-density errer
MixFlow log density error
MixFlow log donsity error
MixFlow log density error
10"
16-3
er.
10
10
10
10
10
10
10 10-
10
10
10
10"
10
4,00400g000000000000
z00
600
g00
0,00,00000000000000000
600y00yd040000000000
ttrnnsformations
stransformations
@transfomatiors
#transformations
(a) Banana
(b) Cross
(c) Lin. Reg.
(d) Log. Reg.
Figure 5: MixFlow relative log-density evaluation error. The lines indicate the median, and error
regions indicate 25th to t5th percentile from 100 evaluations on different positions.
MixFlow
MixFlow
MixFlow
MixFlow ELBO
EU80
550855
EL80
10 10---a
-2209.0
-2209.5
22100.
-250.5
-250.4
-2210.3
-250.6
-2211.0
-250.7
Turwic
300
:00 y00
200
600
ttrandformations
Strandformationn
stranstormations
atranstormations
(a) Banana
(b) Cross
(c) Lin. Reg.
(d) Log. Reg.
Figure 6: MixFlow exact and numerical ELBO estimation over increasing flow length. The lines indicate the averaged ELBO estimates over 200 independent forward orbits. The Monte Carlo error for the given estimates sufficiently small so we omit the error bar.
4.0x10"
40
6.0x30
4.0x10
8.0x
6.80x10
1000 30
3.0x10"
33
2.0x10"
2.0+10
2.00-10
1.0110"
300
600
600
s00
etransfomations
stronsfomatiios:
#transfemations
(a) Banana
(b) Cross
(c) Linear regression
(d) Logistic regression
Figure 7: MixFlow shadowing window size e over increasing flow length. The lines indicate the median, and eиror regions indicate 25th 75 percentile from 10 runs.
Finally, since e depends on Amin (AA7)- which itself depends on N, the shadowing window size may potentially increase with N in practice. Understanding this dependence accurately for general dynamical system is challenging; there are examples where it remains constant, for instance (see Examples [B. and B.2 in Appendix [B.2]. but may grow with N. Our empirical results in next section show that on representative inferential examples with over 500 flow layers, scales roughly linearly with N and is of a similar order of magnitude as the floating point representation error. For further discussion on this topic, we refer readers to Appendix B.2
Experiments
In this section, we verify our error bounds and diagnostic procedure of MixFlow on the banana, cross, and 2 real data targets-dl Bayesian linear regression and logistic regression posterior; detailed

---

model and dataset descriptions can be found in Appendix D.I investigation for Hamiltonian flow [6] on the same synthetic targets in Appendix C.2
We also provide a similar empirical
We begin by assessing the error of trajectory-averagedd estimates based on approximate draws. Fig. displays the relative numerical estimation error compared to the exact estimation based on the same initial draw. Although substantial numerical error was observed in orbit computation (see Fig. [6 in Appendix D.3). the relative sample estimate error was around 1% for all four examples, suggesting that the statistical properties of the forward orbits closely resemble those of the exact orbit.
We then focus the density evaluation. For synthetic examples, assessed densities at 100 evenly distributed points within the target region (Figs.[2)) 2aland (2b). For real data, we evaluated densities at the locations of 100 samples from MCMC method None-U sampler (NUTS); detailed settings for NUTS is described in Appendix D.2 It is evident from the relative error shown in Fig. thaat the numerical density closely matches the exact density evaluation, with the relative numerical error ranging between 0.1% and 1 1,, which is quite small. The absolute error can be found in Fig. 17 in Appendix D.3] showing that the density evaluation error does not grow substantially as the orbit computation error (Fig. 16 in Appendix D.3). which aligns with the bound in Theorem 4.6
Fig. 6 further demonstrates the numerical error for the ELBO estimations (Eq. In each example, both the averaged exact ELBO estimates and numerical ELBO estimates are plotted against an increasing flow length. Each ELBO curve is averaged over 200 independent forward orbits. Across all four examples, the numerical ELBO curve remain closely aligned with the exact ELBO curve, indicating that the numerical error is small in comparison to the scale of the ELBO value. Moreover, the error does not grow with an increasing flow length, and even presents a decreasing trend as N increases in the two synthetic examples and the Logistic regression example. This aligns well with the error bounds provided in Theorem 4.7
Finally, Fig. presents the size of the shadowing window as the flow length N increases. As noted earlier, e depends on the flow map approximation error and potentially N. We evaluated the size of 5 by calculating the approximation error of single F or B for 1000 i..d. draws from the reference distribution. Boxplots of 8 for all examples can be found in Fig. 18 of Appendix D.3 These results show that in the four examples, 5 is close to the floating point representation error (approximately 10-14). Thus, we fixed at 10-14 when computing e. As shown in Fig. 7 e for both the forward and backward orbits grows roughly linearly with the flow length. This growth is significantly less drastic than the exponential growth of the orbit computation error. And crucially. the shadowing window than 10-10 10 for MixFlow of length over 1000, which justifies the validity of the downstream statistical procedures. In contrast, the orbit computation errors in the two synthetic examples rapidly reach a magnitude similar to the scale of the target distributions, as shown in Fig. D
size is very
7 Conclusion
This work delves into the numerical stability of variational flows, drawing insights from shadowing theory and introducing diagnostic tool to assess the shadowing property. Experiments corroborate our theory and demonstrate the effectiveness of the diagnostic procedure. However, the scope of our current analysis is limited to downstream tasks post-training of discrete-time flows. Understanding the impact of numerical instability during the training phase or probing into recent architectures like continuous normalizing flows or neural ODEs (21 37 40), remains an avenue for further exploration.
Additionally, while our theory is centered around error bounds that are proportional to the shadowing window size we have recognized that can grow with N. Further study its theoretical growth rate is needed. Finally, in our experiments, we employed a basic method to calculate the minimum eigenvalue of AAT. Given the sparsity of this matrix, deeper exploration into more efficient computational procedure is merited.
References
[1] Esteban Tabak and Cristina Turner. family of non-parametric density estimation algorithms. Communi- cations on Pure and Applied Mathematics,
2013.
[2] Danilo Rezende and Shakir Mohamed. Variational inference with normalizing flows. In International Conference on Machine Learning, 2015.
10

---

[3] 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.
[4] Ivan Kobyzev, Simon Prince, and Marcus Brubaker. Normalizing flows: introduction and review of current methods. IEEE Transactions on Pattern Analysis and Machine Intelligence, 43(1 43(1 -)33943333397 2021.
[5] Zuheng Xu, Naitong Chen, and Trevor Campbell. MixFlows: principled variational inference via mixed flows. arXiv:2205.07477 077475, 2022.
[6] Naitong Chen, Zubeng Xu, and Trevor Campbell. Bayesian inference via sparse Hamiltonian flows. In Advances in Neural Information Processing Systems, 2022.
[7] Anthony Caterini, Arnaud Doucet, and Dino Sejdinovic. Hamiltonian variational auto-encoder. Advances in Neural Information Processing Systems, 2018.
[8] Solomon Kullback and Richard Leibler. On information and sufficiency. The Annals of Mathematical Statistics, 22(1):79- 86, 1951.
[9] Matthew Hoffmann, David Blei, Chong Wang, and John Paisley. Stochastic variational inference. Journal of Machine Learning Research, 14:1303-134,, 2013.
[10] Alp Kucukelbir, Dustin Tran, Rajesh Ranganath, Andrew Gelman, and David Blei. Automatic Differentia- tion Variational Inference. Journal of Machine Learning Research, 18(14), 2017.
[11] Rajesh Ranganath, Sean Gerrish, and David Blei. Black box variational inference. International Conference on Artificial Intelligence and Statistics, 2014.
[12] Jonas Köhler, Andreas Krämer, and Frank Noe. Smooth normalizing flows. Advances in Neural Information Processing Systems, 2021.
[13] Zhifeng Kong and Kamalika Chaudhuri. The expressive power of class of normalizing flow models. In International Conference on Artificial Intelligence and Statistics, 2020.
[14] Jens Behrmann, Paul Vicol, Kuan-Chieh Wang, Roger Grosse, and Jörn-Henrik Jacobsen. Understanding and mitigating exploding inverses invertible neural networks. International Conference on Artificial Intelligence and Statistics, 2021.
[15] Shifeng Zhang, Chen Zhang, Ning Kang, and Zhenguo Li. ivpf: Numerical invertible volume preserving flow for efficient lossless compression. In The IEEE/CVF Conference on Computer Vision and Pateen Recognition, 2021.
[16] Aidan N Gomez, Mengye Ren, Raquel Urtasun, and Roger Grosse. The reversible residual network: Backpropagation without storing activations. In Advances in Neural Information Processing Systems, 2017.
[17] Brian Coomes, Hüscyin Koçak, and Kenneth Palmer. Shadowing discrete dynamical systems, In Six lectures on dynamical systems, pages 163-211. World Scientifie, 1996.
[18] Kenneth Palmer. Shadowing in dynamical systems: theory and applications, volume 501. Springer Science & Business Media, 2000.
[19] Sergei Yu Pilyugin. Shadowing in dynamical systems. Springer, 2006.
[20] Lucas Backes and Davor Dragičević. Shadowing for nonautonomous dynamies. Advanced Nonlinear Studies, 19(2):4254433 2019.
[21] Bo Chang, Lili Meng, Eldad Haber, Lars Ruthotto, David Begert, and Elliot Holtham. Reversible architectures for arbitrarily deep residual neural networks. In Proceedings of the AAAI Conference on Artificial Intelligence. 2018.
[22] Wacl Bahsoun, Huyi Hu, and Sandro Vaienti. Pscudo-orbits, stationary measures and metastability. dynamical systems, 29(3):322 336, 2014.
[23] Paul Tupper. The relation between approximation in distribution and shadowing in molecular dynamics. SIAM Journal Applied Dynamical Systems, 8(2):734-755, 2009.
[24] Laurent Dinh, Jascha Sohl- Dicks-Dis,,,,, and Samy Bengio. Density estimation using Real NVP. In Interna- tional Conference Learning Representations, 2017.

---

[25] Conor Durkan, Artur Bekasov. lain Murray, and George Papamakarios. Neural spline flows. In Advances in Neural Information Processing Systems, 2019.
[26] Rianne van den Berg. Leonard Hasenclever. Jakub Tomezak, and Max Welling. Sylvester normalizing flows for variational inference. In Conference on Uncertainty in Artificial Intelligence, 2018.
[27] Sameera Ramasinghe, Kasun Fernando, Salman Khan, and Nick Barnes. Robust normalizing flows using Bernstein-yy type polynomials. arXiv:2102.03509 2021.
[28] Sergei Yu Pilyugin. Theory of pscudo-orbit shadowing in dynamical systems. Differential Equations, 47 (13):1929-1938) 2011.
[29] Abbas Fakhari and Helen Ghane. On shadowing: ordinary and ergodic. Journal of Mathematical Analysis and Applications, 364(1):151-15) 2010.
[30] Robert Corless and Sergei Yu Pilyugin. Approximate and real trajectories for generic dynamical systems. Journal of mathematical analysis and applications. 189(2):409-423, 1995.
[31] Jonathan Meddaugh. Shadowing as a structural property of the space of dynamical systems. vnXiv:2106001957 2021.
[32] Sergei Yu Pilyugin and O.B. Plamenevskaya. Shadowing is generic. Topology and its Applications, 97(3): 253 263. 1999.
[33] Ivan Dokmanić and Rémi Gribonval. Beyond moore-penrose part i: generalized inverses that minimize matrix norms. #rXivv/706.000 2017.
[34] Michelle Schatzman. Numerical Analysis: a Mathematical Introduction. Oxford University Press, 2002.
[35] Cornelius Lanczos. iteration method for the solution of the eigenvalue problem of lincar differential and integral operators. 1950.
[36] Ed S. Coakley and Vladimir Rokhlin. A fast divide-and-conquerr algorithm for computing the spectra of real symmetric tridiagonal matrices. Applied and Computational Harmonic Analysis, 2013.
[37] Juntang Zhuang, Nicha C Dvomek, sekhar tatikonda, and James Dunean. MALI: memory efficient and reverse accurate integrator neural ODEs. In International Conference on Learning Representations, 2021.
[38] Juntang Zhuang, Nicha Dvomek, Xiaoxiao Li, Sekhar Tatikonda, Xenophon Papademetris, and James Duncan. Adaptive checkpoint adjoint method for gradient estimation in neural ODE. In International Conference on Machine Learning, 2020.
[39] Ricky T. Q. Chen, Yulia Rubanova. Jesse Bettencourt, and David Duvenaud. Neural ordinary differential cquations. In Advances in Neural Information Processing Systems, 2018
[40] Patrick Kidger. On neural differential equations. urXiv:2202.024333 2022.
[41] Yu Qiao and Nobuaki Minematsu. study on invariance of f-divergence and its application speech recognition. IEEE Transactions on Signal Processing, 58(7):3884 3890, 2010.
[42] George Birkhoff. Proof of the ergodic theorem. Proceedings of the National Academy of Sciences, 17(12): 660, 1931.
[43] Tanja Eisner, Bálint Farkas, Markus Haase, and Rainer Nagel. Operator Theoretic Aspects of Ergodic Theory. Graduate Texts in Mathematics. Springer, 2015.
[44] Brian Coomes, Hüseyin Koçak, and Kenneth Palmer. Rigorous computational shadowing of orbits of ordinary differential equations. Numerische Mathematik, 69:401-421. 1995.
[45] Joseph Frederick Elliott. The characteristic roots of certain real symmetric matrices. Master's thesis, University of Tennessee- Knoxville, 1953.
[46] Sergey Yu Pilyugin. Spaces of dynamical systems. de Gruyter, 2d edition, 2019.
[47] Tomas Geffner and Justin Domke. MCMC variational inference via uncorrected Hamiltonian annealing. In Advances in Neural Information Processing Systems, 2021.
[48] Guodong Zhang, Kyle Hsu, Jianing Li, Chelsea Finn, and Roger Grosse. Differentiable annealed importance sampling and the perils of gradient noise. In Advances in Neural Information Processing Systems, 2021.
12

---

[49] Tomas Geffner and Justin Domke. Langevin diffusion variational inference. In International Conference Artificial Intelligence and Statistics, 2023.
[50] Achille Thin, Nikita Kotelevskii, Alain Durmus, Maxim Panow, Fric Moulines, and Arnaud Doucet, Monte Carlo variational auto-encoders. In International Conference Machine Learning, 2021.
[51] Heikki Haario, Eero Saksman, and Johanna Tamminen, An adaptive Metropolis algorithm. Bernoulli, pages 223-242, 2001.
[52] Athanasios Tsanas, Max A. Little, Patrick F. McSharry, and Lorrane O. Ramig. Accurate telemonitoring of parkinson's ' disease prognession noninvasive speech tests. IEEE Transactions Biomedical Engineering, 57:884-892 19384-893,22 2010. https:///ppp remanti httpa///eerraaooooooo/r////raa.
[53] Sérgio Moro, Paulo Cortez, and Paulo Rita. data-driven approach predict the success of bank celemarketing, Decision Support Systems, 62:22-3 31, 2014.
[54] Radford Neal. MCMC using Hamiltonian dynamics. In Steve Brooks, Andrew Gelman, Galin Jones, and Xiao-Li Meng, editors, Handbook of Markov chain Monte Carla, chapter CRC Press, 2011.
[55] Kai Xu, Hong Ge, Will Tebbutt, Mohamed Tarek, Martin Trapp, and Zoubin Ghahramani. Ad- VancedHMM..j: robust, modular efficient implementation of advanced HMC algorithms. In Symposium on Advances in Approximate Bayesian Inference, 2020.
13

---

A Proofs
A.1 Layer-wise error analysis
The layer-wise error analysis follows by recursing the one-step bound on the difference between the and numerical flow:
FN .... Fix - Fx .... Firl
I FN D... Fix Fx u FN-1 v.. Frail + | FN a FN-1 o... Fix Fi o... Rall
- Lip(Fx)]|xx-1o... Fix Fx-1 .... Fual +
N-1
S
Σ II Lip(E)). N=0 2
A.2 Error bounds for normalizing flows
Proof of Proposition By triangle inequality.
Let S be the random initial state of the orbit that shadows the pseudo-orbit of X qo-
S -B/(x)) | (E/(s))-B((xx)) B{((8)
By e-shadowing and E-Lipshitz continuity of
Bf(5w) - BB((xx) E [(5x) -5((m) -(5x))
s SEB|SSN Exll s lc.
Next
[ff((N) RS((x)]
< sup [B((x) - Ef(Sx)] = 2U (NN,SN) 2U'D (X, AiISO
The last equality
the fact that X Fy o ....o Fi, and the fact that the total variation is invariant under bijections (41 Theorem 1].
are map of S, respectively, under bijection
Proof of Theorem
4.3 Let E X be the initial state of the backward shadowing orbit. By triangle inequality,
log ((x) log g (x)] I log q(x) log g(s)] log ((x) log (())).
For the first term the right-hand side, note that for y on the segment from
llg ((>) log [(8) = v log log - s) se sup
V [l
(7)
For the second term,
log (zx) log q(s) [log (+-xx) log (o x)]
Σ
log JN (a-1--)((--n)) log JN (n- ()(_)()))
We apply the same technique to bound each term as in Eq.
Proof of Theorem By definition, ELBO (qllp) [x))), and by hypothesis the distribution go of the initial flow state equal distribution Eo of the shadowing orbit initial state, so EE [E(X)] = [E()].
Hence ELBO (allp) - E [(x)] < E(S) - B Finally,
E(s) B(x))
(ogp((N) - log (s) + Σ log Jn sm())))
(1og((()) - log 90 (z
Σ log Jn (in)
< [logp(sx) log((xx) + Blog go (s) log (x)) +
log Jn (sn) log Jn (ft))
==1
The proof is completed by bounding each difference with the local Lipschitz constant around the pseudo-orbic times the shadowing window size c, and applying the expectation.
L
14

---

A.3
Error bounds for MixFlows
Proof of Proposition 4.5
Let S be the initial point of the random shadowing orbit. By the definition of
MixFlows,
N
Ef(Y)
Σ
N+1
N
N
1
f(X.) - f(S.)
+E
Σ f(S.) f(8m)
N +
< B
Σ f(x.)-((SA)
Σ f(S.) x(Xa)
-(Xa)-
To bound the second term, apply shadowing and the Lipschitz continuity of f to show that each f(S.)- - -Xx,) < Ee, and hence their average has the same bound. For the first term, the case where N fixed and finite,
sup LiSSSCLLiPCCS
Σ ((x)) ((S.)
Σ
sup AISULLIP(S
[E - f(S.).
N+1
N ++1
Σ
20
sup [B[((x.) f(S.))
N +
f-lsis!
n=D
N
DUD (Xn, S.)
N +1
m=D
N
Σ
2UD (X,S)
N
M=D
2UDry (Eo.0o).
We now consider the large-N limiting case when F is crgrooi and measure- measure-pressrvingg & T and Eo < n. and U f= 1. Let Z Sk, K ~ Unif(0, 1,..., N.. and let W ~ In this case, by the triangle incquality,
N..
sup FUSssUuLiphss
E
- Σ 1(x)) S(S)) S DBL (X, + DBL (2, W).
By S
Theorem 4.1]. both terms on the right-hand side converge to as N + oc.
Lemma A.1. ff [m], a. b, >> 0, then
01 + + an
a2
log
max Isi5
log
+ + bn
Proof. Define A: Σ ac,
+
bi + ... + bn
ag by
log bi + .... bn
- log
log
+
+ an
A as
Jensen's inequality yields that
11
log bi +
log
max I5i5
log
bn
A
b
by
Applying the same technique to the other direction yields the result.
Proof of Theorem
4.6
First, shadowing and the triangle incquality yields
log ((x) - log ((x)] S log ((x) log q(s) + log ((x) - log ((8)
I e ( + log ((x) - log q())
By Lemma A.I have
log 40 (0) II |=1 J (
log ((x) log (())
max 0SwSNN
qu(s- (o(Ssnn/ 3/ 1ز ((-
max N
|log g. log (s))
15

---

where 4n is the density of the length-n normalizing flow with constant forward flow map F, and 4r. is its approximation via the pseudo-orbit. Using the same technique in the proof of Theorem 4.3
max (n (n) S USRSN
(Ten.c(fnn) Σ (.
USaSN
S
max OSPSN
Lg..
Σ (R-n.)
Combining this with earlier bound yields the first stated result. obtain the second result infinite shadowing setting, we repeat the process of bounding maxosagN [log 9n (I) - log gn (a)1, but rather than expressing each supremum around the pseudo-orbit ê-n, we express around the shadowing orbit 8-n. to find that
Σ Lan(e-))
max [log Gn - log om (+)) c OSPSN N
max OSHSN
Lgo (s_.)
We then bound maxogagN with a sum and merge this with the first bound to find that
1 lim sup N
log ((x) - log ((>)] S e lim sup
N [aae(() + Lg (s-n) 1
Σ Lsse(sn))
Since (x) = ((N) as N + oo, the first term decays to O. The latter two terms are orgodic averages under the backward dynamics initialized at if the pointwise ergodic theorem 42). 43
212] holds at then the result follows. However, the pointwise ergodic theorem holds only r-almost surely for each of Lque and L.ee- Denote Z to be the set of m-measure zero for which the theorem does not apply. Suppose there is set A [ X' such that (o(AA) > O, and x E A implies that E Z. But then Eo (2) > O, which contradicts & < m. Therefore the pointwise ergodic theorem applies to s for qo-almost every E X.
Proef of Theorem 4.7
By definition, ELBO (gllp) E. [E(X)] and by hypothesis the distribution of the initial flow state is equal to the distribution εo of the joint shadowing orbit initial state, so E [E(X)] [ES)
Hence ELBO (allp) E
[(x)
$ E(S) - É(X) |]]
So applying the triangle inequality, we find that
N Z5=0 1را J(in-s) yn (0 Il ((n
q0(A,
((a) - ((()) s
1
log
p(sn) p(t...)
+
log
N+1
N
Σj= ==0
N
For the first sum, each term is bounded using the local Lipschitz constant,
log p((n)
N
Σ
logp(sn) p(8m)
N+1
p(ān)
N
W+1
A=0
N
Σ
(fn))
N+1 7=0
For the second, apply Lemma A. and then the local Lipschitz constants,
Σ (0( N+1
N+1 N (0(سn_)))
log
-==0 II: ((n-e)
1
Σ
log
((n
N
( N+1 Σy= mi s(nn-s)
N+1
Σe= Il J(an-s)
go(
== N+1
(2 IP- J(an-1)
N
Σ
max USYSN
log
90 (0( I J(nn-s)
N max OSISN
flog loggo((n--)) log (00(0mn))))) Σ
blo ((En-1) log ((n-
N W1
n=D
max L (fn))(nn)))))))))))))))())) + Σ ( OSIN
6
max - -NSASN
L (2)))(()))
Ls.e Li
N+1
Lon. +
Σ L.se(tn))
- NSS
16

---

Combining this with the previous bound and taking the expectation yields the first result,
ELBO (qllp) - E [E())] c
E (N)) n=0
max Loo Xm -NSAEN
Σ
(xne
N
n=-N
Once again to arrive the second result, we redo the analysis but center the supremum Lipschitz constants around the shadowing orbit points, and replace the maximum with sum, yielding the bound
N
EIBO (gllp) - [(X) E.
E
(Sn) + max -NenSN
(S.)
Σ
N
T==0
n=-N
E..E
(Sn) + (S,) + Σ Ls,e(Sn)
Now since 4o by assumption, S..
d Xn for each N.. so
N Lp (X.) +
N
ELBO (allp) E [{x)] E
L Lm (X.) + E (Xn)
N
n=0
n=-N
n=-N
We divide by N and take the limit supremum:
limsup 1 N-000
ELBO (qllp) - E [Bx)]]
N
N
e semmsup RR.
Σ
(xp) + Lan
(Xx) Tx.((Xa)
Σ
N
N(N +1)
N
n=0
n=-N
n=-N
Consider just the term
= /
dgo (x)m(dz)) Σ La (F"
E
Lxc(Xn)
Lic(F" x)go(dx)
N
N n-u
N
dm
m=0
Theorem 8.10 (vi)] asserts as long as Lif. L (TT) and dae < oo, then
L==² for 1/mm + 1/m²
da
< 1/
lim N-
E (F
dgo dm (z)=(dx)= = (s.e(t))((r)) ( /
dgo (dr) [X)e(X)), X
N
dx
Applying this result each above yields the stated result.
A.4
Norm of matrix right inverse
Notice that the operator A in Theorem 5.1 has following matrix representation:
Dr 1
D2
T
A
R RdNxd(N+1)
where DA VK (8)
-DN 11
Proposition A.2 133
Corollary 4.21). Let as defined in Eq. ( Suppose (Fr) are differentiable.
Then
At = A (AART-- arg min |XII, subject to AX =
(9)
and ||A||| Omin (m)) where aunin (A) denotes the smallest singuaa value of A.
Proof of Proposition A.2 Note that A has full row rank due to the identity blocks. In this case, AAA" is sinvertible, a valid right inverse A. Eq. (9) is then obtained by direct application Corollary
thus A (AA") 4.2].
To obtain the norm At, since A has full row rank.
A= T.VV where U RdNxdN V xx((++))x((++1 are orthonormal matrices,
Σ E RdN xd is reetangular diagonal matrix with full row rank.
17

---

Then
||4 VE"u" (UNV"VS'O)) |
=|<<E>">~11
(A)
Notice that since AA a is an invertible Hermitian matrix, miin (A) is strictly positive.
B Discussion about the shadowing window e
B.1 Numerical error when computing
As mentioned Section due the floating-poin point error involved evaluating Fe. the evaluation perturbed by the numerical error. This consequntely leads to an inaccuracy of the estimation of À = ||4,- I|- Here A, denotes the right inverse of A. In this section, explain how to calibrate the calculation of À to take into account of this error.
This computational challenge. rooted in numerical inaccuracies, previously acknowledged in Coomes al. 17 Section 3.4], with a provided remedy. Suppose we compute terms in (4 with floating-poin error 8. This introduces O(VN') 8) Frobenius norm error into the practically calculated A. Here denote the digital computed A by A, À := ||A," " I|. and = A - AlI. According to Coomes et al. the following upper bound A:
Eq. 45), we can employ
s - - - - (1=0))))))
Substituting values of 6, N, and computed À from experiments (refer Fig. this amounts an inconsequential relative error of around (1 10 10-9)-1
A.
is also worth discussing the the shadowing window computation methods proposed in Coomes al. [17 Section 3.4]. As opposed to our generic eigenvalue computation introduced in Section the computation procedure presented in Coomes et al. 17 Section 3.4] can potentially offer better scalability. However, the
methodology in Coomes et al. [17] operates under the heuristic that dynamical system with shadowing intrinsically hyperbolic,
Specifically, the "hyperbolic splitting" discussed in Coomes al. [44] P.41,, Appendix B] (i.e., choice of &), and the "hyperbolic threshold" choice in Coomes al. P.415, Appendix B] (i.e., the choice of p). both implicitly assume that the dynamics either exhibit hyperbolicity or tend closely to it. Such presuppositions might be overly restrictive for generic variational flows, given that the underlying dynamics general time- time-inhomen inhomogencous systems. When the underlying dynamics are not (nearly) hyperbolic, applying strategies from Coomes et al. (I7] could potentially lead to poor estimation of the shadowing window.
nearly
B.2 Dependence of e on N
One may notice that depends N implicitly through (AA") Recall that can be represented as matrix
D1
I
D2
A
E Ie (
where Dr (. [N].
DN T
And the shadowing window is given by 28Amin ( (AATT) If Amin (AA")-11 scales badly with N, our O(Ne) error bounds become vacuous. Therefore, in this section, we collect some insights about the scaling of over increasing N. It is important to that understanding the explicit quantitative relationship between and general nonautonomous dynamical systems quite challenging. Although current literature dynamical systems focuses on studying conditions so that e does not scale with N, these conditions are typically difficult to verify in practice are only applicable in restrictive settings such as men-hhoogeneouu linear dynamical systems. For more systematic study, [17 20 28 The examples provided below merely serve synthetic illustrations for intuition and should not he taken as definitive statements for all cases.
Example B.1 (Scaling map in R.). Consider the timehomogeneeus dynamics formed by T RR R such that Ya R, T(z) Ca for some > 0. Then the corresponding shadowing window of a ô d-pseudo-orbit
e 26 (C- 1)a 2 2C [1-
(1- cos
N N
18

---

Therefore. for all C. o 1, as N + =x,==0(00 For C I, ((8).. The formula for Amin (AA?) above arises because in this case, AA has the following symmetric tridiagonal form,
C² +1
-C
C
C
-c
AA
+++1 -C
-C
C² +1
whose eigenvalues have the following closed form 45
Page. 19-20]; L=.....,
>x(AA") C² + 2C cos
km
N+1*
Example B.2 (Uniform hyperbolic lincar map[46] Theorem 12.3 and p. 1711). Consider the time- hime-hooogeneuus linear dynamical system in R", T : x + Ax, where A E R gaxd is similar to a block-diagon matrix
diag(A) (A1,(AAA)))) A2) such that [|A || < À and ||A2 1 S X for some E (0,1). Then for all N 14, any ö o-pseudo-
orbit is e-shadowed, with c 1+
1-X 6.
As illustrated in Examples B.1 and B.2 it is possible for the e-shadowing property to hold in way that does not depend on the flow length N and remains in similar magnitude as 6. This phenomenon is referred
the infinite shadowing property 177 L9 32
mentioned in the heginning of Section [33 The infinite
32 meaning that dynamical systems exhibiting this property form a dense subset of the collection of all dynamical systems. However, while genericity offers valuable insights into the typical behavior of dynamical systems, it does not imply that all systems possess the property. As previously explained, only small class of systems has been shown to have the infinite shadowing property, and verifying this property for specific given system remains challenging task.
shadowing property is shown to be generic in the space
dynamical systems
C
Normalizing flow via Hamiltonian dynamics
typically builds flexible flow family by composing numerous parametric flow maps that are agnostic to the target distribution, recent target-informed variational flows facus on designing architectures that take into account the structures the target distribution. common approach for constructing such target-informet flows is to use Langevin and Hamiltonian dynamics as part of the flow transformations 6 7 47.50. In this section, hegin with concise overview of two Hamiltoninn-asse -based normalizing flows: Hamiltonian variational auto-encoder (HVAE) 7 and Hamiltonian flow (HamFlow) 66) Subsequently, we present corollaries derived from Theorems [43 and 4.4 specifically applying these general results to HVAE and HamFlow. Finally, we offer empirical evaluations of our theory for HamFlow on the same synthetic examples as those examined for MixFlow.
While past work normalizing flows 2
3
24 26
Hamiltonian dynamics (Eq. (00) describes the evolution of a particle's 's position 0 and momentum Pe within a physical system. This movement is driven differentiable negative potential energy log p(f,) and kinetic energy 3 Pp px:
dp: dt
db:
= V log क (0c)
= Pr.
(10)
dt
For & E IR, define the mappings He 2d 2dd transforming (es) (Bc+s, Pt+a) according
the dynamics Eq. (10). Intuitively, traversing the trajectory of H,, E allows particles to efficiently explore the distribution of interest's support, thereby motivating the use of He as a flow map. In practice, the
exact Hamiltonian flow H, is replaced by its numerical counterpart- integration- involves K
iterations of the subsequent steps:
ه -- p+ V log (0)
0 + 0+
p Ā + - log (0).
2
Here e denotes the step size of the leapfrog integrator. We use Tk,, sequencing the above three steps K times. There are two key properties of Tk., that make it suitable for normalizing flow constructions. Tk,, is invertible T.. T'k..) and has simple Jacobian determinant of the form: det VTx.el = 1. This enables a straightforward density transformation formula: for any density go(0, p) on 2d pushforward of q under Tk.c . has the density qu (Tk,..(0,))
2d - ue2d to define the map by
Given the above properties, Caterini et al. 7 and Chen et 6 developed flexible normalizing flows flow transformations through alternating compositions of multiple Tk,, and linear transformations. Notably, the momentum variable introduced in the Hamiltonian dynamics Eq. 100 necessitates that both HVAE and HamFlow focus on approximating the augmented target distribution
HVAE and HamFlow, respecivell-uutiiaeee
19

---

p) & =0)) exp(p=pp 20'p) on R This represents the joint distribution of the target and an independent momentum distribution. (0, Ia).
Specifically, flow layers of HVAE alternate between the leapfrog transformation Tk,, and momentum scaling layer:
(), + 7p), 7 E IR'
And those of HamFlow alternate between Tk.c and momentum standardization layer:
p0, - A(p - n)), A igd lower triangular with positive diagonals 1 E I8"
C.1 Error analysis for Hamiltonian-bassd flows
We have provided error analysis for generic normalizing flows in Section 4.2 In this section, specialize the bounds for density evaluation (Theorem 4 and ELBO estimation (Theorem [ for the two Hamiltonian-
based flows. я denotes the target distribution of interest on IR"", and p denotes the augmented target distribution on ad p and p respectively denote the unnormalized density of x and R.
Corollary C.I. Let q be either IVVAE HamFlow of N layers with qo ((0,Izaa) Suppose the backward dynamics has (c,6)-sit --shhdwowing property. Then
Vx E R.2d
log (x) log g(x)) S e (Tq. xx1x + +
Proof of Corollary
Since the Jacobian determinant of each HVAE/HamFlow layer is constant,
LJ,, (2) = 0,
Vn E [M] and æ E ad
(11)
And since qo is standard normal distribution,
Log.c
sup ||v|l S ||a||| +
(12)
I S+
The proof is then complete by directly applying Eqs. and the bound stated in Theorem
Corollary C.2. Let be either HVAE or HamFlow of N layers with = -N0,, Iza). Suppose the forward dynamics has the (c, 5'-hddwwiin property, and &u- qo. Then
ELBO (gllp) E[E(X)] sc (L [e[Lac(ix)] + v2d + c² for X go.
additionally log is L-smooth (i.e., Va, V log Vlog((v)] L||x vlli, then
ELBO (q p - [[L()) s. [1V I|]
+ V2d + (LV1+ 1)e²
for X go-
Proof of Corollary C
Similar to the proof of Fig. 2d
applying Eqs.
and to the error bound stated in
Theorem 4.4 yields that
ELBO (gll)) F(ε(X)]
e E
[xx((Xx)
||||| + e²
for X go-
The first result is then obtained by employing the following bound via Jensen's inequality:
E|||X|| [I] <
X P = 2d for ~ 4o = ((.
To arrive the second result. we apply following bound on Lp., using fact that log p is L-Lisschitz continuous:
L..- W V (())( + sup
|V
[Vlogp(z) V log & logg(x)) + LV l.c.
1lv-
C.2 Empirical results for Hamiltonian flow
In this section, we empirically validate our error bounds and the diagnostic procedure of HamPlow using two synthetic targets: the banana and cross distributions. We employ warm-started HamFlow, detailed in 6 Section 3.4]. This is achieved through single pass of 10,000 samples from the reference distribution. go. without additional optimization. Specifically, the μ
values for each flow standardization layer are determined by the sample mean and the Cholesky factorization of the sample covariance, respectively, of the batched momentum samples processed through the layer. For the banana distribution, we set the leapfrog step K 200 and step size of 0.01. For the cross distribution, set K to 100 with value of 0.02. Figs. Sa
20

---

and showcase comparison between 1000 i.idd samples from the warm-started HamFlow with 400 layers and the true targets. The results illustrate that HamFlow generates samples of decent quality.
We start by examining the numerical orbit computaiton error. As demonstrated in Fig. the pointwise evaluation error grows exponentially with increasing flow depth, and quickly reaches the same order of magnitude as the scale of the target distributions. However, similar to what we have observed in the experiments of MixFlow (described in Section [6) the substantial numerical error of orbit computation does not impact the accuracy of Monte Carlo estimations of expectations and ELBO estimations significantly. Fig. TT shows the relative numerical error of Monte Carlo estimates of the expectations of three functions. The relative samples estimation error remains within 5% for the two bounded L Lipschitz continuous test functions, aligning with our error bound in Proposition 4.2 Moreover, even for the unbounded test function, the error stays within reasonable limits: 10% for the banana and 5% for the cross examples. Fig. [1 then compares the numerical ELBO estimates and exact ELBO estimates for HamFllw over increasing flow length. In both examples. the numerical ELBO curve is almost identical the exact one, notwithstanddnn the considerable sample computation error of the forward orbit.
However, unlike the small log-density evaluation error we observed in the MixFlow experiments, HloFllows log-density shows significant issues because of the numerical errors of backward orbits computation. In Fig. [2] we measure the numerical error for the log-density evaluation of HamFlow at 100 i..d. points from the target distribution. As we add more flow layers. the error grows quickly. For the banana example, the log-density evaluation error can even go above 1000. This raises two main questions: (I) Why is the ELBO estimation error different from the log degsienitty evaluation error in HamFlow? (2) Why is the scaling of the log-densiig evaluation error different between MixFlow and HamFlow? Our theory presented in Section 4 and Appendix C.I offers explanation to these questions.
To address the first question, we compare the bounds for HamFlow on log density evaluation (Corollary C.I and ELBO estimation (Corollary C.2). The main difference hetween these bounds is in the local Lipschitz constants.
The error in log-density evaluation is based on T.ge while the error in ELBO estimation is based on TT[T.p... These two quantities- Ty and E. [[Zg.]-relat to the smoothness of the learned log-density of HamFlow and log oog-disstyy of the target. Figs. [8c to Bijdisplay log and log y side-by-sid for both synthetic examples. Indeed, in both cases, log g markedly deviates from log and presents drastic fluctuations, indicating that Ly can be substantially larger than 18 p... Fig. [3 further provides a quantitative comparison on how Lg and [p.c change as we increase the flow length. Note that since both Lgge and Ly., are intractable in general, we instead focus on the scaling of |V log q| (which is a lower bound of Lg.e) and an upper bound of [B [p.. (derived from case-by-case analysis for both the banana and cross distributions; see Appendix CC33 for detailed derivation).
Results show that in both synthetic examples, Lage can increase exponentially over increasing flow length and eventually larger than [[ by 30 orders of magnitude, while B[L,,e] does not grow and remains within a reasonable scale.
Similarly, the differences in log-density errors between MixFlow and HamFlow come from how smooth the learned log density is in each method, according to the error bounds presented in Theorems 1.3orems/43Jaap and [.6 Fig. [ visualizes the learned log densities of MixFlow, which matches closely to the actual target distribution, and has a smoother output than HamFlow. This matches what was noted Xu et al. [5] Figure 3J. where MixFlow was shown to approximate the target density better than standard normalizing flows.
Finally, Fig. presents the single step evaluation error 5 HamFlow and the size of the shadowing window e as the flow length N increases. Results show that for both synthetic examples, is approximately 10-7, hence we fix 8 to be this value when computing the shadowing window e. Compared to the significant orbit evaluation error when N is large, the size of remains small for both the forward and backward orbits (in the scale of 10 = 400). And more importantly, the size of e grows less drastically than the orbit computation error with the flow
to 10 3 when
quadratically in the banana example, and linearly in the cross example, which is similar to what is observed in the MixFlow experiments.
C.3 Upper bounds of T.p.c for synthetic examples
Recall that we denote p(x) & p(())=2) exp(- [1/3 (4]s4)) the augmented target and denotes the 2d- synthetic target are interested.
Notice that
Lp,e (x)
sup |/ log (.)) [x] || < (2))2))) + [|[]]s=a|| +
where Lp.. requires case-by-cee by-eee analysis. Therefore, in this section, we focus on bounding the local Lipschitz constant for logp. p. Proposition C.3. For the 2-dimensional banana distribution, i.e.,
1 x² 1 (-2 0² 2
(z2 bar -ob))
p(x) x exp
21

---

[1 J2
we have that for all =:
E ER²,
Lp.e(x) S [Vlog((>)] +
2bb/z1| + max 6b'22 2î + 2bx2 2000²²
02 1
1
max {6b|z1l ob(e + ,),
Proof of Proposition
By the mean value thcorem,
(x) - || logp(z)) + sup VVlog(y - v logp(>)]
llv-=lIse
A logp(z)]1 + sup
|| logp(v)l
11v-=lls+
The proof will then complete by showing that
sup S
||
Io)
6b + 2hxz 20 + 1 < 26/ + max {
02 -1
+ 2b max {6b|211 | + 6b(e + e,
(13)
To verify Eq.
15 note that
x6 7 + 2hro 202
+/ 2hz
V2 log
26x1
1
The Gershgorin circle theorem then yiclds that
v logp(x)]] S 26/211| max [66229 2bxx
20262
+ σ² .1
Additionally, we have
sup llls
lyn| S |
and
sup 662b) vî I>v=
2hyz
20²6b²
66b22 . 2hrz 20²6²
2bmax {6b|=1
6b(e + )),
02
The proof is then completed.
Proposition C.4. For the cross distribution, 4-component Gaussian mixture of the form
p(x)
1
1
x
2
[0.15²
0
) +
N
-2]
0
0
0.15
0
1
2
1
0
0
0.15
0
+
4
N
0
0.152
+ v
2
0
we have that for R.?
(x. S 0.15² (2 + |>|| e).
Proof of Proposition C.4 P1, P2,P3.PPPP P4, hence
For notational convenience, denote the density of 4 Gaussian component by
logp(x) = log (pi(z) + p2(x) + pa(x) + pa(x)) - log 4.
Examining V p(x) directly yields that for all æ E R²"
V logg(x)
(1) + +2 +s Pa (x)
conv {Σ, ( - () :i= 1, 2,
Here μi; and V, denotes the Gaussian mean and covariance of ps. Then by Jensen's 's inequality,
|| log < max
sup lly-=lis+
|Σ (#: will ; 1,2,3, 4
|ا ||
[|Hi|] +
|Blll
sup llv-allSe
#4
(2 + ||||||| +e).
0.15a
22

---

Trained HamFlow
Trained HamFlow
0
4
-50
2
A Ae
0
-100
-2
-50
True Distribution Trained Flow
True Distribution Trained
-4
40
-2
0
20
40
-4
-2
0 2 4
(a) Banana: HamFlow sample scatters
(b) Cross: HamFlow sample scatter
(c) True log-pdf
(d) HamFlow log-pdf
(e) True log- log-pdf'
(I) HamFlow log-pdf
Figure 8: Visualization of sample scatters and density of HamFlow. Figure (a)-(b) each shows 1000 i.i.d. draws from HamFlow distribution (red) and target distribution (green). (c) and (c) show respectively the exact log-density of the Banana and the cross distribution, while (d) and (f) show respectively the sliced HamFlow log-density evaluation at p (0, 0) on those two targets. The HamFlow lon-iensity evaluations are computed via 2048-big BigFl oat representation.
HamFllow numerical
HamFFow numerical error
HamFlow numerical
HamFlow numerical erro
10
10.0
Error 1088
10-3
7.5
5.0
10
2.5
0.0
100 200 300
400
100 200
300
400
100
200
300
400
100 200 300 400
#transfrrmations
#transformations
#transformations
@transformatiios
(a) Banana
(b) Banana (log-scale)
(c) Cross
(d) Cross (log-scale)
Figure 9: HamFlow forward (fwd) and backward (bwd) orbit
synthetic examples. (a)
and (c) show the median and upper/lower quartile forward error - f and backward error |Bay y comparing & transformations the forward exact/approximate maps F, F or backward exact/approximate maps B, B. And (b) and (d) respectively display (a) and (c) in a logarithmic scale, intended to provide clearer illustration of the exponential growth of the error. The lines indicate the median, and error regions indicate thth to thth percentile from 100 independent runs, which take an i..dd. draw of from yo and draw y from the actual target distribution p.
HamFlow sampling error
HamFlow sampling error
Hamflow sampling error
HomFlow sampling emror
0.100
10
0.04
0.075
0.03
2/
0.050
10
0.02
10
0.025
10
0.01
10
90142
100
200
300
400
100
200
300
400
100
200
300
400
100
200
300
400
#transtormations
stransformations
stranstormations
stransformations
(a) Banana
(b) Banana (log-scalc)
(c) Cross
(d) Cross (log-scalc)
Figure 10: HamFlow relative sample average computaion error three functions:
Ea Ei (in(() + 1), and [sggood((()) Each line denotes the numerical error computed from Monte Carlo average of the three test functions via 100 independent forward orbits starting at i.i.d. draws from 90-
23

---

HamFlow ELBO est.
hamflow elbo est.
-10
-20
EL
elbo
5
-30
6
-40
numerical exact
-7
numerical exact
-50
0
100
200
300
400
100
200
300
400
#transformations
#transformations
(a) Banana
(b) Cross
Figure 1I: HamFlow exact and numerical ELBO estimation over increasing flow length. The lines indicate the averaged ELBO estimates over 100 independent forward orbits. The Monte Carlo error for the given estimates sufficiently small so we omit the error bar.
Hamfloow low log-density error
Hamfllw log density error
Hamfllow log-density error
Hamf low log-density enror
10
1500
16
error
Encer
1000
Errrr
10
500
10
20
100
200
300
400
100
200
300 400
100
200
400
100 200 300 400
#transformations
*transfarmations
#transformetions
#trensformsstion
(a) Banana
(b) Banana (log-scale)
Cross
(d) Cross (log-scale)
Figure 12: HamFlow log density evaluation error on synthetic examples. Log- Log-densities are assessed positions of 100 independent samples from the target distribution. (b) and (d) respectively display (a) and (c) in logarithmic scale for better illustration of the growth of error over increasing flow length. The lines indicate the median, and error regions indicate 25th th percentile from the 100 evaluations.
HamFlow [|Voog]
HamFlow bound I,,
HamFlow |Voooll
HamFlow bound T.p.
10
260
10
10
240
260
10'0
220
200
100 200
300
400
100
200
300
400
0
100
200
300
400
100
200
300
400
Atransformation
#transformations
#transformations
#transformations
(a) Banana
(c) Banana
(a) Cross
(c) Cross
Figure 13: Comparison of HamFlow's Lg.e and B[Lp.. x)) on the banana and cross distributions.
Figs. 13a and [3c plot |/V log q==)]l against flow length N; the lines indicate the median, and error
regions indicate t5th to 75th percentile from i.i.d. 100 independent draws of z from p. Figs. 13b and plot estimated upper bounds of E E[ (xx) (averaging over 100 independent draws from g) against flow length N. The choice and detailed derivation of the upper bounds of Lø,e can found in Appendix C.3
D Additional experimental details
All experiments were conducted on a machine with an AMD Ryzen 9 3900X and 32G of RAM.
D.1 Model description
The two synthetic distributions tested experiment were
the banana distribution
y
y1
N
0,
100
0
9))
V1 32 + byi
1006
b= 0.1;
24

---

HamFlow single map err.
HamFlow window size
e computation time
0.0015
0.15
Fodd Bwd
Fad 3wd
10
0.0010
0.10
Error
10-9
0.0005
0.05
0.0000
0.00
100 200
300 400
0
100
200
300
400
Fwd er::
Bwd err.
#transformatiios
#transformations
(a) Banana
(b) Banana
(c) Banana
HamFlow single map err.
HamFlow window size
computation time
10
0.000100
0.06
0.000075
Bwd
Bed
Error
10" 5
0.000050
0.04
0.000025
0.02
0.000000
0.00
10
100 200 300 400
100
200
300
400
Fwd err.
Bwd err.
#transformations
#transformations
(d) Cross
(c) Cross
(1) Cross
Figure 14: HamFlow forward (fwd) and backward (bwd) single transformation error 8 on both
synthetic examples (Figure (a) and (d)), shadowing window size over increasing flow length (Figure (b) and (e)), and shadowing window computation time (Figure (c) and (f)). Figure. (a) and (d) show the distribution of forward error {11 F F=-2xoll
I,..., NJ (Fwd err.) and
backward error {1L=TTTO B pm+t1a zoll L...,,,.N (Bwd err These errors computed
along precise forward trajectories that originate at zo, which is sampled 50 times independently and identically from the reference distribution go- N here denotes the flow length for each example. Figure (b) and (e) show c over increasing flow length. Figure (c) and (f) show the wall time of computing shadowing window over increasing flow length. The lines in Figure (b)-(c) and (e)-(f) indicate the median, and error regions indicate 25th to 75th percentile from 50 runs.
a cross-shaped distribution: Gaussian mixture of the form
1
0.15 0
v
1
0
1
0
0.15
0.15
0
0.153
0
)
1
The two real-data experiments are described below.
Linear regression. We consider a Bayesian linear regression problem where the model takes the form
Lin. Reg.: B Lid. (0. log ltd (o,,,),, M/ 1 indep N (=)))
where y is the response and x, C R2 is the feature vector for data point j. For this problem, we use the Oxford Parkinson's 's Disease Telemonitoring Dataset [52]. The original dataset is available at hatt::/arrclive ice uci composed of 16 biomedical voice measurements from 42 patients with staly-stage Parkinson's 's disease. The goal use these voice measurements, as well as subject age and gender information, to predict the total UPDRS scores. We standardize all features and subsample the original dataset down to 500 data points. The posterior dimension of the linear regression inference problems is 21.
Logistic regression. We then consider Bayesian hierarchical logistic regression:
Logis. Reg.: ~ Gam(1,0.01), | NN0 -1r,, yi
indep ~ Bern
53 downsampled to 400 data points. The original dataset is available The goal is to use client information to predict subscription term deposit. We include features from dataset: client age, marital status, balance, housing loan status, duration last contact, number contacts during campaign, number days since last contact, and number contacts before the current campaign. For each of the binary variables (marital status and housing loan status), all unknown entries are removed, All included features are standardized. The resulting posterior dimension of the logistic regression problem is
We use a bank marketing dataset
bteps://archive.. ics. www.uci.e.....................i. htt/al//daaasets/aaaaaaaaa
25

---

inversion error
inversion error
errer
Inversion error
10
10
10
10.8
on
Ener
10 49.9
10
10
10 30
50.0
800
(a) Banana
(b) Cross
(c) Lincar regression
(d) Logistic regression
Figure 15: MixFlow inversion error when using 2048-bit BigFloat representation. Verticle axis
shows the reconstruction error || FN (ro) zoll (Fwd) and || BB (ar) - (Bwd) for
and the two synthetic examples (Figs. 15a and 15b and the two real data examples (Figs. [1c and [5d). F and B are implemented using 2048-bit BigFloat representation, and zo is sampled from the reference distribution qu- The lines indicate the median, and the upper and lower quantiles over 10 independent runs. It can be seen that all four examples, the inversion error of MixFlow ignorable when using 2048-bit computation.
MixFlow numerical error
MixFlow numerical al emor
MixFlow numerical error
Mixflow numericall error
Enor SuNees
Fwd
Fwd
Dwc
0wa
100 200 300 400 500
100 200 300 400
500
100 200 300
400
100 200
400
500
#transformations
#transformatios
*transformetions
#tronsformations
(a) Banana
(b) Cross
(c) Linear regression
(d) Logistic regression
Figure 16: MixFlow forward (fwd) and backward (bwd) orbit errors on both synthetic examples and
real data examples. (a)-(d) show the median and upper/lower quartile forward error || - F* x
and backward error | z BlI comparing k transformations of the forward exact/aprroximate
maps F, f or backward exact/approxiimat maps B, B. Note that Figs. lb le Icaad respectivvy display the initial segments of (a)-(d) in a logarithmic scale, intended to provide a clearer illustration of the exponential growth of the error.
D.2 MixFlow parameter settings
For all four examples. use the uncorrected Hamiltonian MixFlow map F composed of discretized Hamil- tonian dynamic simulation (via lcapfrog integrator) followed by sinusoidal momentum pseudorefreshment as used in Xu et al. 5
(see ( Section 5, Appendices E. and E.2] for details). However, we use a Gaussian momentum distribution, which is standard in Hamiltonian Monte Carlo methods (54) as opposed to the Laplace momentum used previously in MixFlows [3]. In terms of the settings of leapfrog integrators for each target, we used 200 leapfrog steps of size 0.02 and 60 leapfrog steps of size 0.00S for the banana and cross target distributions respectively, and used 40 leapfrog steps size 0.0006 and 50 leapfrog steps of size 0.002 the linear regression and logistic regression examples, respectively. The reference distribution gu for each target is chosen to be the mean- mean-field Gaussian approximation as used in (5)
As for the NUTS used for assessing the log-density evaluation quality two real data examples, we use the Julia package AdvancediMM... j1 (55) with all default settings. NUTS is initialized with the learned mean of the mean-field Gaussian approximation go. and burns in the first 10,000 samples hefore collecting samples; the burn-in samples are used for adapting the hyperparameters NUTS with target acceptance rate set to be 0.7.
D.3 Additional MixFlow results
26

---

MixFlow log-densily enrer
Mixflow log-densiity error
MixFlom log log-density error
low logrdensily etror
0.4
1.00
Erter
0.50
1.0
0.2
0.25
0.5
1000g0--g00 soRgo0gooooo go0
200
0000g00000000 g09g0000ggg00000
0000 y00g000g0 @transtormations
#transtormations
#transtormations
ttranstormations
(a) Banana
(b) Cross
(c) Lin. Reg.
(d) Log. Reg.
Figure 17: MixFlow log-density evaluation error. For synthetic examples (Figs. 17a and 17bp. we assessed densities at 100 evenly distributed points within the target region (Figs, and 2b For real data examples (Figs. and [7d), we evaluated densities at the locations of 100 NUTS samples. The lines indicate the median, and error regions indicate 25th 75th percentile from 100 evaluations on different positions.
MixFlow single map em:.
MixFlow single map err.
MixFlow single map er:
MixFlow single map etr.
10-11
10-11
Bror
10
10
10-12
g 10
10
10"
10
10 15.0
10
10 15
Fwdd ert.
Bwd err.
Pwd err
Bvid ert.
Fwd err.
Bwid
Fwd err.
Bwd err.
(a) Banana
(b) Cross
(c) Linear regression
(d) Logistic regression
Figure 18: MixFlow forward (fwd) and backward (bwd) single transformation error on both synthetic examples and real data examples. (a)-(d) show the distribution of forward eπror |IF"rO F pm-1zo zol I.....N] NN] (Fud err.) and backward error
{1 B p=+1zoll I,...,Nyj (Bud err.). These errors computed along precise forward trajectories that originate at Tp, which is sampled 50 times independently and identically from the reference distribution qo. N here denotes the flow length for each example.
computation time
computation time
computation time
computation time
Woll time in sec
0.3
0.6
Avei
800
Eed
80
Fdd
0.2
Well time
0.4
400
40 2
0.2
0.0
2000000 300 600
0.0
200
400
600
800
Wall
0000g0000000 1000000000000000000 40Pg0P&00
Woll
100 209 100000 40Pg00.00 etrensformations
#tranoformatiors
Atranfformations
(a) Banana
(b) Cross
(c) Linear regression
(d) Logistic regression
Figure 19: MixFlow shadowing window computation time (wall time in second) for both the forward (Fwd) and backward (Bwd) orbits over increasing flow length. The lines indicate the median, and error regions indicate th to 7h percentile from 10 runs. Note that here the Monte Carlo error between different runs is so small that the error bar is too thin to be visualized.
27


Reviews:
[{"review_id": "nAoKA5coNw", "reviewer": "Reviewer_MECE", "summary": "The paper investigates the impact of numerical instability on the reliability of sampling, density evaluation, and evidence lower bound (ELBO) estimation in variational flows. It demonstrates through empirical examples that numerical instability can lead to deviations in the flow map affecting sampling, density, and ELBO computations. However, the paper also finds that despite numerical instability, results from flows can still be accurate enough for practical applications. The paper treats variational flows as dynamical systems and uses shadowing theory to provide theoretical guarantees on the error of downstream tasks. It also develops a diagnostic procedure to validate results produced by numerically unstable flows in practice. The paper concludes by validating its theory and diagnostic procedure on MixFlow with both synthetic and real data examples.", "questions": "* It appears that the exponential scaling in the density evaluation is circumvented by looking at *log*-densities?\n\n* Can one motivate the assumptions (in particular, $\\xi_0=q_0$) for the ELBO result?", "limitations": "Details on limitations are lacking.", "rating": 6, "confidence": 3, "soundness": 2, "presentation": 3, "contribution": 3, "strengths": "* The paper addresses an important and practical problem in the context of variational flows, which are widely used in generative modeling and probabilistic inference.\n* The empirical examples and experimental validation provide concrete evidence for the theoretical claims made in the paper.\n* The use of shadowing theory from dynamical systems provides a rigorous framework to understand the behavior of numerically unstable flows.\n* The diagnostic procedure developed in the paper could be valuable for practitioners using variational flows, allowing them to assess the reliability of their results.", "weaknesses": "* Due to rather strong assumptions, the paper's claims may not be applicable in all scenarios. The experiments focus only on a few datasets as well as on a specific type of flow (MixFlow). Thus, they may not fully capture the behavior of other types of variational flows or other datasets, which limits the generalizability of the findings.\n\n* It should be explained in more detail what novelty the current submission offers compared to the work of [23] (except for an extension to the evaluation of densities and ELBO estimations). ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for reviewing our manuscript! Please see our point-to-point respoonse as follows.\n\n> “Due to rather strong assumptions, the paper's claims may not be applicable in all scenarios.”\n\nPlease see our general response 5 for a comprehensive discussion of the assumptions. \n\n> “The experiments focus only on a few datasets as well as on a specific type\n> of flow (MixFlow). Thus, they may not fully capture the behavior of other\n> types of variational flows or other datasets, which limits the\n> generalizability of the findings.”\n\nThank you for the comment! Although the focus of this paper is theoretical, we are happy to include more experiments in the camera ready; see our general response 4 regarding this comment.\n\n> “It should be explained in more detail what novelty the current submission\n> offers compared to the work of [23] (except for an extension to the\n> evaluation of densities and ELBO estimations).”\n\nPlease see our general response 1 for a detailed comparison to [23].\n\n> “It appears that the exponential scaling in the density evaluation is circumvented by looking at log-densities?”\n\nThis is a very good point! We did not adequately comment on this in the original submission, and will do so in the camera ready. There are two cases to consider here – absolute error and relative error. \n\nIf one cares about the absolute error of the density, then our result suffices without much additional work. Densities tend to be bounded, while log-densities become unbounded in the tails. Generally, if two distributions are pointwise close in log-density, they will be pointwise close in density as well. Things become slightly trickier for unbounded densities, but we will stick with bounded here for simplicity.\n\nIf one cares about the relative error of the density, then indeed our result suggests that exponential growth in the relative error is possible. But in this case, this is again a major improvement on what you would expect; recall that the empirical phenomena we observe is that the *trajectories themselves* diverge exponentially. So if the trajectories diverge exponentially quickly, then the relative error could grow *doubly-exponentially*.\n\nAs a rough example, consider the ratio of two normals: $\\mathcal{N}(0,1)$, and $\\mathcal{N}(x, 1)$. Then if we consider the map $x \\to \\exp(x)$ — a toy model for an exponential divergence in pushforward error—we obtain density ratios like $\\exp(-\\theta^2 + (\\exp(x) - \\theta)^2)$, which grow doubly exponentially in $x$.\n\n> “Can one motivate the assumptions (in particular, $\\xi_0 = q_0$) for the ELBO result?”\n\nThis is a good point, and we will clarify this for the revised manuscript. The assumption $\\xi_0 = q_0$ has been made in past work [23], and was made here primarily for convenience of the analysis (lines 175–177). But it is a reasonable assumption for the following reason. We know that $\\xi_0$ is close to $q_0$ due to shadowing. \nAnd $\\xi_0$ is indeed an implicit function of $q_0$; it’s a fixed point of a twice differentiable function involving the whole trajectory starting at $q_0$, see page 176 of [17] (we do not recommend to solve the fixed point equation in practice). The implicit function theorem yields that $\\xi_0$ is a differentiable function of the numerical orbit (whose initial distribution is $q_0$). So, the relationship between $q_0$ and $\\xi_0$ is not totally pathological. Assuming that $q_0 = \\xi_0$ provides simplicity without compromising the insights of the result.\n\n> Details on limitations are lacking.\n\nPlease see our general response 3 for a comprehensive discussion of limitations.\n\n"}, {"review_id": "nWS5CnEIOz", "reviewer": "Reviewer_c44C", "summary": "The paper discusses the (non)impact of numerical stability concerns when implementing variational flows, specifically, the robustness of sampling, density evaluation or ELBO computation against the chaotic behaviour of numerical implementations of variational flows.\n\nThe results are the following: while small perturbations of the input imply large perturbations of the output of variational flows (Section 3), the quality of samples, density evaluations, and ELBOs are rarely affected. \nThe paper explains this phenomenon with \"shadowing\":\nMore specifically (Section 4), under the assumption of ($\\epsilon$, $\\delta$)-shadowing, the errors of sampling, density evaluation, and ELBO computation behave like $\\epsilon + d_\\text{TV}(q_0, \\xi_0)$, where $\\xi_0$ is the initial condition that shadows $q_0$. \nWhether or not a flow admits shadowing (and how $\\epsilon$ relates to $\\delta$) can be quantified by estimating (and bounding) the operator norm of the inverse of a specific operator that depends on the gradient(s) of the flow (Section 5). Numerical experiments show how the shadowing windows are minimal in typical applications of variational flows, even though numerical inaccuracy in implementing the flow is significant.\n\n", "questions": "- This is not a question but more of a non-technical suggestion for the author(s). There have been discussions about potential ethical problems of the \"Boston Housing\" dataset: https://scikit-learn.org/1.0/modules/generated/sklearn.datasets.load_boston.html. If the corresponding experiment in the paper can be rerun with a different dataset, the longevity of this work may benefit.\n- Out of curiosity: Why are the sampling bounds propositions but the ELBO and density bounds theorems?\n- Please correct the complexity statements as commented in Appendix B.\n", "limitations": "As mentioned under \"weaknesses\" above, one potential weakness of the algorithm derived from Theorem 5.1 could be computational feasibility. From the manuscript, it is unclear how much computational power it takes to implement this algorithm. A reader will benefit from learning about this limitation if the algorithm is prohibitively costly for some applications.\nOther than that, limitations seem to be discussed wherever relevant.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "Analysing the numerical implementation of variational flows through the lens of dynamical systems is exciting because it logically explains a seemingly unlogical phenomenon. On top of that, the paper is very well written: Although some components are rather technical (e.g. shadowing), the manuscript remains easy to follow.\n", "weaknesses": "In general, I like the paper, but Section 5 (\"Computation of the shadowing window size\") leaves some open questions, especially in comparison to the finite-time shadowing Theorem in [33]. In my view, the following points should be addressed before publication:\n\n1. The algorithm for computing $\\lambda = \\|A^{-1}\\|$ is inaccurate: The implementation explained in Appendix B pretends that one can compute $\\nabla F_k$ exactly, but since evaluating $F_k$ suffers from approximation errors, $\\nabla F_k$ must as well. For context, compare Appendix B to Section 3.4 (especially, Equation 45) in [33] (reference in the paper). The inaccuracy is unlikely to be crucial but exists. It can (and should) be corrected. Or have I missed something?\n2. Why does the paper estimate the operator norm the way it does? In other words, why is a call to the default `eigmin` function superior to the algorithm discussed in Section 3.4 in [33]? The paper mentions that the complexity of `eigmin` is $O(d^3N^3)$, whereas the tailored algorithm on p. 187 in Section 3.4 in [33] appears to cost $O(N d^3)$. Please discuss the numerical cost of estimating the shadowing windows in the experiments (for example, by providing estimates of the run-/wall time).\n3. Could the authors please elaborate on how to estimate $M$? Why is estimating $M$ and $\\delta$ problem-specific, and estimating $\\lambda$ not (line 238)? If estimating $M$ is problem-specific, please discuss estimating $M$ for the problems in Section 6. (I have not found this information in Section 6 and Appendix D.)\n4. Sections 5 and 6 seem to focus on the computation of the shadowing window size $\\epsilon$. Does this mean all flows considered in this paper automatically have the shadowing property? Is this true for all variational flows (e.g., those discussed by Papamakarios et al.: neural spline, planar, radial, Hamiltonian, ...)? It would be helpful to discuss the existence of shadowing before computing the window size in the experiments (e.g. by estimating $2M\\lambda^2\\delta$).\n\nAs mentioned above, these points should be addressed before publication. But I expect all questions to be straightforward to answer, and I will increase my score once this is done.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your review and comments! Before we address them, we would like to point out that [33] does not introduce the finite shadowing theorem—[17] does. We assume that when you refer to [33] in your review, you mean [17]. Hopefully we did not misunderstand!\n\n[17]: Brian Coomes, Hüseyin Koçak, and Kenneth Palmer. Shadowing in discrete dynamical systems. In Six lectures on dynamical systems, pages 163–211. World Scientific, 1996.\n\n[33]: Ivan Dokmanić and Rémi Gribonval. Beyond moore-penrose part i: generalized inverses that minimize matrix norms. arXiv:1706.08349, 2017.\n\n> The algorithm for ...\n\nThis is an excellent point! Initially, we felt detailing the errors in computing $\\lambda$ was too nuanced. However, we now agree with your view.  We plan to provide a more comprehensive description of computing $\\lambda$, including handling the small numerical errors that occur there, in the final text.\n\nThe fix is very straightforward, as you suggest – we use Eq. (45) in [17]. And it does not affect the results much: indeed if we compute terms in $\\nabla F_k(x)$ with a floating-point error $\\delta$, this introduces $O(\\sqrt{N}\\delta)$ Frobenius norm error into our representation of $A$, which then translates to an upper bound on $\\lambda$ of roughly \n$$(1- \\delta||\\tilde{A}^{-1}||)^{-1}||\\tilde{A}^{-1}|| = (1 - O(\\sqrt{N}\\delta) \\tilde\\lambda)^{-1} \\tilde\\lambda.$$\nHere $\\tilde{A}$ denotes the digital computed $A$ and $\\tilde\\lambda := ||\\tilde{A}^{-1}||$. \n\nSubstituting values from our experiments, this amounts to an inconsequential relative error of around $(1- 10^{-9})^{-1} \\approx 1$ on our minimum eigenvalues $\\lambda$.\n\n> why does the ...\n\nThis is also a great question!  We didn't claim `eigmin` outperforms the algorithm in Section 3.4 of [17]   (originally in [Coomes95] for general dimensions, so for the rest of the response we refer to [Coomes95]). You're correct that [Coomes95] offers better scaling than a direct call to `eigmin`.\n\nHowever, the method in [Coomes95] operates under the heuristic that a dynamical system with shadowing is hyperbolic, or nearly so. Specifically, the 'hyperbolic splitting' discussed in P.413 Appendix B of [Coomes95] (i.e., the choice of $\\ell$), and the 'hyperbolic threshold' choice in P. 415 Appendix B of [Coomes95] (i.e., the choice of $p$), assumes that the dynamics is hyperbolic or close to it, which are impractical for our variational flows (our dynamics are general time-inhomogeneous systems). Also, as discussed in Appendix B (line 488--492), this method can potentially lead to a significant overestimation of the shadowing window when the dynamics are not (nearly) hyperbolic.\n\nTherefore, we prefer to use a procedure that is simpler and more widely applicable, at the cost of more computation. We found in our experiments that using `eigmin` took just a few seconds in most cases (see our general responses), which was acceptable. We will include the above discussion about the method proposed in [17] in the revised manuscript (likely in the appendix).\n\n[Coomes95] Brian A. Coomes, Hüseyin Koçak, and Kenneth J. Palmer. Rigorous computational shadowing of orbits of ordinary differential equations. Numerische Mathematik, 69:401–421, 1995.\n\n> Could the authors ...\n\nEstimating $M$ can indeed be problem-specific because $M$, by its definition (line 232), requires knowledge of the local Lipschitz smoothness of each flow layer evaluated at the computed orbit. This will necessitate different analyses for distinct flows. A universal approach might involve knowing a uniform upper bound on the 3rd derivative of $F_n$ (see response to Reviewer sEFg). Likewise, determining $\\delta$ demands analyzing the single step error for the given flow.\n\nAs for computing $\\lambda$: we can use standard automatic differentiation tools to differentiate $F_n$ and compute $A$ (up to floating point error). At that point, we just need to compute a eigenspectrum for a symmetric blockwise tri-diagonal matrix. Both computations can be done in a generic way without model-specific derivations.\n\n> Sections 5 and 6 ...\n\nWe appreciate your suggestion, and will include discussion about the existence of shadowing in the revised manuscript! It is not true that all variational flows exhibit the shadowing property. Our Theorem 5.1 states a sufficient condition, specifically, $M\\lambda \\epsilon \\leq 1$ for the shadowing property to exist. The practical verification of the shadowing property requires an estimation of the one-step digital computing error $\\delta$, the local smoothness constant of the flow $M$, and $\\lambda$.\n\nOur manuscript describes the computation of $\\lambda$ (via an eigenvalue problem for a symmetric blockwise tri-diagonal matrix) and the estimation of $\\delta$ (Figure 10). These suffice to compute the shadowing window size. Estimating M can be more difficult, especially for complex flows. While one could analyze the local Lipschitz smoothness of each flow layer, handling the $\\sup$ in the definition of M can be difficult.  One option is to bound  the 3rd derivative of $F_n$ to upper bound the $\\sup$ analytically. However, given the range of results already in our manuscript, we leave this analysis for model-specific applications in future work.\n\n> This is not a question ...\n\nThanks for pointing this out! We were not aware of this concern. Yes, we will re-run our experiments for the camera ready with another dataset.\n\n> Out of curiosity ...\n\nThe sampling bounds were slightly easier to prove compared to the theorems. We can switch them all to “theorem” for consistency’s sake in the camera ready.\n\n> Please correct ...\n\nWill do! \n\n> As mentioned under \"weaknesses\" ...\n\nThanks for your suggestion! We will provide more details regarding computing the shadowing window, with a thorough complexity analysis in the camera ready.  Please also see our general response 2 for a detailed response."}, {"review_id": "AM28UI6V7v", "reviewer": "Reviewer_sFEg", "summary": "The paper \"Embracing the chaos: analysis and diagnosis of numerical instability in variational flows\" investigates the impact of numerical instability on the reliability of sampling, density evaluation, and evidence lower bound (ELBO) estimation in variational flows. The authors treat variational flows as dynamical systems and leverage shadowing theory to elucidate the abnormal behavior: common flows can exhibit a catastrophic accumulation of error, but surprisingly, results produced by flows are often accurate enough for applications despite the presence of serious numerical instability. ", "questions": "1. Does the shaowing property holds for continuous normalizing flow (empirically or theoretically) ?\n\n2. On theorem 5.1, the important constant $M$ depends on $x \\in \\mathcal{X}$, right?  If so, then we cannot derivative the shaowing property (which is a global property independent of $x$) from theorem 5.1.", "limitations": "The paper could benefit from a more detailed discussion of the limitations of the proposed approach.", "rating": 5, "confidence": 3, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "In terms of originality, the paper presents a new approach to understanding the behavior of variational flows in the presence of numerical instability. By treating variational flows as dynamical systems and leveraging shadowing theory, the authors provide theoretical guarantees on the error of sampling, density evaluation, and ELBO estimation. The paper also presents a diagnostic procedure to empirically validate results produced by numerically unstable flows in practice.\n\nIn terms of quality, the paper is well-written. The authors provide a clear motivation for their work and present their results in a coherent manner. The paper also includes several theorems, which adds to the quality of the work.\n\nIn terms of clarity, the paper is easy to follow and understand. The authors provide clear explanations of their methods and results, and the paper includes several figures and examples that help to illustrate the concepts presented.\n\nIn terms of significance, the paper makes mainly two contributions to the field of variational inference. The paper provides a theoretical framework to understand the behavior of variational flows in the presence of numerical instability, which is an important problem in the field. The paper also presents a diagnostic procedure to validate results produced by numerically unstable flows in practice.", "weaknesses": "One potential weakness of the paper \"Embracing the chaos: analysis and diagnosis of numerical instability in variational flows\" is that it focuses primarily on theoretical analysis and does not provide as much empirical validation of the proposed diagnostic procedure as would be desirable. While the authors do provide some empirical results to support their claims, it would be beneficial to see more extensive experiments that demonstrate the effectiveness of the proposed diagnostic procedure on more complex problems datasets like Cifar 10.\n\nAnother potential weakness of the paper is that the theoretical contribution with the concepts of dynamical systems and shadowing theory has an overlap with [1] , which is the main theoretical insight of the paper. Moreover, the theorems in ELBO esitimation and density estimation follows diretcly from the shadowing property.\n\nReference:\n1. Paul Tupper. The relation between approximation and shaowing in molecular dynamics. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for reviewing our manuscript! Please see our point-to-point respoonse as follows.\n> “One potential weakness of the paper \"Embracing the chaos: analysis and\ndiagnosis of numerical instability in variational flows\" is that it focuses\nprimarily on theoretical analysis and does not provide as much empirical\nvalidation of the proposed diagnostic procedure as would be desirable. While\nthe authors do provide some empirical results to support their claims, it\nwould be beneficial to see more extensive experiments that demonstrate the\neffectiveness of the proposed diagnostic procedure on more complex problems\ndatasets like Cifar 10.”\n\nPlease see our general response 4.\n\n> “Another potential weakness of the paper is that the theoretical contribution\n> with the concepts of dynamical systems and shadowing theory has an overlap\n> with [1], which is the main theoretical insight of the paper. Moreover, the\n> theorems in ELBO esitimation and density estimation follows directly from the\n> shadowing property.”\n\nPlease see our general response 1.\n\n> “Does the shadowing property holds for continuous normalizing flow (empirically or theoretically) ?”\n\nThis is a good question, and is also something we are exploring now!  Generally speaking, shadowing theory also applies to numerical discretization of continuous flows such ODE (see e.g. the Finite shadowing theorem in [Coomes95]).  And the downstream tasks of continuous NFs or neural ODEs, e.g., sampling and density evaluation, indeed rely on the numerical discretization of the underlying ODE. Hence, it should be possible to theoretically analyze the error for continuous NFs via shadowing, but we leave a deeper investigation for future work.\n\n[Coomes95] Brian A. Coomes, Hüseyin Koçak, and Kenneth J. Palmer. Rigorous computational shadowing of orbits of ordinary differential equations. Numerische Mathematik, 1995.\n\n\n> “On theorem 5.1, the important constant M\n depends on $x\\in \\mathcal{X}$, right? If so, then we cannot derive the shadowing property (which is a global property independent of $x$) from theorem 5.1.”\n\nThis is a very good point! Indeed, Theorem 5.1 involves a value of $M$ that is valid for a single numerical trajectory. Then as long as we know something about the smoothness of $F_n$ locally around that numerical trajectory (e.g. 3rd derivative bound), we can obtain an upper bound on $M$ for that trajectory. But you are right that in order to prove that shadowing holds globally, we would need to bound $M$ for all starting states $x \\in \\mathcal{X}$. We can obtain such a bound if we make a global smoothness assumption on $F_n$ — e.g., a uniformly bounded 3rd derivative would suffice.\n\nWe will certainly add a discussion of this in the camera ready. Thank you!\n\nOn this note, a very interesting potential way to relax this is to investigate “probabilistic shadowing,” where there is a shadowing trajectory only with high probability under the initial distribution. We will add this to our discussion!"}, {"review_id": "XawPlL56OE", "reviewer": "Reviewer_d5rn", "summary": "This paper investigated the impact of accumulated numerical error on the sampling, density evaluation, and evidence lower bound (ELBO) estimation in variational flows. It demonstrated that the results produced by flows are not destroyed by the serious numerical instability. To explain this phenomenon, it leveraged shadowing theory to theoretically bound the error of sampling, density evaluation, and ELBO estimation. It also developed a diagnostic procedure that can be used to validate results produced by numerically unstable flows in practice.", "questions": "1. What is the definition of \"s_k\" in theorem 4.1.\n2. What is the range of mathcal{X}, is it belong to R or R^N or R\\times[ 0,T]? \n3. I cannot follow the information in Figure 1(b)(d), e.g., what do you want to show through the dots in these figures? ", "limitations": "The limitations are not fully discussed in this paper, e.g., the theory only provide explanations while it does not inspire new algorithm for generative learning. It will narrow the impact of the theory that practitioners may ignore.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 1, "contribution": 3, "strengths": "This paper proposed an interesting question that why the accumulated numerical error does not influence the results of sampling, ELBO, etc. The technique is sound, and the analyses are solid. The derived theoretical results could support its claim.\n", "weaknesses": "The presentation could be improved, and the notations are hard to follow, see details below.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for reviewing our manuscript! Please see our point-to-point respoonse as follows.\n\n> What is the definition of $s_k$ in theorem 4.1.\n\n\"s\" stands for \"shadowing\". $(s_k)_{k = 0}^N$ denotes an exact orbit starting at $s = s_0$ that \"shadows\" the numerical trajectory (as mentioned in line 132--133). Please refer to our definitions and illustrations of orbits, found between lines 117--126 and in Figure 4. We will make sure to clarify this explicitly in the camera ready version!\n\n> What is the range of $\\mathcal{X}$, is it belong to $R$ or $R^N$ or $R\\times[ 0,T]$?\n\nAs outlined between lines 63–68, $\\mathcal{X} \\subseteq R^d$ represents the space on which the variational family $\\\\{q_\\lambda: \\lambda \\in \\Lambda\\\\}$ is supported.\n\n> I cannot follow the information in Figure 1(b)(d), e.g., what do you want to show through the dots in these figures?\n\nFigure 1(b)(d) serves to illustrate the impact of floating-point representation errors during flow computations. Ideally, in the absence of numerical errors, the computed orbit (blue scatters) would coincide perfectly with the exact orbit (red scatters). However, based on the figure, even for a relatively small $N$, there's a significant deviation of the numerical orbit from the exact one.\n\n> The limitations are not fully discussed in this paper, e.g., the theory only\nprovide explanations while it does not inspire new algorithm for generative\nlearning. It will narrow the impact of the theory that practitioners may\nignore.\n\nPlease see our general response 3 for a discussion of the limitations. Note that the purpose of this work is not to introduce a new methodology, but to provide a novel set of theoretical results that explain a counterintuitive, surprising phenomenon related to numerical errors in variational flows that we often see in practice.\n"}, {"review_id": "ydEkUnCgdh", "reviewer": "Reviewer_e3zg", "summary": "The authors present an analytical framework to quantify the effects of numerical errors on sampling, density and ELBO estimation in variational flows. The framework leverages the shadowing theory to derive error bounds based on the shadowing window sizes and the local Lipschitz constants of the dynamical operators. A procedure to verify whether the shadowing property holds in a given system is provided.", "questions": "* The analysis in section 4.3 is predicated on the basis that (ln 181) “when $F$ and $B$ are ergodic and measure-preserving for the target $\\pi$”, which is a pretty strong assumption. Do you have ways to establish this?\n* The horizontal axis in Figure 7 does not have the same range as those in Figures 4 to 6. Is this due to computational constraints, i.e. the proportionality to $N$? \n* It appears that $\\epsilon$ is still growing at the end while the sampling, density and ELBO errors all seem to have converged. Do you have any comments on the qualitative behaviors when extrapolating further? This would imply that the bounds eventually become very loose.\n", "limitations": "There is not much explicit discussion on the limitations. The computational cost for establishing and calculating the shadow window size is an obvious one.", "rating": 5, "confidence": 3, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "* [Significance] Numerical analysis is an important and under-studied aspect in flow-based models. This work takes a step forward in this direction and offers a plausible explanation for why fast accumulation of numerical errors do not necessarily lead to bad results.\n* [Quality] The derivations of the numerical bounds seem to make sense intuitively - although I did not check the math in detail.\n* [Clarity] The material is presented in a clear and logical manner. I find it easy to follow in general.\n* [Novelty] There exists prior work that is based on shadowing theory but this work extends the analysis to the computation of density and ELBO estimations.\n", "weaknesses": "* My primary concern about this work is mainly on the usefulness of the presented framework in practice. All results are predicated on the basis that the shadowing property is already established. The approach presented in section 5 looks cumbersome and does not seem to scale particularly well.\n* I do not find the numerical experiments to support the claims very well. They only demonstrate some qualitative traits - sampling, density and ELBO errors are much lower compared to trajectory errors, which coincides with the observation that the derived shadowing window sizes are relatively low. There is not much numerical evidence showing how tight the derived bounds are.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for reviewing our manuscript! Please see our point-to-point respoonse as follows.\n\n> My primary concern about this work is mainly on the usefulness of the\n> presented framework in practice. All results are predicated on the basis that\n> the shadowing property is already established. \n\nPlease see our general response 5 for a discussion of assumptions and the scope of our work.\n\n\n> The approach presented in section 5 looks cumbersome and does not seem to scale particularly well.\n\nYour concern about the complexity and scalability of the approach presented in Section 5 is valid. However, we wish to highlight that the primary focus of this work is to build a theoretical framework to understand the numerical stability of variational flows, and to explain the counterintuitive phenomenon regarding numerical error in a rigorous manner. The present approach in Section 5 should be considered as a diagnostic tool that enables users to verify when our theory can be applied. It's not positioned as an algorithm that would be run each time variational flows are employed. Even so, as discussed in our general response 2, it is likely that a much more efficient algorithm for computing $\\lambda$ exists.\n\n> \"I do not find the numerical experiments to support the claims very well. They\n> only demonstrate some qualitative traits - sampling, density and ELBO errors\n> are much lower compared to trajectory errors, which coincides with the\n> observation that the derived shadowing window sizes are relatively low. There\n> is not much numerical evidence showing how tight the derived bounds are\"\n\nPlease see our general response 4 for discussion of experiments.\n\n> “It appears that $\\epsilon$ is still growing at the end while the sampling, density\nand ELBO errors all seem to have converged. Do you have any comments on the\nqualitative behaviors when extrapolating further? This would imply that the\nbounds eventually become very loose.”\n\nThis is indeed an interesting phenomenon; as mentioned in our general response 3, understanding how $\\epsilon$ scales with $N$ in more generality is an important direction of future work. Indeed if $\\epsilon$ continues to grow, we expect that eventually the bounds should become loose. Empirically, it seems to grow linearly with a gentle slope, and we suspect the range of $N$ for which these bounds are reasonable is large. But we cannot make a rigorous claim about this at this point.\n\nNonetheless, we believe that our current results offer valuable insights and represent significant contributions to the community. To the best of our knowledge, our theory is the first among the variational flow literature that rigorously explains how numerical error influences the orbit evaluation and downstream tasks in different ways. \n\n> “The analysis in section 4.3 is predicated on the basis that (ln 181) “when \nF and B are ergodic and measure-preserving for the target”, which is a pretty\nstrong assumption. Do you have ways to establish this?”\n\nTo be clear, these assumptions are not required by all analysis in this section; they are just used to obtain stronger results when the dynamical system targets a particular distribution. Since ergodicity and pi-invariance are the two key conditions for convergence of MixFlows paper [5], and are also quite common in ergodic theory [Eisner15], we studied how shadowing would interact with those conditions in variational methods. Our intention was just to provide a new insight that a dynamical system targeting a particular distribution may behave even better numerically than a general system. Given that they are standard conditions, verifying them is beyond the scope of the present work.\n\nHowever, note that there are well-known measure-preserving systems often used in statistics (e.g. Hamiltonian flows, or an ODE satisfying the Liouville conditions [Neklyudov21]). Although one must discretize ODEs in practice, note that our numerically approximated map $\\hat F$ can involve discretization and does not need to be measure-preserving. On the other hand, ensuring ergodicity is challenging in practice, as acknowledged in Section 4.3 of [5]. We agree that verifying ergodicity is a current limitation of MixFlows, but addressing this issue is beyond the scope of our work.\n\n[Eisner15] Eisner, T., Farkas, B., Haase, M., and Nagel, R. Operator\nTheoretic Aspects of Ergodic Theory. Graduate Texts in\nMathematics. Springer, 2015.\n\n[Neklyudov21] Neklyudov, K., Bondesan, R., and Welling, M. Deterministic gibbs sampling via ordinary differential equations.\narXiv:2106.10188, 2021.\n \n\n\n> “The horizontal axis in Figure 7 does not have the same range as those in\n> Figures 4 to 6. Is this due to computational constraints, i.e. the\n> proportionality to N?”\n\nThis was just a minor oversight on our behalf. We originally chose the range of these plots to illustrate the typically linear (or slower) relationship between the shadowing window size and $N$. We can certainly make the range of $N$ match Figure 6 for the camera ready–thank you for pointing this out!\n\n\n> “There is not much explicit discussion on the limitations. The computational\ncost for establishing and calculating the shadow window size is an obvious\none.”\n\nPlease see general response 3 for a discussion of limitations, and general response 2 for a discussion of computational cost."}]

Author rebuttal:
# General response to the reviewers

We thank the reviewers for their valuable feedback. In this response, we address shared comments. Specific responses to each reviewer will follow separately.

## 1. Comparison to [23] (MECE, sFEg)

Our work is not an extension or competitor to the work of [23]; both focus on different objectives. While [23] establishes a relationship between shadowing and the weak distance between inexact and exact trajectories in dynamical systems, our work provides a theoretical basis for the counterintuitive behavior of numerically-implemented variational flows observed in practice. The overlap is minimal; our Proposition 4.2 echoes the main theorem in [23] with some differences (e.g., Levy-Prokhorov vs. bounded Lipschitz distance, our work does not place assumptions on the initial state distribution). We introduce a substantial set of new results in Theorems 4.3, 4.4, 4.6, 4.7, and Proposition 4.5 dedicated to density evaluation, ELBO estimation, and similar and strengthened results for MixFlows. These theorems are not consequences of  [23], employing distinct proof techniques. Our work offers significant contributions even for readers familiar with [23]. 

We'll elaborate on the relationship with [23] in the camera ready.

## 2. Computational cost of the shadowing window (c44C, e3zg). 

In our experiments, this was a fairly minor computational cost. For a flow length of N=500, `eigmin` took just a few seconds. Given that this was acceptable in practice, we didn’t code a specialized method. In the camera ready, we will report the computation time for estimating the shadowing window size.

Even so, we suspect that $\lambda$ can be computed more efficiently by utilizing the sparsity of $AA^T$: specifically, $AA^T$ is a symmetric positive definite block-tridiagonal matrix with bandwidth $d$ (line 476--477 in Appendix B) and so has $O(Nd^2)$ entries. We are currently investigating more efficient methods for calculating $\lambda$ along these lines, e.g., the inverse power method [Schatzman02, Chapter 13.3.3], or tridiagonalization via Lanczos iterations [Lanczos] followed by divide-and-conquer algorithms [Coakley13].

[Schatzman02] Michelle Schatzman. Numerical analysis: a mathematical introduction. Oxford University Press, USA, 2002.

[Lanczos] C. Lanczos: An Iteration Method for the Solution of the Eigenvalue Problem of Linear Differential and Integral Operators, J. Res. Nat. Bur. Stand. 49, 255 (1950)

[Coakley13] Ed S. Coakley and Vladimir Rokhlin. A fast divide-and-conquer algorithm for computing the spectra of real symmetric tridiagonal matrices. Applied and Computational Harmonic Analysis, 2013.



## 3. Discussion of the limitations (MECE, sFEg, d5rn, e3zg)

We agree there was not enough discussion of limitations in the original submission; given the extra page in the camera ready, we are happy to include a thorough discussion of the limitations of our work, based on the following key points:

- Our error bounds mainly address downstream tasks post-training; we do not provide results regarding how numerical instability influences the training process itself.
- While our theory is centered around error bounds that are proportional to the shadowing window size $\epsilon$, we’ve observed that $\epsilon$ can grow with $N$. This suggests further study on its theoretical growth rate; we provided an initial discussion of this in the Appendix, but more work needs to be done here.
- Our Theorem 5.1 verifies shadowing for a particular trajectory, but needs a global smoothness condition on $F_n$ for global shadowing. To help address this, our current analysis may be extended from the global shadowing property to a weaker version that only holds with high probability.  See our response to Reviewer sFEg.
- Our theoretical analysis is focused on standard variational flow methods, and is not adapted to recent architectures like continuous normalizing flows or neural ODEs.
- In our experiments, we employed a basic method to calculate the minimum eigenvalue of $AA^T$. Given its sparsity, a deeper exploration into more efficient techniques is merited.
- Finally, we could include a broader range of variational flows (e.g., RealNVP, GLOW, and Hamiltonian variational flow) in our experiments.

## 4. Lack of numerical experiments (MECE, sFEg, e3zg)

We agree that more experiments that cover other types of variational flow would be a nice addition to the paper. We will plan to include results for more types of flow in the camera ready.

However, we’d like to emphasize that more experimental results are not central to the contributions of this particular work. The key contributions of this work are theoretical: we provide results that use the shadowing property to explain why the error of variational flows in statistical applications remains controlled, despite the catastrophic error growth in their numerical trajectories (Figure 1(b)(d)). Our experiments are designed just to provide empirical illustration of the theory and investigate the typical shadowing window size.
 
## 5. Usefulness of the theory (shadowing is a strong assumption) (MECE, e3zg)

We respectfully disagree that finite shadowing is an overly strong assumption. Indeed, we believe this is precisely the right assumption to explain the counterintuitive error growth phenomena in many common dynamical systems, even though this assumption of course does not hold for all systems. Further, we provide a sufficient condition to verify shadowing in practice via Theorem 5.1, which makes the strength of the assumption a priori less of an issue.

As a side note: even beyond finite shadowing, the *infinite* shadowing property is generic for dynamical systems (informally, there is an infinite shadowing dynamical system “arbitrarily close to” any reasonable dynamical system; see submission Appendix C). This suggests finite shadowing, as a weaker property, should be similarly pervasive.

Meta-review: In this study, the authors delve into the repercussions of numerical instability on the reliability of various processes, including sampling, density evaluation, and evidence lower bound (ELBO) estimation, within the context of variational flows. The investigation begins with empirical examinations and shows that a commonly employed flows can encounter a detrimental accumulation of errors. This manifests in the form of significant deviations between the numerical flow map and its exact counterpart, leading to repercussions on sampling. Similarly, the numerical inverse flow map fails to accurately recover the initial input, thereby influencing the accuracy of density and ELBO computations. They develop theory to illuminate this behavior through theoretical guarantees concerning the errors in sampling, density evaluation, and ELBO estimation.

The reviewers thought the paper addresses an interesting and important problem, the experiment provide good evidence, and use of shadowing theory is interesting and the diagnostic procedure can be of value. They did have some concerns about the assumptions and novelty and overlap with some papers. The authors response seems to have alleviated these concerns and the reviews all recommend acceptance. I concur.
Decision: Accept (poster)