Towards Understanding the Invertibility of Convolutional Neural Networks

TL;DR

This paper investigates the approximate invertibility of CNNs by developing a mathematical model based on sparse signal recovery, connecting it to compressive sensing, and validating it with empirical results on real images.

Contribution

It introduces a theoretical framework linking CNN invertibility to sparse signal recovery and compressive sensing, providing insights into CNN reconstruction capabilities.

Findings

Empirical networks align with the mathematical model

Reasonable image reconstructions are achievable using the model

Discussion of model limitations in practical CNNs

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 re- construction results on real images. We also discuss gaps between our model assumptions and the CNN trained for classification in practical scenarios.

Tables

Table 1. Table 1: 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 Appendix  C.1 for more details about the architectures). The random filters, assumed in our theoretical analysis, perform reasonably well, not far off the learned filters.

Method	1 layer	2 layers	3 layers
Random filters	66.5%	74.6%	74.8%
Learned filters	68.1%	83.3%	89.3%

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

layer	c(1,1)	c(1,2)	p(1)	c(2,1)	c(2,2)	p(2)
% of non-zeros	49.1	69.7	80.8	67.4	49.7	70.7
layer	c(3,1)	c(3,2)	c(3,3)	p(3)	c(4,1)	c(4,2)
% of non-zeros	53.4	51.9	28.7	45.9	35.6	29.6
layer	c(4,3)	p(4)	c(5,1)	c(5,2)	c(5,3)	p(5)
% of non-zeros	12.6	23.1	23.9	20.6	7.3	13.1

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

layer	(1,1)	(1,2)	(2,1)	(2,2)	(3,1)	(3,2)	(3,3)
learned	0.943	0.734	0.644	0.747	0.584	0.484	0.519
random	0.670	0.122	0.155	0.105	0.110	0.090	0.080
layer	(4,1)	(4,2)	(4,3)	(5,1)	(5,2)	(5,3)
learned	0.460	0.457	0.404	0.410	0.410	0.405
random	0.092	0.062	0.062	0.070	0.067	0.067

Table 4. Table 4: Layer-wise relative reconstruction errors by different methods in activation space and image space between reconstructed and original activations. For macro layer i 𝑖 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. 10 10 10 The values in the last column are identical ( = 1.414 absent 1.414 =1.414 ) for all layers because ‖ 𝒇 − 𝒇 ^ ‖ / ‖ 𝒇 ‖ = 2 norm 𝒇 ^ 𝒇 norm 𝒇 2 \|\bm{f}-\hat{\bm{f}}\|/\|\bm{f}\|=\sqrt{2} on average for Gaussian random 𝒇 ^ ^ 𝒇 \hat{\bm{f}} provided ‖ 𝒇 ‖ = ‖ 𝒇 ^ ‖ norm 𝒇 norm ^ 𝒇 \|\bm{f}\|=\|\hat{\bm{f}}\| .

layer	image space			activation space
	relative error			relative error
	learned	random	random	learned	random	random
	filters	filters	activations	filters	filters	activations
1	0.423	0.380	0.610	0.895	0.872	1.414
2	0.692	0.438	0.864	0.961	0.926	1.414
3	0.326	0.345	0.652	0.912	0.862	1.414
4	0.379	0.357	0.436	1.051	0.992	1.414

Table 5. Table 5: Best-performing architecture and classification accuracy of random CNNs on CIFAR-10. “( [ 𝚗 ] delimited-[] 𝚗 \mathtt{[n]} ) [ 𝚔 ] delimited-[] 𝚔 \mathtt{[k]} c” denotes a convolution layer with a stride 1 1 1 , a kernel size [ 𝚔 ] delimited-[] 𝚔 \mathtt{[k]} and [ 𝚗 ] delimited-[] 𝚗 \mathtt{[n]} output channels, “ [ 𝚔 ] delimited-[] 𝚔 \mathtt{[k]} p max max {}_{\textrm{max}} ” denotes a max pooling layer with a kernel size [ 𝚔 ] delimited-[] 𝚔 \mathtt{[k]} and a stride [ 𝚔 ] delimited-[] 𝚔 \mathtt{[k]} , and “ [ 𝚔 ] delimited-[] 𝚔 \mathtt{[k]} p ave ave {}_{\textrm{ave}} ” denotes a average pooling layer. A typical layer consists of four operations, namely convolution, ReLU, batch normalization, and max pooling.

Method	#Layers	1 layers	2 layers	3 layers
Random filters	Best	(2048)5c-2p${}_{\text{max}}$-4p${}_{\text{ave}}$	(2048)3c-2p${}_{\text{max}}$-	(2048)3c-2p${}_{\text{max}}$-(2048)3c-
	param.		(2048)3c-2p${}_{\text{max}}$-2p${}_{\text{ave}}$	2p${}_{\text{max}}$-(1024)3c-2p${}_{\text{max}}$
	Accuracy	66.5%	74.6%	74.8%
Learned filters	Best	(1024)5c-2p${}_{\text{max}}$-4p${}_{\text{ave}}$	(1024)3c-2p${}_{\text{max}}$-	(1024)3c-2p${}_{\text{max}}$-(1024)3c-
	param.		(1024)3c-2p${}_{\text{max}}$-2p${}_{\text{ave}}$	2p${}_{\text{max}}$-(1024)3c-2p${}_{\text{max}}$
	Accuracy	68.1%	83.3%	89.3%

Table 6. Table 6: Comparison of coherence between learned filters in each layer of AlexNet and Gaussian random filters with corresponding sizes.

layer	1	2	3	4	5
coherence of learned filters	0.9172	0.6643	0.6200	0.6382	0.3390
coherence of random filters	0.1996	0.1263	0.0929	0.1073	0.1026

Table 7. Table 7: Layer-wise sparsity of AlexNet on ILSVRC-2012 validation set.

layer	conv1	pool1	conv2	pool2	conv3	conv4	conv5	pool5
% of non-zeros	49.41	87.79	18.97	44.13	31.08	30.95	9.78	28.15

Equations

$$(1-{\delta }_k)\Vert \mathbf{z}{\Vert }_2^2\le \Vert \mathbf{\Phi }\mathbf{z}{\Vert }_2^2\le (1+{\delta }_k)\Vert \mathbf{z}{\Vert }_2^2.$$

$$\hat{\mathbf{z}}=\mathrm{upsample}(max-pool(\mathbf{h}),\mathbf{s}),$$

