{'fraction_non_alphanumeric': 0.04095946485623003, 'fraction_numerical': 0.011905950479233226, 'mean_word_length': 4.465893587994543, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Provably Convergent Algorithms for Solving Inverse Problems Using Generative Models

The traditional approach of hand-crafting priors (such as sparsity) for solving inverse problems is slowly being replaced by the use of richer learned priors (such as those modeled by deep generative networks). In this work, we study the algorithmic aspects of such a learning-based approach from a theoretical perspective. For certain generative network architectures, we establish a simple non-convex algorithmic approach that (a) theoretically enjoys linear convergence guarantees for certain linear and nonlinear inverse problems, and (b) empirically improves upon conventional techniques such as back-propagation. We support our claims with the experimental results for solving various inverse problems. We also propose an extension of our approach that can handle model mismatch (i.e., situations where the generative network prior is not exactly applicable). Together, our contributions serve as building blocks towards a principled use of generative models in inverse problems with more complete algorithmic understanding.


Motivation
Inverse problems arise in a diverse range of application domains including computational imaging, optics, astrophysics, and seismic geo-exploration. In each of these applications, there is a target signal or image (or some other quantity of interest) to be obtained; a device (or some other physical process) records measurements of the target; and the goal is to reconstruct an estimate of the signal from the observations. Let us suppose that x * ∈ R n denotes the signal of interest and y = A(x * ) ∈ R m denotes the observed measurements. The aim is to recover (an estimate of) the unknown signal x * given y and A. Based on the forward measurement operator A, the inverse problem can be defined in two broad categories of linear and nonlinear problems. Many important problems in signal and image processing can be modeled with a linear measurement operator A; examples include compressive sensing, the classical problem of super-resolution or the problem of image inpainting. In case of nonlinear inverse problems, the operator A exhibits a nonlinearity; examples include phase retrieval, blind deconvolution, and de-quantization.
When m < n, the inverse problem is ill-posed, and some kind of prior (or regularizer) is necessary to obtain a meaningful solution. A common technique used to solve ill-posed inverse problems is to seek the minimum of a where z is a k-dimensional latent parameter vector and G is parameterized by the weights and biases of a d-layer neural network. One way to obtain such a model is to train a generative adversarial network [12]. Generative models have found remarkable applications in image analysis [13][14][15][16], and a well-trained generative model closely captures the notion of a signal (or image) being 'natural' [17]. Indeed, generative neural network learning algorithms have been successfully employed to solve inverse problems such as image super-resolution and inpainting [18,19]. However, most of these approaches are heuristic and provable characterization of such algorithms are not readily available.

Contributions
Our goal in this paper is to take some initial steps towards a principled use of generative priors for inverse problems by a) proposing and analyzing the well known projected gradient descent (PGD) algorithm for solving (1) for both linear and nonlinear inverse problems; b) building a general theoretical framework for analyzing performance of such approaches from an algorithmic standpoint. Specifically, apart from providing algorithms to solve inverse problems using generative network models, we also wish to understand the algorithmic costs involved with such algorithms: how computationally challenging they are, whether they provably succeed, and how to make such models robust.
The starting point of our work is the seminal paper by [20], who study the benefits of using generative models in the context of compressive sensing. In this paper, the authors pose the estimated target as the solution to a non-convex optimization problem and establish upper bounds on the statistical complexity of obtaining a "good enough" solution. Specifically, they prove that if the generative network is a mapping G : R k → R n simulated by a d-layer neural network with width ≤ n and with activation functions obeying certain properties, then m = O(kd log n) random observations are sufficient to obtain a good enough reconstruction estimate. However, they do not explicitly discuss an algorithm to perform such non-convex optimization. Moreover, the authors do not study the algorithmic costs of solving the optimization problem, and standard results in non-convex optimization are sufficient to only obtain sublinear convergence rates. In this work, we make several advances towards understanding the convergence properties of gradient descent (and related algorithms).
First: we establish a projected gradient descent (PGD) algorithm with linear convergence rates for the compressive sensing setup identical to [20], and demonstrate its empirical benefits over previous work. This constitutes Contribution I of this paper.
Second: we generalize this to a much wider range of nonlinear inverse problems. Using standard techniques, we propose a generic version of our PGD algorithm named ε−P GD for solving (1) where F is a smooth and strongly-convex objective function. We analyze this algorithm and prove its linear convergence under certain regimes of its smoothness and strong-convexity parameters. We also provide empirical results for solving nonlinear inverse problems. This forms Contribution II of this paper.
Third: we address the challenging inverse problem of phase retrieval [21,22]. This is similar to the compressed sensing setup described above, except with the additional difficulty that only the magnitude of the observations are available. We prove that a natural variation of PGD coupled with an intermediate phase estimation step converges (locally) linearly to the true solution. This forms Contribution III of this paper.  Figure 1: Illustration of our framework. Starting from a zero vector, we perform a gradient descent update step (red arrow) and projection step (blue arrow) alternatively.
Fourth: a drawback of [20] (and our contribution I) is the inability to deal with targets that are outside the range of the generative network model. This is not merely an artifact of their analysis; generative networks are rigid in the sense that once they are learned, they are incapable of reproducing any target outside their range. (This is in contrast with other popular parametric models such as sparsity models that exhibit a "graceful decay" property in the sense that if the sparsity parameter s is large enough, such models capture all possible points in the target space.) We address this gap, and propose an alternative algorithm using our general framework. We call it Myopic ε-PGD algorithm. It is novel, nonlinear extension of our previous work [23,24]. Under (fairly) standard assumptions, this algorithm also can be shown to demonstrate linear convergence. This constitutes Contribution IV of this paper.
In summary: we complement the work of [20] and [25] by providing PGD based algorithms for solving inverse problems, and algorithmic upper bounds for the corresponding problems studied in those works. Together, our contributions serve as further building blocks towards an algorithmic theory of generative models in both linear and nonlinear inverse problems.

