Year: 2021

Abstract:
Training and using modern neural-network based latent-variable generative models (like Variational Autoencoders) often require simultaneously training a generative direction along with an inferential (encoding) direction, which approximates the posterior distribution over the latent variables. Thus, the question arises: how complex does the inferential model need to be, in order to be able to accurately model the posterior distribution of a given generative model? In this paper, we identify an important property of the generative map impacting the required size of the encoder. We show that if the generative map is “strongly invertible" (in a sense we suitably formalize), the inferential model need not be much more complex. Conversely, we prove that there exist non-invertible generative maps, for which the encoding direction needs to be exponentially larger (under standard assumptions in computational complexity). Importantly, we do not require the generative model to be layerwise invertible, which a lot of the related literature assumes and is not satisfied by many architectures used in practice (e.g. convolution and pooling based networks). As a non-invertible generator results in a distribution supported on a low-dimensional manifold, this provides theoretical support for the empirical wisdom that learning deep generative models is harder when data lies on a sufficiently complex low-dimensional manifold.

heading: 1 INTRODUCTION
text: Many modern generative models of choice (e.g. Generative Adversarial Networks (Goodfellow et al., 2014), Variational Autoencoders (Kingma & Welling, 2013)) are modeled as non-linear, possibly stochastic transformations of a simple latent distribution (e.g. a standard Gaussian). A particularly common task is modeling the inferential (encoder) direction: that is, modeling the posterior distribution on the latents z given an observable sample x. Such a task is useful both at train time and at test time. At train time, fitting generative models like variational autoencoders via maximum likelihood often relies on variational methods, which require the joint training of a generative model (i.e. generator/decoder), as well as an inference model (i.e. encoder) which models the posterior distribution of the latent given the observables. At test time, the posterior distribution very often has some practical use, e.g. useful, potentially interpretable feature embeddings for data (Salimans et al., 2016; Berthelot et al., 2018), “intervening” on the latent space to change the sample in some targeted manner (Shen et al., 2020), etc. As such, the question of the “complexity" of the inference model (i.e. number of parameters to represent it using a neural network-based encoder) as a function of the “complexity" of the forward model is of paramount importance, so that training is computationally tractable:
Question: How complex does the inference (encoder) model need to be relative to the complexity of the generative (decoder) model ?
In this paper we identify an important property of the generative direction governing the complexity of the inference direction for variational autoencoders: bijectivity/invertibility of the mean of the generative direction. We prove that when the mean of the generative direction is invertible, the complexity of the inference direction is not much greater than the complexity of the generative direction. Conversely, when the mean of the generative direction is not invertible, modulo standard computational complexity conjectures from cryptography, we can exhibit instances where the inference direction has to be much more complex.
On the mathematical level, our techniques involve a neural simulation of a Langevin random walk to sample from the posterior of the latent variables and uncover novel connections between Langevin diffusions and (hierarchical) deep latent Gaussians. On the lower bound side, we provide a reduction from the existence of one-way Boolean permutations in computational complexity: that is, permutations that are easy to calculate, but hard to invert. We show that the existence of a small encoder for non-invertible generators would allow us to design an invertor for any Boolean permutation, thus violating the existence a one-way permutation. This is the first time such ideas have been applied to generative models.
Note that a non-invertible generator results in a distribution supported on a lower dimensional manifold. Thus, our results show that learning deep generative models is harder when data lies on a sufficiently complex low-dimensional manifold, corroborating similarly flavored empirical observations (Dai & Wipf, 2019; Arjovsky et al., 2017).

heading: 2 OUR RESULTS
text: The Variational Autoencoder (VAE) (Kingma & Welling, 2013) is one of the most commonly used paradigms in generative models. It’s trained by fitting a generator which maps latent variables z to observables x, denoted by pθ (x|z), as well as an encoder which maps the observables to the latent space, denoted by qφ (z|x). Here φ and θ are the encoder parameters and generator parameters respectively. Given n training samples {x(i)}ni=1, the VAE objective is given by
max φ,θ
1
n n∑ i=1 Ez∼qφ(.|x(i)) [ log pθ(x (i)|z) ] −KL ( qφ(z|x(i))||p(z) ) where p(z) is typically chosen to be a standard Gaussian. This loss can be viewed as a variational relaxation of the maximum likelihood objective, where the encoder qφ, in the limit of infinite representational power, is intended to model the posterior distribution pθ(z|x(i)). Setup: We will consider a setting in which the data distribution itself is given by some ground-truth generator G : Rdl → Rdo , and ask how complex (in terms of number of parameters) the encoder needs to be (as a function of the number of parameters of G), s.t. it approximates the posterior distribution p(z|x) of the generator. We will consider two standard probabilistic models for the generator/encoder respectively.
Definition 1 (Latent Gaussian). A latent Gaussian is the conditional distribution given by a stochastic pushforward of a Gaussian distribution. That is, for latent variable z ∈ Rdl and observable x ∈ Rdo , for a neural network G : Rdl → Rdo and noise parameter β2, we have p(x|z) = N (G(z), β2Ido) and p(z) = N (0, Idl).
In other words, a sample from this distribution can be generated by sampling independently z ∼ N (0, Idl) and ξ ∼ N (0, β2Ido) and outputting x = G(z) + ξ. This is a standard neural parametrization of a generator with (scaled) identity covariance matrix, a fairly common choice in practical implementations of VAEs (Kingma & Welling, 2013; Dai & Wipf, 2019).
We will also define a probabilistic model which is a composition of latent Gaussians (i.e. consists of multiple stochastic layers), which is also common, particularly when modeling encoders in VAEs, as they can model potentially non-Gaussian posteriors (Burda et al., 2015; Rezende et al., 2014):
Definition 2 (Deep Latent Gaussian). A deep latent Gaussian is the conditional distribution given by a sequence of stochastic pushforwards of a Gaussian distribution. That is, for observable z0 ∈ Rd0 and latent variables {zi ∈ Rdi}Li=1, for neural networks {Gi : Rdi−1 → Rdi}Li=1 and noise parameters {β2i }Li=1, the conditional distribution p(zL|z0) is a deep latent Gaussian when p(zi|zi−1) = N (Gi(zi−1), β2i Idi),∀i ∈ [L] and p(z0) = N (0, Id0).
In other words, a deep latent Gaussian is a distribution, which can be sampled by ancestral sampling, one layer at a time. Note that this class of distributions is convenient as a choice for an encoder in a VAE, since compositions are amenable to the reparametrization trick of Kingma & Welling (2013)—the randomness for each of the layers can be “presampled” and appropriately transformed (Burda et al., 2015; Rezende et al., 2014). Then, we ask the following:
Question: If a VAE generator is modeled as a latent Gaussian (that is, p(x|z) ≡ N (G(z), β2I)), s.t. the corresponding G has at most N parameters, and we wish to approximate the posterior p(z|x) by a deep latent Gaussian s.t. the total size of the networks in it have at most N ′ parameters, how large must N ′ be as function of N?
We will work in the setting dl = do = d, and prove a dichotomy based on the invertibility of G: namely, if G : Rd → Rd is bijective, and β ≤ O ( 1
d1.5 √ log d/
) , the poste-
rior p(z|x) can be -approximated in total variation distance by a deep latent Gaussian of size N ′ = O (N · poly(d, 1/β, 1/ )). Thus, if the neural network G is invertible, and for a fixed and a small-enough variance term β2, we can approximate the posterior with a deep latent Gaussian polynomially larger than G. On the other hand, if G is not bijective, if one-way-functions exist (a widely believed computational complexity conjecture), we will show there exists a VAE generator G : Rd → Rd of size polynomial in d, for which the posterior p(z|x) cannot be approximated in total variation distance for even an inverse polynomial fraction of inputs x, unless the inferential network is of size exponential in d.
Remark 1: Note that it is in fact sufficient to investigate only the case where dl = do = d. Let us elaborate why:
• For Theorem 1, a function G : Rdl → Rdo can be bijective, Lipschitz and strongly invertible (i.e. satisfies Assumptions 1 and 2) only if dl = do, so Theorem 1 (upper bound) would not make any sense unless the dimensions are equal. (Note: there are bijective maps between, say, R and R2, e.g. space-filling curves, but these cannot be Lipschitz, as Lipschitz maps preserve the Hausdorff dimension, and the Hausdorff dimensions of R and R2 are 1 and 2 respectively).
• Theorem 2 is a lower bound – meaning we aim to exhibit an example of a hard instance when G is not bijective. Having exhibited a hard instance for the dl = do = d case, it is trivially easy to construct a hard instance in the theorem where the output dimension is larger—which is the more common setting in practice. To do so, consider a circuit C̃ : {±1}dl → {±1}do , which equals to C, the one-way-circuit from Conjecture 1 on the first dl coordinates, and is equal to 1 (i.e. is constant) on the last do − dl coordinates. Then C̃ is one-way too—since the last do − dl values are just fixed to 1, if we can invert it, we can invert C too. The reduction in the proof of Theorem 2 then can just be performed with C̃ instead, giving a lower bound instance for the dl < do case.

heading: 2.1 UPPER BOUNDS FOR BIJECTIVE GENERATORS
text: We first lay out the assumptions on the mapG. The first is a quantitative characterization of bijectivity; and the second requires upper bounds on the derivatives of G up to order 3. We also have a centering assumption. We state these below. Assumption 1 (Strong invertibility). We will assume that the latent and observable spaces have the same dimension (denoted d), and G : Rd → Rd is bijective. Moreover, we will assume there exists a positive constant m > 0 such that ∀z1, z2 ∈ Rd, ‖G(z1)−G(z2)‖ ≥ m · ‖z1 − z2‖.
Remark 2: This is a stronger quantitative version of invertibility. Furthermore, the infinitesimal version of this condition (i.e. ‖z1−z2‖ → 0) implies that the smallest magnitude of the singular values of the Jacobian at any point is lower bounded by m, that is ∀z ∈ Rd, min
i∈[d] |σi(JG(z))| ≥ m > 0.
Since m is strictly positive, this in particular means that the Jacobian is full rank everywhere.
Remark 3: Note, we do not require that G is layerwise invertible (i.e. that the each map from one layer to the next is invertible) – if that is the case, at least in the limit β → 0, the existence of an inference decoder of comparable size toG is rather obvious: we simply invert each layer one at a time. This is important, as many architectures based on convolutions perform operations which increase the dimension (i.e. map from a lower to a higher dimensional space), followed by pooling (which decrease the dimension). Nevertheless, it has been observed that these architectures are invertible in practice—Lipton & Tripathi (2017) manage to get almost 100% success at inverting an off-the-shelf trained model—thus justifying this assumption. Assumption 2 (Smoothness). There exists a finite positive constant M > 0 such that :
∀z1, z2 ∈ Rd, ‖G(z1)−G(z2)‖ ≤M · ‖z1 − z2‖
Moreover, we will assume that G has continuous partial derivatives up to order 3 at every z ∈ Rd and the derivatives are bounded by finite positive constants M2 and M3 as
∀z ∈ Rd, ∥∥∇2G(z)∥∥
op ≤M2 <∞, ∥∥∇3G(z)∥∥ op ≤M3 <∞
Remark 4: This is a mild assumption, stating that the map G is smooth to third order. The infinitesimal version of this means that the largest magnitude of the singular values of the Jacobian at any point is upper bounded by M , that is ∀z ∈ Rd, max
i∈[d] |σi(JG(z))| = ‖JG(z)‖op ≤M <∞.
Remark 5: A neural network with activation function σ will satisfy this assumption when σ : R→ R is Lipschitz, and maxa |σ′(a)| & maxa |σ′′(a)| are finite. Assumption 3 (Centering). The map G : Rd → Rd satisfies G(0) = 0.
Remark 6: This assumption is for convenience of stating the bounds—we effectively need the “range” of majority of the samples x under the distribution of the generator. All the results can be easily restated by including a dependence on ‖G(0)‖. Our main result is then stated below. Throughout, the O(.) notation hides dependence on the map constants, namely m,M,M2,M3. We will denote by dTV(p, q) the total variation distance between the distributions p, q. Theorem 1 (Main, invertible generator). Consider a VAE generator given by a latent Gaussian with µ = 0, Σ = I , noise parameter β2 and generator G : Rd → Rd satisfying Assumptions 1 and 2, which has N parameters and a differentiable activation function σ. Then, for
β ≤ O  1 d1.5 √ log d  (1) there exists a deep latent Gaussian withN ′ = O ( N · poly(d, 1β , 1 ) ) parameters and activation functions {σ, σ′, ρ} and a neural network φ with O(N ′) parameters and activation functions {σ, σ′, ρ}, where ρ(x) = x2, such that with probability 1−exp(−O(d)) over a sample x from the VAE generator, φ(x) produces values for the parameters of the deep latent Gaussian such that the distribution q(z|x) it specifies satisfies dTV (q(z|x), p(z|x)) ≤ .
Remark 7: Having an encoder q(z|x) with parameters produced by a neural network taking x as input is fairly standard in practice—it is known as amortized variational inference.(Mnih & Gregor, 2014; Kingma & Welling, 2013; Rezende et al., 2014)
Remark 8: The addition of ρ in the activation functions is for convenience of stating the bound. Using usual techniques in universal approximation it can be simulated using any other smooth activation, see Appendix G.

