Title: COMPRESSED SENSING WITH DEEP IMAGE PRIOR AND LEARNED REGULARIZATION

Abstract: We propose a novel method for compressed sensing recovery using untrained deep generative models. Our method is based on the recently proposed Deep Image Prior (DIP), wherein the convolutional weights of the network are optimized to match the observed measurements. We show that this approach can be applied to solve any differentiable linear inverse problem, outperforming previous unlearned methods. Unlike various learned approaches based on generative models, our method does not require pre-training over large datasets. We further introduce a novel learned regularization technique, which incorporates prior information on the network weights. This reduces reconstruction error, especially for noisy measurements. Finally we prove that, using the DIP optimization approach, moderately overparameterized single-layer networks trained can perfectly fit any signal despite the nonconvex nature of the fitting problem. This theoretical result provides justification for early stopping.

1. INTRODUCTION:

We consider the well-studied compressed sensing problem of recovering an unknown signal x * ∈ R n by observing a set of noisy measurements y ∈ R m of the form y = Ax * + η.

(1)

Here A ∈ R m×n is a known measurement matrix, typically generated with random independent Gaussian entries. Since the number of measurements m is smaller than the dimension n of the unknown vector x * , this is an under-determined system of noisy linear equations and hence ill-posed. There are many solutions, and some structure must be assumed on x * to have any hope of recovery.

Pioneering research (Donoho, 2006; Candès et al., 2006; Candes & Tao, 2005) established that if x * is assumed to be sparse in a known basis, a small number of measurements will be provably sufficient to recover the unknown vector in polynomial time using methods such as Lasso (Tibshirani, 1996) .

Sparsity approaches have proven successful, but more complex models with additional structure have been recently proposed such as model-based compressive sensing (Baraniuk et al., 2010) and manifold models (Hegde et al., 2008; Hegde & Baraniuk, 2012; Eftekhari & Wakin, 2015) . Bora et al. (2017) showed that deep generative models can be used as excellent priors for images. They also showed that backpropagation can be used to solve the signal recovery problem by performing gradient descent in the generative latent space. This method enabled image generation with significantly fewer measurements compared to Lasso for a given reconstruction error. Compressed sensing using deep generative models was further improved in very recent work (Tripathi et al., 2018; Grover & Ermon, 2018a; Kabkab et al., 2018; Shah & Hegde, 2018; Fletcher & Rangan, 2017; Asim et al., 2018) . Additionally a theoretical analysis of the nonconvex gradient descent algorithm (Bora et al., 2017) was proposed by Hand & Voroninski (2017) under some assumptions on the generative model.

Inspired by these impressive benefits of deep generative models, we chose to investigate the application of such methods for medical imaging, a canonical application of compressive sensing. A significant problem, however, is that all these previous methods require the existence of pre-trained models. While this has been achieved for various types of images, e.g. human faces of CelebA (Liu et al., 2015) via DCGAN (Radford et al., 2015) , it remains significantly more challenging for med-ical images (Wolterink et al., 2017; Schlegl et al., 2017; Nie et al., 2017; Schlemper et al., 2017) . Instead of addressing this problem in generative models, we found an easier way to circumvent it.

Surprising recent work by Ulyanov et al. (2017) proposed Deep Image Prior (DIP), which uses untrained convolutional neural networks. In DIP-based schemes, a convolutional neural network generator (e.g. DCGAN) is initialized with random weights; these weights are subsequently optimized to make the network produce an output as close to the target image as possible. This procedure is unlearned, using no prior information from other images. The prior is enforced only by the fixed convolutional structure of the generator network.

Generators used for DIP are typically over-parameterized, i.e. the number of network weights is much larger compared to the output dimension. For this reason DIP has empirically been found to overfit to noise if run for too many iterations (Ulyanov et al., 2017) . In this paper we theoretically prove that this phenomenon occurs with gradient descent and justify the use of early stopping and other regularization methods.

Our Contributions:

• In Section 3 we propose DIP for compressed sensing (CS-DIP). Our basic method is as follows. Initialize a DCGAN generator with random weights; use gradient descent to optimize these weights such that the network produces an output which agrees with the observed measurements as much as possible. This unlearned method can be improved with a novel learned regularization technique, which regularizes the DCGAN weights throughout the optimization process.

• In Section 4 we theoretically prove that DIP will fit any signal to zero error with gradient descent. Our result is established for a network with a single hidden layer and sufficient constant fraction over-parametrization. While it is expected that over-parametrized neural networks can fit any signal, the fact that gradient descent can provably solve this nonconvex problem is interesting and provides theoretical justification for early stopping.

• In Section 5 we empirically show that CS-DIP outperforms previous unlearned methods in many cases. While pre-trained or "learned" methods frequently perform better (Bora et al., 2017) , we have the advantage of not requiring a generative model trained over large datasets. As such, we can apply our method to various medical imaging datasets for which data acquisition is expensive and generative models are difficult to train.

