Projected Gradient Descent Algorithms for Solving Nonlinear Inverse Problems with Generative Priors

1 Introduction

Over the past two decades, the theoretical and algorithmic aspects of high-dimensional linear inverse problems have been studied extensively. The standard compressed sensing (CS) problem, which models low-dimensional structure via the sparsity assumption, is particularly well-understood [Foucart and Rauhut, 2013].

Despite the popularity of linear CS, in many real-world applications, nonlinearities may arise naturally, and it is more desirable to adopt nonlinear measurement models. For example, the semi-parametric single index model (SIM), which is formulated below, is a popular nonlinear measurement model that has long been studied [Han, 1987]:

∗Corresponding author.

where x∗ ∈ Rp is an unknown signal that is close to some structured set K, ai ∈ Rp are the sensing vectors, and fi : R → R are i.i.d. realizations of an unknown (possibly random) function f. In general, f plays the role of a nonlinearity, and it is called a link function. The goal is to estimate x∗ despite this unknown link function. Note that since the norm of x∗ can be absorbed into the unknown f, the signal x∗ is typically assumed to have unit ℓ2-norm.

In addition, inspired by the tremendous success of deep generative models in numerous real-world applications, recently, for the CS problem, it has been of interest to replace the sparsity assumption with the generative model assumption. More specifically, instead of being assumed to be sparse, the signal is assumed to lie near the range of a generative model, typically corresponding to a deep neural network [Bora et al., 2017]. Along with several theoretical developments, the authors of [Bora et al., 2017] perform extensive numerical experiments on image datasets to demonstrate that for a given accuracy, generative priors can reduce the required number of measurements by a factor of 5 to 10. There are a variety of follow-up works of [Bora et al., 2017], including [Heckel and Hand, 2019; Ongie et al., 2020; Jalal et al., 2021], among others.

In this paper, following the developments in both sparsitybased nonlinear inverse problems and inverse problems with generative priors, we provide theoretical guarantees for projected gradient descent (PGD) algorithms devised for nonlinear inverse problems using generative models.

1.1 RelatedWork

The most relevant existing works can roughly be divided into (i) nonlinear inverse problems without generative priors, and (ii) inverse problems with generative priors.

Nonlinear inverse problems without generative priors: The SIM has long been studied in the low-dimensional setting where p ≪ n, based on various assumptions on the sensing vector or link function. For example, the maximum rank correlation estimator has been proposed in [Han, 1987] under the assumption of a monotonic link function. In recent years, the SIM has also been studied in [Plan and Vershynin, 2016; Genzel, 2016; Plan et al., 2017; Oymak and Soltanolkotabi, 2017] in the high-dimensional setting where an accurate estimate can be obtained when n ≪ p, with the sensing vectors being assumed to be Gaussian. In particular, the authors of [Plan and Vershynin, 2016] show that the generalized Lasso approach works for high-dimensional SIM under the assumption that the set of structured signal K is convex, which is in general not satisfied for the range of a generative model with the Lipschitz continuity.

Nonetheless, the generality of the unknown link function in SIM comes at a price. Specifically, as mentioned above, it is necessary for the works studying SIM to assume the distribution of the sensing vector to be Gaussian or symmetric elliptical, and for nonlinear signal estimation problems with general sensing vectors, in order to achieve consistent estimation, knowledge of the link function is required [Zhang et al., 2018]. In addition, when the link function is unknown, since the norm of the signal may be absorbed into this link function, there is an identifiability issue and we are only able to estimate the direction of the signal. In practice, this can be unsatisfactory and may lead to large estimation errors. Moreover, for an unknown nonlinearity, it remains an open problem to handle signals with representation error [Plan and Vershynin, 2016], i.e., the signal (up to a fixed scale factor) is not exactly contained in K.

Based on these issues of SIM and some applications in machine learning such as the activation functions of deep neural networks [Yang et al., 2016], nonlinear measurement models with known and monotonic link functions 1 have been studied in [Yang et al., 2016; Zhang et al., 2018; Soltani and Hegde, 2017]. In [Yang et al., 2016], an ℓ1regularized nonlinear least-squares estimator is proposed, and an iterative soft thresholding algorithm is provided to efficiently approximate this estimator. The authors of [Zhang et al., 2018] propose an iterative hard thresholding (IHT) algorithm to minimize a nonlinear least-squares loss function subject to a combinatorial constraint. In the work [Soltani and Hegde, 2017], the demixing problem is formulated as minimizing a special (not the typical least-squares) loss function under a combinatorial constraint, and a corresponding IHT algorithm is designed to approximately find a minimizer. All the algorithms proposed in [Yang et al., 2016; Zhang et al., 2018; Soltani and Hegde, 2017] can be regarded as special cases of the PGD algorithm.

