Title: DENOISING AND REGULARIZATION VIA EXPLOITING THE STRUCTURAL BIAS OF CONVOLUTIONAL GENER-ATORS

Abstract: Convolutional Neural Networks (CNNs) have emerged as highly successful tools for image generation, recovery, and restoration. A major contributing factor to this success is that convolutional networks impose strong prior assumptions about natural images. A surprising experiment that highlights this architectural bias towards natural images is that one can remove noise and corruptions from a natural image without using any training data, by simply fitting (via gradient descent) a randomly initialized, over-parameterized convolutional generator to the corrupted image. While this over-parameterized network can fit the corrupted image perfectly, surprisingly after a few iterations of gradient descent it generates an almost uncorrupted image. This intriguing phenomena enables state-of-the-art CNNbased denoising and regularization of other inverse problems. In this paper we attribute this effect to a particular architectural choice of convolutional networks, namely convolutions with fixed interpolating filters. We then formally characterize the dynamics of fitting a two layer convolutional generator to a noisy signal and prove that early-stopped gradient descent denoises/regularizes. Our proof relies on showing that convolutional generators fit the structured part of an image significantly faster than the corrupted portion.

1. INTRODUCTION:

Convolutional neural networks are extremely popular for image generation. The majority of image generating networks is convolutional, ranging from deep convolutional Generative Adversarial Networks (DC-GANs) (Radford et al., 2015) to the U-Net (Ronneberger et al., 2015) . It is well known that convolutional neural networks incorporate implicit assumptions about the signals they generate, such as pixels that are close being related. This makes them particularly well suited for modeling distributions of images. It is less well known, however, that those prior assumptions build into the architecture are so strong that convolutional neural networks are useful even without ever being exposed to training data.

The latter was first shown in the Deep Image Prior (DIP) paper (Ulyanov et al., 2018) . Ulyanov et al. (2018) observed that when 'training' a standard convolutional auto-encoder such as the popular U-net (Ronneberger et al., 2015) on a single noisy image and regularizing by early stopping, the network denoises the image with excellent performance. This is based on the empirical observation that un-trained convolutional auto-decoders fit a single natural image faster when optimized with gradient descent than pure noise. Many elements of an auto-encoder, however are irrelevant for this effect: A more recent paper (Heckel & Hand, 2019) proposed a much simpler image generating network, termed the deep decoder. This network can be seen as the relevant part of a convolutional generator architecture to function as an image prior, and can be obtained from a standard convolutional autoencoder by removing the encoder, the skip connections, and perhaps most notably, the trainable convolutional filters of spatial extent larger than one. Thus, the deep decoder does not use learned or trainable convolutional filters like conventional convolutional networks do, and instead only uses convolutions with fixed convolutional kernels to generate an image. In this paper, we study such a simple, untrained convolutional network theoretically. We consider the over-parameterized regime where the network has sufficiently many parameters to represent an arbitrary image (including noise) perfectly and show that: Fitting convolutional generators via early stopped gradient descent provably denoises "natural" images.

To prove this statement, we characterize how the network architecture governs the dynamics of fitting over-parameterzed networks to a single (noisy) image. In particular we prove:

Convolutional generators optimized with gradient descent fit natural images faster than noise.

We depict this phenomena in Figure 1 where we fit a randomly initialized over-parameterized convolutional generator to an image via running gradient descent on the objective L(C) = G(C) -y 2 2 . Here, G(C) is the convolutional generator with weight parameters C, and y is either a noisy image, a clean image, or noise. This experiment demonstrates that an over-parameterized convolutional network fits a natural image (Figure 1b ) much faster than noise (Figure 1c ). Thus, when fitting the noisy image (Figure 1a ), early stopping the optimization enables image denoising. This effect is so strong that it gives state-of-the-art denoising performance (see Figure 2 and Appendix A). Beyond denoising, this effect also enables significant improvements in regularizing a variety of inverse problems such as compressive sensing (Veen et al., 2018; Heckel, 2019) .

1.1. CONTRIBUTIONS AND OVERVIEW OF RESULTS:

In this paper we take a step towards understanding why fitting a convolutional generator via early stopped gradient descent leads to such surprising denoising and regularization capabilities. • We first show experimentally that this denoising phenomena is primarily attributed to convolutions with fixed interpolation kernels, typically implicitly implemented by bi-linear upsampling operations in the convolutional network. • We then show that fitting an over-parameterized convolutional networks (with fixed convolutional filters) via early stopped gradient descent to a signal provably denoises it. Specifically, let x ∈ R n be a smooth signal that can be represented exactly as a linear combination of the p orthogonal trigonometric functions of lowest frequency (defined in equation ( 5), see Figure 3 for a depiction). A smooth signal is a good model for a natural image since natural images are well approximated by low-frequency components. Specifically, Figure 4 in (Simoncelli & Olshausen, 2001) shows that the power spectrum of a natural image (i.e., the energy distribution by frequency) decays rapidly from low frequencies to high frequencies. In contrast, Gaussian noise has a flat power spectrum.

Our goal is to obtain an estimate of this signal from a noisy observation y = x + z where z ∼ N (0, ς 2 n I) is Gaussian noise. Let ŷ = G(C τ ) be the estimate obtained from early stopping the fitting of a two-layer convolutional generator to y. We prove this estimate achieves the denoising rate ŷ -

x 2 2 ≤ c p n ς 2 ,

with c a constant. This rate is optimal up to a constant factor. • Our denoising result follows from a detailed analysis of gradient descent applied to fitting a convolutional generator G(C) to a noisy signal y = x+z with x representing the signal and z the noise. We show that there are weights σ 1 , . . . , σ n associated with the trigonometric basis function that only depend on the convolutional filter used, and which typically decay quickly from the low-frequency to the high-frequency trigonometric basis functions. These weights in turn determine the speed at which the associated components of the signal are fitted. Specifically, we show that the dynamics of gradient descent are approximately given by

G(C τ ) -x ≈ n i=1 w i w i , x (1 -ησ 2 i ) τ error in fitting signal + n i=1 w i w i , z ((1 -ησ 2 i ) τ -1) fit of noise .

The convolutional filters commonly used are such that the weights σ 1 , . . . , σ n decays quickly, implying that low-frequency components in the trigonometric expansion are fitted significantly faster than high frequency ones. So if the signal mostly consists of lowfrequency components, we can choose an early stopping time such that the error in fitting the signal is very low, and thus the signal part is well described, whereas at the same time only a small part of the noise, specifically the part aligned with the low-frequency components has been fitted.

2. CONVOLUTIONAL GENERATORS:

A convolutional generator maps an input tensor B 0 to an image only using upsampling and convolutional operations, followed by channel normalization (a special case of batch normalization) and applications of non-linearities, see Figure 4 . All previously mentioned convolutional generator networks (Radford et al., 2015; Ronneberger et al., 2015) including the networks considered in the DIP paper (Ulyanov et al., 2018) primarily consist of those operations.

For motivating the architecture of the convolutional generators studied in this paper, we first demonstrate in Section 2.1 that convolutions with fixed interpolation filters are critical to the denoising performance with early stopping. Specifically, we empirically show that convolutions with fixed convolutional kernels are critical for convolutional generators to fit natural images faster than noise. In Section 2.2 we formally introduce the class of convolutional generators studied in this paper, and finally, in Section 2.3, we introduce a minimal convolutional architecture which is the focus of our theoretical results.

2.1. THE IMPORTANCE OF FIXED CONVOLUTIONAL FILTERS:

Convolutions with fixed convolutional kernels are critical for denoising with early stopping, because they are critical for the phenomena that natural images are fitted significantly faster than noise. To see this, consider the experiment in Figure 4 in which we fit an image and noise by minimizing the least-squares loss via gradient descent with i) a convolutional generator with only fixed convolutional filters (see Section 2.2 below for a precise description) and ii) a conventional convolutional generator with trainable convolutional filters (essentially the architecture from the popular DC-GAN generators, see Appendix B for details and additional numerical evidence). As illustrated in Figure 4 , the convolutional network with fixed filters fits the natural image much faster than noise, whereas the network with learned convolutional filters, only fits it slightly faster, and this effect vanishes as the network becomes highly overparameterized. Thus, fixed convolutional filters enable un-trained convolutional networks to function as highly effective image priors. We note that the upsampling operation present in most architectures implicitly incorporates a convolution with a fixed convolutional (interpolation) filter.

2.2. ARCHITECTURE OF CONVOLUTIONAL GENERATOR WITH FIXED CONVOLUTIONS:

In this section, we describe the architecture of the deep decoder (Heckel & Hand, 2019) , a convolutional network with fixed convolutional operators only. In this architecture, the channels in the (i + 1)-th layer are given by