heading: 2.2 LOWER BOUNDS FOR NON-BIJECTIVE GENERATORS
text: We now discuss the case when the generative map G is not bijective, showing an instance such that no small encoder corresponding to the posterior exists. The lower bound will be based on a reduction from the existence of one-way functions—a standard complexity assumption in theoretical computer science (more concretely, cryptography). Precisely, we will start with the following form of the one-way-function conjecture: Conjecture 1 (Existence of one-way permutations (Katz & Lindell, 2020)). There exists a bijection f : {−1, 1}d → {−1, 1}d computable by a Boolean circuit C : {−1, 1}d → {−1, 1}d of size poly(d), but for every T (d) = poly(d) and (d) = 1poly(d) and circuit C
′ : {−1, 1}d → {−1, 1}d of size T (d) it holds Prz∼{±1}d [C′(C(z)) = z] ≤ (d).
In other words, there is a circuit of size polynomial in the input, s.t. for every polynomially sized invertor circuit (the two polynomials need not be the same—the invertor can be much larger, so long as it’s polynomial), the invertor circuit succeeds on at most an inverse polynomial fraction of the inputs. Assuming this Conjecture, we show that there exist generators that do not have small encoders that accurately represent the posterior for most points x. Namely:
Theorem 2 (Main, non-invertible generator). If Conjecture 1 holds, there exists a VAE generator G : Rd → Rd with size poly(d) and activation functions {sgn,min,max}, s.t., for every β = o(1/ √ d), every T (d) = poly(d) and every (d) = 1/poly(d), any encoder E that can be represented by a deep latent Gaussian with networks that have total number of parameters bounded by T (d), weights bounded by W , activation functions that are L-Lipschitz, and node outputs bounded by M with probability 1− exp(−d) over a sample x from G and L,M,W = o(exp(poly(d))), we have:
Prx∼G [ dTV(E(z|x), p(z|x)) ≤ 1
10
] ≤ (d)
Thus, we show the existence of a generator for which no encoder of polynomial size reasonably approximates the posterior for even an inverse-polynomial fraction of the samples x (under the distribution of the generator).
Remark 9: The generator G, though mapping from Rd → Rd will be highly non-invertible. Perhaps counterintuitively, Conjecture 1 applies to bijections—though, the point will be that G will be simulating a Boolean circuit, and in the process will give the same output on many inputs (more precisely, it will only depend on the sign of the inputs, rather than their values).
Remark 10: The choice of activation functions {sgn,min,max} is for convenience of stating the theorem. Using standard universal approximation results, similar results can be stated with other activation functions, see Appendix G.
Remark 11: The restrictions on the Lipschitzness of the activations, bounds of the weights and node outputs of E are extremely mild—as they are allowed to be potentially exponential in d—considering that even writing down a natural number in binary requires logarithmic number of digits.

heading: 3 RELATED WORK
text: On the empirical side, the impact of impoverished variational posteriors in VAEs (in particular, modeling the encoder as a Gaussian) has long been conjectured as one of the (several) reasons for the fuzzy nature of samples in trained VAEs. Zhao et al. (2017) provide recent evidence towards this conjecture. Invertibility of generative models in general (VAEs, GANs and normalizing flows), both as it relates to the hardness of fitting the model, and as it relates to the usefulness of having an invertible model, has been studied quite a bit: Lipton & Tripathi (2017) show that for off-theshelf trained GANs, they can invert them with near-100% success rate, despite the model not being encouraged to be invertible during training; Dai & Wipf (2019) propose an alternate training algorithm for VAEs that tries to remedy algorithmic problems during training VAEs when data lies on a lowerdimensional manifold; Behrmann et al. (2020) show that trained normalizing flows, while being by design invertible, are just barely so — the learned models are extremely close to being singular.
On the theoretical side, the most closely relevant work is Lei et al. (2019). They provide an algorithm for inverting GAN generators with random weights and expanding layers. Their algorithm is layerwise — that is to say, each of the layers in their networks is invertible, and they invert the layers one at a time. This is distinctly not satisfied by architectures used in practice, which expand and shrink — a typical example are convolutional architectures based on convolutions and pooling. The same paper also shows NP-hardness of inverting a general GAN, but crucially they assume the network G is part of the input (their proof does not work otherwise). Our lower bound can be viewed as a “non-uniform complexity” (i.e. circuit complexity) analogue of this, since we are looking for a small neural network E, as opposed to an efficient algorithm; crucially, however G is not part of the input (i.e. G can be preprocessed for an unlimited amount of time). Hand & Voroninski (2018) provide similar guarantees for inverting GAN generators with random weights that satisfy layerwise invertibility, albeit via non-convex optimization of a certain objective.

heading: 4 PROOF OVERVIEW
text: Overview of Theorem 1 Recall, given a generator G and noise level β2, we wish to provide a deep latent Gaussian (Definition 2) that approximates a sample from the distribution p(z|x), with high
probability over the draws of x from the underlying density. To show Theorem 1, there are three main ingredients:
1. We show that gradient descent on the function f(z) = 12‖G(z)− x‖ 2, run for O(d/β2) number
of steps, recovers a point zinit close to argmaxz log p(z|x) – i.e. the mode of the posterior.
2. We show that Langevin diffusion, a Markov process that has p(z|x) as it’s stationary distribution, started close to the mode of the posterior via zinit computed by GD above, run for O(log(d/ )) amount of time; returns a sample from a distribution close to p(z|x) in total variation distance.
3. We show that the above two steps can be “simulated" by a deep latent Gaussian of size at most poly(d, 1β , 1 ) larger than the original network, by inductively showing that each next iterate in the
gradient descent and Markov chain can be represented by a neural network of size comparable to the generator G.
Part 1 is formalized in Lemma 1. The main observation is that the gradient of the loss will have the form ∇f(z) = JG(z)T (G(z) − x). By Strong Invertibility (Assumption 1), JG(z) will be an invertible matrix—therefore, any stationary points will necessarily be the (unique) point z s.t. G(z) = x. This in turn, for small β, will be close to argmaxz log p(z|x). Part 2 is formalized in Lemma 2, and is technically most involved. As β → 0, the distribution p(z|x) concentrates more around it’s mode, and approaches a Gaussian centered at it. More precisely, we show that starting at zinit, there is a region D such that the probability of Langevin diffusion leaving this region is small (Lemma 13) and the distribution restricted to D is log-concave (Lemma 4). As a consequence of the latter claim, reflected Langevin diffusion (a “restricted” version of the standard Langevin diffusion) converges rapidly to p(z|x) restricted to D. As a consequence of the former claim, the reflected and standard Langevin diffusion are close in the total variation sense.
Finally, Part 3 is based on the observation that we can approximately simulate gradient descent, as well as the Langevin dynamics using a deep latent Gaussian (Definition 2). Proceeding to gradient descent, the key observation is that if we give a neural network z as the input, the point z − ηJG(z)T (G(z)− x) can be calculated by a network of size comparable to the size of a neural network representing JG; the Jacobian of a neural network. And in turn, the Jacobian network can be represented by a neural network at most a factor O(d) bigger than the neural network representing G—this is by noting that each partial derivative of G can be represented by a neural network of size comparable to G, essentially by the backpropagation algorithm (Lemma 5).
The idea to simulate Langevin dynamics is similar: the Euler-Mariyama discretization of the corresponding stochastic differential equation (with precision h) is a discrete-time Markov chain, s.t. the updates take the form
z ← z − h∇L(z) + √ 2hξt, ξ ∼ N(0, I) (2)
where L the unnormalized log-likelihood of the posterior p(z|x) (equation 3 below). Same as is the case for gradient descent, a sample following the distribution of these updates can be implemented by a latent Gaussian, with a neural net of size at most O(d) bigger than the size of G—the randomness being handled by Gaussian noise that the latent Gaussian has its disposal.
Overview of Theorem 2 The proof of Theorem 2 proceeds by a reduction to Conjecture 1. We consider G(z) = C(sgn(z)), where C is a circuit satisfying the assumptions of Conjecture 1. We show that we can use an encoder q(z|x) that works with good probability over x (wrt to the distribution of G) to construct a circuit C′ that violates Conjecture 1. The main idea is as follows: the samples x will, with high probability, lie near points on the hypercube {±1}d. For such points, the posterior p(z|x) will concentrate around C−1(x̃), so taking the sign of a sample from q(z|x) should, with good probability, produce an inverse of sgn(x).
Simulating this construction using a Boolean circuit, we can get a “randomized” circuit that inverts C. Producing a deterministic circuit follows ideas from complexity theory, where we show that given a randomized circuit, there exists a single string of randomness that works for all inputs, if the success probability of the algorithm is high enough.

heading: 5 MAIN INGREDIENTS: THEOREM 1
text: First, we introduce a few pieces of notation that will be used in multiple places in the proofs. We will also omit subscripting quantities with x to ease notation whenever x is clear from context. For notational convenience, we will denote the unnormalized log-likelihood of p(z|x) as L : Rd → R, s.t.
L(z) := 1
2
( ‖z‖2 + 1
β2 ‖G(z)− x‖2
) (3)
Note, by Bayes rule, p(z|x) ∝ p(x|z) p(z), so p(z|x) ∝ e−L(z). Note this is a well-defined probability distribution since ∀z ∈ Rd, L(z) ≥ 12‖z‖
2, therefore the partition function is finite. Further, we will denote the inverse of x as ẑ, namely:
ẑ := arg min z∈Rd
‖G(z)− x‖2 (4)
Note since G is invertible, ẑ is uniquely defined, and it’s such that G (ẑ) = x. Crucially due to this, the function L(z) is strongly convex in some region around the point ẑ – which can be seen by evaluating∇2L(ẑ) from equation 3. These properties are convenient – that is why we base our analysis on this inverse point.

heading: 5.1 PART 1: CONVERGENCE OF GRADIENT DESCENT
text: First, we show that gradient descent on the function f(z) = 12‖G(z) − x‖ 2, started at the origin, converges to a point close to ẑ in a bounded number of steps. Precisely: Lemma 1 (Convergence of Gradient Descent to the Inverse). When G satisfies Assumptions 1 and 2, with probability 1 − exp(−O(d)) over the choice of x, running gradient descent on f(z) = 1 2‖G(z) − x‖
2 with the starting point z(0) := 0 and learning rate η := 1Q where Q = O(d), for S = O(d 2
δ2 ) steps converges to a point z (S), such that ‖z(S) − ẑ‖ ≤ δ.
As mentioned in Section 4, the idea is that while f is non-convex, its unique stationary point is ẑ, as ∇f(z) = JG(z)T (G(z)− x) and JG(z) is an invertible matrix ∀z ∈ Rd from Assumption 1. The proof essentially proceeds by a slight tweak on the usual descent lemma (Nesterov, 2003), due to the fact that the Hessian’s norm is not bounded in all of Rd and is given in Appendix B.

