Title: TOWARDS UNDERSTANDING THE INVERTIBILITY OF CONVOLUTIONAL NEURAL NETWORKS

Abstract: Several recent works have empirically observed that Convolutional Neural Nets (CNNs) are (approximately) invertible. To understand this approximate invertibility phenomenon and how to leverage it more effectively, we focus on a theoretical explanation and develop a mathematical model of sparse signal recovery that is consistent with CNNs with random weights. We give an exact connection to a particular model of model-based compressive sensing (and its recovery algorithms) and random-weight CNNs. We show empirically that several learned networks are consistent with our mathematical analysis and then demonstrate that with such a simple theoretical framework, we can obtain reasonable reconstruction results on real images. We also discuss gaps between our model assumptions and the CNN trained for classification in practical scenarios.

1. INTRODUCTION:

Deep learning has achieved remarkable success in many technological areas (Bengio et al., 2013; Schmidhuber, 2015) , including computer vision (Krizhevsky et al., 2012; Szegedy et al., 2015; Simonyan and Zisserman, 2015) , automatic speech recognition (Hinton et al., 2012; Hannun et al., 2014) , natural language processing (Collobert et al., 2011; Mikolov et al., 2013; Cho et al., 2014) , bioinformatics (Chicco et al., 2014) , even high energy particle physics (Baldi et al., 2014) . In particular, deep Convolutional Neural Networks (CNNs) (LeCun et al., 1989; Krizhevsky et al., 2012; Simonyan and Zisserman, 2015) have been a critical enabling technique for analyzing images and sequential data.

Following the unprecedented success of deep networks, there has been some theoretical work (e.g., Arora et al. (2014; 2015) ; Paul and Venkatasubramanian (2014) ) that suggest several mathematical models for different deep learning architectures. However, theoretical analysis and understanding lag behind the very rapid evolution and empirical success of deep architectures, and more theoretical analysis is needed to better understand the state-of-the-art deep architectures, and possibly to improve them further.

In this paper, we attempt to address the gap between the empirical success and theoretical understanding of the Convolutional Neural Nets, in particular its invertibility (i.e., reconstructing the input from the hidden activations), by analyzing a simplified mathematical model using random weights. 1This property is intriguing because convolutional neural networks are typically trained with discriminative objectives (i.e., unrelated to reconstruction) with a large amount of labels, such as the ImageNet dataset. For example, Dosovitskiy and Brox (2016) used upsampling-deconvolutional architectures to invert the hidden activations of feedforward CNNs to the input domain. In other related work, Zhao et al. (2016) proposed stacked a what-where network via a (deconvolutional) decoder and demonstrate its promise in unsupervised and semi-supervised settings. Bruna et al. (2014) studied signal discovery from generalized pooling operators using image patches on non-convolutional small scale networks and datasets. Zhang et al. (2016) showed that CNNs discriminately trained for image classification (e.g., VGG Net (Simonyan and Zisserman, 2015) ) are almost fully invertible using pooling switches. Despite these interesting results, there is no clear theoretical explanation as to why CNNs are invertible yet.

We introduce three new concepts that, coupled with the accepted notion that images have sparse representations, guide our understanding of CNNs:

1. we provide a particular model of sparse linear combinations of the learned filters that are consistent with natural images; also, this model of sparsity is itself consistent with the feedforward network;

2. we show that the effective matrices that capture explicitly the convolution of multiple filters exhibit a model-Restricted Isometry Property (model-RIP) (Baraniuk et al., 2010) ; and 3. our model can explain each layer of the feedforward CNN algorithm as one iteration of Iterative Hard Thresholding (IHT) (Blumensath and Davies, 2009) for model-based compressive sensing and, hence, we can reconstruct the input simply and accurately.

In other words, we give a theoretical connection to a particular model of model-based compressive sensing (and its recovery algorithms) and CNNs. We show empirically that large-scale deep convolution networks are consistent with our mathematical analysis. We then demonstrate that with such a simple theoretical framework, we can obtain reasonable reconstruction results on real images, using filters from trained networks. Finally, we observe that it makes a significant difference which filters one uses for encoding and decoding, whether they are trained specifically for reconstruction, or random, or the same for both procedures. This paper explores these properties and elucidate specific empirical aspects that any more sophisticated mathematical model should take into account.2 

2. PRELIMINARIES:

In this section, we set the stage for our mathematical analysis in Section 3. We begin with discussion on the use of random weights in (convolutional) neural networks, and then provide the definitions and models for CNNs. Then, we discuss compressive sensing and sparse signal recovery. We define a particular model of sparsity that we will use throughout our analysis and detail the Iterative Hard Thresholding (IHT) algorithm which is the basis of our reconstruction analysis.

In order to simplify our notation and to make clear our analysis, we focus on a single layer in the analysis instead of multiple layers.3 Also, we assume that all of our input signals are vectors rather than matrices and that any operations we would ordinarily carry out on images (e.g., convolving with a filter bank, dividing into regions over which we pool coefficients), we do on vectors with the appropriate modifications for a simplified structure. While these assumptions ease our exposition, they do not change the nature of our arguments nor their implications for images. Furthermore, we demonstrate the validity of our results in two-dimensional natural images.

2.1. EFFECTIVENESS OF GAUSSIAN RANDOM FILTERS:

We analyze theoretically CNNs with Gaussian random filters, which have been surprisingly effective in unsupervised and supervised deep learning tasks. Jarrett et al. (2009) showed that random filters in 2-layer CNNs work well for image classification. In addition, Saxe et al. (2011) observed that convolutional layer followed by pooling layer is frequency selective and translation invariant, even with random filters, and these properties lead to good performance for object recognition tasks.

On the other hand, Giryes et al. (2016) proved that CNNs with random Gaussian filters have metric preservation property, and they argued that the role of training is to select better hyperplanes discriminating classes by distorting boundary points among classes. According to their observation, random filters are in fact a good choice if training data are initially well-separated. Also, He et al. (2016) empirically showed that random weight CNNs can do image reconstruction well.

To better demonstrate the effectiveness of Gaussian random CNNs, we evaluate their classification performance on CIFAR-10; see Section 4.1 for details. We find that a 3-layer Gaussian random CNN is able to achieve ∼ 75% accuracy on the test set, with only the last classifier layer optimized, (see Table 1 for more details). Even though this number is far from the state-of-the-art results, it is surprisingly good considering the networks are almost untrained. Our theoretical results may provide another new perspective on explaining these phenomena.

n MD x h W = j=0,…,n-1 = Kn D (w 1,1 j ) T t l (w 1,M j ) T (w K,1 j ) T (w K,M j ) T

