Convolutional Neural Networks Analyzed via Convolutional Sparse Coding

Introduction

Deep learning (LeCun et al., 2015), and in particular CNN (LeCun et al., 1990, 1998; Krizhevsky et al., 2012), has gained a copious amount of attention in recent years as it has led to many state-of-the-art results spanning through many fields – including speech recognition (Bengio et al., 2003; Hinton et al., 2012; Mikolov et al., 2013), computer vision (Farabet et al., 2013; Simonyan and Zisserman, 2014; He et al., 2015), signal and image processing (Gatys et al., 2015; Ulyanov et al., 2016; Johnson et al., 2016; Dong et al., 2016), to name a few. In the context of CNN, the forward pass is a multi-layer scheme that provides an end-to-end mapping, from an input signal to some desired output. Each layer of this algorithm consists of three steps. The first convolves the input with a set of learned filters, resulting in a set of feature (or kernel) maps. These then undergo a point wise non-linear function, in a second step, often resulting in a sparse outcome (Glorot et al., 2011). A third (and optional) down-sampling step, termed pooling, is then applied on the result in order to reduce its dimensions. The output of this layer is then fed into another one, thus forming the multi-layer structure, often termed forward pass.

Despite its marvelous empirical success, a clear and profound theoretical understanding of this scheme is still lacking. A few preliminary theoretical results were recently suggested. In (Mallat, 2012; Bruna and Mallat, 2013) the Scattering Transform was proposed, suggesting to replace the learned filters in the CNN with predefined Wavelet functions. Interestingly, the features obtained from this network were shown to be invariant to various transformations such as translations and rotations. Other works have studied the properties of deep and fully connected networks under the assumption of independent identically distributed random weights (Giryes et al., 2015; Saxe et al., 2013; Arora et al., 2014; Dauphin et al., 2014; Choromanska et al., 2015). In particular, in (Giryes et al., 2015) deep neural networks were proven to preserve the metric structure of the input data as it propagates through the layers of the network. This, in turn, was shown to allow a stable recovery of the data from the features obtained from the network.

Another prominent paradigm in data processing is the sparse representation concept, being one of the most popular choices for a prior in the signal and image processing communities, and leading to exceptional results in various applications (Elad and Aharon, 2006; Dong et al., 2011; Zhang and Li, 2010; Jiang et al., 2011; Mairal et al., 2014). In this framework, one assumes that a signal can be represented as a linear combination of a few columns (called atoms) from a matrix termed a dictionary. Put differently, the signal is equal to a multiplication of a dictionary by a sparse vector. The task of retrieving the sparsest representation of a signal over a dictionary is called sparse coding or pursuit. Over the years, various algorithms were proposed to tackle this problem, among of which we mention the thresholding algorithm (Elad, 2010) and its iterative variant (Daubechies et al., 2004). When handling natural signals, this model has been commonly used for modeling local patches extracted from the global data mainly due to the computational difficulties related to the task of learning the dictionary (Elad and Aharon, 2006; Dong et al., 2011; Mairal et al., 2014; Romano and Elad, 2015; Sulam and Elad, 2015). However, in recent years an alternative to this patch-based processing has emerged in the form of the Convolutional Sparse Coding (CSC) model (Bristow et al., 2013; Kong and Fowlkes, 2014; Wohlberg, 2014; Gu et al., 2015; Heide et al., 2015; Papyan et al., 2016a, b). This circumvents the aforementioned limitations by imposing a special structure – a union of banded and Circulant matrices – on the dictionary involved. The traditional sparse model has been extensively studied over the past two decades (Elad, 2010; Foucart and Rauhut, 2013). More recently, the convolutional extension was extensively analyzed in (Papyan et al., 2016a, b), shedding light on its theoretical aspects and prospects of success.

In this work, by leveraging the recent study of CSC, we aim to provide a new perspective on CNN, leading to a clear and profound theoretical understanding of this scheme, along with new insights. Embarking from the classic CSC, our approach builds upon the observation that similar to the original signal, the representation vector itself also admits a convolutional sparse representation. As such, it can be modeled as a superposition of atoms, taken from a different convolutional dictionary. This rationale can be extended to several layers, leading to the definition of our proposed ML-CSC model. Building on the recent analysis of the CSC, we provide a theoretical study of this novel model and its associated pursuits, namely the layered thresholding algorithm and the layered basis pursuit (BP).

Our analysis reveals the relation between the CNN and the ML-CSC model, showing that the forward pass of the CNN is in fact identical to our proposed pursuit – the layered thresholding algorithm. This connection is of significant importance since it gives a clear mathematical meaning, objective and model to the CNN architecture, which in turn can be accompanied by guarantees for the success of the forward pass, studied via the layered thresholding algorithm. Specifically, we show that the forward pass is guaranteed to recover an estimate of the underlying representations of an input signal, assuming these are sparse in a local sense. Moreover, considering a setting where a norm-bounded noise is added to the signal, we show that such a mild corruption in the input results in a bounded perturbation in the output – indicating the stability of the CNN in recovering the underlying representations. Lastly, we exploit the answers to the above questions in order to propose an alternative to the commonly used forward pass algorithm, which is tightly connected to both deconvolutional (Zeiler et al., 2010; Pu et al., 2016) and recurrent networks (Bengio et al., 1994), and also related to residual networks (He et al., 2015). The proposed alternative scheme is accompanied by a thorough theoretical study. Although this and the analysis presented throughout this work focus on CNN, we will show that they also hold for fully connected networks.

This paper is organized as follows. In Section 2 we review the basics of both the CNN and the Sparse-Land model. We then define the proposed ML-CSC model in Section 3, together with its corresponding deep sparse coding problem. In Section 4, we aim to solve this using the layered thresholding algorithm, which is shown to be equivalent to the forward pass of the CNN. Next, having established the relevance of our model to CNN, we proceed to its analysis in Section 5. Standing on these theoretical grounds, we then propose in Section 6 a provably improved pursuit, termed the layered BP, accompanied by its theoretical analysis. We revisit the assumptions of our model in Section 7. First, in Section 7.1 we link the double sparsity model to ours by assuming the dictionaries throughout the layers are sparse. Then, in Section 7.2 we consider an idea typically employed in CNN, termed spatial-stride, showing its benefits from a simple theoretical perspective. Combining our insights from Section 7.1 and 7.2, we move to an experimental phase by constructing a family of signals satisfying the assumptions of our model, which are then used in order to verify our theoretical results. Finally, in Section 9 we conclude the contributions of this paper and present several future directions.

Background

This section is divided into two parts: The first is dedicated to providing a simple mathematical formulation of the CNN and the forward pass, while the second reviews the Sparse-Land model and its various extensions. Readers familiar with these two topics can skip directly to Section 3, which moves to serve the main contribution of this work.

By changing the order of the columns in the convolutional matrix, one can observe that it can be equally viewed as a concatenation of banded and CirculantWe shall assume throughout this paper that boundaries are treated by a periodic continuation, which gives rise to the cyclic structure. matrices, as depicted in Figure 1(b). Using this observation, the above description for one dimensional signals can be extended to images, with the exception that now every Circulant matrix is replaced by a block Circulant with Circulant blocks one.

An illustration of the forward pass algorithm is presented in Figure 2(a) and 2(b). In Figure 2(a) one can observe that ${\mathbf{W}}_{2}$ is not a regular convolutional matrix but a stride one, since it shifts local filters by skipping $m_{1}$ entries at a time. The reason for this becomes apparent once we look at Figure 2(b); the convolutions of the second layer are computed by shifting the filters of ${\mathbf{W}}_{2}$ that are of size $\sqrt{n_{1}}\times\sqrt{n_{1}}\times m_{1}$ across $N$ places, skipping $m_{1}$ indices at a time from the $\sqrt{N}\times\sqrt{N}\times m_{1}$-sized array. A matrix obeying this structure is called a stride convolutional matrix.

Thus far, we have presented the basic structure of CNN. However, oftentimes an additional non-linear function, termed pooling, is employed on the resulting feature map obtained from the ReLU operator. In essence, this step summarizes each $w_{i}$-dimensional spatial neighborhood from the $i$-th kernel map ${\mathbf{Z}}_{i}$ by replacing it with a single value. If the neighborhoods are non-overlapping, for example, this results in the down-sampling of the feature map by a factor of $w_{i}$. The most widely used variant of the above is the max pooling (Krizhevsky et al., 2012; Simonyan and Zisserman, 2014), which picks the maximal value of each neighborhood. In (Springenberg et al., 2014) it was shown that this operator can be replaced by a convolutional layer with increased stride without loss in performance in several image classification tasks. Moreover, the current state-of-the-art in image recognition is obtained by the residual network (He et al., 2015), which does not employ any pooling steps (except for a single layer). As such, we defer the analysis of this operator to a follow-up work.

In the context of classification, for example, the output of the last layer is fed into a simple classifier that attempts to predict the label of the input signal ${\mathbf{X}}$, denoted by $h({\mathbf{X}})$. Given a set of signals $\{{\mathbf{X}}_{j}\}_{j}$, the task of learning the parameters of the CNN – including the filters $\{{\mathbf{W}}_{i}\}_{i=1}^{K}$, the biases $\{{\mathbf{b}}_{i}\}_{i=1}^{K}$ and the parameters of the classifier ${\mathbf{U}}$ – can be formulated as the following minimization problem

In the remainder of this work we shall focus on the feature extraction process and assume that the parameters of the CNN model are pre-trained and fixed. These, for example, could have been obtained by minimizing the above objective via the backpropagation algorithm and the stochastic gradient descent, as in the VGG network (Simonyan and Zisserman, 2014).

2 Sparse-Land

This section presents an overview of the Sparse-Land model and its many extensions. We start with the traditional sparse representation and the core problem it aims to solve, and then proceed to its nonnegative variant. Next, we continue to the dictionary learning task both in the unsupervised and supervised cases. Finally, we describe the recent CSC model, which will lead us in the next section to the proposal of the ML-CSC model. This, in turn, will naturally connect the realm of sparsity to that of the CNN.

For a fixed dictionary, given a signal ${\mathbf{X}}$, the task of recovering its sparsest representation ${\bm{\Gamma}}$ is called sparse coding, or simply pursuit, and it attempts to solve the following problem (Donoho and Elad, 2003; Tropp, 2004; Elad, 2010):

where we have denoted by $\|{\bm{\Gamma}}\|_{0}$ the number of non-zeros in ${\bm{\Gamma}}$. The above has a convex relaxation in the form of the Basis-Pursuit (BP) problem (Chen et al., 2001; Donoho and Elad, 2003; Tropp, 2006), formally defined as

Tighter conditions, relying on sharper characterizations of the dictionary, were also suggested in the literature (Candes et al., 2006; Schnass and Vandergheynst, 2007; Candes et al., 2006; Candes and Tao, 2007). However, at this point, we shall not dwell on these.

One of the simplest approaches for tackling the ${\text{P}_{0}}$ and ${\text{P}_{1}}$ problems is via the hard and soft thresholding algorithms, respectively. These operate by computing the inner products between the signal ${\mathbf{X}}$ and all the atoms in ${\mathbf{D}}$ and then choosing the atoms corresponding to the highest responses. This can be described as solving, for some scalar $\beta$, the following problems:

for the ${\text{P}_{1}}$. The above are simple projection problems that admit a closed-form solution in the formThe curious reader may identify the relation between the notations used here and the ones in the previous subsection, which starts to reveal the relation between CNN and sparsity-inspired models. This connection will be made stringer and clearer as we proceed to CSC. of ${\mathcal{H}}_{\beta}({\mathbf{D}}^{T}{\mathbf{X}})$ or ${\mathcal{S}}_{\beta}({\mathbf{D}}^{T}{\mathbf{X}})$, where we have defined the hard thresholding operator ${\mathcal{H}}_{\beta}(\cdot)$ by