2. BACKGROUND:



2.1. COMPRESSED SENSING: CLASSICAL AND UNLEARNED APPROACHES:

A classical assumption made in compressed sensing is that the vector x * is k-sparse in some basis such as wavelet or discrete cosine transform (DCT). Finding the sparsest solution to an underdetermined linear system of equations is NP-hard in general; however, if the matrix A satisfies conditions such as the Restricted Eigenvalue Condition (REC) or Restricted Isometry Property (RIP) (Candes et al., 2006; Bickel et al., 2009; Donoho, 2006; Tibshirani, 1996) , then x * can be recovered in polynomial time via convex relaxations (Tropp, 2006) or iterative methods. There is extensive compressed sensing literature regarding assumptions on A, numerous recovery algorithms, and variations of RIP and REC (Bickel et al., 2009; Negahban et al., 2009; Agarwal et al., 2010; Bach et al., 2012; Loh & Wainwright, 2011) .

Compressed sensing methods have found many applications in imaging, for example the single-pixel camera (SPC) (Duarte et al., 2008) . Medical tomographic applications include x-ray radiography, microwave imaging, magnetic resonance imaging (MRI) (Winters et al., 2010; Chen et al., 2008; Lustig et al., 2007) . Obtaining measurements for medical imaging can be costly, time-consuming, and in some cases dangerous to the patient (Qaisar et al., 2013) . As such, an important goal is to reduce the number of measurements while maintaining good reconstruction quality.

Aside from the classical use of sparsity, recent work has used other priors to solve linear inverse problems. Plug-and-play priors (Venkatakrishnan et al., 2013; Chan et al., 2017) and Regularization by Denoising (Romano et al., 2017) have shown how image denoisers can be used to solve general linear inverse problems. A key example of this is BM3D-AMP, which applies a Block-Matching and 3D filtering (BM3D) denoiser to an Approximate Message Passing (D-AMP) algorithm (Metzler et al., 2016; 2015) . AMP has also been applied to linear models in other contexts (Schniter et al., 2016) . Another related algorithm is TVAL3 (Zhang et al., 2013; Li et al., 2009) which leverages augmented Lagrangian multipliers to achieve impressive performance on compressed sensing problems. In many different settings, we compare our algorithm to these prior methods: BM3D-AMP, TVAL3, and Lasso.

2.2. COMPRESSED SENSING: LEARNED APPROACHES:

While sparsity in some chosen basis is well-established, recent work has shown better empirical performance when neural networks are used (Bora et al., 2017) . This success is attributed to the fact that neural networks are capable of learning image priors from very large datasets (Goodfellow et al., 2014; Kingma & Welling, 2013) . There is significant recent work on solving linear inverse problems using various learned techniques, e.g. recurrent generative models (Mardani et al., 2017b) and auto-regressive models (Dave et al., 2018) . Additionally approximate message passing (AMP) has been extended to a learned setting by Metzler et al. (2017) . Bora et al. (2017) is the closest to our set-up. In this work the authors assume that the unknown signal is in the range of a pre-trained generative model such as a generative adversarial network (GAN) (Goodfellow et al., 2014) or variational autoencoder (VAE) (Kingma & Welling, 2013) . The recovery of the unknown signal is obtained via gradient descent in the latent space by searching for a signal that satisfies the measurements. This can be directly applied for linear inverse problems and more generally to any differentiable measurement process. Recent work has built upon these methods using new optimization techniques (Chang et al., 2017) , uncertainty autoencoders (Grover & Ermon, 2018b) , and other approaches (Dhar et al., 2018; Kabkab et al., 2018; Mixon & Villar, 2018; Pandit et al., 2019; Rusu et al., 2018; Hand et al., 2018) . The key point is that all this prior work requires pre-trained generative models, in contrast to CS-DIP. Finally, there is significant ongoing work to understand DIP and develop related approaches (Heckel et al., 2018; Dittmer et al., 2018) .

3. PROPOSED ALGORITHM:

Let x * ∈ R n be the signal that we are trying to reconstruct, A ∈ R m×n be the measurement matrix, and η ∈ R m be independent noise. Given the measurement matrix A and the observations y = Ax * + η, we wish to reconstruct an x that is close to x * .

A generative model is a deterministic function G(z; w): R k → R n which takes as input a seed z ∈ R k and is parameterized by "weights" w ∈ R d , producing an output G(z; w) ∈ R n . These models have shown excellent performance generating real-life signals such as images (Goodfellow et al., 2014; Kingma & Welling, 2013) and audio (Van Den Oord et al., 2016) . We investigate deep convolutional generative models, a special case in which the model architecture has multiple cascaded layers of convolutional filters (Krizhevsky et al., 2012) . In this paper we apply a DC-GAN (Radford et al., 2015) model and restrict the signals to be images.

3.1. COMPRESSED SENSING WITH DEEP IMAGE PRIOR (CS-DIP):