Figure 1 : One-dimensional CNN architecture where W ∈ R Kn×M D is the matrix instantiation of convolution over M channels with a filter bank consisting of K different filters. Note that a filter bank has K filters of size l × M , such that there are lM K parameters in this architecture.

2.2. CONVOLUTIONAL NEURAL NETS:

We define a single layer of our CNN as follows. We assume that the input signal x consists of M channels, each of length D, and we write x ∈ R M D . For each of the input channels, m = 1, . . . , M , let w i,m , i = 1, . . . , K denote one of K filters, each of length . Let t be the stride length, the number of indices by which we shift each filter. Note that t can be larger than 1.

We assume that the number of shifts, n = (D -)/t + 1, is an integer. Let w j i,m be a vector of length D that consists of the (i, m)-th filter shifted by jt, j = 0, . . . , n -1 (i.e., w j i,m has at most non-zero entries). We will concatenate over the M channels each of these vectors (as row vectors) to form a large matrix, W , which is the Kn × M D matrix made up of K blocks of the n shifts of each filter in each of M channels. We assume that Kn ≥ M D. We also assume that the Kn row vectors of W span R M D and that we have normalized the rows so that they have unit 2 norm. We assume that the hidden units of the feed-forward CNN are computed by multiplying an input signal x ∈ R M D by the matrix W (i.e., convolving, in each channel, by a filter bank of size K, and summing over the channels to obtain Kn outputs), applying the ReLU function to the Kn outputs, and then selecting the value with maximum absolute value in each of the K blocks; i.e., we perform max pooling over each of the convolved filters and sum over the channels. 4 We use h = W x for the hidden activation computed by a single layer CNN without pooling. Figure 1 illustrates the architecture.

2.3. COMPRESSIVE SENSING:

Let Φ be a i×j matrix with j > i. We say that Φ satisfies the Restricted Isometry Property RIP(k, δ k ) (or, just RIP) if there is a distortion factor δ k > 0 such that for all z ∈ R j with exactly k non-zero entries, (1

-δ k ) z 2 2 ≤ Φz 2 2 ≤ (1 + δ k ) z 2 2 .

If Φ satisfies RIP (for appropriate sparsity level k and sufficiently small δ k ) and if z ∈ R j is k-sparse, then, given the vector x = Φz ∈ R i , we can efficiently recover z (see Candés (2008) for more details)5 . There are many efficient algorithms for doing so, including 1 sparse coding (e.g., 2 minimization with 1 regularization) and greedy, iterative algorithms (such as Iterative Hard Thresholding or IHT).

Model-based compressive sensing. While sparse signals are a natural model for some applications, they are less realistic for CNNs. We consider a vector z ∈ R Kn as the true sparse code for generating the CNN input x with a particular model of sparsity. Rather than permitting k non-zero entries anywhere in the vector z, we divide the support of z into K contiguous blocks of size n and we stipulate that from each block there is at most one non-zero entry in z with a total of k non-zero entries. We call a vector with this sparsity model model-k-sparse and denote the union of all ksparse subspaces with this structure M k . It is clear that M k contains n k K k subspaces. In our analysis, we consider linear combinations of two model-k-sparse signals. To be precise, suppose that z = α 1 z 1 + α 2 z 2 is the linear combination of two elements in M k . Then, we say that z lies in the linear subspace M 2 k that consists of all linear combinations of vectors from M k . We say that a matrix Φ satisfies the model-RIP condition for parameter k if, there is a distortion factor δ k > 0 such that, for all z ∈ M k ,

(1 -δ k ) z 2 2 ≤ Φz 2 2 ≤ (1 + δ k ) z 2 2 .

(1)

See Baraniuk et al. (2010) for the definitions of model sparse and model-RIP, as well as the necessary modifications to account for signal noise and compressible (as opposed to exactly sparse) signals (which we have neglected to consider to keep our analysis simple). Intuitively speaking, a matrix that satisfies the model-RIP is a nearly an orthonormal matrix for a particular set of sparse vectors with a particular sparsity model or pattern.

For our analysis, we also need matrices Φ that satisfy the model-RIP condition for vectors z ∈ M 2 k . We denote the distortion factor δ 2k for such matrices. Note that δ k ≤ δ 2k < 1.

Algorithm 1 Model-based IHT:

Input: model-RIP matrix Φ, measurements x = Φz, structured sparse approximation algorithm M Output: k-sparse approximation z 1: Initialize z 0 = 0, d = x, i = 0 2: while stopping criteria not met do

3: i ← i + 1 4: b ← z i-1 + Φ T d 5: z i ← M(b, k) 6:

d ← x -Φz i 7: end while 8: return z ← z i Many efficient algorithms have been proposed for sparse coding and compressive sensing (Olshausen et al., 1996; Mallat and Zhang, 1993; Beck and Teboulle, 2009) . As with traditional compressive sensing, there are efficient algorithms for recovering model-k-sparse signals from measurements (see Baraniuk et al. (2010) ), assuming the existence of an efficient structured sparse approximation algorithm M, that given an input vector and the sparsity parameter, returns the vector closest to the input with the specified sparsity structure.

In convolutional neural networks, the max pooling operator finds the downsampled activations that are closest to the activations of the original size by retaining the most significant values. The max pooling can be viewed as two steps: 1) zeroing out the locally non-maximum values; 2) downsampling the activations with the locally maximum values retained. To study the pooled activations with sparsity structures, we can recover dimension loss from the second step (downsampling step) by an unsampling operator. This procedure defines our structured sparse approximation algorithm M(z, k), where z is the original (unpooled) code, and k is the sparsity parameter for further sparsification, which guarantees that M(z, k) is a model-k-sparse signal. With the standard layered formulation for neural networks, we have

EQUATION

where s denotes the upsampling switches that indicate where to place the non-zero values in the upsampled activations. Taking the pooling switches known from the max pooling operation as s, we specifically define M as the nesting of the max pooling and the unpooling with known switch. We define this special case as

EQUATION

Alternatively, using the fixed uniform switches as s, we specifically define M as the nesting of the max pooling and the naive unsampling, denoted by M fixed . In the rest of this paper, our theoretical analysis are generic to any type of valid upsampling switches6 , so we use M(z, k) to denote the structured sparse approximation algorithm without worrying about s. The two special cases M known and M fixed are used in the empirical analysis when we need to specify M(z, k) as a fully concrete operator.