and the soft thresholding operator ${\mathcal{S}}_{\beta}(\cdot)$ by

Both of the above, depicted in Figure 3, nullify small entries and thus promote a sparse solution. However, while the hard thresholding operator does not modify large coefficients (in absolute value), the soft thresholding does, by contracting these to zero. This inherent limitation of the soft version will appear later on in our theoretical analysis.

As for the theoretical guarantees for the success of the simple thresholding algorithms; these depend on the properties of ${\mathbf{D}}$ and on the ratio between the minimal and maximal coefficients in absolute value in ${\bm{\Gamma}}$, and thus are weaker when compared to those found for OMP and BP (Donoho and Elad, 2003; Tropp, 2004; Donoho et al., 2006). Still, under some conditions, both algorithms are guaranteed to find the true support of ${\bm{\Gamma}}$ along with an approximation of its true coefficients. Moreover, a better estimation of these can be obtained by projecting the input signal onto the atoms corresponding to the found support (indices of the non-zero entries) by solving a Least-Squares problem. This step, termed debiasing (Elad, 2010), results in a more accurate identification of the non-zero values.

2.2 Nonnegative Sparse Coding

The nonnegative sparse representation model assumes a signal can be decomposed into a multiplication of a dictionary and a nonnegative sparse vector. A natural question arising from this is whether such a modification to the original Sparse-Land model affects its expressiveness. To address this, we hereby provide a simple reduction from the original sparse representation to the nonnegative one.

Consider a signal ${\mathbf{X}}={\mathbf{D}}{\bm{\Gamma}}$, where the signs of the entries in ${\bm{\Gamma}}$ are unrestricted. Notice that this can be equally written as

where we have split the vector ${\bm{\Gamma}}$ to its positive coefficients, ${\bm{\Gamma}}_{P}$, and its negative ones, ${\bm{\Gamma}}_{N}$. Since the coefficients in ${\bm{\Gamma}}_{P}$ and $-{\bm{\Gamma}}_{N}$ are all positive, one can thus assume the signal ${\mathbf{X}}$ admits a non-negative sparse representation over the dictionary $\left[{\mathbf{D}},-{\mathbf{D}}\right]$ with the vector $\left[{\bm{\Gamma}}_{P},-{\bm{\Gamma}}_{N}\right]^{T}$. Thus, restricting the coefficients in the sparsity inspired model to be nonnegative does not change its expressiveness.

Similar to the original model, in the nonnegative case, one could solve the associated pursuit problem by employing a soft thresholding algorithm. However, in this case a constraint must be added to the optimization problem in Equation (7), forcing the outcome to be positive, i.e.,

In other words, the ReLU and the soft nonnegative thresholding operator are equal, a fact that will prove to be important later in our work. We should note that a similar conclusion was reached in (Fawzi et al., 2015). To summarize this discussion, we depict in Figure 3 the hard, soft, and nonnegative soft thresholding operators.

2.3 Unsupervised and Task Driven Dictionary Learning

At first, the dictionaries employed in conjunction with the sparsity inspired model were analytically defined matrices, such as the Wavelet and the Fourier (Daubechies et al., 1992; Mallat and Zhang, 1993; Elad and Bruckstein, 2002; Mallat, 2008). Although the sparse coding problem under these can be done very efficiently, over the years many have shifted to a data driven approach – adapting the dictionary ${\mathbf{D}}$ to a set of training signals at hand via some learning procedure. This was empirically shown to lead to sparser representations and better overall performance, at the cost of complicating the involved pursuit, since the dictionary was usually chosen to be redundant (having more columns than rows).

The task of learning a dictionary for representing a set of signals $\{{\mathbf{X}}_{j}\}_{j}$ can be formulated as follows

The above formulation is an unsupervised learning procedure, and it was later extended to a supervised setting. In this context, given a set of signals $\{{\mathbf{X}}_{j}\}_{j}$, one attempts to predict their corresponding labels $\{h({\mathbf{X}}_{j})\}_{j}$. A common approach for tackling this is first solving a pursuit problem for each signal ${\mathbf{X}}_{j}$ over a dictionary ${\mathbf{D}}$, resulting in

and then feeding these sparse representations into a simple classifier, defined by the parameters ${\mathbf{U}}$. The task of learning jointly the dictionary ${\mathbf{D}}$ and the classifier ${\mathbf{U}}$ was addressed in (Mairal et al., 2012), where the following optimization problem was proposed

Double sparsity – first proposed in (Rubinstein et al., 2010) and later employed in (Sulam et al., 2016) – attempts to benefit from both the computational efficiency of analytically defined matrices, and the adaptability of data driven dictionaries. In this model, one assumes the dictionary ${\mathbf{D}}$ can be factorized into a multiplication of two matrices, ${\mathbf{D}}_{1}$ and ${\mathbf{D}}_{2}$, where ${\mathbf{D}}_{1}$ is an analytic dictionary with fast implementation, and ${\mathbf{D}}_{2}$ is a trained sparse one. As a result, the signal ${\mathbf{X}}$ can be represented as

We propose a different interpretation for the above, which is unrelated to practical aspects. Since both the matrix ${\mathbf{D}}_{2}$ and the vector ${\bm{\Gamma}}_{2}$ are sparse, one would expect their multiplication ${\bm{\Gamma}}_{1}={\mathbf{D}}_{2}{\bm{\Gamma}}_{2}$ to be sparse as well. As such, the double sparsity model implicitly assumes that the signal ${\mathbf{X}}$ can be decomposed into a multiplication of a dictionary ${\mathbf{D}}_{1}$ and sparse vector ${\bm{\Gamma}}_{1}$, which in turn can also be decomposed similarly via ${\bm{\Gamma}}_{1}={\mathbf{D}}_{2}{\bm{\Gamma}}_{2}$.

2.4 Convolutional Sparse Coding Model

In the convolutional model, the classical theoretical guarantees (we are referring to results reported in (Chen et al., 2001; Donoho and Elad, 2003; Tropp, 2006)) for the ${\text{P}_{0}}$ problem, defined in Equation (3), are very pessimistic. In particular, the condition for the uniqueness of the underlying solution and the requirement for the success of the sparse coding algorithms depend on the global number of non-zeros being less than $\frac{1}{2}\left(1+\frac{1}{\mu({\mathbf{D}})}\right)$. Following the Welch bound (Welch, 1974), this expression was shown in (Papyan et al., 2016a) to be impractical, allowing the global number of non-zeros in ${\bm{\Gamma}}$ to be extremely low.

In order to provide a better theoretical understanding of this model, which exploits the inherent structure of the convolutional dictionary, a recent work (Papyan et al., 2016a) suggested to measure the sparsity of ${\bm{\Gamma}}$ in a localized manner. More concretely, consider the $i$-th $n$-dimensional patch of the global system ${\mathbf{X}}={\mathbf{D}}{\bm{\Gamma}}$, given by ${\mathbf{x}}_{i}={\bm{\Omega}}{\bm{\gamma}}_{i}$. The stripe-dictionary ${\bm{\Omega}}$, which is of size $n\times(2n-1)m$, is obtained by extracting the $i$-th patch from the global dictionary ${\mathbf{D}}$ and discarding all the zero columns from it. The stripe vector ${\bm{\gamma}}_{i}$ is the corresponding sparse representation of length $(2n-1)m$, containing all coefficients of atoms contributing to ${\mathbf{x}}_{i}$. This relation is illustrated in Figure 4. Notably, the choice of a convolutional dictionary results in signals such that every patch of length $n$ extracted from them can be sparsely represented using a single shift-invariant local dictionary ${\bm{\Omega}}$ – a common assumption usually employed in signal and image processing.

Intuitively, this seeks for a global vector ${\bm{\Gamma}}$ that can represent sparsely every patch in the signal ${\mathbf{X}}$ using the dictionary ${\bm{\Omega}}$. The advantage of the above problem over the traditional ${\text{P}_{0}}$ becomes apparent as we move to consider its theoretical aspects. Assuming that the number of non-zeros per stripe (and not globally) in ${\bm{\Gamma}}$ is less than $\frac{1}{2}\left(1+\frac{1}{\mu({\mathbf{D}})}\right)$, in (Papyan et al., 2016a) it was proven that the solution for the ${\text{P}_{0,\infty}}$ problem is unique. Furthermore, classical pursuit methods, originally tackling the ${\text{P}_{0}}$ problem, are guaranteed to find this representation.

Similar to the ${\text{P}_{0,\infty}}$ problem, this was also analyzed theoretically, shedding light on the theoretical aspects of the convolutional model in the presence of noise. In particular, a stability claim for the ${\text{P}_{0,\infty}^{\bm{\mathbf{\mathcal{E}}}}}$ problem and guarantees for the success of both the OMP and the BP were provided. Similar to the noiseless case, these assumed that the number of non-zeros per stripe is low.

From Atoms to Molecules: Multi-Layer Convolutional Sparse Model

Convolutional sparsity assumes an inherent structure for natural signals. Similarly, the representations themselves could also be assumed to have such a structure. In what follows, we propose a novel layered model that relies on this rationale.

Intuitively, ${\mathbf{X}}={\mathbf{D}}_{1}{\bm{\Gamma}}_{1}$ assumes that the signal ${\mathbf{X}}$ is a superposition of atoms taken from ${\mathbf{D}}_{1}$. While equation ${\mathbf{X}}={\mathbf{D}}_{1}{\mathbf{D}}_{2}{\bm{\Gamma}}_{2}$ views the signal as a superposition of more complex entities taken from the dictionary ${\mathbf{D}}_{1}{\mathbf{D}}_{2}$, which we coin molecules.

While this proposal can be interpreted as a straightforward fusion between the double sparsity model (Rubinstein et al., 2010) and the convolutional one, it is in fact substantially different. The double sparsity model assumes that ${\mathbf{D}}_{2}$ is sparse, and forces only the deepest representation ${\bm{\Gamma}}_{2}$ to be sparse as well. Here, on the other hand, we replace this constraint by forcing ${\mathbf{D}}_{2}$ to have a stride convolution structure, putting emphasis on the sparsity of both the representations ${\bm{\Gamma}}_{1}$ and ${\bm{\Gamma}}_{2}$. In Section 7.1 we will revisit the double sparsity work and its ties to ours by showing the benefits of injecting the assumption on the sparsity of ${\mathbf{D}}_{2}$ into our proposed model.

Under the above construction the sparse vector ${\bm{\Gamma}}_{1}$ has two roles. In the context of the system of equations ${\mathbf{X}}={\mathbf{D}}_{1}{\bm{\Gamma}}_{1}$, it is the convolutional sparse representation of the signal ${\mathbf{X}}$ over the dictionary ${\mathbf{D}}_{1}$. As such, the vector ${\bm{\Gamma}}_{1}$ is composed from $(2n_{0}-1)m_{1}$-dimensional stripes, $\SS_{1,j}{\bm{\Gamma}}_{1}$, where $\SS_{i,j}$ is the operator that extracts the $j$-th stripe from ${\bm{\Gamma}}_{i}$. From another point of view, ${\bm{\Gamma}}_{1}$ is in itself a signal that admits a sparse representation ${\bm{\Gamma}}_{1}={\mathbf{D}}_{2}{\bm{\Gamma}}_{2}$. Denoting by ${\mathbf{P}}_{i,j}$ the operator that extracts the $j$-th patch from ${\bm{\Gamma}}_{i}$, the signal ${\bm{\Gamma}}_{1}$ is composed of patches ${\mathbf{P}}_{1,j}{\bm{\Gamma}}_{1}$ of length $n_{1}m_{1}$. The above model is depicted in Figure 5, presenting both roles of ${\bm{\Gamma}}_{1}$ and their corresponding constituents – stripes and patches. Clearly, the above construction can be extended to more than two layers, leading to the following definition:

