Year: 2020

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.

text: 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.

heading: 1 INTRODUCTION
text: 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‖22. 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).

heading: 1.1 CONTRIBUTIONS AND OVERVIEW OF RESULTS
text: 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.
1 2 6 21

heading: 2 CONVOLUTIONAL GENERATORS
text: A convolutional generator maps an input tensor B0 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.

heading: 2.1 THE IMPORTANCE OF FIXED CONVOLUTIONAL FILTERS
text: 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.

heading: 2.2 ARCHITECTURE OF CONVOLUTIONAL GENERATOR WITH FIXED CONVOLUTIONS
text: 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
Bi+1 = cn(ReLU(UiBiCi)), i = 0, . . . , d− 1, and finally, the output of the d-layer network is formed as
x = BdCd+1.
Here, the coefficient matrices Ci ∈ Rk×k and Cd+1 ∈ Rk×kout contain the weights of the network. The number of channels, k, determines the number of weight parameters of the network, given by dk2 + koutk. Each column of the tensor BiCi ∈ Rni×k is formed by taking linear combinations of the channels of the tensor Bi 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 Ui ∈ Rni+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.

heading: 2.3 TWO LAYER CONVOLUTIONAL GENERATOR STUDIED THEORETICALLY IN THIS PAPER
text: 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) ∈ Rn. Then, the generator becomes
G(C) = ReLU(UB0C0)c1,
where U ∈ Rn×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 C0 with a non-negativ factor is equivalent to scaling the i-th entry of the output weights c1 with the same factor. We therefore fix the output weights to c1 = 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 B0 is Gaussian, with probability one, it has full rank and thus spans Rn. It follows that optimizing over the parameter C0 in C = B0C0 is equivalent to optimizing over the matrix C ∈ Rn×k. We therefore consider the generator
G(C) = ReLU(UC)v, (1)
where v = [1, . . . , 1,−1, . . . ,−1]/ √ k is fixed and C ∈ Rn×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.

heading: 3 WARMUP: DYNAMICS OF GRADIENT DESCENT ON LEAST SQUARES
text: 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‖22,
and study gradient descent with a constant step size η starting at c0 = 0. The updates are given by
cτ+1 = cτ − η∇L(cτ ), ∇L(c) = JT (Jc− y). The following simple proposition characterizes the trajectory of gradient descent. Proposition 1. Let J ∈ Rn×m be a matrix with left singular vectors w1, . . . ,wn ∈ Rn and corresponding singular values σ1 ≥ σ2 ≥ . . . ≥ σn. Then the residual after τ steps, rτ = y − Jcτ , of gradient descent starting at c0 = 0 is
rτ = n∑ i=1 wi 〈wi,y〉 (1− ησ2i )τ .
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
y = x + z,
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
‖Jcτ − x‖2 ≤ (1− ησ 2 p) τ ‖x‖2 + E(z), E(z) := √√√√ n∑ i=1 ((1− ησ2i )τ − 1)2 〈wi, 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− ησ2i )τ − 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.

