PLUGIn: A simple algorithm for inverting generative models with recovery guarantees
Year: 2021

Abstract:
We consider the problem of recovering an unknown latent code vector under a known generative model. For a d-layer deep generative network G : Rn0 ! Rnd with ReLU activation functions, let the observation be G(x) + ✏ where ✏ is noise. We introduce a simple novel algorithm, Partially Linearized Update for Generative Inversion (PLUGIn), to estimate x (and thus G(x)). We prove that, when weights are Gaussian and layer widths ni & 5n0 (up to log factors), the algorithm converges geometrically to a neighbourhood of x with high probability. Note the inequality on layer widths allows ni > ni+1 when i 1. To our knowledge, this is the first such result for networks with some contractive layers. After a sufficient number of iterations, the estimation errors for both x and G(x) are at most in the order of p 4n0/ndk✏k. Thus, the algorithm can denoise when the expansion ratio nd/n0 is large. Numerical experiments on synthetic data and real data are provided to validate our theoretical results and to illustrate that the algorithm can effectively remove artifacts in an image.

text: p 4dn0/ndk✏k. Thus, the algorithm can denoise when the expansion ratio
nd/n0 is large. Numerical experiments on synthetic data and real data are provided to validate our theoretical results and to illustrate that the algorithm can effectively remove artifacts in an image.