B i+1 = cn(ReLU(U i B i C i )), i = 0, . . . , d -1,

and finally, the output of the d-layer network is formed as

x = B d C d+1 .

Here, the coefficient matrices C i ∈ R k×k and C d+1 ∈ R k×kout contain the weights of the network. The number of channels, k, determines the number of weight parameters of the network, given by dk 2 + k out k. Each column of the tensor B i C i ∈ R ni×k is formed by taking linear combinations of the channels of the tensor B i in a way that is consistent across all pixels, and the ReLU activation function is given by ReLU(t) = max(0, t). Then, cn(•) performs the channel normalization operation which normalizes each channel individually and can be viewed as a special case of the popular batch normalization operation (Ioffe & Szegedy, 2015) .

The operator U i ∈ R ni+1×ni is a tensor implementing an upsampling and most importantly a convolution operation with a fixed kernel. This fixed kernel was chosen in all experiments above as a triangular kernel so that U performs bi-linear 2x upsampling (this is the standard implementation in the popular software packages pytorch and tensorflow). As mentioned earlier this convolution with a fixed kernel is critical for fitting natural images faster than complex ones.

2.3. TWO LAYER CONVOLUTIONAL GENERATOR STUDIED THEORETICALLY IN THIS PAPER:

The simplest model to study the denoising capability of convolutional generators and the phenomena that a natural image is fitted faster than a complex one theoretically is a network with only one hidden layers and one output channel i.e., G(C) ∈ R n . Then, the generator becomes

G(C) = ReLU(UB 0 C 0 )c 1 ,

where U ∈ R n×n is a circulant matrix that implements a convolution with a filter u. In this paper we consider the over-parameterized regime where k ≥ 2n. Note that scaling the i-th column of C 0 with a non-negativ factor is equivalent to scaling the i-th entry of the output weights c 1 with the same factor. We therefore fix the output weights to

c 1 = v, where v = [1, . . . , 1, -1, . . . , -1]/ √ k.

Next, note that both the generators with the fixed output weights and non-fixed output weights have the same range in the regime k ≥ 2n. Next, because the matrix B 0 is Gaussian, with probability one, it has full rank and thus spans R n . It follows that optimizing over the parameter C 0 in C = B 0 C 0 is equivalent to optimizing over the matrix C ∈ R n×k . We therefore consider the generator

EQUATION

where v = [1, . . . , 1, -1, . . . , -1]/ √ k is fixed and C ∈ R n×k is the new coefficient matrix we optimize over. Figure 11 in the appendix shows that even this simple two-layer convolutional network fits a simple image faster than noise. This is the simplest model in which the phenomena that a convolutional networks fits structure faster than noise can reliably be observed. As a consequence, the dynamics of training the model (1) are the focus of the remainder of this paper.

3. WARMUP: DYNAMICS OF GRADIENT DESCENT ON LEAST SQUARES:

As a prelude for studying the dynamics of fitting convolutional generators via a non-linear least square problem, we study the dynamics of gradient descent applied to a linear least squares problem. We demonstrate how early stopping can lead to denoising capabilities even with a simple linear model. We consider a least-squares problem of the form

L(c) = 1 2 y -Jc 2 2

, and study gradient descent with a constant step size η starting at c 0 = 0. The updates are given by

c τ +1 = c τ -η∇L(c τ ), ∇L(c) = J T (Jc -y).

The following simple proposition characterizes the trajectory of gradient descent. Proposition 1. Let J ∈ R n×m be a matrix with left singular vectors w 1 , . . . , w n ∈ R n and corresponding singular values σ 1 ≥ σ 2 ≥ . . . ≥ σ n . Then the residual after τ steps, r τ = y -Jc τ , of gradient descent starting at c 0 = 0 is

r τ = n i=1 w i w i , y (1 -ησ 2 i ) τ .

Suppose that the signal y lies in the column span of J, and that the stepsize is chosen sufficiently small (i.e., η ≤ 1/ J 2 ). Then, by Proposition 1, gradient descent converges to a zero-loss solution and thus fits the signal perfectly. More importantly, gradient descent fits the components of y corresponding to large singular values faster than it fits the components corresponding to small singular values.

To explicitly show how this observation enables regularization via early stopped gradient descent, suppose our goal is to find a good estimate of a signal x from a noisy observation where the signal x lies in a signal subspace that is spanned by the p leading left-singular vectors of J. Then, by Proposition 1, the signal estimate after τ iterations, Jc τ , obeys

y = x + z,

Jc τ -x 2 ≤ (1 -ησ 2 p ) τ x 2 + E(z), E(z) := n i=1 ((1 -ησ 2 i ) τ -1) 2 w i , z 2 .

(2) Thus, after a few iterations most of the signal has been fitted (i.e., (1 -ησ p ) τ is small). Furthermore, if we assume that the ratio σ p /σ p+1 is sufficiently large so that the spread between the two singular values separating the signal subspace from the rest is sufficiently large, most of the noise outside the signal subspace has not been fitted (i.e., ((1

-ησ 2 i ) τ -1) 2 ≈ 0 for i = p + 1, . . . , n).

In particular, suppose the noise vector has a Gaussian distribution given by z ∼ N (0, ς 2 n I). Then E(z) ≈ ς p n so that after order τ = log( )/ log(1 -ησ 2 p ) iterations, with high probability,

Jc τ -x 2 ≤ x 2 + cς p n .

This demonstrates, that provided the signal lies in a subspace spanned by the leading singular vectors, early stoped gradient descent reaches the optimal denoising rate of ς p/n after a few iterations. See Figure 5 for a numerical example demonstrating this phenomena.

4. DYNAMICS OF GRADIENT DESCENT ON CONVOLUTIONAL GENERATORS:

We are now ready to study the implicit bias of gradient descent towards natural/structured images theoretically. Consider the two-layer network (introduced in Section 2.3) of the form

G(C) = ReLU(UC)v, where v = [1, . . . , 1, -1, . . . , -1]/ √ k,

and with weight parameter C ∈ R n×k , and recall that U is a circulant matrix implementing a convolution with a kernel u ∈ R n . We fit the generator to a signal y ∈ R n by minimizing the non-linear least squares objective

L(C) = 1 2 y -G(C) 2 2

(3) with (early-stopped gradient) descent with a constant stepsize η starting at a random initialization C 0 of the weights. The iterates are given by

EQUATION

In our warmup section on linear least squares we saw that the singular vectors and values of the matrix J determine the speed at which different components of the noisy signal y are fitted by gradient descent. The main insight that enables us to extend this intuition to the nonlinear case is that the role of the matrix J can be replaced with the Jacobian of the generator, defined as Contrary to the linear least squares problem, however, in the nonlinear case, the Jacobian is not constant and changes across iterations. Nevertheless, we show that the eigen-values and vectors of the Jacobian at the random initialization govern the dynamics of fitting the network throughout the iterative updates.

J (C) := ∂ ∂C G(C).

For the two-layer convolutional generator that we consider, the left eigenvectors of the Jacobian mapping can be well approximated by the trigonometric basis functions, defined below, throughout the updates. Interestingly, the form of these eigenvectors only depends on the network architecture and not the convolutional kernel used.

Definition 1. The trigonometric basis functions w 1 , . . . , w n are defined as

EQUATION

)

Figure 3 depicts some of these eigenvectors.

In addition to the left eigenvectors we can also approximate the spectrum of the Jacobian throughout the updates by an associated filter/kernel that only depends on the original filter/kernel used in the network architecture.

Definition 2 (Dual kernel). Associated with a kernel u ∈ R n we define the dual kernel σ ∈ R n as

σ = u 2 Fg u u u 2 2 with g(t) = 1 2 1 - cos -1 (t) π t.

Here, for two vectors u, v ∈ R n , u v denotes their circular convolution, the scalar non-linearity g is applied entrywise, and F is the discrete Fourier transform matrix.

In Figure 6 , we depict two commonly used interpolation kernels u, namely a triangular and a Gaussian kernel (recall that the standard upsampling operator is a convolution with a triangle), along with the induced dual kernel σ. The figure shows that the dual kernel σ induced by these kernels has a few large values associated with the low frequency trigonometric functions, and the other values associated with high frequencies are very small.

With these definitions in place we are ready to state our main denoising result. A denoising result requires a signal model-we assume a low-frequency signal x that can be represented as a linear combination of the first p-trigonometric basis functions. Note that this is a good model for a natural image since natural images are well approximated by low-frequency components. Specifically, Figure 4 in (Simoncelli & Olshausen, 2001) shows that the power spectrum of a natural image (i.e., the energy distribution by frequency) decays rapidly from low frequencies to high frequencies. In contrast, Gaussian noise has a flat power spectrum.

