1 Introduction↩︎

Neural networks are the backbone of contemporary applications in machine learning and related fields, having huge influence and significance both in theory and practice. Among the most important and desirable attributes of a trained network are robustness and sparsity. Robustness, is often defined as stability to adversarial perturbations, such as in supervised classification methods. The apparent brittleness of neural networks to adversarial attacks in this context has been considered in the literature for some time, see e.g., [1]–[3] and references therein.

A fundamental question in this regard is how to measure robustness, or more importantly, how to encourage it. One prominent approach supported by theory and practice [4], [5], is to use the Lipschitz constant of the network function to quantize robustness, and regularization to encourage it.

This approach is also supported theoretically with generalization bounds in terms of the layer-wise product of spectral norms [6], [7], which particularly upper-bounds the Lipschitz constant. However, a recent empirical study [8] has found in practice a negative correlation of this measure with generalization. This casts doubts on its usefulness and signals the fact that it is a rather loose upper bound for the Lipschitz constant [9].

Current methods that compute upper bounds on the Lipschitz constant of neural networks can be roughly classified into two classes: (i) the class of product bounds, comprising all upper bounds obtained by the multiplication of layer-wise matrix norms; and, (ii) the class of convex-optimization-based bounds, which addresses the network as a whole entity [4], [9], [10].

A trade-off between computational complexity and quality of the upper bound seems apparent. An ideal bound would achieve a balance between both properties: it should provide a good estimate of the constant while being fast and easy to minimize with iterative first-order algorithms.

Recently, the path-norm of the network [11] has emerged as a complexity measure that is highly-correlated with generalization [8]. Thus, its use as a regularizer holds an increasing interest for researchers in the field.

Despite existing generalization bounds [11], our understanding of the optimization aspects of the path-norm-regularized objective is lacking. [8] refrained from using automatic-differentiation methods in this case because, as they argue, the optimization could fail, thus providing no conclusion about its qualities.

It is then natural to ask: how do we properly optimize the path-norm-regularized objective with theoretical guarantees? What conclusions can we draw about the robustness and sparsity of path-norm-regularized networks? We focus on the 1-path-norm and provide partial answers to those questions, further advancing our understanding of this measure. Let us summarize our main contributions:

Optimization. We show a striking property of the 1-path-norm, that makes it a strong candidate for explicit regularization: despite its non-convexity, it admits an efficient proximal mapping (Algorithm 3). This allows the use of proximal-gradient type methods which are, as of now, the only first-order optimization algorithms to provide guarantees of convergence for composite non-smooth and non-convex problems [12].

Indeed, automatic differentiation modules of popular deep learning frameworks like PyTorch [13] or TensorFlow [14] may not compute the correct gradient for compositions of non-smooth functions, at points where these are differentiable [15], [16]. Our proposed optimization algorithm avoids such issue altogether by using differentiable activation functions like ELU [17] and our novel proximal mapping of the 1-path-norm.

Upper bounds. We show that the 1-path-norm [11] achieves a sweet spot in the computation-quality trade-off observed among upper bounds of the Lipschitz constant: it has a simple closed formula in terms of the weights of the network, and it provides an upper bound on the \((\ell_\infty, \ell_1)\)-Lipschitz constant (cf., 1), which is always better than the product bound.

Sparsity. Neural network regularization schemes promoting sparsity in a principled way are of great interest in the growing field of compression in Deep Learning [18], [19].

Our analysis provides a formula (cf. Lemma 4) for choosing the strength of the regularization, which enforces a desired bound on the sparsity level of the iterates generated by the proximal gradient method. This is a suprising, yet intuitive, result, as the sparsity-inducing properties of non-smooth regularizers have been observed before in convex optimization and signal processing literature, see e.g., [20], [21].

Experiments. In 7, we present numerical evidence that our approach (i) converges faster and to lower values of the objective function, compared to plain SGD; (ii) generates sparse iterates; and, (iii) the magnitude of the regularization parameter of the 1-path-norm allows a better accuracy-robustness trade-off than the common \(\ell_1\) regularization or constraints on layer-wise matrix norms.

2 Problem Setup↩︎

We consider the so-called shallow neural networks with \(n\) hidden neurons and \(p\) outputs \(h:\mathbb{R}^m \rightarrow \mathbb{R}^p\) given by \[\label{eq:nn} h_{V,W}(x) = V^T \sigma(Wx),\tag{1}\] where \(V \in \mathbb{R}^{n \times p}, W\in \mathbb{R}^{n\times m}\) and \(\sigma:\mathbb{R} \rightarrow \mathbb{R}\) is some differentiable activation function with derivative globally bounded between zero and one. This condition is satisfied, for example, by the ELU or softplus activation functions. To control the robustness of the network to perturbations of its input \(x\), we want to regularize training using its Lipschitz constant as a function of the weights \(V\) and \(W\).

To properly define this constant, we utilize the \(\ell_\infty\)-norm for the input space, and the \(\ell_1\)-norm for the output space. Exact computation of such constant is a hard task. A simple and easily computable upper bound can be derived by the product of the layer-wise Lipschitz constants, however, it can be quite loose.

