Chapter 4 Deep Representations as Unrolled Optimization

From the previous chapter, we see that to seek a linear discriminative representation (LDR), mathematically, we are essentially seeking a continuous mapping $f(\cdot):\bm{x}\mapsto\bm{z}$ from the data $\bm{X}=[\bm{x}_{1},\ldots,\bm{x}_{N}]\in\mathbb{R}^{D\times N}$ (or initial features extracted from the data³³3As we will see the necessity of such a feature extraction in the next section.) to an optimal representation $\bm{Z}=[\bm{z}_{1},\ldots,\bm{z}_{N}]\in\mathbb{R}^{d\times N}$ that maximizes the following coding rate reduction objective:

where $\epsilon>0$ is a prescribed quantization error and for simplicity we denote⁴⁴4Notice our use of slightly simplified notation compared to Chapter 3.

The question really boils down to whether there is a constructive way of finding such a continuous mapping $f(\cdot,\bm{\theta})$ from $\bm{x}$ to $\bm{z}$? To this end, let us consider incrementally maximizing the objective $\Delta R_{\epsilon}$ as a function of $\bm{Z}\subseteq\mathbb{S}^{d-1}$. Although there might be many optimization schemes to choose from, for simplicity we first consider the arguably simplest projected gradient ascent (PGA) scheme:⁵⁵5Notice that we use subscript $j$ on $\bm{Z}_{j}$ to indicate features in the $j$-th class and superscript $\ell$ on $\bm{Z}^{\ell}$ to indicate all features at $\ell$-th iteration or layer.

for some step size $\eta>0$ and the iterate starts with the given data $\bm{Z}^{0}=\bm{X}$.⁶⁶6Again, for simplicity, we here first assume the initial features $\bm{Z}^{1}$ are the data themselves. Note that here $\ell$ denotes the number of iterations. Hence, the data and the features have the same dimension $d$. This needs not to be the case though. As we will see in the next section, the initial features can be some (lifted) features of the data to begin with and could in principle have a different (much higher) dimension. All subsequent iterates have the same dimension. This scheme can be interpreted as how one should incrementally adjust locations of the current features $\bm{Z}^{\ell}$, initialized as the input data $\bm{X}$, in order for the resulting $\bm{Z}^{\ell+1}$ to improve the rate reduction $\Delta R_{\epsilon}$, as illustrated in Figure 4.1.

Simple calculation shows that the gradient ${\partial\Delta R_{\epsilon}}/{\partial\bm{Z}}$ entails evaluating the following derivatives of the two terms in $\Delta R_{\epsilon}$:

Notice that in the above, the matrix $\bm{E}^{\ell}$ only depends on $\bm{Z}^{\ell}$ and it aims to expand all the features to increase the overall coding rate; the matrix $\bm{C}^{\ell}_{k}$ depends on features from the $k$-class and aims to compress them to reduce the coding rate of each class. Then the complete gradient $\frac{\partial\Delta R_{\epsilon}}{\partial\bm{Z}}(\bm{Z}^{\ell})\in\mathbb{R}^{d\times N}$ is of the form:

Notice that in the above, the gradient ascent considers all the features $\bm{Z}^{\ell}=[\bm{z}^{\ell}_{1},\dots,\bm{z}^{\ell}_{N}]$ as free variables. The increment $\bm{Z}^{\ell+1}-\bm{Z}^{\ell}=\eta\frac{\partial\Delta R_{\epsilon}}{\partial\bm{Z}}(\bm{Z}^{\ell})$ does not yet give a transformation on the entire feature domain $\bm{z}^{\ell}\in\mathbb{R}^{d}$. According to equation (4.1.5), the gradient cannot be evaluated at a point whose membership is not known, as illustrated in Figure 4.1. Hence, in order to find the optimal $f(\bm{x},\bm{\theta})$ explicitly, we may consider constructing a small increment transform $g(\cdot,\bm{\theta}^{\ell})$ on the $\ell$-th layer feature $\bm{z}^{\ell}$ to emulate the above (projected) gradient scheme:

such that $\big{[}g(\bm{z}_{1}^{\ell},\bm{\theta}^{\ell}),\ldots,g(\bm{z}_{N}^{\ell},\bm{\theta}^{\ell})\big{]}\approx\frac{\partial\Delta R_{\epsilon}}{\partial\bm{Z}}(\bm{Z}^{\ell}).$ That is, we need to approximate the gradient flow $\frac{\partial\Delta R_{\epsilon}}{\partial\bm{Z}}$ that locally deforms all (training) features $\{\bm{z}_{i}^{\ell}\}_{i=1}^{N}$ with a continuous mapping $g(\bm{z},\bm{\theta})$ defined on the entire feature space $\bm{z}^{\ell}\in\mathbb{R}^{d}$. Notice that one may interpret the increment (4.1.8) as a discretized version of a continuous differential equation:

Hence the (deep) network so constructed can be interpreted as a certain neural ODE [CRB+18]. Nevertheless, unlike neural ODE where the flow $g$ is chosen to be some generic structures, here our $g(\bm{z},\bm{\theta})$ is to emulate the gradient flow of the rate reduction on the feature set (as shown in Figure 4.1):

and its structure is entirely derived and fully determined from this objective, without any other priors or heuristics.

By inspecting the structure of the gradient (4.1.5), it suggests that a natural candidate for the increment transform $g(\bm{z}^{\ell},\bm{\theta}^{\ell})$ is of the form:

where $\pi_{k}(\bm{z}^{\ell})\in[0,1]$ indicates the probability of $\bm{z}^{\ell}$ belonging to the $k$-th class. The increment map parameters $\bm{\theta}^{\ell}$ depend on: First, a set of linear maps represented by $\bm{E}^{\ell}$ and $\{\bm{C}^{\ell}_{k}\}_{k=1}^{K}$ that depend only on statistics of features of the training $\bm{Z}^{\ell}$; Second, the membership $\{\pi_{k}(\bm{z}^{\ell})\}_{k=1}^{K}$ of any feature $\bm{z}^{\ell}$. Notice that on the training samples $\bm{Z}^{\ell}$, for which the memberships $\bm{\Pi}_{k}$ are known, the so defined $g(\bm{z}^{\ell},\bm{\theta})$ gives exactly the values for the gradient $\frac{\partial\Delta R_{\epsilon}}{\partial\bm{Z}}(\bm{Z}^{\ell})$.

Since we only have the membership for the training samples, the function $g(\cdot)$ defined in (4.1.10) can only be evaluated on the training. To extrapolate $g(\cdot)$ to the entire feature space, we need to estimate $\pi_{k}(\bm{z}^{\ell})$ in its second term. In conventional deep learning, this map is typically modeled as a deep network and learned from the training data, say via back propagation. Nevertheless, our goal here is not to learn a precise classifier $\pi_{k}(\bm{z}^{\ell})$ already. Instead, we only need a good enough estimate of the class information in order for $g(\cdot)$ to approximate the gradient $\frac{\partial\Delta R_{\epsilon}}{\partial\bm{Z}}$ well.

From the geometric interpretation of the linear maps $\bm{E}^{\ell}$ and $\bm{C}_{k}^{\ell}$ given by Remark 4.1, the term $\bm{p}_{k}^{\ell}\doteq\bm{C}^{\ell}_{k}\bm{z}^{\ell}$ can be viewed as (approximately) the projection of $\bm{z}^{\ell}$ onto the orthogonal complement of each class $j$. Therefore, $\|\bm{p}_{j}^{\ell}\|_{2}$ is small if $\bm{z}^{\ell}$ is in class $j$ and large otherwise. This motivates us to estimate its membership based on the following softmax function:

Hence, the second term of (4.1.10) can be approximated by this estimated membership:

which is denoted as a nonlinear operator $\bm{\sigma}(\cdot)$ on outputs of the feature $\bm{z}^{\ell}$ through $K$ groups of filters: $[\bm{C}^{\ell}_{1},\dots,\bm{C}^{\ell}_{K}]$. Notice that the nonlinearality arises due to a “soft” assignment of class membership based on the feature responses from those filters.

Overall, combining (4.1.8), (4.1.10), and (4.1.12), the increment feature transform from $\bm{z}^{\ell}$ to $\bm{z}^{\ell+1}$ now becomes