For a global signal ${\mathbf{X}}$, a set of convolutional dictionaries $\{{\mathbf{D}}_{i}\}_{i=1}^{K}$, and a vector ${\bm{\lambda}}$, define the deep coding problem ${\text{DCP}_{\bm{\lambda}}}$ as:

where the scalar $\lambda_{i}$ is the $i$-th entry of ${\bm{\lambda}}$.

Denoting ${\bm{\Gamma}}_{0}$ to be the signal ${\mathbf{X}}$, the ${\text{DCP}_{\bm{\lambda}}}$ can be rewritten compactly as

Intuitively, given a signal ${\mathbf{X}}$, this problem seeks for a set of representations, $\{{\bm{\Gamma}}_{i}\}_{i=1}^{K}$, such that each one is locally sparse. As we shall see next, the above can be easily solved using simple algorithms that also enjoy from theoretical justifications. Next, we extend the ${\text{DCP}_{\bm{\lambda}}}$ problem to a noisy regime.

For a global signal ${\mathbf{Y}}$, a set of convolutional dictionaries $\{{\mathbf{D}}_{i}\}_{i=1}^{K}$, and vectors ${\bm{\lambda}}$ and ${\bm{\mathbf{\mathcal{E}}}}$, define the deep coding problem ${\text{DCP}_{\bm{\lambda}}^{\hskip 1.13791pt{\bm{\mathbf{\mathcal{E}}}}}}$ as:

where the scalars $\lambda_{i}$ and $\mathcal{E}_{i}$ are the $i$-th entry of ${\bm{\lambda}}$ and ${\bm{\mathbf{\mathcal{E}}}}$, respectively.

We now move to the task of learning the model parameters. Denote by ${\text{DCP}_{\bm{\lambda}}^{\hskip 1.13791pt\star}}({\mathbf{X}},\{{\mathbf{D}}_{i}\}_{i=1}^{K})$ the representation ${\bm{\Gamma}}_{K}$ obtained by solving the DCP problem (Definition 1, i.e., noiseless) for the signal ${\mathbf{X}}$ and the set of dictionaries $\{{\mathbf{D}}_{i}\}_{i=1}^{K}$. Relying on this, we now extend the dictionary learning problem, as presented in Section 2.2.3, to the multi-layer convolutional sparse representation setting.

A clarification for the chosen name, deep learning problem, will be provided shortly. The solution for the above results in an end-to-end mapping, from a set of input signals to their corresponding labels. Similarly, we can define the ${\text{DLP}_{\bm{\lambda}}^{\hskip 1.13791pt{\bm{\mathbf{\mathcal{E}}}}}}$ problem. However, this is omitted for the sake of brevity. We conclude this section by summarizing, for the convenience of the reader, all notations used throughout this work in Table 1.

Layered Thresholding: The Crux of the Matter

Consider the ML-CSC model defined by the set of dictionaries $\{{\mathbf{D}}_{i}\}_{i=1}^{K}$. Assume we are given a signal

and our goal is to find its underlying representations, $\{{\bm{\Gamma}}_{i}\}_{i=1}^{K}$. Tackling this problem by recovering all the vectors at once might be computationally and conceptually challenging; therefore, we propose the layered thresholding algorithm that gradually computes the sparse vectors one at a time across the different layers. Denoting by ${\mathcal{P}}_{\beta}(\cdot)$ a sparsifying operator that is equal to ${\mathcal{H}}_{\beta}(\cdot)$ in the hard thresholding case and ${\mathcal{S}}_{\beta}(\cdot)$ in the soft one; we commence by computing $\hat{{\bm{\Gamma}}}_{1}={\mathcal{P}}_{\beta_{1}}({\mathbf{D}}_{1}^{T}{\mathbf{X}})$, which is an approximation of ${\bm{\Gamma}}_{1}$. Next, by applying another thresholding algorithm, however this time on $\hat{{\bm{\Gamma}}}_{1}$, an approximation of ${\bm{\Gamma}}_{2}$ is obtained, $\hat{{\bm{\Gamma}}}_{2}={\mathcal{P}}_{\beta_{2}}({\mathbf{D}}_{2}^{T}\hat{{\bm{\Gamma}}}_{1})$. This process, which is iterated until the last representation $\hat{{\bm{\Gamma}}}_{K}$ is acquired, is summarized in Algorithm 1.

One might ponder as to why does the application of the thresholding algorithm on the signal ${\mathbf{X}}$ not result in the true representation ${\bm{\Gamma}}_{1}$, but instead an approximation of it. As previously described in Section 2.2.1, assuming some conditions are met, the result of the thresholding algorithm, $\hat{{\bm{\Gamma}}}_{1}$, is guaranteed to have the correct support. In order to obtain the vector ${\bm{\Gamma}}_{1}$ itself, one should project the signal ${\mathbf{X}}$ onto this obtained support, by solving a Least-Squares problem. For reasons that will become clear shortly, we choose not to employ this step in the layered thresholding algorithm. Despite this algorithm failing in recovering the exact representations in the noiseless setting, as we shall see in Section 5, the estimated sparse vectors and the true ones are close – indicating the stability of this simple algorithm.

Thus far, we have assumed a noiseless setting. However, the same layered thresholding algorithm could be employed for the recovery of the representations of a noisy signal ${\mathbf{Y}}={\mathbf{X}}+{\mathbf{E}}$, with the exception that the threshold constants, $\{\beta_{i}\}_{i=1}^{K}$, would be different and proportional to the noise level.

Assuming two layers for simplicity, Algorithm 1 can be summarized in the following equation

Comparing the above with Equation (1), given by

one can notice a striking similarity between the two. Moreover, by replacing ${\mathcal{P}}_{\beta}(\cdot)$ with the soft nonnegative thresholding, ${\mathcal{S}}_{\beta}^{+}(\cdot)$, we obtain that the aforementioned pursuit and the forward pass of the CNN are equal! Notice that we are relying here on the discussion of Section 2.2.2, where we have shown that the ReLU and the soft nonnegative thresholding are equalA slight difference does exist between the soft nonnegative layered thresholding algorithm and the forward pass of the CNN. While in the former a constant threshold $\beta$ is employed for all entries, the latter uses a bias vector, ${\mathbf{b}}$, that might not be constant in all of its entries. This is of little significance, however, since a similar approach of an entry-based constant could be used in the layered thresholding algorithm as well..

Recall the optimization problem of the training stage of the CNN as shown in Equation (2), given by

and its parallel in the ML-CSC model, the ${\text{DLP}_{\bm{\lambda}}}$ problem, defined by

Notice the remarkable similarity between both objectives, the only difference being in the feature vector on which the classification is done; in the CNN this is the output of the forward pass algorithm, given by $f\left({\mathbf{X}}_{j},\{{\mathbf{W}}_{i}\}_{i=1}^{K},\{{\mathbf{b}}_{i}\}_{i=1}^{K}\right)$, while in the sparsity case this is the result of the ${\text{DCP}_{\bm{\lambda}}}$ problem. In light of the discussion above, the solution for the ${\text{DCP}_{\bm{\lambda}}}$ problem can be approximated using the layered thresholding algorithm, which is in turn equal to the forward pass of the CNN. We can therefore conclude that the problems solved by the training stage of the CNN and the ${\text{DLP}_{\bm{\lambda}}}$ are tightly connected, and in fact are equal once the solution for the ${\text{DLP}_{\bm{\lambda}}}$ is approximated via the layered thresholding algorithm (hence the name ${\text{DLP}_{\bm{\lambda}}}$).

Theoretical Study

Thus far, we have defined the ML-CSC model and its corresponding pursuits – the ${\text{DCP}_{\bm{\lambda}}}$ and ${\text{DCP}_{\bm{\lambda}}^{\hskip 1.13791pt{\bm{\mathbf{\mathcal{E}}}}}}$ problems. We have proposed a method to tackle them, coined the layered thresholding algorithm, which was shown to be equivalent to the forward pass of the CNN. Relying on this, we conclude that the proposed ML-CSC is the global Bayesian model implicitly imposed on the signal ${\mathbf{X}}$ when deploying the forward pass algorithm. Put differently, the ML-CSC answers the question of who are the signals belonging to the model behind the CNN. Having established the importance of our model, we now proceed to its theoretical analysis.

We should emphasize that the following study does not assume any specific form on the network’s parameters, apart from a broad coherence property (as will be shown hereafter). This is in contrast to the work of (Bruna and Mallat, 2013) that assumes that the filters are Wavelets, or the analysis in (Giryes et al., 2015) that considers random weights.

Consider a signal ${\mathbf{X}}$ admitting a multi-layer convolutional sparse representation defined by the sets $\{{\mathbf{D}}_{i}\}_{i=1}^{K}$ and $\{\lambda_{i}\}_{i=1}^{K}$. Can another set of sparse vectors represent the signal ${\mathbf{X}}$? In other words, can we guarantee that, under some conditions, the set $\{{\bm{\Gamma}}_{i}\}_{i=1}^{K}$ is a unique solution to the ${\text{DCP}_{\bm{\lambda}}}$ problem? In the following theorem we provide an answer to this question.

(Uniqueness via the mutual coherence): Consider a signal ${\mathbf{X}}$ satisfying the ${\text{DCP}_{\bm{\lambda}}}$ model,

where $\{{\mathbf{D}}_{i}\}_{i=1}^{K}$ is a set of convolutional dictionaries and $\left\{\mu({\mathbf{D}}_{i})\right\}_{i=1}^{K}$ are their corresponding mutual coherences. If

then the set $\{{\bm{\Gamma}}_{i}\}_{i=1}^{K}$ is the unique solution to the ${\text{DCP}_{\bm{\lambda}}}$ problem, assuming that the thresholds $\left\{\lambda_{i}\right\}_{i=1}^{K}$ are chosen to satisfy

The proof for the above theorem is given in Appendix A. In what follows, we present its importance in the context of CNN. Assume a signal ${\mathbf{X}}$ is fed into a network, resulting in a set of activation values across the different layers. These, in the realm of sparsity, correspond to the set of sparse representations $\{{\bm{\Gamma}}_{i}\}_{i=1}^{K}$, which according to the above theorem are in fact unique representations of the signal ${\mathbf{X}}$.

One might ponder at this point whether there exists an algorithm for obtaining the unique solution guaranteed in this subsection for the ${\text{DCP}_{\bm{\lambda}}}$ problem. As previously mentioned, the layered thresholding algorithm is incapable of finding the exact representations, $\{{\bm{\Gamma}}_{i}\}_{i=1}^{K}$, due to the lack of a Least-Squares step after each layer. One should not despair, however, as we shall see in a following section an alternative algorithm, which manages to overcome this hurdle.

Consider an instance signal ${\mathbf{X}}$ belonging to the ML-CSC model, defined by the sets $\{{\mathbf{D}}_{i}\}_{i=1}^{K}$ and $\{\lambda_{i}\}_{i=1}^{K}$. Assume ${\mathbf{X}}$ is contaminated by a noise vector ${\mathbf{E}}$, generating the perturbed signal ${\mathbf{Y}}={\mathbf{X}}+{\mathbf{E}}$. Suppose we solve the ${\text{DCP}_{\bm{\lambda}}^{\hskip 1.13791pt{\bm{\mathbf{\mathcal{E}}}}}}$ problem and obtain a set of solutions $\{\hat{{\bm{\Gamma}}}_{i}\}_{i=1}^{K}$. How close is every solution in this set, $\hat{{\bm{\Gamma}}}_{i}$, to its corresponding true representation, ${\bm{\Gamma}}_{i}$? In what follows, we provide a theorem addressing this question of stability, the proof of which is deferred to Appendix B.

