Year: 2019

Abstract:
Nonlinearity is crucial to the performance of a deep (neural) network (DN). To date there has been little progress understanding the menagerie of available nonlinearities, but recently progress has been made on understanding the rôle played by piecewise affine and convex nonlinearities like the ReLU and absolute value activation functions and max-pooling. In particular, DN layers constructed from these operations can be interpreted as max-affine spline operators (MASOs) that have an elegant link to vector quantization (VQ) andK-means. While this is good theoretical progress, the entire MASO approach is predicated on the requirement that the nonlinearities be piecewise affine and convex, which precludes important activation functions like the sigmoid, hyperbolic tangent, and softmax. This paper extends the MASO framework to these and an infinitely large class of new nonlinearities by linking deterministic MASOs with probabilistic Gaussian Mixture Models (GMMs). We show that, under a GMM, piecewise affine, convex nonlinearities like ReLU, absolute value, and max-pooling can be interpreted as solutions to certain natural “hard” VQ inference problems, while sigmoid, hyperbolic tangent, and softmax can be interpreted as solutions to corresponding “soft” VQ inference problems. We further extend the framework by hybridizing the hard and soft VQ optimizations to create a β-VQ inference that interpolates between hard, soft, and linear VQ inference. A prime example of a β-VQ DN nonlinearity is the swish nonlinearity, which offers state-of-the-art performance in a range of computer vision tasks but was developed ad hoc by experimentation. Finally, we validate with experiments an important assertion of our theory, namely that DN performance can be significantly improved by enforcing orthogonality in its linear filters.

heading: 1 INTRODUCTION
text: Deep (neural) networks (DNs) have recently come to the fore in a wide range of machine learning tasks, from regression to classification and beyond. A DN is typically constructed by composing a large number of linear/affine transformations interspersed with up/down-sampling operations and simple scalar nonlinearities such as the ReLU, absolute value, sigmoid, hyperbolic tangent, etc. Goodfellow et al. (2016). Scalar nonlinearities are crucial to a DN’s performance. Indeed, without nonlinearity, the entire network would collapse to a simple affine transformation. But to date there has been little progress understanding and unifying the menagerie of nonlinearities, with few reasons to choose one over another other than intuition or experimentation.
Recently, progress has been made on understanding the rôle played by piecewise affine and convex nonlinearities like the ReLU, leaky ReLU, and absolute value activations and downsampling operations like max-, average-, and channel-pooling Balestriero & Baraniuk (2018a;b). In particular, these operations can be interpreted as max-affine spline operators (MASOs) Magnani & Boyd (2009); Hannah & Dunson (2013) that enable a DN to find a locally optimized piecewise affine approximation to the prediction operator given training data. A spline-based prediction is made in
two steps. First, given an input signal x, we determine which region of the spline’s partition of the domain (the input signal space) it falls into. Second, we apply to x the fixed (in this case affine) function that is assigned to that partition region to obtain the prediction ŷ = f(x).
The key result of Balestriero & Baraniuk (2018a;b) is any DN layer constructed from a combination of linear and piecewise affine and convex is a MASO, and hence the entire DN is merely a composition of MASOs.
MASOs have the attractive property that their partition of the signal space (the collection of multidimensional “knots”) is completely determined by their affine parameters (slopes and offsets). This provides an elegant link to vector quantization (VQ) and K-means clustering. That is, during learning, a DN implicitly constructs a hierarchical VQ of the training data that is then used for splinebased prediction.
This is good progress for DNs based on ReLU, absolute value, and max-pooling, but what about DNs based on classical, high-performing nonlinearities that are neither piecewise affine nor convex like the sigmoid, hyperbolic tangent, and softmax or fresh nonlinearities like the swish Ramachandran et al. (2017) that has been shown to outperform others on a range of tasks?
Contributions. In this paper, we address this gap in the DN theory by developing a new framework that unifies a wide range of DN nonlinearities and inspires and supports the development of new ones. The key idea is to leverage the yinyang relationship between deterministic VQ/K-means and probabilistic Gaussian Mixture Models (GMMs) Biernacki et al. (2000). Under a GMM, piecewise affine, convex nonlinearities like ReLU and absolute value can be interpreted as solutions to certain natural hard inference problems, while sigmoid and hyperbolic tangent can be interpreted as solutions to corresponding soft inference problems. We summarize our primary contributions as follows:
Contribution 1: We leverage the well-understood relationship between VQ, K-means, and GMMs to propose the Soft MASO (SMASO) model, a probabilistic GMM that extends the concept of a deterministic MASO DN layer. Under the SMASO model, hard maximum a posteriori (MAP) inference of the VQ parameters corresponds to conventional deterministic MASO DN operations that involve piecewise affine and convex functions, such as fully connected and convolution matrix multiplication; ReLU, leaky-ReLU, and absolute value activation; and max-, average-, and channelpooling. These operations assign the layer’s input signal (feature map) to the VQ partition region corresponding to the closest centroid in terms of the Euclidean distance,
Contribution 2: A hard VQ inference contains no information regarding the confidence of the VQ region selection, which is related to the distance from the input signal to the region boundary. In response, we develop a method for soft MAP inference of the VQ parameters based on the probability that the layer input belongs to a given VQ region. Switching from hard to soft VQ inference recovers several classical and powerful nonlinearities and provides an avenue to derive completely new ones. We illustrate by showing that the soft versions of ReLU and max-pooling are the sigmoid gated linear unit and softmax pooling, respectively. We also find a home for the sigmoid, hyperbolic tangent, and softmax in the framework as a new kind of DN layer where the MASO output is the VQ probability.
Contribution 3: We generalize hard and soft VQ to what we call β-VQ inference, where β ∈ (0, 1) is a free and learnable parameter. This parameter interpolates the VQ from linear (β → 0), to probabilistic SMASO (β = 0.5), to deterministic MASO (β → 1). We show that the β-VQ version of the hard ReLU activation is the swish nonlinearity, which offers state-of-the-art performance in a range of computer vision tasks but was developed ad hoc through experimentation Ramachandran et al. (2017).
Contribution 4: Seen through the MASO lens, current DNs solve a simplistic per-unit (per-neuron), independent VQ optimization problem at each layer. In response, we extend the SMASO GMM to a factorial GMM that that supports jointly optimal VQ across all units in a layer. Since the factorial aspect of the new model would make naı̈ve VQ inference exponentially computationally complex, we develop a simple sufficient condition under which a we can achieve efficient, tractable, jointly optimal VQ inference. The condition is that the linear “filters” feeding into any nonlinearity should be orthogonal. We propose two simple strategies to learn approximately and truly orthogonal weights and show on three different datasets that both offer significant improvements in classification per-
formance. Since orthogonalization can be applied to an arbitrary DN, this result and our theoretical understanding are of independent interest.
This paper is organized as follows. After reviewing the theory of MASOs and VQ for DNs in Section 2, we formulate the GMM-based extension to SMASOs in Section 3. Section 4 develops the hybrid β-VQ inference with a special case study on the swish nonlinearity. Section 5 extends the SMASO to a factorial GMM and shows the power of DN orthogonalization. We wrap up in Section 6 with directions for future research. Proofs of the various results appear in several appendices in the Supplementary Material.