Theorem 1 (Denoising with early stopping). Let x ∈ R n be a signal in the span of the first p trigonometric basis functions, and consider a noisy observation

y = x + z,

where z is Gaussian noise with distribution N (0, ς 2 n I), for some variance ς 2 ≥ 0. To denoise this signal, we fit a two layer generator network

EQUATION

channels, for some > 0, and with convolutional kernel u of the convolutional operator U and associated dual kernel σ, to the noisy signal y. Here, C u a constant only depending on the convolutional kernel u. Then, with probability at least 1-e -k 2 -1 n 2 , the reconstruction error obtained after

EQUATION

)

Note that for this choice of stopping time, provided that dual kernel decays sharply around the p-th singular value, the first term in the bound (7) (i.e., (1 -ησ 2 p ) τ ≈ 0) essentially vanishes and the error bound becomes O(ς p n ). The dual kernel decays sharply around the leading eigenvalues provided the kernel is for example a sufficiently wide triangular or Gaussian kernel (see Figure 6 ). This result demonstrates that when the noiseless signal x is sufficiently structured (e.g. contains only the p lowest frequency components in the trigonometric basis) and the convolutional generator has sufficiently many channels, then early stopped gradient descent achieves a near optimal denoising performance proportional to ς p n . This theorem is obtained from a more general result stated in the appendix which characterizes the evolution of the reconstruction error obtained by the convolutional generator. Theorem 2 (Reconstruction dynamics of convolutional generators). Consider the setting and assumptions of Theorem 1 but now with a fixed noise vector z, and without an explicit stopping time. Then, for all iterates τ obeying τ ≤ 100 ησ 2 p and provided that k ≥ C u n/ 8 , for some ∈ (0, σp σ1 ], with probability at least 1 -e -k 2 -1 n 2 , the reconstruction error obeys

G(C τ ) -x 2 ≤ (1 -ησ 2 p ) τ x 2 + n i=1 ((1 -ησ 2 i ) τ -1) 2 w i , z 2 + y 2 .

This theorem characterizes the reconstruction dynamics of convolutional generators throughout the updates. In particular, it explains why convolutional generators fit a natural signal significantly faster than noise, and thus early stopping enables denoising and regularization. To see this, note that as mentioned previously each of the basis functions w i has a (positive) weight σ i > 0 associated with it that only depends on the convolutional kernel used in the architecture (through the definition of the dual kernel). These weights determine how fast the different components of the noisy signal are fitted by gradient descent. As we demonstrated earlier in Figure 6 , for typical convolutional filters those weights decay very quickly from low to high frequency basis functions. As a result, when the signal x is sufficiently structured (i.e. lies in the range of the p trigonometric functions with lowest frequencies), after a few iterations most of the signal is fitted (i.e., (1 -ησ 2 p ) τ is small), while most of the noise has not been fitted (i.e., ((1 -ησ 2 i ) τ -1) 2 ≈ 0 for i = p + 1, . . . , n). Thus, early stopping achieves denoising.

The proof, provided in the appendix, is based on associating the following linear least-squares problem with the non-linear least squares problem (3):

L lin (c) = 1 2 G(C 0 ) + J(c -vect(C 0 )) -y 2 2 .

Here, vect(C 0 ) ∈ R kn is a vectorized version of the matrix C 0 , and J ∈ R n×nk approximates the Jacobian of G(C) around the initialization C 0 . The linear least-squares problem is obtained by linearizing the non-linear problem around the initialization. In the proof we show that if i) the network is sufficiently over-parameterized and if ii) the network is randomly initialized, then the linearized problem is a good approximation of the non-linear problem. We then show that the singular values of the Jacobian J are the trigonometric basis functions, and the singular values are the values of the dual kernel. The proof is then concluded by characterizing the trajectory of gradient descent applied to the linear least-squares problem above.

4.1. MULTILAYER NETWORKS AND MODERATE OVER-PARAMETERIZATION:

Our theoretical results rely on the finding that for over-parameterized single hidden-layer networks, the leading singular vectors of the Jacobian are the trigonometric functions throughout all iterations, and that the associated weights (i.e., the singular values) are concentrated towards the low frequency components. In this regime, the dynamics are well approximated by a linear model. This general strategy can be extended to multi-layer networks, however if the network is not highly overparameterized, an associated linear model might not be a good approximation.

In order to understand whether our finding of low-frequency components being being fitted faster than high-frequency ones carries over to multi-layer networks and to the moderately overparameterized regime, in this section, we study muli-layer networks in the moderately overparameterized regime numerically. We show that for a multilayer network, the spectrum of the Jacobian is concentrated towards singular vectors/functions that are similar to the low-frequency components.

We also show that throughout training those functions do vary, albeit the low frequency components do not change significantly and the spectrum remains concentrated towards the low frequency components. This shows that the implications of our theory continue to apply to muli-layer networks.

In more detail, we take a standard one dimensional deep decoder with d = 4 layers with output in R 512 and with k = 64 channels in each layer. Recall that the standard one dimensional decoder obtains layer i + 1 from layer i by linearly combining the channels of layer i with learnable coefficients followed by linear upsampling (which involves convolution with the triangular kernel [1/2, 1, 1/2]).

The number of parameters is d × k 2 = 32 • 512, so the network is over-parameterized by a factor of 32. In Figure 7 , we display the singular values as well as the leading singular vectors/function of the Jacobian at initialization (t = 1) and after t = 50 and t = 3500 iterations of gradient descent. As can be seen the leading singular vectors (s = 1-5) are close to the trigonometric basis functions and do not change dramatically throughout training. The singular vectors corresponding to increasingly smaller singular values (s = 20, 50, 100, 150) contain increasingly higher frequency components but are far from the high-frequency trigonometric basis functions.

5. RELATED LITERATURE:

As mentioned before, the DIP paper (Ulyanov et al., 2018) was the first to show that overparameterized convolutional networks enable solving denoising, inpainting, and super-resolution problems well even without any training data. Subsequently, the paper (Heckel & Hand, 2019) proposed a much simpler image generating network, termed the deep decoder. The papers (Veen et al., 2018; Heckel, 2019; Jagatap & Hegde, 2019; Bostan et al., 2020) have shown that the DIP and the deep decoder also enable solving or regularizing compressive sensing problems and other inverse problems.

Since the convolutional generators considered here are image-generating deep networks, our work is also related to methods that rely on trained deep image models. Deep learning based methods are either trained end-to-end for tasks ranging from compression (Toderici et al., 2016; Agustsson et al., 2017; Theis et al., 2017; Burger et al., 2012; Zhang et al., 2017) to denoising (Burger et al., 2012; Zhang et al., 2017) , or are based on learning a generative image model (by training an autoencoder or GAN (Hinton & Salakhutdinov, 2006; Goodfellow et al., 2014) ) and then using the resulting model to solve inverse problems such as compressed sensing (Bora et al., 2017; Hand & Voroninski, 2018) , denoising (Heckel et al., 2018) , or phase retrieval (Hand et al., 2018) , by minimizing an associated loss. In contrast to the method studied here, where the optimization is over the weights A large body of work focuses on understanding the optimization landscape of the simple nonlinearities or neural networks (Candes et al., 2015; Soltanolkotabi, 2017; Alon Brutzkus & Globerson, 2017; Zhong et al., 2017; Oymak, 2018; Fu et al., 2018; Tu et al., 2016) when the labels are created according to a planted model. Our proofs rely on showing that the dynamics of gradient descent on an over-parameterized network can be related to that of a linear network or a kernel problem. This proof technique has been utilized in a variety of recent publication (Soltanolkotabi et al., 2018; Venturi et al., 2019; Du et al., 2018; Oymak & Soltanolkotabi, 2019a; b; Arora et al., 2019; Oymak et al., 2019) . Two recent publication have used this proof technique to show that functions are learned at different rates: A DENOISING PERFORMANCE OF UNTRAINED CONVOLUTIONAL

GENERATORS:

In this section, we provide further details on the denoising performance of deep neural networks. We compare the denoising performance of four methods: i) the BM3D algorithm (Dabov et al., 2007) as a standard baseline for denoising, ii) the deep image prior (Ulyanov et al., 2018) , which is a U-net like encoder-decoder convolutional architecture applied with exactly the parameters as proposed for denoising in the paper (Ulyanov et al., 2018) , and early stopped after 1900 iterations, again with the same stopping time as proposed in the original paper, iii) an under-parameterized deep decoder with five layers and k = 128 channels in each layers, trained for 3000 iterations (which is close to convergence), iv) an over-parameterized deep decoder with five layers and k = 512 channels in each layer, early stopped at 1500 iterations. We compared the performance of those four methods on denoising 100 randomly chosen images from the ImageNet validation set. The code to reproduce the results contains a list of the images we considered. Each image has 512 × 512 pixels and three color channels. We added the same random Gaussian noise to each color channel, because we are interested in evaluating the performance by imposing structural assumptions on the image. 