(Stability of the solution to the ${\text{DCP}_{\bm{\lambda}}^{\hskip 1.13791pt{\bm{\mathbf{\mathcal{E}}}}}}$ problem): Suppose a signal ${\mathbf{X}}$ that has a decomposition

is contaminated with noise ${\mathbf{E}}$, where $\|{\mathbf{E}}\|_{2}\leq\mathcal{E}_{0}$, resulting in ${\mathbf{Y}}={\mathbf{X}}+{\mathbf{E}}$. For all $1\leq i\leq K$, if

$\|{\bm{\Gamma}}_{i}\|_{0,\infty}^{\scriptscriptstyle{\SS}}\leq\lambda_{i}<\frac{1}{2}\left(1+\frac{1}{\mu({\mathbf{D}}_{i})}\right)$; and

${\mathcal{E}_{i}}^{2}=\frac{4\mathcal{E}_{i-1}^{2}}{1-(2\|{\bm{\Gamma}}_{i}\|_{0,\infty}^{\scriptscriptstyle{\SS}}-1)\mu({\mathbf{D}}_{i})}$,

where the set $\{\hat{{\bm{\Gamma}}}_{i}\}_{i=1}^{K}$ is the solution for the ${\text{DCP}_{\bm{\lambda}}^{\hskip 1.13791pt{\bm{\mathbf{\mathcal{E}}}}}}$ problem.

Is this necessarily the true behavior of a deep network? Perhaps the answer to this resides in the choice we made above of considering the noise as adversary. A similar, yet somewhat more involved, analysis with a random noise assumption should be done, with the hope to see a better controlled noise propagation in this system. We leave this for our future work.

Another important remark is that the above bounds the absolute error between the estimated and the true representation. In practice, however, the relative error is of more importance. This is measured in terms of the signal to noise ratio (SNR), which we shall define in Section 8.

Having established the stability of the ${\text{DCP}_{\bm{\lambda}}^{\hskip 1.13791pt{\bm{\mathbf{\mathcal{E}}}}}}$ problem, we now turn to the stability of the algorithms attempting to solve it, the chief one being the forward pass of CNN.

3 Stability of the Layered Hard Thresholding

Define the $\|\cdot\|_{2,\infty}^{\scriptscriptstyle{{\mathbf{P}}}}$ and $\|\cdot\|_{0,\infty}^{\scriptscriptstyle{{\mathbf{P}}}}$ norm of ${\bm{\Gamma}}_{i}$ to be

respectively. The operator ${\mathbf{P}}_{i,j}$ extracts the $j$-th patch of length $n_{i}m_{i}$ from the $i$-th sparse vector ${\bm{\Gamma}}_{i}$.

In the above definition, the letter p emphasizes that the norms are computed by sweeping over all patches, rather than stripes. Recall that we have defined $m_{0}=1$, since the number of channels in the input signal ${\mathbf{X}}={\bm{\Gamma}}_{0}$ is equal to one.

(Stable recovery of hard thresholding in the presence of noise): Suppose a clean signal ${\mathbf{X}}$ has a convolutional sparse representation ${\mathbf{D}}_{1}{\bm{\Gamma}}_{1}$, and that it is contaminated with noise ${\mathbf{E}}$ to create the signal ${\mathbf{Y}}={\mathbf{X}}+{\mathbf{E}}$, such that $\|{\mathbf{E}}\|_{2,\infty}^{\scriptscriptstyle{{\mathbf{P}}}}\leq\epsilon_{0}$. Denote by $|\Gamma_{1}^{\text{min}}|$ and $|\Gamma_{1}^{\text{max}}|$ the lowest and highest entries in absolute value in ${\bm{\Gamma}}_{1}$, respectively. Denote further by $\hat{{\bm{\Gamma}}}_{1}$ the solution obtained by running the hard thresholding algorithm on ${\mathbf{Y}}$ with a constant $\beta_{1}$, i.e. $\hat{{\bm{\Gamma}}}_{1}={\mathcal{H}}_{\beta_{1}}({\mathbf{D}}_{1}^{T}{\mathbf{Y}})$. Assuming that

$\|{\bm{\Gamma}}_{1}\|_{0,\infty}^{\scriptscriptstyle{\SS}}<\frac{1}{2}\left(1+\frac{1}{\mu({\mathbf{D}}_{1})}\frac{|\Gamma_{1}^{\text{min}}|}{|\Gamma_{1}^{\text{max}}|}\right)-\frac{1}{\mu({\mathbf{D}}_{1})}\frac{\epsilon_{0}}{|\Gamma_{1}^{\text{max}}|}$; and

The threshold $\beta_{1}$ is chosen according to Equation (87) (see below),

The support of the solution $\hat{{\bm{\Gamma}}}_{1}$ is equal to that of ${\bm{\Gamma}}_{1}$; and

\|{\bm{\Gamma}}_{1}-\hat{{\bm{\Gamma}}}_{1}\|_{2,\infty}^{\scriptscriptstyle{{\mathbf{P}}}}\leq\sqrt{\|{\bm{\Gamma}}_{1}\|_{0,\infty}^{\scriptscriptstyle{{\mathbf{P}}}}}\Big{(}\epsilon_{0}+\mu({\mathbf{D}}_{1})\left(\|{\bm{\Gamma}}_{1}\|_{0,\infty}^{\scriptscriptstyle{\SS}}-1\right)|\Gamma_{1}^{\text{max}}|\Big{)}.

Notice that by plugging $\epsilon_{0}=0$ the above theorem covers the noiseless scenario. Notably, even in such a case, we obtain a deviation from the true representation due to the lack of a Least-Squares step.

(Stability of layered hard thresholding in the presence of noise): Suppose a clean signal ${\mathbf{X}}$ has a decomposition

and that it is contaminated with noise ${\mathbf{E}}$ to create the signal ${\mathbf{Y}}={\mathbf{X}}+{\mathbf{E}}$, such that $\|{\mathbf{E}}\|_{2,\infty}^{\scriptscriptstyle{{\mathbf{P}}}}\leq\epsilon_{0}$. Denote by $|\Gamma_{i}^{\text{min}}|$ and $|\Gamma_{i}^{\text{max}}|$ the lowest and highest entries in absolute value in the vector ${\bm{\Gamma}}_{i}$, respectively. Let $\{\hat{{\bm{\Gamma}}}_{i}\}_{i=1}^{K}$ be the set of solutions obtained by running the layered hard thresholding algorithm with thresholds $\{\beta_{i}\}_{i=1}^{K}$, i.e. $\hat{{\bm{\Gamma}}}_{i}={\mathcal{H}}_{\beta_{i}}({\mathbf{D}}_{i}^{T}\hat{{\bm{\Gamma}}}_{i-1})$ where $\hat{{\bm{\Gamma}}}_{0}={\mathbf{Y}}$. Assuming that $\forall\ 1\leq i\leq K$

$\|{\bm{\Gamma}}_{i}\|_{0,\infty}^{\scriptscriptstyle{\SS}}<\frac{1}{2}\left(1+\frac{1}{\mu({\mathbf{D}}_{i})}\frac{|\Gamma_{i}^{\text{min}}|}{|\Gamma_{i}^{\text{max}}|}\right)-\frac{1}{\mu({\mathbf{D}}_{i})}\frac{\epsilon_{i-1}}{|\Gamma_{i}^{\text{max}}|}$; and

The threshold $\beta_{i}$ is chosen according to Equation (109),

The support of the solution $\hat{{\bm{\Gamma}}}_{i}$ is equal to that of ${\bm{\Gamma}}_{i}$; and

$\|{\bm{\Gamma}}_{i}-\hat{{\bm{\Gamma}}}_{i}\|_{2,\infty}^{\scriptscriptstyle{{\mathbf{P}}}}\leq\epsilon_{i}$,

where \epsilon_{i}=\sqrt{\|{\bm{\Gamma}}_{i}\|_{0,\infty}^{\scriptscriptstyle{{\mathbf{P}}}}}\ \Big{(}\epsilon_{i-1}+\mu({\mathbf{D}}_{i})\left(\|{\bm{\Gamma}}_{i}\|_{0,\infty}^{\scriptscriptstyle{\SS}}-1\right)|\Gamma_{i}^{\text{max}}|\Big{)}.

The proof for the above is given in Appendix D. We now turn to an analogous theorem for the forward pass of the CNN, prior to discussing the surprising implications of these theorems.

4 Stability of the Forward Pass (Layered Soft Thresholding)

In light of the discussion in Section 4, the equivalence between the layered thresholding algorithm and the forward pass of the CNN is achieved assuming that the operator employed is the nonnegative soft thresholding ${\mathcal{S}}_{\beta}^{+}(\cdot)$. However, thus far, we have analyzed the closely related hard version ${\mathcal{H}}_{\beta}(\cdot)$ instead. In what follows, we show how the stability theorem presented in the previous subsection can be modified to the soft version, ${\mathcal{S}}_{\beta}(\cdot)$. For simplicity, and in order to stay in line with the vast sparse representation theory, herein we choose not to assume the nonnegative assumption. This implies that we are proposing a slightly different CNN architecture in which the ReLU function is two sided (Kavukcuoglu et al., 2010). We now move to the stable recovery of the soft thresholding algorithm.

(Stable recovery of soft thresholding in the presence of noise): Suppose a clean signal ${\mathbf{X}}$ has a convolutional sparse representation ${\mathbf{D}}_{1}{\bm{\Gamma}}_{1}$, and that it is contaminated with noise ${\mathbf{E}}$ to create the signal ${\mathbf{Y}}={\mathbf{X}}+{\mathbf{E}}$, such that $\|{\mathbf{E}}\|_{2,\infty}^{\scriptscriptstyle{{\mathbf{P}}}}\leq\epsilon_{0}$. Denote by $|\Gamma_{1}^{\text{min}}|$ and $|\Gamma_{1}^{\text{max}}|$ the lowest and highest entries in absolute value in ${\bm{\Gamma}}_{1}$, respectively. Denote further by $\hat{{\bm{\Gamma}}}_{1}$ the solution obtained by running the soft thresholding algorithm on ${\mathbf{Y}}$ with a constant $\beta_{1}$, i.e. $\hat{{\bm{\Gamma}}}_{1}={\mathcal{S}}_{\beta_{1}}({\mathbf{D}}_{1}^{T}{\mathbf{Y}})$. Assuming that

$\|{\bm{\Gamma}}_{1}\|_{0,\infty}^{\scriptscriptstyle{\SS}}<\frac{1}{2}\left(1+\frac{1}{\mu({\mathbf{D}}_{1})}\frac{|\Gamma_{1}^{\text{min}}|}{|\Gamma_{1}^{\text{max}}|}\right)-\frac{1}{\mu({\mathbf{D}}_{1})}\frac{\epsilon_{0}}{|\Gamma_{1}^{\text{max}}|}$; and

The threshold $\beta_{1}$ is chosen according to Equation (87),

The support of the solution $\hat{{\bm{\Gamma}}}_{1}$ is equal to that of ${\bm{\Gamma}}_{1}$; and

\|{\bm{\Gamma}}_{1}-\hat{{\bm{\Gamma}}}_{1}\|_{2,\infty}^{\scriptscriptstyle{{\mathbf{P}}}}\leq\sqrt{\|{\bm{\Gamma}}_{1}\|_{0,\infty}^{\scriptscriptstyle{{\mathbf{P}}}}}\Big{(}\epsilon_{0}+\mu({\mathbf{D}}_{1})\left(\|{\bm{\Gamma}}_{1}\|_{0,\infty}^{\scriptscriptstyle{\SS}}-1\right)|\Gamma_{1}^{\text{max}}|+\beta_{1}\Big{)}.

