Provable Compressed Sensing with Generative Priors via Langevin Dynamics

We focus on the problem of recovering a signal $x^{*}\in\mathbb{R}^{n}$ from a set of linear measurements $y\in\mathbb{R}^{m}$ where

To keep our analysis and results simple, we consider zero measurement noise, i.e., $\varepsilon=0$²²2We not in passing that our analysis techniques succeed for any vector $\varepsilon$ with bounded $\ell_{2}$ norm.. Here, $A\in\mathbb{R}^{m\times n}$ is a matrix populated with i.i.d. Gaussian entries with mean 0 and variance $1/m$. We assume that $x^{*}$ belongs to the range of a known generative model $G\mathrel{\mathop{\mathchar 58\relax}}\mathcal{D}\subset\mathbb{R}^{d}\rightarrow\mathbb{R}^{n}$; that is,

Following Bora et al. (2017b), we restrict $z$ to belong to a $d$-dimensional Euclidean ball, i.e., $\mathcal{D}=\mathcal{B}(0,R)$. Then, given the measurements $y$, our goal is to recover $x^{*}$. Again following Bora et al. (2017b), we do so by solving the usual optimization problem:

Hereon and otherwise stated, $\|\cdot\|$ denotes the $\ell_{2}$-norm. The most popular approach to solving (3.1) is to use gradient descent Bora et al. (2017b). For generative models $G(z)$ defined by deep neural networks, the function $F(z)$ is highly non-convex, and as such, it is impossible to guarantee global signal recovery using regular (projected) gradient descent.

We adopt a slightly more refined approach. Starting from an initial point $z_{0}\sim\mu_{0}$, our algorithm computes stochastic gradient updates of the form:

where $\xi_{k}$ is a unit Gaussian random vector in $\mathbb{R}^{d}$, $\eta$ is the step size and $\beta$ is an inverse temperature parameter. This update rule is known as stochastic gradient Langevin dynamics (SGLD) Welling & Teh (2011) and has been widely studied both in theory and practice Raginsky et al. (2017); Zhang et al. (2017). Intuitively, (3.2) is an Euler discretization of the continuous-time diffusion equation:

where $Z(0)\sim\mu_{0}$. Under standard regularity conditions on $F(z)$, one can show that the above diffusion has a unique invariant Gibbs measure.

We refine the standard SGLD to account for the boundedness of $z$. Specifically, we require an additional Metropolis-like accept/reject step to ensure that $z_{k+1}$ always belongs to the support $\mathcal{D}$, and also is not too far from $z_{k}$ of the previous iteration. We study this variant for theoretical analysis; in practice we have found that this is not necessary. Algorithm 1 (CS-SGLD) shows the detailed algorithm. Note that we can use stochastic (mini-batch) gradient instead of the full gradient $\nabla_{z}F(z)$.

We wish to derive sufficient conditions on the convergence (in distribution) of the random process in Algorithm 1 to the target distribution $\pi$, denoted by:

and study its consequence in recovering the true signal $x^{*}$. This leads to the first guarantees of a stochastic gradient-like method for compressed sensing with generative priors. In order to do so, we make the following three assumptions on the generator network $G(z)$.

1. (A.1)

   Boundedness. For all $z\in\mathcal{D}$, we have that $\|G(z)\|\leq B$ for some $B>0$.

2. (A.2)

   Near-isometry. $G(z)$ is a near-isometric mapping if there exist $0<\iota_{G}\leq\kappa_{G}$ such that the following holds for any $z,z^{\prime}\in\mathcal{D}$:

   $\displaystyle\iota_{G}\|z-z^{\prime}\|\leq\|G(z)-G(z^{\prime})\|\leq\kappa_{G}\|z-z^{\prime}\|.$

3. (A.3)

   Lipschitz gradients. The Jacobian of $G(z)$ is $M$-Lipschitz, i.e., for any $z,z^{\prime}\in\mathcal{D}$, we have

   $$\|\nabla_{z}G(z)-\nabla_{z}G(z^{\prime})\|\leq M\|z-z^{\prime}\|,$$

   where $\nabla_{z}G(z)=\frac{\partial G(z)}{\partial z}$ is the Jacobian of the mapping $G(\cdot)$ with respect to $z$.

All three assumptions are justifiable. Assumption (A.1) is reasonable due to the bounded domain $K$ and for well-trained generative models $G(z)$ whose target data distribution is normalized. Assumption (A.2) is reminiscent of the ubiquitous restricted isometry property (RIP) used for compressed sensing analysis (Candes & Tao, 2005) and is recently adopted in (Latorre et al., 2019). Finally, Assumption (A.3) is needed so that the loss function $F(z)$ is smooth, following typical analyses of Markov processes.