Techniques
At a high level, our algorithms are standard. The primary novelty is in their applications to generative network models, and some aspects of their theoretical analysis.
Suppose that G : R k → R n is the generative network model under consideration. The cornerstone of our analysis is the assumption of an ε-approximate (Euclidean) projection oracle onto the range of G. We pre-suppose the availability of a computational routine P G that, given any vector x ∈ R n , can return a vector x ∈ Range(G) that approximately minimizes x − x 2 2 . The availability of this oracle, of course, depends on the nature of G. Some further comments on how to heuristically approximate this oracle are in Section 5.
For a special case of linear inverse problems (and compressive sensing in particular), we assume such oracle to be simply a gradient descent routine minimizing the x − x 2 2 over the latent variable z with x = G(z). Though this loss function is highly non-convex due to the presence of G, we find empirically that the gradient descent (implemented via back-propagation) works very well, and can be used as a projection oracle. Our procedure is depicted in Fig. 1. We choose a zero vector as our initial estimate (x 0 ), and in each iteration, we update our estimate by following the standard gradient descent update rule (red arrow in Fig. 1), followed by projection of the output onto the span of generator (G) (blue arrow in Fig. 1).
We support this specific PGD algorithm via a rigorous theoretical analysis. We show that the final estimate at the end of T iterations is an approximate reconstruction of the original signal x * , with very small reconstruction error; moreover, under certain sufficiency conditions on the linear operator A, PGD demonstrates linear convergence, meaning that T = log(1/δ) is sufficient to achieve δ-accuracy. Further, we present a series of numerical results as validation of our approach.
We also provide a direct generalization of the above approach for nonlinear inverse problems, that we call ε − P GD. We analyze this generic algorithm to show a linear convergence by assuming that the objective function in (1) obeys the Restricted Strong Convexity/Smoothness assumptions [26]. With this assumption, the proof of convergence follows from a straightforward modification of the proof given in [27]. Through our analysis, it indeed can be seen that the PGD algorithm for linear inverse problems is in fact a special case of ε-PGD.
The fourth algorithm (Myopic ε-PGD) is a novel approach for handling model mismatch in the target. The main idea (following the lead of [25]) is to pose the target x * as the superposition of two components: x * = G(z) + ν, where ν can be viewed as an "innovation" term that is s-sparse in some fixed, known basis B. The goal is now to recover both G(z) and ν. This is reminiscent of the problem of source separation or signal demixing [28], and in our previous work [24,29] we proposed greedy iterative algorithms for solving such demixing problems. We extend this work by proving a nonlinear extension, together with a new analysis, of the algorithm proposed in [24].