B A NUMERICAL STUDY OF THE IMPLICIT BIAS OF CONVOLUTIONAL:



NETWORKS:

In this section, we empirically demonstrate that convolutions with fixed convolutional kernels are critical for convolutional generators to fit natural images faster than noise. Towards this goal, we study numerically the following four closely related architectural choices, which differ in the upsampling/no-upsampling and convolutional operations which generate the activations in the (i + 1)-st layer, B i+1 , from the activations in the i-th layer, B i : i) Bilinear upsampling and linear combinations. Layer i + 1 is obtained by linearly combining the channels of layer i with learnable coefficients (i.e., performing one-times-one convolutions), followed by bi-linear upsampling. This is the deep decoder architecture from (Heckel et al., 2018) . ii) Fixed interpolation kernels and linear combinations. Layer i + 1 is obtained by linearly combining the channels of layer i with learnable coefficients followed by convolving each channel with the same 4x4 interpolation kernel that is used in the linear upsampling operator. iii) Parameterized convolutions: Layer i + 1 is obtained from layer i though a convolutional layer. iv) Deconvolutional network: Layer i + 1 is obtained from layer i though a deconvolution layer. This is essentially the DC-GAN (Radford et al., 2015) generator architecture.  and 16 ). Gradient descent on convolutional generators involving fixed convolutional matrixes fit an image significantly faster than noise.

To emphasize that architectures i)-iv) are structurally very similar operations, we recall that each operation consists only of upsampling and convolutional operations. Let T(c) : R n → R n be the convolutional operator with kernel c, let u the linear upsampling kernel (equal to u = [0.5, 1, 0.5] in the one-dimensional case), and let U : R n → R 2n be an upsampling operator, that in the one dimensional case transforms [x 1 , x 2 , . . . , x n ] to [x 1 , 0, x 2 , 0, . . . , x n , 0]. In each of the architectures i)-iv), the -th channel of layer i + 1 is obtained from the channels in the i-th layer as: b i+1, = ReLU k j=1 M(c ij )b i , where the linear operator M is defined as follows for the four architectures i)

M(c) = cT(u)U, ii) M(c) = cT(u), iii) M(c) = T(c), iv) M(c) = T(c)U.

The coefficients associated with the i-th layer are given by C i = {c ij }, and all coefficients of the networks are C = {c ij }. Note that here, the coefficients or parameters of the networks are the weights and not the input to the network.

B.1 DEMONSTRATING IMPLICIT BIAS OF CONVOLUTIONAL GENERATORS:

We next show that convolutional generators with fixed convolutional operations fit natural or simple images significantly faster than complex images or noise. Throughout this section, for each image or signal x * we fit weights by minimizing the loss

L(C) = G(C) -x * 2 2

with respect to the network parameters C using plain gradient descent with a fixed stepsize.

In order to exemplify the effect, we fit the phantom MRI image as well as noise for each of the architectures above for a 5-layer network. We choose the number of channels, k, such that the over-parameterization factor (i.e., the ratio of number of parameters of the network over the output dimensionality) is 1, 4, and 16, respectively. The results in Figure 8 show that for architectures i) and ii) involving fixed convolutional operations, gradient descent requires more than one order of magnitude fewer iterations to obtain a good fit of the phantom MRI image relative to noise. For architectures iii) and iv), with trainable convolutional filters, we see a smaller effect, but the effect essentially vanishes when the network is highly over-parameterized. This effect continues to exist for natural images in general, as demonstrated by Figure 10 which depicts the average and standard deviation of the loss curves of 100 randomly chosen images from the imagenet dataset.

We also note that the effect continues to exist in the sense that highly structured images with a large number of discontinuities are difficult to fit. An example is the checkerboard image in which each pixel alternates between 1 and 0; this image leads to the same loss curves as noise.

In our next experiment, we highlight that the distance between final and initial network weights is a key feature that determines the difference of fitting a natural image and noise. Towards this goal, we again fit the phantom MRI image and noise for the architecture i) and an over-parameterization factor of 4 and record, for each layer i the relative distance

C (t) i -C (0) i / C (0) i

, where C (0) i are the weights at initialization (we initialize randomly), and C (t) i are the weights at the optimizer step t. The results, plotted in Figure 9 , show that to fit the noise, the weights have to change significantly, while for fitting a natural image they only change slightly.

C THE SPECTRUM OF THE JACOBIAN OF THE DEEP DECODER AND DEEP IMAGE PRIOR:

Our theoretical results predict that for over-parameterized networks, the parts of the signal that is aligned with the leading singular vectors of the Jacobian at initialization is fitted fastest. In this section we briefly show that natural images are much more aligned with the leading singular vectors than with Gaussian noise, which is equally aligned with each of the singular vectors. For both deep decoder and deep image prior, this quantity is significantly smaller for noise than for a natural image. Thus, a natural image is better aligned with the leading singular vectors of the Jacobian than noise. This demonstrates that the Jacobian of the networks is approximately low-rank, with natural images lying in the space spanned by the leading singularvectors.

Towards this goal, we compute the norm of the product of the Jacobian at a random initalization, C 0 , with a signal y as this measures the extent to which the signal is aligned with the leading singular vectors, due to

J T (C 0 )y 2 = VΣW T y 2 = i σ 2 i w i , y 2 ,

where J (C 0 ) = WΣV T is the singular value decomposition of the Jacobian.

Figure 12 depicts the distribution of the norm of the product of the Jacobian at initialization, J (C 0 ), with an image y * or noise z of equal norm ( y * = z ). Since for both the deep decoder and the deep image prior, the norm of product of the Jacobian and the noise (i.e., J T (C 0 )z ) is significantly smaller than that with a natural image (i.e., J T (C 0 )y * ), it follows that a structured image is much better aligned with the leading singular vectors of the Jacobian than the Gaussian noise, which is approximately equally aligned with any of the singular vectors. Thus, the Jacobian at random initialization has an approximate low-rank structure, with natural images lying in the space spanned by the leading singularvectors.

D PROOFS AND FORMAL STATEMENT OF RESULTS:

The results stated in the main text are obtained from a slightly more general result which applies beyond convolutional networks. Consider neural network generators of the form

EQUATION

with C ∈ R n×k , and U ∈ R n×n an arbitrary fixed matrix, and v ∈ R k , with half of the entries equal to +1/ √ k and the other half equal to -1/ √ k.

The (transposed) Jacobian of ReLU(Uc) is U T diag(ReLU (Uc)). Thus the Jacobian of G(C) is given by

EQUATION

)

where ReLU is the derivative of the activation function. Next we define a notion of expected Jacobian. Towards this goal, we first define the matrix

Σ(U) := E J (C)J T (C) ∈ R n×n associated with the generator G(C) = ReLU(UC)v.

Here, expectation is over C with iid N (0, ω 2 ) entries. Note that since the derivative of the ReLU non-linearity is invariant to the scaling of C, so is the Jacobian J (C), but to be consistent throughout we take expectation over C with iid N (0, ω 2 ) entries. Consider the eigenvalue decomposition of Σ(U) given by

Σ(U) = n i=1 σ 2 i w i w T i .

Our results depend on the largest and smallest eigenvalue of Σ(U), defined throughout as

α = σ 2 n , β = σ 2 1 = U .

With these definitions in place we are ready to state our result about neural generators. Theorem 3. Consider a noisy signal y ∈ R n given by y = x + z, where x ∈ R n is assumed to lie in the signal subspace spanned by the p leading singular vectors w 1 , . . . , w n of Σ(U), and z ∈ R n is an arbitrary noise vector. Suppose that the number of channels obeys

EQUATION

)

where c is a fixed numerical constant and ξ an error tolerance parameter obeying 0 < ξ ≤ 1/ 32 log 2n δ . Moreover, T is the maximal number of gradient descent steps, and we assume it obeys 1 ≤ T ≤ 2 5 β 2 ηξ 2 α 2 . We fit the neural generator G(C) to the noisy signal y ∈ R n by minimizing the loss

L(C) = 1 2 G(C) -y 2 2 (11)