We derive an improved upper bound which is still easy to compute. In the following, we denote with \(\|W\|_\infty\) the operator norm of a matrix \(W\) with respect to the \(\ell_\infty\) norm for both input and output space; it is equal to the maximum \(\ell_1\)-norm of its rows. We denote with \(\|V\|_{\infty, 1}\) the operator norm of the matrix \(V\) with respect to the \(\ell_\infty\) norm in input space and \(\ell_1\)-norm in output space; it is equal to the sum of the \(\ell_1\) norm of its columns.

The proof is provided in appendix 9. The term in the middle of inequality ?? belongs to the family of path-norms, introduced in [11]. Throughout, we refer to it as the 1-path-norm.

Notice that although the path-norm and layer wise product bounds can be equal, this only happens in the following worst case: For the weight matrix in the first layer, the 1-norms of every row are equal. Thus, in practice the bounds can differ drastically.

We now turn to the task of minimizing an empirical risk functional regularized by the improved upper bound on the Lipschitz constant given in (?? ): \[\label{eq:objective} \min_{V,W} \mathbb{E}_{(x, y)}[\ell(h_{V, W}(x), y)] + \lambda \sum_{i=1}^n\sum_{j=1}^m\sum_{k=1}^p |W_{ij}V_{ik}|\tag{2}\] The objective function in problem 2 is composed of an expectation of a nonconvex smooth loss, and a nonconvex nonsmooth regularizer, meaning that it is essentially a composite problem (cf. [22]). That is, the objective function 2 can be cast as use these notation hereafter) \[\label{eq:objective95fn} \min_{V, W} \mathcal{F} (V, W) \equiv f(V, W) + \lambda g(V, W),\tag{3}\] where \(f\) is a nonconvex continuously differentiable function, and \(g(V, W)= \sum_{i=1}^n\sum_{j=1}^m\sum_{k=1}^p |W_{ij}V_{ik}|\) is a continuous, nonconvex, nonsmooth, function. We assume that the objective function is bounded below, i.e., \(\displaystyle\mathop{\rm inf}\mathcal{F} :=\mathcal{F}_*>0\).

A natural choice for a scheme to obtain critical points for (3 ) is the proximal-gradient framework. However, for a nonconvex \(g\), solving the proximal gradient problem is a hard problem in general. In Section 4 we develop a method that computes the proximal gradient with respect to \(g\) efficiently.

To streamline our approach and techniques in a compact and user-friendly manner, we will illustrate the majority of our results and proofs via the particular single-output scenario in which \(h\) and \(g\) are reduced to \[h_{v, W}(x) = v^T \sigma(Wx), \;g (v,W) \equiv \lambda \| \mathop{\mathrm{vec}}(\mathop{\mathrm{Diag}}(v)W) \|_1.\] The multi-output case follows from the same techniques and insights, however, requires more tedious computations and arguments, on which we elaborate in Section 5, and detail in the appendix.

3 The Prox-Grad Method↩︎

Assume that \(f\) has a Lipschitz continuous gradient with Lipschitz constant \(L>0\), that is \[\| \nabla f(z) - \nabla f(u) \| \leq L \| z - u\|, \qquad \forall z,u \in\mathbb{R}^n.\] The prox-grad method is described by Algorithm 1; since \(g\) is nonconvex, the prox in 4 can be a set of solutions.

Theoretical guarantees for the prox-grad method with respect to a nonconvex regularizer were established by [12] (for a more general prox-grad type scheme).

4 Computing the Proximal Mapping↩︎

Throughout this section we assume the single-output setting. The path-norm regularizer we propose is a nonconvex nonsmooth function, suggesting that the prox-grad scheme in Algorithm 1 is intractable.

In this section we will not only prove that in fact it is tractable in the single output case, but that it can also be implemented efficiently with complexity of \(O(m\log(m))\); we prove the stated in detail in Section 11, and provide here a concise version.

Denote the given pair \((x, Y)\) by \(z\). The proximal mapping with respect to \(\lambda g\) at \(z\) is defined as \[\label{eq:proximal} \mathop{\mathrm{prox}}_{\lambda g}(z) = \displaystyle\mathop{\rm argmin}_{u} \lambda g(u) + \frac{1}{2} \|{\mathop{\mathrm{vec}}(u - z)}\|_2^2.\tag{4}\] By the choice of \(g\), the objective function in 4 is coercive and lower bounded, implying that there exists an optimal solution (cf. [24]).

Additionally, we have that 4 is separable with respect to the \(i\)-th entry of the vector \(v\) and the \(i\)-th row of the matrix \(W\), meaning that problem 4 can be solved in a distributed manner by applying the same solution procedure coordinate-wise for \(v\) and row-wise for \(W\). In light of this, let us consider the \(i\)-th row related problem \[\label{eq:7b} \min_{v, w\in\mathbb{R}\times\mathbb{R}^m} \frac{1}{2} (v-x)^2 + \frac{1}{2} \sum_{j=1}^m (w_j - y_j)^2 + \lambda |v| \sum_{j=1}^m |w_j|.\tag{5}\] The signs of the elements of the decision variables in (5 ) are determined by the signs of \((x, y)\), and consequently, the problem in (5 ) is equivalent to \[\label{eq:cw95prox95prob} \min_{v, w\in\mathbb{R}_+\times\mathbb{R}_+^m} \frac{1}{2} (v-|x|)^2 + \frac{1}{2} \sum_{j=1}^m (w_j - |y_j|)^2 + \lambda v \sum_{j=1}^m w_j.\tag{6}\]