heading: 1 Introduction
text: We consider the inverse problem of recovering an unknown vector x⇤ 2 Rn0 from a noisy observation y 2 Rnd of the form y = G(x⇤) + ✏, (1) where ✏ 2 Rnd is noise and G : Rn0 ! Rnd is a known d-layer feed-forward neural network with ReLU activation functions (in fact, it is straightforward to generalize to positively homogeneous activation functions, which we will elaborate later). Precisely, G has the form
G(x) = (Ad (Ad 1 . . . (A1x) . . .)), (2)
where (·) = max(·, 0) is the ReLU activation function and Ai 2 Rni⇥ni 1 is the weight matrix in the i-th layer.
The inverse problem (1) has applications in, for example, signal denoising and signal compression [1, 2, 3]. In the case of denoising the typical goal is to recover a clean signal y⇤ from its noisy observation y⇤ + ✏ and, in the case of compression, the goal is to find an efficient low dimensional representation of y⇤. Traditional approaches for these problems often use priors on signals, for example, a sparsity prior with respect to a fixed basis or dictionary [4, 5, 6]. An emerging viewpoint is to use a generative prior that assumes the unknown signal y⇤ is in the range of a deep generative model G, i.e., y⇤ = G(x⇤), and develop algorithms that can estimate the latent code vector x⇤ from y ⇤ + ✏, thus recovering y⇤.
Recent advancements in training deep neural networks have shown that priors in the form of generative models can effectively map low dimensional vectors to the space of natural image classes
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
[7, 8, 9]. Learned generative models can then be used as priors to solve various inverse problems including denoising [1, 10], compressive sensing [11, 12, 13, 14, 15, 16], phase retrieval [17], blind deconvolution [18, 19], low-rank matrix recovery [20] and have been shown to perform on par or outperform classical sparsity based approaches for these inverse problems. For example, in [1] the authors empirically showed that an end-to-end approach for denoising using a neural network that maps noisy patches in an image to noise-free ones achieves state-of-the-art performance and is on par with BM3D. Similarly, in [11] the authors empirically showed that for compressive sensing using generative prior, optimization of the empirical risk objective over the latent code space (of the generative prior) can recover a vector that effectively estimates the uncompressed signal with 5-10 times less measurements compared to Lasso in some cases.
Given that y equals G(x⇤), with possibly some additive noise, a standard way to estimate x⇤, i.e., to invert G(x⇤), would be to look for a minimizer of the program
min x2Rn0
ky G(x)k2. (3)
Unfortunately, this program is non-convex and to our knowledge there is no known efficient method that can achieve its global minimum in general. On the other hand, for generative networks with random (Gaussian) weight matrices, a line of papers showed that gradient-based algorithms can provably avoid local minima with high probability [10, 14, 17, 20]. In particular, [10] considers a random Gaussian noise model, and random Gaussian weight matrices Ai which are highly expansive at each layer. Under these conditions, the authors show that the latent code vector x⇤ can be accurately estimated using a gradient-based method that uses the (sub-)gradient updates given by
x k+1 = xk ⌘(D1A1)|(D2A2)| · · · (DdAd)| G(xk) y , (4)
where xk is the k-th estimate, ⌘ 2 R is step size, and Dj is a diagonal matrix with entries that are either zero or one. Each Dj zeros out the inactive rows of Aj with respect to the estimate xk and so it is a function of xk (and Ap for p < j). Thus, at each iteration all Dj need to be updated.
In this paper, we show that we can drop the Dj and still recover an accurate estimate of x⇤. We propose the following iterative algorithm, Partially Linearized Update for Generative Inversion (PLUGIn), to estimate x⇤:
PLUGIn: x k+1 = xk ⌘A|1A | 2 · · ·A | d G(xk) y . (5)
This algorithm was inspired by previous work showing that latent vectors for non-linear single-index function can be approximately estimated by treating the function as linear [21, 22]. Similar to [10, 11, 12, 14, 17, 19, 23], for theoretical analysis, we assume the weight matrices are Gaussian. To show that the algorithm works more broadly, we conduct real data simulations. Applying the ideas in [21, 22] one can show that for any fixed x0, the first iteration of (5) provides an unbiased estimate of x⇤ with ⌘ = 2d, which is generally not the case for the gradient descent estimates given by (4), see Section 3 for further details. Additionally, each iteration of PLUGIn maps the difference G(xk) y to the low dimensional latent code space using a static matrix A|1A | 2 . . . A | d , which can be pre-multiplied and reused in subsequent iterations.
In addition to the randomness assumption on the weight matrices, when the width of the each layer satisfies ni & 5in0 (up to log factors), noise ✏ is independent of weight matrices, and the step size ⌘ is fixed, we show that the estimates provided by PLUGIn converge to a neighbourhood of x⇤ geometrically with high probability. This result allows ni > ni+1 for i 1 and thus can provide theoretical guarantees for estimating x⇤ from (1) even when G has some contractive layers. To our knowledge, this is the first such result. Additionally, we show that the k-th recovery and reconstruction errors satisfy
kxk x⇤k . s
4dn0 nd
k✏k and kG(xk) G(x⇤)k . s
4dn0 nd k✏k
with high probability when k is large. Thus, PLUGIn can effectively denoise noisy observations satisfying (1) provided that the expansion ratio nd/n0 is large. These contributions are highlighted below:
1. We introduce a simple novel algorithm (5), called Partially Linearized Update for Generative Inversion (PLUGIn), for inverting a generative model. Note that each update of PLUGIn
can be implemented with two blackboxes. The first blackbox provides the G(xk) for each iterate and the second blackbox provides the product of the weight matrices A|1A | 2 . . . A | d .
Thus, given these two blackboxes the individual matrices are not needed. 2. Under relatively mild conditions, we show that for generative priors with possibly contractive
layers, PLUGIn can effectively estimate the unknown latent code x⇤ and the unknown signal G(x⇤) from observation y satisfying (1). The convergence is geometric and holds for a range of step sizes.
3. In the presence of arbitrary noise ✏, we show that the recovery and reconstruction error corresponding to the estimate xk, for sufficiently large k, is approximately p 4dn0/ndk✏k.
In contrast to typical random noise assumptions, we only require the noise to be independent of the weight matrices. This allows stochastic noise (independent of weight matrices) as a special case if one can bound the norm of the stochastic noise with high probability.
Organization of the paper: In Section 1.1, we introduce notations used through the paper. In Section 2, we state the main results. In Section 3, we provide an outline of the proof. In Section 4, we provide numerics on synthetic and real data. In Section 5, we briefly discuss a generalization of activation functions, as well as limitations in our results.