The main recovery algorithm that we focus on is a model-sparse version of Iterative Hard Thresholding (IHT) (see Blumensath and Davies (2009) ), not because we are interested in recovering modelsparse signals, per se, but because one iteration of IHT for our model of sparsity captures exactly a feedforward CNN. 7 Algorithm 1 describes the model-based IHT algorithm. In particular, the sequence of steps 4-6 in the middle IHT (without the outer iterative loop) is exactly one layer of a feedforward CNN. As a result, the theoretical analysis of IHT for model-based sparse signal recovery serves as a guide for how to analyze the approximation activations of a CNN.

3. ANALYSIS:

To motivate our more formal analysis, we begin with a simple example. Suppose that the matrix W is an orthonormal basis for R M D and define

Ψ = W T -W T .

Proposition 1. A one-layer CNN using the matrix Ψ T , with no pooling, gives perfect reconstruction (with the matrix Ψ) for any input vector x ∈ R M D .

Proof. Because we have both the positive and the negative dot products of the signal with the basis

vectors in ReLU(Ψ T x) = ReLU W x -W x

, we have positive and negative versions of the hidden

units h + = ReLU(W x) and h -= ReLU(-W x)

where we decompose h = W x = h + -h - into the difference of two non-negative vectors, the positive and the negative entries of h. From this decomposition, we can easily reconstruct the original signal via

Ψ h + h - = W T -W T h + h - = W T (h + -h -) = W T h = W T W x = x.

In the example above, we have pairs of vectors (w, -w) in our matrix Ψ. This settings allow us to turn what would ordinarily be a nonlinear function, ReLU, into a linear one. In fact, the assumption that trained CNN filters come in positive and negative is validated by Shang et al. (2016) , which makes a CNN much easier to analyze within the model compressed sensing framework.

Suppose that we have a vector z that we split into positive and negative components, z = z + -z -, and that we synthesize (or construct) a signal x from z using the matrix W T -W T . Then, we have

W T -W T z + z - = W T (z + -z -) = W T z = x.

Next, suppose that we multiply x = W T z by the transpose of the same matrix, we find

W -W x = W W T z -W W T z

and, if we apply ReLU to this vector, we produce

(W W T z) + (W W T z) - a vector that is split

into its positive and negative components. To determine whether or not we have "reconstructed" the vector z, the structure of the product W W T is crucial. In addition, this calculation shows that if we have both positive and negative pairs of filters or vectors, then the ReLU function applied to both the positive and negative dot products simply splits the vector into the positive and negative components. These components are then reassembled in the next computation. For this reason, in the analysis in the following sections, it is sufficient to consider W T z = x and W x = h with max pooling alone applied to h, assuming that all of the entries in the vectors are real numbers, rather than only non-negative.

3.1. MODEL-RIP AND RANDOM FILTERS:

Our first main result says that if we use Gaussian random filters in our CNN, then, with high probability, the transpose of the matrix W formed by the convolutions with these filters has the model-RIP property. In other words, Gaussian random filters generate a matrix whose transpose W T is almost an orthonormal transform for sparse signals with a particular sparsity pattern (that is consistent with our pooling procedure). The bounds in the theorem tell us that we must balance the size of the filters and the number of channels M against the sparsity of the hidden units k, the number of the filter banks K, the number of shifts n, the distortion parameter δ k , and the failure probability . The proof is in Appendix A.

Theorem 3.1. Assume that we have M K vectors w i,m of length in which each entry is a scaled i.i.d. (sub-)Gaussian random variable with mean zero and variance 1 (the scaling factor is 1/ √ M ). Let t be the stride length (where n = (D -)/t + 1) and build the structured random matrix W as the weight matrix in a single layer CNN for M -channel input dimension D. If

M 2 D ≥ Ck δ 2 k log(K) + log(n) -log( ) ,

then, with probability 1 -, the M D × Kn matrix W T satisfies the model-RIP for model M k with parameter δ k .

We also note that the same analysis can be applied to the sum of two model-k-sparse signals, with changes in the constants (that we do not track here). Corollary 3.2. Random matrices with the CNN structure have, with high probability, the model-RIP property for

M 2 k .

Other examples of matrices that satisfy model-RIP (both empirically and via a less sophisticated analysis on the dot products between any two columns) include wavelets and localized Fourier bases; both examples that can be easily and efficiently implemented via convolutions in a CNN.

3.2. RECONSTRUCTION BOUNDS:

To distinguish the true sparse code z and its reconstruction, we use ẑ = M(h, k) = M(W x, k) for the reconstruction by CNN. Our next result tells us that if we compute the hidden units h from an input signal x using a weight matrix W whose transpose has the model-RIP and using max pooling over each filter ( ẑ), then we can reconstruct (approximately) the input signal x simply by multiplying the hidden units by W . This result bounds the relative error between the approximate reconstruction x and the input as a function of the distortion for the model-RIP. In our analysis, we assume that the input signal x = W T z is a sparse linear combination of hidden activations, captured approximately by the filters in W . See Appendix B for the detailed proofs. Part of our analysis also shows that the hidden units ẑ are approximately the putative coefficient vector z in the sparse linear representation for the input signal. Theorem 3.3. We assume that W T satisfies the M 2 k -RIP with constant δ k ≤ δ 2k < 1. If we use W in a single layer CNN both to compute the hidden units ẑ and to reconstruct the input x from these hidden units as x so that x = W T M(W x, k), the error in our reconstruction is

x -x 2 ≤ 5δ 2k 1 -δ k √ 1 + δ 2k √ 1 -δ 2k x 2 .

Recall that the structured sparsity approximation algorithm M includes the downsampling caused by pooling and an unsampling operator. Theorem 3.3 is applicable to any type of upsampling switches, so our reconstruction bound is generic to the particular design choice on how to recover the activation size in a decoding neural network.

4. EXPERIMENTAL EVIDENCE AND ANALYSIS:

In this section, we provide experimental validation of our theoretical model and analysis. We first validate experimentally the relevance of our assumption by examining the effectiveness of random filter CNNs. We then provide an experimental validation of our theoretical analysis on the synthetic 1D case, then we provide experimental results on more realistic scenarios. In particular, we study popular deep neural networks trained for image classification on the ImageNet ILSVRC 2012 dataset (Deng et al., 2009) . We calculate empirical model-RIP bounds for W T , showing that they are consistent with theory. Our results are also consistent with a long line of research shows that it is reasonable to model real, natural images as sparse linear combinations over learned dictionaries (e.g., Boureau et al. (2008) ; Le et al. (2013) ; Lee et al. (2008) ; Olshausen et al. (1996) ; Ranzato et al. (2007) ; Yang et al. (2010) ). In addition, we verify our theoretical bounds for the reconstruction error x -W T ẑ 2 / x 2 on real images. (This is the relative 2 distance between the original image and the reconstruction.) We investigate both randomly sampled filters and empirically learned filters in these experiments. Our implementation is based on the Caffe (Jia et al., 2014) and MatConvNet (Vedaldi and Lenc, 2015) toolboxes.