via running gradient descent starting from C 0 with i.i.d. N (0, ω 2 ) entries, ω = y 2 √ nβ ξ α 2 β 2 , and step size obeying η ≤ 1/β 2 . Then, with probability at least 1 -e -k 2 -δ, for all iterations τ ≤ T ,

EQUATION

)

D.1 PROOF OF THEOREM 2

Theorem 2 stated in the main text follow directly from Theorem 3 above as follows. We first note that for a convolutional generator (where U implements a convolution and thus is circulant) the eigenvectors of the matrix Σ(U) are given by the trigonometric basis functions per Definition 1 and the eigenvalues are the square of the entries of the dual kernel (Definition 2). To see this, we note that as detailed in Section H,

EQUATION

)

Because the matrix U implements a convolution with a kernel u that is equal to its first column, the matrix Σ(U) is again a circulant matrix and is also Hermitian. Thus, its spectrum is given by the Fourier transform of the first column of the circulant matrix, and its left-singular vectors are given by the trigonometric basis functions defined in equation ( 5) and depicted in Figure 3 .

Furthermore, using the fact that the eigenvalues of a circulant matrix are given by its discrete Fourier transform we can substitute β = U = Fu ∞ and α = σ n > 0. With the assumption τ ≤ T = 100 ησ 2 p , and using that

1 + ξ 4 α β T β 2 ≤ 1 + 32ξ β σ p ≤ 33,

where we used the assumption ξ ≤ σp β , the condition ( 10) is implied by the assumption k ≥ C u n/ξ 8 , with C u ∝ β 18 α 18 (see equation ( 6), where this assumption is made). This yields Theorem 2.

D.2 PROOF OF THEOREM 1:

Finally, we note that to obtain the simplified final expression in Theorem 1 from Theorem 2 we also used the fact that for a Gaussian vector z, the vector W T z is also Gaussian. Furthermore, by the concentration of Lipschitz functions of Gaussians with high probability we have

n i=1 ((1 -ησ 2 i ) τ -1) 2 w i , z 2 ≈ E n i=1 ((1 -ησ 2 i ) τ -1) 2 w i , z 2 (i) = ς 2 n n i=1 ((1 -ησ 2 i ) τ -1) 2 (ii) ≤ ς 2 2p n .

Here, equation (i) follows from w i , z being zero mean Gaussian with variance ς 2 /n (since z ∼ N (0, (ς 2 /n)I), and w i 2 = 1). Finally, (ii) follows by choosing the early stopping time so that it obeys τ ≤ log(1

-p/n)/ log(1 -ησ 2 p+1 ) which in turn implies that (1 -ησ 2 i ) τ ≥ 1 -p/n, for all i > p, yielding that ((1 -ησ 2 i ) τ -1) 2 ≤ p/n, for all i > p.

E THE DYNAMICS OF LINEAR AND NONLINEAR LEAST-SQUARES:

Theorem 3 builds on a more general result on the dynamics of a non-linear least squares problem which applies beyond convolutional networks, and that is stated and discussed in this section. Consider a nonlinear least-squares fitting problem of the form

L(θ) = 1 2 f (θ) -y 2 2 .

Here, f : R N → R n is a non-linear model with parameters θ ∈ R N . To solve this problem, we run gradient descent with a fixed stepsize η, starting from an initial point θ 0 , with updates of the form

EQUATION

Here, J (θ) ∈ R n×N is the Jacobian associated with the nonlinear map f with entries given by [J (θ)] i,j = ∂fi(θ) ∂θj . In order to study the properties of the gradient descent iterates ( 14), we relate the non-linear least squares problem to a linearized one in a ball around the initialization θ 0 . We note that this general strategy has been utilized in a variety of recent publications (Du et al., 2018; Arora et al., 2019; Oymak & Soltanolkotabi, 2019b; Oymak et al., 2019) . The associated linearized least-squares problem is defined as

EQUATION

Here, J ∈ R n×p , refered to as the reference Jacobian, is a fixed matrix independent of θ that approximates the Jacobian mapping at initialization, J (θ 0 ). Starting from the same initial point θ 0 , the gradient descent updates of the linearized problem are

EQUATION

To show that the non-linear updates ( 14) are close to the linearized iterates ( 16), we make the following assumptions: Assumption 1 (Bounded spectrum). We assume the singular values of the reference Jacobian obey

EQUATION

Furthermore, we assume that the Jacobian mapping associated with the nonlinear model f obeys

EQUATION

)

Assumption 2 (Closeness of the reference and initialization Jacobians). We assume the reference Jacobian and the Jacobian of the nonlinearity at initialization J (θ 0 ) are 0 -close in the sense that

EQUATION

Assumption 3 (Bounded variation of Jacobian around initialization). We assume that within a radius R around the initialization, the Jacobian varies by no more than in the sense that

EQUATION

where B R (θ 0 ) := {θ : θθ 0 ≤ R} is the ball with radius R around θ 0 .

Our first result shows that under these assumptions the nonlinear iterative updates ( 14) are intimately related to the linear iterative updates (16). Specifically, we show that the residuals associated with these two problems defined below nonlinear residual:

EQUATION

linear residual:

EQUATION

are close in the proximity of the initialization. Theorem 4 (Closeness of linear and nonlinear least-squares problems). Assume the Jacobian mapping J (θ) ∈ R n×N associated with the function f (θ) obeys Assumptions 1, 2, and 3 around an initial point θ 0 ∈ R N with respect to a reference Jacobian J ∈ R n×N and with parameters α, β, 0 , , and R. Furthermore, assume the radius R is given by R

EQUATION

with T a constant obeying 1 ≤ T ≤ 1 2η 2 , and J † the pseudo-inverse of J. We run gradient descent with stepsize η ≤ 1 β 2 on the linear and non-linear least squares problem, starting from the same initialization θ 0 . Then, for all τ ≤ T the iterates of the original and the linearized problems and their corresponding residuals obey

EQUATION

)

Moreover, for all iterates τ ≤ T ,

EQUATION

)

The above theorem formalizes that in a (small) radius around the initialization, the non-linear problem behaves similarly as its linearization. Thus, to characterize the dynamics of the nonlinear problem, it suffices to characterize the dynamics of the linearized problem. This is the subject of our next theorem.

Theorem 5. Consider a linear least squares problem (15) and let J = WΣV T ∈ R n×N = n i=1 σ i w i v T i be the singular value decomposition of the matrix J. Then the residual r τ after τ iterations of gradient descent with updates (16) is

EQUATION

Moreover, using a step size satisfying η ≤ 1

σ 2 1

, the linearized iterates (16) obey

EQUATION

)

In the next section we combine these two general theorems to provide guarantees for denoising using general neural networks.

E.1 PROOF OF THEOREM 4 (CLOSENESS OF LINEAR AND NON-LINEAR LEAST-SQUARES):

The proof is by induction. We suppose the statement, in particular the bounds ( 24), ( 25), and ( 26) hold for iterations t ≤ τ -1. We then show that they continue to hold for iteration τ in four steps.

In Step I, we show that a weaker version of (26) holds, specifically that θ τ -θ 0 2 ≤ R. Using this result, in Steps II and III we show that the bounds ( 24) and ( 25) hold, respectively. Finally, in

Step IV we utilize Steps I-III to complete the proof of equation ( 26).

Step I: Next iterate obeys θ τ ∈ B R (θ 0 ). To prove θ τ ∈ B R (θ 0 ), first note that by the triangle inequality and the induction assumption (26) we have

θ τ -θ 0 2 ≤ θ τ -θ τ -1 2 + θ τ -1 -θ 0 2 , ≤ θ τ -θ τ -1 2 + R 2 .

So to prove θ τ -θ 0 2 ≤ R it suffices to show that θ τ -θ τ -1 2 ≤ R/2. To this aim note that

1 η θ τ -θ τ -1 2 = ∇L(θ τ -1 ) 2 = J T (θ τ -1 )r τ -1 2 ≤ J T (θ τ -1 ) r τ -1 2 + J (θ τ -1 ) r τ -1 -r τ -1 2 ≤ J T r τ -1 2 + J (θ τ -1 ) -J r τ -1 2 + J (θ τ -1 ) r τ -1 -r τ -1 2 (i) ≤ J T r τ -1 2 + ( + 0 ) r τ -1 2 + β r τ -1 -r τ -1 2 (ii) ≤ 1 η J † r 0 2 + ( + 0 ) r 0 2 + 2 β 2 α 2 ( 0 + ) r 0 2 ≤ R 2 .

In the above (i) follows from Assumptions 1, 2, and 3. For (ii), we bounded the last term with the induction hypothesis (24), the middle term with r τ -1 2 ≤ r 0 2 , and the first term with the bound