heading: 5.2 PART 2: FAST MIXING OF LANGEVIN DIFFUSION
text: We first review some basic definitions and notation. Recall, Langevin diffusion is described by a stochastic differential equation (SDE): Definition 3 (Langevin Diffusion). Langevin Diffusion is a continuous-time Markov process {zt : t ∈ R+} defined by the stochastic differential equation
dzt = −∇g(zt)dt+ √
2dBt t ≥ 0 (5) where g : Rd → R, and {Bt : t ∈ R+} is a Brownian motion in Rd with covariance matrix Id. Under mild regularity conditions, the stationary distribution of this process is P : Rd → R+ such that P (z) ∝ e−g(z).
For the reader unfamiliar with SDEs, the above process can be interpreted as the limit as h→ 0 of the Markov chain
z ← z − h∇L(z)dt+ √
2hξ, ξ ∼ N(0, I) (6) The convergence time of the chain to stationarity is closely tied to the convexity properties of L, through functional-analytic tools like Poincaré inequalities and Log-Sobolev inequalities (Bakry & Émery, 1985). We show that the distribution p(z|x) can be sampled -approximately in total variation distance by running Langevin dynamics for O(log(d/ )) number of steps. Namely:
Lemma 2 (Sampling Complexity for general G). Let β = O (
1 d1.5 √ log d
) . Then, with probability
1 − exp(−O(d)) over the choice of x, if we initialize Langevin diffusion equation 6 at z0 which satisfies ‖z0−ẑ‖ ≤ Ω ( β · √ d log d ) , for some T = O( 1d2 log 1 ) we have dTV (PT , p(z|x)) ≤ /2
where Pt is the density corresponding to zt.
Note, the lemma requires the variance of the latent Gaussian β to be sufficiently small. This intuitively ensures the posterior p(z|x) is somewhat concentrated around its mode — and we will show that when this happens, in some region around the mode, the posterior is close to a Gaussian. The proof of this has two main ingredients: (1) zt stays in a region D around ẑ with high probability; (2) The distribution p(z|x), when restricted to D is log-concave, so a “projected” version of Langevin diffusion (formally called reflected Langevin diffusion, see Definition 4) mixes fast to this restricted version of p(z|x). The first part is formalized as Lemma 13 in Appendix C. The main mathematical tool is a stochastic differential equation characterizing the change of the distance from ẑ. Namely, we show: Lemma 3 (Change of distance from ẑ). Let η(z) := 12‖z − ẑ‖
2. If z follows the Langevin SDE corresponding to zt, then for z ∈ D as defined in Lemma 13 and β ≤ β0 as defined above, it holds that
dη(z) ≤ −m 2
β2 · η(z)dt+
( ‖ẑ‖2 + d ) dt+ 2 √ η(z)dBt
This SDE is analyzed through the theory of Cox-Ingersoll-Ross (Cox et al., 2005) processes, which are SDEs of the type dMt = (a−bMt)dt+c √ MtdBt, where a, b, c are positive constants. Intuitively, this SDE includes an ‘attraction term’ −bMtdt which becomes stronger as we increase Mt (in our case, as we move away from the center of the region ẑ), and a diffusion term c √ MtdBt due to the injection of noise. The term (‖ẑ‖2 + d)dt comes from the heuristic Itó calculus equality dB2t = dt. The “attraction term” −m 2
β2 η(z)dt comes from the fact that L(z) is approximately strongly convex near ẑ. (Intuitively, this means the Langevin process is attracted towards ẑ.)
To handle the second part, we show the distribution p(z|x) restricted to D is log-concave: Lemma 4 (Strong convexity of L over D). For all z ∈ D, ∇2L(z) I .
Since the set D is a l2 ball (which is of course convex), we can show an alternate Markov process which always remains inside D, called the reflected Langevin diffusion mixes fast. Definition 4 (Reflected Langevin diffusion, (Lions & Sznitman, 1984; Saisho, 1987)). For a sufficiently regular region D ⊆ Rd, there exists a measure U supported on ∂D, such that the continuoustime diffusion process {yt : t ∈ R+} starting with y0 ∈ D defined by the stochastic differential equation:
dyt = −∇g(yt)dt+ √
2dBt + νtU(yt)dt t ≥ 0 (7) has as stationary measure P̃ : D → R+ such that P̃ (y) ∝ e−g(y). Here g : Rd → R, {Bt : t ∈ R+} is a Brownian motion in Rd with covariance matrix Id, and νt is an outer normal unit vector to D.
For readers unfamiliar with SDEs, the above process can be interpreted as the limit as h→ 0 of ∀k ≥ 0, yk+1 = ΠD ( yk − h∇g(yk) + √ 2hξk ) ξk ∼ N (0, Id) (8)
where ΠD is the projection onto the set D. For total variationa distance, the mixing time in total variation distance of this diffusion process was characterized by Bubeck et al. (2018): Theorem 3 (Mixing for log-concave distributions over convex set, Bubeck et al. (2018): Proposition 2). Suppose a measure p : Rd → R+ of the form p(x) ∝ e−g(x) is supported over S ⊆ Rd of diameter R which is convex, and ∀x ∈ S,∇2g & 0. Furthermore, let x0 ∈ S. Then, if pt is the distribution after running equation 7 for time t starting with a Dirac distribution at x0, we have
TV(pt, p) . e − t 2R2
Since the set D in Lemma 13 is an l2 ball (which is of course convex), it suffices to show that L is convex on D—which is a straightforward calculation. The details of this are again in Appendix C.

heading: 5.3 PART 3: SIMULATION BY NEURAL NETWORKS
text: Finally, we show that the operations needed for gradient descent as well as Langevin sampling can be implemented by a deep latent Gaussian. Turning the gradient descent first, we need to “simulate” by a
neural net the operation z ← z− ηJG(z)T (G(z)− x) roughlyO(d/β2) number of times. We prove that there is a neural network that given an input z, outputs z− ηJG(z)T (G(z)− x). The main ideas for this are that: (1) If G is a neural network with N parameters, each coordinate of JG(z) can be represented as a neural network with O(N) parameters. (2) If f, g are two neural networks with N and M parameters respectively, f + g and fg can be written as neural networks with O(M + N) parameters.
The first part essentially follows by implementing backpropagation with a neural network: Lemma 5 (Backpropagation Lemma, Rumelhart et al. (1986)). Given a neural network G : Rd → R of depth l and size N , there is a neural network of size O(N + l) which calculates the gradient ∂G/∂zi for i ∈ [d].
For the second part, addition is trivially implemented by “joining” the two neural networks and adding one extra node to calculate the sum of their output. For multiplication, the claim follows by noting that xy = 14 (x+ y) 2 − 14 (x− y) 2. Using the square activation function, this can be easily computed by a 2-layer network.
Proceeding to the part about simulating Langevin diffusion, we first proceed to create a discrete-time random process that closely tracks the diffusion. The intuition for this was already mentioned in Section 4— the SDE describing Langevin can be viewed as the limit h→ 0 of
z ← z − h [ z + 1
β2 JG(z)
T (G(z)− x) ] + √ 2hξ, ξ ∼ N(0, I)
It can be readily seen that the drift part, that is 1β2 JG(z) T (G(z)− x) can be calculated by a neural net via the same token as in the case of the gradient updates. Thus, the only difference is in the injection of Gaussian noise. However, each layer Gi in a deep latent Gaussian exactly takes a Gaussian sample ξ and outputs Gi(z) + γξ – so each step of the above discretization of Langevin can be exactly implement by a layer in a deep latent Gaussian. The details are in Appendix D.

heading: 6 MAIN INGREDIENTS: THEOREM 2
text: We prove Theorem 2 proceeds by a reduction to Conjecture 1. Namely, if C is a Boolean circuit satisfying the one-way-function assumptions in Conjecture 1, we considerG(z) = C(sgn(z)). Suppose there is a deep latent Gaussian of size T (n) with corresponding distribution q(z|x), s.t. with probability at least (d) over the choice of x (wrt to the distribution of G), it satisfies dTV(p(z|x), q(z|x)) ≤ 1/10. We will how to construct from the encoder a circuit C′ that violates Conjecture 1. Namely, consider a x̃ constructed by sampling z̃ ∼ {±1}d uniformly at random, and outputting x̃ = C(z̃). To violate Conjecture 1, we’d like to find C−1(x̃), using the encoder E. We first note that x = x̃ + ξ, ξ ∼ N(0, β2Id) is distributed according to the distribution of the generator G. Moreover, we can show that with probability at least (d)/2, it will be the case that x is s.t. both dTV(p(z|x), q(z|x)) ≤ 1/10 and x is very close to a point on the discrete hypercube {±1}d. The reason for this is that the samples x themselves for small enough β with high probability are close to a point on the hypercube. For such points, it will be the case that q(z|x) is highly concentrated around C−1(x̃), thus returning the sign of a sample will with high probability invert x̃. By simulating this process using a Boolean circuit, we get a “randomized” circuit (that is, a circuit using randomness) that succeeds at inverting C with inverse polynomial success. To finally produce a deterministic Boolean circuit that violates Conjecture 1, we use a standard trick in complexity theory: boosting the success probability by generating multiple samples together with the assumption on the weights and Lipschitzness of E. For details, refer to Appendix E.

heading: 7 CONCLUSION
text: In this paper we initiated the first formal study of the representational complexity of the good encoders in VAEs. We proved the following dichotomy: generators with invertible means give rise to posteriors that can be approximated by an encoders not much larger than the generator. On the other hand, for non-invertible generators, the corresponding encoder may need to be exponentially larger. We hope our work will stimulate research of other data distribution properties that render inference tractable.