Denote \[h_\lambda(v,w;x,y)= \frac{1}{2} (v-|x|)^2 + \frac{1}{2} \sum_{j=1}^m (w_j - |y_j|)^2 + \lambda v \sum_{j=1}^m w_j.\] Although \(h_\lambda\) is nonconvex, we will show that a global optimum to 6 can be obtained efficiently by utilizing several tools, the first being the first-order optimality conditions of 6 (cf. [24]) given below.

A key insight following Lemma 3 is that: the elements of any solution to 6 , satisfy a monotonic relation in magnitude, correlated with the magnitude of the elements of \(y\); this is formulated by the next corollary.

To supplement the results above, we now show that we can actually upper-bound the sparsity level of the prox-grad output by adjusting the value of \(\lambda\).

Moreover, the function \(h_\lambda\) is monotonically decreasing in the sparsity level, which implies that instead of exhaustively checking the value of \(h_\lambda\) for any sparsity level, we can employ a binary search. Denote for any \(s\in \{0, \ldots, m\}\) the \(m+1\) possible solutions: \[\begin{align} v^{(s)} &= \dfrac{1}{1-s \lambda^2}\left(|x| - \lambda\sum_{j=1}^s|y_j|\right) \\ w^{(s)}_j &= |y_j| - \lambda v^{(s)} \;\text{ for } j\in [s], \;\text{ and } w^{(s)}_j = 0 \;\text{otherwise}. \end{align}\]

Lemma 5 follows from algebraic considerations, and thus its proof is deferred to Section 11. Its substantial implication is the following.

Note that since, by definition, the \(s\) first entries of the vector \(w^{(s)}\) are ordered in decreasing order, the constrained \(w^{(s)}_s > 0\) ensures that the full vector \(w^{(s)}\) has exactly \(s\) nonzero entries, which are all strictly positive.

The final ingredient required for designing an efficient algorithm is the following monotone property of the feasibility criterion in problem ?? :

This property, whose proof is also deferred to Section 11, implies that the optimal sparsity parameter \(s^*\) can be efficiently found using a binary search approach.

We conclude this section by combining all the ingredients above to develop Algorithm 2, and to prove that it yields a solution to 4 .

4.0.0.1 Time complexity of Algorithm 2

In the worst case where \(m \leq \lambda^{-2}\), the number of searches for finding \(s^*\) is at most \(\log_2(m)\). Each search requires to compute \(v^{(s)}\), and in particular \(\sum_{j=1}^s |y_j|\), as well as \(w_j^{(s)}\), \(j=1,\ldots,s\), each taking \(\mathcal{O}(s)\) steps. Thus, the overall loop complexity is \(\mathcal{O}(m)\).

Moreover, this algorithm assumes that the input vector \(y\) is already sorted in decreasing magnitude order. This can easily be achieved by a sorting procedure in time \(\mathcal{O}(m \log m)\).

5 Multi-Output↩︎

The efficient computation of the robust-sparse proximal mapping we derived for the single-output scenario will now be generalized to the multi-output case. Although we use similar arguments and insights, the analysis is much more complicated and requires more delicate and advanced treatment. Due to the tedious computations that accompany the analysis, the proofs are deferred to appendix 12.

When the network has multiple-output, the proximal operator \(\text{prox}_{\lambda g}(X,Y)\) can be written as the solution set of \[\label{eq:multi-prox} \begin{align} \min_{V,W} \|V - X\|_F + \|W - Y\|_F + 2\lambda \sum_{i=1}^n\sum_{j=1}^m\sum_{k=1}^p |W_{ij}V_{ik}|, \end{align}\tag{7}\] where \(V\in \mathbb{R}^{n\times p}\) and \(W\in \mathbb{R}^{n\times m}\). As in the single-output case, we observe that the proximal mapping 7 is separable with respect to the \(i\)-th rows of the matrices \(V\) and \(W\), and that the signs of the decision variables are determined by the signs of \((X,Y)\). Therefore, it is enough to consider the problem related to the \(i\)-th row of \(V\), denoted as \(x\), and \(i\)-th row of \(W\), denoted as \(y\), i.e., \[\label{eq:5} \begin{align} \min_{v, w \in \mathbb{R}_+^p \times \mathbb{R}_+^m} h_\lambda(v,w;x,y), \end{align}\tag{8}\] where we redefine \(h_\lambda(v,w;x,y)\) to include the multi-output case: \(h_\lambda(v,w;x,y) = \frac{1}{2} \sum_{k=1}^p (v_k-|x_k|)^2 + \frac{1}{2} \sum_{j=1}^m (w_j - |y_j|)^2 + \lambda \sum_{k=1}^p v_k \sum_{j=1}^m w_j.\) To improve readability, we will abuse notation and just write \(h_\lambda(v,w)\), assuming that \((x,y)\) are understood from context.