Armed with the above lemma, which is proven in Appendix E, we now proceed to the stability of the forward pass of the CNN.

(Stability of the forward pass (layered soft thresholding algorithm) in the presence of noise): Suppose a clean signal ${\mathbf{X}}$ has a decomposition

and that it is contaminated with noise ${\mathbf{E}}$ to create the signal ${\mathbf{Y}}={\mathbf{X}}+{\mathbf{E}}$, such that $\|{\mathbf{E}}\|_{2,\infty}^{\scriptscriptstyle{{\mathbf{P}}}}\leq\epsilon_{0}$. Denote by $|\Gamma_{i}^{\text{min}}|$ and $|\Gamma_{i}^{\text{max}}|$ the lowest and highest entries in absolute value in the vector ${\bm{\Gamma}}_{i}$, respectively. Let $\{\hat{{\bm{\Gamma}}}_{i}\}_{i=1}^{K}$ be the set of solutions obtained by running the layered soft thresholding algorithm with thresholds $\{\beta_{i}\}_{i=1}^{K}$, i.e. $\hat{{\bm{\Gamma}}}_{i}={\mathcal{S}}_{\beta_{i}}({\mathbf{D}}_{i}^{T}\hat{{\bm{\Gamma}}}_{i-1})$ where $\hat{{\bm{\Gamma}}}_{0}={\mathbf{Y}}$. Assuming that $\forall\ 1\leq i\leq K$

$\|{\bm{\Gamma}}_{i}\|_{0,\infty}^{\scriptscriptstyle{\SS}}<\frac{1}{2}\left(1+\frac{1}{\mu({\mathbf{D}}_{i})}\frac{|\Gamma_{i}^{\text{min}}|}{|\Gamma_{i}^{\text{max}}|}\right)-\frac{1}{\mu({\mathbf{D}}_{i})}\frac{\epsilon_{i-1}}{|\Gamma_{i}^{\text{max}}|}$; and

The threshold $\beta_{i}$ is chosen according to Equation (109) (with the $\epsilon_{i}$ defined below),

The support of the solution $\hat{{\bm{\Gamma}}}_{i}$ is equal to that of ${\bm{\Gamma}}_{i}$; and

$\|{\bm{\Gamma}}_{i}-\hat{{\bm{\Gamma}}}_{i}\|_{2,\infty}^{\scriptscriptstyle{{\mathbf{P}}}}\leq\epsilon_{i}$,

where \epsilon_{i}=\sqrt{\|{\bm{\Gamma}}_{i}\|_{0,\infty}^{\scriptscriptstyle{{\mathbf{P}}}}}\ \Big{(}\epsilon_{i-1}+\mu({\mathbf{D}}_{i})\left(\|{\bm{\Gamma}}_{i}\|_{0,\infty}^{\scriptscriptstyle{\SS}}-1\right)|\Gamma_{i}^{\text{max}}|+\beta_{i}\Big{)}.

The above theorem guarantees that the distances between the original representations and the ones obtained from the CNN are bounded. Even if we set $\epsilon_{0}=0$, the recovered activations deviate from the true ones, simply because the layered thresholding algorithm does not do a perfect job, even on a noiseless signal. When the signal is noisy, these deviations are strengthened, but still in a controlled way.

This, by itself, might not be surprising. After all, the CNN is a deterministic system of linear operations (convolutions), followed by simple non-linearities that are non-expanding. If we feed a slightly perturbed signal to such a system, it is clear that the activations all along the network will be perturbed as well with a bounded effect. However, the above theorem shows far more than that. There are, in fact, two types of stabilities, the trivial one that considers the sensitivity of the whole feed-forward network to perturbations in its input, and the more intricate one that shows that this system enables a rather accurate recovery of the generating representations. The second option is the stability we prove here.

5 Guarantees for Fully Connected Networks

One should note that the convolutional structure imposed on the dictionaries in our model could be removed, and the theoretical guarantees we have provided above would still hold. The reason being is that the unconstrained dictionary can be regarded as a convolutional one, constructed from a single shift of a local matrix with no circular boundary. In the context of CNN, this is analogous to a fully connected layer. As such, the theoretical analysis provided here sheds light on both convolutional and fully connected networks. A different point of view on the same matter can also be proposed; fully connected layers can be viewed as convolutional ones with filters that cover their entire input (Long et al., 2015).

Layered Basis Pursuit – The Future of Deep Learning?

The stability analysis presented above unveils two significant limitations of the forward pass of the CNN. First, this algorithm is incapable of recovering the unique solution for the ${\text{DCP}_{\bm{\lambda}}}$ problem, the existence of which is guaranteed from Theorem 4. This acts against our expectations, since in the traditional sparsity inspired model it is a well known fact that such a unique representation can be retrieved, assuming certain conditions are met.

A solution for the first problem, already presented throughout this work, is a two-stage approach. First, run the thresholding operator in order to recover the correct support. Then, once the atoms are chosen, their corresponding coefficients can be obtained by solving a linear system of equations. In addition to retrieving the true representation in the noiseless case, this step can also be beneficial in the noisy scenario, resulting in a solution closer to the underlying one. However, since no such step exists in current CNN architectures, we refrain from further analyzing its theoretical implications.

Next, we present an alternative to the layered soft thresholding algorithm, which will tackle both of the aforementioned problems. Recall that the result of the soft thresholding is a simple approximation of the solution for the ${\text{P}_{1}}$ problem, previously defined in Equation (4). In every layer, instead of applying a simple thresholding operator that estimates the sparse vector by computing $\hat{{\bm{\Gamma}}}_{i}={\mathcal{S}}_{\beta_{i}}({\mathbf{D}}_{i}^{T}\hat{{\bm{\Gamma}}}_{i-1})$; we propose to tackle the full pursuit, i.e. to minimize

Notice that one could readily obtain the nonnegative sparse coding problem by simply adding an extra constraint in the above equation, forcing the coefficients in ${\bm{\Gamma}}_{i}$ to be nonnegative. More generally, Equation (53) can be written in its Lagrangian formulation

where the constant $\xi_{i}$ is proportional to the noise level and should tend to zero in the noiseless scenario. We name the above the layered basis pursuit (BP) algorithm. In practice, one possible method for solving it is the iterative soft thresholding (IST). Formally, this obtains the minimizer of Equation (54) by repeating the following recursive formula

where $\hat{{\bm{\Gamma}}}_{i}^{t}$ is the estimate of ${\bm{\Gamma}}_{i}$ at iteration $t$. The above can be interpreted as a simple projected gradient descent algorithm, where the constant $c_{i}$ is inversely proportional to its step size. As a result, if $c_{i}$ is chosen to be large enoughThe constant $c_{i}$ should satisfy $c_{i}>0.5\lambda_{\text{max}}\left({\mathbf{D}}_{i}^{T}{\mathbf{D}}_{i}\right)$, where $\lambda_{\text{max}}\left({\mathbf{D}}_{i}^{T}{\mathbf{D}}_{i}\right)$ is the maximal eigenvalue of the gram matrix ${\mathbf{D}}_{i}^{T}{\mathbf{D}}_{i}$ (Combettes and Wajs, 2005)., the above algorithm is guaranteed to converge to its global minimum that is the solution of (54), as was shown in (Daubechies et al., 2004). The method obtained by gradually computing the set of sparse representations, $\{{\bm{\Gamma}}_{i}\}_{i=1}^{K}$, via the IST is summarized in Algorithm 2 and named layered iterative soft thresholding. Notice that this algorithm coincides with the simple layered soft thresholding if it is run for a single iteration with $c_{i}=1$ and initialized with $\hat{{\bm{\Gamma}}}_{i}^{0}=\mathbf{0}$. This implies that the above algorithm is a natural extension to the forward pass of the CNN.

With respect to the computational aspects of the IST algorithm, the work of (Gregor and LeCun, 2010) proposed the LISTA method, showing how the number of iterations required by the IST to convergence can be reduced using neural networks. Analogously, the work of (Xin et al., 2016) presented a generalization of the iterative hard thresholding (IHT), which was shown both theoretically and empirically to be superior to the original IHT.

The original motivation for the layered IST was its theoretical superiority over the forward pass algorithm – one that will be explored in detail in the next subsection. Yet more can be said about this algorithm and the CNN architecture it induces. In (Gregor and LeCun, 2010) it was shown that the IST algorithm can be formulated as a simple recurrent neural network. As such, the same can be said regarding the layered IST algorithm proposed here, with the exception that the induced recurrent network is much deeper. The reader can therefore interpret this part of the work as a theoretical study of a special case of recurrent neural networks.

From another perspective, the underlying architecture of the layered IST algorithm is a cascade of $K$ blocks. Each of these corresponds to a fixed number of unfolded iterations, $T_{i}$, of a single IST algorithm. As such, it contains several convolutional layers with shared weights, as well as skip connections in order to compute the residual, $\hat{{\bm{\Gamma}}}_{i-1}-{\mathbf{D}}_{i}\hat{{\bm{\Gamma}}}_{i}^{t-1}$, as defined in Equation (55). Interestingly, the above description is reminiscent (though not exact) of residual networks (He et al., 2015), which have recently led to state-of-the-art results in image recognition.

where $\{{\mathbf{D}}_{i}\}_{i=1}^{K}$ is a set of convolutional dictionaries and $\left\{\mu({\mathbf{D}}_{i})\right\}_{i=1}^{K}$ are their corresponding mutual coherences. Assuming that $\forall\ 1\leq i\leq K$

then the layered BP algorithm is guaranteed to recover the set $\{{\bm{\Gamma}}_{i}\}_{i=1}^{K}$.

2 Stability of Layered BP Algorithm

Having established the guarantee for the success of the layered BP algorithm, we now move to its stability analysis. In particular, in a noisy scenario where obtaining the true underlying representations is impossible, does this algorithm remain stable? If so, how do its guarantees compare to those of the layered thresholding algorithm? The following theorem, which we prove in Appendix F, aims to answer these questions.

(Stability of the layered BP algorithm in the presence of noise): Suppose a clean signal ${\mathbf{X}}$ has a decomposition

and that it is contaminated with noise ${\mathbf{E}}$ to create the signal ${\mathbf{Y}}={\mathbf{X}}+{\mathbf{E}}$, such that $\|{\mathbf{E}}\|_{2,\infty}^{\scriptscriptstyle{{\mathbf{P}}}}\leq\epsilon_{0}$. Let $\{\hat{{\bm{\Gamma}}}_{i}\}_{i=1}^{K}$ be the set of solutions obtained by running the layered BP algorithm with parameters $\{\xi_{i}\}_{i=1}^{K}$. Assuming that $\ \forall\ 1\leq i\leq K$

$\|{\bm{\Gamma}}_{i}\|_{0,\infty}^{\scriptscriptstyle{\SS}}<\frac{1}{3}\left(1+\frac{1}{\mu({\mathbf{D}}_{i})}\right)$; and

The support of the solution $\hat{{\bm{\Gamma}}}_{i}$ is contained in that of ${\bm{\Gamma}}_{i}$;

$\|{\bm{\Gamma}}_{i}-\hat{{\bm{\Gamma}}}_{i}\|_{2,\infty}^{\scriptscriptstyle{{\mathbf{P}}}}\leq\epsilon_{i}$;