heading: A TECHNICAL DEFINITIONS
text: In this section, we state several definitions and standard lemmas, to cover some technical prerequisites for the proof.
We first recall a few definitions and lemmas about two notions of distance between probabilities we will use extensively, total variation and χ2 distance. Definition 5 (Total variation distance). The total variation distance between two distributions P,Q is defined as
dTV(P,Q) := sup A |P (A)−Q(A)|
Definition 6 (χ2 distance). The χ2 distance between two distributions P,Q, s.t. P is absolutely continuous with respect to Q, is defined as
χ2(P,Q) := EQ ( P Q − 1 )2
The following two lemmas about the total variation and χ2 distance are standard: Lemma 6 (Coupling Lemma). Let P,Q : Ω → R be two distributions, and c : Ω⊗2 → R be any coupling of P,Q. Then, if (X,X ′) are random variables following the distribution c, we have:
dTV(P,Q) ≤ Pr[X 6= X ′]
Lemma 7 (Inequality between TV and χ2). Let P,Q be probability measures, s.t. P is absolutely continuous with respect to Q. We then have:
dTV(P,Q) ≤ 1
2
√ χ2(P,Q)
We will be heavily using notions from continuous-time Markov Processes. In particular: Definition 7 (Markov semigroup). We say that a family of functions {Pt(x, y)}t≥0 on a state space Ω is a Markov semigroup if Pt(x, ·) is a distribution on Ω and
Prt+s(x, y) = ∫ Ω Pt(x, z)Ps(z, y)dz
for all x, y ∈ Ω and s, t ≥ 0. Definition 8 (Continuous time Markov processes). A continuous time Markov process (Xt)t≥0 on state space Ω is defined by a Markov semigroup {Pt(x, y)}t≥0 as follows. For any measurable A ⊆ Ω
Pr(Xs+t ∈ A) = ∫ A Pt(x, y)dy := Pt(x,A)
Moreover Pt can be thought of as acting on a function g as (Ptg)(x) = EPt(x,·)[g(y)] = ∫
Ω
g(y)Pt(x, y)dy
Finally we say that p(x) is a stationary distribution if X0 ∼ p implies that Xt ∼ p for all t.
With this in place, we can finally define the Langevin diffusion Markov process: Definition 9 (Langevin Diffusion, Definition 3 restated). Langevin Diffusion is a continuous-time Markov process {zt : t ∈ R+} defined by the stochastic differential equation
dzt = −∇g(zt)dt+ √ 2dBt t ≥ 0 (9)
where g : Rd → R, and {Bt : t ∈ R+} is a Brownian motion in Rd with covariance matrix Id. Under mild regularity conditions, the stationary distribution of this process is P : Rd → R+ such that P (z) ∝ e−g(z).
It will also be useful to have a definition of a diffusion process that only stays within some support D, and thus has a stationary distribution supported on D. Specifically:
Definition 10 (Reflected Langevin diffusion, (Lions & Sznitman, 1984; Saisho, 1987), Definition 4 restated). For a sufficiently regular region D ⊆ Rd, there exists a measure U supported on ∂D, such that the continuous-time diffusion process {yt : t ∈ R+} starting with y0 ∈ D defined by the stochastic differential equation:
dyt = −∇g(yt)dt+ √ 2dBt + νtU(yt)dt t ≥ 0 (10)
has as stationary measure P̃ : D → R+ such that P̃ (y) ∝ e−g(y). Here g : Rd → R, {Bt : t ∈ R+} is a Brownian motion in Rd with covariance matrix Id, and νt is an outer normal unit vector to D.
Towards introducing quantities governing the mixing time of a continuous-time Markov process, we introduce Definition 11. The generator L of the Markov Process is defined (for appropriately restricted functionals g) as
Lg = lim t→0 Ptg − g t .
Moreover if p is the unique stationary distribution, the corresponding Dirichlet form and variance are
EM (g, h) = −Ep〈g,Lh〉 and Varp(g) = Ep(g − Epg)2
respectively. We will use the shorthand E(g) := E(g, g).
A quantity of particular interest when analyzing the mixing time of Markov processes is the Poincaré constant: Definition 12 (Poincaré Inequality). A continuous Markov process satisfies a Poincaré inequality with constant C if for all functions g such that EM (g) is defined (finite), it holds that:
EM (g) ≥ 1
C · VarP (g)
where P (in the subscript) is the stationary distribution of the Markov process. We will abuse notation, and for such a Markov process with stationary distribution P , denote by Cpc(P ) the smallest C such that above Poincaré inequality is satisfied; and call it the Poincaré constant of P .
The relationship between the Poincaré constant and the mixing time of a random walk (in particular, the one in Definition 4) is given by the following (standard) lemma: Lemma 8 (Mixing in χ2 from Poincaré). Let z̃t follow the stochastic differential equation of reflected Langevin diffusion (Definition 4) with stationary measure P̃ . Let P̃t be the density of z̃t and let P̃0 be absolutely continuous with respect to the Lebesgue measure. Then:
(1) If P̃0 is supported on D, P̃t is supported on D, ∀t > 0
(2) χ2(P̃t, P̃ ) ≤ e−t/Cpcχ2(P̃0, P̃ )
Finally, we will also need a discretized version of the continuous-time Langevin diffusion. Definition 13 (Euler-Mariyama Discretization for the Langevin Diffusion). For a potential g : Rd → R, the Euler-Mariyama discretization of the continuous Langevin diffusion with discretization level h is the discrete-time random walk
∀k ≥ 0, zk+1 = zk − h∇g(zk) + √
2hξk ξk ∼ N (0, Id) (11) We will also associate with this discrete-time process a stochastic differential equation
dẑt = −∇g(ẑbt/hch)dt+ √ 2dBt (12)
which has the property that ẑkh and zk follow the same distribution for all k ≥ 0, if z0 and ẑ0 follow the same distribution.
Note the equivalence between equation 11 and equation 12 follows by noting that the drift in the interval [th, (t+ 1)h] in equation 12 is constant – so we only need to integrate the Brownian motion for time h.
We can also consider the discretized version of the reflected Langevin diffusion, which looks like projected gradient descent with Gaussian noise Bubeck et al. (2018).
Definition 14 (Discretization for the Reflected Langevin Diffusion). For a potential g : Rd → R, the Euler-Mariyama discretization of the reflected Langevin diffusion with discretization level h is the discrete-time random walk
∀k ≥ 0, yk+1 = ΠD ( yk − h∇g(yk) + √ 2hξk ) ξk ∼ N (0, Id) (13)
where ΠD is the projection onto the set D. We can also associate with this discrete-time process a stochastic differential equation
dŷt = −∇g(ŷbt/hch)dt+ √ 2dBt + νtU(ŷt)dt (14)
for a measure U supported on ∂D, which has the property that ŷkh and yk follow the same distribution for all k ≥ 0, if y0 and ŷ0 follow the same distribution.
B PROOF OF LEMMA 1: CONVERGENCE OF GRADIENT DESCENT TO ẑ
First, we state the tweak on the usual descent lemma (Nesterov, 2003), due to the fact that the Hessian’s norm is not bounded in all of Rd. We only need to identify a region A containing the starting point z(0) such that the function’s value at points outside A are guaranteed to be more than the function value at z(0). Namely: Lemma 9 (Adjusted Descent Lemma). Let f : Rd → R be a twice differentiable function. For a regionA, let it hold that ∀z ∈ A, ‖∇2f(z)‖op ≤ Q <∞. Further, let ∀z ∈ Rd\A, f(z) ≥ f(z(0)) for a known point z(0) ∈ A. Then the gradient descent iterates starting with z(0) with the learning rate η := 1Q satisfy
f(z(s+1)) ≤ f(z(s))− 1 2Q ‖∇f(z(s))‖2 s ≥ 0
We first state the well-known descent lemma. Lemma 10 (Descent Lemma, Eq 1.2.13 in (Nesterov, 2003)). Let f be a twice differentiable function, and at a point z, the hessian satisfies ‖∇2f(z)‖op ≤ Q <∞. Then with the learning rate η := 1Q , the next gradient descent iterate z+ satisfies
f(z+) ≤ f(z)− 1 2Q ‖∇f(z)‖2
Using this, we prove the Adjusted Descent Lemma.
Proof of Lemma 9. We will prove the Lemma by induction. To that end, let us call by I(s),∀s ≥ 0 the proposition in the Lemma statement, i.e.
I(s) : f(z(s+1)) ≤ f(z(s))− 1 2Q ‖∇f(z(s))‖2
We first write two more propositions J(s) and K(s):
J(s) : z(s) ∈ A K(s) : f(z(s)) ≤ f(z) ∀z ∈ Rd \ A
We will use induction to show that I(s) ∧ J(s) ∧K(s) hold true ∀s ≥ 0. J(0) holds true, since z(0) ∈ A is given in the statement. K(0) is also given in the statement of the lemma. Since ∀z ∈ A, ‖∇2f(z)‖op ≤ Q; from Lemma 10, we conclude that
J(s) =⇒ I(s) ∀s ≥ 0 (15) Hence I(0) also holds true. So the base case of the induction (s = 0) is done, i.e. I(0)∧J(0)∧K(0). Since I(s) implies descent, and K(s) implies that all points outside region A have larger function value, we also have
I(s) ∧K(s) =⇒ J(s+ 1) ∧K(s+ 1) ∀s ≥ 0 (16)
Combining implications equation 15 and equation 16 and Lemma 10, we get
I(s) ∧ J(s) ∧K(s) =⇒ I(s+ 1) ∧ J(s+ 1) ∧K(s+ 1) ∀s ≥ 0 (17)
Hence the induction step also goes through.
We can now use this to prove convergence of gradient descent.
Proof of Lemma 1. To use the adjusted descent lemma (Lemma 9), let us define the appropriate region A.
A := { z : ‖z‖ ≤ ( M
m + 1 ) · ‖x‖ m } (18)
Obviously z(0) = 0 ∈ A. We also have ∀z ∈ Rd \ A, the following holds
f(z)− f(z(0)) = f(z)− f(0)
= 1
2
( ‖G(z)− x‖2 − ‖G(0)− x‖2 ) = 1
2
( ‖G(z)−G(ẑ)‖2 − ‖G(0)−G(ẑ)‖2 ) ≥(0) 1
2
( m2‖z − ẑ‖2 −M2‖ẑ‖2 ) ≥ C+ (m‖z − ẑ‖ −M‖ẑ‖)
≥ mC+ ( ‖z − ẑ‖ − M
m · ‖ẑ‖ ) ≥ mC+ ( ‖z‖ − ( M
m + 1
) · ‖ẑ‖ ) ≥(1) mC+ (( M
m + 1 ) · ‖x‖ m − ( M m + 1 ) · ‖ẑ‖ ) ≥ ( M
m + 1
) ·mC+ · ( ‖x‖ m − ‖ẑ‖ ) ≥(2) 0
C+ = 1 2 (m‖z − ẑ‖+M‖ẑ‖) denotes a non-negative quantity, which does not affect the sign of these inequalities. In the above calculation, (0) follows from Assumptions 1 and 2, (1) uses the definition of A from equation equation 18, and (2) uses Lemma 23. In A, the norm of the Hessian can be bounded as below. The expression of∇2f can be inferred from the expression of∇2L given in equation equation 56. We then get∥∥∇2f(z)∥∥
op ≤ ‖JG(z)‖2op + ∑ i∈[d] |Gi(z)− xi| · ∥∥∇2Gi(z)∥∥op
≤(0) M2 + ∑ i∈[d] |Gi(z)− xi|  ·max i∈[d]
∥∥∇2Gi(z)∥∥op ≤(1) M2 + √ d · ‖G(z)− x‖ ·M2 ≤M2 + √ dMM2 · ‖z − ẑ‖ ≤(2) M2 + √ dMM2 · 2 ( M
m + 1 ) · ‖x‖ m =: Q z ∈ A (19)
where (0) follows from Assumption 2, (1) follows since ∀i ∈ [d], ∥∥∇2Gi(z)∥∥op ≤ ‖∇2G(z)‖op; which follows from the sup-definition of the `2 norm (Lim, 2005). In (1), we also used the fact that ‖v‖1 ≤ √ dim(v) · ‖v‖2. And in (2) we used the diameter of A of equation equation 18.
With A, Q defined through equations equation 18 and equation 19, the conditions of Lemma 9 are satisfied, so we can use the descent equation. We can then use the standard argument below
∀s ∈ [0, S], ∥∥∥∇f(z(s))∥∥∥ ≥ α
=⇒ f(z(s)) ≤ f(z(0))− S · 1 2Q · α2
Since f is a non-negative function, it is bounded below by zero. This means that
min s∈[0,S]
‖∇f(z(s))‖ ≤ √
2Q · f(z(0)) S
We note that ‖∇f(z)‖ = ‖JG(z)T (G(z)− x) ‖
≥ min i∈[d] |σi(JG(z))| ‖G(z)− x‖
≥ min i∈[d] |σi(JG(z))|
√ 2f(z)
≥ m √ 2f(z)
By Lemma 9 we then have
f(z(S)) ≤ 1 2
( mins∈[0,S] ‖∇f(z(s))‖
m
)2
≤ 1 m2
Qf(z(0))
S Furthermore, by Assumption 1 we have
‖z(S) − ẑ‖ ≤ 1 m ‖G(z(S))− x‖
= 1
m
√ 2f(z(S))
≤ 1 m2
√ 2Qf(z(0))
S Hence, setting δ = 1m2 · √ 2Qf(z(0)) S gives us the required number of gradient descent steps. We will finally simplify the expression for S. Since f(z(0)) = f(0) = 12‖G(0)− x‖ 2 = 12‖x‖
2 (from Assumption 3), we get
S = Q · ‖x‖2 m4 · 1 δ2
(20)
Plugging in the expression for Q from equation 19, and using ‖x‖ = O( √ d) from Lemma 24, we get S = O(d 2
δ2 ) as desired.

heading: C PROOF OF LEMMA 2: FAST MIXING OF LANGEVIN DIFFUSION
text: Notation: We first set up some notation. We denote by zt the time t iterate of Langevin diffusion (Definition 3) with g := L, having as stationary distribution P (z) ≡ p(z|x). Further, we denote by Pt the density of zt.
Similarly, we denote by z̃t the iterates of the reflected Langevin chain (Definition 4) with g := L, with stationary distribution P̃ (z) ∝ e−L(z) supported over z ∈ D. We denote by P̃t the density of z̃t. Lemma 11. rad(D) defined in Lemma 13 along with the restriction on β from the same lemma satisfies
rad(D) ≤ min
{ m2
6dMM2 , m√ 2 √ dMM3 } This in particular means that rad(D) ≤ O(1/d).
Proof. Follows by straightforward substitution.
Organization: This section will be organized as follows. In Subsection C.1, we prove a lemma quantifying the difference between the distributions of the regular and reflected Langevin Chains (Lemma 12). To use this lemma, we require a property of the Langevin diffusion for our setting (Lemma 13) — namely, that the chain with probability 1 − remains in the region D. We tackle this in Subsection C.2. In subsection C.3, we prove Lemma 4, namely that p(z|x) restricted to D is log-concave. Finally, using the above three pieces, we will prove prove the main lemma of this section (Lemma 2) in subsection C.4.
We note that some of the techniques here have overlap with techniques used in (Moitra & Risteski, 2020), but we provide all the proofs in detail for completeness.

heading: C.1 RELATING RESTRICTED AND UNRESTRICTED LANGEVIN CHAINS
text: The main lemma of this section is as follows: Lemma 12 (Comparing with reflected Langevin diffusion). Let z̃t follow the reflected Langevin SDE (Definition 4) with g := L), with stationary measure P̃ and with domain D. Let zt follow the Langevin diffusion SDE (Definition 3) with g := L, with stationary measure P . Let P̃t be the density of z̃t and let Pt be the density of zt and let P̃0 be absolutely continuous with respect to the Lebesgue measure.
Then if ∀t ≥ 0,Pr [∃a ∈ [0, t], za /∈ D] ≤ α, it holds that:
dTV(Pt, P ) . 2α+ e − t 2rad2(D) ∀t ≥ 0
Proof. Consider the total variation distance between Pt and P . Using triangle inequality, we have
dTV(Pt, P ) ≤ dTV(Pt, P̃t) + dTV(P̃t, P̃ ) + dTV(P̃ , P ) (21) We first bound the second term using Theorem 3. Since D is convex and L is convex over D by Lemma 4, we have
dTV(P̃t, P̃ ) . e − t
2rad2(D) (22) For the first term, consider the coupling of zt and z̃t such that the Brownian motion Bt is the same. We then have for any t ≥ 0
dTV(Pt, P̃t) ≤ Pr[zt 6= z̃t] ≤ Pr[∃s ∈ [0, t], zt /∈ D] ≤ α (23) where the first inequality follows by Lemma 6, and the second inequality follows from a rewrite that uses the fact that if the unrestricted chain zt stays inside the region D, then the two diffusions won’t differ at all due to the coupling of the Brownian motion. The term α bounds the probability of the unrestricted chain zt leaving the region D.
For the third term, we need to bound dTV(P̃ , P ). Note that P : Rd → R+, P (z) ∝ e−L(z) and P̃ : D → R+, P̃ (z) ∝ e−L(z). This means that ∀z ∈ D, P̃ (z) > P (z) and ∀z ∈ Rd \ D, P (z) > P̃ (z). So we can write
dTV(P̃ , P ) = Pr [z∞ /∈ D] (24) where z∞ is the follows the stationary distribution of the Langevin diffusion P .
From the drift bound in the lemma statement, we can take the limit of t→∞ to get lim t→∞ ∀t ≥ 0,Pr [∃a ∈ [0, t], za /∈ D] ≤ α
=⇒ Pr [∃a ∈ [0,∞), za /∈ D] ≤ α =⇒ Pr [z∞ /∈ D] ≤ α (25)
Using equation equation 24 and equation 25, we get
dTV(P̃ , P ) ≤ α (26)
The lemma statement follows from equations equation 23, equation 22 and equation 26.

heading: C.2 EFFECTIVE REGION FOR LANGEVIN DIFFUSION
text: In this section, we prove Langevin mostly stays in an effective region D, namely: Lemma 13 (Effective region for Langevin diffusion). Let D := {z : ‖z − ẑ‖ ≤ rad(D)} with
rad(D) = 4βm √√√√(2d+ ‖x‖2m2 ) log ( 4 ( 2d+ ‖x‖2 m2 ) ) . Let z0 satisfy ‖z0 − ẑ‖ ≤ 12rad(D). Then with
β ≤ β0, we have that ∀T > 0, Pr [∃t ∈ [0, T ], zt /∈ D] ≤ /4, where
β0 = 4
d √√√√(2d+ ‖x‖2m2 ) log ( 4 ( 2d+ ‖x‖2 m2 )
) ·min{ m36MM2 , d 0.75m2√ 2MM3 }
Towards proving Lemma 13, we will first prove Lemma 3.
Proof of Lemma 3. Using dz = −∇L(z)dt+ √ 2dBt along with Itó’s Lemma, we can compute
dη(z) = 〈∇η(z),−∇L(z)〉dt+ 1 2 · ( √ 2)2∆η(z)dt+ 〈∇η(z), √ 2dBt〉 (27)
We will bound each term on the right hand side. First note that∇η(z) = z − ẑ, and so ‖∇η(z)‖ =√ 2η(z) for all z ∈ Rd. We proceed to the second and third term since they’re straightforward. For
the third term, we have 〈∇η(z), dBt〉 = ‖η(z)‖dBt = √
2η(z)dBt, using a slight abuse of notation to denote the 1-dimensional Brownian motion (in the RHS) as also Bt. For the second term, note that ∆η(z) = d since ∇2η(z) = I . So:
dη(z) = −〈∇η(z),∇L(z)〉dt+ d · dt+ 2 √ η(z)dBt (28)
For the first term, from Lemma 26 we get
〈∇η(z),∇L(z)〉 = 〈z − ẑ,∇L(z)〉 = 〈z − ẑ, ẑ〉+ 〈z − ẑ, ( I + 1
β2 JG(ẑ)
TJG(ẑ) ) (z − ẑ)〉+ 〈z − ẑ, R∇L(z)〉
= 〈z − ẑ, ẑ〉+ ‖z − ẑ‖2 + 1 β2 ‖JG(ẑ)(z − ẑ)‖2 + 〈z − ẑ, R∇L(z)〉 (29)
First, we analyze 〈z − ẑ, ẑ〉. Using Cauchy-Schwartz we can write 〈z − ẑ, ẑ〉 ≥ −‖z − ẑ‖‖ẑ‖. Case 1. ‖z − ẑ‖ ≤ ‖ẑ‖. We can write that 〈z − ẑ, ẑ〉 ≥ −‖ẑ‖2. Case 2. ‖z − ẑ‖ > ‖ẑ‖. We can write that 〈z − ẑ, ẑ〉 ≥ −‖z − ẑ‖2. In either case, it holds that
〈z − ẑ, ẑ〉 ≥ − ( ‖z − ẑ‖2 + ‖ẑ‖2 ) (30)
Second, by Assumption 1, we can say that
‖JG(ẑ)(z − ẑ)‖2 ≥ m2‖z − ẑ‖2 (31) Third, by the expression of the norm of the remainder from Lemma 26, we get
〈z − ẑ, R∇L(z)〉 ≥ −‖z − ẑ‖ · ‖R∇L(z)‖ ≥ − 1 2β2 · ( 3dMM2 · ‖z − ẑ‖+ √ dMM3 · ‖z − ẑ‖2 ) · ‖z − ẑ‖2
For z ∈ D we have ‖z − ẑ‖ ≤ rad(D) so
〈z − ẑ, R∇L(z)〉 ≥ − 1 2β2 · ( 3dMM2 · rad(D) + √ dMM3 · rad(D)2 ) · ‖z − ẑ‖2 z ∈ D
The inequality on rad(D) from Lemma 11 implies that 3dMM2 · rad(D) ≤ m 2 2 and √ dMM3 · rad(D)2 ≤ m 2
2 , and so
〈z − ẑ, R∇L(z)〉 ≥ − m2
2β2 ‖z − ẑ‖2 z ∈ D (32)
Plugging in equation equation 30, equation 31, equation 32 in equation equation 29, we get
〈∇η(z),∇L(z)〉 ≥ −‖ẑ‖2 + m 2 β2 ‖z − ẑ‖2 − m 2 2β2 ‖z − ẑ‖2 = −‖ẑ‖2 + m 2 2β2 ‖z − ẑ‖2 z ∈ D
Using this in equation equation 28, and substituting the definition of η(z), i.e. ‖z − ẑ‖2 = 2η(z), finishes the proof.
Next, we prove high probability deviation bounds for Cox-Ingersoll-Ross processes. Precisely, we show:
Lemma 14 (Concentration Bounds on Cox-Ingersoll-Ross Process). Consider the SDE dXt = −wXtdt+ Ñdt+ 2 √ XtdBt
for Ñ ∈ N and w > 0. Then for any > 0
∀T > 0, Pr [ ∀t ∈ [0, T ], s.t. Xt ≤ 2X0 + 4Ñ
w log
( 4Ñ )] ≥ 1−
Proof. The SDE describes a Cox-Ingersoll-Ross process of dimension Ñ , which equals in distribution to the random variable Yt as defined below, and with matching initial condition Y0 := X0.
Yt = Ñ∑ i=1 ( V (i) t )2 (33)
Here V (i)t follow the Ornstein-Uhlenbeck equation dV (i) t = −w2 V (i) t dt+dB (i) t , and V (i) 0 = √ Y0/Ñ for all i ∈ [Ñ ]. dB(i)t are (1-d) Brownian motions, and we have referenced them separately because they are independent for each i ∈ [Ñ ]. Indeed, applying Itó’s Lemma, we recover the same SDE as:
dYt = Ñ∑ i=1 [( 2V (i) t ( −w 2 V (i) t ) + 1 2 · 2 ) dt+ 2V (i) t dB (i) t ]
= −w Ñ∑ i=1 ( V (i) t )2 dt+ Ñdt+ 2 Ñ∑ i=1 V (i) t dB (i) t
= −wYtdt+ Ñdt+ 2 √√√√ Ñ∑ i=1 ( V (i) t )2 dBt
= −wYtdt+ Ñdt+ 2 √ YtdBt
Here in the RHS of the first expression, the term 12 · 2 is the term 1 2 ·
∂2 ( V
(i) t )2 ∂B (i)2
t
from Itó’s Lemma.
And in the third expression, notice that we have replaced a sum of brownian motions with a single one of matching variance.
Now we can write an explicit solution for the SDE. Since we know the solution of an OrnsteinUhlenbeck process, each V (i)t can be solved as (now we drop the indexing i since the solution is same for each i ∈ [Ñ ])
Vt = V0e −w2 t + ∫ t 0
e− w 2 (t−s)dBs︸ ︷︷ ︸ F (t)
(34)
Applying the reflection principle to F (t), we get
∀T > 0,∀a > 0 Pr [∃t ∈ [0, T ], F (t) ≥ a] = 2Pr[F (T ) ≥ a] (35)
Further, writing the distribution of F (t), we have F (t) ∼ N ( 0, 1
w
( 1− e−wt )) From standard Gaussian tail bounds, we have
∀T > 0,∀a > 0 Pr [F (T ) ≥ a] ≤ exp ( − a 2
2 w (1− e−wT ) ) From equation 35, substituting a := r √ 2 w (1− e−wT ) we have
∀T > 0,∀r > 0 Pr [ ∃t ∈ [0, T ], F (t) ≥ r √ 2
w (1− e−wT )
] ≤ 2e−r 2
Since the quantity F (t) has a symmetric density about zero for all t, we can write the two sided bound as
∀T > 0,∀r > 0 Pr [ ∃t ∈ [0, T ], |F (t)| ≥ r √ 2
w (1− e−wT )
] ≤ 4e−r 2
(36)
Choosing r such that /Ñ = 4e−r 2
in equation 36, and using a union bound over the indices i ∈ [Ñ ], we have ∀T > 0 w.p. (1− ) : ∀t ∈ [0, T ],∀i ∈ [Ñ ], ( V
(i) t − V (i) 0 e
−w2 t )2 ≤ log
( 4Ñ ) · 2 w ( 1− e−wT ) ≤ 2 w log ( 4Ñ )
Plugging in V (i)0 = √ Y0/Ñ for all i ∈ [Ñ ], we get
∀T > 0 w.p. (1− ) : ∀t ∈ [0, T ],∀i ∈ [Ñ ], V (i)t ≤ √ Y0/Ñe −w2 t + √√√√ 2 w log ( 4Ñ )
≤ √ Y0/Ñ + √√√√ 2 w log ( 4Ñ )
Finally, using (a+ b)2 ≤ 2(a2 + b2) and summing over all indices i ∈ [Ñ ], we get
∀T > 0 w.p. (1− ) : ∀t ∈ [0, T ], Ñ∑ i=1 ( V (i) t )2 ≤ 2Y0 + 4Ñ w log
( 4Ñ ) (37)
Replacing the sum by Yt using equation 33 and replacing Xt, X0 in place of Yt, Y0 gives the result.
Next, since all of our SDEs so far have been actually inequalities, we will need a standard comparison theorem for SDEs:
Lemma 15 (Comparison theorem, (Ikeda & Watanabe, 1977)). LetUt,Xt be two stochastic processes following the SDEs
dUt = f(Ut)dt+ h(Ut)dBt
dXt = g(Xt)dt+ h(Xt)dBt
driven by the same Brownian motion Bt. Let at least one of these SDEs have a pathwise unique solution. If ∀x ∈ Rd, f(x) ≤ g(x), and if U0 = X0, then with probability 1,
Ut ≤ Xt ∀t ≥ 0
Lemma 15 shows domination of stochastic variables based on an inequality in their respective SDEs. This is almost what we need, since we showed an inequality in SDEs in Lemma 3 and showed that the upper bound is concentrated in Lemma 14. We need a slight adjustment on top of this though. Note that the inequality between SDEs in Lemma 3 holds only for z ∈ D, however the inequality that the above Lemma requires is on all points α, f(α) ≤ g(α). We can work through this by simply noting that conditioned on Xt (the upper bound SDE) being concentrated in a region, Ut will also be concentrated in the region. We state and show this below.
Lemma 16. Let Ut, Xt be two stochastic processes whose SDEs follow
dUt = f(Ut)dt+ h(Ut)dBt
dXt = g(Xt)dt+ h(Xt)dBt
driven by the same Brownian motion Bt. Let at least one of these SDEs have a pathwise unique solution. If f(α) ≤ g(α) for all α ≤ α0, and if U0 = X0 ≤ α0; then
∀T > 0,∀X̄ ≤ α0 Pr [ ∀t ∈ [0, T ], Ut ≤ X̄ ] ≥ Pr [ ∀t ∈ [0, T ], Xt ≤ X̄ ] Proof. Consider any T > 0 and X̄ ≤ α0.
Pr [ ∀t ∈ [0, T ], Ut ≤ X̄ ] ≥ Pr [ ∀t ∈ [0, T ], Ut ≤ Xt and ∀t ∈ [0, T ], Xt ≤ X̄ ] ≥ Pr [ ∀t ∈ [0, T ], Xt ≤ X̄ ] ·Pr [ ∀t ∈ [0, T ], Ut ≤ Xt | ∀t ∈ [0, T ], Xt ≤ X̄
]︸ ︷︷ ︸ 1
Where the last probability is one, because both Ut and Xt are continuous-time processes, and so Ut > Xt given Xt ≤ X̄ ≤ α0 is not possible because Ut cannot cross over Xt in the f ≤ g region (using comparison Lemma 15).
Given the above preliminaries, we can finally complete the proof of Lemma 13.
Proof of Lemma 13. Instantiate Lemma 16 with Ut := η(zt) and instantiate another SDE on Xt with the same Brownian motion as Ut, as below
dXt := − m2
β2 Xtdt+ d‖ẑ‖2 + de · dt+ 2
√ XtdBt
such that X0 := U0 = η(z0) = 12‖z0 − ẑ‖ 2. From Lemma 14 (plugging in Ñ := d‖ẑ‖2 + de and w := m 2
β2 ), we can conclude that for any > 0
∀T > 0, Pr [ ∀t ∈ [0, T ], s.t. Xt ≤ 2X0 + 4Ñ
w log
( 4Ñ )] ≥ 1− (38)
Using the Lemma 16 on Ut and Xt, with X̄ := 2X0 + 4Ñw log ( 4Ñ ) , and with α0 := 12rad(D) 2 (since the inequality between SDEs of Ut and Xt holds in the region Ut ≤ 12rad(D) 2 from lemma 3), we get
X̄ ≤ α0 ⇒ ∀T > 0, Pr [∀t ∈ [0, T ], Ut ≤ α0] ≥ Pr [ ∀t ∈ [0, T ], Ut ≤ X̄ ] ≥ Pr [ ∀t ∈ [0, T ], Xt ≤ X̄
] ≥ 1−
where the last inequality follows from equation 38. Note that Ut ≤ α0 ≡ ‖zt − ẑ‖ ≤ rad(D) ≡ zt ∈ D. This is exactly the event whose probability we need to lower bound. We now check that the condition for the above is met, namely:
X̄ ≤ α0 ≡ 2X0 + 4Ñ
w log
( 4Ñ ) ≤ 1
2 rad(D)2
Where X0 = 12‖z0 − ẑ‖ 2 ≤ 12 · ( 1 2rad(D) )2 from the lemma statement. This then becomes
≡ 4Ñ w log
( 4Ñ ) ≤ 1
4 rad(D)2, i.e.
rad(D) ≥ 4 √√√√Ñ w log ( 4Ñ )
Putting in the values Ñ , w from above, and using Lemma 23 to bound ‖ẑ‖ (since the RHS above is an increasing function of Ñ ), and bounding the ceiling trivially by dκe ≤ κ+ d (since d ≥ 1) means it suffices that
rad(D) ≥ 4β m √√√√√(2d+ ‖x‖2 m2 ) log 4 ( 2d+ ‖x‖ 2 m2 )  Using the bound on β the above inequality is indeed satisfied for our choice of rad(D).

heading: C.3 PROOF OF LEMMA 4
text: Proof. Using the expression of ∇2L from equation equation 56, we get that
∇2L(z) 1 + 1 β2 m2 − ∥∥∥∥∥∥ ∑ i∈[d] ∇2Gi(z) (Gi(z)− xi) ∥∥∥∥∥∥ op  I (39) And we have ∥∥∥∥∥∥ ∑ i∈[d] ∇2Gi(z) (Gi(z)− xi) ∥∥∥∥∥∥ op ≤ ∑ i∈[d] |Gi(z)− xi| · ∥∥∇2Gi(z)∥∥op
≤ ∑ i∈[d] |Gi(z)− xi| · ∥∥∇2Gi(z)∥∥op ≤(0) √ d · ‖G(z)− x‖ · ∥∥∇2G(z)∥∥ op ≤(1) √ d ·M‖z − ẑ‖ ·M2 ≤ √ dMM2 · rad(D) z ∈ D
In (0), we used the fact that ∀i ∈ [d], ∥∥∇2Gi(z)∥∥op ≤ ‖∇2G(z)‖op; which can be derived from the sup-definition of the tensor `2 norm from (Lim, 2005). In (1), we’ve used Assumption 2. Using the expression of rad(D) from Lemma 11, we get∥∥∥∥∥∥ ∑ i∈[d] ∇2Gi(z) (Gi(z)− xi) ∥∥∥∥∥∥ op ≤ m 2 6 √ d z ∈ D
≤ m2 z ∈ D (since d ≥ 1)
Using this in equation 39, we get∇2L(z) I for z ∈ D.