Using the same observations we exploited to enumerated all stationary points of the proximal mapping in the single-output setup, we can identify the stationary points depending on the number of non zero elements of \(v\) and \(w\).

From the two first points in Lemma 7, the argument in Remark 5 is also valid in the multi-output case, and so we assume hereafter that the input vectors \(x, y\) are sorted in decreasing magnitude order.

Using the second order stationary conditions, we can generalize our sparsity bound in the single-output scenario, given in Lemma 4, to an upper bound on the product of the sparsities of the solutions based on the value of \(\lambda\); indeed, \(s_v = 1\) yields the bound in Lemma 4.

A possible algorithm for computing this proximal mapping would thus be to compute the value of \(h_\lambda \left(v^{(s_v, s_w)},w^{(s_v, s_w)}\right)\) for each pair of sparsities \((s_v, s_w) \in \{0,\ldots,p\} \times \{0,\ldots,m\}\) satisfying \(s_v s_w \leq \lambda^{-2}\) and return the pair achieving the smallest value.

However, such an approach would be computationally inefficient. In order to avoid computing the value of \(h_\lambda\) at each pair, we show the following monotonicity property of \(h_\lambda\) in the sparsity levels, which generalizes the same property in the single-output case.

Moreover, the feasibility criterion \(v \geq 0, w \geq 0\) also has a monotonic property:

To properly address the complications arising from handling two intertwining sparsity levels at the same time, we introduce the notion of maximal feasibility boundary (MFB) which acts a frontier of possible sparsity levels.

The efficient computation of the multi-output robust-sparse proximal mapping is based on the fact that we only need to compute the value of \(h_\lambda\) for sparsity levels that are at the frontier of the MFB. This allows us to find the optimal sparsity in time \(\mathcal{O}(p + m)\), improving upon the \(\mathcal{O}(pm)\) complexity of the exhaustive search. Algorithm 3 implements the above by employing a binary search type procedure defined in Algorithm 5 to calculate the MFB.

5.0.0.1 Time complexity of Algorithm 3

It is easy to see that the maximal feasibility boundary contains at most \(\min(m,p)\) pairs, and Algorithm 5 finds them all in time \(\mathcal{O}(m+p)\). Then, for each such pair \((s_v, s_w)\), we must compute \(v^{(s_v, s_w)}\) and \(w^{(s_v, s_w)}\) and \(h_\lambda(v^{(s_v, s_w)}, w^{(s_v, s_w)})\), which takes time \(\mathcal{O}(m+p)\). The total complexity of Algorithm 3 is thus \(\mathcal{O}(\min(m,p)(m+p))\). In most practical application, the output layer size \(p\) can be considered \(\mathcal{O}(1)\), so that the complexity of computing this proximal mapping is comparable to the complexity of computing one stochastic gradient.

6 Related Work↩︎

The path regularization approach to train neural networks can be traced back to the seminal paper by [11], who introduced the \(p\)-path-norm as a heuristic proxy to control the capacity of the network.

In this paper (cf. 1), we took a step forward by moving from heuristic explanations to rigorous arguments by establishing a new connection between the \(1\)-path-norm and the Lipschitz constant of the network. This result also reads as a relation between the 1-path-norm and the product bound, of which variants have been found to be useful in deriving generalization bounds [6].

Generalization bounds in terms of the \(p\)-path-norm were also derived in [11], but the question of how to methodologically exploit these as regularizers remains open; our algorithmic contribution is a first step in this direction. Additionally, issues regarding optimization with path-norm regularization were reported by [27], which examined a conditional gradient method in the context of path-norm regularization.

A growing collection of works have focused on the task of network compression, doing so via sparsity-inducing regularizers [28]–[31]. They have achieved a great level of success by setting the regularization term in an ad-hoc manner. In contrast, we follow a principled regularization approach with theoretical properties of generalization and robustness, and as a consequence, we are able to quantify the sparsity of the resulting networks.

Moreover, the aforementioned works only use convex regularizers for which efficient proximal mappings are available (see [20] and references therein). We tackle the much harder non-convex regularization task, and derive a new method to compute the proximal mapping in this case. The merits of non-convex non-smooth regularization, and difficulties regarding their optimization, have been extensively studied in the imaging and signal sciences, see e.g., [32] and the recent survey [33].

Layer-wise constraints or regularization with matrix-norms, which are also motivated by the product bound, have been used for robust training [5], [34] and generative models [7]. These focus on robustness with respect to the \(\ell_2\)-norm, which requires a careful handling of operations on the singular values of the weight matrices, and does not have the extra benefit of inducing sparsity.

In 7 we compare to this class of methods for the \(\ell_\infty\)-norm case, in which a simple rescaling of the rows in the weight matrices yields a numerically stable procedure [35], [36].

7 Experiments↩︎

We empirically evaluate shallow neural networks trained by regularized empirical risk minimization 3 using cross-entropy loss. In terms of the weight matrices \(V\) and \(W\) of the network 1 , the following regularizers are considered:

\(\ell_1\) regularization. We penalize the \(\ell_1\)-norm of the parameters of the network, i.e., \(g(V, W) = \| \mathop{\mathrm{vec}}(V)\|_1 + \| \mathop{\mathrm{vec}}(W) \|_1\) in the objective function 3 .

1-path-norm regularization. We set \(g(V, W)\) as \(\sum_{i=1}^n\sum_{j=1}^m\sum_{k=1}^p |W_{ij}V_{ki}|\) in the objective function 3 .

Layer-wise regularization (Parseval Networks). we minimize the cross-entropy loss with a hard constrain on the \(\ell_\infty\)-operator-norm of the weight matrices i.e., \(\|W\|_\infty \leq \lambda^{-1}\) and \(\|V\|_\infty \leq \lambda^{-1}\), as described by [5]. The projection on such set is achieved by projecting each row of the matrices onto an \(\ell_1\)-ball using efficient algorithms [35], [36].

Remark. We will refer (incorrectly) to the training loop defined by PyTorch’s SGD optimizer as Stochastic gradient descent (SGD) (see the discussion in 1).

Experimental setup. Our benchmarks are the MNIST [37], Fashion-MNIST [38] and Kuzushiji-MNIST [39]. For a wide range of learning rates, number of hidden neurons and regularization parameters \(\lambda\), we train networks with SGD and Proximal-SGD (with constant learning rate). We do so for 20 epochs and with batch size set to 100. For each combination of parameters we train 6 networks with the default random initialization. Details and further experiments are reported in appendix 13.

7.1 Convergence of SGD vs Proximal-SGD↩︎

Due to the non-differentiability of the \(\ell_1\)- and path-norm regularizers, we expect Proximal-SGD to converge faster, and to lower values of the regularized loss, when compared to SGD. This is examined in Figure 4, where we plot the value of the loss function across iterations. For both SGD and Proximal-SGD, the loss function decays rapidly in the first few epochs. We then enter a second regime where SGD suffers from slow convergence, whereas Proximal-SGD continues to reduce the loss at a fast rate. At the end of the 20-epochs, Proximal-SGD consistently achieves a lower value of the loss.

An advantage of Proximal-SGD over plain SGD is that the proximal mappings of both the \(\ell_1\)- and path-norm regularizers can set many weights to exactly zero. In [fig:nnz95vs95iter] we plot the average error and robust test error obtained, as a fuction of the sparsity of the network. Compared to \(\ell_1\) regularization, the sparsity pattern induced by the 1-path-norm correlates with the robustness to a higher degree. As a drawback, it appears that in more difficult datasets like KMNIST, the 1-path-norm struggles to obtain good accuracy and sparsity simultaneously.

7.2 The robustness-accuracy trade-off↩︎

The relation between the Lipschitz constant of a network and its robustness to adversarial perturbations has been extensively studied in the literature. In 1 we have shown that the 1-path-norm of a single-output network is a tighter upper bound of its Lipschitz constant, compared to the corresponding product bound.

To the best of our knowledge, the \(\ell_1\)-norm regularizer only provides an upper boud on the already loose product bound [11], which makes it less attractive as a regularizer, despite its sparsity-inducing properties. Hence, the 1-path-norm regularizer is, in theory, a better proxy for robustness than the other regularization schemes.

[fig:error95vs95robust] shows the misclassification error on clean and adversarial examples as a function of \(\lambda\), and corresponds to the learning rate minimizing the error on clean samples. The adversarial perturbations were obtained by PGD [3].

Any training procedure which promotes robustness of a classifier may decrease its accuracy, and this effect is consistently observed in practice [40]. Hence, the merits of a regularizer should be measured by how efficiently it can trade-off accuracy for robustness. We observe that for all three regularization schemes, there exists choices of \(\lambda\) that attain the best possible error on clean samples.

On the other hand, the error obtained by the \(\ell_1\) regularization degrades significantly. The layer-wise and 1-path-norm regularization achieve a noticeably low error on adversarial examples. Comparing the latter schemes, the 1-path-norm regularization shows only a slight advantage over the layer-wise methods, which merits further investigation.

8 Future Work: Multilayer Extension↩︎

A natural extension of our approach is to apply path regularization to multi-layered networks. Since the number of paths is potentially huge, this scenario requires more sophisticated treatment, and hence left for future research. Nonetheless, a trivial extension of our approach is to divide a multi-layered network into pairs of consecutive layers, and apply our method in a sensible manner. We now describe this approach, to complement the theory.

Precisely, assume that the network contains an even number of layers. For some lists of matrices \(\mathbf{V}=(V^1,\ldots, V^k)\) and \(\mathbf{W}=(W^1, \ldots, W^k)\) of appropriate sizes, the network can be written as a composition of activation functions and shallow networks \[h_{\mathbf{V}, \mathbf{W}}:=h_{W^l, V^l} \circ \sigma \circ h_{W^{l-1}, V^{l-1}} \circ \ldots \circ \sigma \circ h_{V^1, W^1}\]