4.1. EVALUATION OF GAUSSIAN RANDOM CNNS ON CIFAR-10:

To show the practical relevance of our theoretical assumptions on using random filters for CNNs as stated in Section 2.1, we evaluate simple CNNs with Gaussian random filters (with i.i.d. zeromean unit-variance entries) on the CIFAR-10 dataset. The goal of this experiment is not to achieve state-of-the-art results, but to examine practical relevance of our assumption on random filter CNNs.

Once the CNNs weights are initialized (randomly), they are fixed during the training of the classifiers. Specifically, we test random CNNs with 1, 2, and 3 convolutional layers, where we use ReLU as the activation. A 2×2 max pooling layer follows each convolutional layer to down-sample the feature map. 8 We experiment with different filter sizes (3, 5, 7) and numbers of channels (64, 128, 256, 1024, 2048) and report the classification accuracy of the best-performing architectures based on cross-validation in Table 1 . We also report the best performance using learnable filters for comparison. More details about the architectures can be found in Section C.1 of the supplementary materials. We observe the CNNs with Gaussian random filters achieve surprisingly good classification performance (implying that they serve as reasonable representation of input data), although fully learnable CNN counterparts perform better. Our experimental results are also consistent with the observations made by Jarrett et al. (2009) and Saxe et al. (2011) . Overall, these results seem to suggest that the CNNs with Gaussian random filters might be a reasonable setup which is amenable to mathematical analysis while not being too far off in terms of practical relevance. 1) and its corollary from Lemma B.1 in the supplementary materials. We also empirically show the reconstruction bound in Theorem 3.3 on synthetic vectors x = W T z (Figure 2c ). The reconstruction error is concentrated at around 0.1-0.2 and bound under 0.5. Results in Figure 2 suggests the practical validity of our theory when the model assumptions hold.

4.3. ARCHITECTURES AND DATASET:

We conduct the rest of our experimental evaluations on the 16-layer VGGNet (Model D in Simonyan and Zisserman (2015)), 9 where the computation is carried out on images; e.g., convolution with a 2-D filter bank and pooling on square regions. In contrast to the theory, the realistic network does not pool activations over all the possible shifts for each filter, but rather on non-overlapping patches. The networks are trained for the large-scale ImageNet classification task, which is important for extending to other supervised tasks in vision. The main findings on VGGNet are presented in the rest of this section; we also provide some analysis on AlexNet (Krizhevsky et al., 2012) in the supplementary materials.

8 Implementation detail: We add a batch normalization layer together with a learnable scale and bias before the activation so that we do not need to tune the scale of the filters. The filter weights of the intermediate layers in the CNNs are not trained after random initialization. On top of the network, we use an optional average pooling layer to reduce the feature map size to 4 × 4 and a dropout layer for better regularization before feeding the feature to a learnable soft-max classifier for image classification. 9 VGGNet is practically important as it is popularly used in the community and is one of the best-performing "single-pathway" networks (i.e., no skip connections). We expect that the ResNet (e.g., trained from ImagNet) can also reconstruct images from its activations well in practice. However, the ResNet architectures are too complicated to be in the scope of our theory without further nontrivial customization. 1); supposed to be concentrated at 1), (b) W W T z 2/ z 2 (model-RIP corollary from Lemma B.1 in the supplementary materials; supposed to be concentrated at 1), and (c) xx 2/ x 2 (reconstruction bound in Theorem 3.3, supposed to be small), where z is a M k sparse signal that generates the vector x and x = W T M fixed (W x, k) is the reconstruction of x, where we use the naive unsampling to recover the reduced dimension due to pooling (see Section 2.3). VGGNet contains five groups of convolution and pooling layers, each group has 2~3 convolutional layers followed by a pooling layer. We denote the j-th convolutional layer in the i-th group "conv(i, j)," and the pooling layer "pool(i)." When we say the activations/features are from i-th layer, we mean they are the output of pool(i). Our analysis is for single convolutional layers. When evaluating the i-th layer, we take the activations from the (i -1)-th layer, and investigate the filters and output of conv(i, 1).

layer c(1,1) c(1,2) p(1) c(2,1) c(2,2) p(2) c(3,1) c(3,2) c(3,3) p(3) % of

4.4. 2D MODEL-RIP:

The key to our reconstruction bound is Theorem 3.3 is the model-RIP condition for our particular model of sparsity in Equation (1). We empirically evaluate the model-RIP property, i.e., W T z / z , for real CNN filters of the pretrained VGGNet. We use two-dimensional coefficients (or hidden units) z (each block of coefficients is of size D × D), K filters of size × , and pool the coefficients over smaller pooling regions (i.e., not over all possible shifts of each filter). The following experimental evidence suggest that the sparsity model and the model-RIP property of the filters are consistent with what we conclude from the mathematical analysis on the simpler one-dimensional case.