J T r τ -1 2 = J T (I -ηJJ T ) τ -1 r 0 2 = Σ(I -ηΣ 2 ) τ -1 W T r 0 2 ≤ n j=1 σ 2 j w j , r 0 2 ≤ β 2 n j=1 1 σ 2 j w j , r 0 2 = β 2 J † r 0 2 ≤ 1 η J † r 0 2 .

Finally, the last inequality follows by definition of R in ( 23) together with the fact that T ≥ 1.

Step II: Original and linearized residuals are close: In this step, we bound the deviation of the residuals of the original and linearized problem defined as

e τ := r τ -r τ .

This step relies on the following lemma, which is a variant of (Oymak et al., 2019, Lem. 6.7) .

Lemma 1 (Bound on growth of perturbations). Suppose that Assumptions 1, 2, and 3 hold and that θ τ , θ τ +1 ∈ B R (θ 0 ). Then, provided the stepsize obeys η ≤ 1/β 2 , the deviation of the residuals obeys

EQUATION

By the previous step, θ τ , θ τ +1 ∈ B R (θ 0 ). We next bound the two terms on the right hand side of (29). Regarding the first term, we note that an immediate consequence of Theorem 5 is the following bound on the residual of the linearized iterates:

EQUATION

)

In order to bound the second term in (29), namely, e τ 2 , we used the following lemma, proven later in Section E.1.1. Lemma 2. Suppose that for positive scalars α, η, ρ, ξ > 0, η ≤ 1/α 2 , the sequences r τ and e τ obey

EQUATION

Then, for all τ ≤ 1 2η 2 , we have that

EQUATION

)

With these lemmas in place we now have all the tools to prove that the original and linear residuals are close. In particular, from Step I, we know that θ τ ∈ B R (θ 0 ) so that the assumptions of Lemma 1 are satisfied. Lemma 1 implies that the assumption (32) of Lemma 2 is satisfied with ξ = β( 0 + ) and the bound (30) implies that the assumption (31) of Lemma 2 is satisfied with ρ = r 0 2 . Thus, Lemma 2 implies that for all τ ≤ 1 2η 2

EQUATION

)

This concludes the proof of (24).

Step III: Original and linearized parameters are close: First note that by the triangle inequality and Assumptions 2 and 3 we have

EQUATION

)

The difference between the parameter of the original iterate θ and the linearized iterate θ obey

1 η θ τ -θ τ 2 = τ -1 t=0 ∇L(θ t ) -∇L lin ( θ t ) 2 = τ -1 t=0 J T (θ t )r t -J T r t 2 ≤ τ -1 t=0 (J T (θ t ) -J T ) r t 2 + J T (θ t )(r t -r t ) 2 (i) ≤ τ -1 t=0 ( 0 + ) r t 2 + β e t 2 (ii) ≤ ( 0 + ) τ -1 t=0 1 -ηα 2 τ -1 r 0 2 + β τ -1 t=0 e t 2 = ( 0 + ) 1 -(1 -ηα 2 ) τ ηα 2 r 0 2 + β τ -1 t=0 e t 2 (iii) ≤ ( 0 + ) ηα 2 r 0 2 + 2τ β 2 α 2 ( 0 + ) r 0 2 = 1 ηα 2 ( 0 + ) 1 + 2ητ β 2 r 0 2 .

Here, inequality (i) follows from ( 35) and an application of Assumption 1, inequality (ii) from (30), inequality (iii) from η ≤ 1/β 2 which implies (1 -ηα 2 ) ≥ 0 and (34). This concludes the proof of the bound (25).

Step IV: Completing the proof of (26): By the triangle inequality

θ τ -θ 0 2 ≤ θ τ -θ 0 2 + θ τ -θ τ 2 (i) ≤ J † r 0 2 + 1 α 2 ( 0 + ) 1 + 2ητ β 2 r 0 2 (ii) ≤ R/2,

where inequality (i) follows from the bound (25), which we just proved, and the fact that, from equation (36) in Theorem 5,

θ τ -θ 0 2 2 = n i=1 w i , r 0 2 (1 -(1 -ησ 2 i ) τ ) 2 σ 2 i ≤ n i=1 w i , r 0 2 /σ 2 i = J † r 0 2 .

Finally, inequality (ii) follows from the definition of R in equation ( 23) along with τ ≤ T , by assumption.

E.1.1 PROOF OF LEMMA 2:

We prove the result by induction. Assume equation ( 33) holds true for some τ . We prove that then it also holds true for τ + 1. By the two assumptions in the lemma,

e τ +1 -e τ ≤ η 2 e τ + ηξ r τ ≤ η 2 e τ + ηξ(1 -ηα 2 ) τ ρ (i) ≤ ηξ 2 2ρ α 2 + (1 -ηα 2 ) τ ρ ,

where (i) follows from the induction assumption (33). Summing up the difference of the errors gives

e τ ξ = τ -1 t=0 (e t+1 -e t ) ξ ≤ τ η 2 2ρ α 2 + ηρ τ -1 t=0 (1 -ηα 2 ) t = τ η 2 2ρ α 2 + ηρ 1 -(1 -ηα 2 ) τ ηα 2 ≤ ρ α 2 (η2τ 2 + 1) ≤ 2 ρ α 2 ,

where the last inequality follows from the assumption that τ ≤ 1 2η 2 .

E.2 PROOF OF THEOREM 5:

The proof of identity ( 27) is equivalent to the proof of Proposition (1). Regarding inequality (28), note that

EQUATION

This in turn implies that

θ τ -θ 0 2 2 = n i=1 v i , θ τ -θ 0 2 = n i=1 w i , r 0 1 -(1 -ησ 2 i ) τ σ i 2 .

F PROOFS FOR NEURAL NETWORK GENERATORS (PROOF OF THEOREM 3):

The proof of Theorem 3 relies on the fact that, in the overparameterized regime, the non-linear least squares problem is well approximated by an associated linearized least squares problem. Studying the associated linear problem enables us to prove the result. We apply Theorem 4, which ensures that the associated linear problem is a good approximation of the non-linear least squarest problem, with the nonlinearity G(C) = ReLU(UC)v and with the parameter given by θ = C. Recall that v is a fixed vector with half of the entries 1/ √ k, and the other half -1/ √ k.

As the reference Jacobian in the associated linear problem, we choose a matrix J ∈ R n×nk , that is very close to the original Jacobian at initialization, J (C 0 ), and that obeys JJ T = Bound on initial residual: For what follows, we require a bound on the initial residual. To prove this bound, we apply the following lemma. Lemma 6 (Initial residual). Consider G(C) = ReLU(UC)v, and let C ∈ R n×k be generated at random with i.i.d. N (0, ω 2 ) entries. Suppose half of the entries of v are ν/ √ k and the other half are -ν/ √ k, for some constant ν > 0. Then, with probability at least 1 -δ, G(C) 2 ≤ νω 8 log(2n/δ) U F .

Now, the initial residual can be upper bounded as

EQUATION

where we used that, by Lemma 6,

EQUATION

Here, the second inequality follows from U F ≤ √ n U and our choice of the parameter ω, and the last inequality follows from ξ ≤ 1/ 32 log(2n/δ) and α/β ≤ 1.

Verifying Assumption 3: Verification of the assumption requires us to control the perturbation of the Jacobian matrix around a random initialization. Lemma 7 (Jacobian perturbation around initialization). Let C 0 be a matrix with i.i.d. N (0, ω 2 ) entries. Then, for all C obeying

C -C 0 ≤ ω R with R ≤ 1 2 √ k,

the Jacobian mapping (9)associated with the generator G(C) = ReLU(UC)v obeys

J (C) -J (C 0 ) ≤ v ∞ 2(k R) 1/3 U , with probability at least 1 -ne -1 2 R 4/3 k 7/3 .

In order to verify Assumption 3, first note that the radius in the theorem, defined in equation ( 23), obeys

R = 2 J † r 0 2 + 2 α 2 ( 0 + ) 1 + 2ηT β 2 r 0 2 (i) ≤ 2 α + 2 α 2 ( 0 + ) 1 + 2ηT β 2 r 0 2 (ii) ≤ 2 α 1 + 2 α 1 + 2ηT β 2 2ω √ nξ -1 β β 2 α 2 = 4 1 + 2 α 1 + 2ηT β 2 ω √ nξ -1 β 3 α 3 (iii) ≤ 4ω 2 + ξ 2 α β ηT β 2 √ nξ -1 β 3 α 3 (iv) ≤ ω ξ 16 α 2 β 2 3 √ k := ω R.