We build the regularized objective as \[\label{eq:objective95multilayer} \min_{\mathbf{V}, \mathbf{W}} \mathbb{E}_{(x, y)}[\ell(h_{\mathbf{V}, \mathbf{W}}(x), y)] + \lambda \sum_{l=1}^L P_1(V^l, W^l).\tag{9}\] where we have introduced the shorthand \[P_1(V^l, W^l):= \sum_{i=1}^{n_l} \sum_{j=1}^{m_l} \sum_{k=1}^{p_l} |W^l_{ij} V^l_{ik}|\] For the 1-path-norm of the \(l\)-th subnetwork \(h_{W^{l}, V^l}\).

Because the nonsmooth nonconvex regularizer in 9 is separable in the variables \(\{(V_i, W_i): i=1,\ldots,k\}\), its proximal mapping is indeed nothing but our multioutput prox algorithm (section 5), applied independently to each shallow subnetwork component. Thus, the proximal gradient methods can be applied efficiently to optimize 9 .

Acknowledgements↩︎

This work is funded (in part) through a PhD fellowship of the Swiss Data Science Center, a joint venture between EPFL and ETH Zurich. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement n. 725594 - time-data). This work was supported by the Swiss National Science Foundation (SNSF) under grant number 200021_178865 / 1.

9 Proof of Theorem 1↩︎

We will first prove a particular case of 1, the single-output case (\(p=1\)).

First, note that because the output space is \(\mathbb{R}\), the \(\ell_1\)-norm is just the absolute value of the output. In this case the Lipschitz constant of the single-output function \(h\) is equal to the supremum of the \(\ell_1\)-norm of its gradient, over its domain (c.f., [9]).

Proof of 1. We now proceed with the general case where \(V \in \mathbb{R}^{n \times p}\), \(W \in \mathbb{R}^{n \times m}\) and \(h_{V,W}(x)=V^T \sigma(Wx)\).

10 Proof of Theorem 2↩︎

In this section we prove the theoretical guarantees stated in Theorem 2 of the prox-grad method described by Algorithm 1. The first and second parts of Theorem 2 follow immediately from the results establish by [12]. Part two in Theorem 2 states that Algorithm 1 is globally convergent under the celebrated Kurdyka–Lojasiewicz (KL) property [41]. The broad classes of semi-algebraic and subanalytic functions, widely used in optimization, satisfy the KL property (see e.g. [12]), and in particular, most convex functions encountered in finite dimensional applications satisfy it (see [12]). We refer the reader to the works [12], [41], [42], in particular to [12] for additional information and results.

For Part three we require the sufficient decrease property stated next.

11 Single output proximal map computation↩︎

This section provides the theoretical background and the required intermediate results to prove Theorem 3.

11.1 Moving to an Equivalent Easier Problem↩︎

We are interested in minimizing the nonconvex twice continuously differentiable function \[\label{eq:7} \min_{v, w\in\mathbb{R}\times\mathbb{R}^m} \frac{1}{2} (v-x)^2 + \frac{1}{2} \sum_{j=1}^m (w_j - y_j)^2 + \lambda |v| \sum_{j=1}^m |w_j|.\tag{11}\] The signs of the elements of the decision variables in (11 ) are determined by the signs of \((x, y)\), and consequently, the problem in (11 ) is equivalent to problem 12 ; this is (partly) formulated by Lemma 12. \[\label{eq:8} \min_{v, w\in\mathbb{R}_+\times\mathbb{R}_+^m} h_\lambda(v,w;x,y) \equiv \frac{1}{2} (v-|x|)^2 + \frac{1}{2} \sum_{j=1}^m (w_j - |y_j|)^2 + \lambda v \sum_{j=1}^m w_j.\tag{12}\]

To summarize, we have established that, finding an optimal solution to 12 and then changing signs accordingly, yields an optimal solution to 11 . We will now focus on obtaining an optimal solution for (11 ).

11.2 Solving the Prox Problem↩︎

First we note that problem 12 is well-posed.

In light of Lemma 13, and due to the fact that in 12 we minimize a continuously differentiable function over a closed convex set, the set of optimal solutions of 12 is a nonempty subset of the set of stationary points of 12 . These satisfy the following conditions (cf. [24]).

The stationarity conditions given in Lemma 14 imply a solution form that we exploit in Algorithm 2; this is described by Corollary 3, where we use the convention that \(\sum_{j=1}^{0} a_j \equiv 0\) for any \(\{a_j\}\subseteq\mathbb{R}\).

In our analysis, we strictly distinguish between the trivial solution \((v^*,w^*) = (0,y)\), and the non-trivial solution in which \(v^*>0\). A practical point-of-view suggests that if \(v^* = 0\), then the corresponding succeeding weights should also be zero, even though the optimality conditions imply otherwise. However, to avoid hindering the training process, this observation is considered only in the end of the training.

Recall that our analysis so-far implies that the magnitude order of \(y\) determines the order magnitude of \(w\), effectively implying on set of non-zero entries in \(w\) (cf. Remark 5). For clarity of indices, and without loss of generality, we assume throughout this section that the vector \(y\) is already sorted in decreasing order of magnitude, that is \(y \equiv \bar{y}\). We will use, without confusion, both notation to describe the same entity in order to maintain coherence with our procedures and results.