Our approach is to find a set of weights for the convolutional network such that the measurement matrix applied to the network output, i.e. AG(z; w), matches the measurements y we are given. Hence we initialize an untrained network G(z; w) with some fixed z and solve the following:

w * = arg min w y -AG(z; w) 2 . (2)

This is, of course, a non-convex problem because G(z; w) is a complex feed-forward neural network.

Still we can use gradient-based optimizers for any generative model and measurement process that is differentiable. Generator networks such as DCGAN are biased toward smooth, natural images due to their convolutional structure; thus the network structure alone provides a good prior for reconstructing images in problems such as inpainting and denoising (Ulyanov et al., 2017) . Our finding is that this applies to general linear measurement processes. Furthermore, our method also directly applies to any differentiable forward operator A. We restrict our solution to lie in the span of a convolutional neural network. If a sufficient number of measurements m is given, we obtain an output such that x * ≈ G(z; w * ).

Note that this method uses an untrained generative model and optimizes over the network weights w. In contrast previous methods, such as that of Bora et al. (2017) , use a trained model and optimize over the latent z-space, solving z * = arg min z y -AG(z; w) 2 . We instead initialize a random z with Gaussian i.i.d. entries and keep this fixed throughout the optimization process.

In our algorithm we leverage the well-established total variation regularization (Rudin et al., 1992; Wang et al., 2008; Liu et al., 2018) , denoted as T V (G(z; w)). We also propose an additional learned regularization technique, LR(w); note that without this technique, i.e. when λ L = 0, our method is completely unlearned. Lastly we use early stopping, a phenomena that will be analyzed theoretically in Section 4.

Thus the final optimization problem becomes

EQUATION

The regularization term contains hyperparameters λ T and λ L for total variation and learned regularization: R(w; λ T , λ L ) = λ T T V (G(z; w)) + λ L LR(w). Next we discuss this LR(w) term.

3.2. LEARNED REGULARIZATION:

Without learned regularization CS-DIP relies only on linear measurements taken from one unknown image. We now introduce a novel method which leverages a small amount of training data to optimize regularization. In this case training data refers to measurements from additional ground truth of a similar type, e.g. other x-ray images.

To leverage this additional information, we pose Eqn. 3 as a Maximum a Posteriori (MAP) estimation problem and propose a novel prior on the weights of the generative model. This prior then acts as a regularization term, penalizing the model toward an optimal set of weights w * .

For a set of weights w ∈ R d , we model the likelihood of the measurements y = Ax, y ∈ R m , and the prior on the weights w as Gaussian distributions given by

p(y|w) = exp - y -AG(z; w) 2 2λ L (2πλ L ) m ; p(w) = exp -1 2 (w -µ) T Σ -1 (w -µ) (2π) d |Σ| , where µ ∈ R d and Σ ∈ R d×d .

In this setting we want to find a set of weights w * that maximizes the posterior on w given y, i.e.,

w * = arg max w p(w|y) ≡ arg min w y -AG(z; w) 2 + λ L (w -µ) T Σ -1 (w -µ) . (4)

This gives us the learned regularization term

EQUATION

where the coefficient λ L in Eqn. 4 controls the strength of the prior.

Our motivation for assuming a Gaussian distribution on the weights is to build upon the proven success of 2 regularization, which also makes this assumption. Notice that when µ = 0 and Σ = I d×d , this regularization term is equivalent to 2 -regularization. Thus this method can be thought of as an adaptive version of standard weight decay. Further, because the network weights are initialized Gaussian i.i.d., we assumed the optimized weights would also be Gaussian. Previous work has shown evidence that the convolutional weights in a trained network do indeed follow a Gaussian distribution (Ma et al., 2018) .

3.2.1. LEARNING THE PRIOR PARAMETERS:

In the previous section, we introduced the learned regularization term LR(w) defined in Eqn. 5. However we have not yet learned values for parameters (µ, Σ) that incorporate prior knowledge of the network weights. We now propose a way to estimate these parameters.

Assume we have a set of measurements

S Y = {y 1 , y 2 , • • • , y Q } from Q different images S X = {x 1 , x 2 , • • • , x Q },

each obtained with a different measurement matrix A. For each measurement y q , q ∈ {1, 2, ..., Q}, we run CS-DIP to solve the optimization problem in Eqn. 3 and obtain an optimal set of weights

W * = {w * 1 , w * 2 , • • • , w * Q }.

Note that when optimizing for the weights W * , we only have access to the measurements S Y , not the ground truth S X .

The number of weights d in deep networks tends to be very large. As such, learning a distribution over each weight, i.e. estimating µ ∈ R d and Σ ∈ R d×d , becomes intractable. We instead use a layer-wise approach: with L network layers, we have µ ∈ R L and Σ ∈ R L×L . Thus each weight within layer l ∈ {1, 2, ..., L} is modeled according to the same N (µ l , Σ ll ) distribution. For simplicity we assume Σ ij = 0 ∀ i = j, i.e. that network weights are independent across layers. The process of estimating statistics (µ, Σ) from W * is described in Algorithm 1 of the appendix.