In particular, every entry of ${\bm{\Gamma}}_{i}$ greater in absolute value than $\frac{\epsilon_{i}}{\sqrt{\|{\bm{\Gamma}}_{i}\|_{0,\infty}^{\scriptscriptstyle{{\mathbf{P}}}}}}$ is guaranteed to be recovered; and

The solution $\hat{{\bm{\Gamma}}}_{i}$ is the unique minimizer of the Lagrangian BP problem (Equation (54)),

where $\epsilon_{i}=\|{\mathbf{E}}\|_{2,\infty}^{\scriptscriptstyle{{\mathbf{P}}}}\ 7.5^{i}\ \prod_{j=1}^{i}\sqrt{\|{\bm{\Gamma}}_{j}\|_{0,\infty}^{\scriptscriptstyle{{\mathbf{P}}}}}$.

Several remarks are due at this point. The condition for the stability of the layered thresholding algorithm, given by

is expected to be more strict than that of the theorem presented above, which is

In addition, similar to the stability analysis presented in Section 5.2, the above shows the growth (as a function of the depth) of the distance between the recovered representations and the true ones.

A Closer Look at the Proposed Model

In this section, we revisit the assumptions of our model by imposing additional constraints on the dictionaries involved and showing their theoretical benefits. These additional assumptions originate from the current common practice of both CNN and sparsity.

Consider the representation ${\bm{\Gamma}}_{K-1}$, given by

where ${\bm{\Omega}}_{K}$ is the stripe-dictionary of ${\mathbf{D}}_{K}$, the vector ${\mathbf{P}}_{K-1,i}{\bm{\Gamma}}_{K-1}$ is the $i$-th patch in ${\bm{\Gamma}}_{K-1}$ and ${\bm{\gamma}}_{K,i}$ is its corresponding stripe. Recalling the definition of the $\|\cdot\|_{0,\infty}^{\scriptscriptstyle{{\mathbf{P}}}}$ norm (Definition 6 in Section 5.3), we have that

The multiplication ${\bm{\Omega}}_{K}\hskip 1.42271pt{\bm{\gamma}}_{K,i}$ can be seen as a linear combination of at most $\|{\bm{\gamma}}_{K,i}\|_{0}$ atoms, each contributing no more than $\|{\bm{\Omega}}_{K}\|_{0}$ non-zeros. As such

Noticing that $\|{\bm{\Omega}}_{K}\|_{0}=\|{\mathbf{D}}_{K}\|_{0}$ (as can be seen in Figure 4), and using the definition of the $\|\cdot\|_{0,\infty}^{\scriptscriptstyle{\SS}}$ norm, we conclude that

In other words, given $\|{\bm{\Gamma}}_{K}\|_{0,\infty}^{\scriptscriptstyle{\SS}}$ and $\|{\mathbf{D}}_{K}\|_{0}$, we can bound the maximal number of non-zeros in a patch from ${\bm{\Gamma}}_{K-1}$.

The claims in Section 5 and 6 are given in terms of not only $\|{\bm{\Gamma}}_{K-1}\|_{0,\infty}^{\scriptscriptstyle{{\mathbf{P}}}}$, but also $\|{\bm{\Gamma}}_{K-1}\|_{0,\infty}^{\scriptscriptstyle{\SS}}$. According to Table 1, the length of a patch in ${\bm{\Gamma}}_{K-1}$ is $n_{K-1}m_{K-1}$, while the size of a stripe is $(2n_{K-2}-1)m_{K-1}$. As such, we can fit $(2n_{K-2}-1)/n_{K-1}$ patches in a stripe. Assume for simplicity that this ratio is equal to one. As a result, we obtain that a patch in the signal ${\bm{\Gamma}}_{K-1}$ extracted from the system

is also a stripe in the representation ${\bm{\Gamma}}_{K-1}$ when considering

hence the name of this subsection. Leveraging this assumption, we return to Equation (71) and obtain that

Using the same rationale for the remaining layers, and assuming that once again the patches become stripes, we conclude that

2 On the Role of the Spatial-Stride

A common step among practitioners of CNN (Krizhevsky et al., 2012; Simonyan and Zisserman, 2014; He et al., 2015) is to convolve the input to each layer with a set of filters, skipping a fixed number of spatial locations in a regular pattern. One of the primary motivations for this is to reduce the dimensions of the kernel maps throughout the layers, leading to computational benefits. In this subsection we unveil some theoretical benefits of this common practice, which we coin spatial-stride.

Note that while the original $\|{\bm{\Gamma}}_{i}\|_{0,\infty}^{\scriptscriptstyle{\SS}}$ is equal to the maximal number of non-zeros in a stripe of length $(2n_{i-1}-1)m_{i}$ in ${\bm{\Gamma}}_{i}$, the term $\|{\mathbf{Q}}_{i}{\bm{\Gamma}}_{i}\|_{0,\infty}^{\scriptscriptstyle{\SS}}$ counts the same quantity but for stripes of length $(2\ \left\lceil n_{i-1}/s_{i-1}\right\rceil-1)m_{i}$ in ${\mathbf{Q}}_{i}{\bm{\Gamma}}_{i}$.

According to the study in Section 5 and 6, the theoretical advantage of the spatial-stride is twofold. First, consider the mutual coherence of the stride convolutional dictionary ${\mathbf{D}}_{i}$. Due to the locality of the filters and their restriction to certain spatial shifts, the mutual coherence of ${\mathbf{D}}_{i}{\mathbf{Q}}_{i}^{T}$ is expected to be lower than that of ${\mathbf{D}}_{i}$, thus leading to more non-zeros allowed per stripe. Second, the length of a stripe in ${\mathbf{Q}}_{i}{\bm{\Gamma}}_{i}$ is equal to $(2\ \left\lceil n_{i-1}/s_{i-1}\right\rceil-1)m_{i}$, while that of ${\bm{\Gamma}}_{i}$ is $(2n_{i-1}-1)m_{i}$. As such, our analysis allows a larger number of non-zeros per a smaller-sized stripe. From another perspective, notice that imposing a spatial-stride on the dictionary ${\mathbf{D}}_{i}$ is equivalent to forcing a portion of the entries in ${\bm{\Gamma}}_{i}$ to be zero. As such, the spatial-stride encourages sparser solutions.

Experiments: The Generator Behind the CNN

In this section we combine the above notions in order to achieve our goal – generate a set of signals that will satisfy the ML-CSC assumptions. These will then serve as a playground for several experiments, which will compare both theoretically and practically the different pursuits presented in this paper.

We commence by describing the design of the dictionaries, and in the next subsection continue to the actual generation of the signals. In our experiments, the signal is one dimensional and therefore $m_{0}=1$. Moreover, for simplicity, the dictionary in every layer contains a single atom with its shifts and thus $m_{i}=1\ \forall 1\leq i\leq K$. We should note that the choice of a single atom simplifies the involved pursuit problem, but as we will see, even in such a case the suggested layered pursuits (including the forward pass) may fail. This is because the mutual coherence and the amount of non-zeros are still non-trivial.

In the first layer we choose this filter to be the analytically defined discrete Meyer Wavelet of length $n_{0}=29$. In order to obtain sparser representations and improve the coherence of the global dictionary ${\mathbf{D}}_{1}$, we employ a stride of $s_{0}=6$, resulting in $\mu({\mathbf{D}}_{1})=2.44\times 10^{-4}$. As a consequence of our choice of ${\mathbf{D}}_{1}$, the signals resulting from our model are a superposition of shifted versions of discrete Meyer Wavelets, multiplied by different coefficients.

2 Noiseless Experiments

Given the signals, we attempt to retrieve their underlying representations using the layered pursuits presented in this work. Recall that our analysis in Section 5 and 6 indicates that the layered hard thresholding is superior to its soft counterpart, which is equivalent to the forward pass, and that the layered BP is even better than both of these algorithms. We now turn to asserting this claim empirically. While doing so, we aim to study the gap between the theoretical guarantees presented throughout our paper and the empirical performance obtained in practice.

Several remarks are due here. First and foremost, the theoretical bounds indeed hold, since the blue points are above their corresponding green ones and the correct supports are always recovered. Second, our analysis predicts that the distance between the estimated sparse representation, $\hat{{\bm{\Gamma}}}_{i}$, and the true ones, ${\bm{\Gamma}}_{i}$, should increase with the layer. This is evident by the decrease in the values of the green points with the layers. The empirical results presented here (blue dots) corroborate this prognosis, as the error in both algorithms is lowest in the first layer and highest in the lastInterestingly, the error in the layered hard thresholding algorithm is approximately equal in the second and third layers.. Third, our analysis suggests that the layered hard thresholding algorithm should be superior to its soft counterpart. Once again, this can be deduced from the figure by comparing the values of the green points in both of the algorithms. The empirical results presented in Figure 6 confirm this behavior, as can be clearly seen by comparing the errors (blue points) obtained by both algorithms in the $i$-th layer. One should note that the performance gap exhibited here is due to the constant $\beta_{i}$ being subtracted from every entry in the soft thresholding algorithm.

The implications of the above discussion might be troubling in the context of CNN, as what this experiment shows is a deterioration of the empirical SNR throughout the layers of the network. Is this truly the behavior of CNN? Recall that in practice the biases of the different layers (thresholds) are learned in order to achieve the best possible performance in solving a certain task. As such, it might be possible that the decline in SNR presented here is alleviated when better thresholds are employed in lieu of the theoretical ones used thus far. We demonstrate this by running the layered soft thresholding algorithm with an oracle parameter, chosen to be the minimal threshold that leads to $\|{\bm{\Gamma}}_{i}\|_{0}$ non-zeros being chosen in the estimated sparse representation $\hat{{\bm{\Gamma}}}_{i}$. The results for this are presented in Figure 6 and colored in red. Indeed, we observe that this better choice of parameters improves the empirical performance of the layered soft thresholding algorithm and leads to a slower decline in SNR. Still, the performance of the layered soft thresholding is inferior to that of its hard variantNote that in the layered hard thresholding, as long as the correct support is chosen, the threshold does not affect the error and as such the oracle version for it is meaningless., as can be seen by comparing the red points with the blue ones in the subplots below.

Next, we proceed our experiments by running the layered BP algorithm, as defined in Section 6, on the same set of signals. Recall that one of the prime motivations for proposing this algorithm was its ability to retrieve the exact underlying representations, as justified theoretically in Theorem 11. In our experiments, we validate this claim by checking that its conditions hold for each signal and that the underlying representations are indeed retrieved. We omit showing a plot for this and comparing it to the layered thresholding algorithms since the errors obtained are simply zeros.

3 Noisy Experiments

Having established the stability of our proposed algorithms in a perfect scenario, where $\epsilon_{0}=0$, we now turn to a noisy setting. Naturally, the estimation task becomes now even more challenging – not only does the SNR drop with each layer, as demonstrated previously, but also the input SNR is no longer infinity. In order to facilitate the success of our algorithms, in this section we demonstrate the empirical performance and theoretical bounds on $K=2$ layers and a small noise level.

We present the obtained results in terms of the local SNR in Figure 7, showing the stability of the different algorithms that is in accordance with our theoretical bounds. Similar to the noiseless experiment, we observe that for all the algorithms the error increases both theoretically (green points) and empirically (blue points) with the layer depth. As previously discussed, a performance gap exists between the soft and hard layered thresholding algorithms. To mitigate this, we run the layered soft thresholding with an oracle parameter and compare the obtained errors (red points) to those of the other algorithms. The results, depicted in the same figure, show a clear improvement in the performance.