with the nonlinear function $\bm{\sigma}(\cdot)$ defined above and $\bm{\theta}^{\ell}$ collecting all the layer-wise parameters. That is $\bm{\theta}^{\ell}=\left\{\bm{E}^{\ell},\bm{C}^{\ell}_{1},\dots,\bm{C}^{\ell}_{K},\gamma_{k},\lambda\right\}$. Note features at each layer are always “normalized” by projecting onto the unit sphere $\mathbb{S}^{d-1}$, denoted as $\mathcal{P}_{\mathbb{S}^{d-1}}$. The form of increment in (4.1.13) can be illustrated by a diagram in Figure 4.3(a).

Notice that the increment is constructed to emulate the gradient ascent for the rate reduction $\Delta R_{\epsilon}$. Hence by transforming the features iteratively via the above process, we expect the rate reduction to increase, as we will see in the experimental section. This iterative process, once converged say after $L$ iterations, gives the desired feature map $f(\bm{x},\bm{\theta})$ on the input $\bm{x}=\bm{z}^{0}$, precisely in the form of a deep network, in which each layer has the structure shown in Figure 4.3 left:

As this deep network is derived from maximizing the rate reducation, we call it the ReduNet. By comparing the architecture of ReduNet with those of popular empirically designed networks, ResNet and ResNeXt shown in Figure 4.3, the similarity is somewhat uncanny. Conceptually, ReduNet could also be used to justify the popular mixture of experts (MoE) architecture [SMM+17] as each parallel channel, $\bm{C}^{\ell}_{k}$, can be viewed as an “expert” trained for each class of objects.

We summarize the training and evaluation of ReduNet in Algorithm 4.1 and Algorithm 4.2, respectively. Notice that all parameters of the network are explicitly constructed layer by layer in a forward propagation fashion. The construction does not need any back propagation! The so-learned features can be directly used for classification, say via a nearest subspace classifier.

To classify one-dimensional data $\bm{x}=[x(0),x(1),\ldots,x(D-1)]\in\mathbb{R}^{D}$ invariant under shifting, we take $\mathbb{G}$ to be the group of all circular shifts. Each observation $\bm{x}_{i}$ generates a family $\{\bm{x}_{i}\circ\mathfrak{g}\,|\,\mathfrak{g}\in\mathbb{G}\}$ of shifted copies, which are the columns of the circulant matrix $\mathsf{circ}(\bm{x}_{i})\in\mathbb{R}^{D\times D}$ given by

We refer the reader to [KS12] for properties of circulant matrices. For simplicity, let $\bm{Z}^{1}\doteq[\bm{z}_{1}^{1},\dots,\bm{z}_{N}^{1}]=\bm{X}\in\mathbb{R}^{d\times N}$.⁷⁷7Again, to simplify discussion, we assume for now that the initial features $\bm{Z}^{1}$ are $\bm{X}$ themselves hence have the same dimension $d$, i.e., $D=d$. But that does not need to be the case as we will soon see that we need to lift $\bm{X}$ to a higher dimension. Then what happens if we construct the ReduNet from their circulant families $\mathsf{circ}(\bm{Z}^{1})=\left[\mathsf{circ}(\bm{z}_{1}^{1}),\dots,\mathsf{circ}(\bm{z}_{N}^{1})\right]\in\mathbb{R}^{d\times(dN)}$? That is, we want to compress and map all these into the same subspace by the ReduNet.

Notice that now the data covariance matrix:

associated with this family of samples is automatically a (symmetric) circulant matrix. Moreover, because the circulant property is preserved under sums, inverses, and products, the matrices $\bm{E}^{1}$ and $\bm{C}^{1}_{k}$ are also automatically circulant matrices, whose application to a feature vector $\bm{z}\in\mathbb{R}^{d}$ can be implemented using circular convolution “$\circledast$”. Specifically, we have the following proposition.

Not only do the first-layer parameters $\bm{E}^{1}$ and $\bm{C}^{1}_{k}$ of the ReduNet become circulant convolutions but also the next-layer features remain circulant matrices. That is, the incremental feature transform in (4.1.13) applied to all shifted versions of a $\bm{z}^{1}\in\mathbb{R}^{d}$, given by

is a circulant matrix. This implies that there is no need to construct circulant families from the second layer features as we did for the first layer. By denoting