heading: C.4 PROOF OF LEMMA 2
text: We can now put together the proof of Lemma 2.
Proof of Lemma 2. Using Lemma 12 with D defined as in Lemma 13, with α := /4 in Lemma 13 and setting T = rad2(D) log 1 we get e
−T/2rad2(D) ≤ 2 . Using rad(D) = O(1/d) from lemma 11, the result follows.

heading: D PROOF OF THEOREM 1: EULER-MARIYAMA DISCRETIZATION AND NEURAL SIMULATION
text: In this section, we use the mixing time bounds in the previous sections, along with a discretization error analysis when running a disretized (in the simplest manner, via a Euler-Mariyama scheme) version of the Langevin diffusion chain. This is done in Subsection D.1. Subsequently, we show that each of these steps can be simulated by a “layer” of a neural network. This is done in Subsection D.2. Finally, we put everything together to prove Theorem 1 in Subsection D.3.

heading: D.1 BOUNDING THE ERROR IN THE EULER-MARIYAMA DISCRETIZATION
text: We finally wish to show that with a polynomially small discretization level (i.e. h), the EulerMariyama scheme approximates appropriately for our purposes the Langevin diffusion process. We recall from equation 2 that this Markov process is described by
z ← z − h∇L(z) + √
2hξt, ξ ∼ N(0, I) (40) There is a rich field in the literature on discretizing Langevin diffusion (Dalalyan, 2016; Raginsky et al., 2017; Dalalyan, 2017) primarily used for designing efficient sampling algorithms — a large fraction of them analyzing the above discretization. We will mostly rely on tools from Bubeck et al. (2018). We show: Lemma 17 (Euler-Mariyama discretization). Let zt follow the Langevin diffusion dzt = −∇L(zt)dt+ √ 2dBt, and let Pt denote the distribution of zt.
Let ẑt follow the continuous-time Markov process determined by the Euler-Mariyama discretization scheme with discretization level h given by equation 12, namely
dẑt = −∇L(ẑbt/hch)dt+ √ 2dBt (41)
Let us denote by P̂t the distribution of ẑt. Finally, let z0 = ẑ0 ∈ D satisfying ‖z0 − ẑ‖ ≤ 12rad(D), with D defined as in Lemma 13.
Then, if T = (2rad(D))2 log(1/ ) and h = poly( , 1/d, 1/M,m, β), we have dTV(P̂T , PT ) ≤
Proof. Let yt follow the Reflected Langevin SDE (Definition 4), with g := L and region D, and with the initialization y0 = z0. Let ŷt follow the continuous-time extension of the discretized reflected Langevin SDE given by equation 14 with g := L and region D), with the initial point ŷ0 = z0. Let Qt and Q̂t denote the respective distributions. Let Q denote the stationary distribution of yt.
From triangle inequality, we can decompose our error as dTV ( P̂T , PT ) ≤ dTV (QT , PT ) + dTV ( Q̂T , QT ) + dTV ( P̂T , Q̂T ) (42)
We will bound each of the terms in turn.
For the first term, consider a coupling of yt and zt that uses the same Brownian motion for both processes, so long as neither has left D, and chooses an arbitrary coupling henceforth. By Lemma 6, we have
dTV(QT , PT ) ≤ Pr[zT 6= yT ] ≤ Pr[zT leaves D] + Pr[yT leaves D] ≤ /2
where the last inequality follows by Lemma 13.
For the second term, we can further break it up, again by using the triangle inequality as dTV ( Q̂T , QT ) ≤ dTV ( Q̂T , Q ) + dTV (Q,QT ) (43)
We will use a result about the discretization of reflected Langevin diffusion from (Bubeck et al., 2018). Towards that, let us denote
L1D = max z∈D ‖∇L(z)‖ (44)
L2D = max z∈D ‖∇2L(z)‖op (45)
Using Theorem 1 from (Bubeck et al., 2018), if T = (2rad(D))2 log(1/ ) and h = poly( , 1/d, 1/L1D, 1/L 2 D), both terms in equation 43 are upper bounded by /4. On the other hand, by equations equation 55 and equation 56, we have
‖∇L(z)‖ ≤ ‖z‖+ 1 β2 ‖JG(z)T ‖2‖G(z)− x‖
≤ (‖ẑ‖+ rad(D)) + 1 β2 M2rad(D)2 z ∈ D
Using Lemma 23 and the first term of the RHS from Lemma 11, we get
‖∇L(z)‖ ≤ ‖x‖ m +
( 1 + M2
β2
) m4
36d2M2M22 z ∈ D
≤ O( √ d) +O
( 1
β2d2
) z ∈ D
Which follows using Lemma 24. For the second order bound, a similar calculation as the proof of Lemma 4 we get
‖∇2L(z)‖op ≤ 1 + 1
β2
( ‖JG(z)‖2op + ∑ i ‖∇2Gi‖2|Gi(z)− xi| )
≤ 1 + 1 β2 (M2 +m2)
= O ( 1
β2
)
Plugging the bounds above, we get that if T = (2rad(D))2 log(1/ ) and h = poly( , 1/d, 1/M,m, β), then dTV ( Q̂T , QT ) ≤ /2.
Finally, for the third term, we will show that dTV ( P̂T , Q̂T ) ≤ /4 (46)
Towards that, by Lemma 6, coupling the Gaussian noise to be the same for {P̂t}t∈[T ] and {Q̂t}t∈[T ], we have
dTV ( P̂T , Q̂T ) ≤ Pr[P̂T 6= Q̂T ]
≤ Pr[∃t ∈ [T ], ẑt /∈ D]
By the definition of total variation distance, and Lemma 6 again, we have Pr[∃t ∈ [T ], ẑt /∈ D] ≤ dTV ( {ẑt}t∈[T ], {zt}t∈[T ] ) + Pr[∃t ∈ [T ], zt /∈ D]
By Lemma 13 and Lemma 6, we have Pr[∃t ∈ [T ], zt /∈ D] ≤ /4 and dTV ( {ẑt}t∈[T ], {zt}t∈[T ] ) ≤
/4 which proves equation 46.