Denote \[\label{eq:vsws} \begin{align} v^{(s)} &= \dfrac{1}{1-s \lambda^2}\left(|x| - \lambda\sum_{j=1}^s|y_j|\right) \\ w^{(s)}_j &= |y_j| - \lambda v^{(s)} \;\text{ for } j=1,2,\ldots,s, \;\text{ and } w^{(s)}_j = 0 \;\text{otherwise}. \end{align}\tag{13}\] Lemma 5 which states the monotonicity property \[h_\lambda(v^{(s)}, w^{(s)}; x, y) < h_\lambda(v^{(s-1)}, w^{(s-1)}; x, y).\] is proved next.

We can now prove our key result formulated in Corollary 2, that states that \((v^{(s^*)}, w^{(s^*)} )\) is an optimal solution of 6 for \[s^* = \max \left\lbrace s \in [\bar{s}] : v^{(s)}, w^{(s)} > 0\right\rbrace, \;\text{ where } \;\bar{s} = \min(\lfloor \lambda^{-2}\rfloor, m).\]

Finally, we will show that the problem of finding \(s^*\) can be easily solved using binary search. To this end, we show that the feasibility criterion (i.e., \(v^{(s)} > 0\) and \(w^{(s)} > 0\)) satisfies that \[(v^{(k)}, w^{(k)}) \text{ is feasible } \Rightarrow (v^{(i)}, w^{(i)}) \text{ is feasible } \forall i < k\]

12 Multi-output proximal map computation↩︎

12.1 Solving the prox problem↩︎

Returning to the multi-output setting, recall that \(h_{V,W}(x) = V^T \sigma(Wx)\) with \(V \in \mathbb{R}^{p\times n}, W\in \mathbb{R}^{n\times m}\) and \[g(V, W)= \sum_{i=1}^n\sum_{j=1}^m\sum_{k=1}^p W_{ij}V_{ki}.\] The proximal mapping can then be written as: \[\begin{align} \text{prox}_{\lambda g}(\bar{V}, \bar{W}) &= \displaystyle\mathop{\rm argmin}_{V, W} \frac{1}{2} \|V - \bar{V}\|_F + \frac{1}{2} \|W - \bar{W}\|_F + \lambda \sum_{i=1}^n\sum_{j=1}^m\sum_{k=1}^p W_{ij}V_{ki} \\ &= \displaystyle\mathop{\rm argmin}_{V, W} \sum_{i=1}^n \left( \frac{1}{2} \sum_{k=1}^p (V_{ki} - \bar{V}_{ki})^2 + \sum_{j=1}^p (W_{ij} - \bar{W}_{ij})^2 + \sum_{j=1}^m\sum_{k=1}^p W_{ij}V_{ki}\right). \end{align}\] Noting that the proximal mapping is separable with respect to the columns of \(W\) and the rows of \(V\), and using the same sign trick applied in the single-output case, it is enough to solve for any \(i=1,\ldots, n\), \[\min_{v, w\in \mathbb{R}_+^p \times \mathbb{R}_+^m} h_\lambda(v,w;x,y) \equiv \frac{1}{2} \sum_{k=1}^p (v_k - |x_k|)^2 + \frac{1}{2} \sum_{j=1}^m (w_j - |y_j|)^2 + \lambda \sum_{j=1}^m\sum_{k=1}^p v_k w_j, \label{eq:multi-ouptut-prox}\tag{15}\] where \(x\) denotes the \(i\)-th row of \(V\) and \(y\) the \(i\)-th column of \(W\), in order to compute the prox operator.

The stationarity conditions for 15 are stated next; the arguments are the same as in the single-output case.

The next lemma restates the result in Lemma 7 which expands on the monotonic relation in magnitude originally established for single-output networks in Corollary 1.

Finally, we show, as in the single-output case, that second order stationarity condition constraints the ranges of sparsities of \(v^*\) and \(w^*\); this relation is given by Lemma 8, and is proved next.

A possible way of solving this proximal problem is thus to exhaustively compute the value of \(h_\lambda\) at each stationary point associated with sparsities \(s_v = 1, \ldots, p\), \(s_w = 1, \ldots, m\) such that \(s_v s_w \leq \lambda^{-2}\). However, trying all possible pairs of sparsities \((s_v, s_w)\) is computationally costly. Similarly as is the single output case, we can exploit some structure of the objective function \(h_\lambda\) in order to reduce the possible candidate pairs of sparsities.

Without loss of generality, we assume hereafter that the vectors \(x, y\) are already sorted in decreasing order of magnitude.