We use this learned (µ, Σ) in the regularization term LR(w) from Eqn. 5 for reconstructing measurements of images. We refer to this technique as learned regularization. While this may seem analogous to batch normalization (Ioffe & Szegedy, 2015) , note that we only use (µ, Σ) to penalize the 2 -norm of the weights and do not normalize the layer outputs themselves.

3.2.2. DISCUSSION OF LEARNED REGULARIZATION:

The proposed CS-DIP does not require training if no learned regularization is used, i.e. if λ L = 0 in Eqn. 3. This means that CS-DIP can be applied only with measurements from a single image and no prior information of similar images in a dataset.

Our next idea, learned regularization, utilizes a small amount of prior information, requiring access to measurements from a small number of similar images (roughly 5 -10). In contrast, other pretrained models such as that of Bora et al. (2017) require access to ground truth from a massive number of similar images (tens of thousands for CelebA). If such a large dataset is available, and if a good generative model can be trained on that dataset, we expect that pre-trained models would outperform our method. Our approach is instead more suitable for reconstructing problems where large amounts of data or good generative models are not readily available.

4. THEORETICAL RESULTS:

In this section we provide theoretical evidence to highlight the importance of early stopping for DIPbased approaches. Here we focus on denoising a noisy signal y ∈ R n by optimizing over network weights. This problem takes the form:

EQUATION

)

This is a special instance of Eqn. 2 with the measurement matrix A = I corresponding to denoising. We focus on generators consisting of a single hidden-layer ReLU network with k inputs, d hidden units, and n outputs. Using w = (W, V ) the generator model in this case is given by

EQUATION

where z ∈ R k is the input, W ∈ R d×k the input-to-hidden weights, and V ∈ R n×d the hiddento-output weights. We assume V is fixed at random and train over W using gradient descent. With these formulations in place, we are now ready to state our theoretical result. Theorem 4.1. Consider fitting a generator of the form

W → G(z; W, V ) = V • ReLU (W z) to a signal y ∈ R n with z ∈ R k , W ∈ R d×k , V ∈ R n×d ,

and ReLU (z) = max(0, z). Furthermore, assume V is a random matrix with i.i.d. N (0, ν 2 ) entries with ν = 1 √ dn y z . Starting from an initial weight matrix W 0 selected at random with i.i.d. N (0, 1) entries, we run gradient descent updates of the form W τ +1 = W τ -η∇L(W τ ) on the loss

L(W ) = 1 2 V • ReLU (W z) -y 2 ,

with step size η = η y 2 8n 4n+d where η ≤ 1. Assuming that d ≥ Cn, with C a fixed numerical constant, then

V • ReLU (W τ z) -y ≤ 3 1 - η 8(4n + d) τ y

holds for all τ with probability at least 1 -5e -n/2 -e -d/2 -e -4d 2 3 n 1 3 .

Our theoretical result shows that after many iterative updates, gradient descent will solve this nonconvex optimization problem and fit any signal y, if the generator network is sufficiently wide. This occurs as soon as the number of hidden units d exceeds the signal size n by a constant factor. In our theorem we have focused on the case where the measurement matrix is the identity (A = I). We note however that our theorem directly applies to many compressed sensing measurement matrices, in particular any matrix obtained by subsampling the rows of an orthonormal matrix (e.g. subsampling a Fourier matrix). This is possible as for any such matrix, AV has the same distribution as a Gaussian matrix with i.i.d. entries. Therefore, our theorem directly applies in this setting by replacing V with AV . This result demonstrates that early stopping is necessary for DIP-based methods to be successful; otherwise the network can fit any signal, including one that is noisy.

Our proof builds on theoretical ideas from Oymak & Soltanolkotabi (2019) which provide a general framework for establishing global convergence guarantees for overparameterized nonlinear learning problems based on various properties of the Jacobian mapping along the gradient descent trajectory. Related literature can be found in Du et al. ( 2018); Oymak & Soltanolkotabi (2018) . We combine tools from empirical process theory, random matrix theory, and matrix algebra to show that, starting from a random initialization, the Jacobian mapping across all iterates has favorable properties with high probability, hence facilitating convergence to a global optima.

5. EXPERIMENTS:



5.1. EXPERIMENTAL SETUP:

Measurements: We evaluate our algorithm using two different measurements processes, i.e. matrices A ∈ R m×n . First we set the entries of A to be Gaussian i.i.d. such that A i,j ∼ N (0, 1 m ). Recall m is the number of measurements, and n is the number of pixels in the ground truth image. This measurement process is standard practice in compressed sensing literature; hence we use it on each dataset. Additionally we use a Fourier measurement process common in MRI applications (Mardani et al., 2018; 2017a; Hammernik et al., 2018; Lehtinen et al., 2018; Lustig et al., 2008) .