heading: D.2 NEURAL SIMULATION OF ARITHMETIC OPERATIONS
text: The deep latent Gaussian will be an overall map from observables x ∈ Rd to latents z ∈ Rd, such that z is an approximate sample from the posterior p(z|x). In this section, we prove “simulation” lemmas that show how to combine neural networks for two functions f, g to create a neural network for some arithmetic operation applied to f, g (e.g. f + g, f × g, etc.). Recall, ρ : R→ R denotes the square activation function, i.e. ρ(x) = x2.
Lemma 18. If f : Rd → R is a neural network with M parameters and g : Rd → R is a neural network with N parameters, both with set of activation functions Σ, then w1f + w2g can be represented by a neural network with O(M +N) parameters and set of activation functions Σ for any w1, w2 ∈ R.
Proof. We simply note that we can make a network with one more layer than the maximum of the depth of f and g that adds up the scaled outputs of f , g to produce w1f + w2g. No additional nonlinearity is needed.
Lemma 19. If f : Rd → R is a neural network with M parameters and g : Rd → R is a neural network with N parameters, both with set of activation functions Σ, then for any w ∈ R, v = wfg can be represented by a neural network with O(M +N) parameters with set of activations Σ ∪ {ρ}. Proof. Noting xy = 14 ( (x+ y)2 − (x− y)2 ) for all x, y ∈ R, we can implement the multiplication of scalars using a two-layer neural network with O(1) parameters using the ρ activation function. We then proceed same as for addition, adding one layer to multiply the outputs of f, g.
Lemma 20. IfG : Rd → Rd is represented by a neural network withN parameters and differentiable activation function set Σ, given an x ∈ Rd and z ∈ Rd, one can compute for any c1, c2 ∈ R:
c1z + c2JG(z) T (G(z)− x)
via a neural network with O(Nd2) parameters and activation functions Σ ∪ Σ′ ∪ {ρ}, where Σ′ contains the derivatives of the activations in Σ.
Proof. Each coordinate of JG : Rd → Rd×d can be written as a neural network with O (N) parameters by Lemma 5, with activation functions Σ ∪ Σ′. Each coordinate of G(z) − x can be computed in O(N) parameters using lemma 18. Performing the multiplications needed for JG(z)T (G(z)− x) via Lemma 19 results in a network of size O(Nd2) due to the d2 size of the Jacobian, and together with the addition of c1z (via O(d) parameters) results in the final bound of O(Nd2) on the parameters.