Next, we introduce a new concept of smoothness for generative networks. This concept is a weaker version of a condition on $G(\cdot)$ introduced in Latorre et al. (2019).

Following Latorre et al. (2019) (Assumption 2), we call this property “strong smoothness”. However, our definition of strong smoothness requires two parameters instead of one, and is weaker since we allow for an additive slack parameter $\gamma\geq 0$.

Definition 3.1 can be closely linked to the following property of the loss function $F(z)$ that turns out to be crucial in establishing convergence results for CS-SGLD.

It is straightforward to see that (3.6) essentially recovers the strong smoothness condition (3.5) if the measurement matrix $A$ is assumed to be the identity matrix. In compressed sensing, it is often the case that $A$ is a (sub)Gaussian matrix and that given a sufficient number of measurements as well as Assumptions (A.1), (A.2) and (A.3), the dissipativity of $F(z)$ for such an $A$ can still be established.

Once $F$ is shown to be dissipative, the machinery of Raginsky et al. (2017); Zhang et al. (2017); Zou et al. (2020) can be adapted to show that the convergence of CS-SGLD. The majority of the remainder of the paper is devoted to proving this series of technical claims.

We first show that a very broad class of generator networks satisfies the assumptions made above. The following proposition is an extension of a result in Latorre et al. (2019).

This proposition merits a thorough discussion. First, architectures with increasing layer sizes are common; many generative models (such as GANs) assume architectures of this sort. Observe that the non-decreasing layer size condition is much milder than the expansivity ratios of successive layers assumed in related work Hand & Voroninski (2018); Asim et al. (2019).

Second, the compactness assumption of the domain of $G$ is mild, and traces its provenance to earlier related works Bora et al. (2017b); Latorre et al. (2019). Moreover, common empirical techniques for training generative models (such as GANs) indeed assume that the latent vectors $z$ lie on the surface of a sphere White (2016).

Third, common activation functions such as the sigmoid, or the Exponential Linear Unit (ELU) are continuously differentiable and monotonic. Note that the standard Rectified Linear Unit (ReLU) activation does not satisfy these conditions, and establishing similar results for ReLU networks is deferred to future work.

The key for our theoretical analysis, as discussed above, is Definition 3.1, and establishing this requires Proposition 3.1. Interestingly however, in Section 5 below we provide empirical evidence that strong smoothness holds for generative adversarial networks with ReLU activation trained on the MNIST and CIFAR-10 image datasets.

We now obtain a measurement complexity result by deriving a bound on the number of measurements required for $F$ to be dissipative.

The above result can be derived using covering number arguments, similar to the treatment in Bora et al. (2017b). Observe that the number of measurements scales linearly with the dimension of the latent vector $z$ instead of the ambient dimension, keeping in line with the flavor of results in standard compressed sensing. Recent lower bounds reported Liu & Scarlett (2020) also have shown that the scaling of $m$ with respect to $d$ and $\log L$ might be tight for compressed sensing recovery in several natural parameter regimes.

We need two more quantities to readily state our convergence guarantee. Both definitions are widely used in the convergence analysis of MCMC methods. The first quantity defines the goodness of an initial distribution $\mu_{0}$ with respect to the target distribution $\pi$.

The next quantity is the Cheeger constant that connects the geometry of the objective function and the hitting time of SGLD to a particular set in the domain Zhang et al. (2017).

Putting all the above ingredients together, our main theoretical result describing the convergence of Algorithm 1 (CS-SGLD) for compressed sensing recovery is given as follows.

In words, if we choose a high enough inverse temperature and appropriate step size, CS-SGLD converges (in expectation) to a signal estimate with very low loss within a polynomial number of iterations.

Let us parse the above result further. First, observe that the right hand side of (3.7) consists of two terms. The first term can be made arbitrarily small (at the cost of greater computational cost since $\eta$ decreases ). The second term represents the irreducible expected error of the exact sampling algorithm on the Gibbs measure $\pi({\mathrm{d}}z)$, which is worse than the optimal loss obtained at $z=z^{*}$.

Second, suppose the right hand side of (3.7) is upper bounded by $\epsilon^{\prime}$. Once SGLD finds an $\epsilon^{\prime}$-approximate minimizer of the loss, in the regime of sufficient compressed sensing measurements (as specified by Lemma 3.1), we can invoke Theorem 1.1 in Bora et al. (2017b) along with Jensen’s inequality to immediately obtain a recovery guarantee, i.e.,