heading: 1.1 Notations
text: For a positive integer n, let [n] = {1, 2, . . . , n}. For a vector x, let kxk be its Euclidean norm; for a matrix A, let kAk be its operator norm; for a matrix A and a set T , let kAkT := supx2T \{0} kAxk kxk . Let N (0, t2) be the normal distribution with mean zero and standard deviation t. Let B(x, r) be the Euclidean ball of radius r centered at x and let Bn(0, r) be the Euclidean ball in Rn with radius r, centered at origin. We use C and c to denote absolute constants (often c for small ones and C for large ones) which may vary from line to line. We also use C0, C1, etc., to denote particular absolute constants, which do not change throughout the paper.

heading: 2 Main Results
text: Let the generative model G : Rn0 ! Rnd be defined as in Equation (2). We consider the inverse problem:
Let: x⇤ 2 Rn0 , ✏ 2 Rnd , y = G(x⇤) + ✏, Given: generative model G and observation y, Estimate: latent code vector x⇤ and G(x⇤).
We assume:
A1. All activation functions are ReLU. A2. Each Ai 2 Rni⇥ni 1 has i.i.d. N (0, 1/ni) entries and {Ai}i2[d] are independent. A3. Layer widths (number of nodes in each layer) satisfy
ni C05in0 log
0 @ i 1Y
j=0
enj
n0
1
A , i 2 [d] (6)
for some (sufficiently large) absolute constant C0. A4. The noise ✏ does not depend on {Ai}i2[d]. (The noise may be deterministic or random.)
Let x0 2 Rn0 and let xk be the result of applying k iterations of PLUGIn (5). Under assumptions A1-A4, xk geometrically converges to a neighbourhood of x⇤ (and also G(xk) to a neighbourhood of G(x⇤)) for a range of step sizes near 2d. Precisely, we have the following theorem. Theorem 1. Let ✓ 2 (0, 43 ) and let ↵ = |1 ✓|+ 1 2✓ 2 (0, 1). Let R be a positive number such that kx0 x⇤k  R. Under assumptions A1-A4, the k-th estimate xk given by PLUGIn algorithm (5) with constant step size ⌘ = ✓2d satisfies
kxk x⇤k  ↵kR+ 15✓ 1 ↵2
d p
n0/ndk✏k, and
kG(xk) G(x⇤)k  3↵kR+ 45✓ 1 ↵2
d p n0/ndk✏k
with probability at least 1 2(k + 3)e 10n0 .
When ✓ = 1, Theorem 1 reduces to the following corollary. Corollary 1. Let R be a positive number such that kx0 x⇤k  R. Under assumptions A1-A4, the k-th estimate x k given by PLUGIn algorithm (5) with constant step size ⌘ = 2d satisfies
kxk x⇤k  2 kR+ 30 · 2d p n0/ndk✏k, and
kG(xk) G(x⇤)k  2 k(3R) + 90 · 2d p n0/ndk✏k
with probability at least 1 2(k + 3)e 10n0 . Remark 1 (Allowed number of iterations and comparison of error bound summands). The probability in Theorem 1 (and also Corollary 1) decreases as the number of iterations increases. However, this is generally not a problem because we can take up to exponentially many iterations before this probability becomes trivial. In fact, if we require a probability of at least 99% on the error bounds, then we can take any k  k?(n0) := 1200e
10n0 3. Note that k? grows exponentially in n0 and it is already large when n0 is relatively small1.
The noise term in the error bounds (proportional to k✏k) will dominate the other term (long) before k ? many iterations are reached, unless the noise is virtually zero. Indeed, in Corollary 1, if 2 k ? R p 4dn0/ndk✏k, then the noise must satisfy
k✏k < R q
nd 4dn0 1 2
k? = 8R
q nd
4dn0 1 2
exp(10n0)/200 .
Note that 1 2 exp(10n0)/200 is a double exponential and extremely small2, so we will consider the noise as being virtually zero in this case. On the other hand, if 2 k ? R < p 4dn0/ndk✏k, then after at most k? iterations, with probability at least 99%, the estimation errors for x⇤ and G(x⇤) will be at most C p 4dn0/ndk✏k. The same argument can be made for Theorem 1 when ↵ is close to 12 (i.e., ✓ is close to 1) so that ↵k decays fast and 11 ↵ is not close to being singular.
Remark 2 (Contractive layers). In A3, (6) states a lower bound on ni with respect to the latent code dimension n0 (up to log factors). While this bound strictly increases with layer depth i, it is not necessary for ni to always increase with i (except in the first layer). For example, consider ni = C05dn0d(2d i) where is any fixed number such that C0 2 N and 4 + logC0. It is easy to see n1 > n2 > · · · > nd, and we can also verify (see Appendix D) that such ni satisfy (6). In this case, the network is contractive in each layer after the first, and Theorem 1 still applies. Remark 3 (Initialization may depend on random weight matrices). The results of the theorem can still hold when x0 is chosen randomly, dependent on the weight matrices Ai. In this case, suppose that kx0 x⇤k  R with probability at least 1 . Then, the error bounds hold with probability at least 1 2(k + 3)e 10n0 . This does not follow directly from the theorem as stated (which fixes x 0, then takes random weight matrices), but follows from the proof.
Remark 4 (Comparison to guarantees for gradient-based method). Here we compare our results to the ones in [10], which uses (4) for iterations and considers a random noise model.
When the noise ✏ ⇠ N (0, 1nd 2 Ind), the error bounds on the k-th PLUGIn iterate in Corollary 1, with high probability, reduces to
max{kxk x⇤k, kG(xk) G(x⇤)k}  C ⇣ 2 kR+ 2d p n0/nd ⌘ .
Thus, we get geometric convergence with a rate of 1/2 and after sufficiently many iterations, the errors are at most C2d p n0/nd .
A similar result was shown in [10], which studies a gradient-based method for solving the inverse problem (1). They show that when the weight matrices are random Gaussian and sufficiently
1For example, k?(5) ⇡ 2.59⇥ 1019 and k?(10) ⇡ 1.34⇥ 1041. 2For example, 2 exp(10)/200 ⇡ 7⇥ 10 34.
expansive at each layer, the iterates of the gradient-based method converge to a neighborhood of the target signal x⇤. The result holds for a fixed step size, dependent on d. After sufficiently many iterations, the iterates converge geometrically to a neighborhood of x⇤ of radius at most on the order of p n0/nd , up to log factors. This rate of convergence takes the form (1 C/d2), thus giving slower convergence for deeper nets. On the other hand, we note that dependence on d is of relatively minor concern. Generative models usually have small depth in practice, our MNIST experiments (below) work well with depth 3, and typical applications use depth less than 8.