heading: D.3 PROOF OF THEOREM 1
text: Finally, we can put all the ingredients together to prove Theorem 1.
Proof of Theorem 1. The deep latent Gaussian we construct will have all latent dimensions d as well as the output dimension d. The input will be x ∈ Rd, the observable on which we want to carry out the inference. We will first create a part for the GD simulation, and then one for the Langevin chain.
Using Lemma 1 with δ := 14rad(D), if we run Gradient Descent for S = O( d2 rad(D)2 ) = O( d β2 ) steps (using the bound rad(D) = Ω(β √ d) from Lemma 13), we can get a zinit with ‖zinit − ẑ‖ ≤ 1 4rad(D).
Using Lemma 20 , we can construct a network Done with O ( Nd2 ) parameters that executes a single step of gradient descent. In the deep latent Gaussian we construct, we repeat Done with zero variance O( dβ2 ) times to simulate all the steps of gradient descent and construct zinit. The overall number
of parameters of the deep latent Gaussian so far is then O ( poly(d, 1β ) ·N )
. Note, though the first variable z0 has unit variance, we can easily transform it to a zero variance Gaussian to simulate the first step by multiplying by 0.
Starting with the random point z0, if we run the discretized Langevin dynamics given by equation 11 for T/h steps with T and h as given by Lemma 17, we get by the error bounds of Lemma 17 and Lemma 2 and the triangle inequality:
dTV ( P̂T , p(z | x) ) ≤
which is the required overall bound.
We further note that each step of the discretized Langevin dynamics can be simulated by a layer of the deep latent Gaussian: by Lemma 20 we can construct a network Sone with O ( Nd2 ) parameters that computes the next Langevin iterate. Choosing the variance to be 2h, then, accomplishes exactly the distribution of a single step of the discretized dynamics. Given the choice of T and h, the number of steps we need to run is O (poly(1/ , d, 1/β)). Hence, the overall number of parameters required for simulating the discretized dynamics is O (poly(1/ , d, 1/β) ·N). Finally, it’s clear that the dependence on x in the layers comes through the JG(z)T (G(z)− x) terms in the gradient and Langevin updates—so the weights for the deep latent Gaussian can be produced by a network of size O (poly(1/ , d, 1/β) ·N) as well. The multiplication by x can be implemented via the square activation ρ. Thus, the proof is finished

heading: E PROOF OF THEOREM 2: LOWER BOUND
text: First, we observe the following proposition, which shows that the samples x with high probability lie near a vertex of the hypercube {±1}d
Lemma 21 (Closeness to hypercube). A sample x from the VAE withG(z) = C(sgn(z)) and variance β2I satisfies, for all c > 1
Pr[‖x−G(z)‖ ≤ 6cβ √ d] ≥ 1− exp(−c2d)
Proof. Since x−G(z) is a Gaussian with variance β2I , the claim follows via standard concentration of the norm of a Gaussian.
Furthermore, we also have the following lemma, that shows that when x is near a vertex of the hypercube, the posterior p(z|x) concentrates near the C−1(sgn(x)), that is, the pre-image of the nearest point to x on the hypercube. Lemma 22 (Concentration near pre-image). Let x be such that ‖x− sgn(x)‖ ≤ 6β √ d. Then,
Prz∼p(z|x)[sgn(z) 6= C−1(sgn(x))] ≤ exp(−2d)
Proof. By Bayes rule, since p(sgn(z)) is uniform, we have
p(sgn(z)|x) ∝ p(x|sgn(z))
Consider then any b ∈ {±1}d, b 6= sgn(z). We have:
p(x|b) p(x|sgn(z)) = e−1/β 2‖x−b‖2 e1/β 2‖x−sgn(z)‖2
= e−1/β 2(‖x−b‖2−‖x−sgn(z)‖2)
≤ e−10d
where the last inequality follows for β = o(1/ √ d) and a small enough constant in the o(·) notation. Hence,
p(sgn(z)|x) = 1− ∑
b 6=sgn(z) p(b|x)
≥ 1− 2de−10dp(x|sgn(z)) ≥ 1− e−2d
Proof of Theorem 2. The proof will proceed by a reduction to Conjecture 1. Namely, C be a Boolean circuit satisfying the one-way-function assumptions in Conjecture 1. Let G be defined as G(z) = C(sgn(z))
Suppose that there is an encoder E of size T (n), s.t. with probability at least (d) over the choice of x (wrt to the distribution of G), we have TV(p(z|x), q(z|x)) ≤ 1/10, where q(z|x) is the distribution of the encoder when x is supplied as input.
We will first show how to construct a neural sequential Gaussian C′′ which violates Conjecture 1 in an “average sense”: namely, for 1/poly(d) fraction of inputs x̃ ∈ {±1}d, with probability 1− exp(−d) it outputs the z̃, s.t. x̃ = C(z̃). Subsequently we will remove the randomness — i.e. we will produce from this a deterministic and a neural network, rather than a circuit – subsequently we will show how to remove both the randomness, and use simulate the neural networks using a Boolean circuit.
This deep latent Gaussian C′′ receives a x̃ ∈ {±1}d and constructs (randomly) a multiple samples x ∈ Rd as x = x̃+ ξ, ξ ∼ N(0, β2Id) (47) Subsequently, it samples (multiple) z from E(z|x) and checks whether x̃ = C(sgn(z)). If yes, it outputs sgn(z). If for none of the samples x and z, sgn(z) = x̃, it outputs (1, 1, . . . , 1).
We will show that for 1/poly(d) fraction of inputs x̃ ∈ {±1}d, with probability 1− exp(−d), C̃ ′′ outputs the z̃, s.t. x̃ = C(z̃). Note that x is distributed according to the distribution of the generator G. (Since x̃ is uniformly distributed on the hypercube.) Let E1 be the event that dTV(p(z|x), E(z)) ≤ 1 10 , E2 the event that
‖x− sgn(x)‖ ≤ 6β √ d. By Assumption, Pr[Ē1] ≤ 1− (d); by Lemma 21, for a large enough d, Pr[Ē2] ≤ (d)2 . Hence, Pr[Ē1 ∨ Ē2] ≤ 1− (d) 2 ., i.e. Pr[E1 ∧ E2] ≥ (d) 2 .
Consider an x for which both events E1 and E2 attain. From the the definition of TV distance and Lemma 22, we have
Prz∼E(z|x)[sgn(z) 6= C−1(sgn(x))] ≤ Prz∼p(z|x)[sgn(z) 6= C−1(sgn(x))] + 1
10
≤ exp(−2d) + 1 10 ≤ 1 5
Hence, we have
Prx [ 1 ( Prz∼E(z|x)[sgn(z) 6= C−1(sgn(x))] ≤ 1/5 )] ≥ (d)
2
By the tower rule of probability, we can rewrite this as Eẑ∼N(0,Id) [ Prx∼p(x|ẑ) [ 1 ( Prz∼E(z|x)[sgn(z) 6= C−1(sgn(x))] ≤ 1/5 )]] ≥ (d)/2
Let us denote the random variable A(ẑ) := Prx∼p(x|ẑ) [ 1 ( Prz∼E(z|x)[sgn(z) 6= C−1(sgn(x))] ≤ 1/5 )] such that we have Eẑ[A(ẑ)] ≥ (d)/2 By the reverse Markov inequality, it follows that
Pr[A(z) ≤ a] ≤ 1− (d)/2 1− a
Taking a = (d)/2 we get Pr[A(z) ≤ (d)/2] ≤ 1 − 2(d)/4. Hence, with probability at least 2(d)/4 over the choice of ẑ, we have
Prx∼p(x|ẑ) [ 1 ( Prz∼E(z|x)[sgn(z) 6= C−1(sgn(x))] ≤ 1/5 )] ≥ (d)/2
As we indicated, C′′ will repeat sampling x and z|x multiple times. Let us denote by ξ1 the Gaussian noise sampled in equation 47 and {ξi2}Li=1 the Gaussian samples used in sampling z ∼ q(z|x). Consider sampling {ξi1, i ∈ [1,M ]}, where M = Ω ( 1
(d) log d
) and {ξi ′,j 2 , i ′ ∈ [L], j ∈ [M ′]},
where M ′ = Ω(log d), and for each i ∈ [M ], j ∈ [M ′] outputs x = x̂+ ξi1, and samples from q(z|x) using the randomness in {ξi ′,j 2 }Li′=1. By Chernoff bounds, with probability 1− exp(−Ω(d)), at least one pair of indices (i, j) is such that the corresponding choice of random vectors results in a choice of z, such that C(sgn(z)) = x̂. Since the total number of possible strings x̂ is 2d, taking the constants in M,M ′ sufficiently large, by union bound, there is a pair (i, j), s.t. the corresponding vectors ξi1, {ξ i′,j 2 }Li′=1 works for all x̂. Hard-coding these into the weights of a neural network, we get a deterministic neural network of size O(|C|+ T (d)), that with probability 2(d)/4 inverts C. The only leftover issue is that E is a neural network (hence, a continuous function), whereas we are trying to produce a Boolean circuit. Conversions between neural networks applied to Boolean inputs and Boolean circuits are fairly standard (see Siegelmann (2012), Chapter 4, Theorem 6), but we include a proof sketch for completeness nevertheless.
First, we show that having logarithmic in L,W,M precision for the weights and activations in E suffices. For notational convenience, let us denote by N the maximum degree of any node in E. The proof is essentially the same as Lemma 4.2.1 in Siegelmann (2012) – the only difference is that our inputs are close to binary, but not exactly binary. Namely, for some constants δw, δa to be chosen, consider a network Ê : Rd → Rd, which has the same architecture as E, and additionally:
• The edge weights {w̃i,j} are produced by truncating the weights wi,j to log(δw) significant bits, s.t. ∀i, j : |wi,j − w̃i,j | ≤ δw.
• If a node in E has activation σ, the corresponding node in E′ has activation σ′ that truncates the value of σ to log(δa) significant bits, s.t. ∀x, |σ(x)− σ′(x)| ≤ δa.
Then, we claim that for nodes at depth t from the input, denoting v(x) and v̂(x) the function calculated by the node in E and Ê respectively, we will show by induction:
|v(x)− v̂(x)| ≤ L(LNW )t−1(LNMδw + δa) Let us denote by t the maximum of |v(x)− v̂(x)| for any nodes v, v̂ at depth t. Suppose the claim is true for nodes at depth t and consider a node at depth t+ 1. Denoting by v1, v2, . . . , vn the inputs to v, we have
|v(x)− v̂(x)| = ∣∣∣∣∣σ (∑
i
wivi(x)
) − σ′ (∑ i ŵiv̂i(x) )∣∣∣∣∣ (48) ≤ ∣∣∣∣∣σ (∑
i
wivi(x)
) − σ (∑ i ŵiv̂i(x) )∣∣∣∣∣+ ∣∣∣∣∣σ′ (∑ i ŵiv̂i(x) ) − σ′ (∑ i ŵiv̂i(x) )∣∣∣∣∣ (49)
≤ L ∑ i |wivi(x)− ŵiv̂i(x)|+ δa (50)
≤ L ∑ i (|wivi(x)− wiv̂i(x)|+ |ŵiv̂i(x)− wiv̂i(x)|) + δa (51)
≤ L ∑ i (W t +Mδw) + δa (52)
≤ LNW t + LNMδw + δa (53) where equation 50 follows by Lipschitzness of σ and equation 52 follows from the fact that the weights are bounded by W and the outputs v, v̂ are bounded by M .
Unfolding the recursion, we get
t+1 ≤ t∑ i=0 (LNW )i(LNMδw + δa) ≤ (LNW )t(LNMδw + δa)
Given this estimate, we can choose δw, δa, s.t. the values at each output coordinate of E and Ê match. It suffices to have (LNW )d(LNMδw + δa) < 1, which obtains when δw, δa = O ( 1 M (LNW ) −d) – i.e. it suffices to only keep O (d log(LNW ) + log(M)) significant bits.
Next, we show how to simulate Ê by a Boolean circuit – this is essentially part of the proof of Lemma 4.2.2 in Siegelmann (2012). Namely, each pre-activation can be computed by a subcircuit of size O ( N (d log(LNW ) + log(M)) 2 )
by the classic carry-lookahead algorithm for addition. Since we only require log(δa) bits of accuracy for the activation function, and each activation is at most M (so we need O(logM) bits for the integer part), we need O(logM + log(δa)) nodes for the activation function.
Hence, we can simulate Ê by a Boolean circuit which is at most a factor of poly (N, log(L), log(W ), log(M)). Since L,W,M are o(exp(poly(d))), the size of the size of the resulting Boolean circuit is still poly(d). Thus, we’ve constructed a polynomial-sized circuit which inverts the one-way-permutation, thus violating Conjecture 1.