Datasets: We use our algorithm to reconstruct both grayscale and RGB images. For grayscale we use the first 100 images in the test set of MNIST (LeCun et al., 1998) and also 60 random images from the Shenzhen Chest X-Ray Dataset (Jaeger et al., 2014) , downsampling a 512x512 crop to 256x256 pixels. For RGB we use retinopathy images from the STARE dataset (Hoover et al., 2000) with 512x512 crops downsized to 128x128 pixels. Baselines: We compare our algorithm to state-of-the-art unlearned methods such as BM3D-AMP (Metzler et al., 2016; 2015) , TVAL3 (Li, 2011; Li et al., 2009; Zhang et al., 2013) , and Lasso in a DCT basis (Ahmed et al., 1974) . We also evaluated the performance of Lasso in a Daubechies wavelet basis (Daubechies, 1988; Wasilewski, 2010) but found this performed worse than Lasso -DCT on all datasets. Thus hereon we refer to Lasso -DCT as "Lasso" and do not include results of Lasso -Wavelet.

Metrics: To quantitatively evaluate the performance of our algorithm, we use per-pixel meansquared error (MSE) between the reconstruction x and true image x * , i.e. x-x * 2 n . Note that because these pixels are over the range [-1, 1], it's possible for MSE to be greater than 1.

Implementation: To find a set of weights w * that minimize Eqn. 3, we use PyTorch (Paszke et al., 2017) with a DCGAN architecture. For baselines BM3D-AMP and TVAL3, we use the repositories provided by the authors Metzler et al. (2018) and Li et al., respectively. For baseline reconstructions Lasso, we use scikit-learn (Pedregosa et al., 2011) . Section A in the appendix provides further details on our experimental procedures, e.g. choosing hyperparameters. The supplementary material also contains our code repository for these experiments.

5.2. EXPERIMENTAL RESULTS:



5.2.1. RESULTS: LEARNED REGULARIZATION:

We first evaluate the benefits of learned regularization by comparing our algorithm with and without learned regularization, i.e. λ L = 100 and λ L = 0, respectively. The latter setting is an unlearned method, as we are not leveraging (µ, Σ) from a specific dataset. In the former setting we first learn (µ, Σ) from a particular set of x-ray images; we then evaluate on a different set of x-ray images. We compare these two settings with varying noise and across different number of measurements.

Our results in Table 1 show that learned regularization does indeed provide benefit. This benefit tends to increase with more noise or fewer measurements. Thus we can infer that assuming a learned Gaussian distribution over weights is useful, especially when the original signal is noisy or significantly compressed.

5.2.2. RESULTS: UNLEARNED CS-DIP:

For the remainder of this section, we evaluate our algorithm in the noiseless case without learned regularization, i.e. when η = 0 in Eqn. 1 and λ L = 0 in Eqn. 3. Hence CS-DIP is completely unlearned; as such, we compare it to other state-of-the-art unlearned algorithms on various datasets and with different measurement matrices. 

MNIST::

In Figure 1b we plot reconstruction error with varying number of measurements m of n = 784. This demonstrates that our algorithm outperforms baselines in almost all cases. Figure 2b shows reconstructions for 75 measurements, while remaining reconstructions are in the appendix.

Chest X-Rays: In Figure 1a we plot reconstruction error with varying number of measurements m of n = 65536. Figure 2a shows reconstructions for 2000 measurements; remaining reconstructions are in the appendix. On this dataset we outperform all baselines except BM3D-AMP for higher m. However for lower m, e.g. when the ratio m n ≤ 3%, BM3D-AMP often doesn't converge. This finding seems to support the work of Metzler et al. (2015) : BM3D-AMP performs well on higher m, e.g. m n ≥ 10%, but recovery at lower sampling rates is not demonstrated. Retinopathy: We plot reconstruction error with varying number of measurements m of n = 49152 in Figure 3a of the appendix. On this RGB dataset we quantitatively outperform all baselines except BM3D-AMP on higher m; however, even at these higher m, patches of green and purple pixels corrupt the image reconstructions as seen in Figure 10 . Similar to x-ray for lower m, BM3D-AMP often fails to produce anything sensible. All retinopathy reconstructions are located in the appendix.

Fourier Measurement Process: All previous experiments used a measurement matrix A containing Gaussian i.i.d. entries. We now consider the case where the measurement matrix is a subsampled Fourier matrix, as discussed in Section A of the appendix. Our algorithm outperforms baselines on the x-ray dataset, as shown in the appendix Figure 3b .

Additional Experiments, Appendix A: In Figure 4a we compare our algorithm against the learned approach of Bora et al. (2017) . In Figure 4b we compare our algorithm against baselines in the presence of noise, i.e. when η = 0 in Eqn. 1. We also demonstrate runtimes for all algorithms on the x-ray dataset in Table 2 .