Inverse problems
The study of solving inverse problems has a long history. As discussed above, the general approach to solve an ill-posed inverse problem is to assumes that the target signal/image obeys a prior. Classical methods mainly used hand-crafted signal priors to distinguish 'natural' signals from the infinite set of feasible solutions. The prior can be encoded in the form of either a constraint set (as in Eq. (1)) or an extra regularization penalty. Several methods (including [30][31][32]) employ sparsity priors to solve linear inverse problems such as denoising, super-resolution, and inpainting. Despite their successful practical and theoretical results, all such hand-designed priors often fail to restrict the solution space only to natural images, and it is easily possible to generate signals satisfying the prior but do not resemble natural data.

Neural network models
The last few years have witnessed the emergence of trained neural networks for solving such problems. The main idea is to eschew hand-crafting any priors, and instead learn an end-to-end mapping from the measurement space to the image space. This mapping is simulated via a deep neural network, whose weights are learned from a large dataset of input-output training examples [33]. The works [34][35][36][37][38][39][40][41] have used this approach to solve several types of inverse problems, and has met with considerable success. However, the major limitations are that a new network has to be trained for each new linear inverse problem; moreover, most of these methods lack concrete theoretical guarantees. An exception of this line of work is the powerful framework of [42], which does not require retraining for each new problem; however, this too is not accompanied by theoretical analysis of statistical and computational costs.