heading: 4 DYNAMICS OF GRADIENT DESCENT ON CONVOLUTIONAL GENERATORS
text: 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 ∈ Rn×k, and recall that U is a circulant matrix implementing a convolution with a kernel u ∈ Rn. We fit the generator to a signal y ∈ Rn by minimizing the non-linear least squares objective
L(C) = 1 2 ‖y −G(C)‖22 (3)
with (early-stopped gradient) descent with a constant stepsize η starting at a random initialization C0 of the weights. The iterates are given by
Cτ+1 = Cτ − η∇L(Cτ ). (4)
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 J (C) := ∂∂CG(C).
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.
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 w1, . . . ,wn are defined as
[wi+1]j = 1√ n  1 i = 0√ 2 cos(2πji/n) i = 1, . . . , n/2− 1 (−1)j i = n/2√
2 sin(2πji/n) i = n/2 + 1, . . . , n− 1
. (5)
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 ∈ Rn we define the dual kernel σ ∈ Rn as
σ = ‖u‖2 √√√√∣∣∣∣∣Fg ( u ~ u
‖u‖22
)∣∣∣∣∣ with g(t) = 12 ( 1− cos −1 (t) π ) t.
Here, for two vectors u,v ∈ Rn, 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 ∈ Rn 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 G(C) = ReLU(UC)v with,
k ≥ Cun/ 8, (6)
channels, for some > 0, and with convolutional kernel u of the convolutional operator U and associated dual kernel σ, to the noisy signal y. Here, Cu a constant only depending on the convolutional kernel u. Then, with probability at least 1−e−k2− 1n2 , the reconstruction error obtained after τ = log(1 − √ p/n)/ log(1 − ησ2p+1) iterations of gradient descent (4) with step size η ≤ 1‖Fu‖2∞ (Fu is the Fourier transform of u) starting from C0 with i.i.d. N (0, ω2), entries, ω ∝ ‖y‖2√ n
, is bounded by
‖G(Cτ )− x‖2 ≤ (1− ησ 2 p) τ ‖x‖2 + ς
√ 2p
n + ‖y‖2. (7)
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−ησ2p)τ ≈ 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ησ2p and provided that k ≥ Cun/ 8, for some ∈ (0, σpσ1 ], with probability at least 1− e−k2 − 1n2 , the reconstruction error obeys
‖G(Cτ )− x‖2 ≤ (1− ησ 2 p) τ ‖x‖2 + √√√√ n∑ i=1 ((1− ησ2i )τ − 1)2 〈wi, 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 wi 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− ησ2p)τ is small), while most of the noise has not been fitted (i.e., ((1 − ησ2i )τ − 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):
Llin(c) = 1
2 ‖G(C0) + J(c− vect(C0))− y‖22.
Here, vect(C0) ∈ Rkn is a vectorized version of the matrix C0, and J ∈ Rn×nk approximates the Jacobian of G(C) around the initialization C0. 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.

heading: 4.1 MULTILAYER NETWORKS AND MODERATE OVER-PARAMETERIZATION
text: 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 R512 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× k2 = 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.

heading: 5 RELATED LITERATURE
text: 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
of the network, in all the aforementioned methods, the weights are adjusted only during training and then are fixed upon solving the inverse problem.
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: Basri et al. (2019) have shown that functions of different frequencies are learned at different speeds, and Arora 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.

heading: CODE AND ACKNOWLEDGEMENTS
text: Code to reproduce the experiments is available at https://github.com/MLI-lab/ overparameterized_convolutional_generators.
R. Heckel is partially supported by NSF award IIS-1816986, acknowledges support of the NVIDIA Corporation in form of a GPU, and would like to thank Tobit Klug for proofreading a previous version of this manuscript. M. Soltanolkotabi is supported by the Packard Fellowship in Science and Engineering, a Sloan Research Fellowship in Mathematics, an NSF-CAREER under award #1846369, the Air Force Office of Scientific Research Young Investigator Program (AFOSR-YIP) under award #FA9550-18-1-0078, an NSF-CIF award #1813877, DARPA under the Learning with Less Labels (LwLL) program, and a Google faculty research award.

heading: A DENOISING PERFORMANCE OF UNTRAINED CONVOLUTIONAL GENERATORS
text: 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. Table 1 below show that the overparameterized deep decoder with early stopping performs best for this task.

heading: B A NUMERICAL STUDY OF THE IMPLICIT BIAS OF CONVOLUTIONAL NETWORKS
text: 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, Bi+1, from the activations in the i-th layer, Bi:
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.
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) : Rn → Rn 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 : Rn → R2n be an upsampling operator, that in the one dimensional case transforms [x1, x2, . . . , xn] to [x1, 0, x2, 0, . . . , xn, 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: bi+1,` = ReLU (∑k j=1 M(cij`)bi ) , 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 Ci = {cij`}, and all coefficients of the networks are C = {cij`}. Note that here, the coefficients or parameters of the networks are the weights and not the input to the network.

heading: B.1 DEMONSTRATING IMPLICIT BIAS OF CONVOLUTIONAL GENERATORS
text: 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∗‖22 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,
(a) fitting images
(b) fitting pure noise
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.