Inverse problems with generative priors: Bora et al. show that when the generative model is L-Lipschitz continuous with bounded k-dimensional inputs, roughly O(k logL) random Gaussian linear measurements are sufficient to attain accurate estimates [Bora et al., 2017]. Their analysis is based on minimizing a linear least-squares loss function, and the objective function is minimized directly over the latent variable in Rk using gradient descent. A PGD algorithm in the ambient space in Rp has been proposed in [Shah and Hegde, 2018; Peng et al., 2020] for noiseless and noisy Gaussian linear measurements respectively. It has been empirically demonstrated that this PGD algorithm leads to superior reconstruction performance over the algorithm used in [Bora et al., 2017]. Various nonlinear measurement models with known nonlinearity have also been studied for generative priors. Specifically, near-optimal sample complexity bounds for 1bit measurement models have been presented in [Qiu et al., 2020; Liu et al., 2020]. Furthermore, the works [Wei et al., 2019; Liu and Scarlett, 2020a] have provided near-optimal non-uniform recovery guarantees for nonlinear compressed sensing with an unknown nonlinearity. More specifically, the authors of [Wei et al., 2019] assume that the link function is differentiable and propose estimators via score functions based on the first and second order Steins identity. The differentiability assumption fails to hold for 1-bit and other quantized measurement models. To take such measurement models into consideration, the work [Liu and Scarlett, 2020a] instead makes the assumption that the (uncorrupted) observations are sub-Gaussian, and proposes to use a simple linear least-squares estimator despite the unknown nonlinearity. While obtaining these estimators is practically hard due to the typical non-convexity of the range of a generative model, both works are primarily theoretical, and no practical algorithm is provided to approximately find the estimators.

1 For this case, a natural idea is to apply approaches for linear measurement models to the inverted data {f -1 (yi), ai}ni=1. Unfortunately, such a simple idea works well only in the noiseless setting. See [Yang et al., 2016, Section 2] for a discussion.

1.2 Contributions

Throughout this paper, we make the assumption that the generative model is L-Lipschitz continuous with bounded kdimensional inputs (see, e.g., [Bora et al., 2017]). The main contributions of this paper are as follows:

• For the scenario of SIMwith unknown nonlinearity, we assume that the sensing vector is Gaussian and the signal is exactly contained in the range of the generative model, and propose a PGD algorithm for a linear least-squares estimator. We show that roughly O(k logL) samples suffice to ensure that this PGD algorithm converges linearly and yields an estimator with optimal statistical rate, which is roughly of order √k logL/n. While this PGD algorithm is identical to the PGD algorithm for solving linear inverse problems using generative models as proposed in [Shah and Hegde, 2018; Peng et al., 2020], the corresponding analysis is significantly different since we consider the SIM with an unknown nonlinear function f, instead of the simple linear measurement model. Moreover, unlike [Shah and Hegde, 2018; Peng et al., 2020], we have provided a neat theoretical guarantee for choosing the step size.
• For the scenario where the link function is known, we make much weaker assumptions for sensing vectors and allow for representation error, i.e., the signal do not quite reside in the range of the generative model, and propose a nonlinear least-squares estimator. We prove that the estimator enjoys optimal statistical rate, and show that a corresponding PGD algorithm converges linearly to this estimator. To the best of our knowledge, the corresponding PGD algorithm (cf. (18)) is novel.
• We perform various numerical experiments on image datasets to back up our theoretical results.

Remark 1. Generative model based phase retrieval has been studied in [Hand et al., 2018; Jagatap and Hegde, 2019; Hyder et al., 2019; Liu et al., 2021b]. However, for the scenario of SIM with unknown nonlinearity, we follow the settings in [Liu and Scarlett, 2020a], and as mentioned therein, phase retrieval is beyond the scope of this setup. For the case of a known link function, phase retrievel is also not applicable since its corresponding nonlinear functions are not monotonic. Moreover, it is typically unavoidable for phase retrieval with generative priors to require the strong assumption about the existence of a good initial vector, whereas for the nonlinear function (whether it is unknown or known) and the corresponding PGD algorithm considered in our work, the initial vector can be arbitrary.