$$\mathbf{\Psi }\left[\begin{array}{c}{\mathbf{h}}_{+}\\ {\mathbf{h}}_{-}\end{array}\right]$$

$$={\mathbf{W}}^T\mathbf{h}={\mathbf{W}}^T\mathbf{W}\mathbf{x}=\mathbf{x}.$$

$$\left[\begin{array}{cc}{\mathbf{W}}^T& -{\mathbf{W}}^T\end{array}\right]\left[\begin{array}{c}{\mathbf{z}}_{+}\\ {\mathbf{z}}_{-}\end{array}\right]={\mathbf{W}}^T({\mathbf{z}}_{+}-{\mathbf{z}}_{-})={\mathbf{W}}^T\mathbf{z}=\mathbf{x}.$$

$$\frac{M\ell^{2}}{D}\geq\frac{C}{\delta_{k}^{2}}\Big{(}k(\log(K)+\log(n))-\log(\epsilon)\Big{)}$$

$$\Vert \hat{\mathbf{x}}-\mathbf{x}{\Vert }_2\le \frac{5{\delta }_{2k}}{1-{\delta }_k}\frac{\sqrt{1+{\delta }_{2k}}}{\sqrt{1-{\delta }_{2k}}}\Vert \mathbf{x}{\Vert }_2.$$

$$\frac{M\ell^{2}}{D}\geq\frac{C}{\delta_{k}^{2}}\Big{(}k(\log(K)+\log(n))-\log(\epsilon)\Big{)}$$

$$\mathbb{E}(\Vert {\mathbf{W}}^T\mathbf{z}{\Vert }_2^2)=\Vert \mathbf{z}{\Vert }_2^2$$

$$\mathbb{E}\left(\Vert {\mathbf{W}}^T\mathbf{z}{\Vert }_2^2\right)$$

$$\Big{|}\|\bm{y}\|_{2}^{2}-\|\bm{z}\|^{2}_{2}\Big{|}>\delta\|\bm{z}\|^{2}_{2}.$$

$$\mathbb{E}(\exp (tZ))\le \exp (t^2{C}^2)$$

$$\|Z\|_{\psi_{2}}=\sup_{p\geq 1}\frac{1}{\sqrt{p}}\Big{(}\mathbb{E}|Z|^{p}\Big{)}^{1/p}.$$

$$\mathbb{P}\Big{(}|X|>t\Big{)}\leq\exp(1-t/C)$$

$${\mathbf{y}}_i=\sum _{j=1}^{Kn}{\mathbf{W}}_{i,j}{\mathbf{z}}_j$$

$$\mathbb{P}\Bigg{(}\Big{|}\sum_{i=1}^{MD}\bm{a}_{i}(Y_{i}^{2}-\mathbb{E}Y_{i}^{2})\Big{|}\geq t\Bigg{)}\leq 2\exp\Bigg{(}-C\min\Big{(}\frac{t^{2}}{T^{2}\|\bm{a}\|_{2}^{2}},\frac{t}{T\|\bm{a}\|_{\infty}}\Big{)}\Bigg{)}.$$

$$\mathbb{P}\Big{(}\Big{|}\|\bm{y}\|_{2}^{2}-\|\bm{z}\|_{2}^{2}\Big{|}>\delta\|\bm{z}\|_{2}^{2}\Big{)}=\mathbb{P}\Bigg{(}\Big{|}\sum_{i=1}^{MD}\|\bm{y}_{i}\|^{2}_{\psi_{2}}(\tilde{\bm{y}}^{2}_{i}-\mathbb{E}\tilde{\bm{y}}^{2}_{i})\Big{|}>\delta\|\bm{z}\|_{2}^{2}\Bigg{)}.$$

$$\Vert \mathbf{a}{\Vert }_2^2=\sum _{i=1}^{MD}\Vert {\mathbf{y}}_i{\Vert }_{{\psi }_2}^4\le \frac{CD\Vert w{\Vert }_{{\psi }_2}^4\Vert \mathbf{z}{\Vert }_2^4}{M{\ell }^2}\text{and}\Vert \mathbf{a}{\Vert }_{\infty }\le \frac{C\Vert w{\Vert }_{{\psi }_2}^2\Vert \mathbf{z}{\Vert }_2^2}{M\ell }.$$

$$\mathbb{P}\Big{(}\Big{|}\|\bm{y}\|_{2}^{2}-\|\bm{z}\|_{2}^{2}\Big{|}>\delta\|\bm{z}\|_{2}^{2}\Big{)}\leq 2\exp\Bigg{(}-C\min\Big{(}\frac{\delta^{2}M\ell^{2}}{D\|w\|_{\psi_{2}}^{4}},\frac{\delta M\ell}{\|w\|_{\psi_{2}}^{2}}\Big{)}\Bigg{)}.$$

$$\exp\Big{(}-C\frac{M\ell^{2}\delta^{2}}{D\|w\|_{\psi_{2}}^{2}}\Big{)}\exp\Big{(}k(\log(K)+\log(n))\Big{)}\leq\exp(\log(\epsilon))$$

$$\frac{M\ell^{2}}{D}\geq\frac{\|w\|_{\psi_{2}}^{2}}{\delta^{2}}\Bigg{(}k(\log(K)+\log(n))-\log(\epsilon)\Bigg{)}=\frac{C}{\delta^{2}}\Big{(}k(\log(K)+\log(n))-\log(\epsilon)\Big{)}.$$

$$D\ge \frac{C}{{\delta }^2}\left(k(\log (K)+\log (n))-\log (ϵ)\right).$$

$$\Vert \mathbf{z}{\Vert }_2\le \frac{1}{\sqrt{1-\delta }}\Vert \mathbf{x}{\Vert }_2.$$

$$\mathbf{h}=\mathbf{W}\mathbf{x}=\mathbf{W}{\mathbf{W}}^T\mathbf{z},$$

$$\Vert \hat{\mathbf{x}}-\mathbf{x}{\Vert }_2\le \frac{5{\delta }_{2k}}{1-{\delta }_k}\frac{\sqrt{1+{\delta }_{2k}}}{\sqrt{1-{\delta }_{2k}}}\Vert \mathbf{x}{\Vert }_2.$$

$$\Vert {\mathbf{W}}_{\Omega }\mathbf{x}{\Vert }_2$$