heading: C THE SPECTRUM OF THE JACOBIAN OF THE DEEP DECODER AND DEEP IMAGE PRIOR
text: 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.
Towards this goal, we compute the norm of the product of the Jacobian at a random initalization, C0, with a signal y as this measures the extent to which the signal is aligned with the leading singular vectors, due to ∥∥J T (C0)y∥∥2 = ∥∥VΣWTy∥∥2 = ∑
i
σ2i 〈wi,y〉 2 ,
where J (C0) = WΣVT 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 (C0), 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 (C0)z∥∥) is significantly smaller than that with a natural image (i.e.,
∥∥J T (C0)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.

heading: D PROOFS AND FORMAL STATEMENT OF RESULTS
text: 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
G(C) = ReLU(UC)v, (8)
with C ∈ Rn×k, and U ∈ Rn×n an arbitrary fixed matrix, and v ∈ Rk, with half of the entries equal to +1/ √ k and the other half equal to −1/ √ k.
The (transposed) Jacobian of ReLU(Uc) is UTdiag(ReLU′(Uc)). Thus the Jacobian of G(C) is given by
J T (C) = v1U Tdiag(ReLU′(Uc1))
... vkU Tdiag(ReLU′(Uck))  ∈ Rnk×n, (9) 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) ] ∈ Rn×n
associated with the generatorG(C) = ReLU(UC)v. Here, expectation is over C with iidN (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 iidN (0, ω2) entries. Consider the eigenvalue decomposition of Σ(U) given by
Σ(U) = n∑ i=1 σ2iwiw T i .
Our results depend on the largest and smallest eigenvalue of Σ(U), defined throughout as
α = σ2n, β = σ 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 ∈ Rn given by
y = x + z,
where x ∈ Rn is assumed to lie in the signal subspace spanned by the p leading singular vectors w1, . . . ,wn of Σ(U), and z ∈ Rn is an arbitrary noise vector. Suppose that the number of channels obeys
k ≥ cnξ−8 ( 1 + ξ
4
α β ηTβ2
)2 β18
α18 (10)
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 minimiz-
ing the loss
L(C) = 1 2 ‖G(C)− y‖22 (11)
via running gradient descent starting from C0 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−k2 − δ, for all iterations τ ≤ T ,
‖x−G(Cτ )‖2 ≤ (1− ησ 2 p) τ‖x‖2 + √√√√ n∑ i=1 ((1− ησ2i )τ − 1)2 〈wi, z〉 2 + ξ‖y‖2. (12)

heading: D.1 PROOF OF THEOREM 2
text: 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,
[Σ(U)]ij = 1
2
( 1− cos−1 ( 〈ui,uj〉 ‖ui‖‖uj‖ ) /π ) 〈ui,uj〉 . (13)
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 ησ2p , 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 ≥ Cun/ξ 8, with Cu ∝ β 18
α18 (see equation (6), where this assumption is made). This yields Theorem 2.

heading: D.2 PROOF OF THEOREM 1
text: 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 WT z is also Gaussian. Furthermore, by the concentration of Lipschitz functions of Gaussians with high probability we have
n∑ i=1 ((1− ησ2i )τ − 1)2 〈wi, z〉 2 ≈ E [ n∑ i=1 ((1− ησ2i )τ − 1)2 〈wi, z〉 2 ]
(i) = ς2
n n∑ i=1 ((1− ησ2i )τ − 1)2
(ii) ≤ ς2 2p
n .
Here, equation (i) follows from 〈wi, z〉 being zero mean Gaussian with variance ς2/n (since z ∼ N (0, (ς2/n)I), and ‖wi‖2 = 1). Finally, (ii) follows by choosing the early stopping time so that it obeys τ ≤ log(1 − √ p/n)/ log(1 − ησ2p+1) which in turn implies that (1 − ησ2i )τ ≥ 1 − √ p/n, for all i > p, yielding that ((1− ησ2i )τ − 1)2 ≤ p/n, for all i > p.

heading: E THE DYNAMICS OF LINEAR AND NONLINEAR LEAST-SQUARES
text: 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‖22.
Here, f : RN → Rn is a non-linear model with parameters θ ∈ RN . To solve this problem, we run gradient descent with a fixed stepsize η, starting from an initial point θ0, with updates of the form
θτ+1 = θτ − η∇L(θτ ) where ∇L(θ) = J T (θ)(f(θ)− y). (14)
Here, J (θ) ∈ Rn×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
Llin(θ) = 1
2 ‖f(θ0) + J(θ − θ0)− y‖22. (15)
Here, J ∈ Rn×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
θ̃τ+1 = θ̃τ − ηJT ( f(θ0) + J(θ̃τ − θ0)− y ) = θ̃τ − ηJTJ ( θ̃τ − θ0 ) − ηJT (f(θ0)− y) . (16)
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
α ≤ σn ≤ σ1 ≤ β. (17) Furthermore, we assume that the Jacobian mapping associated with the nonlinear model f obeys
‖J (θ)‖ ≤ β for all θ ∈ RN . (18) 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
‖J (θ0)− J‖ ≤ 0. (19) 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
‖J (θ)− J (θ0)‖ ≤ 2 , for all θ ∈ BR(θ0), (20)
where BR(θ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: rτ := f(θτ )− y (21) linear residual: r̃τ := ( I− ηJJT )τ r0, (22)
are close in the proximity of the initialization. Theorem 4 (Closeness of linear and nonlinear least-squares problems). Assume the Jacobian mapping J (θ) ∈ Rn×N associated with the function f(θ) obeys Assumptions 1, 2, and 3 around an initial point θ0 ∈ RN with respect to a reference Jacobian J ∈ Rn×N and with parameters α, β, 0, , and R. Furthermore, assume the radius R is given by
R 2 := ∥∥J†r0∥∥2 + 1α2 ( 0 + ) (1 + 2ηTβ2) ‖r0‖2, (23)
with T a constant obeying 1 ≤ T ≤ 12η 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
‖rτ − r̃τ‖2 ≤ 2 β
α2 ( 0 + )‖r0‖2 (24)∥∥∥θτ − θ̃τ∥∥∥
2 ≤ 1 α2
( 0 + ) ( 1 + 2ητβ2 ) ‖r0‖2. (25)

heading: Moreover, for all iterates τ ≤ T ,
text: ‖θτ − θ0‖2 ≤ R
2 . (26)
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ΣVT ∈ Rn×N =∑n i=1 σiwiv T i be the singular value decomposition of the matrix J. Then the residual r̃τ after τ iterations of gradient descent with updates (16) is
r̃τ = n∑ i=1 ( 1− ησ2i )τ wi 〈wi, r0〉 . (27)
Moreover, using a step size satisfying η ≤ 1 σ21 , the linearized iterates (16) obey
∥∥∥θ̃τ − θ0∥∥∥2 2 = n∑ i=1 ( 〈wi, r0〉 1− (1− ησ2i )τ σi )2 . (28)
In the next section we combine these two general theorems to provide guarantees for denoising using general neural networks.

heading: E.1 PROOF OF THEOREM 4 (CLOSENESS OF LINEAR AND NON-LINEAR LEAST-SQUARES)
text: 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 θτ ∈ BR(θ0). To prove θτ ∈ BR(θ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
≤ ∥∥JT r̃τ−1∥∥2 + ‖J (θτ−1)− J‖ ‖r̃τ−1‖2 + ‖J (θτ−1)‖ ‖rτ−1 − r̃τ−1‖2
(i) ≤ ∥∥JT r̃τ−1∥∥2 + ( + 0)‖r̃τ−1‖2 + β‖rτ−1 − r̃τ−1‖2
(ii) ≤ 1 η ∥∥J†r0∥∥2 + ( + 0)‖r0‖2 + 2β2α2 ( 0 + )‖r0‖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 ≤ ‖r0‖2, and the first term with the bound∥∥JT r̃τ−1∥∥2 = ∥∥JT (I− ηJJT )τ−1r0∥∥2
= ∥∥Σ(I− ηΣ2)τ−1WT r0∥∥2
≤ √√√√ n∑ j=1 σ2j 〈wj , r0〉2
≤ β2 √√√√ n∑
j=1
1
σ2j 〈wj , r0〉2
= β2 ∥∥J†r0∥∥2 ≤ 1 η
∥∥J†r0∥∥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 ∈ BR(θ0). Then, provided the stepsize obeys η ≤ 1/β2, the deviation of the residuals obeys
‖eτ+1‖2 ≤ ηβ ( 0 + ) ‖r̃τ‖2 + ( 1 + η 2 ) ‖eτ‖2. (29)
By the previous step, θτ ,θτ+1 ∈ BR(θ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:
‖r̃τ‖2 ≤ ( 1− ηα2 )τ ‖r0‖2. (30) 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 r̃τ ≤ ( 1− ηα2 )τ ρ (31)
eτ ≤ ( 1 + η 2 ) eτ−1 + ηξr̃τ−1 (32)
Then, for all τ ≤ 12η 2 , we have that
eτ ≤ 2ξ ρ
α2 . (33)
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 θτ ∈ BR(θ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 ρ = ‖r0‖2. Thus, Lemma 2 implies that for all τ ≤ 12η 2
‖eτ‖2 ≤ 2 β
α2 ( 0 + )‖r0‖2. (34)
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
‖J (θτ )− J‖ ≤ ‖J (θτ )− J (θ0)‖ + ‖J (θ0)− J‖ ≤ 0 + . (35)
The difference between the parameter of the original iterate θ and the linearized iterate θ̃ obey
1
η ∥∥∥θτ − θ̃τ∥∥∥ 2 = ∥∥∥∥∥ τ−1∑ t=0 ∇L(θt)−∇Llin(θ̃t) ∥∥∥∥∥ 2
= ∥∥∥∥∥ τ−1∑ t=0 J T (θt)rt − JT r̃t ∥∥∥∥∥ 2
≤ τ−1∑ t=0 ∥∥(J T (θt)− JT )r̃t∥∥2 + ∥∥J T (θt)(rt − r̃t)∥∥2 (i) ≤ τ−1∑ t=0 ( 0 + )‖r̃t‖2 + β‖et‖2
(ii) ≤ ( 0 + ) τ−1∑ t=0 ( 1− ηα2 )τ−1 ‖r0‖2 + β τ−1∑ t=0 ‖et‖2
= ( 0 + ) 1− (1− ηα2)τ
ηα2 ‖r0‖2 + β τ−1∑ t=0 ‖et‖2
(iii) ≤ ( 0 + )
ηα2 ‖r0‖2 + 2τ
β2 α2 ( 0 + )‖r0‖2
= 1
ηα2 ( 0 + )
( 1 + 2ητβ2 ) ‖r0‖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†r0∥∥2 + 1α2 ( 0 + ) (1 + 2ητβ2) ‖r0‖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 〈wi, r0〉2 (1− (1− ησ2i )τ )2 σ2i
≤ n∑ i=1 〈wi, r0〉2 /σ2i
= ∥∥J†r0∥∥2.
Finally, inequality (ii) follows from the definition of R in equation (23) along with τ ≤ T , by assumption.

heading: E.1.1 PROOF OF LEMMA 2
text: 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τ ≤ η 2eτ + ηξr̃τ ≤ η 2eτ + ηξ(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 (et+1 − et) ξ
≤ τη 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 τ ≤ 12η 2 .

heading: E.2 PROOF OF THEOREM 5
text: The proof of identity (27) is equivalent to the proof of Proposition (1). Regarding inequality (28), note that
θ̃τ − θ̃0 = −η τ−1∑ t=0 ∇Llin(θ̃t) = −η τ−1∑ t=0 JT r̃t = −ηV ( τ−1∑ t=0 Σ ( I− ηΣ2 )t)〈w1, r0〉... 〈wn, r0〉  . Thus〈
vi, θ̃τ − θ̃0 〉 = −ησi 〈wi, r0〉 ( τ−1∑ t=0 (1− ησ2i )t ) = −ησi 〈wi, r0〉 1− (1− ησ2i )τ ησ2i . (36)
This in turn implies that∥∥∥θ̃τ − θ̃0∥∥∥2 2 = n∑ i=1 〈 vi, θ̃τ − θ̃0 〉2 = n∑ i=1 ( 〈wi, r0〉 1− (1− ησ2i )τ σi )2 .

heading: F PROOFS FOR NEURAL NETWORK GENERATORS (PROOF OF THEOREM 3)
text: 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 ∈ Rn×nk, that is very close to the original Jacobian at initialization, J (C0), and that obeys JJT =
E [ J (C)J T (C) ] . The concrete choice of J is discussed in the next paragraph. We apply Theorem 4 with the following choices of parameters:
α = σn (Σ(U)) , β = ‖U‖ , 0 = β ( 4 log( 2nδ )
k
)1/4 , = ξ
8 β α2 β2 , ω = ‖y‖2√ nβ ξ α2 β2 .
A few comments on the intuition behind those choices: As shown below, α and β are chosen so that Assumption 1 is satisfied. Regarding the choices of , 0, and ω, recall that the difference between the linear and non-linear residual is bounded by 2 βα2 ( 0 + )‖r0‖2 (by Theorem 4). We want this error to be bounded by ξ‖y‖2. This bound is guaranteed by our choice of , by choosing 0 such that it obeys 0 ≤ , and by choosing ω so that the initial residual is bounded by ‖r0‖2 ≤ 2‖y‖2. Note that 0 ≤ holds by assumption (10). We next verify by applying a series of Lemmas proven in Appendix I that the conditions of Theorem 4 are satisfied (specifically, Assumptions 1, 2, 3 on the Jacobians of the non-linear map and the associated linearized map, to be bounded and sufficiently close to each other).
Choice of reference Jacobian and verification of assumption 2: We chose the reference Jacobian J so that it is 0 close to the random initialization J (C0) (with high probability), and thus satisfies assumption 2. Towards this goal, we require the following two lemmas.
Lemma 3 (Concentration lemma). Consider G(C) = ReLU(UC)v with v ∈ Rk and U ∈ Rn×k and associated Jacobian J (C) (9). Let C ∈ Rn×k be generated at random with i.i.d. N (0, ω2) entries. Then,
∥∥J (C)J T (C)−Σ (U)∥∥ ≤ ‖U‖2 √√√√log(2n
δ ) k∑ `=1 v4` ,
holds with probability at least 1− δ.
To see that Lemma 3 implies the condition (19), we use the following lemma.
Lemma 4 ((Oymak et al., 2019, Lem. 6.4)). Let X ∈ Rn×N , N ≥ n and let Σ be n×n psd matrix obeying ∥∥XXT −Σ∥∥ ≤ 20/4, for a scalar ≥ 0. Then there exists a matrix J ∈ Rn×N obeying Σ = JJT such that
‖J−X‖ ≤ 0.
In order to verify the condition (19), note that using the fact that ∑k ` v 4 ` = 1 k by Lemma 3, with probability at least 1− δ, we have
∥∥J (C0)J T (C0)−Σ(U)∥∥ ≤ ‖U‖2 √ log ( 2n δ ) k = ( 0/2) 2, (37)
Combining inequality (37) with Lemma 4 (applied with X = J (C0)), it follows that condition (19) holds for the chosen value of 0, concluding the proof of Assumption 2 being satisfied. This proof specifies the reference Jacobian J as the matrix that is 0-close to J (C0), and that exists by Lemma 4.
Verifying Assumption 1: To verify Assumption 1, note that by definition JJT = Σ(U) thus trivially σn (Σ(U)) ≥ α holds. Furthermore, Lemma 5 below combined with the fact that ‖v‖2 = 1 implies that ‖J‖ ≤ β and ‖J (C)‖ ≤ β for all C. This completes the verification of Assumption 1. It remains to show that the Jacobian has bounded spectrum:
Lemma 5 (Spectral norm of Jacobian). Consider G(C) = ReLU(UC)v with v ∈ Rk and U ∈ Rn×k and associated Jacobian J (C) (9), and let J be any matrix obeying JJT = E [ J (C)J T (C) ] , where the expectation is over a matrix C with iid N (0, ω2) entries. Then
‖J (C)‖ ≤ ‖v‖2 ‖U‖ and ‖J‖ ≤ ‖v‖2 ‖U‖ .
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 ∈ Rn×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
‖r0‖2 ≤ ‖y‖2 + ‖G(C0)‖2 ≤ 3
2 ‖y‖2, (38)
where we used that, by Lemma 6,
‖G(C0)‖2 ≤ ω √ 8 log(2n/δ) ‖U‖F ≤ ‖y‖2ξ α2
β2
√ 8 log(2n/δ) ≤ 1
2 ‖y‖2. (39)
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 C0 be a matrix with i.i.d. N(0, ω2) entries. Then, for all C obeying
‖C−C0‖ ≤ ωR̃ with R̃ ≤ 1
2
√ k,
the Jacobian mapping (9)associated with the generator G(C) = ReLU(UC)v obeys
‖J (C)− J (C0)‖ ≤ ‖v‖∞ 2(kR̃) 1/3 ‖U‖ ,
with probability at least 1− ne− 12 R̃4/3k7/3 .
In order to verify Assumption 3, first note that the radius in the theorem, defined in equation (23), obeys
R = 2 ∥∥J†r0∥∥2 + 2α2 ( 0 + ) (1 + 2ηTβ2) ‖r0‖2
(i) ≤ ( 2
α +
2
α2 ( 0 + )
( 1 + 2ηTβ2 )) ‖r0‖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†Sr0∥∥∥
2 ≤ 1α‖r0‖2, (ii) from the bound on the initial resid-
ual (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 (C0)‖ ≤ ‖v‖∞ 2(kR̃) 1/3 ‖U‖ = 1√ k 2 ξ 16
α2
β2 (k · k1/2)1/3β = ξ 8
α2
β = .
holds with probability at least
1− ne− 12 R̃k 7/3 = 1− ne−2 −17ξ4 α
8 β8 k3 (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 ≤ 12η 2 , from Theorem 4 is satisfied under the number of constraints of Theorem 3 by 12η 2 = 25β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)τ WT r0 + z∥∥∥
2 +
1 2 ξ‖y‖2
(iv) = ∥∥∥W (I− ηΣ2)τ WT (G(C0)− x)− (W (I− ηΣ2)τ WT − I) z∥∥∥
2 +
1 2 ξ‖r0‖2
(v) ≤ ∥∥∥(I− ηΣ2)τ WTx∥∥∥ 2 + ∥∥∥((I− ηΣ2)τ − I)WT z∥∥∥ 2 + ‖G(C0)‖2 + 1 2 ξ‖y‖2
(vi) ≤ (1− ησ2p)τ‖x‖2 + √√√√ n∑ i=1 ((1− ησ2i )τ − 1)2 〈wi, z〉 2 + ξ‖y‖2.
Here, (i) follows from the triangular inequality, (ii) from Theorem 4 equation (24) combined with our choice for and ‖r0‖2 ≤ 1 2‖y‖2, shown above. Moreover, (iii) follows from Theorem 5, (iv) from r0 = G(C0)−y = G(C0)−x−z, (v) from the triangular inequality. Finally, for (vi) we used the fact that x ∈ span{w1,w2, . . . ,wp} combined with the fact that
∥∥W (I− ηΣ2)τ WT∥∥ ≤ 1, as well as using the bound ‖G(C0)‖2 ≤ 1 2‖y‖2 proven in (39). This proves the final bound (12), as desired.

heading: G PROOF OF PROPOSITION 1 AND EQUATION (2)
text: The residual of gradient descent at iteration τ is rτ = y − Jcτ
= y − J(cτ−1 − ηJT (Jcτ−1 − y)) = (I− ηJJT )(y − Jcτ−1) = (I− ηJJT )τ (y − Jc0) = (I− ηJJT )τy = (I− ηWΣ2WT )τy
where we used that c0 = 0 and the SVD J = WΣVT . Expanding y in terms of the singular vectors wi (i.e., the columns of W), as y = ∑ i wi 〈wi,y〉, and noting that (I − ηWΣ2WT )τ =∑
i(1− ησ2i )τwiwTi we get rτ = ∑ i (1− ησ2i )τwi 〈wi,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 wi(1− ησ2i )τ 〈wi,x〉+ n∑ i=1 wi((1− ησ2i )τ − 1) 〈wi, z〉 .
By the triangle inequality,
‖x− Jcτ‖2 ≤ ∥∥∥∥∥ p∑ i=1 wi(1− ησ2i )τ 〈wi,x〉 ∥∥∥∥∥ 2 + ∥∥∥∥∥ n∑ i=1 wi((1− ησ2i )τ − 1) 〈wi, z〉 ∥∥∥∥∥ 2
≤ (1− ησ2p)τ‖x‖2 + √√√√ n∑ i=1 ((1− ησ2i )τ − 1)2 〈wi, z〉 2 ,
where the second inequality follows by using orthogonality of the wi and by using (1− ησ2i ) ≤ 1, from η ≤ 1/σ2max. This concludes the proof of equation (2).

heading: H THE EXPECTED JACOBIAN FOR CONVOLUTIONAL GENERATORS
text: 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
J (C)J T (C) = k∑ `=1 v2`diag(σ ′(Uc`))UU Tdiag(σ′(Uc`))
= k∑ `=1 v2`σ ′(Uc`)σ ′(Uc`) T UUT
= σ′(UC)diag(v21 , . . . , v 2 k)σ ′(UC) T UUT , where denotes the entrywise product of the two matrices. Then,
E [ J (C)J T (C) ] = k∑ `=1 v2`E [ σ′(Uc`)σ ′(Uc`) T ] UUT . (40)
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 ( 〈ui,uj〉 ‖ui‖2‖uj‖2 ) /π ) .
Using that ‖v‖2 = 1, we get[ E [ J (C)J T (C) ]] ij = 1
2
( 1− cos−1 ( 〈ui,uj〉 ‖ui‖2‖uj‖2 ) /π ) 〈ui,uj〉 ,
where ui 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 ∈ Rn×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−1ÛF,
where F ∈ Cn×n is the DFT matrix with entries [F]jk = e
i2πjk/n, j, k = 0, . . . , n− 1, and Û 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 û = Fu is conjugate symmetric (since u is real) so that
[û]n−k+2 = û ∗, k = 2, . . . , n/2.
Let U = WΣVT 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 û.

heading: I PROOFS OF LEMMAS FOR NEURAL NETWORK DENOISERS (PROOFS OF AUXILIARY LEMMAS IN SECTION F)

heading: I.1 PROOF OF LEMMA 7: JACOBIAN PERTURBATION AROUND INITIALIZATION
text: 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
‖J (C)− J (C′)‖ ≤ ‖v‖∞‖U‖maxj ∥∥σ′(uTj C)− σ′(uTj C′)∥∥2. (41)
To see this, first note that, by (9),
J (C)− J (C′) = [. . . vj(diag(σ′(Ucj))− diag(σ′(Ucj))′)U . . .].
This in turn implies that ‖J (C)− J (C′)‖2 = ∥∥∥(J (C)− J (C′))(J (C)− J (C′))T∥∥∥
= ∥∥∥(σ′(UC)− σ′(UC′))diag(v21 , . . . , v2k)(σ′(UC)− σ′(UC′))T UUT∥∥∥2
(i) ≤ ‖U‖2 max
j ∥∥(σ′(uTj C)− σ′(uTj C′))diag(v)∥∥2 ≤ ‖v‖2∞‖U‖
2 max j ∥∥σ′(uTj C)− σ′(uTj C′)∥∥2, where for (i) we used that for two positive semidefinite matrices A,B, λmax(A B) ≤ λmax(A) maxi Bii. This concludes the proof of equation (41).
Step 2: Step one implies that we only need to control σ′(Ucj) 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 |uTi C′|π(q) ‖ui‖ .
Then max i ∥∥σ′(uTi C)− σ′(uTi C′)∥∥ ≤√2q Proof. Suppose that σ′(uTi C) and σ ′(uTi 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 ≥ ∥∥uTi (C−C′)∥∥2/‖ui‖2
≥ q |uTi C′|2π(q) ‖ui‖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−kq2/2, the q-th smallest entry of uTi C ′ ∈ Rk obeys
|uTi C′|π(q) ‖ui‖ ≥ q 2k ν for all i = 1, . . . , n. (42)
We note that it is sufficient to prove this result for ν = 1. This follows from anti-concentration of Gaussian random variables. Specifically, with the entries of C′ being iid N (0, 1) distributed, the entries of g = uTi C
′/‖ui‖ ∈ Rk are iid standard Gaussian random variables as well. We show that with probability at least 1 − e−kq2/2, at most q entries are larger than q2k . Let γδ be the number for which P [|g`| ≤ γδ] ≤ δ, where g is a standard Gaussian random variable. Note that γδ ≥ √ π/2δ ≥ δ. Define the random variable
δ` = { 1 if |g`| ≤ γδ, 0 otherwise.
with δ = q2k . With E [δ`] = δ, by Hoeffding’s inequality,
P [ k∑ `=1 δ` ≥ m ] = P [ k∑ `=1 δ` − E [δ`] ≥ m/2 ] ≤ e−2k(m/2) 2 = e−km 2/2. (43)
Thus, with probability at least 1− ke−kq2/2 no more than m entries are smaller than γδ ≥ δ = q2k . The results now follows by taking the union bound over all i = 1, . . . , n.
We are now ready to conclude the proof of the lemma. By equation (41), ‖J (C)− J (C′)‖ ≤ ‖v‖∞‖U‖maxj ∥∥σ′(uTj C)− σ′(uTj C′)∥∥
≤ ‖v‖∞‖U‖ √ 2q
provided that ‖C−C′‖ ≤ √q q
2k ν,
with probability at least 1 − ne−kq2/2. Setting q = (2kR)2/3 concludes the proof (note that the assumption R ≤ 12 √ k ensures q ≤ k).

heading: I.2 PROOF OF LEMMA 5: BOUNDED JACOBIAN
text: By the expression of the Jacobian in equation (9), ‖J (C)‖2 = ∥∥∥J (C)J (C)T∥∥∥
= ∥∥∥σ′(UC)diag(v21 , . . . , v2k)σ′(UC)T UUT∥∥∥2
(i) ≤ ‖U‖2 max
j ∥∥σ′(uTj C)diag(v)∥∥22 ≤ ‖v‖22 ‖U‖ 2 ,
where for (i) we used that for two positive semidefinite matrices A,B, λmax(A B) ≤ λmax(A) maxi Bii. To prove the second inequality note that∥∥JJT∥∥ = ‖Σ(U)‖
=‖v‖22 ∥∥∥Ec∼N (0,I) [ReLU′(Uc)ReLU′(Uc)T ] (UUT )∥∥∥
(i) ≤‖v‖22 ∥∥UUT∥∥ . Here, (i) follows from the fact that for two positive semidefinite matrices A,B, λmax(A B) ≤ λmax(A) maxi Bii.

heading: I.3 PROOF OF LEMMA 3: CONCENTRATION LEMMA
text: We begin by defining the zero-mean random matrices
S` = v 2 ` ( σ′(Ũc`)σ ′(Ũc`) T − E [ σ′(Ũc`)σ ′(Ũc`) T ]) ( ŨŨT ) .
With this notation,
J (C)J T (C)−Σ(Ũ) = k∑ `=1 v2` ( σ′(Ũc`)σ ′(Ũc`) T − E [ σ′(Ũc`)σ ′(Ũc`) 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( σ′(Ũc`)σ ′(Ũc`) T − E [ σ′(Ũc`)σ ′(Ũc`) T ]) ( ŨŨT ) ( σ′(Ũc`)σ ′(Ũc`) T ) ( ŨŨT ) =diag ( σ′(Ũc`) T ) ŨŨT diag ( σ′(Ũc`)
) B2ŨŨT
Similarly,( σ′(Ũc`)σ ′(Ũc`) T − E [ σ′(Ũc`)σ ′(Ũc`) T ]) ( ŨŨT ) − E [ σ′(Ũc`)σ ′(Ũc`) T ] ( ŨŨT ) −B2ŨŨT .
Thus,
−v2`B2ŨŨT S` v2`B2ŨŨT .
Therefore, using A` := v2`B 2ŨŨT we have
S2` A2` ,
and
σ2 := ∥∥∥∥∥ k∑ `=1 A2` ∥∥∥∥∥ ≤ B4 ( k∑ `=1 v4` )∥∥∥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 S2` A2` almost surely.

heading: Then, for all t ≥ 0,
text: P [∥∥∥∥∥ k∑ `=1 S` ∥∥∥∥∥ ≥ t ] ≤ 2ne− t2 8σ2 where σ2 := ∥∥∥∥∥ k∑ `=1 A2` ∥∥∥∥∥ . Therefore, applying matrix Hoeffding inequality we get that
P [∥∥∥∥∥ k∑ `=1 S` ∥∥∥∥∥ ≥ t ] ≤ 2ne − t2 B4‖v‖44‖Ũ‖ 4 ,
which concludes the proof.

heading: I.4 PROOF OF LEMMA 6 (BOUND ON INITIAL RESIDUAL)
text: Without loss of generality we prove the result for ν = 1. First note that by the triangle inequality
‖r0‖2 = ‖σ(UC)v − y‖2 ≤ ‖σ(UC)v‖2 + ‖y‖2.
We next bound ‖σ(UC)v‖2. Consider the i-th entry of the vector σ(UC)v ∈ Rn, given by σ(uTi C)v, and note that qj = (σ(u T i cj) − σ(uTi cn−j))/‖ui‖2 is sub-Gaussian with parameter 2, i.e., P [|qj | ≥ t] ≤ 2e−2t 2 . It follows that
P ∣∣∣∣∣∣ k/2∑ j=1 qj ∣∣∣∣∣∣ ≥ β√k  ≤ 2e− β28 .
Thus, P [∣∣σ(uTi C)v∣∣ ≥ ‖ui‖2ξβ] ≤ 2e− β28 , where we used that |vj | = ξ/ √ k. Taking a union bound over all n entries,
P [ ‖σ(UC)v‖22 ≥ ‖U‖ 2 F ξ 2β2 ] ≤ 2ne− β2 8 .
Choosing β = √ 8 log(2n/δ) concludes the proof.

Review summaries:
This paper studies the situation in which a two-layer CNN with RELU nonlinearity is fit to a single image and the observation that it is able to fit a "natural" image in fewer iterations than a "noisy" image. Theorems on the convergence of this fitting are discussed and proven in the appendices. Intermediate results study the convergence of fitting a linear model to an image plus noise and fitting a single-layer CNN. Denoising is demonstrated on two images with additive white Gaussian noise and the approach under study is shown to provide a better signal-to-noise ratio than BM3D, another untrained denoising approach. The main result is that the use of upsampling via a fixed interpolation filter provides an inductive bias towards "natural" images.