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):
== After reading authors* response == The authors didn*t really offer a response. Their new PDF just mentions that their proof for the case can be extended to any orthogonal matrix (since if is a iid Gaussian matrix with unit variance, so is for any orthogonal matrix ). Hence, I don*t change my response, and still support *weak accept* == Original review == The paper proposes an extension to the 2017 *Deep Image Prior* (DIP) idea. DIP is a method to use untrained neural nets for compressed sensing, but uses the structure of the neural net to implicitly regularize and constrain the solution (though, as a major result of this paper shows, it doesn*t regularize -- it can learn noise). It*s an interesting combination of compressed sensing and deep learning. This particular paper adds the following to the DIP idea: (1) it adds a total-variation (TV) semi-norm regularization term, (2) it adds a weighted Tikhonov/ridge-regression style regularization term, with learned weight-matrix, and some overall lambda_L hyper-parameter, (3) they provide a theorem, based on a 2019 preprint, that a very simplified version of the unregularized DIP model (1 layer, with one set of weights fixed and assumed drawn Gaussian) can be fully trained by gradient descent, and hence it can overfit to noise, and thus necessarily must be stopped early. They provide some numerical experiments to demonstrate the performance of their approach. == I have some high-level criticism of the approach (things that cannot be fixed): -- The idea of an algorithm that must be stopped early is distasteful to me and perhaps to others as well. It*s vague, and it means the model is wrong. It doesn*t feel like math, it feels like sorcery. You haven*t actually provided a result that says when to stop, or that stopping early at all will help, only that running until convergence is a bad idea. -- Compressed sensing is a bit old these days, and part of the greatness of CS was what it inspired, even if single-pixel cameras have not overturned the camera industry. But pure CS, with a Gaussian measurement matrix, and noiseless measurements, is quite academic and not applicable to the real-world. In particular, your method seems to do best when there are a lot of measurements (for low measurement numbers, the Bora et al approach does well, as does TVAL on the MNIST data -- see Fig 1b). The experiments were also noiseless, and they had Gaussian or the Fourier sampling (which is not random, but rather heavily biased toward DC). I*m not that convinced this is useful for medical imaging. On the plus side, any approach for medical imaging is welcome, even if the connection is a bit academic. The writing was overall good. Furthermore, the proof of the theorem seems like it combines quite a bit of good math. It wasn*t clear if empirical process theory was used in their proof, or if that*s from the 2019 result they cite. == Some major comments (things that *can* be addressed) -- Experimentally, I wasn*t convinced that the learned regularization (LR) was useful. Table 1 appears to show it is helpful, but this is just saying that some kind of l2 regularization is useful. I would be convinced if you showed that vanilla Tikhonov regularization is not as helpful as your learned l2 regularization. -- Furthermore, for the benefit of your LR approach, you did a grid search on the lambda_L parameter to choose the best value. Since you never mentioned trainining/validation/testing data, I assume you did this on the final set of data? In that case, you*ve tuned your hyperparameter to your testing data, and so I don*t find the results that convincing. I think this is a big deal. Or maybe I don*t understand what you*ve done, because the language *.. and selected the best one* was unclear what *best* meant, e.g., if you look at a residual, or if you use ground-truth data. -- For both your LR and vanilla Tikhonov, you can try selecting lambda_L based on the residual or other non-oracle method. Standard methods include generalized cross-validation, Morozov’s discrepancy principle, and the unbiased predictive risk estimator method. -- For the results comparing with Bora et al., how much of your performance advantage was due to the TV term? Please re-do Fig 4b with variants of your method with lambda_T=0, lambda_L=0, and both lambdas=0 added to the plot. -- The handling of the TV term was not discussed. It seems you fed this into PyTorch and hoped for the best. Actually minimizing the usual isotropic 2D TV is tricky, and there are many optimization papers written about it. It*s not differentiable, and if you actually converge, then you usually land at points of non-differentiability, so it cannot just be ignored. Sometimes it is smoothed, or a separable but non-isotropic version is used. I think your response to this point is, *but we don*t actually run the algorithm to convergence* (since it will overfit), which I still find bizarre (imagine saying that in a paper submitted to an optimization journal!). == Summary of my opinion == I liked the idea, but I think I liked the core idea of DIP more than the extensions that this paper provides. The experiments show enough promising results to make this interesting, but they have not explored the details of the extensions. The theorem is very complex but quite limited in scope, and maybe would be better as its own (math) result rather than as an add-on to this paper. The application of CS is a bit limited, and the biggest benefits seem to be in the high-measurement and low or zero-noise case. Overall, I lean toward weak reject, as I think the regularization parts of the paper would benefit from a major revision of the paper. == Minor comments == - In section 3.1, the description of G is confusing. Writing G(z;w): R^k --> R^n is unclear. You mean G(.;w): R^k --> R^n. And calling z a *seed* confused me, until I realized you meant this is the *input* to the layer. - statements like *we use early stopping, a phenomena that will be analyzed theoretically in Section 4* are misleading, since you don*t analyze early stopping at all. You show that your model doesn*t regularize, so that it will fit noise (in a special simplified case), so that running to convergence is a bad idea. This doesn*t mean you know when to stop. - Section 3.2, for p(y|w) and p(w), why not just say these are Gaussian and give the parameters? In the sentence above eq (4), missing a *log* before *posterior*. The justification for the Gaussian weights (last paragraph on page 4) felt rather weak to me. Like everyone for the past 60 years, you are adding a l2 regularization term because it*s simple and it works fairly well and there are not that many alternatives. - It took me a while to figure out that Sigma was L x L. The description here was confusing (in section 3.2.1), and it was unclear of Sigma_U was a sub-matrix (and if it was diagonal or a scaled identity) or a scalar. - Thm 4.1: -- *V is fixed at random* sounds like it would make a proper probabilist squirm! The statement of the theorem is that *... holds ... with probability at least ...*. You mean that this probability is due to the randomness in V? -- You have a 1-layer network with 2 sets of weights, but in fact fix one set of weights; and it*s not compressed sensing but denoising since A=I. I understand that results in this area are hard, but make sure not to oversell your result earlier in the paper, since this limited setting doesn*t actually apply to the main problem you*re interested in. -- *While our proof is for the case of A=I, a similar result can be shown for other measurement matrices.* If that*s true, then please show it! - Table 1: for with and without LR, did you have a TV term for both? - Section 5.2.1, *thus we can infer that assuming a learned Gaussian distribution over weights is useful.* No, I don*t think this is an appropriate assumption. Maybe you can infer that shrinking the estimator toward a fixed vector (e.g., 0), and training to find how much to shrink, makes your solution more robust (a bias-variance tradeoff). In particular, in Table 1, I*d like to see the improvements for one of these entries as a box plot, to see if the improvement is due to controlling bad outliers, or if actually the bulk is really improved. - In bibliography, capitalize acronyms like *admm*, *amp*, *mri*, etc.