Generative networks
A special class of neural networks that attempt to directly model the distribution of the input training samples are known as generative adversarial training networks, or GANs [12]. GANs have been shown to provide visually striking results [17,[43][44][45]. The use of GANs to solve linear inverse problems was advocated in [20]. Specifically, given (noisy) linear observations y = Ax * + e of a signal x * ∈ R n , assuming that x * belongs to the range of a generative network G : R n → R k , this approach constructs the reconstructed estimatex as follows: Gaussian measurements, then together with regularity assumptions on the generative network, they prove that the solutionx satisfies: for some constant C that can be reliably upper-bounded. In particular, in the absence of noise the recovery of x * is exact. However, there is no discussion of how computationally expensive this procedure is. Observe that the above minimization is highly non-convex (since for any reasonable neural network, G is a non-convex function) and possibly also non-smooth (if the neural network contains non-smooth activation functions, such as rectified linear units, or ReLUs). More recently, [46] improve upon the approach in [20] for solving more general nonlinear inverse problems (in particular, any inverse problem that has a computable derivative). Also, [47] have analyzed the above problem for the untrained case where the weights of the generative model G obeys certain randomness assumptions. See, also, [48].
Under similar statistical assumptions as [20], the work of [1] provably establishes a linear convergence rate, provided that a projection oracle (on to the range of G) is available, but only for the special case of compressive sensing. Our generalized result (Contribution II) extends this algorithm (and analysis) to more general nonlinear inverse problems.
More recently, [49] proposed a method that learns the a network-based projector for use in the PGD algorithm, making the projection step faster computationally. However, their theoretical result assumes the learned projector to be δ−approximate, indicating that the effective training of the projector is crucial for the success of their method posing an additional challenge.

Phase retrieval
The phase retrieval problem has been extensively studied over the last few decades [50][51][52] and it appears in several applications, including optical imaging [51,53], microscopy [54], and X-ray crystallography [55]. Phase retrieval is a non-convex problem and classical solution methods rely on alternating projection heuristics; examples include Gerchberg-Saxton [50] and Fienup [51]. In recent years, lifting-based methods were introduced that reformulate phase retrieval as a semidefinite program. [52]. Subsequently, non-convex methods have been proposed for solving phase retrieval problem with theoretical performance guarantees [56][57][58][59][60][61][62]. Most of the non-convex methods rely on estimating a good initial solution via the so-called spectral initialization method. A number of methods for solving phase retrieval using trained neural networks have been recently proposed [41,63,64].

Model mismatch
A limitation of most generative network models is that they can only reproduce estimates that are within their range; adding more observations or tweaking algorithmic parameters are completely ineffective if a generative network model is presented with a target that is far away from the range of the model. To resolve this type of model mismatch, the authors of [25] propose to model the signal x * as the superposition of two components: a "base" signal u = G(z), and an "innovation" signal v = Bν, where B is a known ortho-basis and ν is an l-sparse vector. In the context of compressive sensing, the authors of [25] solve a sparsity-regularized loss minimization problem: and prove that the reconstructed estimatex = G(ẑ) +v is close enough to x provided m = O((k + l)d log n) measurements are sufficient. However, as before, the algorithmic costs of solving the above problem are not discussed. Our third main result (Contribution III) proposes a new algorithm for dealing with model mismatches in generative network modeling, together with an analysis of its convergence and iteration complexity.

Main Algorithms and Analysis
Let us first establish some notational conventions. Below, · will denote the Euclidean norm unless explicitly specified. We use O(·)-notation in several places in order to avoid duplication of constants. We use F (·) to denote a (scalar) objective function.

Contribution I: Solving linear inverse problems
Let S ⊆ R n be the set of 'natural' images in data space with a vector x * ∈ S. We consider an ill-posed linear inverse problem (2) with the linear operator A(x) = Ax, where A is a Gaussian random matrix. For simplicity, we do not consider the additive noise term.
To solve for x (estimate of x * ), we choose Euclidean measurement error as the loss function F (·) in Eqn. (1). Therefore, given y and A, we seek

Algorithm
Our algorithm is described in Alg. 1. We assume that our trained generator network (G) well approximates the highdimensional probability distribution of the set S. With this assumption, we limit our search for x only to the range of Algorithm 1 PGD 1: Inputs: y, A, G, T , Output: t ← t + 1 7: x ← x T the generator function (G(z)). The function G is assumed to be differentiable, and hence we use back-propagation for calculating the gradients of the loss functions involving G for gradient descent updates.
The optimization problem in Eqn. 3 is similar to a least squares estimation problem, and a typical approach to solve such problems is to use gradient descent. However, the candidate solutions obtained after each gradient descent update need not represent a 'natural' image and may not belong to set S. We solve this limitation by projecting the candidate solution on the range of the generator function after each gradient descent update.
Thus, in each iteration of our proposed algorithm 1, two steps are performed in alternation: a gradient descent update step and a projection step. The first step is simply an application of a gradient descent update rule on the loss function F (·) with the learning rate η. In projection step, we minimize the projection loss by gradient descent updates with learning rate η in : Though the projection loss function is highly non-convex due to the presence of G, we find empirically that the gradient descent (implemented via back-propagation) works very well. Thus, the gradient descent based minimization serves as a projection oracle in this case. In each of the T iterations, we run T in gradient descent updates for calculating the projection. Therefore, T × T in is the total number of gradient descent updates required in our approach.

Analysis
Drawing parallels with standard compressive sensing theory, in our case, we need to ensure that the difference vector of any two signals in the set S lies away from the nullspace of the matrix A. This condition is encoded via the S-REC (Set Restricted Eigenvalue Condition) as defined and established in [20]. We slightly modify this condition and present it in the form of squared 2 -norm : Further, based on [65,66], we propose the following theorem about the convergence of our algorithm: Then, for every vector x * ∈ S, the sequence (x t ) defined by the algorithm PGD [1] exhibits linear convergence for a carefully chosen range of stepsizes η.
Proof: Suppose F (·) is the squared error loss function as defined above. Then, we have: Substituting y = Ax * and rearranging yields, Define: Then, by definition of the projection operator P G , the vector x t+1 is a better (or equally good) approximation to w as the true image x * . Therefore, we have: Substituting for w t and expanding both sides, we get: Substituting y = Ax * and rearranging yields, We now use 4 and 5 to obtain, Now, since A is a random Gaussian with sufficiently many rows, it satisfies the S-REC [20]: As x * , x t and x t+1 are 'natural' vectors, Substituting 7 in 6, Morever, as shown in [20], we have Av ≤ ρ v with high probability for any v. Therefore, we write: Let us choose learning rate(η) such that 1 2γ < η < 1 γ . We also have ρ 2 ≤ γ. Combining both, we get ρ 2 < 1 η , which makes the L.H.S. in the above equation negative. Therefore, where δ is inversely proportional to the number of measurements m [20]. If the slack parameter δ is small enough then it can be ignored. Also, 1 2γ < η < 1 γ yields, Hence, i.e., the sequence F (x t ) converges linearly upto a neighborhood of radius O(δ). In the noiseless case, this also implies that x t+1 converges to a small neighborhood around x * .

Contribution II: Solving nonlinear inverse problems
We now present the generic version of PGD algorithm suitable for large class of nonlinear inverse problems by generalizing for the loss function F (·) and the projection oracle P G .
We denote the F (·) to be a scalar function with continuous gradient, and assume that F has a continuous gradient which can be evaluated at any point x ∈ R n .
Recall that we wish to solve the problem where G is a generative network. To do so, we now employ a generalized version of projected gradient descent algorithm using the ε-approximate projection oracle for G. The algorithm is described in Alg. 2.
We define the ε-approximate projection oracle P G as, Definition 2 (Approximate projection) A function P G : R n → Range(G) is an ε-approximate projection oracle if for all x ∈ R n , P G (x) obeys: We will assume that for any given generative network G of interest, such a function P G exists and is computationally tractable 5 . Here, ε > 0 is a parameter that is known a priori.
In contrast to our previous analysis, here we introduce more general restriction conditions on the F (·): This assumption is by now standard; see [26,27] for in-depth discussions. This means that the objective function is strongly convex / strongly smooth along certain directions in the parameter space (in particular, those restricted to the set S of interest). The parameter α > 0 is called the restricted strong convexity (RSC) constant, while the parameter β > 0 is called the restricted strong smoothness (RSS) constant. Clearly, β ≥ α. In fact, throughout the paper, we assume that 1 ≤ β α < 2, which is a fairly stringent assumption but again, one that we do not know at the moment how to relax.
Definition 4 (Incoherence) A basis B and Range(G) are called µ-incoherent if for all u, u ∈ Range(G) and all v, v ∈ Span(B), we have: for some parameter 0 < µ < 1.
Remark 1 In addition to the above, we will make the following assumptions in order to aid the analysis. Below, γ and ∆ are positive constants.
Some comments about these assumptions may be warranted. The first says that the gradient of the loss at at the minimizer is small. The second says that the range of G is compact. The third links the parameters from the first two assumptions.
Proof: The proof is a minor modification of that in [27]. For simplicity we will assume that · refers to the Euclidean norm. Let us suppose that the step size η = 1 β . Define By invoking RSS, we get: where the last few steps are consequences of straightforward algebraic manipulation. Now, since x t+1 is an ε-approximate projection of w t onto Range(G) and x * ∈ Range(G), we have: Therefore, we get: However, due to RSC, we have: Therefore, where the last inequality follows from Cauchy-Schwartz and the assumptions on ∇F (x * ) and the diameter of Range(G). Further, by assumption, γ∆ ≤ O(ε). Rearranging terms, we get: for some constant C > 0.
This theorem asserts that the distance between the objective function at any iteration to the optimum decreases by a constant factor in every iteration. (The decay factor is β α − 1, which by assumption is a number between 0 and 1). Therefore, we immediately obtain linear convergence of ε-PGD up to a ball of radius O(ε): Therefore, the overall running time can be bounded as follows: It is noticeable that the analysis of PGD algorithm for linear problem is a special case of the generalized analysis given by Theorem 2. That is because once we set the F (·) as defined in Eq. (3), the RSC for F (·) can be obtained through the S-REC condition from Eq. (4). Similarly, we can use the upper bound of the spectral norm for the Gaussian matrix A to obtain RSS for F (·).
In Sec. 4, we provide empirical results for solving nonlinear inverse problems using Alg. 2. Specifically, we consider two nonlinear forward models: a sinusoidal model with A(x * ) = Ax * + sin(Ax * ); and a sigmoid model with A(x * ) = sigmoid(Ax * ) = 1 1+exp(−Ax * ) . While we use the L 2 -loss as a loss function in the case of sinusoidal model, for the sigmoid nonlinearity, we use a loss function specified as: where, Θ(·) is integral of A(·), and a i represents the rows of the measurement matrix A. The gradient of the loss can be calculated in closed form: Such choice of the loss function is inspired by the problem of single-index model (SIM) estimation [68]. [69] also advocates the usage of such a loss function.

Contribution III: Phase retrieval
We now propose a method to solve a different category of nonlinear inverse problems using our overall broad approach. Specifically, we will solve inverse problems where the observations are of the form: i.e., only magnitude information of the phaseless measurements are retained. This is the well-known phase retrieval setting [21,22]; here the added twist is that x is assumed to obey a generative prior.
One approach to solving this problem is to cast the recovery as a nonlinear optimization problem based on a suitable loss function F (x) that measures the (squared) difference between the measurements and the purported phaseless observations. We will slightly depart from this and instead propose a new algorithm (Phase-PGD) for solving the above problem. See Alg. 3. In each iteration of the Phase-PGD algorithm , three steps are performed: a phase update step, a gradient descent update step, and a projection step.
The first step is to calculate the phase of Ax. For real A and x, at the t th iteration, we update the phase estimate: After calculating the phase vector p, we can use an element-wise product between p and y as an estimate of linear measurements and convert the phase retrieval problem into a linear inverse problem.
The second step is simply an application of a gradient descent update rule on the loss function f (·) which is given as: Thus, the gradient descent update at the t th iteration is given by: where η is the learning rate.
The third steps is the projection step, in which we aim to find an image from the span of the generator, M which is closest to our current estimate w t . We define the projection operator P G as follows: where L in is the inner loss function defined as, We solve the inner optimization problem by running gradient descent with T in number of updates on L in (z).
We now analyze this algorithm. This part of the algorithm is described in Lines 3-7 of Algorithm 3. We can prove that provided a good initial estimate (x 0 ), Phase-PGD provably converges to x * . The high level intuition is as follows. Ignoring the noise, the observation model for phase retrieval can be restated as follows: for all i = {1, 2, . . . , m}. To ease notation, denote the phase vector p ∈ R m as a vector that contains the unknown signs of the measurements, i.e., p i = sign ( a i , x ) for all i = {1, 2, . . . , m}. Let p * denote the true phase vector and let P denote the set of all phase vectors, i.e. P = {p : p i = ±1, ∀i}. Then our measurement model gets modified as: Therefore, the recovery of x * can be posed as a (non-convex) optimization problem: To solve this problem, we alternate between estimating p and x. We perform two estimation steps: (a) if we fix the signal estimate x, then the minimizer p ∈ P is given in closed form as: (b) and if we fix the phase vector p, the signal vector x ∈ M can be obtained by solving: We now analyze our proposed descent scheme. We obtain: Theorem 3 Suppose we have an initialization x 0 ∈ M satisfying dist (x 0 , x * ) ≤ δ 0 x * 2 , for 0 < δ 0 < 1, and suppose the number of (Gaussian) measurements, m > C (kd log n) , for some large enough constant C. Then with high probability the iterates x t+1 of Algorithm 3, satisfy: where x t , x t+1 , x * ∈ M, and 0 < ρ < 1 is a constant.

Proof sketch:
The high level idea behind the proof is that with a δ-ball around the true signal x * , the "phase noise" can be suitably bounded in terms of a constant times the signal estimation error. To be more precise, suppose that z * = Ax * = p * y. Then, at any iteration t, we have: where e t can be viewed as the "phase noise". Now, examining Line 6 of the above algorithm, we have that x t is the output of Phase-PGD after t iterations. An "unpacking" argument similar to the one in [70,71] indicates that: where α is a small enough constant. We will show that e t can be also bounded in terms of x t−1 − x * , via Lemma 4 below. Consequently: where ρ is a small enough constant.
We therefore achieve a per-step error reduction scheme if the initial estimate x 0 satisfies x 0 − x * ≤ δ 0 x * . This result can be trivially extended to the case where the initial estimate x 0 satisfies x 0 + x * ≤ δ 0 x * , hence giving the convergence criterion of the form (for ρ < 1): We now state Lemma 4 without proof; this is a straightfoward adaptation of the covering number analysis of analysis of [20].
Theorem 4 Suppose that the generator network model G(·) is comprised of d layers of neurons with ReLU activation functions and weight matrices with bounded operator norms. As long as the initial estimate is a small distance away from the true signal x * ∈ M ( i.e. dist (x 0 , x * ) ≤ δ x * ) and subsequently, dist (x t , x * ) ≤ δ x * , where x t is the t th update of Algorithm 3, then the following bound holds for any t ≥ 0: with high probability, as long as m > C(kd log n) and ρ 1 < 1 is a constant.

Contribution IV: Addressing signal model mismatch
A key assumption in all the above discussion has been that the underlying signal is well-modeled by the generator G. This, of course, is idealistic; natural signals (and images) exhibit a wide variety of features, some of which may not be modeled by (even) a well-trained generator.
We now generalize the ε-PGD algorithm to handle situations involving signal model mismatch. Assume that the target signal can be decomposed as: where B T v 0 ≤ l n for some ortho-basis B.
Theorem 5 Let ⊕ denote the Minkowski sum. If F satisfies RSC/RSS over Range(G) ⊕ Span(B) with constants α and β, and if we assume µ-incoherence between B and Range(G), we have: Proof: We will generalize the proof technique of [24]. We first define some auxiliary variables that help us with the proof. Let: . and let x * = u * + v * be the minimizer that we seek. As above, by invoking RSS and with some algebra, we obtain: However, by definition, Therefore, . 6 The algorithm is a variant of block-coordinate descent, except that the block updates share the same gradient term.
But u t+1 is an ε-projection of w u t and u * is in the range of G, we have: Similarly, since v t+1 is an l-sparse thresholded version of w v t , we have: Plugging in these two upper bounds, we get: . Expanding squares and cancelling (several) terms, the right hand side of the above inequality can be simplified to obtain: Plugging this into (17), we get: We already know how to bound the first term T 1 , using an identical argument as in the proof of Theorem 2. We get: The second term T 2 can be bounded as follows. First, observe that This gives us the following inequalities: Similarly, Combining, we get: Moreover, by invoking RSC and Cauchy-Schwartz (similar to the proof of Theorem 2), we have: Therefore we obtain the upper bound on T 2 : Plugging in the upper bounds on T 1 and T 2 and re-arranging terms, we get: which leads to the desired result.

Experimental Results
Our primary focus in this paper is theoretical. However, we supplement our theory with representative numerical experiments that show the promise of our proposed algorithms. in this section, we describe our experimental setup and report algorithm performance.

Experiments with compressed sensing
We use two different GAN architectures and two different datasets in our experiments to show that our approach [Alg. 1] can work with variety of GAN architectures and datasets. We provide comparisons with the CSGM algorithm proposed in [20] as well as the Lasso.
In our experiments, we choose the entries of the matrix A independently from a Gaussian distribution with zero mean and variance 1/m. In these experiments, we ignore the presence of noise; however, our experiments can be replicated in the presence of additive Gaussian noise. For comparison with CSGM [20], we use a gradient descent optimizer keeping the total number of update steps (T × T in ) fixed for both algorithms to produce fair comparisons.
In the first experiment, we use a very simple GAN model trained on the MNIST dataset, which is collection of 60, 000 handwritten digit images, each of size 28 × 28 [72]. In our GAN, both the generator and the discriminator are fullyconnected neural networks with only one hidden layer. The generator consists of 20 input neurons, 200 hidden-layer neurons and 784 output neurons, while the discriminator consists of 784 input neurons, 128 hidden layer neurons and 1 output neuron. The size of the latent space is set to k = 20, i.e., the input to our generator is a standard normal vector z ∈ R 20 . We train the GAN using the method described in [12]. We use the Adam optimizer [73] with learning rate 0.001 and mini-batch size 128 for the training.
We test the MNIST GAN with 10 test images taken from the span of the generator to avoid model mismatch issues, and provide both quantitative and qualitative results. For PGD-GAN, because of the zero initialization, a high learning rate is required to get a meaningful output before passing it to the projection step. Therefore, we choose η = 0.5. The parameter η in is set to 0.01 with T = 15 and T in = 200. Thus, the total number of update steps is fixed to 3000. Similarly, the algorithm of [20] is tested with 3000 updates and η = 0.01. For reporting purposes, we use the reconstruction error = x − x * 2 . In Fig. 2(a), we show the reconstruction error comparisons for increasing values of number of measurements. It can be seen that our algorithm performs better than the other two methods. Also, as the input images are chosen from the span of the generator itself, it is possible to get close to zero error with only 100 measurements. Fig. 2(b) depicts reconstruction results for selected MNIST images.
The second set of our experiments are performed on a Deep Convolutional GAN (DCGAN) trained on the celebA dataset, which contains more than 200, 000 face images of celebrities [74]. We use a pre-trained DCGAN model, which was made available by [20]. The dimension of the latent space for the DCGAN model is k = 100. We report the results on a held out test dataset, unseen by the GAN at the time of training. The total number of updates is set to 1000, with T = 10 and T in = 100. Learning rates for PGD-GAN are set as η = 0.5 and η in = 0.1. The algorithm of [20] is run with η = 0.1 and 1000 update steps. Image reconstruction results from m = 1000 measurements with our algorithm are displayed in Fig. 2(c). We observe that our algorithm produces better reconstructions compared to the other baselines.

Experiments with nonlinear inverse problems
We extend our experiments to nonlinear models to depict the performance of our algorithm as described in Sec. 3.2. We present image reconstructions from the measurements obtained using two non-linear forward models: a sinusoidal model with A(x * ) = Ax * + sin(Ax * ), and a sigmoidal model with A(x * ) = sigmoid(Ax * ). Similar to the linear case, these experiments are performed using a DCGAN trained on celebA. It is evident that our algorithm produces superior reconstructions as depicted in Fig. 3.

Experiments with phase retrieval
In this section, we describe our experimental setup and report the performance comparisons of our proposed Phase-PGD and deep phase retrieval (DPR) method proposed in [75]. The DPR method estimates the signal as which requires solving an optimization problem directly over the latent space z of the generative model G(·).
In our experiments, we choose the entries of the matrix A independently from the N (0, 1 m ) distribution. Although we ignore the presence of noise, it is possible to replicate our experiments with additive Gaussian noise.
We use two different generative models for the MNIST and CelebA datasets. The generative model for CelebA follows the DCGAN framework [76] except that we do not use any batchnorm layer since the gradient for this layer is dependent on batch size and the distribution of the batch. We train our generators by jointly optimizing generator parameters and the latent code z using SGD by following the procedure in [77]. We use the squared-loss function, x − x 2 to train the generators. We choose z from the standard normal distribution on R k and normalize it to unit norm. We project z back to the unit norm after each gradient update.
In our first set of experiments, we use a generator trained over the MNIST training dataset resized to 32 × 32 pixel. We test two approaches on 10 images from the test set of MNIST dataset and provide both quantitative and qualitative results. For Phase-PGD, we choose gradient descent step size η = 0.9 and T = 50. For fair comparison, we use 2500 iterations for DPR. We show the reconstruction error comparison in Fig. 4b and SSIM comparison in Fig. 4c for increasing number of measurements. Since the input images are not chosen from the span of the generator, the reconstruction error does not necessarily reduce to zero. Nevertheless, we observe in Fig. 4b that Phase-PGD error gets close to zero with nearly 60 measurements, which is significantly smaller than those of DPR method. Fig. 4a depicts reconstruction results for some of the selected MNIST images for the two approaches.
For our second set of experiments, we train a generator for the CelebA dataset. For training, we resize 202,599 color images of celebrity faces in celebA dataset to 64 × 64 × 3 and kept 1 32 of the images for testing. We do not use the aligned and cropped version, which includes only the faces in the images.
We experiment on a subset of 10 images from the test dataset and report reconstruction results. We set the total number of updates to 1500 with T = 50 for Phase-PGD. Image reconstruction results from m = 1000 measurements using Phase-PGD algorithm are presented in Fig. 5a. A comparison between our method and DPR in terms of reconstruction error and SSIM is shown in Fig. 5b and Fig. 5c, respectively. We observe that Phase-PGD can achieve good reconstruction with significantly fewer measurements compared to DPR method.

DISCUSSION
Our contributions in this paper are primarily theoretical. We also explored the practical benefits of our approach in the context of inverse problems such as compressive sensing and phase retrieval. The algorithms proposed in this paper are generic and can be potentially used to solve other nonlinear inverse problems as well.