Interestingly, although theoretically superior, the layered BP leads to similar performance to that of the layered soft thresholding and worse performance than that of the layered hard thresholding (when comparing the blue points). We attribute this phenomenon to the suboptimal choice of the parameter $\xi_{i}$, which was chosen thus far according to our theoretical analysis. To validate this suspicion, we run the layered BP with hand-picked $\xi_{i}$ and plot the obtained SNR in red in Figure 7. Not only are the correct supports retrieved for all the signals, but we can also see a clear improvement in terms of the SNR. In the first layer, the layered BP outperforms the layered soft thresholding and leads to similar results to those of the layered hard thresholding, while in the second, the layered BP significantly outperforms both of the other pursuit algorithms.

Next, each signal ${\mathbf{X}}$ is contaminated with a zero-mean white additive Gaussian noise ${\mathbf{E}}$, resulting in a noisy signal ${\mathbf{Y}}={\mathbf{X}}+{\mathbf{E}}$. The average SNR of the noisy signals obtained is $124.43$ dB. Note that this is a weak noise, chosen due to the deterioration of the SNR throughout the layers (one that is worsened when the theoretical parameters are employed). The signals are then fed into the layered BP algorithm, resulting in a set of estimated sparse representations, $\hat{{\bm{\Gamma}}}_{i}$. The parameters $\xi_{i}$ employed are the theoretically justified ones, $\xi_{i}=4\epsilon_{i-1}$. We should note that in our experiments we attempted to run the layered thresholding algorithms, however, as our theory predicts these failed in recovering the correct supports.

Given the estimated representations, we compute the errors $\|\hat{{\bm{\Gamma}}}_{i}-{\bm{\Gamma}}_{i}\|_{2,\infty}^{\scriptscriptstyle{{\mathbf{P}}}}$ and compare these to their corresponding theoretical bounds, obtained from Theorem 12. In addition, we verify that the retrieved supports are contained in the true one, as the theorem guarantees. In practice, we obtain that the layered BP always finds the full support. The obtained results are depicted in Figure 8 in terms of the local SNR. For comparison, we run the layered BP with hand-picked $\xi_{i}$ and present the obtained results in the same figure. We conclude that the layered BP remains stable despite the poor coefficient ratio, unlike the layered thresholding algorithms. Moreover, tuning the $\xi_{i}$ results in a much better performance, similar to what we have seen in the previous experiment.

At this point, one might ponder as to whether the hurdle of poor coefficient ratio is one that the layered soft thresholding (forward pass) can not overcome. We believe that several ideas currently used in CNN, such as Batch Normalization (Ioffe and Szegedy, 2015) or Local Response Normalization (Krizhevsky et al., 2012), are tightly connected to this problem. However, their exact relation to this issue and its theoretical analysis is a matter of future work.

Conclusion

Definition: “A guiding question is the fundamental query that directs the search for understanding” (Traver, 1998). In this work our guiding question was who are the signals that the CNN architecture is designed for? To answer this we have defined the ML-CSC model, for which the thresholding pursuit is nothing but the forward pass of the CNN. Although nothing promises that the forward pass will lead to the original representation of a signal emerging from the ML-CSC model, we have shown this is indeed the case. Having established the relevance of our model to CNN, we then turned to its theoretical analysis. In particular, we provided guarantees for the uniqueness of the feature maps CNN aims to recover, and the stability of the problem CNN aims to solve.

Inspired by the evolution of the pursuit methods in the theory of Sparse-Land, we continued our work by proposing the layered BP algorithm. In the noiseless case, this was theoretically shown to be capable of finding the unique solution of the deep coding problem, the existence of which has been also guaranteed; while in the noisy setting, we have proved the stability of this algorithm.

We analyzed the theoretical benefits of two popular ideas employed in the CNN community, namely the use of sparse filters and the spatial-stride. Leveraging those, we then generated signals satisfying the ML-CSC assumptions and demonstrated the performance of the pursuits presented throughout this work.

We conclude this work by presenting our ongoing research directions:

Through this paper we have assumed the worst – an adversary noise. Can our theoretical analysis be extended to a setting where the noise is random?

Thus far in tackling the deep coding problem, we have restricted ourself to existing methods, such as the forward pass of the CNN or deconvolutional networks (Zeiler et al., 2010). Can we suggest better approximations for the solution of this problem?

Clearly a relation exists between our proposed layered iterative thresholding algorithm and the current throne holder in the task of image recognition – residual networks (He et al., 2015). Can our theory reveal the benefits of introducing skip connections to a CNN?

What is the role of common tricks currently employed in CNN in the context of the ML-CSC model? These include but are not limited to, Batch Normalization (Ioffe and Szegedy, 2015), Local Response Normalization (Krizhevsky et al., 2012), Dropout (Srivastava et al., 2014) and Pooling (LeCun et al., 1990; Krizhevsky et al., 2012; Simonyan and Zisserman, 2014).

The research leading to these results has received funding from the European Research Council under European Union’s Seventh Framework Programme, ERC Grant agreement no. 320649. The authors would like to thank Jeremias Sulam for the inspiring discussions and creative advice.

A Uniqueness via the Mutual Coherence (Proof of Theorem 4)

Let $\{\hat{{\bm{\Gamma}}}_{i}\}_{i=1}^{K}$ be a set of representations of the signal ${\mathbf{X}}$, obtained by solving the ${\text{DCP}_{\bm{\lambda}}}$ problem. According to our assumptions, $\|{\bm{\Gamma}}_{1}\|_{0,\infty}^{\scriptscriptstyle{\SS}}<\frac{1}{2}\left(1+\frac{1}{\mu({\mathbf{D}}_{1})}\right)$. Moreover, since the set $\{\hat{{\bm{\Gamma}}}_{i}\}_{i=1}^{K}$ is a solution of the ${\text{DCP}_{\bm{\lambda}}}$ problem, we also have that $\|\hat{{\bm{\Gamma}}}_{1}\|_{0,\infty}^{\scriptscriptstyle{\SS}}\leq\lambda_{1}<\frac{1}{2}\left(1+\frac{1}{\mu({\mathbf{D}}_{1})}\right)$. As such, in light of the aforementioned uniqueness theorem, both representations are equal. Once we have concluded that ${\bm{\Gamma}}_{1}=\hat{{\bm{\Gamma}}}_{1}$, we would also like to show that the representations ${\bm{\Gamma}}_{2}$ and $\hat{{\bm{\Gamma}}}_{2}$ are identical. Similarly, the assumptions $\|{\bm{\Gamma}}_{2}\|_{0,\infty}^{\scriptscriptstyle{\SS}}<\frac{1}{2}\left(1+\frac{1}{\mu({\mathbf{D}}_{i})}\right)$ and $\|\hat{{\bm{\Gamma}}}_{2}\|_{0,\infty}^{\scriptscriptstyle{\SS}}\leq\lambda_{2}<\frac{1}{2}\left(1+\frac{1}{\mu({\mathbf{D}}_{2})}\right)$ guarantee that ${\bm{\Gamma}}_{2}=\hat{{\bm{\Gamma}}}_{2}$. The same set of steps can be applied for all $1\leq i\leq K$, leading to the fact that both sets of representations are identical.

Proof In (Papyan et al., 2016b), for a signal ${\mathbf{Y}}={\mathbf{X}}+{\mathbf{E}}={\mathbf{D}}_{1}{\bm{\Gamma}}_{1}+{\mathbf{E}}$, it was shown that if the following hold:

$\|{\bm{\Gamma}}_{1}\|_{0,\infty}^{\scriptscriptstyle{\SS}}<\frac{1}{2}\left(1+\frac{1}{\mu({\mathbf{D}}_{1})}\right)$ and $\|{\mathbf{E}}\|_{2}=\|{\mathbf{Y}}-{\mathbf{D}}_{1}{\bm{\Gamma}}_{1}\|_{2}\leq\mathcal{E}_{0}$,

$\|\hat{{\bm{\Gamma}}}_{1}\|_{0,\infty}^{\scriptscriptstyle{\SS}}<\frac{1}{2}\left(1+\frac{1}{\mu({\mathbf{D}}_{1})}\right)$ and $\|{\mathbf{Y}}-{\mathbf{D}}_{1}\hat{{\bm{\Gamma}}}_{1}\|_{2}\leq\mathcal{E}_{0}$,

In the above, we have defined ${\bm{\Delta}}_{1}$ as the difference between the true sparse vector, ${\bm{\Gamma}}_{1}$, and the corresponding representation obtained by solving the ${\text{DCP}_{\bm{\lambda}}^{\hskip 1.13791pt{\bm{\mathbf{\mathcal{E}}}}}}$ problem, $\hat{{\bm{\Gamma}}}_{1}$. In item 2 we have used the fact that the solution for the ${\text{DCP}_{\bm{\lambda}}^{\hskip 1.13791pt{\bm{\mathbf{\mathcal{E}}}}}}$ problem, $\hat{{\bm{\Gamma}}}_{1}$, must satisfy $\|{\mathbf{Y}}-{\mathbf{D}}_{1}\hat{{\bm{\Gamma}}}_{1}\|_{2}\leq\mathcal{E}_{0}$ and $\|\hat{{\bm{\Gamma}}}_{1}\|_{0,\infty}^{\scriptscriptstyle{\SS}}\leq\lambda_{1}$; and our assumption that $\lambda_{1}<\frac{1}{2}\left(1+\frac{1}{\mu({\mathbf{D}}_{1})}\right)$. Next, notice that $\hat{{\bm{\Gamma}}}_{1}={\bm{\Gamma}}_{1}+{\bm{\Delta}}_{1}={\mathbf{D}}_{2}{\bm{\Gamma}}_{2}+{\bm{\Delta}}_{1}$, and that the following hold:

$\|{\bm{\Gamma}}_{2}\|_{0,\infty}^{\scriptscriptstyle{\SS}}<\frac{1}{2}\left(1+\frac{1}{\mu({\mathbf{D}}_{2})}\right)$ and $\|{\bm{\Delta}}_{1}\|_{2}=\|\hat{{\bm{\Gamma}}}_{1}-{\mathbf{D}}_{2}{\bm{\Gamma}}_{2}\|_{2}\leq\mathcal{E}_{1}$,

$\|\hat{{\bm{\Gamma}}}_{2}\|_{0,\infty}^{\scriptscriptstyle{\SS}}<\frac{1}{2}\left(1+\frac{1}{\mu({\mathbf{D}}_{2})}\right)$ and $\|\hat{{\bm{\Gamma}}}_{1}-{\mathbf{D}}_{2}\hat{{\bm{\Gamma}}}_{2}\|_{2}\leq\mathcal{E}_{1}$.

The second item relies on the fact that both $\hat{{\bm{\Gamma}}}_{1}$ and $\hat{{\bm{\Gamma}}}_{2}$, obtained by solving the ${\text{DCP}_{\bm{\lambda}}^{\hskip 1.13791pt{\bm{\mathbf{\mathcal{E}}}}}}$ problem, must satisfy $\|\hat{{\bm{\Gamma}}}_{1}-{\mathbf{D}}_{2}\hat{{\bm{\Gamma}}}_{2}\|_{2}\leq\mathcal{E}_{1}$ and $\|\hat{{\bm{\Gamma}}}_{2}\|_{0,\infty}^{\scriptscriptstyle{\SS}}\leq\lambda_{2}$. In addition, the second expression uses the assumption that $\lambda_{2}<\frac{1}{2}\left(1+\frac{1}{\mu({\mathbf{D}}_{2})}\right)$. Employing once again the aforementioned stability theorem, we are guaranteed that

Using the same set of steps presented above, we conclude that

C Stable Recovery of Hard Thresholding in the Presence of Noise (Proof of Lemma 7)

Proof Denote by $\mathcal{T}_{1}$ the support of ${\bm{\Gamma}}_{1}$. Denote further the $i$-th atom from ${\mathbf{D}}_{1}$ by ${\mathbf{d}}_{1,i}$. The success of the hard thresholding algorithm with threshold $\beta_{1}$ in recovering the correct support is guaranteed if the following holds