the features at the next level can be written as

Continuing inductively, we see that all matrices $\bm{E}^{\ell}$ and $\bm{C}^{\ell}_{k}$ based on such $\mathsf{circ}(\bm{Z}^{\ell})$ are circulant, and so are all features. By virtue of the properties of the data, ReduNet has taken the form of a convolutional network, with no need to explicitly choose this structure!

There is one problem though: In general, the set of all circular permutations of a vector $\bm{z}$ gives a full-rank matrix. That is, the $d$ “augmented” features associated with each sample (hence each class) typically already span the entire space $\mathbb{R}^{d}$. For instance, all shifted versions of a delta function $\delta(d)$ can generate any other signal as their (dense) weighted superposition. The MCR² objective (3.4.12) will not be able to distinguish classes as different subspaces.

One natural remedy is to improve the separability of the data by “lifting” the original signal to a higher dimensional space, e.g., by taking their responses to multiple, filters $\bm{k}_{1},\ldots,\bm{k}_{C}\in\mathbb{R}^{d}$:

The filters can be pre-designed invariance-promoting filters,⁸⁸8For 1D signals like audio, one may consider the conventional short-time Fourier transform (STFT); for 2D images, one may consider 2D wavelets as in the ScatteringNet [BM13]. or adaptively learned from the data,⁹⁹9For learned filters, one can learn filters as the principal components of samples as in the PCANet [CJG+15] or from convolution dictionary learning [LB19, QLZ19]. or randomly selected as we do in our experiments. This operation lifts each original signal $\bm{x}\in\mathbb{R}^{d}$ to a $C$-channel feature, denoted as $\bar{\bm{z}}\doteq[\bm{z}[1],\ldots,\bm{z}[C]]^{\top}\in\mathbb{R}^{C\times d}$. Then, we may construct the ReduNet on vector representations of $\bar{\bm{z}}$, denoted as $\vec{(}\bar{\bm{z}})\doteq[\bm{z}[1]^{\top},\ldots,\bm{z}[C]^{\top}]\in\mathbb{R}^{dC}$. The associated circulant version $\mathsf{circ}(\bar{\bm{z}})$ and its data covariance matrix, denoted as $\bar{\bm{\Sigma}}(\bar{\bm{z}})$, for all its shifted versions are given as:

where $\mathsf{circ}(\bm{z}[c])\in\mathbb{R}^{d\times d}$ with $c\in[C]$ is the circulant version of the $c$-th channel of the feature $\bar{\bm{z}}$. Then the columns of $\mathsf{circ}(\bar{\bm{z}})$ will only span at most a $d$-dimensional proper subspace in $\mathbb{R}^{dC}$. However, this simple lifting operation (if linear) is not sufficient to render the classes separable yet—features associated with other classes will span the same $d$-dimensional subspace. This reflects a fundamental conflict between invariance and linear (subspace) modeling: one cannot hope for arbitrarily shifted and superposed signals to belong to the same class.

One way of resolving this conflict is to leverage additional structure within each class, in the form of sparsity: signals within each class are not generated as an arbitrary linear superposition of some base atoms (or motifs), but only sparse combinations of them and their shifted versions, as shown in Figure 4.7. More precisely, let $\bm{D}_{k}=[\bm{d}_{k,1},\ldots,\bm{d}_{k,c}]$ denote a matrix with a collection of atoms associated for class $k$, also known as a dictionary, then each signal $\bm{x}$ in this class is sparsely generated as:

for some sparse vector $\bm{z}$. Signals in different classes are then generated by different dictionaries whose atoms (or motifs) are incoherent from one another. Due to incoherence, signals in one class are unlikely to be sparsely represented by atoms in any other class. Hence all signals can be represented as

where $\bar{\bm{z}}$ is sparse.¹⁰¹⁰10Notice that similar sparse representation models have long been proposed and used for classification purposes in applications such a face recognition, demonstrating excellent effectiveness [WYG+09, WWG+12]. Recently, the convolution sparse coding model has been proposed by [PRE17] as a framework for interpreting the structures of deep convolution networks. There is a vast literature on how to learn the most compact and optimal sparsifying dictionaries from sample data, e.g. [LB19, QLZ19] and subsequently solve the inverse problem and compute the associated sparse code $\bm{z}$ or $\bar{\bm{z}}$. Recent studies of [QLZ20, QZL+20] even show that under broad conditions the convolution dictionary learning problem can be solved effectively and efficiently.