Here, (i) follows from the fact that J † S r 0 2 ≤ 1 α r 0 2 , (ii) from the bound on the initial residual (38), using that 0 ≤ , as well as our choice of ω. Moreover (iii) follows from the definition of and (iv) from the lower bound on k in assumption (10).

Note that since ξ 2 16 α 2 β 2 ≤ 1 2 for this choice of radius R, by Lemma 7 we have, with

v ∞ = 1/ √ k, that J (C) -J (C 0 ) ≤ v ∞ 2(k R) 1/3 U = 1 √ k 2 ξ 16 α 2 β 2 (k • k 1/2 ) 1/3 β = ξ 8 α 2 β = .

holds with probability at least

1 -ne -1 2 Rk 7/3 = 1 -ne -2 -17 ξ 4 α 8 β 8 k 3 (i) ≥ 1 -δ,

where in (i) we used (10) together with ξ ≤ 8 log 2n δ . Therefore, Assumption 3 holds with high probability by our choice of = ξ 8 α 2 β .

Verifying the bound on number of iterations: Finally we note that the constraint on the number of iterations, T , T ≤ 1 2η 2 , from Theorem 4 is satisfied under the number of constraints of Theorem 3 by 1 2η 2 = 2 5 β 2 ηξ 2 α 4 .

Concluding the proof of Theorem 3: Now that we have verified the conditions of Theorem 4 we can apply the theorem. This allows us to conclude that

G(C τ ) -x 2 = G(C τ ) + z -y 2 = r τ + z 2 = r τ + z + r τ -r τ 2 (i) ≤ r τ + z 2 + r τ -r τ 2 (ii) ≤ r τ + z 2 + ξ y 2 (iii) = W I -ηΣ 2 τ W T r 0 + z 2 + 1 2 ξ y 2 (iv) = W I -ηΣ 2 τ W T (G(C 0 ) -x) -W I -ηΣ 2 τ W T -I z 2 + 1 2 ξ r 0 2 (v) ≤ I -ηΣ 2 τ W T x 2 + I -ηΣ 2 τ -I W T z 2 + G(C 0 ) 2 + 1 2 ξ y 2 (vi) ≤ (1 -ησ 2 p ) τ x 2 + n i=1 ((1 -ησ 2 i ) τ -1) 2 w i , z 2 + ξ y 2 .

Here, (i) follows from the triangular inequality, (ii) from Theorem 4 equation ( 24) combined with our choice for and r 0 2 ≤ 1 2 y 2 , shown above. Moreover, (iii) follows from Theorem 5, (iv) from r 0 = G(C 0 ) -y = G(C 0 ) -xz, (v) from the triangular inequality. Finally, for (vi) we used the fact that x ∈ span{w 1 , w 2 , . . . , w p } combined with the fact that W I -ηΣ 2 τ W T ≤ 1, as well as using the bound G(C 0 ) 2 ≤ 1 2 y 2 proven in (39). This proves the final bound (12), as desired.

G PROOF OF PROPOSITION 1 AND EQUATION (2):

The residual of gradient descent at iteration τ is

r τ = y -Jc τ = y -J(c τ -1 -ηJ T (Jc τ -1 -y)) = (I -ηJJ T )(y -Jc τ -1 ) = (I -ηJJ T ) τ (y -Jc 0 ) = (I -ηJJ T ) τ y = (I -ηWΣ 2 W T ) τ y

where we used that c 0 = 0 and the SVD J = WΣV T . Expanding y in terms of the singular vectors w i (i.e., the columns of W), as y = i w i w i , y , and noting that (I -

ηWΣ 2 W T ) τ = i (1 -ησ 2 i ) τ w i w T i we get r τ = i (1 -ησ 2 i ) τ w i w i , y , as desired.

Proof of equation (2): By proposition 1 and using that x and lies in the signal subspace

x -Jc τ = p i=1 w i (1 -ησ 2 i ) τ w i , x + n i=1 w i ((1 -ησ 2 i ) τ -1) w i , z .

By the triangle inequality,

x -Jc τ 2 ≤ p i=1 w i (1 -ησ 2 i ) τ w i , x 2 + n i=1 w i ((1 -ησ 2 i ) τ -1) w i , z 2 ≤ (1 -ησ 2 p ) τ x 2 + n i=1 ((1 -ησ 2 i ) τ -1) 2 w i , z 2 ,

where the second inequality follows by using orthogonality of the w i and by using (1 -ησ 2 i ) ≤ 1, from η ≤ 1/σ 2 max . This concludes the proof of equation ( 2).

H THE EXPECTED JACOBIAN FOR CONVOLUTIONAL GENERATORS:

We first prove the closed form expression of the expected Jacobian, or more precisely of the matrix Σ(U) defined in (13). By the expression for the Jacobian in equation ( 9), we have that (40)

J (C)J T (C) =

Next, we have with (Daniely et al., 2016, Sec. 4 .2) and using that the derivative of the ReLU function is the step function, E σ (Uc )σ (Uc )

T ij = 1 2 1 -cos -1 u i , u j u i 2 u j 2 /π .

Using that v 2 = 1, we get

E J (C)J T (C) ij = 1 2 1 -cos -1 u i , u j u i 2 u j 2 /π u i , u j ,

where u i are the rows of U. This concludes the proof of equation ( 13).

We next briefly comment on the singular value decomposition of a circulant matrix and explain that the singular vectors are given by Definition 1. Recall, that U ∈ R n×n is a circulant matrix, implementing the convolution with a filter. Assume for simplicity that n is even. It is well known that the discrete Fourier transform diagonalizes U, i.e., U = F -1 ÛF, where F ∈ C n×n is the DFT matrix with entries

[F] jk = e i2πjk/n , j, k = 0, . . . , n -1, and Û is a diagonal matrix with diagonal Fu, where u is the first column of the circulant matrix U. From this, we can compute the singular value decomposition of U by using that û = Fu is conjugate symmetric (since u is real) so that [û] n-k+2 = û * , k = 2, . . . , n/2. Let U = WΣV T be the singular value decomposition of U. The entries of the left singular vectors are given by the trigonometric basis functions defined in (5), and the singular values are given by the absolute values of û. The proof follows that of (Oymak & Soltanolkotabi, 2019b, Lem. 6.9) .

Step 1: We start by relating the perturbation of the Jacobian to a perturbation of the activation patterns. For any C, C , we have that

EQUATION

)

To see this, first note that, by (9), 

≤ U 2 max j (σ (u T j C) -σ (u T j C ))diag(v) 2 ≤ v 2 ∞ U 2 max j σ (u T j C) -σ (u T j C ) 2 ,

where for (i) we used that for two positive semidefinite matrices A, B, λ max (A B) ≤ λ max (A) max i B ii . This concludes the proof of equation (41).

Step 2:

Step one implies that we only need to control σ (Uc j ) around a neighborhood of c j . Since σ is the step function, we need to count the number of sign flips between the matrices UC and UC . Let |v| π(q) be the q-th smallest entry of v in absolute value. Lemma 8. Suppose that, for all i, and q ≤ k,

C -C ≤ √ q |u T i C | π(q) u i .

Then

max i σ (u T i C) -σ (u T i C ) ≤ 2q

Proof. Suppose that σ (u T i C) and σ (u T i C ) have 2q many different entries, then the conclusion of the statement would be violated. We show that this implies the assumption is violated as well, proving the statement by contraction. By the contradiction hypothesis,

C -C 2 ≥ u T i (C -C ) 2 / u i 2 ≥ q |u T i C | 2 π(q) u i 2 ,

where the last inequality follows by noting that at least 2q many entries have different signs, thus their difference is larger than their individual magnitudes, and at least q many individual magnitudes are lower bounded by the q-th smallest one.

Step 3: Next, we note that, with probability at least 1 -ne -kq 2 /2 , the q-th smallest entry of

u T i C ∈ R k obeys |u T i C | π(q) u i ≥ q 2k

ν for all i = 1, . . . , n.

With this notation,

J (C)J T (C) -Σ( U) = k =1 v 2 σ ( Uc )σ ( Uc ) T -E σ ( Uc )σ ( Uc ) T U U T = k =1 S .

To show concentration we use the matrix Hoeffding inequality. To this aim note that the summands are centered in the sense that E S = 0. Next note that

σ ( Uc )σ ( Uc ) T -E σ ( Uc )σ ( Uc ) T U U T σ ( Uc )σ ( Uc ) T U U T =diag σ ( Uc ) T U U T diag σ ( Uc ) B 2 U U T Similarly, σ ( Uc )σ ( Uc ) T -E σ ( Uc )σ ( Uc ) T U U T -E σ ( Uc )σ ( Uc ) T U U T -B 2 U U T .