Using the same set of steps as those used in proving Theorem 4 in (Papyan et al., 2016b), we can lower bound the left-hand-side by

we ensure the success of the thresholding algorithm. This condition can be equally written as

Equation (82) also implies that the threshold $\beta_{1}$ that should be employed must satisfy

Thus far, we have considered the successful recovery of the support of ${\bm{\Gamma}}_{1}$. Next, assuming this correct support was recovered, we shall dwell on the deviation of the thresholding result, $\hat{{\bm{\Gamma}}}_{1}$, from the true ${\bm{\Gamma}}_{1}$. Denote by ${\bm{\Gamma}}_{1,\mathcal{T}_{1}}$ and $\hat{{\bm{\Gamma}}}_{1,\mathcal{T}_{1}}$ the vectors ${\bm{\Gamma}}_{1}$ and $\hat{{\bm{\Gamma}}}_{1}$ restricted to the support $\mathcal{T}_{1}$, respectively. We have that

In what follows, we shall upper bound both of the expressions in the right hand side of the inequality.

In the remainder of this proof we will localize the above bound into one that is posed in terms of patch-errors. Note that $\|{\bm{\Gamma}}_{1}-\hat{{\bm{\Gamma}}}_{1}\|_{2,\infty}^{\scriptscriptstyle{{\mathbf{P}}}}$ is equal to the maximal energy of an $n_{1}m_{1}$-dimensional patch taken from it, where the $i$-th patch can be extracted using the operator ${\mathbf{P}}_{1,i}$. Relying on this and the relation $\|{\mathbf{V}}\|_{2}\leq\sqrt{\|{\mathbf{V}}\|_{0}}\ \|{\mathbf{V}}\|_{\infty}$, we have that

Recalling that, based on Definition 6, $\|{\bm{\Gamma}}_{1}-\hat{{\bm{\Gamma}}}_{1}\|_{0,\infty}^{\scriptscriptstyle{{\mathbf{P}}}}$ denotes the maximal number of non-zeros in a patch of length $n_{1}m_{1}$ extracted from this vector, we obtain that

In the last inequality we have used the success of the first stage in recovering the correct support, resulting in $\|{\bm{\Gamma}}_{1}-\hat{{\bm{\Gamma}}}_{1}\|_{0,\infty}^{\scriptscriptstyle{{\mathbf{P}}}}\leq\|{\bm{\Gamma}}_{1}\|_{0,\infty}^{\scriptscriptstyle{{\mathbf{P}}}}$. Plugging inequality (98) into the above equation, we conclude that

D Stability of the Layered Hard Thresholding in the Presence of Noise (Proof of Theorem 8)

Proof The stability of the first stage of the layered hard thresholding algorithm is obtained from Lemma 7. Denoting by ${\bm{\Delta}}_{1}=\hat{{\bm{\Gamma}}}_{1}-{\bm{\Gamma}}_{1}$, notice that $\hat{{\bm{\Gamma}}}_{1}={\bm{\Gamma}}_{1}+{\bm{\Delta}}_{1}={\mathbf{D}}_{2}{\bm{\Gamma}}_{2}+{\bm{\Delta}}_{1}$. In other words, ${\bm{\Gamma}}_{1}$ is a signal that admits a convolutional sparse representation ${\mathbf{D}}_{2}{\bm{\Gamma}}_{2}$, which is contaminated with noise ${\bm{\Delta}}_{1}$, resulting in $\hat{{\bm{\Gamma}}}_{1}$. Next, we would like to employ Lemma 7 for the signal $\hat{{\bm{\Gamma}}}_{1}={\bm{\Gamma}}_{1}+{\bm{\Delta}}_{1}={\mathbf{D}}_{2}{\bm{\Gamma}}_{2}+{\bm{\Delta}}_{1}$, with the local noise level

Assuming the above hold, Lemma 7 guarantees that the support of $\hat{{\bm{\Gamma}}}_{2}$ is equal to that of ${\bm{\Gamma}}_{2}$, and also that

Using the same steps as above, we obtain the desired claim for all the remaining layers, assuming that

and that the thresholds $\beta_{i}$ are chosen to satisfy

E Stable Recovery of Soft Thresholding in the Presence of Noise (Proof of Lemma 9)

Proof The success of the soft thresholding algorithm with threshold $\beta_{1}$ in recovering the correct support is guaranteed if the following holds

Since the soft thresholding operator chooses all atoms with correlations greater than $\beta_{1}$, the above implies that the true support $\mathcal{T}_{1}$ will be chosen. This condition is equal to that of the hard thresholding algorithm, and thus using the same steps as in Lemma 7, we are guaranteed that the correct support will be chosen under Assumptions (a) and (b).

The difference between the hard thresholding algorithm and its soft counterpart becomes apparent once we consider the estimated sparse vector. While the former estimates the non-zero entries in $\hat{{\bm{\Gamma}}}_{1}$ by computing ${\mathbf{D}}_{1,\mathcal{T}_{1}}^{T}{\mathbf{Y}}$, the latter subtracts or adds a constant $\beta_{1}$ from these, obtaining ${\mathbf{D}}_{1,\mathcal{T}_{1}}^{T}{\mathbf{Y}}-\beta_{1}{\mathbf{B}}$, where ${\mathbf{B}}$ is a vector of $\pm 1$. As a result, the distance between the true sparse vector and the estimated one is given by

F Stability of the Layered BP Algorithm in the Presence of Noise (Proof of Theorem 12)

Proof In (Papyan et al., 2016b), for a signal ${\mathbf{Y}}={\mathbf{X}}+{\mathbf{E}}={\mathbf{D}}_{1}{\bm{\Gamma}}_{1}+{\mathbf{E}}$, it was shown that if the following hold:

$\|{\mathbf{Y}}-{\mathbf{X}}\|_{2,\infty}^{\scriptscriptstyle{{\mathbf{P}}}}\leq\epsilon_{0}$; and

$\|{\bm{\Gamma}}_{1}\|_{0,\infty}^{\scriptscriptstyle{\SS}}<\frac{1}{3}\left(1+\frac{1}{\mu({\mathbf{D}}_{1})}\right)$,

then the solution $\hat{{\bm{\Gamma}}}_{1}$ for the Lagrangian formulation of the BP problem with $\beta_{1}=4\epsilon_{0}$ (see Equation (54)) satisfies that

The support of the solution $\hat{{\bm{\Gamma}}}_{1}$ is contained in that of ${\bm{\Gamma}}_{1}$;

$\|{\bm{\Delta}}_{1}\|_{\infty}=\|{\bm{\Gamma}}_{1}-\hat{{\bm{\Gamma}}}_{1}\|_{\infty}\leq 7.5\ \epsilon_{0}$;

In particular, every entry of ${\bm{\Gamma}}_{1}$ greater in absolute value than $7.5\ \epsilon_{0}$ is guaranteed to be recovered; and

The solution $\hat{{\bm{\Gamma}}}_{1}$ is the unique minimizer of the Lagrangian BP problem (Equation (54)).

Using similar steps to those employed in the proof of Theorem 8, we obtain that

Plugging above the inequality $\|{\bm{\Delta}}_{1}\|_{\infty}\leq 7.5\ \epsilon_{0}$, we get

Since the support of $\hat{{\bm{\Gamma}}}_{1}$ is contained in that of ${\bm{\Gamma}}_{1}$, we have that $\|{\bm{\Delta}}_{1}\|_{0,\infty}^{\scriptscriptstyle{{\mathbf{P}}}}=\|{\bm{\Gamma}}_{1}-\hat{{\bm{\Gamma}}}_{1}\|_{0,\infty}^{\scriptscriptstyle{{\mathbf{P}}}}\leq\|{\bm{\Gamma}}_{1}\|_{0,\infty}^{\scriptscriptstyle{{\mathbf{P}}}}$, leading to

We conclude that the first stage of the layered BP is stable and the following must hold

The support of the solution $\hat{{\bm{\Gamma}}}_{1}$ is contained in that of ${\bm{\Gamma}}_{1}$;

$\|{\bm{\Delta}}_{1}\|_{2,\infty}^{\scriptscriptstyle{{\mathbf{P}}}}\leq\ \epsilon_{1}$;

In particular, every entry of ${\bm{\Gamma}}_{1}$ greater in absolute value than $\frac{\epsilon_{1}}{\sqrt{\|{\bm{\Gamma}}_{1}\|_{0,\infty}^{\scriptscriptstyle{{\mathbf{P}}}}}}$ is guaranteed to be recovered; and

The solution $\hat{{\bm{\Gamma}}}_{1}$ is the unique minimizer of the Lagrangian BP problem (Equation (54)).

Next, we turn to the stability of the second stage of the layered BP algorithm. Notice that $\hat{{\bm{\Gamma}}}_{1}={\bm{\Gamma}}_{1}+{\bm{\Delta}}_{1}={\mathbf{D}}_{2}{\bm{\Gamma}}_{2}+{\bm{\Delta}}_{1}$. Put differently, ${\bm{\Gamma}}_{1}$ is a signal that admits a convolutional sparse representation ${\mathbf{D}}_{2}{\bm{\Gamma}}_{2}$ that is perturbed by ${\bm{\Delta}}_{1}$, resulting in $\hat{{\bm{\Gamma}}}_{1}$. As such, we can invoke once again the same theorem from (Papyan et al., 2016b). Since we have that

$\|{\bm{\Delta}}_{1}\|_{2,\infty}^{\scriptscriptstyle{{\mathbf{P}}}}\leq\epsilon_{1}$; and

$\|{\bm{\Gamma}}_{2}\|_{0,\infty}^{\scriptscriptstyle{\SS}}<\frac{1}{3}\left(1+\frac{1}{\mu({\mathbf{D}}_{2})}\right)$,

we are guaranteed that the solution $\hat{{\bm{\Gamma}}}_{2}$ for the Lagrangian formulation of the BP problem with parameter $\beta_{2}=4\epsilon_{1}$ satisfies

The support of the solution $\hat{{\bm{\Gamma}}}_{2}$ is contained in that of ${\bm{\Gamma}}_{2}$;

$\|{\bm{\Delta}}_{2}\|_{\infty}=\|{\bm{\Gamma}}_{2}-\hat{{\bm{\Gamma}}}_{2}\|_{\infty}\leq 7.5\ \epsilon_{1}$ ;

In particular, every entry of ${\bm{\Gamma}}_{2}$ greater in absolute value than $7.5\ \epsilon_{1}$ is guaranteed to be recovered; and

The solution $\hat{{\bm{\Gamma}}}_{2}$ is the unique minimizer of the Lagrangian BP problem (Equation (54)).

We conclude that, similar to the first one, the second stage of the layered BP is stable and the following must hold

The support of the solution $\hat{{\bm{\Gamma}}}_{2}$ is contained in that of ${\bm{\Gamma}}_{2}$;

$\|{\bm{\Delta}}_{2}\|_{2,\infty}^{\scriptscriptstyle{{\mathbf{P}}}}\leq\ \epsilon_{2}$;

In particular, every entry of ${\bm{\Gamma}}_{2}$ greater in absolute value than $\frac{\epsilon_{2}}{\sqrt{\|{\bm{\Gamma}}_{2}\|_{0,\infty}^{\scriptscriptstyle{{\mathbf{P}}}}}}$ is guaranteed to be recovered; and

The solution $\hat{{\bm{\Gamma}}}_{2}$ is the unique minimizer of the Lagrangian BP problem (Equation (54)).

Using the same set of steps, we obtain similarly the stability of the remaining layers.

References