6. CONCLUSION:

We demonstrate how Deep Image Prior (DIP) can be generalized to solve any differentiable linear inverse problem. We further propose a learned regularization method which, with a small amount of prior information, reduces reconstruction error for noisy or compressed measurements. Lastly we prove that the DIP optimization technique can fit any signal given a sufficiently wide single-layer network. This provides theoretical justification for regularization methods such as early stopping.

A ADDITIONAL EXPERIMENTS, DETAILS, AND INSIGHTS:

Hyperparameters: Our algorithm CS-DIP is implemented in PyTorch using the RMSProp optimizer (Tieleman & Hinton, 2012) with learning rate 10 -3 , momentum 0.9, and 1000 update steps for every set of measurements. These parameters are the same across all datasets. We initialize one random measurement matrix A for each image in all experiments. Further we set the TV hyperparameter λ T = 0.01, and found that high frequency components of the image were not reconstructed as sharply with λ T = 0.

Dataset-specific design choices: On larger images such as x-ray (n = 65536) and retinopathy (n = 49152), we found no difference using random restarts of the initial seed z. However for smaller vectors such as MNIST (n = 784), restarts did provide some benefit. As such our experiments utilize 5 random restarts for MNIST and one initial seed (no restarts) for x-ray and retinopathy images. For choosing hyperparameter λ L in Eqn. 3, we used a standard grid search and selected the best one. We used a similar grid search procedure for choosing dataset-specific hyperparameters in baseline algorithms BM3D-AMP, TVAL3, and Lasso.

Learned Regularization: For learned regularization, 10 x-ray images are used to learn µ and Σ. The results are then evaluated, i.e. averaged over x * , on 50 additional x-ray images.

Network Architecture: Our network has depth 7 and uses convolutional layers with ReLU activations. The initial seed z in Eqn. 3 is initialized with random Gaussian i.i.d. entries and then held fixed as we optimize over network weights w. We found negligible difference when varying the dimension of z (within reason), as this only affects the number of channels in the network's first layer. Hence we set the dimension of z to be 128, a standard choice for DCGAN architectures.

Stopping Criterion: We stop after 1000 iterations in all experiments. We did not optimize this hyperparameter, so admittedly there could be room for improvement. We further note that the "Error vs. Iterations" curve of CS-DIP with RMSProp did not monotonically decrease for some learning rates, even though error gradually decreased in all cases. As such we implemented a stopping condition which chooses the reconstruction with least error over the last 20 iterations. Note we choose this reconstruction based off measurement loss and do not look at the ground truth image.

Fourier Sampling: The measurements are obtained by sampling Fourier coefficients along radial lines as demonstrated in Figure 13 . For a 2D image x and a set of indices Ω, the measurements we receive are given by y (i,j) = [F(x)] (i,j) , (i, j) ∈ Ω, where F is the 2D Fourier transform. We choose Ω to be indices along radial lines, as shown in Figure 13 of the appendix; this choice of Ω is common in literature (Candès et al., 2006) and MRI applications (Mardani et al., 2017a; Lustig et al., 2008; Eksioglu & Tanc, 2018) . While Fourier subsampling is common in MRI applications, we use it here on images of x-rays simply to demonstrate that our algorithm performs well with different measurement processes.

We compare our algorithm to baselines on the x-ray dataset for {3, 5, 10, 20} radial lines in the Fourier domain, which corresponds to {381, 634, 1260, 2500} Fourier coefficients, respectively. We plot reconstruction error with varying number of Fourier coefficients in Figure 3b , outperforming baselines. Reconstructions be found in Figure 14 .

Table 2: Runtime (seconds) for each algorithm with varying number of measurements. Note that our algorithm was run on a NVIDIA GTX 1080-Ti GPU, while baselines only leverage CPU. Thus our goal here is not to compare runtimes; instead we demonstrate that our algorithm can run in a reasonable amount of time, which is an issue with other DIP methods. Unfortunately an RGB version of TVAL3 does not currently exist, although related TV algorithms such as FTVd perform similar denoising tasks (Wang et al., 2008) . Get v ∈ R S , a vector of S uniformly sampled weights from the l th layer of w * q 5:

M t [l, :] ← v T where M t [l, :] is the l th row of matrix M t ∈ R L×S 6: µ t [l] ← 1 S S i=1 v i 7:

end for 8:

Σ t ← 1 S M t M T t -µ t µ T t 9: end for 10: µ ← 1 T T t=1 µ t 11: Σ ← 1 T T t=1 Σ t

B PROOF OF SECTION 4: THEORETICAL JUSTIFICATION FOR EARLY STOPPING:

In this section we prove our theoretical result in Theorem 4.1. We begin with a summary of some notations we use throughout in Section B.1. Next, we state some preliminary calculations in Section B.2. Then, we state a few key lemmas in Section B.3 with the proofs deferred to Appendix C. Finally, we complete the proof of Theorem 4.1 in Section B.4.