Thus,

-v 2 B 2 U U T S v 2 B 2 U U T .

Therefore, using A := v 2 B 2 U U T we have

S 2 A 2 ,

and

σ 2 := k =1 A 2 ≤ B 4 k =1 v 4 U 4

To continue we will apply matrix Hoeffding inequality stated below. Theorem 6 (Matrix Hoeffding inequality, Theorem 1.3 (Tropp, 2011) ). Consider a finite sequence S of independent, random, self-adjoint matrices with dimension n, and let {A } be a sequence of fixed self-adjoint matrices. Assume that each random matrix satisfies E[S ] = 0 and S 2 A 2 almost surely.

Then, for all t ≥ 0,

P k =1 S ≥ t ≤ 2ne -t 2 8σ 2

where

σ 2 := k =1 A 2 .

Therefore, applying matrix Hoeffding inequality we get that Without loss of generality we prove the result for ν = 1. First note that by the triangle inequality r 0 2 = σ(UC)v -y 2 ≤ σ(UC)v 2 + y 2 .

:



Figures and Tables:

Figure 1: Fitting an over-parameterized Deep Decoder (DD-O) and the deep image prior (DIP) to a (a) noisy image, (b) clean image, and (c) pure noise. Here, MSE denotes Mean Square Error of the network output with respect to the clean image in (a) and fitted images in (b) and (c). While the network can fit the noise due to over-parameterization, it fits natural images in significantly fewer iterations than noise. Hence, when fitting a noisy image, the image component is fitted faster than the noise component which enables denoising via early stopping.

Figure 2: Denoising with BM3D and various convolutional generators. The relative ranking of algorithms in this picture is representative and maintained on a larger set of test images (see Appendix A). DD-U is an under-parameterized deep decoder. DIP and DD-O are over-parameterized convolutional generators and with early stopping outperform the BM3D algorithm, the next best method that does not require training data.

Figure 3: The 1st, 2nt, 6th, and 21th trigonometric basis functions in dimension n = 300.

Figure 4: Left panel: Convolutional generators. The output is generated through repeated convolutional layers, channel normalization, and applying ReLU non-linearities. Right panel: Fitting the phantom MRI and noise with different architectures of depth d = 5, for different number of over-parameterization factors (1,4, and 16). Gradient descent on convolutional generators involving fixed convolutional matrixes fit an image significantly faster than noise.

Figure 5: Gradient descent on the least squares problem of minimizing y -Jc τ 2 , where J ∈ R 100×100 has decaying singular values (left panel) and the observation is the sum of a signal component, equal to the leading singular vector w 1 of J, and a noisy component z ∼ N (0, (1/n)I), i.e., y = w 1 + z. The signal component w 1 is fitted significantly faster than the other components (right panel), thus early stopping enables denoising.

Figure 6: Triangular and Gaussian kernels and the weights associated to low-frequency trigonometric functions they induce, for a generator network of output dimension n = 300. The wider the kernels are, the more the weights are concentrated towards the low-frequency components of the signal.

Figure 7: The Singular value distribution of the Jacobian of a four-layer deep decoder after t = 50 and t = 3500 iterations of gradient descent (panel (a)), along with the corresponding singular vectors/function (b-i). The singular functions corresponding to the large singular vectors are close to the low-frequency Fourier modes and do not change significantly through training.

Figure 8: Fitting the phantom MRI and noise with different architectures of depth d = 5, for different number of over-parameterization factors (1,4, and 16). Gradient descent on convolutional generators involving fixed convolutional matrixes fit an image significantly faster than noise.

Figure9: The relative distances of the weights in each layer from its random initialization. The weights need to change significantly more to fit the noise, compared to an image, thus a natural image lies closer to a random initialization than noise.

Figure 11: Fitting a step function and noise with a two-layer deep decoder: Even for a two-layer network, the simple image (step function) is fitted significantly faster than the noise.

(σ (Uc ))UU T diag(σ (Uc )) where denotes the entrywise product of the two matrices. Then,E J (C)J T (C) = k =1 v 2 E σ (Uc )σ (Uc ) T UU T .

PROOFS OF LEMMAS FOR NEURAL NETWORK DENOISERS (PROOFS OF AUXILIARY LEMMAS IN SECTION F) I.1 PROOF OF LEMMA 7: JACOBIAN PERTURBATION AROUND INITIALIZATION

(C) -J (C ) = [. . . v j (diag(σ (Uc j )) -diag(σ (Uc j )) )U . . .].This in turn implies thatJ (C) -J (C ) 2 = (J (C) -J (C ))(J (C) -J (C )) T = (σ (UC) -σ (UC ))diag(v 2 1 , . . . , v 2 k )(σ (UC) -σ (UC ))

PROOF OF LEMMA 6 (BOUND ON INITIAL RESIDUAL)

Basri et al. (2019) have shown that functions of different frequencies are learned at different speeds, andArora et al. (2019) has provided a theoretical explanation of the empirical observation that a simple 2-layer network fits random labels slower than actual labels in the context of classification. A recent publication(Li et al., 2019) focuses on demonstrating how early stopping leads to robust classification in the presence of label corruption under a cluster model for the input data. Neither of the aforementioned publication however, does address denoising in a regression setting or fitting convolutional generators of the form studied in this paper. K. Zhang, W. Zuo, Y. Chen, D. Meng, and L. Zhang. Beyond a Gaussian denoiser: Residual learning of deep CNN for image denoising. IEEE Transactions on Image Processing, 26(7):3142-3155, 2017. Kai Zhong, Zhao Song, Prateek Jain, Peter L Bartlett, and Inderjit S Dhillon. Recovery guarantees for one-hidden-layer neural networks. In International Conference on Machine Learning, 2017.

Table content: <table/>

Table1below show that the overparameterized deep decoder with early stopping performs best for this task. Average performance for denoising images color image with the same Gaussian noise added to each color channel with BM3D and un-trained convolutional neural networks. The overparameterized deep decoder performs best for this task.

Table content: <table><tr><td>Noise level</td><td>20.60</td></tr><tr><td>Deep decoder with k=512 and early stopping</td><td>28.08dB</td></tr><tr><td colspan="2">deep decoder with k=128 without early stopping 27.84dB</td></tr><tr><td>DIP</td><td>25.76 dB</td></tr><tr><td>BM3D</td><td>25.52 dB</td></tr></table>

Peer review (iclr2020):
The paper provides a theoretical study of regularization capabilities of over-parameterized convolutional generators trained via gradient descent, in the context of denoising with an approach similar to the *deep image prior*. The authors show that when using appropriate upsampling operators, gradient descent biases the reconstructed images towards low-frequency components, while the noise components, which typically consist of high-frequency patterns, take longer to fit, so that early stopping can provide a useful bias for denoising. The proofs are based on recent results for *lazy* training of over-parameterized networks, namely neural tangent kernels. I find the paper interesting and novel, as it is the first to my knowledge to study the implicit bias of optimization for such generator-based denoising approaches, and provides interesting theoretical and empirical denoising results. While the current results are for a simple architecture with one hidden layer, they may pave the way for a formal study of more complex architectures. Nevertheless, the paper presents some limitations which should be addressed further in the paper: * in terms of denoising capabilities, it seems that the results for the non-linear, over-parameterized network are essentially similar to the linear case (with an appropriately designed operator J in Section 3). Are there any benefits to the non-linear case that I am missing? Whether or not this is a limitation of the study, it should be discussed further in the paper. * Some comments on the empirical results reported in Figure 2: (i) how are hyper-parameters chosen? this seems to be crucial for the effective use of such denoising strategies, and it is unclear how robust these methods are, e.g., to the early stopping time; (ii) there are other methods that only learn on the given image (thus without the need for a training set) which outperform BM3D [e.g. 1,2], it is unclear how the proposed methods compare to those. other minor comments: - Figure 2 *maintained on a large test set of images* more details are welcome, including on hyperparameter selection - Section 1.1: * *demystifying* -> understanding? * different notations are used for the (same?) iterates, C^t and C_tau - Section 2.3 * *it follow that optimizing...*: please clarify. I agree that the optimum is the same, but is the implicit bias from initialization the same? - end of Section 3: early stoped -> stopped - Section 4 * after Definition 2, *this in turn implies that the jacobian spectrum ...rapidly*: clarify, is it a consequence of the analysis? * after Thm 1: maybe give some insight about the proof and the use of NTK? also, discuss similarities or differences with the linear case [1] Dabov et al. (2009) BM3D Image Denoising with Shape-Adaptive Principal Component Analysis [2] Mairal et al. (2009) Non-local Sparse Models for Image Restoration