1.3 Notation

We use upper and lower case boldface letters to denote matrices and vectors respectively. For any positive integer N, we write [N] = {1, 2, . . . ,N} and we use IN to represent an identity matrix in RN×N . A generative model is a function G : D → Rp , with latent dimension k, ambient dimension p, and input domain D ⊆ Rk . We focus on the setting where k ≪ p. For a set S ⊆ Rk and a generative model G : Rk → Rp , we write G(S) = {G(z) : z ∈ S}. We use ‖X‖2→2 to denote the spectral norm of a matrix X. We define the ℓq-ball Bk q (r) := {z ∈ Rk : ‖z‖q ≤ r} for q ∈ [0,+∞]. The symbols c,C are absolute constants whose values may be different per appearance.

2 Preliminaries

We present the definition for a sub-Gaussian random variable. Definition 1. A random variable X is said to be subGaussian if there exists a positive constant C such that (E [|X|q ])1/q ≤ C√q for all q ≥ 1. The sub-Gaussian norm of a sub-Gaussian random variable X is defined as ‖X‖ψ2 := supq≥1 q-1/2 (E [|X|q ])1/q .

Throughout this paper, we make the assumption that the generative model G : Bk 2 (r) → Rp is L-Lipschitz continuous, and we fix the structured set K to be the range of G, i.e., K := G(Bk 2 (r)).

In the following, we state the definition of the Twosided Set-Restricted Eigenvalue Condition (TS-REC), which is adapted from the S-REC proposed in [Bora et al., 2017].

Definition 2. Let S ⊆ Rp . For parameters ϵ ∈ (0, 1), δ ≥ 0, a matrix ˜A ∈ Rn×p is said to satisfy the TS-REC(S, ϵ, δ) if, for every x1, x2 ∈ S, it holds that

Suppose that B ∈ Rn×p has i.i.d. N (0, 1) entries. We have the following lemma, which says that 1 √n B satisfies TS-REC for the set K = G(Bk 2 (r)) with high probability.

Lemma 1. (Adapted from [Bora et al., 2017, Lemma 4.1]) For ϵ ∈ (0, 1) and δ > 0, if n = Ω ( k ϵ2 log Lr δ ), 2 then a random matrix 1 √n B ∈ Rn×p with bij i.i.d. ∼ N (0, 1) satisfies the TS-REC(K, ϵ, δ) with probability 1 - e-Ω(ϵ2 n) .

2 Here and in subsequent statements of lemmas and theorems, the implied constant is assumed to be sufficiently large.

3 PGD for Unknown Nonlinearity

In this section, we provide theoretical guarantees for a PGD algorithm in the case that the nonlinear link function f is unknown. For this case, we follow [Liu and Scarlett, 2020a] and make the assumptions:

• Let x∗ ∈ Rp be the signal to estimate. We assume that µx∗ is contained in the set K = G(Bk 2 (r)), where 3

is a fixed parameter depending solely on f. Since the norm of x∗ may be absorbed into the unknown f, for simplicity of presentation, we assume that ∥∥x∗∥∥2 = 1.

• ai are i.i.d. realizations of a random vector a ∼ N(0, Ip), with a being independent of f. We write the sensing matrix as A = [a1, . . . , an]T ∈ Rn×p .
• We assume the SIM for the (unknown) uncorrupted measurements y1, y2, . . . , yn as in (1).
• Similarly to [Plan et al., 2017; Liu and Scarlett, 2020a], the random variable y := f(aT x∗ ) is assumed to be sub-Gaussian with sub-Gaussian norm ψ, i.e.,

where g ∼ N(0, 1). Such an assumption will be satisfied, e.g., when f does not grow faster than linearly, i.e., for any x ∈ R, |f(x)| ≤ a + b|x| for some a and b. Hence, various noisy 1-bit measurement models and non-binary quantization schemes satisfy this assumption [Liu and Scarlett, 2020a].

• In addition to possible random noise in f, we allow for adversarial noise that may depend on a. In particular, instead of observing y directly, we only assume access to (corrupted) measurements ˜y = [˜ y1, . . . , ˜ yn]T ∈ Rn satisfying

for some τ ≥ 0, where y = [y1, y2, . . . , yn]T ∈ Rn .

• To derive an estimate of the signal x∗ (up to constant scaling), we minimize the linear ℓ2 loss over K:

The above optimization problem is referred to as the generalized Lasso or K-Lasso. The idea behind using the K-Lasso to derive an accurate estimate even for nonlinear observations is that the nonlinearity is regarded as noise and the nonlinear observation model can be converted into a scaled linear model with unconventional noise [Plan and Vershynin, 2016].

The authors of [Liu and Scarlett, 2020a] provide recovery guarantees with respect to globally optimal solutions of (6), but they have not designed practical algorithms to find an optimal solution. Solving (6) may be practically difficult since

3µ is important for analyzing the recovery performance, but note that the knowledge of µ cannot be assumed since f is unknown.

Algorithm 1 A PGD algorithm for approximately solving (6) (PGD-GLasso)

Input: A, ˜y, ν > 0, number of iterations T , generative model

G, arbitrary initial vector x(0)

Procedure: Iterate as in (8) for t = 0, . . . , T -1; return x(T)

in general, K = G(Bk 2 (r)) is not a convex set. In this section, we consider using the following iterative procedure to approximately solve (6):

where PK(·) is the projection function onto K and ν > 0 is a tuning parameter. For convenience, the corresponding algorithm is described in Algorithm 1.

Remark 2. We will implicitly assume the exact projection in analysis. Our proof technique does not require PK to be unique, but only requires it to be a retraction onto the manifold of the generative prior. The exact projection assumption is also made in relevant works including [Hyder et al., 2019; Liu et al., 2022; Peng et al., 2020; Shah and Hegde, 2018]. In practice approximate methods might be needed, and both gradient-based projection [Shah and Hegde, 2018] and GANbased projection [Raj et al., 2019] have been shown to be highly effective. Compared to exact projection, in practice global optima of optimization problems like (6) are typically much more difficult to approximate, and projectionbased methods may serve as powerful tools for approximating the global optima. For example, for the simple linear measurement model, it has been numerically verified that performing gradient descent over the latent variable (without using projection) leads to inferior performance and cannot approximate the global optima of (6) well, whereas a projection-based gradient descent method gives better reconstruction [Shah and Hegde, 2018].

Algorithm 1 is identical to the PGD algorithm for solving linear inverse problems using generative models as proposed in [Shah and Hegde, 2018; Peng et al., 2020]. However, the corresponding analysis is significantly different, since we consider the SIMwith an unknown nonlinear function f, instead of the simple linear measurement model. In particular, we have the following theorem showing that if 2µ1 < 1, Algorithm 1 converges linearly and achieves optimal statistical rate, which is roughly of order √k logL/n. The proof of Theorem 1 is provided in the supplementary material.

Theorem 1. Recall that µ and ψ are defined in (3) and (4) respectively. For any ϵ ∈ (0, 1), letting

For any δ > 0 satisfying Lr = Ω(δp), if n = Ω( k ϵ2 log Lr δ ) and 2µ1 < 1 with 1 - 2µ1 = Θ(1), then for any t ∈ N, with probability 1 - e-Ω(ϵ2 n) , we have

To ensure that 2µ1 < 1, if ϵ is chosen to be a sufficiently small positive constant, we should select the parameter ν from the interval (0.5, 1.5), and a good choice of ν is ν = 1. In addition, a d-layer neural network generative model typically has Lipschitz constant L = pΘ(d) [Bora et al., 2017], and thus we may set r = pΘ(d) and δ = 1 pΘ(d) without affecting the scaling of the term log Lr δ (and the assumption Lr = Ω(δp) is certainly satisfied for fixed (r, δ)). Hence, if there is no adversarial noise, i.e., τ = 0, and ψ is a fixed constant, we see that after a sufficient number of iterations, Algorithm 1 will return a point x(T) satisfying ∥∥x(T ) - µx∗∥∥2 = O(√ k log Lr δ n ). By the analysis of sample complexity lower bounds for noisy linear CS using generative models [Liu and Scarlett, 2020b; Kamath et al., 2020], this statistical rate is optimal and cannot be improved without extra assumptions.

4 PGD for Known Nonlinearity

In this section, we provide theoretical guarantees for the case when the nonlinear link function f is known. Throughout this section, we make the following assumptions:

• Unlike in the case of unknown nonlinearity, we now allow for representation error and assume that the signal x∗ lies near (but does not need to be exactly contained in) K = G(Bk 2 (r)). Note that since we have precise knowledge of f, we do not need to make any assumption on the norm of x∗ .
• The sensing matrix A ∈ Rn×p satisfies the following two assumptions with high probability:
• Johnson-Lindenstrauss embeddings (JLE): For any ϵ ∈ (0, 1) and any finite set E ⊆ Rp satisfying n = Ω( 1 ϵc1 · logc2 |E|) for some absolute constants c1, c2, we have for all x ∈ E that