$$\Vert {\mathbf{W}}_{\Omega }{\mathbf{W}}_{\Omega }^T\mathbf{z}{\Vert }_2$$

$$\Vert {\mathbf{W}}_{\Omega }{\mathbf{W}}^T\mathbf{z}{|}_{{\Omega }^c}{\Vert }_2\le {\delta }_{2k}\Vert \mathbf{z}{|}_{{\Omega }^c}{\Vert }_2.$$

$$\Vert \mathbf{z}{|}_{{\Omega }^c}{\Vert }_2\le \frac{2{\delta }_{2k}}{1-{\delta }_k}\Vert \mathbf{z}{\Vert }_2.$$

$$\Vert \mathbf{h}-{\mathbf{h}}_{\Omega }{\Vert }_2\le \Vert \mathbf{h}-{\mathbf{h}}_{\Pi }{\Vert }_2,$$

$$\Vert \mathbf{h}{|}_{\Omega \setminus \Pi }{\Vert }_2^2\ge \Vert \mathbf{h}{|}_{\Pi \setminus \Omega }{\Vert }_2^2.$$

Full text

1 Introduction

Deep learning has achieved remarkable success in many technological areas, including 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), and computer vision, in particular with deep Convolutional Neural Networks (CNNs) (LeCun et al., 1989; Krizhevsky et al., 2012; Simonyan and Zisserman, 2015; Szegedy et al., 2015).