Third, the convergence rate of CS-SGLD can be slow. In particular, SGLD may require a polynomial number of iterations to recover the true signal, while linearized ADMM Latorre et al. (2019) converges within a logarithmic number of iterations up to a neighborhood of the true signal. Obtaining an improved characterization of CS-SGLD convergence (or perhaps devising a new linearly convergent algorithm) is an important direction for future work.

Fourth, the above result is for noiseless measurements. A rather similar result can be derived with noisy measurements of bounded noise (says, $\|\varepsilon\|\leq\sigma)$. This quantity (times a constant depending on $A$) will affect (3.7) up to an additive term that scales with $\sigma$. This is precisely in line with most compressed sensing recovery results and for simplicity we omit such a derivation.

Let $g(z)=\nabla_{z}F(z)$, $u$ and $w$ respectively be the points before and after one iteration of Algorithm 1; the Markov chain is written as $u\rightarrow v\rightarrow w$, where $v\sim\mathcal{N}(u-\eta g(u),\frac{2\eta}{\beta}I)$ with the following density:

Without the correction step, $P(v|u)$ is exactly the transition probability of the standard Langevin dynamics. Note also that one can construct a similar density with a stochastic (mini-batch) gradient. The process of $v\rightarrow w$ is

Let $p(u)=\mathbb{P}_{v\sim P(\cdot|u)}[v\in\mathcal{B}(u,r)\cap\mathcal{D}]$ be the probability of accepting $v$. The conditional density $Q(w|u)$ is

where $\delta_{u}(\cdot)$ is the Dirac-delta function at $u$. Similar to Zou et al. (2020); Zhang et al. (2017), we consider the $1/2$-lazy version of the above Markov process, with the transition distribution

and construct an auxiliary Markov process by adding an extra Metropolis accept/reject step. While proving the ergodicity of the Markov process with transition distribution $\mathcal{T}_{u}(w)$ is difficult, the auxiliary chain does indeed converge to a unique stationary distribution $\pi\propto e^{-\beta F(z)}\cdot{\mathbf{1}}(z\in\mathcal{D})$ due to the Metropolis-Hastings correction step.

The auxiliary Markov chain is given as follows: starting from $u$, let $w$ be the state generated from $\mathcal{T}_{u}(\cdot)$. The Metropolis-Hasting SGLD accepts $w$ with probability,

Let $\mathcal{T}^{\star}_{u}(\cdot)$ denote the transition distribution of the auxiliary Markov process, such that

Below, we establish the connection between $\mathcal{T}_{u}(\cdot)$ and $\mathcal{T}^{\star}_{u}(\cdot)$, as well as the convergence of the original chain in Algorithm 1 through a conductance analysis on $\mathcal{T}^{\star}_{u}(\cdot)$.

In Appendix B, we show that $F(z)$ is $L$-smooth with $L=(MB+\kappa_{G}^{2})$ and its gradient is bounded by $D=\kappa_{G}^{2}\|A^{\top}A\|$.

One can verify that $\mathcal{T}^{\star}_{u}(\cdot)$ is time-reversible (Zhang et al., 2017). Moreover, following Lovász et al. (1993); Lovász & Vempala (2007), the convergence of a time-reversible Markov chain to its stationary distribution depends on its conductance, which is defined as follows:

Using the conductance parameter $\phi$ and the closeness $\delta$ between $\mathcal{T}_{u}(\cdot)$ and $\mathcal{T}^{\star}_{u}(\cdot)$, we can derive the convergence of $\mathcal{T}_{u}(\cdot)$ in total variation distance.

We will further give a lower bound on $\delta$ in order to establish an explicit convergence rate.

Armed with these tools, we formally establish the first step of the proof.

The convergence rate is polynomial in the Cheeger constant $\rho$ whose lower bound is difficult to obtain generally. A rough bound $\rho=e^{-\tilde{O}(d)}$ can be derived using the Poincaré constant of the distribution $\pi$, under the smoothness assumption. See (Bakry et al., 2008) for details.

Apart from the step size and the number of iterations, the convergence depends on $\lambda$, the goodness of the initial distribution $\mu_{0}$. In this part, we specify a particular choice of $\mu_{0}$ in establish this.

As mentioned above, our theory relies on the assumption that the following condition holds for some constants $\alpha>0,\gamma\geq 0$ and $\forall z,z^{\prime}\in\mathcal{D}$ for a domain $\mathcal{D}$. Here, we take $\mathcal{D}=\mathbb{R}^{d}$.