heading: 3 Proof Outline
text: Here we give a sketch for the proof of Theorem 1. For simplicity, we will only focus on analyzing one iteration of PLUGIn with step size ⌘ = 2d. The complete proof can be found in appendices.

heading: The Special Case
text: Let us first look at the special case where d = 1 and ✏ = 0. The analysis here highlights some of the key ideas in our proofs, while its result Lemma 1 serves as a building block for proof in the general case. In this special case, PLUGIn with ⌘ = 2d reduces to
x k+1 = xk 2A| ⇥ (Axk) (Ax⇤) ⇤
where = ReLU and A 2 Rm⇥n is random with i.i.d. N 0, 1m entries.
In fact, the first iterate provides an unbiased estimate of x⇤ when x0 does not depend on A. Indeed, the rotation invariance property of the Gaussian distribution may be leveraged to show [21, 22], for any fixed x, EA| (Ax) = 12x. (7) For completeness, we also include a proof for (7) in Appendix A, Lemma 2. Applying (7) to the first iteration gives
Ex1 = x0 2EA| (Ax0) + 2EA| (Ax⇤) = x0 x0 + x⇤ = x⇤
Thus, even the first iterate can be shown to be a good estimate by showing that x1 concentrates around its mean. Further iterates are generally no longer unbiased estimators because they pick up complex dependence on the random matrix A. We overcome this by developing a series of uniform deviation inequalities, as below.
Let us suppose we have shown that, with high probability, kxk x⇤k  r for some (small) constant r > 0. Then we wish to show that kxk+1 x⇤k  r/2 with high probability. Notice that
(xk+1 x⇤) = 2A| ⇥ (Axk) (Ax⇤) ⇤ (xk x⇤)
kxk+1 x⇤k = sup u2Sn 1
2 ⌦ Au, (Axk) (Ax⇤) ↵ ⌦ u, x k x⇤ ↵
= 2 sup u2Sn 1
Z(u, xk;x⇤)
where Z(u, v;x⇤) := hAu, (Av) (Ax⇤)i 12 hu, v x
⇤i . We wish to bound the supremum of random process Z(u, xk;x⇤) over u 2 Sn 1. However, this process is challenging to analyze since xk depends on A when k 1. To alleviate this dependency, we bound by the supremum of Z(u, v;x⇤) over (u, v) 2 T 0 := Bn(0, 1) ⇥ B(x⇤, r) instead. It is worth noting that Z(u, v;x⇤) is centred, namely EZ(u, v;x⇤) = 0 for any fixed (u, v). We now arrive at the estimate
kxk+1 x⇤k  2 sup (u,v)2T 0 Z(u, v;x⇤) if kxk x⇤k  r. (8)
The following Lemma 1 provides a bound on supT 0 Z(u, v;x⇤). In fact, it is slightly more general because we replaced T 0 with T1 ⇥ T2 (this replacement is helpful when studying the general case d > 1). The proof of this lemma can be found in Appendix B. The proof idea is to first establish
that Z(u, v;x⇤) has mixed (sub-Gaussian and sub-exponential) tail increments through Bernstein’s inequality, and then apply the result from [24], which provides a general bound for the supremum of random processes with mixed tail increments. Lemma 1. Let = ReLU. Fix w 2 Rn and let A 2 Rm⇥n have i.i.d. N
0, 1m entries. Define
Z(u, v;w) := hAu, (Av) (Aw)i 12 hu, v wi . Suppose T1, T2 are sets (not depending on A) such that