heading: F TECHNICAL LEMMAS
text: This section contains several equations and calculational lemmas that we repeatedly use across the proofs of Theorem 1.
F.1 BOUND ON ‖ẑ‖
Recall that we defined the inverse point ẑ for a given x ∈ Rd in equation equation 4. As mentioned earlier, due to bijectivity, it is unique and satisfied G(ẑ) = x. Further, we can easily prove a bound on the norm of ẑ using our assumptions, which will be useful in various places in the proof.
Lemma 23. ẑ defined as in equation 4 with G that satisfies Assumption 1, is such that
‖ẑ‖ ≤ ‖x‖ m
Proof. Using strong invertibility from Assumption 1 and centering from Assumption 3, we can write
‖x−G(0)‖ = ‖G(ẑ)−G(0)‖ ≥ m · ‖ẑ − 0‖
=⇒ ‖ẑ‖ ≤ ‖x‖ m
(54)
F.2 HIGH PROBABILITY BOUND ON ‖x‖
Lemma 24 (High-probability bound on ‖x‖). With x following the density from Definition 1 (with dl = do = d), s.t. G satisfies Assumption 2, we have that
‖x‖ ≤ 12(M + β) √ d w.p. 1− 2 exp(−4d)
Proof. The density of x is defined by the latent Gaussian. Let’s denote the error term as e ∼ N (0, β2I). We have
x = G(z) + e
=⇒ ‖x‖ ≤ ‖G(z)‖+ ‖e‖ =⇒ ‖x‖ ≤M‖z‖+ ‖e‖
(Where we used assumption 2). Now z and e follow zero-mean Gaussian densities, with variances I and β2I respectively. Using standard gaussian tail bounds (similar to lemma 21 with c = 2), we get the required claim.
Remark: Since we work under the setting β ≤ O(1) from equation 1, we can use ‖x‖ ≤ O( √ d) with probability 1− exp(Ω(d)).
F.3 CHARACTERIZING LOSS NEAR ẑ
We need to understand how L and ∇L behaves near ẑ for several claims in the proofs. We will use a Taylor expansion, bounding the contribution of higher order terms.

heading: F.3.1 EXPRESSIONS FOR DERIVATIVES
text: Using equation equation 3, we first calculate the expressions of the derivatives as
∇L(z) = z + 1 β2 JG(z) T (G(z)− x) (55)
∇2L(z) = I + 1 β2 JG(z)TJG(z) + ∑ i∈[d] ∇2Gi(z) (Gi(z)− xi)  (56) At the inverse point ẑ, since G(ẑ) = x, we have: L(ẑ) = 12‖ẑ‖ 2
∇L(ẑ) = ẑ ∇2L(ẑ) = I + 1β2 JG(ẑ) TJG(ẑ) I (57)
Since∇2L(ẑ) I , the function L(z) is strongly convex in the neighbourhood of ẑ.

heading: F.3.2 BOUND ON 3RD DERIVATIVE OF L
text: To show the behaviour of∇L around ẑ, we will first bound the norm of the 3rd order derivatives of L below. Lemma 25 (Bounding the 3rd derivative of L). For all z ∈ Rd, it holds that∥∥∇3L(z)∥∥ op ≤ √ dM β2 · ( 3 √ dM2 +M3‖z − ẑ‖
) Proof. An elementary calculation shows that
∇3L(z) = 1 β2 3 ∑ i∈[d] ∇Gi(z)⊗T ∇2Gi(z) + ∑ i∈[d] ∇3Gi(z) (Gi(z)− xi)  (58) Here, the tensor product ⊗T multiplies elements from Rd and Rd×d to give elements in Rd×d×d. We subscript with T to avoid confusion with the kronecker product.
Then, we have:∥∥∇3L(z)∥∥ op ≤ 1 β2 3 · ∥∥∥∥∥∥ ∑ i∈[d] ∇Gi(z)⊗T ∇2Gi(z) ∥∥∥∥∥∥ op + ∥∥∥∥∥∥ ∑ i∈[d] ∇3Gi(z) (Gi(z)− xi) ∥∥∥∥∥∥ op  (triangle inequality) ≤ 1 β2 3 ·∑ i∈[d] ∥∥∇Gi(z)⊗T ∇2Gi(z)∥∥op + ∑ i∈[d] |Gi(z)− xi| · ∥∥∇3Gi(z)∥∥op
 ≤ 1 β2 3 · ∑ i∈[d] ‖∇Gi(z)‖  ·max i∈[d] ∥∥∇2Gi(z)∥∥op + ‖G(z)− x‖1 ·maxi∈[d] ∥∥∇3Gi(z)∥∥op 
≤(1) 1 β2
( 3 · √ d · ‖JG(z)‖F · ∥∥∇2G(z)∥∥ op + √ d · ‖G(z)− x‖ · ∥∥∇3G(z)∥∥ op ) ≤(2) 1
β2
( 3d ·MM2 + √ d ·MM3 · ‖z − ẑ‖ ) In (1) we use ∑ i∈[k] |αi| ≤ √ k · √∑ i∈[k] α 2 i on the first part of both the terms; and the fact that
∀i ∈ [d], ‖T (i)‖op ≤ ‖T‖op for any tensor T (where T (i) denotes a sub-tensor) on the second part of both terms. This fact follows from the sup-definition of `2 norm on tensors from (Lim, 2005). In (2), we use Assumption 2, as well as ‖A‖F ≤ √ α · ‖A‖op for A being an α× α matrix.

heading: F.3.3 PERTURBATION BOUND OF ∇L
text: Using the above, we show the Taylor expansion result for the gradient of L: Lemma 26 (Perturbation of ∇L). For all z ∈ Rd, it holds that
∇L(z) = ẑ + ( I + 1
β2 JG(ẑ)
TJG(ẑ) ) (z − ẑ) +R∇L(z)
where ‖R∇L(z)‖ ≤ √ dM 2β2 · ( 3 √ dM2 +M3‖z − ẑ‖ ) · ‖z − ẑ‖2.
Proof. We Taylor expand ∇L around ẑ. Using equation equation 57 we get : ∇L(z) = ∇L(ẑ) +∇2L(ẑ) (z − ẑ) +R∇L(z)
= ẑ + ( I + 1
β2 JG(ẑ)
TJG(ẑ) ) (z − ẑ) +R∇L(z) (59)
To bound the norm of the remainder term, we note that
‖R∇L(z)‖ = 1 2! · ∥∥∥∇3L(z)∣∣ z=zmid ∥∥∥ op · ‖z − ẑ‖2
For zmid = tmidz + (1− tmid) ẑ for some tmid ∈ [0, 1]. And for the third derivative, we use lemma 25 directly to get
‖R∇L(z)‖ ≤ √ dM 2β2 · ( 3 √ dM2 +M3‖z − ẑ‖ ) · ‖z − ẑ‖2 (60)
This finishes the proof.

heading: G REMARK ON ACTIVATION FUNCTIONS
text: The choice of activation functions in our results is largely a matter of convenience. Using standard techniques from approximation theory, e.g. Hornik (1991); Yarotsky (2017), one can approximate a neural network with one choice of nonlinearity via a (comparably sized) neural network with another choice of nonlinearity, under very mild conditions on the nonlinearities. Crucially, this simulation only increases the size by a dimension-independent factor. This result frees us (for purposes of deriving an expressibility result) to work with activation functions chosen for mathematical convenience and produce results that hold without loss of generality.
For instance, if we wished to use ReLU activations, we can use the following lemma. We note that this proof is almost verbatim the same as the proof of Lemma 1.3 in Telgarsky (2017) and proofs for other activations like sigmoid or tanh can be written completely analogously. Lemma 27. Let Ω ⊆ [−M,M ]d and let G1 : [−M,M ]d → R be a neural network with at most l layers and n parameters, s.t., the weights W (i) for each layer i and node j in G1 are bounded as, i.e., ∀i, j : ∑ k |W (i) j,k | ≤ B. Furthermore, assume that the activation functions used in G1 belong to the set Ξ, s.t. all functions σ : R→ R satisfy supx∈[−B·M,B·M ] σ < M . Then there exists a neural network G2 with ReLU activation and O(n (LB)llBM log( (LB)llBM )) parameters, such that for all supx∈[−M,M ]d |G1(x)−G2(x)| ≤ .
Proof. For any σ ∈ Ξ, from Theorem 1 in Yarotsky (2017) it follows that there exists a neural network R with ReLU activations and O(LBM ′ log( LBM ′ )) parameters such that supx∈[−B·M,B·M ]|σ(x)− R(x)| ≤ ′. We will construct the network G2 by replacing each activation in G1 with the corresponding network R as given by the result above with ′ = /l. Note, this network is at most a factor of O( (LB) llBM log( (LB)llBM )) bigger than G1, as the lemma requires.
We will prove the claim of the lemma by induction on l. More precisely, we will show (by induction) that for each node at layer i, the network G2 calculates a function that is (LB)ii ′ away in l∞ norm from the corresponding node in G1, and the inputs to the node are in [−BM,BM ].
For the base case i = 1, since the input x ∈ [−M,M ]d, the result follows by Theorem 1 in [Yarotsky (2017)].
We proceed to the inductive claim. Let H(x) denote the vector valued mapping computed by the nodes at layer i, and let HR(x) be the corresponding vector in G2. As inductive hypothesis, we assume that ‖H(x) − HR(x)‖∞ ≤ (LB)ii ′ for all x ∈ [−M,M ]d and ‖H(x)‖∞ ≤ M as well as ‖HR(x)‖∞ ≤ M . Therefore, for the jth node in layer (i + 1) in network G1 we have |WTj H(x)| ≤ ‖Wj‖1‖H(x)‖∞ ≤ BM and σ is bounded by M on this interval, so we have ‖σ1(WTj H(x))‖∞ ≤ M . Along with the bound on the activations, the part of the inductive hypothesis about the size of the input is proven. To prove the error bound, we have:
|σ1(WTj H(x))−R(WTj HR(x))| ≤ |σ1(WTj H(x))− σ1(WTj HR(x))|+ |σ1(WTj HR(x))−R(WTj HR(x))| ≤ L|WTj (H(x)−HR(x))|+ ′
≤ L‖Wj‖1‖H(x)−HR(x)‖∞ + ′
≤ (LB)i+1(i+ 1) ′
This finishes the proof of the inductive step, and thus the lemma.

Review summaries:
This purely mathematical paper investigates the important question of how does the necessary level of complexity of an inference subnetwork depend on the the complexity of the corresponding generative subnetwork, in a VAE model. The paper introduces a specific measure of invertibility, and uses it to show that a generative subnetwork that is sufficiently strongly invertible only requires an inference model at similar complexity, while some "non-invertible" generative subnetworks will require an exponentially more complex inference subnetwork. The paper provides a lengthy and comprehensive proof of both.

The paper answers the question how complex inference models need to be to accurately estimate posterior distributions. The conclusion is when a latent Gaussian model with N parameters satisfies (1) strong invertibility, (2) 3-th smoothness, the posterior can be approximated by a deep latent Gaussian model with O(N) parameters. A special case is also given to show that strong invertibility is necessary to the conclusion.