B.1 NOTATION:

In this section we gather some notation used throughout the proofs. We use φ(z) =ReLU(z) = max(0, z) with φ (z) = I {z≥0} . For two matrices/vectors x and y of the same size we use x y to denote the entrywise Hadamard product of these two matrices/vectors. We also use x ⊗ y to denote their Kronecker product. For two matrices B ∈ R n×d1 and C ∈ R n×d2 , we use the Khatrio-Rao product as the matrix A = B * C ∈ R n×d1d2 with rows A i given by A i = B i ⊗ C i . For a matrix M ∈ R m×n we use vect(M ) ∈ R mn to denote a vector obtained by aggregating the rows of the matrix M into a vector, i.e. vect(M ) = [M 1 M 2 . . . M m ] T . For a matrix X we use σ min (X) and X denotes the minimum singular value and spectral norm of X. Similarly, for a symetric matrix M we use λ min (M ) to denote its smallest eigenvalue.

B.2 PRELIMINARIES:

In this section we carryout some simple calculations yielding simple formulas for the gradient and Jacobian mappings. We begin by noting we can rewrite the gradient descent iterations in the form

vect (W τ +1 ) = vect (W τ ) -ηvect (∇L (W τ )) .

Here, vect (∇L (W τ )) = J T (W τ )r (W τ ) , where

J (W ) = ∂ ∂vect (W ) f (W ) and

is the Jacobian mapping associated to the network and

r (W ) = φ (V φ (W z)) -y.

is the misfit or residual vector. Note that

∂ ∂vect (W ) v T φ (W z) = v 1 φ w T 1 z x T v 2 φ w T 2 z x T . . . v k φ w T k z x T = (v φ (W x)) T ⊗ x T Thus J (W ) = (V diag (φ (W z))) * 1z T ,

This in turn yields

EQUATION

B.3 LEMMAS FOR CONTROLLING THE SPECTRUM OF THE JACOBIAN AND INITIAL MISFIT:

In this section we state a few lemmas concerning the spectral properties of the Jacobian mapping, its perturbation and initial misfit of the model with the proofs deferred to Appendix C.

Theorem C.1. Assume A 1 , . . . , A d ∈ R n×n are i.i.d. random positive semidefinite matrices whose coordinates have bounded second moments. Define Σ := E[A 1 ] (this is an entry-wise expectation) and

Σ d = 1 d d =1

A .

Let h ∈ (1, +∞) be such that E (u T A 1 u) 2 ≤ hu T Σu for all u ∈ R n . Then for any δ ∈ (0, 1)

we have

P ∀u ∈ R n : u T Σ k u ≥ 1 -7h n + 2 ln(2/δ) d u T Σu ≥ 1 -δ

We shall apply this theorem with A := a a T . To do this we need to calculate the various parameters in the theorem. We begin with Σ and note that for ReLU we have

Σ :=E[A 1 ] =E a 1 a T 1 =E w∼N (0,I k ) (φ ( w, z )) 2 E v∼N (0,In) [vv T ] =E w∼N (0,I k ) φ (w T z) 2 I n =E w∼N (0,I k ) I {w T z≥0} I n = 1 2 I n .

To calculate h we have

E (u T A 1 u) 2 ≤ E a T 1 u 4 ≤ E w∼N (0,I k ) I {w T z≥0} • E v∼N (0,In) (v T u) 4 ≤ 3 2 u 4 ≤ √ 3 √ 2 u 2 = √ 6u T 1 2 I n u = √ 6 • u T Σu.

Thus we can take h = √ 6. Therefore, using Theorem C.1 with δ = 2e -n we can conclude that

λ min 1 d d =1 a a T ≥ 1 4

holds with probability at least 1 -2e -n as long as

d ≥ 3528 • n.

Plugging this into equation 11 we conclude that with probability at least 1 -2e

-n σ min (J (W )) ≥ 1 2 √ d z . C.2 PROOF OF LEMMA B.2

We prove the result for ν = 1, the general result follows from a simple rescaling. Based on equation 8 we have

EQUATION

where S(W ) ⊂ {1, 2, . . . , d} is the set of indices where W z and W 0 z have different signs i.e. S(W ) := { : sgn(e T W z) = sgn(e T W 0 z)} and V S(W ) is a submatrix V obtained by picking the columns corresponding to S(W ).

To continue further note that by Gordon's lemma we have

sup |S|≤s V S ≤ √ n + 2s log(d/s) + t,

with probability at least 1 -e -t 2 /2 . In particular using t = √ n we conclude that

EQUATION

with probability at least 1 -e -n/2 . To continue further we state a lemma controlling the size of |S(W )| based on the size of the radius R. To prove this result we utilize two lemmas from Oymak & Soltanolkotabi (2019) . In these lemmas we use |v| m-to denote the mth smallest entry of v after sorting its entries in terms of absolute value. 

:



Figures and Tables:

Figure 1: Per-pixel reconstruction error (MSE) vs. number of measurements. Vertical bars indicate 95% confidence intervals. BM3D-AMP frequently fails to converge for fewer than 4000 measurements on x-ray images, as denoted by error values far above the vertical axis.

Figure 2: Reconstruction results on x-ray images for m = 2000 measurements (of n = 65536 pixels) and MNIST for m = 75 measurements (of n = 784 pixels). From top to bottom row: original image, reconstructions by our algorithm, then reconstructions by baselines BM3D-AMP, TVAL3, and Lasso. For x-ray images the number of measurements obtained are 3% the number of pixels (i.e. m n = .03), for which BM3D-AMP often fails to converge.

Figure 3: Per-pixel reconstruction error (MSE) vs. number of measurements. Vertical bars indicate 95% confidence intervals.Unfortunately an RGB version of TVAL3 does not currently exist, although related TV algorithms such as FTVd perform similar denoising tasks(Wang et al., 2008).

Reconstruction error (MSE) on MNIST for varying number of measurements. As expected, the trained algorithm of Bora et al. (CSGM) outperforms our method for fewer measurements; however, CSGM saturates after 75 measurements, as its output is constrained to the range of the generator. This saturation is discussed inBora et al., Section  6.1.1.

error (MSE) on x-ray images for varying amounts of noise; number of measurements m fixed at 2000. The term σ 2 corresponds to variance of the noise vector η in y = Ax + η, i.e. each entry of η is drawn independently N (0, σ 2 m ). Other baselines have error far above the vertical axis and are thus not visible in this plot.

Figure 4

Lemma C.2 (sign changes in local neighborhood). Let W 0 ∈ R d×k be a matrix with i.i.d. N (0, 1) entries. Also for a matrix W ∈ R d×k define S(W) := { : sgn(e T W z) = sgn(e T W 0 z)}. Then for any W ∈ R d×k obeying W -W 0 ≤ R |S(W )| ≤ 2 (2dR) at least 1 -e -n/2 -e -

Lemma C.3.(Oymak & Soltanolkotabi, 2019, Lemma C.2) Given an integer m, supposeW -W 0 ≤ √ m |W 0 z| m- z , then |S(W )| ≤ 2m.Lemma C.4.(Oymak & Soltanolkotabi, 2019, Lemma C.3) Let z ∈ R k . Also let W 0 ∈ R d×k be a matrix with i.i.d. N (0, 1) entries. Then, with probability at least 1 -e -m 6 , |W 0 z| m-

Figure 5: Reconstruction results on MNIST for m = 25, 50 measurements respectively (of n = 784 pixels). From top to bottom row: original image, reconstructions by our algorithm, then reconstructions by baselines BM3D-AMP, TVAL3, and Lasso.

Figure 6: Reconstruction results on MNIST for m = 100, 200 measurements respectively (of n = 784 pixels). From top to bottom row: original image, reconstructions by our algorithm, then reconstructions by baselines BM3D-AMP, TVAL3, and Lasso.

Figure 7: Reconstruction results on x-ray images for m = 500, 1000 measurements respectively (of n = 65536 pixels). From top to bottom row: original image, reconstructions by our algorithm, then reconstructions by baselines BM3D-AMP, TVAL3, and Lasso.

Figure 13: A radial sampling pattern of coefficients Ω in the Fourier domain. The measurements are obtained by sampling Fourier coefficients along these radial lines.

Evaluating the benefits of learned regularization (LR) on x-ray images with varying levels of noise and number of measurements. Table values are percent decrease in error, e.g. at σ 2 η = 0 and m = 500, LR reduces MSE by 9.9%. The term σ 2 η corresponds to variance of the noise vector η in Eqn. 1, i.e. each entry of η is drawn independently N (0,

Table content: <table><tr><td>σ 2 η m ). These results indicate that LR</td></tr></table>

Peer review (iclr2020):
This paper proposes a deep learning based compressed sensing method, CS-DIP which employs deep image prior algorithm to recovering signals especially for medical images from noisy measurements without pre-learning over large dataset. For further optimization, a learned regularization is leveraged over small training dataset. Experimental results show that the proposed methods outperformed the others. The theoretical analysis of early stopping is also given. Overall, the writing of this paper is good, the idea is novel and the contribution is sufficient. Here are some of my comments: 1. The authors claim that CS-DIP is *unlearned*, which makes me confusing. Although the latent representation in Eq.(2) is not learned, the weights is still learned from training over measurements. I think maybe the description is not accurate of saying that the network is unlearned. 2. It would be better if the authors can compare CS-DIP with more deep-learning based CS methods. However, the proposed method is only compared with three traditional CS recovery method over several datasets. In Figure 4(a), only one deep-learning based method, CSGM, was compared with over MNIST dataset. More comparison should be conducted over all the datasets to demonstrate its efficiency, especially over medical image datasets.