heading: T1 = S1 \ Bn(0,↵) and T2 = S2 \ B(w,↵r)
text: for some q-dimensional (affine) subspaces S1,S2 ✓ Rn and real numbers ↵, r > 0. Then for any t 1,
sup u2T1 v2T2
|Z(u, v;w)|  C1↵2r r q
m +
q m +
r t
m +
t
m
!
with probability at least 1 e t. Here C1 > 0 is an absolute constant.
We can apply Lemma 1 to estimate (8) (with S1 = S2 = Rn) and get, for example,
kxk+1 x⇤k  2C1r ⇣q n m + n m + q n m + n m ⌘  12r
with probability at least 1 e n, provided that m (16C1)2n.

heading: The General Case
text: Let us illustrate the proof idea with d = 2 (the extension to d > 2 is straightforward). Denote x k 1 = (A1x k) and x⇤1 = (A1x⇤). By adding and subtracting 2A | 1(x k 1 x⇤1) we can write (5) with ⌘ = 2d as
x k+1 x⇤ = xk x⇤ 22A|1A | 2 [G(xk) G(x⇤) ✏]
= (xk x⇤) 2A|1 (A1x k) (A1x⇤)
+ 2A|1 ⇥ (xk1 x⇤1) 2A | 2 (A2x k 1) (A2x⇤1) ⇤ + 22A|1A | 2✏.
Similar to the special case above, we can get
kxk+1 x⇤k  sup u2Sn0 1 2Z1(u, x k) + sup u2Sn0 1 22Z2(A1u, x k 1) + 2 2kA|1A | 2✏k (9)
where (denote x⇤0 = x⇤)
Zj(u, v) := ⌦ Aju, (Ajv) (Ajx⇤j 1) ↵ 12 ⌦ u, v x⇤j 1 ↵ , j = 1, 2.
Also assume that kxk x⇤k  r, it remains to bound each term on the right hand side of (9). The first term can be bounded directly through Lemma 1 (with t = 10n0). The last term is also easy to bound by the randomness of Ai (see Appendix C, Lemma 6), in which case we have kA|1A | 2✏k  15 p n0/ndk✏k with high probability.
Denote G1(x) = (A1x) for x 2 Rn0 . For the second term, first notice that range(A1) is a n0dimensional subspace in Rn1 . Using the ideas from [11, 25], we can also show that range(G1) is contained in a union of N many n0-dimensional (affine) subspaces, where N  (en1/n0)n0 . Furthermore, let E be the event such that mappings A1,G1,G all have Lipschitz constants being at most 3, then we can show (Appendix C, Lemma 8) that P(E) 1 3e 10n0 . Also on event E (note that kA1k  3 and kxk1 x⇤1k  3r), we have
A1Sn0 1 ✓ range(A1) \ Bn1(0, 3) = S1 \ Bn1(0, 3) =: T1 x k 1 2 range(G1) \ B(x⇤1, 3r) ✓ [j2[N ] (S1,j \ B(x⇤1, 3r)) =: [j2[N ]T2,j
where S1 and S2,j are n0-dimensional (affine) subspaces. Applying Lemma 1 on each T1 ⇥ T2,j , followed by a union bound over j 2 [N ], we get (denote T2 = [j2[N ]T2,j)
sup T1⇥T2
Z2(u, v)  C1(9r) ⇣q
n0 n2 + n0n2 + q t n2 + tn2 ⌘
with probability (over A2 and conditioning on A1) at least 1 Ne t. By choosing t = 2n0 log(en1/n0), we obtain a high probability bound for supu2Sn0 1 Z2(A1u, xk1).
Finally, if C0 is sufficiently large, we can thus show from (9) that, with high probability,
kxk+1 x⇤k  1 2
⇣ r + 30 · 22 p n0/ndk✏k ⌘ .

heading: 4 Numerical Experiments
text: In this section, we provide numerical experiments on synthetic random data and MNIST images where the oberservations follow the model in (1). The first experiment illustrates a phase portrait that verifies Theorem 1 in the noiseless case. The second experiment illustrates that the PLUGIn algorithm (5) is stable to small dense noise and compares its performance to gradient descent. The third experiment illustrates decay rate of recovery error and reconstruction error of PLUGIn in the presence of noise and compares its performance to gradient descent. Lastly, the final experiments illustrates that PLUGIn can effectively remove dense unstructured random noise introduced to MNIST images. All experiments were conducted using Google Colaboratory.
In the synthetic experiments, we let the generative prior be a 3-layer neural network G : Rn0 ! Rn3 of the form G(z) = (A3 (A2 (A1x))), where the entries of weight matrix Ai 2 Rni⇥ni 1 are sampled from N (0, 1/ni). We sample the target latent code x⇤ uniformly from Sn0 1, set the noise level as ↵ 2 R, and set the noise to be ↵⌫ where ⌫ is sampled uniformly from Sn3 1. Then we set y = G(x⇤) + ↵⌫. We run PLUGIn and gradient descent each for 10,000 iterations, or until the difference in norm between consecutive iterates is less than 10 13 times the norm of the current iterate, whichever comes first, and set x̂ to be the output. We use a fixed step size of ⌘ = 3 for PLUGIn; for gradient descent, we use a parameter of 1000 to move along the gradient descent direction computed using PyTorch [26].
For the first experiment, we fix n0 2 {1, 2, . . . , 19}, n1 2 {1, 21, 41, . . . , 441}, n3 = n2 = n1, and ↵ = 0. We randomly initialize PLUGIn, setting x0 to have independent standard normal entries. Note that the norm of x0 concentrates around p n0, and thus is generally far from the norm of the target latent vector. For each trial, we say PLUGIn successfully recovers the target latent code if
kx̂ x⇤k  10 5. Figure 1 shows the fraction of successful recoveries from 20 independent trials using the from observation as described above. Black squares corresponds to no successful recovery and white squares correspond to 100% successful recovery.
For the second experiment, we fix n0 = 20, n1 = 600, n2 = 500, n3 = 500, and sample the noise level ↵ uniformly in the interval [0, 1]. In figures 2a and 2b, the solid line corresponds to the performance of PLUGIn and the dotted line represents the performance of gradient descent. Figure 2a shows the empirical dependence of the the relative recovery error kx̂ x⇤k/kx⇤k on the noise-to-signal ratio given, given by ↵, from 20 independent trials. Similarly, figure 2b shows the empirical dependence of the the relative reconstruction error kG(x̂) G(x⇤)k/kG(x⇤)k from 20 independent trials. The figures show that PLUGIn can be used to stably denoise signals using a generative prior with contractive layers.
For the third synthetic experiment, we fix n0 = 20, noise level ↵ = 0.1 and sample n1 = n2 = n3 in the interval [200, 5000]. Figures 2c and 2d show the empirical dependence of relative recovery error kx̂ x⇤k/kx⇤k and relative reconstruction error kG(x̂) G(x⇤)k/kG(x⇤)k, respectively, on the expansion ratio n3/n0 from 20 independent trials. In figures 2c and 2d, the solid line corresponds to PLUGIn and the dotted line corresponds to gradient descent. The figures confirm the linear dependence in the recovery and reconstruction errors of PLUGIn, in log scale, with respect to nd/n0, provided ni, for all i 1, are sufficiently large compared to n0. We now empirically show that PLUGIn can effectively remove noise synthetically introduced to MNIST images and compare its performance to gradient descent. We trained a VAE [27] using Adam optimizer [28] with a learning rate of 0.001 and mini-batch size 100 on the MNIST dataset consisting of 60,000 number of 28⇥ 28 images of handwritten digits [29]. The latent code dimension of the trained network G is 20 and the decoder network in the VAE is a fully connected network with parameters 20 500 500 784. Finally, a noisy image y is generated via pixelwise addition of an image y⇤ from the MNIST database and an noise vector ⌫ 2 R784, i.e. y = y⇤ + ⌫. In all MNIST
experiments, we use a fixed step size of ⌘ = 1/ for PLUGIn, where is the product of the operator norms of the weight matrices; for gradient descent, we use a parameter of 1000 to move along the gradient descent direction computed using PyTorch. Similar to the synthetic experiment, we run PLUGIn and gradient descent each for 10,000 iterations, or until the difference in norm between consecutive iterations is less than 10 13 times the norm of the current iterate, whichever comes first.
In figure 3, the images in the top row are the observations, the images in the second row and third row are the recovered images corresponding to PLUGIn and gradient descent, respectively. Figure 3 shows the result of using PLUGIn and gradient descent to remove noise ⌫ sampled uniformly over the sphere S783 from observation y. The noise level ↵ was adjusted so that the signal-to-noise ratio is in the interval [50, 200], with 50 corresponding to the leftmost column and 200 corresponding to the rightmost one.

heading: 5 Final Remarks
text: Extension to positively homogeneous activation functions
Our proofs and results can be extended to generative model G with positively homogeneous activation functions3 (such as Leaky ReLU). To see this, we modify the definition for Z(u, v;w) in Lemma 1 to
Z(u, v;w) := hAu, (Av) (Aw)i hu, v wi ,
where := E g · (g) with g ⇠ N (0, 1), and change the step size from ⌘ = 2d to ⌘ = d. The random process Z(u, v;w) remains centred (Appendix A, Lemma 2) and one can prove a similar version of Lemma 1. Moreover, in the multi-layer case, the range at each layer is still contained in a union of n0-dimensional (affine) subspaces in Rni (Appendix C, Lemma 7), so one can similarly prove a version of Theorem 1 in this case.

heading: Limitations
text: In Theorem 1 we assumed that weight matrices are independent and Gaussian. Although this is a commonly used assumption in literature (e.g. [10, 11, 12, 14, 17, 19, 23]) for theoretical analysis, practitioners may wonder whether theory based upon this assumption is instructive towards predicting convergence when the algorithm is applied to generative models with learned weights. Our experiments with MNIST give some evidence that the algorithm does indeed generalize. Without the Gaussian assumption, one needs a heuristic for choosing the step size. We find that the inverse of the product of the norm of the weight matrices works well. Secondly, in numerical experiments, comparison between PLUGIn and gradient descent over synthetic data seems to suggest that PLUGIn has error with similar dependence on parameters to gradient descent, but with a larger constant.
3Activation functions satisfies (↵t) = ↵ (t) for all ↵ 0 and t 2 R.

heading: Acknowledgments and Disclosure of Funding
text: Y. Plan is partially supported by an NSERC Discovery Grant (GR009284), an NSERC Discovery Accelerator Supplement (GR007657), and a Tier II Canada Research Chair in Data Science (GR009243). O. Yılmaz is partially supported by an NSERC Discovery Grant (22R82411) and Pacific Institute for the Mathematical Sciences (PIMS) CRG 33: High-Dimensional Data Analysis. B. Joshi is partially supported by the Pacific Institute for the Mathematical Sciences (PIMS). (The research and findings may not reflect those of the Institute.)

Review summaries:
This paper proposes an algorithm to solve the problem of inverting random generative models from noisy measurements where the weights of the model are drawn from $\mathcal{N}(0, 1/n)$ where n is the width of the corresponding layer, and the activation function is a ReLU unit.   Their proposed algorithm "Partially Linearized Update for Generative Inversion" (PLUGIn) is a simple iterative algorithm where each iterate is given by $x^{k+1} = x^k - \eta A_1^T \dots A_d^T (\mathcal{G}(x^k) - y)$. The authors demonstrate that   1. PLUGIn requires 2 blackboxes, (1) $\mathcal{G}(x^k)$ for every $x^k$ and (2) $ A_1^T \dots A_d^T $.  2. They show that for generative priors with possibly contractive layers, PLUGIn can effectively estimate the unknown signal from noisy observations, at a geometric rate.  3. In the presence of noise independent of the weight vector, the reconstruction error corresponding to $x^k$ for large enough $k$ goes as $\sqrt{4^d n_0/n_d}\| \epsilon \|$.   The proof of the result relies on simple properties of Gaussian matrices. In the case of the first layer it is enough to see that $E[A^T \sigma(Ax)] = 1/2 x$ and so with $\eta = 2$ we are done, for subsequent layers one needs to account for bias that creeps in and increase $\eta$ layer-wise. The authors overcome this by developing a series of uniform deviation inequalities, and bounding the resulting error. 

This paper introduces an algorithm for recovering an unknown latent code vector under a specific generative model: a deep generative network with ReLU activation functions. Innovation of paper compared with [10] is to show that one drops the $D_j$ and still recover an accurate estimate of $x^*$. Theoretical results on geometrical convergence to a neighborhood of $x^{*}$ are confirmed by experiments on synthetic random data and MNIST images.

This paper studies the question of inverting a generative model; formally, given relu network G and a noisy observation y = G(x^*) + eps, approximately recover x. Similar to existing works in this literature, they assume the weight matrices in G are all Gaussian, suitably scaled. But whereas previous works needed to assume expansivity, namely that at each layer the ratio between output and input dimensions is at least some universal constant (or more), this paper only requires that the ratio between the output dimension of any layer L and the input dimension of *the network itself* scales at least exponentially in L. In particular, they can handle networks where some of the intermediate layers are contractive, as long as the network is expansive "on average" across layers.  In this setting, they analyze the update rule x' <- x - eta M^T(G(x) - y), where M is the product of the weight matrices, and show that it converges exponentially quickly to a neighborhood of x^*. Their assumption on eps is very mild: it just needs to be independent of the randomness of the weight matrices.  The motivation for this update rule is that when eps = 0, after a single step from an arbitrary initialization, the next iterate is equal to x^* in expectation with respect to the randomness of the weight matrices. This is not true for subsequent iterates as they depend on this randomness, but one can nevertheless look at the norm of the distance between any iterate and x^*. To break the dependence, the author(s) upper bound this by the sum of suprema of certain random processes corresponding to the different layers of the network. They give a high probability bound for these suprema via a chaining argument. The key fact that lets them avoid the need for expansivity is that these processes are over points that lie in a union of subspaces whose dimension is given by the *input dimension of the network*.