heading: 2 BACKGROUND ON MAX-AFFINE SPLINES AND DEEP NETWORKS
text: We first briefly review max-affine spline operators (MASOs) in the context of understanding the inner workings of DNs Balestriero & Baraniuk (2018a;b). A MASO is an operator S[A,B] : RD → RK that maps an input vector of length D into an output vector of length K by concatenating K independent max-affine splines Magnani & Boyd (2009); Hannah & Dunson (2013), with each spline formed from R piecewise affine and convex mappings. The MASO parameters consist of the “slopes” A ∈ RK×R×D and the “offsets/biases” B ∈ RK×R. See Appendix A for the precise definition. Given the input x ∈ RD and parameters A,B, a MASO produces the output z ∈ RK via
[z]k = [S[A,B](x)]k = maxr=1,...,R (〈[A]k,r,·,x〉+ [B]k,r) , (1)
where [z]k denotes the kth dimension of z. The three subscripts of the slopes tensor [A]k,r,d correspond to output k, partition region r, and input signal index d. The two subscripts of the offsets/biases tensor [B]k,r correspond to output k and partition region r.
An important consequence of (1) is that a MASO is completely determined by its slope and offset parameters without needing to specify the partition of the input space (the “knots” when D = 1). Indeed, solving (1) automatically computes an optimized partition of the input space RD that is equivalent to a vector quantization (VQ) Nasrabadi & King (1988); Gersho & Gray (2012). We can make the VQ aspect explicit by rewriting (1) in terms of the Hard-VQ (HVQ) matrix TH ∈ RK×R. that contains K stacked one-hot row vectors, each with the one-hot position at index [t]k ∈ {1, . . . , R} corresponding to the arg max over r = 1, . . . , R of (1). Given the HVQ matrix, (or equivalently, a region of the input space), the input-output mapping is affine and fully determined by
[z]k = R∑ r=1 [TH ]k,r (〈[A]k,r,·,x〉+ [B]k,r) . (2)
We retrieve (1) from (2) by noting that [t]k = arg maxr=1,...,R(〈[A]k,r,·,x〉+ [B]k,r). The key background result for this paper is that the layers of a very large class of DN are MASOs. Hence, such a DN is a composition of MASOs, where each layer MASO has as input the feature map z(`−1) ∈ RD(`−1) and produces z(`) ∈ RD(`) , with ` corresponding to the layer. Each MASO has thus specific parameters A(`), B(`). Theorem 1. Any DN layer comprising a linear operator (e.g., fully connected or convolution) composed with a convex and piecewise affine operator (such as a ReLU, leaky-ReLU, or absolute value activation; max/average/channel-pooling; maxout; all with or without skip connections) is a MASO Balestriero & Baraniuk (2018a;b).
Appendix A provides the parameters A(`), B(`) for the MASO corresponding to the `th layer of any DN constructed from linear plus piecewise affine and convex components. Given this connection, we will identify z(`−1) above as the input (feature map) to the MASO DN layer and z(`) as the output (feature map). We also identify [z(`)]k in (1) and (2) as the output of the kth unit (aka neuron) of the `th layer. MASOs for higher-dimensional tensor inputs/outputs are easily developed by flattening.

heading: 3 MAX-AFFINE SPLINES MEET GAUSSIAN MIXTURE MODELS
text: The MASO/HVQ connection provides deep insights into how a DN clusters and organizes signals layer by layer in a hierarchical fashion Balestriero & Baraniuk (2018a;b). However, the entire ap-
proach requires that the nonlinearities be piecewise affine and convex, which precludes important activation functions like the sigmoid, hyperbolic tangent, and softmax. The goal of this paper is to extend the MASO analysis framework of Section 2 to these and an infinitely large class of other nonlinearities by linking deterministic MASOs with probabilistic Gaussian Mixture Models (GMMs).
3.1 FROM MASO TO GMM VIA K-MEANS
For now, we focus on a single unit k from layer ` of a MASO DN, which contains both linear and nonlinear operators; we generalize below in Section 5. The key to the MASO mechanism lies in the VQ variables [t(`)]k ∀k, since they fully determine the output via (2). For a special choice of bias, the VQ variable computation is equivalent to the K-means algorithm Balestriero & Baraniuk (2018a;b).
Proposition 1. Given − 12 ∥∥[A(`)] k,r,· ∥∥2 = [B(`)] k,r , the MASO VQ partition corresponds to a K-
means clustering1 with centroids [ A(`) ] k,r,· computed via [t̂
(`)]k =arg min r=1,...,R
∥∥[A(`)] k,r,·−z (`−1) ∥∥2.
For example, consider a layer ` using a ReLU activation function. Unit k of that layer partitions its input space using a K-means model with R(`) = 2 centroids: the origin of the input space and the unit layer parameter [A(`)]k,1,·. The input is mapped to the partition region corresponding to the closest centroid in terms of the Euclidean distance, and the corresponding affine mapping for that region is used to project the input and produce the layer output as in (2).
We now leverage the well-known relationship between K-means and Gaussian Mixture Models (GMMs) Bishop (2006) to GMM-ize the deterministic VQ process of max-affine splines. As we will see, the constraint on the value of [ B(`) ] k,r
in Proposition 1 will be relaxed thanks to the GMM’s ability to work with a nonuniform prior over the regions (in contrast to K-means).
To move from a deterministic MASO model to a probabilistic GMM, we reformulate the HVQ selection variable [t(`)]k as an unobserved categorical variable [t(`)]k ∼ Cat([π(`)]k,·) with parameter [π(`)]k,· ∈ 4R(`) and 4R(`) the simplex of dimension R(`). Armed with this, we define the following generative model for the layer input z(`−1) as a mixture of R(`) Gaussians with mean [A(`)]k,r,· ∈ RD (`−1) and identical isotropic covariance with parameter σ2
z(`−1) = R(`)∑ r=1 1 ([ t(`) ] k = r ) [ A(`) ] k,r,· + , (3)
with ∼ N(0, Iσ2). Note that this GMM generates an independent vector input z(`−1) for every unit k = 1, . . . , D(`) in layer `. For reasons that will become clear below in Section 3.3, we will refer to the GMM model (3) as the Soft MASO (SMASO) model. We develop a joint, factorial model for the entire MASO layer (and not just one unit) in Section 5.

heading: 3.2 HARD VQ INFERENCE
text: Given the GMM (3) and an input z(`−1), we can compute a hard inference of the optimal VQ selection variable [t(`)]k via the maximum a posteriori (MAP) principle
[t̂(`)]k = arg max t=1,...,R(`)
p ( t|z(`−1) ) . (4)
The following result is proved in Appendix E.1.
Theorem 2. Given a GMM with parameters σ2 = 1 and [π(`)]k,t = exp([B(`)]k,r+ 1 2‖[A (`)]k,r,·‖2)∑ r exp([B (`)]k,r+ 1 2‖[A(`)]k,r,·‖2) , t = 1, . . . , R(`), the MAP inference of the latent selection variable [t̂(`)]k given in (4) can be computed via the MASO HVQ (1)
[t̂(`)]k = arg max r=1,...,R(`)
〈 [A(`)]k,t,·, z (`−1) 〉 + [B(`)]k,t, ∀A(`) ∀B(`). (5)
1It would be more accurate to call this R(`)-means clustering in this case.
The optimal HVQ selection matrix is given by [T̂ (`)H ]k,r = 1 ( r = [t̂(`)]k ) .
Note in Theorem 2 that the bias constraint of Proposition 1 (which can be interpreted as imposing a uniform prior [π(`)]k,·) is completely relaxed.
HVQ inference of the selection matrix sheds light on some of the drawbacks that affect any DN employing piecewise affine, convex activation functions. First, during gradient-based learning, the gradient will propagate back only through the activated VQ regions that correspond to the few 1-hot entries in T (`)H . The parameters of other regions will not be updated; this is known as the “dying neurons phenomenon” Trottier et al. (2017); Agarap (2018). Second, the overall MASO mapping is continuous but not differentiable, which leads to unexpected gradient jumps during learning. Third, the HVQ inference contains no information regarding the confidence of the VQ region selection, which is related to the distance of the query point to the region boundary. As we will now see, this extra information can be very useful and gives rise to a range of classical and new activation functions.

heading: 3.3 SOFT VQ INFERENCE
text: We can overcome many of the limitations of HVQ inference in DNs by replacing the 1-hot entries of the HVQ selection matrix with the probability that the layer input belongs to a given VQ region
[T̂ (`) S ]k,r = p ( [t(`)]k = r | z(`−1) ) = exp (〈 [A(`)]k,r,·, z (`−1)〉+ [B(`)]k,r)∑
r exp (〈 [A(`)]k,r,·, z(`−1) 〉 + [B(`)]k,r ) , (6)
which follows from the simple structure of the GMM. This corresponds to a soft inference of the categorical variable [t(`)]k. Note that T (`) S → T (`) H as the noise variance in (3)→ 0. Given the SVQ selection matrix, the MASO output is still computed via (2). The SVQ matrix can be computed indirectly from an entropy-penalized MASO optimization; the following is reproved in Appendix E.2 for completeness.
Proposition 2. The entries of the SVQ selection matrix [T̂ (`)S ]k,· from (6) solve the following entropypenalized maximization, where H(·) is the Shannon entropy2
[T̂ (`) S ]k,· = arg max
t∈4 R
(`) k
R (`) k∑
r=1
[t]r (〈 [A(`)]k,r,·, z (`−1) 〉 + [B(`)]k,r ) +H(t). (7)
Proposition 2, which was first established in Manning & Klein (2003); Mount (2011), unifies HVQ and SVQ in a single optimization problem. The transition from HVQ (5) to SVQ (7) is obtained simply by adding the entropy regularization H(t). Notice that removing the Entropy regularization from (7) leads to the same VQ as (5). We summarize this finding in Table. 1.

heading: 3.4 SOFT VQ MASO NONLINEARITIES
text: Remarkably, switching from HVQ to SVQ MASO inference recovers several classical and powerful nonlinearities and provides an avenue to derive completely new ones. Given a set of MASO parameters A(`), B(`) for calculating the layer-` output of a DN via (1), we can derive two distinctly different DNs: one based on the HVQ inference of (5) and one based on the SVQ inference of (6). The following results are proved in Appendix E.5.
Proposition 3. The MASO parameters A(`), B(`) that induce the ReLU activation under HVQ induce the sigmoid gated linear unit Elfwing et al. (2018) under SVQ.
Proposition 4. The MASO parameters A(`), B(`) that induce the max-pooling nonlinearity under HVQ induce softmax-pooling Boureau et al. (2010) under SVQ.
Appendix C discusses how the GMM and SVQ formulations shed new light on the impact of parameter initialization in DC learning plus how these formulations can be extended further.
2The observant reader will recognize this as the E-step of the GMM’s EM learning algorithm.

heading: 3.5 ADDITIONAL NONLINEARITIES AS SOFT DN LAYERS
text: Changing viewpoint slightly, we can also derive classical nonlinearities like the sigmoid, tanh, and softmax Goodfellow et al. (2016) from the soft inference perspective. Consider a new soft DN layer whose unit output [z(`)]k is not the piecewise affine spline of (2) but rather the probability [z(`)]k = p([t(`)]k = 1|z(`−1)) that the input z(`) falls into each VQ region. The following propositions are proved in Appendix E.6.
Proposition 5. The MASO parameters A(`), B(`) that induce the ReLU activation under HVQ induce the sigmoid activation in the corresponding soft DN layer.3
A similar train of thought recovers the softmax nonlinearity typically used at the DN output for classification problems.
Proposition 6. The MASO parametersA(`), B(`) that induce a fully-connected-pooling layer under HVQ (with output dimensionD(L) equal to the number of classesC) induce the softmax nonlinearity in the corresponding soft DN layer.

heading: 4 HYBRID HARD/SOFT INFERENCE VIA ENTROPY REGULARIZATION
text: Combining (5) and (6) yields a hybrid optimization for a new β-VQ that recovers hard, soft, and linear VQ inference as special cases
[T̂ (`) β ]k = arg max
t∈4 R
(`) k
[β(`)]k
R (`) k∑
r=1
[t]r (〈 [A(`)]k,r,., z (`−1) 〉 + [B]k,r ) + ( 1− [β(`)]k ) H(t), (8)
with the new hyper-parameter [β(`)]k ∈ (0, 1). The β-VQ obtained from the above optimization problem utilizes [β(`)]k to balance the impact of the regularization term (introduced in the SVQ derivation (7)), allowing to recover and interpolate the VQ between linear, soft and hard (see Table. 1). The following is proved in Appendix E.3. Theorem 3. The unique global optimum of (8) is given by
[T̂ (`) β ]k,r =
exp (
[β(`)]k 1−[β(`)]k
(〈 [A(`)]k,r, z (`−1)〉+ [B(`)]k,r))∑R(`) j=1 exp ( [β(`)]k
1−[β(`)]k
(〈 [A(`)]k,j,., z(`−1) 〉 + [B(`)]k,j )) . (9) The β-VQ covers all of the theory developed above as special cases: β = 1 yields HVQ, β = 12 yields SVQ, and β = 0 yields a linear MASO with [T̂ (`)0 ]k,r = 1 R(`)
. See Figure 1 for examples of how the β parameter interacts with three example activation functions. Note also the attractive property that (9) is differentiable with respect to [β(`)]k.
The β-VQ supports the development of new, high-performance DN nonlinearities. For example, the swish activation σswish(u) = σsig([η(`)]ku)u extends the sigmoid gated linear unit Elfwing et al. (2018) with the learnable parameter [η(`)]k Ramachandran et al. (2017). Numerous experimental studies have shown that DNs equipped with a learned swish activation significantly outperform those with more classical activations like ReLU and sigmoid.4
3The tanh activation is obtained similarly by reparametrizing A(`) and B(`); see Appendix E.6. 4Best performance was usually achieved with [η(`)]k ∈ (0, 1) Ramachandran et al. (2017).
Proposition 7. The MASO A(`), B(`) parameters that induce the ReLU nonlinearity under HVQ induce the swish nonlinearity under β-VQ, with [η(`)]k =
[β(`)]k 1−[β(`)]k .
Table 1 summarizes some of the many nonlinearities that are within reach of the β-VQ.

heading: 5 OPTIMAL JOINT VQ INFERENCE VIA ORTHOGONALIZATION
text: The GMM (3) models the impact of only a single layer unit on the layer-` input z(`−1). We can easily extend this model to a factorial model for z(`−1) that enables all D(`) units at layer ` to combine their syntheses:
z(`−1) = D(`)∑ k=1 R(`)∑ r=1 1 ( [t(`)]k = r ) [A(`)]k,r,· + , (10)
with ∼ N(0, Iσ2). This new model is a mixture of R(`) Gaussians with means [A(`)]k,r,· ∈ RD(`−1) and identical isotropic covariances with variance σ2. The factorial aspect of the model means that the number of possible combinations of the t(`) values grow exponentially with the number of units. Hence, inferring the latent variables t(`) quickly becomes intractable.
However, we can break this combinatorial barrier and achieve efficient, tractable VQ inference by constraining the MASO slope parameters A(`) to be orthogonal〈
[A(`)]k,r,., [A (`)]k′,r′,. 〉 = 0 ∀k 6= k′ ∀r, r′. (11)
Orthogonality is achieved in a fully connected layer (multiplication by the dense matrix W (`) composed with activation or pooling) when the rows of W (`) are orthogonal. Orthogonality is achieved in a convolution layer (multiplication by the convolution matrix C(`) composed with activation or pooling) when the rows of C(`) are either non-overlapping or properly apodized; see Appendix E.4 for the details plus the proof of the following result. Theorem 4. If the slope parameters A(`) of a MASO are orthogonal in the sense of (11), then the random variables [t(`)]1|z(`−1), . . . , [t(`)]1|z(`−1) of the model (10) are independent and hence p ( [t(`)]1, . . . , [t (`)]D(`) |z(`−1) ) = ∏D(`) k=1 p ( [t(`)]k|z(`−1) ) .
In an orthogonal, factorial MASO, optimal inference can be performed independently per factor, as opposed to jointly over all of the factors. Orthogonality renders the joint MAP inference of the factorial model’s VQs tractable. The following result is proved in Appendix E.4.
Practically, this not only lowers the computational complexity tremendously but also imparts the benefit of “uncorrelated unit firing,” which has been shown to be advantageous in DNs Srivastava et al. (2014). Beyond the scope of this paper, such an orthogonalization strategy can also be applied to more general factorial models such as factorial GMMs Zemel (1994); Ghahramani (1995) and factorial HMMs Ghahramani & Jordan (1996).
Corollary 1. When the conditions of Theorem 4 are fulfilled, the joint MAP estimate for the VQs of the factorial model (10)
t̂ (`) f = arg max t∈{1,...,R(`)}×···×{1,...,R(`)} p ( t|z(`−1) ) = [ [t̂(`)]1, . . . , [t̂(`)]D(`) ]ᵀ (12)
and thus can be computed with linear complexity in the number of units.
The advantages of orthogonal or near-orthogonal filters have been explored empirically in various settings, from GANs Brock et al. (2016) to RNNs Huang et al. (2017), typically demonstrating improved performance. Table 2 tabulates the results of a simple confirmation experiment with the largeCNN architecture described in Appendix D. We added to the standard cross-entropy loss a term λ ∑ k ∑ k′ 6=k ∑ r,r′〈[A(`)]k,r,·, [A(`)]k′,r′,·〉2 that penalizes non-orthogonality (recall (11)). We did not cross-validate the penalty coefficient λ but instead set it equal to 1. The tabulated results show clearly that favoring orthogonal filters improves accuracy across both different datasets and different learning settings.
Since the orthogonality penalty does not guarantee true orthogonality but simply favors it, we performed one additional experiment where we reparametrized the fully-connected and convolution matrices using the Gram-Schmidt (GS) process Daniel et al. (1976) so that they were truly orthogonal. Thanks to the differentiability of all of the operations involved in the GS process, we can backpropagate the loss to the orthogonalized filters in order to update them in learning. We also used the swish activation, which we showed to be a β-VQ nonlinearity in Section 4. Since the GS process adds significant computational overhead to the learning algorithm, we conducted only one experiment on the largest dataset (CIFAR100). The exactly orthogonalized largeCNN achieved a classification accuracy of 61.2%, which is a major improvement over all of the results in the bottom (CIFAR100) cell of Table 2. This indicates that there are good reasons to try to improve on the simple orthogonality-penalty-based approach.

heading: 6 FUTURE WORK
text: Our development of the SMASO model opens the door to several new research questions. First, we have merely scratched the surface in the exploration of new nonlinear activation functions and pooling operators based on the SVQ and β-VQ. For example, the soft- or β-VQ versions of leakyReLU, absolute value, and other piecewise affine and convex nonlinearities could outperform the new swish nonlinearity. Second, replacing the entropy penalty in the (7) and (8) with a different penalty will create entirely new classes of nonlinearities that inherit the rich analytical properties of MASO DNs. Third, orthogonal DN filters will enable new analysis techniques and DN probing methods, since from a signal processing point of view problems such as denoising, reconstruction, compression have been extensively studied in terms of orthogonal filters. This work was partially supported by NSF grants IIS-17-30574 and IIS-18-38177, AFOSR grant FA9550-18-1-0478, ARO grant W911NF-15-1-0316, ONR grants N00014-17-1-2551 and N00014-18-12571, DARPA grant G001534-7500, and a DOD Vannevar Bush Faculty Fellowship (NSSEFF) grant N00014-18-1-2047.

heading: SUPPLEMENTARY MATERIALS

heading: A BACKGROUND
text: A Deep Network (DN) is an operator fΘ : RD → RC that maps an input signal x ∈ RD to an output prediction y ∈ RC . All current DNs can be written as a composition of L intermediate mappings called layers
fΘ(x) = ( f (L)
θ(L) ◦ · · · ◦ f (1) θ(1)
) (x), (13)
where Θ = { θ(1), . . . , θ(L) } is the collection of the network’s parameters from each layer. The DN layer
at level ` is an operator f (`) θ(`)
that takes as input the vector-valued signal z(`−1)(x) ∈ RD (`−1) and produces
the vector-valued output z(`)(x) ∈ RD (`)
with D(L) = C. The signals z(`)(x), ` > 1 are typically called feature maps an the input is denoted as z(0)(x) = x. For concreteness, we will focus here on processing multichannel images x but adjusting the appropriate dimensionalities can be used to adapt our results. We will use two equivalent representations for the signal and feature maps, one based on tensors and one based on flattened vectors. In the tensor representation, z(`) contains C(`) channels of size ( I(`) × J(`) ) pixels. In the vector
representation, [z(`)(x)]k represents the entry of the kth dimension of the flattened, vector version z(`)(x) of z(`)(x). Hence, D(`) = C(`)I(`)J(`), C(L) = C, I(L) = 1, and J(L) = 1. For conciseness we will often denote z(`)(x) as z(`). When using nonlinearities and pooling which are piecewise affine and convex, the layers and whole DN fall under the analysis of max-affine spline operators (MASOs) developed in Balestriero & Baraniuk (2018a). In this framework, a max-affine spline operator with parametersA(`) ∈ RD (`)×R×D(`−1) and B(`) ∈ RD (`)×R is defined as
z(`) = S [ A(`), B(`) ] (z(`−1)) =  maxr=1,...,R〈[A (`)]1,r,.,z (`−1)〉+ [B(`)]1,r ...
maxr=1,...,R〈[A(`)]K,r,.,z(`−1)〉+ [B(`)]K,r  . (14) Any DN layer made of convex and piecewise affine nonlinearities or pooling can be rewritten exactly as a MASO. Hence, such operators take place of the layer mappings of (13) We first proceed by modifying (14) to highlight the internal inference problem. We first introduce the VQ-matrix T (`) ∈ RD
(`)×R which will be used to make the mapping region specific, as in
A(`)[T (`)] =  ( ∑R r=1[T (`)]1,r[A (`)]1,r,.) T
... ( ∑R r=1[T (`)]K,r[A (`)]K,r,.) T
 , B(`)[T (`)] =  ( ∑R r=1[T (`)]1,r[B (`)]1,r,.) T
... ( ∑R r=1[T (`)]K,r[B (`)]K,r,.) T  , (15) effectively making A(`)[T (`)] a matrix of shape (D(`), D(`−1)) and B(`)[T (`)] a vector of length D(`). Hence the VQ-matrix is used to combined the per region parameters. In a standard MASO, each row of T (`) is a onehot vector at position corresponding to the region in which the input falls into. Due to the one-hot encoding present in T (`) we refer to this inference as a hard-VQ.
Proposition 8. For a MASO, the VQ-matrix is denoted as T (`)H and is obtained via the internal maximization process of (14). It corresponds to the (hard-)VQ of the input. Once computed the output is a simple affine transform of the input as
z(`) = A(`)[T (`) H ]z (`−1) +B(`)[T (`) H ]. (16)
with [T (`)H ]k,r = 1{r=arg maxr=1,...,R〈[A(`)]k,r,.,z(`−1)〉+[B(`)]k,r}.
The VQ matrix T (`)H always belongs to the set of all matrices with different one-hot positions (from 1 to R) for each of the output dimensions k = 1, . . . , D(`). We denote this VQ-matrix space as T(`)H = {[a1, . . . , aD(`) ] T , ak ∈ {e1, . . . , eR}} with er = δr, dim(er) = R.

heading: B ORTHOGONAL FILTERS DETAILS
text: The developed results on orthogonality induce orthogonality of the case of fully-connected layers. For the case on convolutional layer it implies orthogonality as well as non overlapping patches. This is not practical as it considerably reduces the spatial dimensions making very deep network unsuitable. As such we now
propose a brief approximation result. Due to the specificity of the convolution operator we are able to provide a tractable inference coupled with an apodization scheme. To demonstrate this, we first highlight that any input can be represented as a direct sum of its apodized patches. Then, we see that filtering apodized patches with a filter is equivalent to convolving the input with apodized filters. We first need to introduce the patch notation. We define a patch P[z(`−1)](pi, pj) ∈ {1, . . . , I(`)} × {1, . . . , J(`)} as the slice of the input with indices c = 1, . . . ,K(`), i = (all channels) and (i, j) ∈ {pi, . . . , pi + I(`)C } × {pi, . . . , pi + J (`) C }, hence a patch starting at position (pi, pj) and of same shape as the filters.
Apodizing a signal in general corresponds to applying an apodization function (or windowing function)Weisstein (2002) h onto it via an Hadamard product. Let define the 2D apodized functions h : Ω(I
(`) C , J (`) C )→ R + with Ω(I(`)C , J (`) C ) = {1, . . . , I (`) C }×{1, . . . , J (`) C } and where we remind that (I (`) C , J (`) C ) is the spatial shape of the convolutional filters. Given a function h such that ∑ u∈Ω(I(`)
C ,J (`) C
) h(u) = 1 one can
represent an input by summing the apodized patches as in [z(`)]k,i,j = ∑
(pi,pj)∈{i−I (`) C ,...,i}×{j−J(`) C ,...,j}
P[z(`−1)](pi, pj) h. (17)
The above highlights the ability to treat an input via its collection of patches with the condition to apply the defined apodization function. With the above, we can demonstrate how minimizing the per patch reconstruction loss leads to minimizing the overall input modeling
0 ≤ ‖ ∑ i,j (h P[z(`)](i, j)− [W (`)]t(`)(i,j))‖ 2 ≤ ∑ i,j ‖h P[z(`)](i, j)− [W (`)]t(`)(i,j)‖ 2, (18)
which represents the internal modeling of the factorial model applied across filters and patches. As a result, when performing the per position minimization one minimizes an upper bound which ultimately reaches the global minimum as
‖P[z(`−1)](pi, pj)− P[ẑ(`−1)](pi, pj)‖2 → 0 =⇒ ‖z(`−1) − ∑
(pi,pj)
P[z(`−1)](pi, pj)‖2 = 0. (19)
C INTERPRETATION: INITIALIZATION AND INPUT SPACE PARTITIONING
The GMM formulation and related inference also allows interpretation of the internal layer parameters. First we demonstrate how the region prior π(`) is affected by the layer parameters especially at initialization. Then we highlight how our result allows to generalize the input space partitioning results from Balestriero & Baraniuk (2018b;a).
Region Prior. The region prior of the GMM-MASO model [π(`)]k,. (recall Thm. 2) depends on the bias and norm of the layer weight as [π(`)]k,. ∝ e[B (`)]k,r+ 1 2 ‖[A(`)]k,r,.‖2 . We can study how this region prior looks like at initialization. At initialization, common practice uses [B(`)]k,r = 0, ∀k, r and [A(`)]k,r,d ∼ N(0, (v(`))2). This bias initialization leads to a cluster prior probability proportional to the norm of the weights. For example, the case of absolute value leads to E(‖[A(`)]k,1,.‖2) = E(‖[A(`)]k,2,.‖2) and thus uniform prior as E([π(`)]k,.) = (0.5, 0.5)T for any initialization standard deviation v(`). On the other hand, ReLU has
always ‖[A(`)]k,2,.‖2 = 0 and E ( ‖[A(`)]k,r,.‖2 ) = D(`)(v(`))2. If one uses Xavier initialization Glorot &
Bengio (2010) then D(`)(v(`))2 = 1 and we thus have as prior probability [π(`)]k,. ≈ (0.62, 0.38)T . THe latter slightly favors the inactive state of the ReLU and thus sparser activations. In general, the smaller v(`) is, the more the region prior will favor inactive state of the ReLU.
Input Space Partitioning. We now generalize the ability to study the input space partitioning which was before limited to the special case of [B(`)]k,r = − 12‖[A
(`)]k,r,.‖2 (recall Prop. 1). Studying the input space partition is crucial as the MASO property implies that for each input region, an observation is transformed via a simple linear transformation. However, deriving insights on that is the actual partition is cumbersome as analytical formula are impractical and one thus has to probe the input space and record the observed VQ for each point to estimate the input space partitioning. We are now able to derive some clear links between the MASO partition and standard models which will allows much more efficient computation of the input space partitions.
Corollary 2. A MASO with arbitrary parameters [A(`)]k,r,., [B(`)]k,r has an input space partitioning being the same as a GMM with parameters from Thm. 2.
This augments previous study of the MASO input space partitioning only related to k-mean (recall Prop. 1) which required specific bias values.

heading: D DEEP NETWORK TOPOLOGIES AND DATASETS
text: We first present the topologies used in the experiments except for the notation ResNetD-W which is the standard wide ResNet based topology with depth D and width W . We thus have the following network architectures for smallCNN and largeCNN:
largeCNN
Conv2DLayer(layers[-1],96,3,pad=’same’) Conv2DLayer(layers[-1],96,3,pad=’same’) Conv2DLayer(layers[-1],96,3,pad=’same’,stride=2) Conv2DLayer(layers[-1],192,3,pad=’same’) Conv2DLayer(layers[-1],192,3,pad=’same’) Conv2DLayer(layers[-1],192,3,pad=’same’,stride=2) Conv2DLayer(layers[-1],192,3,pad=’valid’) Conv2DLayer(layers[-1],192,1) Conv2DLayer(layers[-1],10,1)

heading: GlobalPoolLayer(layers[-1],2)
text: where the Conv2DLayer(layers[-1],192,3,pad=’valid’) denotes a standard 2D convolution with 192 filters of spatial size (3, 3) and with valid padding (no padding).

heading: E PROOFS

heading: E.1 THEOREM 2
text: Proof. The log-probability of the model corresponds to
[t(`)]k = arg max r 〈[A(`)]k,r,., z(`−1)〉+ [B(`)]k,r
= arg max r 〈[A(`)]k,r,., z(`−1)〉+ [B(`)]k,r +
1 2 ‖[A(`)]k,r,.‖2 − 1 2 ‖[A(`)]k,r,.‖2
= arg max r 〈[A(`)]k,r,., z(`−1)〉+ [B(`)]k,r +
1 2 ‖[A(`)]k,r,.‖2
− log (∑ r e[B (`)]k,r+ 1 2‖[A (`)]k,r,.‖2 ) − 1 2 ‖[A(`)]k,r,.‖2
= arg max r 〈[A(`)]k,r,., z(`−1)〉+ log
( e[B (`)]k,r+ 1 2‖[A (`)]k,r,.‖2 )
− log (∑ r e[B (`)]k,r+ 1 2‖[A (`)]k,r,.‖2 ) − 1 2 ‖[A(`)]k,r,.‖2
= arg max r 〈[A(`)]k,r,., z(`−1)〉+ log
( e[B (`)]k,r+ 1 2‖[A
(`)]k,r,.‖2∑ r e [B(`)]k,r+ 1 2‖[A(`)]k,r,.‖2
) − 1
2 ‖[A(`)]k,r,.‖2
= arg max r log (p(x|r)p(r))− 1 2 ‖z(`−1)‖2
= arg max r
p(x|r)p(r)
We also remind the reader that arg maxr p(z (`−1)|r)p(r) = arg maxr log(p(z(`−1)|r)p(r)). Based on the above it is straightforward to derive (5) from the above.

heading: E.2 ENTROPY REGULARIZED OPTIMIZATION
text: Proof. We are interested into the following optimization problem:
[t(`)∗]k = arg max q[`,k] F (q[`, k],Θ) = arg max q[`,k] Eq[log(p(z (`−1)|[t(`)]k)p([t(`)]k))] +H([t(`)]k)
= arg max u(`)∈4R (∑ r [u(`)]k,r[ −1 2σ2 ‖z(`−1) − µr‖2 + log(πr)]− ∑ r [u(`)]k,r log([u (`)]k,r) ) .
We now use the KKT and Lagrange multiplier to optimize the new loss function (per k) including the equality constraint
L(u) = ∑ r [u]r[ −1 2σ2 ‖z(`−1) − µr‖2 + log(πr)]− ∑ r [u]r log([u (`)]r) + λ( ∑ r [u]r − 1)
Due to the strong duality we can directly optimize the primal and dual problems and solve jointly all the partial derivatives to 0. We thus obtain by denoting Ar := [ −12σ2 ‖z (`−1) − µr‖2 + log(πr)]
∂L
∂[u]p = Ap − log([u]p)− 1 + λ,∀p
∂L ∂λ = ∑ r [u]r − 1
we can now set the derivatives to 0 and see that this leads to [u]p = eAp−1+λ,∀p. We can now sum over p to obtain
[u]p = e Ap−1+λ,∀p =⇒ ∑ p [u]p = ∑ p eAp−1+λ
=⇒ 1 = ∑ p eAp−1+λ
=⇒ 1 = eλ ∑ p eAp−1
=⇒ 0 = λ+ log( ∑ p eAp−1)
which leads to λ = − log( ∑ p e Ap−1). Plugging this back into the above equation we obtain
[u]p =e Ap−1+λ = eAp−1∑ p e Ap−1 = eAp∑ p e Ap

heading: E.3 THEOREM 3
text: For the proof of Theorem 3 please refer to the proof in E.2 by applying the convex combination with coefficients β.

heading: E.4 THEOREM 4
text: Proof. The proof to demonstrate this inference and VQ equality is essentially the same as the one of GMM-MASO (E.1) with addition of the following first step:
‖z(`−1) − D(`)∑ k=1 [W (`)]k,r,.‖2 = ‖z(`−1)‖2 − 2 D(`)∑ k=1 R(`)∑ r=1 [W (`)]k,[rk],. + D(`)∑ k=1 ‖[W (`)]k,[rk],.‖ 2
for any configuration r ∈ {1, . . . , R(`)}D(`) . Using the same results we can re-write the independent joint optimization as multiple independent optimization problems.

heading: E.5 PROPOSITIONS 3 AND 4
text: For Proposition 4 using the developed formula one can extend the following proof for max-pooling.
Proof.
[z(`)(x)]k = 〈eA(`)[k,2],z(`−1)(x)〉+B(`)[k,2]
1 + e〈A(`)[k,2],z(`−1)(x)〉+B(`)[k,2] × (〈A(`)[k, 2], z(`−1)(x)〉+B(`)[k, 2])
= σsigmoid(〈A(`)[k, 2], z(`−1)(x)〉+B(`)[k, 2])(〈A(`)[k, 2], z(`−1)(x)〉+B(`)[k, 2])
= σsigmoid(〈[C(`)]k,., z(`−1)(x)〉+ [b(`)C ]k)× (〈[C (`)]k,., z (`−1)(x)〉+ [b(`)C ]k) (20)
with 1 for the first region exponential as e〈A (`)[k,2],z(`−1)(x)〉+B(`)[k,1] = e0 = 1 and the last line demonstrating the case where ReLU activation and convolution was the internal layer configuration for illustrative purposes.

heading: E.6 PROPOSITIONS 5 AND 6
text: Proof.
p([t(`)]k = 1|z(`−1)) = p(z(`−1)|[t(`)]k = 1)p([t(`)]k = 1)
p(z(`−1))
= p(z(`−1)|[t(`)]k = 1)p([t(`)]k = 1)
p(z(`−1)|[t(`)]k = 0)p([t(`)]k = 0) + p(z(`−1)|[t(`)]k = 1)p([t(`)]k = 1)
= e(−‖z (`−1)−[C(`)]k,.‖2 2 ) e( 12‖[C (`)]k,.‖2+[B(`)]k,1) 1+e( 12‖[C(`)]k,.‖2+[B(`)]k,1)
e(−‖z (`−1)‖2 2 )
1 1+e( 12‖[C(`)]k,.‖2+[B(`)]k,1)
+ e(−‖z (`−1)−[C(`)]k,.‖2 2 ) e( 12‖[C(`)]k,.‖2+[B(`)]k,1)
1+e( 12‖[C(`)]k,.‖2+[B(`)]k,1)
= e(−‖z (`−1)−[C(`)]k,.‖2 2 )e( 1 2‖[C (`)]k,.‖2 + [B(`)]k,1)
e(−‖z (`−1)‖2 2 ) + e(− ‖z(`−1)−[C(`)]k,.‖2 2 )e( 1 2‖[C(`)]k,.‖2 + [B(`)]k,1)
= e(〈z(`−1), [C(`)]k,.〉+ [B(`)]k,1)
1 + e(〈z(`−1), [C(`)]k,.〉+ [B(`)]k,1) = σsigmoid([u (`)]k).
While this is direct for sigmoid DNs, the use of hyperbolic tangent requires to reparametrize the current and following layer weights and biases to represent the shifting scaling as in C(`) := 2C(`) and C(`+1) := 2C(`+1), b(`+1)C := b (`+1) C − 1 with C replaced by W for fully connected operators.

Review summaries:
At the core of this paper is the insight from [1] that a neural network layer constructed from a combination of linear, piecewise affine and convex operators can be interpreted as a max-affine spline operator (MASO). MASOs are directly connected to vector quantization (VQ) and K-means clustering, which means that a deep network implicitly constructs a hierarchical clustering of the training data during learning. This paper now substitutes VQ with probabilistic clustering models (GMMs) and extends the MASO interpretation of a wider range of possible operations in deep neural networks (sigmoidal activation functions, etc.).

This work extends the applicability of the spline theory of deep networks explored in previous works of Balestriero/ Baraniuk. The previous works setup DNs as layer-wise max-affine spline operators (MASOs) and recovers several non-linearities practically used as special cases of these MASOs. The previous works already recover RELU variants and some downsampling operators that the current submission characterizes as "hard" quantization.