To estimate these constants, we generate samples $z$ and $z^{\prime}$ from $\mathcal{N}(0,\mathbb{I})$. To establish $\alpha$ and $\gamma$, we perform experiments on two different datasets (i) MNIST (Net1) and (ii) CIFAR10 (Net2). For both datasets, we compute the terms $u(z,z^{\prime})=\nabla_{z}G(z)^{\top}(G(z)-G(z^{\prime})),z-z^{\prime}\rangle$ and $v(z,z^{\prime})=\|z-z^{\prime}\|^{2}$ for 500 different instantiations of $z$ and $z^{\prime}$. We then plot these pairs of $(\alpha v-\gamma,u)$ samples for different $z$’s and $z^{\prime}$’s and compute the values of $\alpha$ and $\gamma$ by a simple linear program. We do this experiment for two DCGAN generators trained on MNIST (Figure 1 (a)) as well as on CIFAR10 (Figure 1 (c)).

Similarly for the compressed sensing case, we also derive values $\alpha_{A}$ and $\gamma_{A}$, where a compressive matrix $A$ acts on the output of the generator $G$. Here, we have picked the number of measurements $m=0.1n$ where $n$ is the signal dimension. This is encapsulated in the following equation:

for all possible Gaussian matrices $A$ and different instantiations of $z$ and $z^{\prime}$. Here, we capture the left side of the inequality in $u(z,z^{\prime})=\langle\nabla_{z}(AG(z))^{\top}(AG(z)-AG(z^{\prime})),z-z^{\prime}\rangle$. We similarly plot points $(\alpha_{A}v(z,z^{\prime})-\gamma_{A},u(z,z^{\prime}))$ for all pairs $(z,z^{\prime})$. The scatter plot is generated for 500 different instantiations of $z$ and $z^{\prime}$ and $5$ different instantiations of $A$. We do this experiment for two DCGAN generators, one trained on MNIST (Figure 1 (b)) and the other trained on CIFAR10 (Figure 1 (d)). These experiments indicate that the dissipativity constant $\alpha$ is positive in all cases.

We test the SGLD reconstruction by using the update rule in (3.2) and compare against the optimizing the updates of $z$ using standard gradient descent as in Bora et al. (2017b). For all experiments, we use a pre-trained DCGAN generator, with network configuration described as follows: the generator consists of four different layers consisting of transposed convolutions, batch normalization and RELU activation; this is followed by a final layer with a transposed convolution and $\tanh$ activation Radford et al. (2015).

We display the reconstructions on MNIST in Figure 2. Note that the implementation in Bora et al. (2017b) requires 10 random restarts for CS reconstruction and they report the results corresponding to the best reconstruction. This likely suggests that the standard implementation is likely to get stuck in bad local minima or saddle points. For the sake of fair comparison, we fix the same random initialization of latent vector $z$ for both GD and SGLD with no restarts. We select $m=0.2n$. In Figure 2 we show reconstructions for the 16 different examples, which were all reconstructed at once using same $k=2000$ steps, learning rate of $\eta=0.02$ and the inverse temperature $\beta=1$ for both approaches. The only difference is the additional noise term in SGLD (Figure 2 part (d)). Notice that this additional noise component helps achieve better reconstruction performance overall as compared to simple gradient descent.

Phase transition plots scanning a range of compression ratios $m/n$ as well as example reconstructions on CIFAR-10 images can be found in the supplement. More thorough empirical comparisons with PGD-based approaches Shah & Hegde (2018); Raj et al. (2019) are deferred to future work.

We display the reconstructions on CIFAR10 in Figure 3. As with the implementation for MNIST, for the sake of fair comparison, we fix the same random initialization of latent vector $z$ for both GD and SGLD with no restarts. We select $m=0.3n$. In Figure 3 we show reconstructions for the 16 different examples from MNIST, which were all reconstructed at once using same $k=2000$ steps, learning rate of $\eta=0.05$ and the inverse temperature $\beta=1$ for both approaches. The only difference is the additional noise term in SGLD (Figure 2 part (d)). Similar to our experiments on MNIST we notice that this additional noise component helps achieve better reconstruction performance overall as compared to simple gradient descent.

Next, we plot phase transition diagrams by scanning the compression ratio $f=m/n=[0.2,0.4,0.6,0.8,1.0]$ for the MNIST dataset in Figure 4. For this experiment, we have chosen 5 different instantiations of the sampling matrix $A$ for each compression ratio $f$. In Figure 4 we report the average Mean Square Error (MSE) of reconstruction $\|\hat{x}-x\|^{2}$ over 5 different instances of $A$.

We conclude that SGLD gives improved reconstruction quality as compared to GD.