Nevertheless, for tasks such as classification, we are not necessarily interested in the precise optimal dictionary nor the precise sparse code for each individual signal. We are mainly interested if collectively the set of sparse codes for each class are adequately separable from those of other classes. Under the assumption of the sparse generative model, if the convolution kernels $\{\bm{k}_{c}\}_{c=1}^{C}$ match well with the “transpose” or “inverse” of the above sparsifying dictionaries $\bm{D}=[\bm{D}_{1},\ldots,\bm{D}_{K}]$, also known as the analysis filters [NDE+13, RE14], signals in one class will only have high responses to a small subset of those filters and low responses to others (due to the incoherence assumption). Nevertheless, in practice, often a sufficiently large number of, say $C$, random filters $\{\bm{k}_{c}\}_{c=1}^{C}$ suffices to ensure that the extracted $C$-channel features

for different classes have different response patterns to different filters hence make different classes separable [CJG+15].

Therefore, in our framework, to a large extent the number of channels (or the width of the network) truly plays the role as the statistical resource whereas the number of layers (the depth of the network) plays the role as the computational resource. The theory of compressive sensing precisely characterizes how many measurements are needed in order to preserve the intrinsic low-dimensional structures (including separability) of the data [WM21].

The multi-channel responses $\bar{\bm{z}}$ should be sparse. So to approximate the sparse code $\bar{\bm{z}}$, we may take an entry-wise sparsity-promoting nonlinear thresholding, say $\bm{\tau}(\cdot)$, on the above filter outputs by setting low (say absolute value below $\epsilon$) or negative responses to be zero:

Figure 4.8 illustrates the basic ideas. One may refer to [RE14] for a more systematical study on the design of the sparsifying thresholding operator. Nevertheless, here we are not so interested in obtaining the best sparse codes as long as the codes are sufficiently separable. Hence the nonlinear operator $\bm{\tau}$ can be simply chosen to be a soft thresholding or a ReLU. These presumably sparse features $\bar{\bm{z}}$ can be assumed to lie on a lower-dimensional (nonlinear) submanifold of $\mathbb{R}^{dC}$, which can be linearized and separated from the other classes by subsequent ReduNet layers, as illustrated later in Figure 4.9.

The ReduNet constructed from circulant version of these multi-channel features $\bar{\bm{Z}}\doteq[\bar{\bm{z}}_{1},\ldots,\bar{\bm{z}}_{N}]\in\mathbb{R}^{C\times d\times N}$, i.e., $\mathsf{circ}(\bar{\bm{Z}})\doteq[\mathsf{circ}(\bar{\bm{z}}_{1}),\dots,\mathsf{circ}(\bar{\bm{z}}_{N})]\in\mathbb{R}^{dC\times dN}$, retains the good invariance properties described above: the linear operators, now denoted as $\bar{\bm{E}}$ and $\bar{\bm{C}}_{k}$, remain block circulant, and represent multi-channel 1D circular convolutions. Specifically, we have the following result.

From Proposition 4.2, shift invariant ReduNet is a deep convolutional network for multi-channel 1D signals by construction. Notice that even if the initial lifting kernels are separated (4.1.26), the matrix inverse in (4.1.27) for computing $\bar{\bm{E}}$ (similarly for $\bar{\bm{C}_{k}}$) introduces “cross talk” among all $C$ channels. Hence, these multi-channel convolutions in general are not depth-wise separable, unlike the Xception nets [Cho17] that were once suggested to simplify multi-channel convolutional neural networks.¹¹¹¹11It remains open what additional structures on the data would lead to depth-wise separable convolutions.

Following the above derivation, we see that in order to find a linear discriminative representation (LDR) for multiple classes of signals/images that is invariant to translation, sparse coding, a multi-layer architecture with multi-channel convolutions, different nonlinear activation, and spectrum computing all become necessary components for achieving the objective effectively and efficiently. Figure 4.9 illustrates the overall process of learning such a representation via invariant rate reduction on the input sparse codes.