1. Bounded spectral norm: For some absolute constant a, it holds that

When A has independent isotropic sub-Gaussian rows, 4 from Lemma [Vershynin, 2010, Proposition 5.16], we have that when n = Ω( 1 ϵ2 log |E|) (thus c1 = 2 and c2 = 1), the event corresponding to (11) occurs with probability 1 - e-Ω(nϵ2 ) . In addition, similarly to [Vershynin, 2010, Corollary 5.35], we have that with probability 1 - e-Ω(n) , the assumption about bounded spectral norm is satisfied with a = 0.5. Moreover, from the theoretical results concerning JLE in [Krahmer and Ward, 2011] and the standard inequality ‖A‖2→2 ≤ max{‖A‖1→1, ‖A‖∞→∞}, we know that when A is a subsampled Fourier matrix or a partial Gaussian circulant matrix with random column sign flips and isotropic rows, these two assumptions are also satisfied with high probability for appropriate absolute constants c1, c2 and a. Hence, these assumptions on A are significantly more generalized than the i.i.d. Gaussian assumption made for the case of unknown nonlinearity. Notably, when the two assumptions are satisfied, by a chaining argument [Bora et al., 2017; Liu et al., 2021a], the random matrix 1 √n A satisfies the TS-REC for K.

4 A random vector v is said to be isotropic if E[vvT ] = Ip.

• The (unknown) uncorrupted measurements are generated from the following measurement model:

where f : R → R is a known (deterministic) nonlinear function, and ηi are additive noise terms. Similarly to [Yang et al., 2016; Soltani and Hegde, 2017; Zhang et al., 2018], we assume that f is monotonic, differentiable, and for all x ∈ R, f ′ (x) ∈ [l, u] with u ≥ l > 0 being fixed constants. 5 In addition, we assume that ηi are independent realizations of zeromean sub-Gaussian random variables with maximum sub-Gaussian norm σ.

• We also allow for adversarial noise and assume that for some τ ≥ 0, the observed (corrupted) vector ˜y satisfies

• To estimate the signal x∗ , we utilize the knowledge of f , and consider minimizing the nonlinear ℓ2 loss over K:

Under the preceding assumptions, we have the following theorem that gives a recovery guarantee for optimal solutions to (15). The proof is placed in the supplementary material.

Theorem 2. Let ¯x = argminx∈K ‖x-x∗ ‖2. For any δ > 0, we have that any solution ˆx to (15) satisfies

In (16), the term ‖¯x - x∗ ‖2 corresponds to the representation error. Similarly to the discussion after Theorem 1, we see that when there is no representation error or adversarial noise, and considering δ being sufficiently small and σ being a fixed constant, we obtain the optimal statistical rate, i.e., ‖ˆx - x∗ ‖2 = O(√ k log Lr δ ).

Theorem 2 is concerned with globally optimal solutions of (15), which are intractable to obtain due to the nonconvexity of the corresponding objective function. In the

5 The case that l ≤ u < 0 can be similarly handled.

Algorithm 2 A PGD algorithm for approximately solving (15) (PGD-NLasso)

Input: A, ˜y, ζ > 0, number of iterations T , generative model G, arbitrary initial vector x(0)

Procedure: Iterate as (18) for t = 0, . . . , T - 1; return x(T)

following, we use a PGD algorithm to approximately minimize (15). In particular, we select an initial vector x(0) arbitrarily, and for any non-negative integer t, letting

where ζ > 0 is the step size and "⊙" represents elementwise product. The corresponding algorithm is described in Algorithm 2 for convenience.

Next, we present the following theorem, which establishes a theoretical guarantee similar to that of Theorem 1, except that there is an extra ‖¯x-x∗ ‖2 term corresponding to representation error in the upper bound. The proof of Theorem 3 can be found in the supplementary material.

Theorem 3. Let ¯x = argminx∈K ‖x-x∗ ‖2. For any δ > 0 and ϵ ∈ (0, 1) that is sufficiently small, letting

If 2µ2 < 1 with 1 - 2µ2 = Θ(1), we have for all t ∈ N that

5 Experiments