To check the significance of the model-RIP property (i.e., how close W T z / z is to 1) in controlled settings, we first synthesize the hidden activations z with sparse uniform random variables, which fully agree with our model assumptions. The sparsity of z is constrained to the average level of the real CNN activations (refer to Table 2 ). Given the filters of a certain convolutional layer, we use the synthetic z (in equal position to this layer's output activations) to get statistics for the model-RIP property. To be consistent with succeeding experiments, we choose conv(5, 2), while other layers (c) z is recovered by Algorithm 2 from the conv(5,1) activations after applying ReLU. The learned filters admits similar model-RIP value distributions to the random filters except for a bit larger bandwidth, which means the model-RIP condition in Equation ( 1) can empirically hold even when the filters do not necessarily subject to the i.i.d Gaussian random assumption.

show similar results. Figure 3 (a) summarizes the distribution of empirical model-RIP values, which is clearly centered around 1 and satisfies Equation ( 1) with a short tail roughly bounded by δ k < 1.

For more details of the algorithm, we normalize the filters from the conv(5, 2) layer, which are × ( = 3). All K = 512 filters with M = 512 input channels are used. 10 We set D = 15 (the same as the output activations of conv(5, 2)) and use 2 × 2 pooling regions11 (commonly used in recent deep networks). We generate 1000 M k randomly sampled sparse activation (z) maps by first sampling their non-zero supports and then filling elements on the supports uniformly from [-1, 1]. The sparsity is the same as that in conv(5, 1) activations.

Algorithm 2 Sparse hidden activation recovery:

Input: convolution matrix W , input activation/image x Output: hidden code z, satisfying our model-RIP assumption with M k and reconstructing x with W 1:

z init = argmin z x -W T z 2 2 + λ z 1 2: z model = M known (z init , k) 3: z = argmin z x -W T z 2 2 + λ z 1 , s.t. z i = 0 if z model i = 0

To conduct more realistic experiments, we observe the actual conv(5, 2) activations from VGGNet are not necessarily drawn from a model-sparse uniform distribution. This motivates us to evaluate the empirical model-RIP property on the hidden activations z that reconstruct the actual input activations x from conv(5, 1) by W T z. Per theory, the x is given by a max pooling layer, so we constrain the sparsity (i.e., the size of the support set is no more than 1 in a pooling region for a single channel). We use a simple and efficient algorithm to recover z from x in Algorithm 2. The algorithm is inspired by " 1 heuristic" method that are commonly used in practice (e.g. Boyd (2015) ). As shown in Algorithm 2, we first do 1 -regularized least squares without constraining the support set. Max pooling is then applied to figure out the support set for each pooling region. In particular, we use M known , defined in (3), to zero out the locally non-maximum values without messing up the support structures. We perform 1 -regularized least squares again on the fixed support set to recover the hidden activations satisfying the model sparsity. As shown in Figures 3 (b)-(c ), the empirical model-RIP property values for visual activations x from conv(5, 1) with/without ReLU are both close to 1. The center offset to 1 is less than 0.05 and the range bound δ k is rough less then 0.05, which agrees with the theoretical bound (1) quite well.

To gain more insight, we summarize the learned filter coherence in Table 4 for all the convolutional layers in VGGNet.12 This measures the correlation or similarity between the columns of W T and is a proxy for the value of the model-RIP parameter δ k (which we can only estimate computationally). The smaller the coherence, the smaller δ k is, and the better the reconstruction. The coherence of the learned filters is not low, which is inconsistent with our theoretical assumptions. However, the model-RIP property turns out to be robust to this mismatch. It also demonstrates the strong invertibility of CNN in practice.

With model-RIP as a sufficient condition, Theorem 3.3 provides theoretical bounds for layer-wise reconstruction via x = W T M(W x, k). This operator consists of the projection and reconstruction in one IHT iteration. Without confusion, we refer to it as IHT for notational convenience. We investigate the practical reconstruction errors on Layer 1~4 activations (i.e., pool(1)~( 4)) of VGGNet.

To encode and reconstruct intermediate activations of CNNs, we employ IHT with sparsity estimated from the real CNN activations on ILSVRC-2012 validation set (see Table 2 ). We also reconstruct input images, since CNN inversion is not limited to a single layer, and images are easier to visualize than hidden activations. To implement image reconstruction, we project the reconstructed activations into the image space via a pretrained decoding network as in (Zhang et al., 2016) , which extends a similar autoencoder architecture as in (Dosovitskiy and Brox, 2016) to a stacked "what-where" autoencoder (Zhao et al., 2016) . The reconstructed activations were scaled to have the same norm as the original activations so that we can feed them into the decoding network. As an example, Figure 4 illustrates the image reconstruction results for the hidden activations of the 4-th layer, the ground truth of which is obtained by feeding natural images to the CNNs. Interestingly, the decoding network itself is powerful, since it can reconstruct the glimpse of images with Gaussian random input, as shown in Figure 4 (e). Object shapes are recovered by using the pooling switches only in the "what-where" autoencoder. This result suggests that it is important to determine which pooling units are active and then to estimate these values accurately. These steps are consistent with the steps in the inner loop of any iterative sparse signal reconstruction algorithm.

In Figure 4 (c), we take the pretrained conv(5, 1) filters for IHT. The images recovered from the IHT reconstructed 4-th layer activations are reasonable and the reconstruction quality is significantly better than the random input baseline. We also try Gaussian random filters (Figure 4 (d) ), which agree more with the model assumptions (e.g., lower coherence, see Table 4 ). The learned filters from VGGNet perform equally well visually. IHT ties the encoder and decoder weights (no filter learning for the decoder), so it does not perform as well as the decoding network trained with a huge batch of data (Figure 4 (b) ). Nevertheless, we show both theoretically and experimentally decent reconstruction bounds for these simple reconstruction methods on real CNNs. More visualization results for more layers are in the supplementary materials (Figure 5 in Section C.3).

In Table 3 , we summarize reconstruction performance for all 4 layers. With random filters, the model assumptions hold and the IHT reconstruction is the best quantitatively. IHT with real CNN filters performs comparable to the best case and much better than the baseline established by the randomly sampled activations.

Additionally, reconstruction performance of IHT is strongly related to the filter coherence, summarized in Table 4 . Lower coherence agrees more closely with the model assumptions and leads to higher reconstruction quality. Higher coherence yields worse recovery of the hidden activation (i.e., large ẑz , where ẑ is the hidden activations recovered by IHT, z is the true activation). Compared to Algorithm 2, (one-step) IHT is not so robust to high coherence.

In summary, when the assumption of i.i.d Gaussian randomness of the CNN filters holds, our theoretical reconstruction bound strictly match with the empirical observations. More importantly, we demonstrate that the bound can still reasonably hold in practice for discriminatively learned CNN layers, which is particularly true for layers with relatively lower coherence.

5. CONCLUSION:

We introduce three concepts that tie together a particular model of compressive sensing (and the associated recovery algorithms), the properties of learned filters, and the empirical observation that CNNs are (approximately) invertible. Our experiments show that filters in trained CNNs are consistent with the mathematical properties we present while the hidden units exhibit a much richer structure than mathematical analysis suggests. Perhaps simply moving towards a compressive, rather than exactly sparse, model for the hidden units will capture the sophisticated structure in these layers of a CNN or, perhaps, we need a more sophisticated model. Our experiments also demonstrate that there is considerable information captured in the switch units (or the identities of the non-zeros in the hidden units after pooling) that no mathematical model has yet expressed or explored thoroughly.

for all t ≥ 0.

Let y i denote the ith entry of the vector y = W T z. We can write

y i = Kn j=1

W T i,j z j and observe that y i is a linear combination of i.i.d. sub-Gaussian random variables (or it is identically equal to 0) and, as such, is itself a sub-Gaussian random variable with mean zero and sub-Gaussian norm y i ψ2 ≤ C/ √ M w ψ2 z 2 (see Vershynin (2010) , Lemma 5.9). The structure of the random matrix and how many non-zero entries are in row i of W do enter the more refined bound on the sub-Gaussian norm of y i ψ2 (again, see Vershynin (2010) , Lemma 5.9 for details) but we ignore such details for this estimate as they are not necessary for the next estimate.

To obtain a concentration bound for y i 2 2 , we recall from Park et al. (2011) ; Vershynin (2010) that the sum of squares of sub-Gaussian random variables tightly concentrate.

Theorem A.2. Let Y 1 , . . . , Y M D be independent sub-Gaussian random variables with sub-Gaussian norms Y i ψ2 for all i = 1, . . . , M D. Let T = max i Y i ψ2 . The for every t ≥ 0 and every a

∈ R M D , P M D i=1 a i (Y i -EY 2 i ) ≥ t ≤ 2 exp -C min Ct 2 T 2 a 2 2 , Ct T a ∞ .

We note that although some entries y i may be identically zero, depending on the sparsity pattern of z, not all entries are. Let us define ỹi = yi yi ψ 2 so that ỹi ψ2 = 1 and observe that

P y 2 2 -z 2 2 > δ z 2 2 = P M D i=1 y i 2 ψ2 ( ỹ2 i -E ỹ2 i ) > δ z 2 2 .

We apply Theorem A.2 to the sub-Gaussian random variables ỹi with the weights y i 2 ψ2 . We have

a 2 2 = M D i=1 y i 4 ψ2 ≤ CD w 4 ψ2 z 4 2 M 2 and a ∞ ≤ C w 2 ψ2 z 2 2 M .

If we set T = 1, t = δ z 2 2 , and use the above estimates for the norms of a, we have

EQUATION

Finally, we use the concentration of measure result in a crude union bound to bound the failure probability over all vectors z ∈ M k . We take n k K k ≈ (nK) k and for a desired constant failure probability. Using the smaller term in Equation ( 4), (note that δ < 1, /D < 1, and w ψ2 ≥ 1) we have

exp -C M 2 δ 2 D w 4 ψ2 exp k(log(K) + log(n)) ≤ exp(log( )) which implies M 2 D ≥ k δ 2 w 4 ψ2 log(K) + log(n) -log( ) = C k δ 2 log(K) + log(n) -log( ) .

Therefore, if design our matrix W as described and with the parameter relationship as above, the matrix W T with satisfy the model-RIP for M k and parameter δ with probability 1 -.

Let us discuss the relationship amongst the parameters in our result. First, if we have only one channel M = 1 and the filter length = D, then our bound on the number of measurements D matches those of traditional (model-based) compressive sensing; namely, D ≥ C k δ 2 (log(K) + log(n) -log( )) . If < D (i.e., the filters are much shorter than the length of the input signal as in a CNN), then we can compensate by adding more channels; i.e., the filter length needs to be larger than √ D, or, if add more channels, D/M .

B MATHEMATICAL ANALYSIS: RECONSTRUCTION BOUNDS:

The consequences of having model-RIP are two-fold. The first is that if we assume that an input image is the structured sparse linear combination of filters, x = W T z (where z ∈ M k and W T satisfies the model-RIP property), then we know an upper and lower bound on the norm of x in terms of the norm of its sparse coefficients, x 2 ≤ (1 ± δ) z 2 . Additionally,

z 2 ≤ 1 √ 1 -δ x 2 .

More importantly, when we calculate the hidden units of x,

h = ReLU(W x) = ReLU(W W T z)

we can see that the computation of h is nothing other than the first step of a reconstruction algorithm analogous to that of model-based compressed sensing. As a result, we have a bound on the error between h and z and we see that we can analyze the approximation properties of a feedfoward CNN and its linear reconstruction algorithm. In particular, we can conclude that a feedforward CNN and a linear reconstruction algorithm provide a good approximation to the original input image.

Theorem 3.3(Restated) We assume that W T satisfies the M 2 k -RIP with constant δ k ≤ δ 2k < 1. If we use W in a single layer CNN both to compute the hidden units ẑ and to reconstruct the input x from these hidden units as x so that x = W T M(W x, k), the error in our reconstruction is

x -x 2 ≤ 5δ 2k 1 -δ k √ 1 + δ 2k √ 1 -δ 2k x 2 .

Proof. To show this result, we recall the two following lemmas from Baraniuk et al. (2010) and rephrase them in the setting of a feedforward CNN.

Lemma B.1. Suppose W T has M k -RIP with constant δ k . Let Ω be a support corresponding to a subspace in M k . Then we have the following bounds:

W Ω x 2 ≤ 1 + δ k x 2 (5)

W Ω W T Ω z 2 ≤ (1 + δ k ) z 2 (6) W Ω W T Ω z 2 ≥ (1 -δ k ) z 2 (7)

Lemma B.2. Suppose that W T has M 2 k -RIP with constant δ 2k . Let Ω be a support corresponding to a subspace of M k and suppose that z ∈ M k (not necessarily supported on Ω). Then

W Ω W T z| Ω c 2 ≤ δ 2k z| Ω c 2 .

Let Π denote the support of the M k sparse vector z. Set h = W x and set ẑ to be the result of max pooling applied to the vector h, or the best fit (with respect to the 2 norm) to h in the model M k .

Let Ω denote the support set of ẑ ∈ M k . For simplicity, we assume |Π| = k = |Ω|. Lemma B.3 (Identification). The support set, Ω, of the switch units captures a significant fraction of the total energy in the coefficient vector

z z| Ω c 2 ≤ 2δ 2k 1 -δ k z 2 .

Proof. Let h Ω and h Π be the vector h restricted to the support sets Ω and Π, respectively. Since both are support sets for M k and since Ω is the best support set for h,

h -h Ω 2 ≤ h -h Π 2 ,

and, after several calculations, we have

h| Ω\Π 2 2 ≥ h| Π\Ω 2 2 .

Using Lemma B.2 and the size |(Ω \ Π) Π| ≤ 2k, we have

h Ω\Π 2 = W Ω\Π W T z 2 ≤ δ 2k z 2 .

We can bound the other side of the inequality as

h Π\Ω 2 ≥ W Π\Ω (W T z| Π\Ω ) 2 -W Π\Ω (W T z| Ω ) 2 ≥ (1 -δ k ) z| Π\Ω 2 -δ 2k z| Ω 2 .

Since the support of z is the set Π, Π \ Ω = Ω c and we can conclude that

δ 2k z 2 ≥ (1 -δ k ) z| Ω c 2 -δ 2k z| Ω 2 ,

and with some rearrangement, we have

z| Ω c 2 ≤ 2δ 2k 1 -δ k z 2 .

To set the value of ẑ on its support set Ω, we simply set ẑ = h| Ω and ẑ| Ω c = 0. Then Lemma B.4 (Estimation).

z -ẑ 2 ≤ 5δ 2k 1 -δ k z 2 Proof. First, note that I -W Ω W T Ω 2 ≤ δ k since (1 -δ k ) ≤ sup z =0 W T Ω z 2 2 z 2 2 = σ 2 max (W T Ω ) = σ max (W Ω W T Ω ) ≤ (1 + δ k ),

where σ max is the maximum singular value. Therefore,

z -ẑ 2 ≤ z| Ω c 2 + z| Ω -ẑ| Ω 2 = z| Ω c 2 + z| Ω -W Ω (W T z| Ω + W T z| Ω c ) 2 ≤ z| Ω c 2 + (I -W Ω W T Ω )z| Ω 2 + W Ω W T z| Ω c 2 ≤ z| Ω c 2 + I -W Ω W T Ω 2 z| Ω 2 + δ 2k z| Ω c 2 ≤ z| Ω c 2 + δ k z| Ω 2 + δ 2k z| Ω c 2 ≤ (1 + δ 2k ) 2δ 2k 1 -δ k + δ k z 2 ≤ 5δ 2k 1 -δ k z 2 .

Finally, if we use the autoencoder formulation to reconstruct the original image x by setting x = W T ẑ, we can estimate the reconstruction error. We note that ẑ is M k -sparse by construction and remind the reader that W T satisfies M 2 k -model-RIP with constants δ k ≤ δ 2k 1. Then, using Lemma B.4 as well as the M 2 k -sparse properties of W T ,

x -x 2 = W T (z -ẑ) 2 ≤ 1 + δ 2k z -ẑ 2 ≤ 5δ 2k 1 -δ k 1 + δ 2k z 2 ≤ 5δ 2k 1 -δ k √ 1 + δ 2k √ 1 -δ 2k x 2 .

This proves that a feedforward CNN with a linear reconstruction algorithm is an approximate autoencoder and bounds the reconstruction error of the input image in terms of the geometric properties of the filters. 

:

For more discussion about random filters, see Sections

2.1 and 4.1.

We note that our model may not be an exact replica of a real setting, but for mathematical analysis, it is a simplified but representative abstraction of practical settings. A number of works show that random weight CNNs still achieve surprisingly good classification accuracy although they may not match the state-of-the-art results; see Sections 2.1 and 4.1 for more discussion.

We can extend the equivalency on a single layer of CNNs to multiple layer CNNs simply by using the output on one layer as the input to another, still using the steps of the inner loop of IHT.

The convolution can be computed more efficiently than a straight-forward matrix multiplication, but they are mathematically equivalent.

We note that this is a sufficient condition and that there are other, less restrictive sufficient conditions, as well as more complicated necessary conditions. Furthermore, we have not given the exact, quantitative relations amongst the parameters. For simplicity, we stick with this definition.

Valid switches should place a non-zero value at exactly one location.

Multiple iterations of IHT can improve the quality of signal recovery. However, it is rather equivalent to the recurrent version of CNNs and does not fit to the scope of this work.

We do not remove any filters including those in approximate positive/negative pairs (refer to 3).

In VGGNet, no pooling layer follows conv(5, 2). Here, we just use it in this way to analyze the convolutionpooling pair targeted by our theory.

12 The coherence is defined as the maximum (in absolute value) dot product between distinct pairs of columns of the matrix W T , i.e. µ = max i =j |WiW T j |, where Wi denote the i-th row of matrix W .

The relative error in activation space of random activations (the last column) are identical (1.414) for all layers because ff / f = √ 2 on average for Gaussian random f provided f = f .

Figures and Tables:

EXPERIMENTAL VALIDATION OF THE ANALYSIS IN 1D SYNTHETIC DATA We use 1-D synthetic data to empirically show the basic validity of our theory in terms of the model-RIP condition in Equation (1) and reconstruction bound in Theorem 3.3. We plot the histograms of the empirical model-RIP values of 1D Gaussian random filters W ( scaled by 1/ √ lM ) with size l × 1 × M × K = 5 × 1 × 32 × 96 on 1D M k sparse signal z with size D = 32 and sparsity k = 10, whose non-zero elements are drawn from a uniform distribution on [-1, 1]. The histograms in Figure 2a and 2b are tightly centered around 1, suggesting that W T satisfies the model-RIP condition in Equation (

Figure 2: For 1D scaled Gaussian random filters W , we plot the histogram of ratios (a) W T z 2/ z 2 (model-RIP condition in Equation (1); supposed to be concentrated at 1), (b) W W T z 2/ z 2 (model-RIP corollary from Lemma B.1 in the supplementary materials; supposed to be concentrated at 1), and (c) xx 2/ x 2 (reconstruction bound in Theorem 3.3, supposed to be small), where z is a M k sparse signal that generates the vector x and x = W T M fixed (W x, k) is the reconstruction of x, where we use the naive unsampling to recover the reduced dimension due to pooling (see Section 2.3).

Figure3: For VGGNet's conv(5, 2) filters W , we plot the histogram of ratios W T z 2/ z 2 (the model-RIP value derived from Equation (1); supposed to be concentrated at 1) where z is a M k sparse signal. (a) z is randomly generated with the same sparsity as the conv(5, 2) activations and from a uniform distribution for the non-zero magnitude. (b) z is recovered by Algorithm 2 from the conv(5,1) activations before applying ReLU. (c) z is recovered by Algorithm 2 from the conv(5,1) activations after applying ReLU. The learned filters admits similar model-RIP value distributions to the random filters except for a bit larger bandwidth, which means the model-RIP condition in Equation (1) can empirically hold even when the filters do not necessarily subject to the i.i.d Gaussian random assumption.

Figure 4: Visualization of images reconstructed by a pretrained decoding network with VGGNet's pool(4) activation reconstructed using different methods: (a) original image, (b) output of the 5-layer decoding network with original activation, (c) output of the decoding net with reconstructed activation by IHT with learned filters, (d) output of the decoding net with reconstructed activation by IHT with Gaussian random filters, (e) output of the decoding net with Gaussian random activation.

Figure 5: Visualization of images reconstructed by a pretrained decoding network with VGGNet's pool(4) activation reconstructed using different methods: (a) original image, (b) output of the 5-layer decoding network with original activation, (c) output of the decoding net with reconstructed activation by IHT with learned filters, (d) output of the decoding net with reconstructed activation by IHT with Gaussian random filters, (e) output of the decoding net with Gaussian random activation.

Classification accuracy of CNNs with random and learnable filters on CIFAR-10. A typical layer consists of four operators: convolution, ReLU, batch normalization and max pooling. Networks with optimal filter size and numbers of output channels are used (see Section C.1 in the supplementary materials for the architecture details). The random filters, assumed in our theoretical analysis, perform reasonably well, not far off the learned filters.

Table content: <table><tr><td>Method</td><td colspan="3">1 layer 2 layers 3 layers</td></tr><tr><td colspan="2">Random filters 66.5%</td><td>74.6%</td><td>74.8%</td></tr><tr><td>Learned filters</td><td>68.1%</td><td>83.3%</td><td>89.3%</td></tr></table>

Layer-wise sparsity of VGGNet on ILSVRC-2012 validation set. "c" stands for convolutional layers while "p" represents pooling layers. CNN with random filters in Section 4.4 can be simulated with the same sparsity.

Table content: <table><tr><td>non-zeros</td><td>49.1</td><td>69.7</td><td>80.8</td><td>67.4</td><td>49.7</td><td>70.7</td><td>53.4</td><td>51.9</td><td>28.7</td><td>45.9</td></tr><tr><td>layer</td><td colspan="3">c(4,1) c(4,2) c(4,3)</td><td>p(4)</td><td colspan="3">c(5,1) c(5,2) c(5,3)</td><td>p(5)</td><td/><td/></tr><tr><td>% of non-zeros</td><td>35.6</td><td>29.6</td><td>12.6</td><td>23.1</td><td>23.9</td><td>20.6</td><td>7.3</td><td>13.1</td><td/><td/></tr></table>

Layer-wise relative reconstruction errors by different methods in activation space and image space between reconstructed and original activations. For layer i, we take its activation after pooling from that layer and reconstruct it with different methods (using learned filters from the layer above or scaled Gaussian random filters) and feed the reconstructed activation to a pretrained corresponding decoding network.13

Table content: <table><tr><td/><td colspan="3">image space relative error</td><td colspan="4">activation space relative error</td></tr><tr><td>layer</td><td colspan="2">learned random</td><td>random</td><td colspan="3">learned random</td><td>random</td></tr><tr><td/><td>filters</td><td>filters</td><td>activations</td><td>filters</td><td>filters</td><td colspan="2">activations</td></tr><tr><td>1</td><td>0.423</td><td>0.380</td><td>0.610</td><td>0.895</td><td>0.872</td><td/><td>1.414</td></tr><tr><td>2</td><td>0.692</td><td>0.438</td><td>0.864</td><td>0.961</td><td>0.926</td><td/><td>1.414</td></tr><tr><td>3</td><td>0.326</td><td>0.345</td><td>0.652</td><td>0.912</td><td>0.862</td><td/><td>1.414</td></tr><tr><td>4</td><td>0.379</td><td>0.357</td><td>0.436</td><td>1.051</td><td>0.992</td><td/><td>1.414</td></tr><tr><td>layer</td><td/><td>(1,1)</td><td>(1,2)</td><td>(2,1)</td><td>(2,2)</td><td>(3,1)</td><td>(3,2)</td><td>(3,3)</td></tr><tr><td colspan="2">coherence of learned filters</td><td colspan="7">0.9427 0.7340 0.6435 0.7465 0.5838 0.4844 0.5194</td></tr><tr><td colspan="2">coherence of random filters</td><td colspan="7">0.6701 0.1218 0.1546 0.1053 0.1099 0.0895 0.0802</td></tr><tr><td>layer</td><td/><td>(4,1)</td><td>(4,2)</td><td>(4,3)</td><td>(5,1)</td><td>(5,2)</td><td>(5,3)</td></tr><tr><td colspan="2">coherence of learned filters</td><td colspan="6">0.4596 0.4574 0.4043 0.4099 0.4099 0.4046</td></tr><tr><td colspan="2">coherence of random filters</td><td colspan="6">0.0920 0.0619 0.0617 0.0696 0.0674 0.0674</td></tr></table>

Comparison of coherence between learned filters in each convolutional layer of VGGNet and Gaussian random filters with corresponding sizes.

Table content: <table/>

Peer review (iclr2017):
Summary of the paper The paper studies the invertiblity of convolutional neural network in the random model. A reconstruction algorithm similar to IHT is proposed for layer-wise inversion of the network. Clarity: - The paper is confusing wrt to standard notations in deep learning. Comments: The paper makes two simplifications in the analysis of a CNN, that makes it map to a model based compressive sensing framework: 1- The non linearity (RELU) is dropped. This is a big simplification, for random gaussian weights for instance we know by JL that we can preserve L_2 distance, when RELU is applied the metric changes (see for instance the kernel for n=1 in http://cseweb.ucsd.edu/~saul/papers/nips09_kernel.pdf). 2- The pooling operation is modified to be shrinkage operator, that keeps the maximum value in a block and sets to the zero other values, hence the dimensionality of the pooled representation is the same of the un-pooled representations. Note that in that cases the locations of where the maximum happens is known. This is not the case in the ’standard’ max pooling definition. IHT does not map to a forward of a CNN as described in the paper. (see next point) 3- It maybe that the notations used in the paper are implying some confusions wrt to the standard notations in deep learning. - Let z be the standard pooled representation of dimension ‘k’. - Define U as the unpooling operation U(z, locations of maximum) := M(Wx,k). Hence your model of inversion is assuming the knowledge of the ’standard’ pooling representation and the switches (max locations). Referring to M as a pooling operation is misleading and confusing, it is the ‘standard’ unpooling operation. - Under this notations W^{	op} U(z, locations of maximum) is a sensible reconstruction algorithm , with the simplification of ignoring the RELU. - Under these notations, IHT is similar to a backward of the encoding neural network not a forward. It is known if you take the derivative with respect to the input in a neural network, you get the transpose convolution also known as ‘deconvolution’ (which is not a correct naming) , hence this IHT iteration is nothing else then the well known transposed convolution. 4- I suggest the authors to rewrite the paper taking into account standard definitions and notation in deep learning as discussed in point 3, and to clearly state that the recovery is done under the knowledge of the switches. ===== After reading the authors rebuttal and the revisions , the reviewer maintains his main concerns with the paper. Many improvements are still needed for the paper to be ready to be published in ICLR, at the notation , presentation and the theory levels. =====