Lemma 16 shows that for each pair \((s_v, s_w)\), \(s_v = 0,1,\ldots,p\), \(s_w = 0,1,\ldots,m\), there exists a stationary point \((v^{(s_v, s_w)}, w^{(s_v, s_w)})\) of \(h_\lambda(\cdot, \cdot; x,y)\) such that \(|\{k: v^{(s_v, s_w)}_k > 0\}| = s_v\), \(|\{j: w^{(s_v, s_w)}_j > 0\}| = s_w\), given by \[\label{eq:vsws-multi} \begin{align} v^{(s_v, s_w)}_k &= |x_k| + \dfrac{1}{1 - s_v s_w \lambda^2}\left(\lambda^2 s_w \sum_{l=1}^{s_v} |x_l| - \lambda\sum_{j=1}^{s_w} |y_j|\right) \;\text{ for } k=1,2,\ldots,s_v, \;\text{ and } v^{(s_v, s_w)}_k = 0 \;\text{otherwise} \\ w^{(s_v, s_w)}_j &= |y_j| + \dfrac{1}{1 - s_v s_w \lambda^2}\left(\lambda^2 s_v \sum_{l=1}^{s_w} |y_l| - \lambda\sum_{k=1}^{s_v} |x_k|\right) \;\text{ for } j=1,2,\ldots,s_w, \;\text{ and } w^{(s_v, s_w)}_j = 0 \;\text{otherwise}. \end{align}\tag{16}\]

We now move to prove the monotonicity property stated in Lemma 9.

In order to derive an efficient algorithm , we will again exploit the monotone property of the feasibility criterion \(v^{(s_v, s_w)} > 0\), \(w^{(s_v, s_w)} > 0\) restated from Lemma 10:

We now introduce the efficient procedure to compute the maximal feasibility boundary (MFB), and prove that it indeed delivers, as promised, all sparsity pairs in the MFB set.

12.1.0.3 Time complexity of Algorithm 5

At each iteration of the loop, either \(s_v\) is incremented by \(1\) or \(s_w\) is decremented by \(1\). Since \(s_v\) starts from \(0\) and \(s_w\) from \(m\), and that the stopping criterion is \(s_v > p\) or \(s_w < 0\), it follows that the maximal number of iterations inside the loop is \(m+p\). At each iteration, we must compute \(v_{s_v}^{(s_v, s_w)}\) and \(w_{s_w}^{(s_v, s_w)}\), which requires in particular to compute \(\sum_{k=1}^{s_v} |x_k|\) and \(\sum_{j=1}^{s_w} |y_j|\). However, these cumulative sums can be efficiently computed before the loop in time \(\mathcal{O}(m+p)\), so that computing \(v_{s_v}^{(s_v, s_w)}\) and \(w_{s_w}^{(s_v, s_w)}\) inside the loop can be done in constant time. The overall complexity of this algorithm is thus \(\mathcal{O}(m+p)\).

Moreover, we can see that each time we add a pair to \(S\), we must both decrement \(s_w\) by \(1\) (just after adding the element in the algorithm), and increment \(s_v\) by \(1\) (in order for the boolean \(maximal\) to become true again). Therefore, there can be at most \(\min(m,p)\) pairs in the final set \(s\) at the end of the algorithm.

Merging all previous results, we can finally prove Theorem 4.

13 Experimental details and other plots↩︎

We consider the following values for the parameters that determine the training loop:

• batch size: 100

• epochs: 20

• learning rate: 1e-1, 1e-2, 1e-3, 1e-4, 5e-1, 5e-2, 5e-3, 5e-4

• dataset: mnist, fmnist, kmnist

• hidden neurons: 200

• lambda (\(\lambda\)): 0., 1e-5, 1e-4, 1e-3, 1e-2, 1e-1, 1e0, 1e1, 1e2, 2e-5, 2e-4, 2e-3, 2e-2, 2e-1, 2e0, 2e1, 2e2, 3e-5, 3e-4, 3e-3, 3e-2, 3e-1, 3e0, 3e1, 3e2, 4e-5, 4e-4, 4e-3, 4e-2, 4e-1, 4e0, 4e1, 4e2, 5e-5, 5e-4, 5e-3, 5e-2, 5e-1, 5e0, 5e1, 5e2

The \(\ell_\infty\)-bounded adversarial examples used to evaluate the robustness of the networks were generated using the PGD method described in [3] and implemented in the advertorch toolbox (https://github.com/BorealisAI/advertorch) using the following parameters:

• epsilon: 0.05, 0.1, 0.15, 0.2, 0.25, 0.3

• iterations: 40

• step size: epsilon / 20

• random initialization: True

13.1 sparsity per iteration↩︎

One advantage of the proximal mapping of the 1-path-norm and the \(\ell_1\)-norm is that they can set many weights to exactly zero. This has the effect of providing sparse networks from early iterations. This is in contrast to SGD with a constant stepsize which does not generate sparse iterates. In Figures 4, 5, 6 and 7 we plot the percentage of nonzero weights as a function of the iteration count, for both plain SGD and proximal SGD. We observe that in fact this is the case, and that the sparsity of the \(\ell_1\) and 1-path-norm regularized network can be controlled with the regularization parameter \(\lambda\).

13.2 Robustness vs accuracy tradeoff↩︎

For all possible values of \(\lambda\), in Figure 8 we plot the data corresponding to the lerning rate with least error. We plot the value of the error on clean samples and the error on adversarial examples. This allows us to understand the tradeoff between accuracy and robustness that is controlled by the regularization paramter \(\lambda\).

References↩︎