In this section, we empirically evaluate the performance of Algorithm 1 (PGD-GLasso; abbreviated to PGD-G) for the linear least-squares estimator (6) and Algorithm 2

(PGD-NLasso; abbreviated to PGD-N) for the nonlinear least-squares estimator (15) on the MNIST [LeCun et al., 1998] and CelebA [Liu et al., 2015] datasets. For both datasets, we use the generative models pre-trained by the authors of [Bora et al., 2017]. For PGD-G, the step size ν is set to be 1. For PGD-N, we set the step size ζ to be 0.2. For both algorithms, the total number of iterations T is set to be 30. On the MNIST dataset, we do 5 random restarts, and pick the best estimate among these random restarts. The generative model G is set to be a variational autoencoder (VAE) model with a latent dimension of k = 20. The projection step PK(·) with K being the range of G is approximated using gradient descent, performed using the Adam optimizer with a learning rate of 0.03 and 200 steps. On the CelebA dataset, we use a Deep Convolutional Generative Adversarial Networks (DCGAN) generative model with a latent dimension of k = 100. We select the best estimate among 2 random restarts. An Adam optimizer with 100 steps and a learning rate of 0.1 is used for the projection operator PK(·). Throughout this section, for simplicity, we consider the case that there is no adversarial noise, i.e., τ = 0. Since for unknown nonlinearity, the signal is recovered up to a scalar ambiguity, to compare performance across algorithms, we use a scale-invariant metric named Co-

where x∗ is the signal vector to estimate, and x(T) denotes the output vector of an algorithm. We use Python 2.7 and TensorFlow 1.0.1, with a NVIDIA Tesla K80 24GB GPU.

We follow the measurement model (13) to generate the observations, with f(x) = 2x+0.5 cos(x) and the noises being independent zero-mean Gaussian with standard deviation σ. We observe that f is smooth and monotonically increasing with l = 1.5 and u = 2.5. Thus, we perform PGD-N on the generated data. In addition, since the assumption made for unknown nonlinearity that f(g) for g ∼ N(0, 1) is subGaussian is satisfied, we also compare with PGD-G. The standard deviation σ is set to be 0.1 for the MNIST dataset, and 0.01 for the CelebA dataset. The baseline is the Lasso using 2D Discrete Cosine Transform (2D-DCT) basis [Tibshirani, 1996] and the method for linear inverse problem with generative models proposed in [Bora et al., 2017] (denoted by CSGM). For Lasso, CSGM and PGD-G, the sensing matrix A is assumed to contain i.i.d. standard Gaussian entries. Moreover, since for a known and monotonic nonlinear link function, we allow for a wide variety of distributions on the sensing vectors, for PGD-N, we consider both cases where A is a standard Gaussian matrix or a partial Gaussian circulant matrix similar to that in [Liu et al., 2021b]. The corresponding Algorithm 2 is denoted by PGD-N(G) or PGD-N(C) respectively, where "G" refers to the standard Gaussian matrix and "C" refers to the partial Gaussian circulant matrix.

We perform experiments to compare the performance of these algorithms, and the reconstructed results are reported in Figures 1, 2, and 3. We observe that for our settings, the sparsity-based method Lasso always attains poor reconstructions, while all three generative model based PGD methods attain accurate reconstructions even when the number of measurements n is small compared to the ambient dimension p. In addition, from Figure 3(a), we observe that on the MINST dataset, PGD-N(G) and PGD-N(C) lead to similar Cosine Similarity, and they are clearly better than CSGM and PGD-G when n < 300. This is not surprising since both CSGM and PGD-G does not make use of the knowledge of the nonlinear link function. From Figures 2 and 3(b), we see that for the CelebA dataset, three generative prior based PGD methods give similar reconstructed images, with the Cosine Similarity corresponding to PGD-N(G) and PGD-N(C) being slightly higher than that of PGD-G, and all of them lead to clearly higher Cosine Similarity compared to that of CSGM.

Numerical results for noisy 1-bit measurements are presented in the supplementary material.

6 Conclusion

We have proposed PGD algorithms to solve generative model based nonlinear inverse problems, and we have provided theoretical guarantees for these algorithms for both unknown and known link functions. In particular, these algorithms are guaranteed to converge linearly to points achieving optimal statistical rate in spite of the model nonlinearity.

Acknowledgments

We sincerely thank the three anonymous reviewers for their careful reading and insightful comments. We are extremely grateful to Dr. Jonathan Scarlett for proofreading the manuscript and giving valuable suggestions.