Following the unprecedented success of deep networks, there have been some theoretical works (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 address the gap between the empirical success and theoretical understanding of the CNNs, in particular its invertibility (i.e., reconstructing the input from the hidden activations), by analyzing a simplified mathematical model using random weights (See Section 2.1 and 4.1 for the practical relevance of the assumption).

This property is intriguing because CNNs are typically trained with discriminative objectives (i.e., unrelated to reconstruction) with a large amount of labels, such as the ImageNet dataset (Deng et al., 2009). Bruna et al. (2014) studied signal discovery from generalized pooling operators using image patches on non-convolutional small scale networks and datasets. Dosovitskiy and Brox (2016) used upsampling-deconvolutional architectures to invert the hidden activations of feedforward CNNs to the input domain. In another related work, Zhao et al. (2016) proposed a stacked what-where autoencoder network and demonstrated its promise in unsupervised and semi-supervised settings. Zhang et al. (2016) showed that CNNs discriminately trained for image classification (e.g., VGGNet (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:

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. 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. 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 version of model-based compressive sensing (and its recovery algorithms) and CNNs. Using the connection, we give a reconstruction bound for a single layer in CNNs, which can be possibly extended to multiple layers. In the experimental sections, we show empirically that large-scale CNNs are consistent with our mathematical analysis. This paper explores these properties and elucidates specific empirical aspects that further mathematical models might need to take into account.

2 Preliminaries

In this section, we begin with discussion on the effectiveness of random weights in CNNs, and then we provide the notations for CNNs, compressive sensing, and sparse signal recovery.

2.1 Effectiveness of Gaussian Random Filters

CNNs with Gaussian random filters have been shown to be 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. Also,  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 (Krizhevsky, 2009) in Section 4.1. Although the performance is not the state-of-the-art, it is surprisingly good considering that the networks are almost untrained. Our theoretical results may provide a new perspective on explaining these phenomena.

2.2 Convolutional Neural Nets

For simplicity, we vectorize input signals to 1-d signal; for any operations we would ordinarily carry out on images, we do on vectors with the appropriate modifications. We define a single layer of our CNN as follows. We assume that the input signal $\bm{x}$ consists of $M$ channels, each of length $D$, and we write $\bm{x}\in\mathbb{R}^{MD}$. For each of the input channels, $m=1,\ldots,M$, let $\bm{w}_{i,m}$, $i=1,\ldots,K$ denote one of $K$ filters, each of length $\ell$. 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-\ell)/t+1$, is an integer. Let $\bm{w}_{i,m}^{j}$ be a vector of length $D$ that consists of the $(i,m)$-th filter shifted by $jt$, $j=0,\ldots,n-1$ (i.e., $\bm{w}_{i,m}^{j}$ has at most $\ell$ non-zero entries). We concatenate over the $M$ channels each of these vectors (as row vectors) to form a large matrix $\bm{W}$, which is the $Kn\times MD$ matrix made up of $K$ blocks of the $n$ shifts of each filter in each of $M$ channels. We assume that $Kn\geq MD$ and the $Kn$ row vectors of $\bm{W}$ span $\mathbb{R}^{MD}$ and that we have normalized the rows so that they have unit $\ell_{2}$ norm. The hidden units of the feed-forward CNN are computed by multiplying an input signal $\bm{x}\in\mathbb{R}^{MD}$ by the matrix $\bm{W}$ (i.e., convolving, in each channel, by a filter bank of size $K$, and summing over the channels to obtain $Kn$ outputs).111Convolution can be computed more efficiently than matrix multiplication, but they are mathematically equivalent. We use $\bm{h}=\bm{W}\bm{x}$ for the hidden activation computed by a single layer CNN without pooling. Figure 1 illustrates the architecture. As a nonlinear activation, we apply the ${\rm 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.

2.3 Compressive Sensing

In compressive sensing, we assume that there is a latent sparse code $\bm{z}$ that generates the visible signal $\bm{x}$. We say that a $p\times q$ matrix $\bm{\Phi}$ with $q>p$ satisfies the Restricted Isometry Property RIP$(k,\delta_{k})$ if there is a distortion factor $\delta_{k}>0$ such that for all $\bm{z}\in\mathbb{R}^{q}$ with exactly $k$ non-zero entries, $(1-\delta_{k})\|\bm{z}\|^{2}_{2}\leq\|\bm{\Phi}\bm{z}\|^{2}_{2}\leq(1+\delta_{k})\|\bm{z}\|^{2}_{2}$. If $\bm{\Phi}$ satisfies RIP with sufficiently small $\delta_{k}$ and if $\bm{z}$ is $k$-sparse, then given the vector $\bm{x}=\bm{\Phi}\bm{z}\in\mathbb{R}^{p}$, we can efficiently recover $\bm{z}$ (see Candés (2008) for more details)222We note that this is a sufficient condition and that there are other, less restrictive sufficient conditions, as well as more complicated necessary conditions. . There are many efficient algorithms, including $\ell_{1}$ sparse coding (e.g., $\ell_{2}$ minimization with $\ell_{1}$ regularization) and greedy and iterative algorithms, such as Iterative Hard Thresholding (IHT) (Blumensath and Davies, 2009).

Model-based compressive sensing. While sparse signals are a natural model for some applications, they are less realistic for CNNs. We consider a vector $\bm{z}\in\mathbb{R}^{Kn}$ as the true sparse code generating the CNN input $\bm{x}$ with a particular model of sparsity. Rather than permitting $k$ non-zero entries anywhere in the vector $\bm{z}$, we divide the support of $\bm{z}$ into $K$ contiguous blocks of size $n$ and we stipulate that from each block there is at most one non-zero entry in $\bm{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 $k$-sparse subspaces with this structure ${\cal M}_{k}$. It is clear that ${\cal M}_{k}$ contains $n^{k}\binom{K}{k}$ subspaces. In our analysis, we consider linear combinations of two model-$k$-sparse signals. To be precise, suppose that $\bm{z}=\alpha_{1}\bm{z}_{1}+\alpha_{2}\bm{z}_{2}$ is the linear combination of two elements in ${\cal M}_{k}$. Then, we say that $\bm{z}$ lies in the linear subspace ${\cal M}^{2}_{k}$ that consists of all linear combinations of vectors from ${\cal M}_{k}$.333Intuitively, ${\cal M}^{2}_{k}$ is a subspace where the error signal $\hat{\bm{z}}-\bm{z}$ lies in and used for reconstruction bound derivation; see Appendix A. We say that a matrix $\bm{\Phi}$ satisfies the model-RIP if there is a distortion factor $\delta_{k}>0$ such that, for all $\bm{z}\in{\cal M}_{k}$,

[TABLE]

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 don’t consider in this paper to keep our analysis simple. Intuitively, a matrix satisfying the model-RIP is a nearly orthonormal matrix of a particular set of sparse vectors with a particular sparsity model or pattern.

For our analysis, we also need matrices $\bm{\Phi}$ that satisfy the model-RIP for vectors $\bm{z}\in{\cal M}^{2}_{k}$. We denote the distortion factor $\delta_{2k}$ for such matrices; note that $\delta_{k}\leq\delta_{2k}<1$.

Many efficient algorithms have been proposed for sparse coding and compressive sensing (Olshausen and others, 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 (Baraniuk et al., 2010), assuming the existence of an efficient structured sparse approximation algorithm $\mathbb{M}$, that given an input vector and the sparsity parameter, returns the vector closest to the input with the specified sparsity structure.

In CNNs, 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 downsampling step by an unsampling operator. This procedure defines our structured sparse approximation algorithm $\hat{\bm{z}}=\mathbb{M}(\bm{h},k)$, where $\bm{h}$ is the original (unpooled) response, $k$ is the sparsity parameter for block-sparsification, and $\hat{\bm{z}}$ is the sparsified response after pooling but without shrinking the length (i.e., the locally non-maximum values are zeroed out such that $\hat{\bm{z}}$ has the same dimension as $\bm{h}$). Note that $\hat{\bm{z}}$ is a model-$k$-sparse signal by construction. On the other hand, without considering the block-sparsification, we actually apply the following max pooling and upsampling operations:

[TABLE]

where $\hat{\bm{z}}$ is the pooled response, $\bm{h}$ is the filter response of CNN given input before max pooling (see Section 2.2), and $\bm{s}$ denotes the upsampling switches that indicate whereto place the non-zero values in the upsampled activations. Since our theoretical analysis does not depend on $\bm{s}$ but depends on $k$, any type of valid upsampling switches will be consistent with the block-sparsification (model-$k$-sparse) assumption, thus we will use $\mathbb{M}(\bm{h},k)$ to denote the structured sparse approximation algorithm (2) without worrying about $\bm{s}$.

We use model-sparse version of IHT (Blumensath and Davies, 2009) as our recovery algorithm, as one iteration of IHT for our model of sparsity captures exactly a feedforward CNN.444Multiple 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. Algorithm 1 describes the model-based IHT algorithm. In particular, the sequence of steps 4–6 in the middle IHT 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

Following the idea of compressive sensing in Section 2.3, we assume that the input $\bm{x}$ is generated from a latent model-$k$-sparse signal $\bm{z}$ with basis vectors $\bm{\Phi}$, which turns out to be $\bm{W}^{T}$ by Theorem 3.1 (i.e., $\bm{x}=\bm{W}^{T}\bm{z}$). Therefore, our analysis views the output of CNN (with pooling) is a reconstruction of $\bm{z}$ (i.e., $\hat{\bm{z}}=\mathbb{M}(\bm{W}\bm{x},k)$), and $\bm{W}^{T}$ can be used to reconstruct $\bm{x}$ from $\hat{\bm{z}}$: that is, $\hat{\bm{x}}=\bm{W}^{T}\hat{\bm{z}}$.

3.1 CNN Filters with Positive and Negative Pairs

Here we assume that all of the entries in the vectors are real numbers rather than only non-negative like when using ${\rm ReLU}$. This setup is equivalent to using Concatenated ${\rm ReLU}$ (CReLU) (Shang et al., 2016) as an activation function (i.e., keeping the positive and negative activations as separate hidden units) with tied decoding weights. The CReLU activation scheme is justified by the fact that trained CNN filters come in positive and negative pairs and that it achieves superior classification performance in several benchmarks. This setting makes a CNN much easier to analyze within the model compressed sensing framework.

To motivate the setting, we begin with a simple example. Suppose that the matrix $\bm{W}$ is an orthonormal basis for $\mathbb{R}^{MD}$ and define $\bm{\Psi}=\begin{bmatrix}\bm{W}^{T}&-\bm{W}^{T}\end{bmatrix}$.

Proposition 1**.**

A one-layer CNN using the matrix $\bm{\Psi}^{T}$, with no pooling, gives perfect reconstruction (with the matrix $\bm{\Psi}$) for any input vector $\bm{x}\in\mathbb{R}^{MD}$.

Proof.

Because we have both the positive and the negative dot products of the signal with the basis vectors in ${\rm ReLU}(\bm{\Psi}^{T}\bm{x})={\rm ReLU}\left(\begin{bmatrix}\bm{W}\bm{x}\\ -\bm{W}\bm{x}\end{bmatrix}\right)$, we have positive and negative versions of the hidden units $\bm{h}_{+}={\rm ReLU}(\bm{W}\bm{x})$ and $\bm{h}_{-}={\rm ReLU}(-\bm{W}\bm{x})$ where we decompose $\bm{h}=\bm{W}\bm{x}=\bm{h}_{+}-\bm{h}_{-}$ into the difference of two non-negative vectors, the positive and the negative entries of $\bm{h}$. From this decomposition, we can easily reconstruct the original signal via

[TABLE]

∎

In the example above, we have pairs of vectors $(\bm{w},-\bm{w})$ in our matrix $\bm{\Psi}$. Now suppose that we have a vector $\bm{z}$ where its positive and negative components can be split into $\bm{z}=\bm{z}_{+}-\bm{z}_{-}$, and that we synthesize a signal $\bm{x}$ from $\bm{z}$ using the matrix $\begin{bmatrix}\bm{W}^{T}&-\bm{W}^{T}\end{bmatrix}$. Then, we have

[TABLE]

Next, we multiply $\bm{x}=\bm{W}^{T}\bm{z}$ by a concatenation of positive and negative $\bm{W}$, then we get $\begin{bmatrix}\bm{W}\\ -\bm{W}\\ \end{bmatrix}\bm{x}=\begin{bmatrix}\bm{W}\bm{W}^{T}\bm{z}\\ -\bm{W}\bm{W}^{T}\bm{z}\end{bmatrix}$ and if we apply ${\rm ReLU}$ to this vector, we get $\begin{bmatrix}(\bm{W}\bm{W}^{T}\bm{z})_{+}\\ (\bm{W}\bm{W}^{T}\bm{z})_{-}\end{bmatrix}$, which is a vector $\bm{W}\bm{W}^{T}\bm{z}$ that is split into its positive and negative components. The structure of the product $\bm{W}\bm{W}^{T}$ is crucial to the reconstruction quality of the vector $\bm{z}$. In addition, this calculation shows that if we have both positive and negative pairs of filters or vectors, then the ${\rm 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 assume that all of the entries in the vectors are real numbers, rather than only non-negative.

3.2 Model-RIP and Random Filters

Our first main result shows that if we use Gaussian random filters in our CNN, then, with high probability, $\bm{W}^{T}$, the transpose of a large matrix formed by the convolution filters satisfies the model-RIP. In other words, Gaussian random filters generate a matrix whose transpose $\bm{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 $\ell$ 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 $\delta_{k}$, and the failure probability $\epsilon$. The proof is in Appendix A.

Theorem 3.1**.**

Assume that we have $MK$ vectors $\bm{w}_{i,m}$ of length $\ell$ in which each entry is a scaled i.i.d. (sub-)Gaussian random variable with zero mean and unit variance (the scaling factor is $1/\sqrt{M\ell})$. Let $t$ be the stride length (where $n=(D-\ell)/t+1$) and $\bm{W}$ be a structured random matrix, which is the weight matrix of a single layer CNN with $M$ channels and input length $D$. If

[TABLE]

for a positive constant $C$, then with probability $1-\epsilon$, the $MD\times Kn$ matrix $\bm{W}^{T}$ satisfies the model-RIP for model ${\cal M}_{k}$ with parameter $\delta_{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 satisfy, with high probability, the model-RIP for ${\cal M}^{2}_{k}$.

Other examples of matrices that satisfy the model-RIP include wavelets and localized Fourier bases; both examples can be easily and efficiently implemented via convolutions.

3.3 Reconstruction Bounds

Suppose $\bm{W}^{T}$ satisfies the model-RIP and $\hat{\bm{z}}$ is the reconstruction of true sparse code $\bm{z}$ through a CNN layer followed by pooling, i.e., $\hat{\bm{z}}=\mathbb{M}(\bm{W}\bm{x},k)$. Then, Theorem 3.3 shows that $\hat{\bm{x}}=\bm{W}^{T}\hat{\bm{z}}$ is an approximate reconstruction of the input signal, and the relative error is bounded on a function of the distortion parameters of the model-RIP.

Theorem 3.3**.**

We assume that $\bm{W}^{T}$ satisfies the ${\cal M}_{k}^{2}$-RIP with constant $\delta_{k}\leq\delta_{2k}<1$. If we use $\bm{W}$ in a single layer CNN both to compute the hidden units $\hat{\bm{z}}$ and to reconstruct the input $\bm{x}$ from these hidden units as $\hat{\bm{x}}$ so that $\hat{\bm{x}}=\bm{W}^{T}\mathbb{M}(\bm{W}\bm{x},k)$, the error in our reconstruction is

[TABLE]

See Appendix B for the detailed proofs. Part of our analysis also shows that the hidden units $\hat{\bm{z}}$ are approximately the putative coefficient vector $\bm{z}$ in the sparse linear representation for the input signal. Recall that the structured sparsity approximation algorithm $\mathbb{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. We can extend the analysis for a single layer CNN to multiple layer CNN by using the output on one layer as the input to another, following the proof in Appendix B. We leave further investigation of this idea as future work.

4 Experimental Evidence and Analysis

In this section, we provide experimental validation of our theoretical model and analysis. We first validate the practical relevance of our assumption by examining the effectiveness of random filter CNNs, and then provide results on more realistic scenarios. In particular, we study popular deep CNNs trained for image classification on ILSVRC 2012 dataset (Deng et al., 2009). We calculate empirical model-RIP bounds for $\bm{W}^{T}$, showing that they are consistent with our theory. Our results are also consistent with a long line of research shows that it is reasonable to model real and natural images as sparse linear combinations overcomplete dictionaries (Boureau et al., 2008; Le et al., 2013; Lee et al., 2008; Olshausen and others, 1996; Ranzato et al., 2007; Yang et al., 2010). In addition, we verify our theoretical bounds for the reconstruction error $\|\bm{x}-\bm{W}^{T}\hat{\bm{z}}\|_{2}/\|\bm{x}\|_{2}$ on real images. We investigate both randomly sampled filters and empirically learned filters in these experiments. Our implementation is based on Caffe (Jia et al., 2014) and MatConvNet (Vedaldi and Lenc, 2015).

Recall that our theoretical analysis is generic to any upsampling switches in (2) for reconstruction. In the experiments, we specifically use the naive upsampling to reverse max-pool activations to its original size, where only the first element in a pooling region is assigned with the pooled activation, and the rest elements are all zero. Thus, no extra information other than the pooled activation values are taken into account.

4.1 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. zero mean unit variance entries on the CIFAR-10 (Krizhevsky, 2009). Note that 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.555Implementation 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. See Appendix C.1 for more details. Specifically, we test random CNNs with 1, 2, and 3 convolutional layers followed by ReLU activation and $2\times 2$ max pooling layer. We tested different filter sizes ($3,5,7$) and numbers of channels ($64,128,256,1024,2048$) and report the best classification accuracy by 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 Appendix C.1. We observe that CNNs with Gaussian random filters achieve good classification performance (implying that they serve as reasonable representation of input data), which is not too far off the learned filters. Our experimental results are also consistent with the observations made by Jarrett et al. (2009) and Saxe et al. (2011). In conclusion, those results suggest that 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.

4.2 1-d Model-RIP

We use 1-d synthetic data to empirically show the basic validity of our theory in terms of the model-RIP in (1) and reconstruction bound in Theorem 3.3. We plot the histograms of the empirical model-RIP values of 1-d Gaussian random filters $\bm{W}$ ( scaled by $1/\sqrt{lM}$ ) with size $l\times 1\times M\times K=5\times 1\times 32\times 96$ on 1-d ${\cal M}_{k}$ sparse signal $\bm{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 2 (a)–(b) are tightly centered around $1$, suggesting that $\bm{W}^{T}$ satisfies the model-RIP in (1) and its corollary from Lemma B.1, respectively. We also empirically show the reconstruction bound in Theorem 3.3 on synthetic vectors $\bm{x}=\bm{W}^{T}\bm{z}$ (Figure 2 (c)). The reconstruction error is concentrated at around $0.1$–$0.2$ and bound under $0.5$. Results in Figure 2 suggest the practical validity of our theory when the model assumptions hold.

4.3 Architectures for 2-d Model-RIP

We conduct the rest of our experimental evaluations on the 16-layer VGGNet (Model D in Simonyan and Zisserman (2015)), 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 image 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 Appendix C.2.

VGGNet contains five macro layers of convolution and pooling layers, and each macro layer has 2 or 3 convolutional layers followed by a pooling layer. We denote the $j$-th convolutional layer in the $i$-th macro layer “conv$(i,j)$,” and the pooling layer “pool$(i)$.” The activations/features from $i$-th macro layer are the output of pool$(i)$. Our analysis is for single convolutional layers.

4.4 2-d Model-RIP

The key to our reconstruction bound is Theorem 3.3. We empirically evaluate the model-RIP, i.e., $\|\bm{W}^{T}\bm{z}\|/\|\bm{z}\|$, for real CNN filters of the pretrained VGGNet. We use two-dimensional coefficients $\bm{z}$ (each block of coefficients is of size $D\times D$), $K$ filters of size $\ell\times\ell$, and pool the coefficients over smaller pooling regions (i.e., not over all possible shifts of each filter). The following experimental evidence suggests that the sparsity model and the model-RIP of the filters are consistent with our mathematical analysis on the simpler one-dimensional case.

To check the significance of the model-RIP (i.e., how close $\|\bm{W}^{T}\bm{z}\|/\|\bm{z}\|$ is to $1$) in controlled settings, we first synthesize the hidden activations $\bm{z}$ with sparse uniform random variables, which fully agree with our model assumptions.

The sparsity of $\bm{z}$ is constrained to the average level of the real CNN activations, which is reported in Table 2. Given the filters of a certain convolutional layer, we use the synthetic $\bm{z}$ (in equal position to this layer’s output activations) to get statistics for the model-RIP. To be consistent with succeeding experiments, we choose conv$(5,2)$, while other layers show similar results. Figure 3 (a) summarizes the distribution of empirical model-RIP values, which is clearly centered around $1$ and satisfies (1) with a short tail roughly bounded by $\delta_{k}<1$. For more details of the algorithm, we normalize the filters from the conv$(5,2)$ layer, which are $\ell\times\ell$ ($\ell=3$). All $K=512$ filters with $M=512$ input channels are used.666We do not remove any filters including those in approximate positive/negative pairs (see Section 3.) We set $D=15$, which is the same as the output activations of conv$(5,2)$, and use $2\times 2$ pooling regions777No pooling layer follows conv$(5,2)$ in VGGNet. However, we use it in this way to analyze the convolution-pooling pair per theory., which is commonly used in recent CNNs. We generate 1000 ${\cal M}_{k}$ randomly sampled sparse activation maps $\bm{z}$ 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.

More realistically, we observe that 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 on the hidden activations $\bm{z}$ that reconstruct the actual input activations $\bm{x}$ from conv$(5,1)$ by $\bm{W}^{T}\bm{z}$. Per theory, the $\bm{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 $\bm{z}$ from $\bm{x}$ in Algorithm 2. The algorithm is inspired by “$\ell_{1}$ heuristic" method that are commonly used in practice (e.g.,  Boyd (2015)). As shown in Algorithm 2, we first do $\ell_{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 max pooling and unpooling with known switches (line 2) to zero out the locally non-maximum values without messing up the support structures. We perform $\ell_{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 values for visual activations $\bm{x}$ from conv$(5,1)$ with/without ${\rm ReLU}$ are both close to $1$. The center offset to $1$ is less than $0.05$ and the range bound $\delta_{k}$ is roughly less than $0.05$, which agrees with the theoretical bound in (1). To gain more insight, we summarize the learned filter coherence in Table 3 for all convolutional layers in VGGNet.888The coherence is defined as the maximum (in absolute value) dot product between distinct pairs of columns of the matrix $W^{T}$, i.e. $\mu=\max_{i\not=j}|W_{i}W_{j}^{T}|$, where $W_{i}$ denote the $i$-th row of matrix $W$. This measures the correlation or similarity between the columns of $W^{T}$ and is a proxy for the value of the model-RIP parameter $\delta_{k}$ (which we can only estimate computationally). The smaller the coherence, the smaller $\delta_{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 turns out to be robust to this mismatch. It demonstrates the strong practical invertibility of CNN.

4.5 Reconstruction Bounds

With model-RIP as a sufficient condition, Theorem 3.3 provides a theoretical bound for layer-wise reconstruction via $\hat{\bm{x}}=\bm{W}^{T}\mathbb{M}(\bm{W}\bm{x})$, which 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 pool$(1)$ to $(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 pool$(4)$. Interestingly, the decoding network itself is quite powerful, since it can reconstruct the rough (although very noisy) glimpse of images with Gaussian random input, as shown in Figure 4 (e). Object shapes are recovered up to some extent 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 pool$(4)$ 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 3). 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 Appendix C.3.

In Table 10, we summarize reconstruction performance for all 4 macro 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.

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. We leave such explorations as future work.

Acknowledgments

This work was supported in part by ONR N00014-16-1-2928, NSF CAREER IIS-1453651, and Sloan Research Fellowship. We would like to thank Michael Wakin for helpful discussions about concentration of measure for structured random matrices.

Appendix A Mathematical Analysis: Model-RIP and Random Filters

Theorem 3.1(Restated) Assume that we have $MK$ vectors $\bm{w}_{i,m}$ of length $\ell$ in which each entry is a scaled i.i.d. (sub-)Gaussian random variable with zero mean and unit variance (the scaling factor is $1/\sqrt{M\ell})$. Let $t$ be the stride length (where $n=(D-\ell)/t+1$) and $\bm{W}$ be a structured random matrix, which is the weight matrix of a single layer CNN with $M$ channels and input length $D$. If

[TABLE]

*for a positive constant $C$, then with probability $1-\epsilon$, the $MD\times Kn$ matrix $\bm{W}^{T}$ satisfies the model-RIP for model ${\cal M}_{k}$ with parameter $\delta_{k}$. *

Proof.

We note that the proof follows the same structure of those in other papers such as Park et al. (2011) and Vershynin (2010), though we make minor tweaks to account for the particular structure of $\bm{W}^{T}$.

Suppose that $\bm{z}\in{\cal M}_{k}$, i.e., $\bm{z}$ consists of at most $k$ non-zero entries that each appears in a distinct block of size $n$ (there are a total of $K$ blocks). First, Lemma A.1 shows that the expectation of the norm of $\bm{W}^{T}\bm{z}$ is preserved.

Lemma A.1**.**

[TABLE]

Proof.

Note that each entry of $\bm{W}^{T}$ is either zero or Gaussian random variable $w\sim N(0,1)$ before scaling. Therefore, it is obvious that $\mathbb{E}(\bm{W}\bm{W}^{T})=\bm{I}$ since each row of $\bm{W}$ satisfies $\mathbb{E}\left(\left(\bm{w}_{i_{1},m}^{j_{1}}\right)^{T}\left(\bm{w}_{i_{2},m}^{j_{2}}\right)\right)=0$ if $j_{1}\neq j_{2}$ or $i_{1}\neq i_{2}$ for any $m=1\dots M$, and we normalized the random variables so that $\mathbb{E}\left(\left\|\left[\left(\bm{w}_{i,1}^{j}\right)^{T},\ldots,\left(\bm{w}_{i,M}^{j}\right)^{T}\right]\right\|_{2}\right)=1$ for all $i,j$’s. Finally, we have

[TABLE]

∎

Let $\bm{y}=\bm{W}^{T}\bm{z}$. We aim to show that the square norm of the random variable $\|\bm{y}\|_{2}^{2}$ concentrates tightly about its mean; i.e., with exceedingly low probability

[TABLE]

To do so, we need several properties of sub-Gaussian and sub-exponential random variables. A mean-zero sub-Gaussian random variable $Z$ has a moment generating function that satisfies

[TABLE]

for all $t\in\mathbb{R}$ and some constant $C$. The sub-Gaussian norm of $Z$, denoted $\|Z\|_{\psi_{2}}$ is

[TABLE]

If $Z\sim N(0,\sigma^{2})$, then $\|Z\|_{\psi_{2}}\leq c\sigma$ where $c$ is a positive constant (following Definition 5.7 in Vershynin (2010)).

A sub-exponential random variable $X$ satisfies111111There are two other equivalent properties. See Vershynin (2010) for details.

[TABLE]

for all $t\geq 0$.

Let $\bm{y}_{i}$ denote the $i$th entry of the vector $\bm{y}=\bm{W}^{T}\bm{z}$. We can write

[TABLE]

and observe that $\bm{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 zero mean and sub-Gaussian norm $\|\bm{y}_{i}\|_{\psi_{2}}\leq C/\sqrt{M\ell}\|w\|_{\psi_{2}}\|\bm{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 $\bm{W}$ do enter the more refined bound on the sub-Gaussian norm of $\|\bm{y}_{i}\|_{\psi_{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 $\|\bm{y}_{i}\|^{2}_{2}$, we recall from Park et al. (2011) and Vershynin (2010) that the sum of squares of sub-Gaussian random variables tightly concentrate.

Theorem A.2**.**

Let $Y_{1},\ldots,Y_{MD}$ be independent sub-Gaussian random variables with sub-Gaussian norms $\|Y_{i}\|_{\psi_{2}}$ for all $i=1,\ldots,MD$. Let $T=\max_{i}\|Y_{i}\|_{\psi_{2}}$. For every $t\geq 0$ and every $\bm{a}\in\mathbb{R}^{MD}$ and a positive constant $C$,

[TABLE]

We note that although some entries $\bm{y}_{i}$ may be identically zero, depending on the sparsity pattern of $\bm{z}$, not all entries are. Let us define $\tilde{\bm{y}}_{i}=\frac{\bm{y}_{i}}{\|\bm{y}_{i}\|_{\psi_{2}}}$ so that $\|\tilde{\bm{y}}_{i}\|_{\psi_{2}}=1$.

From Lemma A.1 and the relation $\bm{y}=\bm{W}^{T}\bm{z}$, we have

[TABLE]

See Proposition 5.16 in Vershynin (2010) for the proof of Theorem A.2. We apply Theorem A.2 to the sub-Gaussian random variables $\tilde{\bm{y}}_{i}$ with the weights $\|\bm{y}_{i}\|^{2}_{\psi_{2}}$. We have

[TABLE]

If we set $T=1$, $t=\delta\|\bm{z}\|_{2}^{2}$, and use the above estimates for the norms of $\bm{a}$, we have

[TABLE]

Finally, we use the concentration of measure result in a crude union bound to bound the failure probability over all vectors $\bm{z}\in{\cal M}_{k}$. We take $n^{k}\binom{K}{k}\approx(nK)^{k}$ and $\epsilon$ for a desired constant failure probability. Using the smaller term in (3), (note that $\delta<1$, $\ell/D<1$, and $\|w\|_{\psi_{2}}\leq 1$) we have

[TABLE]

which implies

[TABLE]

Therefore, if design our matrix $\bm{W}$ as described and with the parameter relationship as above, the matrix $\bm{W}^{T}$ satisfies the model-RIP for ${\cal M}_{k}$ and parameter $\delta$ with probability $1-\epsilon$. ∎

Let us discuss the relationship amongst the parameters in our result. First, if we have only one channel $M=1$ and the filter length $\ell=D$ ; namely,

[TABLE]

If $\ell<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 $\ell$ needs to be larger than $\sqrt{D}$, or, if add more channels, $\sqrt{D/M}$.

Appendix B Mathematical Analysis: Reconstruction Bounds

The consequences of having the model-RIP are two-fold. The first is that if we assume that an input image is the structured sparse linear combination of filters, $\bm{x}=\bm{W}^{T}\bm{z}$ where $\bm{z}\in{\cal M}_{k}$ and $\bm{W}^{T}$ satisfies the model-RIP, then we know an upper and lower bound on the norm of $\bm{x}$ in terms of the norm of its sparse coefficients, $\|\bm{x}\|_{2}\leq(1\pm\delta)\|\bm{z}\|_{2}$. Additionally,

[TABLE]

More importantly, when we calculate the hidden units of $\bm{x}$,

[TABLE]

then we can see that the computation of $\bm{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 $\bm{h}$ and $\bm{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 $\bm{W}^{T}$ satisfies the ${\cal M}_{k}^{2}$-RIP with constant $\delta_{k}\leq\delta_{2k}<1$. If we use $\bm{W}$ in a single layer CNN both to compute the hidden units $\hat{\bm{z}}$ and to reconstruct the input $\bm{x}$ from these hidden units as $\hat{\bm{x}}$ so that $\hat{\bm{x}}=\bm{W}^{T}\mathbb{M}(\bm{W}\bm{x},k)$, the error in our reconstruction is

[TABLE]

Proof.

To show this result, we recall the following lemmas from Baraniuk et al. (2010) and rephrase them in the setting of a feedforward CNN. Note that Lemma B.1 and  B.2 are the same as Lemma 1 and 2 in Baraniuk et al. (2010), respectively.

Lemma B.1**.**

Suppose $\bm{W}^{T}$ has ${\cal M}_{k}$-RIP with constant $\delta_{k}$. Let $\Omega$ be a support corresponding to a subspace in ${\cal M}_{k}$. Then we have the following bounds:

[TABLE]

Lemma B.2**.**

Suppose that $\bm{W}^{T}$ has ${\cal M}_{k}^{2}$-RIP with constant $\delta_{2k}$. Let $\Omega$ be a support corresponding to a subspace of ${\cal M}_{k}$ and suppose that $\bm{z}\in{\cal M}_{k}$ (not necessarily supported on $\Omega$). Then

[TABLE]

Let $\Pi$ denote the support of the ${\cal M}_{k}$ sparse vector $\bm{z}$. Set $\bm{h}=\bm{W}\bm{x}$ and set $\hat{\bm{z}}$ to be the result of max pooling applied to the vector $\bm{h}$, or the best fit (with respect to the $\ell_{2}$ norm) to $\bm{h}$ in the model ${\cal M}_{k}$. Let $\Omega$ denote the support set of $\hat{\bm{z}}\in{\cal M}_{k}$. For simplicity, we assume $|\Pi|=k=|\Omega|$.

Lemma B.3** (Identification).**

The support set, $\Omega$, of the switch units captures a significant fraction of the total energy in the coefficient vector $\bm{z}$

[TABLE]

Proof.

Let $\bm{h}_{\Omega}$ and $\bm{h}_{\Pi}$ be the vector $\bm{h}$ restricted to the support sets $\Omega$ and $\Pi$, respectively. Since both are support sets for ${\cal M}_{k}$ and since $\Omega$ is the best support set for $\bm{h}$,

[TABLE]

and, after several calculations, which is identical to those in the proof of Lemma 3 in Baraniuk et al. (2010), we have

[TABLE]

Using Lemma B.2 and the size $|(\Omega\setminus\Pi)\bigcup\Pi|\leq 2k$, we have

[TABLE]

Using (6) and Lemma B.2, we can bound the other side of the inequality as

[TABLE]

Since the support of $\bm{z}$ is the set $\Pi$, $\Pi\setminus\Omega=\Omega^{c}$ for $\bm{z}$, so we can conclude that

[TABLE]

and with some rearrangement, we have

[TABLE]

∎

To set the value of $\hat{\bm{z}}$ on its support set $\Omega$, we simply set $\hat{\bm{z}}=\bm{h}|_{\Omega}$ and $\hat{\bm{z}}|_{\Omega^{c}}=0$. Then

Lemma B.4** (Estimation).**

[TABLE]

Proof.

First, note that $\|\bm{I}-\bm{W}_{\Omega}\bm{W}_{\Omega}^{T}\|_{2}\leq\max\{(1+\delta_{k})-1,1-(1-\delta_{k})\}=\delta_{k}$ since

[TABLE]

where $\sigma_{\max}$ is the maximum singular value. Therefore,

[TABLE]

∎

Finally, if we use the autoencoder formulation to reconstruct the original image $\bm{x}$ by setting $\hat{\bm{x}}=\bm{W}^{T}\hat{\bm{z}}$, we can estimate the reconstruction error. We note that $\hat{\bm{z}}$ is ${\cal M}_{k}$-sparse by construction and remind the reader that $\bm{W}^{T}$ satisfies ${\cal M}_{k}^{2}$-model-RIP with constants $\delta_{k}\leq\delta_{2k}\ll 1$. Then, using Lemma B.4 as well as the ${\cal M}_{k}^{2}$-sparse properties of $\bm{W}^{T}$,

[TABLE]

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. ∎

Appendix C More Experimental Results

C.1 More Details on Evaluation of CNNs with Gaussian Random Filters

In this section, we provide more details on the network architectures that we used in Table 1. For the network architecture, 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\times 4$ and a dropout layer for better regularization before feeding the feature to a learnable soft-max classifier for image classification. We describe the best performing architectures for all cases in Table 5.

C.2 Layer-wise Coherence and Sparsity for AlexNet

We present coherence (see Table 6) and sparsity level (see Table 7) for each layer in AlexNet.

C.3 Visualization of Image Reconstruction for VGGNet

In Figure 5, we show reconstructed images from each layer using different reconstruction methods via a pretrained decoding network.