# Structured Sparsity Inducing Adaptive Optimizers for Deep Learning

Tristan Deleu  Yoshua Bengio

###### Abstract

The parameters of a neural network are naturally organized in groups, some of
which might not contribute to its overall performance. To prune out
unimportant groups of parameters, we can include some non-differentiable
penalty to the objective function, and minimize it using proximal gradient
methods. In this paper, we derive the weighted proximal operator, which is a
necessary component of these proximal methods, of two structured sparsity
inducing penalties. Moreover, they can be approximated efficiently with a
numerical solver, and despite this approximation, we prove that existing
convergence guarantees are preserved when these operators are integrated as
part of a generic adaptive proximal method. Finally, we show that this
adaptive method, together with the weighted proximal operators derived here,
is indeed capable of finding solutions with structure in their sparsity
patterns, on representative examples from computer vision and natural language
processing.

[width=]figures/convnet-sparsity

Figure 1: Structured sparsity in neural networks. (Left) An example of the
effect of unstructured sparsity inducing penalties, e.g. the $\ell_{1}$
penalty, on a convolutional layer. The zeros of each filter are represented in
blue. (Middle) An example of the effect of structured sparsity inducing
penalties, e.g. the mixed $\ell_{1}/\ell_{2}$ norm, on a similar convolutional
layer. The filters are grouped together, depending on their input channel, and
whole groups are encouraged to be equal to zero: this is a _channel-wise
structure_ (Wen et al., 2016). Here, the sparsity pattern effectively removes
the effect of the second input channel from the network’s predictions. (Right)
Graphs of the two structured sparsity inducing penalties (in 1-D, for
visualization) studied in this paper. The $\ell_{1}$ norm applies a larger
penalty for large values of $x$, whereas the Minimax Concave Penalty (MCP)
saturates after a certain point.

##  1 Introduction

In machine learning, a large majority of problems involve the minimization of
a composite loss function of the form

$\min_{{\bm{x}}\in{\mathbb{R}}^{N}}f({\bm{x}})+h({\bm{x}}),$ (1)

where $f$ is a differentiable loss function, and $h$ is a penalty function (or
_regularizer_). Some standard penalties used in the context of training neural
networks include the squared $\ell_{2}$ norm (also called _ridge regression_
in statistics, or _weight decay_ in machine learning) and elastic net (Zou &
Hastie, 2005). In this work, we focus on non-differentiable, and possibly non-
convex, penalties leading to sparse solutions of this optimization problem. In
particular, we are interested in penalties which encourage structure in the
sparsity patterns of these solutions (Kyrillidis et al., 2015), such as the
mixed $\ell_{1}/\ell_{2}$ norm, also known as _group Lasso_ when $f$ is the
squared loss (Yuan & Lin, 2006). These penalties tend to produce solutions
where complete groups of variables are zeroed-out, as opposed to the more
frequently used $\ell_{1}$ norm which yields unstructured sparsity. This class
of penalties is called _structured sparsity inducing_ penalties (Bach et al.,
2012a, b).

This family of regularizers is particularly appropriate for neural networks,
where the parameters are naturally organized in layers and groups of filters
(Krizhevsky et al., 2012; Ioannou et al., 2017). If some of these groups of
parameters were exactly equal to zero, this would effectively be equivalent to
deactivating certain neurons in the network (Zhou et al., 2016). As a
consequence, it may have an impact on the size of the model, reducing its
memory footprint and computational cost at inference time. However in these
structured sparsity inducing penalties, the groups of variables can be
arbitrary, and they do not necessarily have to follow the layered nature of
the neural network if we have some prior knowledge about the structure of the
problem, such as a causal structure (Germain et al., 2015; Lachapelle et al.,
2020; Ke et al., 2019), or for feature selection.

To solve an optimization problem like the one in Eq. 1, it is standard
practice in Deep Learning to use first-order methods, and in particular
adaptive methods such as Adagrad (Duchi et al., 2011), RMSprop (Tieleman &
Hinton, 2012), or Adam (Kingma & Ba, 2015). Even when the penalty function $h$
is non-differentiable, as is the case for structured sparsity inducing
penalties, it is still possible to optimize this composite loss function using
stochastic subgradient methods (Wen et al., 2016). But even if these penalties
are meant to introduce sparsity, the solutions found by subgradient methods
are typically not sparse (Bach et al., 2012a), and they often require a post-
processing step where the parameters with a small magnitude are pruned out. An
alternative approach is to use _proximal gradient methods_ , which are
specifically designed for the optimization of composite loss functions. These
methods avoid choosing an arbitrary subgradient of $h$ at non-differentiable
points by applying a _proximal operator_ (Moreau, 1962), which can often be
efficiently computed, to a standard gradient descent update rule.

While proximal gradient methods are often presented as derived from vanilla
gradient descent, it has been shown that they can also be applied to adaptive
first-order methods (Yang et al., 2020; Yun et al., 2020), using a weighted
version of the proximal operator (Hiriart-Urruty & Lemaréchal, 1993). For many
penalties, including the $\ell_{1}$ penalty (Becker et al., 2019), the
weighted proximal operator is a straightforward extension of its unweighted
counterpart, and can often be computed in closed-form. However for structured
sparsity inducing penalties, computing their weighted proximal operator
becomes more challenging than their unweighted equivalent, due to their
grouping property.

In this paper, we derive the _weighted proximal operator_ of two different
structured sparsity inducing penalties: the mixed $\ell_{1}/\ell_{2}$ norm,
and the group Minimax Concave Penalty (group MCP; Breheny & Huang, 2009). We
show that unlike their (unweighted) proximal operator, no closed-form solution
exists for the weighted proximal operator of these penalties, although they
can be approximated efficiently using a root finding algorithm. When this
approximation is integrated as part of a stochastic proximal gradient method
(Yun et al., 2020), we prove that the resulting inexact method retains the
same convergence guarantees as if the exact weighted proximal operator was
applied. We further support empirically that this algorithm can find solutions
with structure in their sparsity patterns, with examples from computer vision
and natural language processing.

##  2 Related work

#### Pruning

In addition to having a smaller memory footprint, pruning groups of variables
that only have a small impact on the network’s performance, such as filters in
a convolutional neural network, can also significantly speedup inference on
standard architectures (Li et al., 2017). Closely related to our work, Wen et
al. (2016) used the mixed $\ell_{1}/\ell_{2}$ norm in order to identify these
groups, and minimized the objective with subgradient methods. However this
required a post-processing step at the end of training because subgradient
methods may fail at finding truly sparse solutions (Bach et al., 2012a). In
contrast, proximal gradient methods tend to not suffer from this problem, and
they can find solutions with exact zeros. Lin et al. (2019) also proposed to
use structured sparsity regularization for pruning, with an approach inspired
by ADMM.

#### Proximal gradient methods

Proximal gradient methods (also known as _forward-backward splitting_ methods;
Duchi & Singer, 2009) is a class of first-order methods that are particularly
designed for the optimization of composite objectives of the form of Eq. 1,
where the non-differentiable penalty $h$ admits a proximal operator that can
be efficiently computed. See Parikh & Boyd (2014) for a general introduction
to proximal algorithms, and Bach et al. (2012a) for an overview of proximal
gradient methods for sparsity-inducing penalties (including structured
sparsity).

Proximal gradient methods have also been embedded into adaptive first-order
methods, such as Adam or RMSProp (Lee & Lee, 2019). Yang et al. (2020)
introduced Prox-SGD, an algorithm based on a weighted proximal operator to
correct for the rescaling introduced by these adaptive methods. This algorithm
was further refined in ProxGen (Yun et al., 2020), to treat the update as an
exact proximal gradient step. Both of these methods are generic, and they can
be applied to any penalty, provided the weighted proximal operator can be
computed efficiently. But so far, their applications have been limited to
sparsity-inducing penalties such as $\ell_{q}^{q}$, with $q\in[0,1]$, whose
weighted proximal operators can be obtained in closed-form thanks to the
separability of these penalties (Becker et al., 2019). Our work is
complementary, since we give an efficient algorithm to compute the weighted
proximal operator of some structured sparsity inducing penalties, for which
there is no closed-form solution available. Melchior et al. (2019) also
propose an alternative algorithm where the weighted proximal operator is
approximated by iterating over the (unweighted) proximal operator.

Finally, although we are focusing our attention here on first-order methods,
which are more appropriate for Deep Learning applications, it is important to
note that there also exists proximal quasi-Newton methods (Becker et al.,
2019), also based on weighted proximal operators, which are adapted from
limited-memory quasi-Newton methods such as SR1 and L-BFGS.

##  3 Structured sparsity inducing penalties

Throughout this paper, we assume that we have a collection ${\mathcal{G}}$ of
(disjoint) groups of variables; for example, these groups might correspond to
individual filters in a convolutional neural network (Wen et al., 2016).
Structured sparsity inducing penalties are penalty functions that encourage
whole groups of variables as defined by ${\mathcal{G}}$ to be ignored (see
Fig. 1 for an illustration). In this section, we will first recall the
definition of two of these penalties, namely the mixed $\ell_{1}/\ell_{2}$
norm and group MCP, and how they have been used in the past in conjunction
with proximal gradient descent.

###  3.1 Mixed $\ell_{1}/\ell_{2}$ norm

Taking inspiration from the $\ell_{1}$ penalty, which is widely adopted for
unstructured sparsity, the mixed $\ell_{1}/\ell_{2}$ norm can be thought of as
the $\ell_{1}$ norm of a vector consisting of the $\ell_{2}$ norms over each
group in ${\mathcal{G}}$. More precisely, this penalty is defined as

$h({\bm{x}})\triangleq\sum_{g\in{\mathcal{G}}}\lambda_{g}\|{\bm{x}}_{g}\|_{2},$
(2)

where $\|{\bm{x}}_{g}\|_{2}=\sqrt{\sum_{j\in g}x_{j}^{2}}$ is the $\ell_{2}$
norm of the vector ${\bm{x}}$, restricted to group $g$. To account for groups
of varying sizes, the weighting $\lambda_{g}$ is typically set to
$\lambda_{g}=\lambda\sqrt{|g|}$, where $|g|$ is the number of elements in
group $g$ (Murphy, 2012).

###  3.2 Bias reduction with non-convex penalties

Although applying the $\ell_{1}$ penalty may result in sparse solutions, this
penalty is also known to suffer from a _shrinkage_ effect (Copas, 1983),
introducing bias in the model. The mixed $\ell_{1}/\ell_{2}$ norm inherits
from this bias, but this time at the group level: groups of variables get
pushed invariably towards $0$, regardless of their $\ell_{2}$ norm.

In order to reduce this bias, Zhang (2010) introduced a non-convex penalty
called the _Minimax Concave Penalty_ (MCP). MCP operates under two regimes:
similarly to the $\ell_{1}$ penalty, small values are shrunk towards $0$, but
contrary to $\ell_{1}$, sufficiently large values are not penalized as much
anymore. Formally, this penalty function is defined as

$\operatorname{MCP}(x;\beta,\lambda)\triangleq\left\\{\begin{array}[]{cl}\lambda|x|-\dfrac{x^{2}}{2\beta}&\textrm{if
$|x|\leq\beta\lambda$}\\\
\dfrac{\beta\lambda^{2}}{2}&\textrm{otherwise,}\end{array}\right.$ (3)

with $\beta>1$ controlling the magnitude beyond which values are no longer
penalized. Fig. 1 shows a comparison of $\operatorname{MCP}$ against the
$\ell_{1}$ penalty. Similar to the mixed $\ell_{1}/\ell_{2}$ norm, it is also
possible to apply MCP to the $\ell_{2}$ norm of groups of variables, in order
to encourage structured sparsity while retaining the benefits of MCP: this is
called group MCP (Breheny & Huang, 2009)

$h({\bm{x}})=\sum_{g\in{\mathcal{G}}}\operatorname{MCP}(\|{\bm{x}}_{g}\|_{2};\beta,\lambda_{g}).$
(4)

###  3.3 Proximal gradient descent

In order to minimize the composite objective in Eq. 1, the proximal gradient
descent algorithm updates an iterate ${\bm{x}}_{t}$ with an update similar to
standard gradient descent, based on the gradient of the differentiable part
$f$

${\bm{x}}_{t+1}\leftarrow\operatorname{prox}_{\alpha
h}\\!\big{(}{\bm{x}}_{t}-\alpha\nabla_{{\bm{x}}}f({\bm{x}}_{t})\big{)},$

where $\alpha$ is the learning rate, and $\operatorname{prox}_{\alpha h}$ is
the _proximal operator_ (Moreau, 1962) of the penalty $h$, defined as

$\operatorname{prox}_{\alpha
h}({\bm{x}})=\operatorname*{argmin}_{{\bm{z}}\in{\mathbb{R}}^{n}}\frac{1}{2}\|{\bm{z}}-{\bm{x}}\|_{2}^{2}+\alpha
h({\bm{z}}).$ (5)

Interestingly, proximal gradient descent generalizes both standard gradient
descent (when $h({\bm{x}})\equiv 0$), and projected gradient descent (when
$h({\bm{x}})={\bm{I}}_{C}(x)$ is the indicator function, equal to $0$ if
${\bm{x}}\in C$, and $+\infty$ otherwise). For many penalty functions of
interest, even non-differentiable and possibly non-convex ones, the proximal
operator can be computed efficiently, often in closed-form. This is true in
particular for the structured sparsity inducing penalties presented in this
paper.

Let’s take the example of the mixed $\ell_{1}/\ell_{2}$ norm. First note that
since the penalty in Eq. 2 decomposes along the (disjoint) groups
$g\in{\mathcal{G}}$, the proximal operator in Eq. 5 benefits from the same
decomposition, involving only the proximal operators of the $\ell_{2}$ norm
evaluated on each ${\bm{x}}_{g}$. Therefore, it is sufficient to know the
proximal operator of the $\ell_{2}$ norm $h({\bm{x}})=\lambda\|{\bm{x}}\|_{2}$
in order to construct the update for proximal gradient descent, which is given
in closed-form (Combettes & Wajs, 2005) by

$\operatorname{prox}_{\alpha
h}({\bm{x}})=\left[1-\frac{\alpha\lambda}{\|{\bm{x}}\|_{2}}\right]_{+}\\!{\bm{x}},$
(6)

where $[z]_{+}=\max(0,z)$. Similarly, the group MCP enjoys the same
decomposition property along groups, and there exists a closed-form expression
for the proximal operator of the MCP / $\ell_{2}$ penalty (i.e.
$h({\bm{x}})=\operatorname{MCP}(\|{\bm{x}}\|_{2})$; Breheny & Huang, 2009).

##  4 Adaptive proximal optimizers for structured sparsity inducing penalties

In practice in Deep Learning, the minimization of the objective is typically
performed using stochastic and adaptive first-order methods, such as Adagrad,
RMSprop, or Adam, where the update direction is rescaled with a diagonal
preconditioning matrix ${\bm{D}}_{t}$

${\bm{x}}_{t+1}\leftarrow{\bm{x}}_{t}-\alpha{\bm{D}}_{t}^{-1}{\bm{m}}_{t},$

where ${\bm{m}}_{t}$ is an estimate of the mean of the gradients of $f$; see
also Table 3 in Section C.1 for examples of matrices ${\bm{D}}$. Although it
appears to be similar to the standard gradient descent update, the
preconditioning ${\bm{D}}$ makes composing these methods with a proximal
operator more challenging. Yun et al. (2020) proposed a general proximal
method that is capable of dealing with this preconditioning, with an update of
the form

${\bm{x}}_{t+1}\leftarrow\operatorname{prox}_{\alpha
h}^{{\bm{D}}_{t}}\\!\big{(}{\bm{x}}_{t}-\alpha{\bm{D}}_{t}^{-1}{\bm{m}}_{t}\big{)},$
(7)

where $\operatorname{prox}_{\alpha h}^{{\bm{D}}}$ is the _weighted proximal
operator_ of $h$ (Hiriart-Urruty & Lemaréchal, 1993; Lee et al., 2014)

$\operatorname{prox}_{\alpha
h}^{{\bm{D}}}({\bm{x}})=\operatorname*{argmin}_{{\bm{z}}\in{\mathbb{R}}^{n}}\frac{1}{2}\|{\bm{z}}-{\bm{x}}\|_{{\bm{D}}}^{2}+\alpha
h({\bm{z}}),$ (8)

and $\|{\bm{z}}\|_{{\bm{D}}}^{2}=\langle{\bm{z}},{\bm{D}}{\bm{z}}\rangle$.
This naturally extends to the case where the update is stochastic, see
Appendix C, and Section 5 for details and additional theoretical guarantees.
Contrary to the proximal operator, for which proximity is computed using the
Euclidean distance, the weighted proximal operator captures the reweighting by
${\bm{D}}$ by changing the geometry of the space over which proximity is
measured.

However unlike their proximal operators (see Section 3.3), the weighted
proximal operators of the $\ell_{2}$ norm and MCP / $\ell_{2}$ cannot be
evaluated in closed-form in general. Becker et al. (2019); Yang et al. (2020)
show that we can get a closed-form expression for the mixed
$\ell_{1}/\ell_{2}$ norm when the preconditioning matrix ${\bm{D}}$ decomposes
along the groups into spherical matrices ${\bm{D}}_{g}=m{\bm{I}}_{g}$, which
is not practical for Deep Learning applications: while being diagonal, the
preconditioning matrix usually does not have the same value for all the
variables belonging to the same group (e.g. it might be based on aggregated
statistics from past gradients). In Theorem 1, we give an expression of the
weighted proximal operator for the $\ell_{2}$ norm, based on an implicit
expression. Recall from Section 3.3 that due to group decomposition, this is
sufficient to compute the weighted proximal operator of the mixed
$\ell_{1}/\ell_{2}$ norm.

######  Theorem 1 (Weighted proximal operator of $\ell_{2}$).

Let ${\bm{D}}=\operatorname{diag}(d_{1},\ldots,d_{n})$ be a positive definite
diagonal matrix (i.e. $d_{i}>0$ for all $i$), and
$h({\bm{x}})=\lambda\|{\bm{x}}\|_{2}$ the $\ell_{2}$ penalty. The weighted
proximal operator of $h$ is given by

$\big{[}\operatorname{prox}_{\alpha
h}^{{\bm{D}}}({\bm{x}})\big{]}_{i}=\left\\{\begin{array}[]{cl}\dfrac{d_{i}\theta^{\star}x_{i}}{d_{i}\theta^{\star}+\alpha\lambda}&\textrm{if
$\|{\bm{D}}{\bm{x}}\|_{2}>\alpha\lambda$}\\\\[10.00002pt]
0&\textrm{otherwise},\end{array}\right.$ (9)

where $\theta^{\star}>0$ is the unique positive solution of

$\sum_{i=1}^{n}\left[\frac{d_{i}x_{i}}{d_{i}\theta^{\star}+\alpha\lambda}\right]^{2}=1.$
(10)

The proof of Theorem 1 is given in Section A.1. Unlike the (unweighted)
proximal operators for many standard penalties, including for the $\ell_{2}$
norm itself, the weighted proximal operator of $\ell_{2}$ cannot be expressed
in closed-form. Instead in practice, we need to solve Eq. 10 numerically, for
example using the Newton-Raphson algorithm. We note that Duchi et al. (2011)
also provides a bisection-based procedure to solve a somewhat different
minimization problem, related to Eq. 8 up to a change of variable. Theorem 1,
on the other hand, gives a more explicit formulation of the weighted proximal
operator when the preconditioning matrix is diagonal, more appropriate for the
adaptive proximal gradient method presented above (Yun et al., 2020).

One consequence of this approximation of $\theta^{\star}$ is that we can only
obtain an approximation of the weighted proximal operator of $\ell_{2}$, as
opposed to an exact expression; the resulting proximal gradient method is then
called _inexact_ (Schmidt et al., 2011). Eq. 11 shows that this error on the
weighted proximal operator is controlled by the approximation error induced by
the root finding algorithm used to find $\theta^{\star}$.

######  Proposition 1.

Let $\varepsilon>0$ and $\tilde{\theta}>0$ such that
$|\tilde{\theta}-\theta^{\star}|\leq\varepsilon$. Let
$\widetilde{\operatorname{prox}}^{{\bm{D}}}_{\alpha h}({\bm{x}})$ be an
approximation of the weighted proximal operator of $\ell_{2}$ (i.e.
$h({\bm{x}})=\lambda\|{\bm{x}}\|_{2}$), replacing $\theta^{\star}$ by
$\tilde{\theta}$ in Eq. 9. Then we have for all ${\bm{x}}\in{\mathbb{R}}^{n}$:

$\big{\|}\widetilde{\operatorname{prox}}_{\alpha
h}^{{\bm{D}}}({\bm{x}})-\operatorname{prox}_{\alpha
h}^{{\bm{D}}}({\bm{x}})\big{\|}_{2}\leq\varepsilon.$ (11)

The proof of Eq. 11 is provided in Section A.2, along with bounds on
$\theta^{\star}$ to reduce the search space as much as possible, and to ensure
fast convergence of the numerical solver.

Likewise, we also give in Theorem 2 an (implicit) expression for the weighted
proximal operator of MCP / $\ell_{2}$. Again, this is sufficient to compute
the weighted proximal operator of group MCP.

######  Theorem 2 (Weighted proximal operator of MCP / $\ell_{2}$).

Let ${\bm{D}}=\operatorname{diag}(d_{1},\ldots,d_{n})$ be a positive definite
diagonal matrix (i.e. $d_{i}>0$ for all $i$), and
$h({\bm{x}})=\operatorname{MCP}(\|{\bm{x}}\|_{2})$ the MCP / $\ell_{2}$
penalty. Suppose that $\alpha$ and $\beta$ satisfy $\alpha<\beta d_{\min}$,
where $d_{\min}$ is the smallest value of the diagonal of ${\bm{D}}$. The
weighted proximal operator of $h$ is given by

$\big{[}\operatorname{prox}_{\alpha
h}^{{\bm{D}}}({\bm{x}})\big{]}_{i}=\left\\{\begin{array}[]{cl}x_{i}&\textrm{\scalebox{0.83}{if
$\|{\bm{x}}\|_{2}>\beta\lambda$}}\\\\[8.00003pt]
\dfrac{d_{i}\beta\theta^{\star}x_{i}}{(d_{i}\beta-\alpha)\theta^{\star}+\alpha\beta\lambda}&\textrm{\scalebox{0.62}{\begin{tabular}[]{@{}l}if
$\|{\bm{x}}\|_{2}\leq\beta\lambda$\\\\[3.99994pt] and
$\|{\bm{D}}{\bm{x}}\|_{2}>\alpha\lambda$\end{tabular}}}\\\\[15.00002pt]
0&\textrm{\scalebox{0.62}{\begin{tabular}[]{@{}l}if
$\|{\bm{x}}\|_{2}\leq\beta\lambda$\\\\[3.99994pt] and
$\|{\bm{D}}{\bm{x}}\|_{2}\leq\alpha\lambda$,\end{tabular}}}\end{array}\right.$

(12)

where $\theta^{\star}>0$ is the unique positive solution of

$\beta^{2}\sum_{i=1}^{n}\left[\frac{d_{i}x_{i}}{(d_{i}\beta-\alpha)\theta^{\star}+\alpha\beta\lambda}\right]^{2}=1.$
(13)

The proof is provided in Section A.3. Moreover, guarantees on the
approximation error induced by the numerical solver to determine
$\theta^{\star}$ for MCP / $\ell_{2}$, similar to Eq. 18, are given in Section
A.4.

##  5 Convergence analysis

In this section, we will make stochasticity more explicit, and we will
consider the following composite objective function

$\min_{{\bm{x}}\in{\mathbb{R}}^{N}}\big{(}F({\bm{x}})\triangleq\mathbb{E}_{\xi}[f({\bm{x}};\xi)]+h({\bm{x}})\big{)}.$
(14)

Prior work studied the non-asymptotic convergence properties of stochastic
proximal gradient descent on this kind of objective (Xu et al., 2019; Yun et
al., 2020). However, their results depend on the fact that the (weighted)
proximal operator of $h$ can be computed exactly, which is unfortunately
impossible in practice for structured sparsity inducing penalties (their
weighted version at least). We will see that these results remain valid,
despite the approximation of the weighted proximal operator. Note that related
to our work, convergence guarantees of (deterministic) inexact proximal
gradient methods also exist, although with stronger conditions on the
convexity of $f$ and $h$ (Schmidt et al., 2011).

Because the overall objective function $F$ can be non-differentiable and non-
convex, convergence is proven in terms of the expected distance of the
_Frechet subdifferential_ $\widehat{\partial}F({\bm{x}}_{t})$ to zero111This
reduces to $\|\nabla F({\bm{x}}_{t})\|_{2}$ when $F$ is differentiable; see
Eq. 31 for the definition of the Frechet subdifferential. In particular at a
stationary point ${\bm{x}}^{\star}$, we have
${\bm{0}}\in\widehat{\partial}F({\bm{x}}^{\star})$. (Rockafellar, 1976). To
derive the convergence bound, we borrow the assumptions made by Yun et al.
(2020), which we recall here:

######  Assumption 1 (Yun et al., 2020).

x

  1. 1.

L-smoothness The loss function $f$ is $L$-smooth and lower-bounded:
$f({\bm{x}}^{\star})>-\infty$ for the optimal solution ${\bm{x}}^{\star}$, and
$\forall{\bm{x}},{\bm{y}},\|\nabla f({\bm{x}})-\nabla f({\bm{y}})\|_{2}\leq
L\|{\bm{x}}-{\bm{y}}\|_{2}$.

  2. 2.

Bounded variance The stochastic gradient ${\bm{g}}_{t}=\nabla
f({\bm{x}}_{t};\xi_{t})$ is unbiased, and has bounded variance:
$\mathbb{E}_{\xi}[{\bm{g}}_{t}]=\nabla f({\bm{x}}_{t})$, and
$\mathbb{E}_{\xi}[\|{\bm{g}}_{t}-\nabla
f({\bm{x}}_{t})\|^{2}_{2}]\leq\sigma^{2}$.

  3. 3.

(i) The update and (ii) the stochastic gradient are bounded, and (iii) the
momentum parameter is exponentially decaying: (i)
$\|{\bm{x}}_{t+1}-{\bm{x}}_{t}\|_{2}\leq D$, (ii) $\|{\bm{g}}_{t}\|_{2}\leq
G$, (iii) $\rho_{t}=\rho_{0}\mu^{t}$, with $\mu\in[0,1)$.

  4. 4.

Sufficiently positive-definite For all $t$,
${\bm{D}}_{t}\succeq\delta{\bm{I}}$ (i.e. ${\bm{D}}_{t}-\delta{\bm{I}}$ is
positive semi-definite), and
$\alpha_{t}{\bm{D}}_{t}^{-1}\succeq\gamma{\bm{I}}$, with $\delta,\gamma>0$.

In addition to these assumptions, we also need to account for the possible
approximation of the weighted proximal operator of the structured sparsity
inducing penalties. Here, we add the following assumption

######  Assumption 2.

Let ${\bm{x}}_{t+1}^{\star}=\operatorname{prox}_{\alpha
h}^{{\bm{D}}_{t}}\\!\big{(}{\bm{x}}_{t}-\alpha{\bm{D}}_{t}^{-1}{\bm{m}}_{t}\big{)}$
be the exact proximal update. For all $t\geq 0$, we have either

  1. 1.

${\bm{x}}_{t+1}$ is the exact proximal update:
${\bm{x}}_{t+1}\leftarrow{\bm{x}}_{t+1}^{\star}$;

  2. 2.

or ${\bm{x}}_{t+1}$ is an $\varepsilon_{t+1}$-approximation of
${\bm{x}}_{t+1}^{\star}$,
$\|{\bm{x}}_{t+1}-{\bm{x}}_{t+1}^{\star}\|_{2}\leq\varepsilon_{t+1}$, and $h$
is $L^{\prime}$-smooth in the $\varepsilon_{t+1}$-ball ${\mathcal{B}}$ around
${\bm{x}}_{t+1}^{\star}$ (i.e. ${\bm{x}}_{t+1}\in{\mathcal{B}}$).

Moreover, we assume that the approximation of the weighted proximal operator
becomes more accurate as $t$ grows. Specifically, with the convention
$\varepsilon_{t+1}=0$ if ${\bm{x}}_{t+1}={\bm{x}}_{t+1}^{\star}$,

$\sum_{t=0}^{+\infty}\varepsilon_{t+1}^{2}=K<+\infty.$

2 is a slightly weaker assumption than that of Yun et al. (2020): instead of
always relying on an exact proximal update
(${\bm{x}}_{t+1}={\bm{x}}_{t+1}^{\star}$), we allow the update to be
approximated at points where the penalty is smooth. This condition is in
particular satisfied for both structured sparsity inducing penalties studied
in this paper, since the value of $\theta^{\star}$, which is the only
approximated quantity, is only used to define the weighted proximal operator
in a regime where the corresponding penalty is smooth (at least if
$\|{\bm{D}}{\bm{x}}\|_{2}>\alpha\lambda$). The decreasing sequence of
approximation errors is also a well accepted assumption in the context of
inexact proximal methods (Rockafellar, 1976). Under these conditions, we can
derive a bound on the convergence of stochastic proximal methods, similar to
the one given in Xu et al. (2019). Note that Algorithm 2 is given in Appendix
C.

######  Theorem 3.

Suppose that 1 and 2 are satisfied. If we run Algorithm 2 with a non-
increasing step-size $\alpha_{t}$, such that $\alpha_{0}<\delta/2L$, then the
output ${\bm{x}}_{R}$ of Algorithm 2, where $R$ is sampled uniformly in
$\\{1,\ldots,T\\}$, satisfies

$\displaystyle\mathbb{E}_{R}\big{[}\operatorname{dist}\\!\big{(}{\bm{0}},\widehat{\partial}F({\bm{x}}_{R})\big{)}^{2}\big{]}\leq\frac{C_{1}}{T}\sum_{t=0}^{T-1}\|{\bm{g}}_{t}-\nabla
f({\bm{x}}_{t})\|_{2}^{2}+\frac{C_{2}\Delta}{T}+\frac{C_{3}}{T},$

with $C_{1}$, $C_{2}$, and $C_{3}$ positive constants independent of $T$. Here
$\operatorname{dist}({\bm{z}},S)$ is the distance of a set $S$ to a point
${\bm{z}}$, defined as the minimal distance of any point in $S$ to ${\bm{z}}$.

The proof of Theorem 3 is given in Appendix B. We note that this matches
previous (non-asymptotic) convergence guarantees of adaptive proximal gradient
methods from Yun et al. (2020), albeit with the slightly weaker 2. Moreover,
the corollary results from Xu et al. (2019) of the convergence of mini-batch
stochastic proximal methods can be directly transposed here; in particular, in
the case where the mini-batch size is fixed:

######  Corollary 1 (Fixed mini-batch size).

If the assumptions of Theorem 3 are satisfied, with
$T=2(C_{2}\Delta+C_{3})/\varepsilon^{2}$ and with a fixed mini-batch size
$m_{t}$ with $m_{t}=2C_{1}\sigma^{2}/\varepsilon^{2}$, then the output
${\bm{x}}_{R}$ of Algorithm 2 satisfies

$\mathbb{E}_{R}\big{[}\operatorname{dist}\\!\big{(}{\bm{0}},\widehat{\partial}F({\bm{x}}_{R})\big{)}^{2}\big{]}\leq\varepsilon^{2},$

where $C_{1}$, $C_{2}$, and $C_{3}$ are the constants from Theorem 3. To have
$\mathbb{E}[\operatorname{dist}({\bm{0}},\widehat{\partial}F({\bm{x}}_{R}))]\leq\varepsilon$,
it is then sufficient to have $T=O(1/\varepsilon^{2})$, making the total
complexity $O(1/\varepsilon^{4})$.

##  6 Experimental results

To validate that proximal gradient methods with structured sparsity inducing
penalties are indeed capable of finding solutions with structure in their
sparsity patterns, we experimented with two families of representative
architectures: convolutional neural networks from computer vision, and
transformers from natural language processing. In all experiments with
proximal gradient methods, we used the ProxGen algorithm (Yun et al., 2020)
presented in Section 4, and the weighted proximal operators given in Theorems
1 and 2; the adaptive optimizer (i.e. the form of the preconditioning matrix)
is Adam (Kingma & Ba, 2015). In order to approximate the weighted proximal
operators, we used the Newton-Raphson algorithm for its fast convergence
properties; see Algorithm 3 for details about this procedure. We will get back
to the choice of this algorithm in Section 6.4. The code is available here:
https://github.com/tristandeleu/pytorch-structured-sparsity

###  6.1 Convolutional Neural Networks

To show the advantage of structured sparsity in convolutional architectures,
we trained a VGG-16 (Simonyan & Zisserman, 2015) on CIFAR-10 using both the
mixed $\ell_{1}/\ell_{2}$ norm and group MCP. The penalty is only applied to
the weights of the neural network, leaving the biases unpenalized. Following
Wen et al. (2016), we study _channel-wise_ structured sparsity for the
convolutional layers, where the groups correspond to all the outgoing weights
of a single channel; see Fig. 1 for an illustration. Likewise, we use a _row-
wise_ structure for the weights of the unique fully-connected layer. This
structure is motivated by finding a network where intermediate representations
are only influenced by a subset of channels from the previous layer; we will
return to the advantages of this choice in Section 6.3. Overall, there are 4k
groups of variables, in a model containing about 15M parameters.

Table 1: Performance of VGG-16 trained on CIFAR-10, with different structured sparsity inducing penalties, reported as the mean and standard deviation over 3 runs. For each penalty, the first line corresponds to training with subgradient methods (Wen et al., 2016), and “+ prox.” with proximal gradient methods. Here, group sparsity is the proportion of groups (out of 4k) with non-zero norm. ⋆Required a thresholding step, see the text for details. | Group Sparsity | Test accuracy
---|---|---
Baseline | – | $90.76\pm 0.29\%\phantom{{}^{\star}}$
$\ell_{1}/\ell_{2}$ | $88.81\pm 0.09\%^{\star}$ | $88.56\pm 0.39\%^{\star}$
\+ prox. | $22.03\pm 0.33\%\phantom{{}^{\star}}$ | $89.55\pm 0.14\%\phantom{{}^{\star}}$
Group MCP | $88.98\pm 0.14\%^{\star}$ | $89.38\pm 0.67\%^{\star}$
\+ prox. | $22.63\pm 0.15\%\phantom{{}^{\star}}$ | $89.80\pm 0.06\%\phantom{{}^{\star}}$

With this experiment, we also want to show the effectiveness of proximal
gradient methods with the (non-differentiable) penalties studied in this
paper, as opposed to subgradient methods (Wen et al., 2016), where the
composite objective was simply trained with Adam. Table 1 compares the
performance for both training methods, with both structured sparsity inducing
penalties. In the case of subgradient methods, the solution found for both
penalties was not sparse, i.e. no group had zero norm, meaning a group
sparsity of $100\%$. Following Wen et al. (2016), we applied a post-processing
thresholding step to these networks, where groups with small enough $\ell_{2}$
norm were zeroed-out. We found that this step was highly sensitive to the
choice of the threshold, jumping from $90\%$ accuracy down to $10\%$ (i.e. a
random predictor) over a small range of values; the values reported in Table 1
are trading-off accuracy for group sparsity.

Despite tuning the value of the threshold, the amount of group sparsity for
subgradient methods remains high, showing little advantage in terms of the
number of non-zero parameters left. On the other hand, the solutions found
using proximal gradient methods are capable of reaching a lower group
sparsity, without any post-processing necessary (because the proximal operator
can set groups of parameters to $0$), while maintaining a good accuracy
compared to our baseline model.

Figure 2: Comparison of the size of the network before (light) and after
(dark) pruning, for each layer of VGG-16 (left) and BERT (right). For VGG-16,
convolutional layers (with their corresponding batch-normalization) are
represented in blue, and the linear layer in green; note that since there is
only a single linear layer at the end, this network has only 14 layers. The
label above each bar represents the size of the layer after pruning (all
convolutions have a kernel of size $3\times 3$).

###  6.2 Large-scale Transformers

Unlike in the large-scale computer vision community, which has mostly moved
away from adaptive optimizers (Wilson et al., 2017), Adam is still a popular
choice for training and fine-tuning large-scale transformers for language
modeling. To test our inexact proximal gradient method on a large-scale
problem, we fine-tuned a BERT model (Devlin et al., 2018) on SQuAD 1.1, a
question-answering benchmark (Rajpurkar et al., 2016); the language model was
initially pre-trained on a large corpus of articles from Wikipedia and the
BookCorpus dataset. Taking inspiration from our experiment with VGG-16, we use
a row-wise structure for the weight matrices of all the attention layers in
the network. In the context of BERT, Guo et al. (2019) found that this kind of
structure was already emerging in some weight matrices with unstructured
sparsity inducing penalties. Overall, there are 86k groups of variables, in a
model containing 108M parameters. To be consistent with prior work on pruning
language models (Sanh et al., 2020), the embeddings are kept fixed during
fine-tuning.

Table 2: Performance of BERT fine-tuned on SQuAD, with different structured sparsity inducing penalties, and different levels of group sparsity. Here, group sparsity is the proportion of groups (out of 86k) with non-zero norm. | Group Sparsity | EM | F1
---|---|---|---
Baseline | – | $81.01$ | $88.27$
$\ell_{1}/\ell_{2}$ | $84.34\%$ | $72.29$ | $81.96$
$76.14\%$ | $66.16$ | $77.50$
Group MCP | $81.49\%$ | $75.40$ | $84.50$
$61.33\%$ | $69.47$ | $79.95$

Table 2 shows the performance of BERT, both in terms of F1 score and Exact
Match (EM), for two choices of hyperparameters ($\lambda$ and $\beta$),
leading to two values of group sparsity. In contrast to our experiments in
Section 6.1, the levels of group sparsity reached using the structured
sparsity inducing penalties are higher, meaning that fewer parameters are
zeroed-out. The gap in performance between the sparse models and the baseline
matches the gap found in prior work using unstructured sparsity (Sanh et al.,
2020), although with more limited levels of sparsity here. We can also observe
that at comparable levels of group sparsity, the solutions found with group
MCP tend to perform better than the ones found with the mixed
$\ell_{1}/\ell_{2}$ norm. These results are encouraging, and show that
structured sparsity can also be effective on large-scale models, and remains
practical even with a very large number of groups.

###  6.3 Pruning with indirect sparsity

Interestingly, as observed by Li et al. (2017), enforcing channel-wise and
row-wise structured sparsity at the level of one layer has an effect on the
computational efficiency of that one layer, but it has consequences on the
neighboring layers as well: if one channel sees its outgoing weights being
zeroed-out, then it does not contribute to the output of the network anymore,
and can therefore be ignored in upstream computations as well. This has a
significant impact on the effective size of the network, beyond group
sparsity, where connections can be further pruned _indirectly_ thanks to these
structured sparsity inducing penalties without affecting predictions. The
details of this pruning procedure are given in Section C.3.

Fig. 2 shows the proportion of parameters preserved after applying this
procedure on two networks: VGG-19 trained with $\ell_{1}/\ell_{2}$ (group
sparsity $22.03\%)$, and BERT fine-tuned with group MCP (group sparsity
$61.33\%$). The effective sparsity (i.e. the size of the pruned network,
divided by the size of the original network) is about $3\%$ for VGG-16, and
$48\%$ for BERT, with significant gains on some layers (e.g. 3 orders of
magnitude fewer parameters for Conv9-12 in VGG-16). For BERT, we are still
below the high (unstructured) sparsity of Sanh et al. (2020), where the same
level of performance is reached with about $3\%$ of the size of the original
network. To achieve better sparsity with structured sparsity inducing
penalties, Transformers in general could benefit from better choices of
groups, e.g. spanning over multiple weight matrices.

Figure 3: Number of iterations, normalized by the total number of groups (here
4k), to approximate the weighted proximal operator with 3 different
algorithms, while training a VGG-16 network on CIFAR-10. Note that Newton’s
algorithm reaches a number of iterations per group close to 1 at the end of
training (about 1.65 iterations / group), in part thanks to the amortized cost
of running the algorithm in a vectorized way.

However, a big advantage of structured sparsity, with our choice of groups, is
that the pruned models are smaller _dense_ networks: in addition to having a
smaller memory footprint, this can accelerate inference, without any special
requirement in terms of hardware (see Sanh et al. (2020) for discussions).

###  6.4 Approximation of the weighted proximal operator

Although our formulation of the weighted proximal operators for the mixed
$\ell_{1}/\ell_{2}$ norm and group MCP in Theorems 1 and 2 is independent of
the choice of the numerical solver to find $\theta^{\star}$, in practice we
chose the Newton-Raphson algorithm. In this section, we evaluate empirically
how efficient this algorithm is over the course of training, compared to other
approximations of the weighted proximal operator. This also helps us to
estimate the overhead induced by this routine over standard first-order
methods.

As a direct replacement of Newton’s algorithm, we consider a bisection-based
method, similar to the one described in Duchi et al. (2011). We also compare
it to AdaProx (Melchior et al., 2019), a general method that approximates the
weighted proximal operator by iterating over the (unweighted) proximal
operator. Since all three methods are iterative methods, we evaluate their
efficiency as the number of iterations per gradient update. Similar to Section
6.1, we trained a VGG-16 network on CIFAR-10 with the mixed
$\ell_{1}/\ell_{2}$ norm as our penalty. In order to control for variations
during the learning process, at each gradient update we approximate the
weighted proximal operator with all three methods, but we only apply the one
found by Newton’s algorithm. We also use the same tolerance
$\varepsilon=10^{-6}$ for all algorithms.

Fig. 3 shows the evolution of the number of iterations per group, and per
gradient update, with the number of training steps. Thanks to the increasing
group sparsity, the cost induced by both Newton’s algorithm and bisection
decreases over time: in the case of $\ell_{1}/\ell_{2}$, only non-zero
$\operatorname{prox}({\bm{x}})$ require an approximation of $\theta^{\star}$
(i.e. if $\|{\bm{D}}{\bm{x}}\|_{2}>\alpha\lambda$; see Theorem 1). AdaProx, on
the other hand, iterates over the proximal operator of $\ell_{2}$ regardless
of our knowledge of $\|{\bm{D}}{\bm{x}}\|_{2}$, and therefore has a non-
decreasing cost over the course of training. This cost is also significantly
higher than the other two methods, because AdaProx was often reaching the
maximum number of iterations, set at $100$ per gradient update. Finally,
Newton’s algorithm is about 3 times as fast as the bisection-based method,
showing a clear advantage in terms of efficiency, and supporting our choice of
this algorithm as part of our inexact proximal gradient method.

##  7 Conclusion

The weighted proximal operator is a quantity that plays a fundamental role in
adaptive proximal gradient methods, a class of methods derived from adaptive
optimizers in Deep Learning, specifically adapted for optimizing composite
objectives with non-differentiable penalties. In this work, we derived the
weighted proximal operators of two structured sparsity inducing penalties: the
mixed $\ell_{1}/\ell_{2}$ norm, and group MCP. We saw that unlike their
unweighted counterparts, there exists no closed-form expression for the
weighted proximal operators of these penalties. Nevertheless, we found that
they can be efficiently approximated using the Newton-Raphson algorithm, and
they can be incorporated as part of a general proximal gradient method called
ProxGen (Yun et al., 2020). We also proved that the convergence guarantees of
this algorithm were maintained, despite the resulting method being inexact. We
showed experimentally that this proximal gradient method applied to structured
sparsity inducing penalties was capable of finding sparse solutions, with
structure in their sparsity patterns, on examples from computer vision with
VGG-16, and natural language processing with BERT, with a moderate loss of
performance. Finally, we leveraged this structure to prune these models even
further, and obtain much smaller models while keeping them functionally
identical.

Going beyond pruning neural networks, proximal gradient methods, together with
the weighted proximal operators derived here, can be applied to any problem
where structure is known or expected a priori. In future work, we would like
to evaluate this method on other choices of groups, for example motivated by
the causal structure of the problem.

## Acknowledgement

We would like to thank Emmanuel Bengio, Marwa El Halabi, and Salem Lahlou for
discussions and feedback on the paper. We would like to also thank Olexa
Bilaniuk for his help with Torchscript to accelerate the approximation of the
weighted proximal operators.

## References

  * Bach et al. (2012a) Bach, F., Jenatton, R., Mairal, J., and Obozinski, G.  Optimization with Sparsity-Inducing Penalties.  _Foundations and Trends in Machine Learning_ , 2012a.
  * Bach et al. (2012b) Bach, F., Jenatton, R., Mairal, J., Obozinski, G., et al.  Structured Sparsity through Convex Optimization.  _Statistical Science_ , 2012b.
  * Becker et al. (2019) Becker, S., Fadili, J., and Ochs, P.  On Quasi-Newton Forward-Backward Splitting: Proximal Calculus and Convergence.  _SIAM Journal on Optimization_ , 2019.
  * Breheny & Huang (2009) Breheny, P. and Huang, J.  Penalized methods for bi-level variable selection.  _Statistics and its Interface_ , 2009.
  * Combettes & Wajs (2005) Combettes, P. L. and Wajs, V. R.  Signal recovery by proximal forward-backward splitting.  _Multiscale Modeling & Simulation_, 2005.
  * Copas (1983) Copas, J. B.  Regression, Prediction and Shrinkage.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 1983.
  * Devlin et al. (2018) Devlin, J., Chang, M.-W., Lee, K., and Toutanova, K.  BERT: Pre-training of Deep Bidirectional Transformers for Language Understanding.  _arXiv preprint_ , 2018.
  * Duchi et al. (2011) Duchi, J., Hazan, E., and Singer, Y.  Adaptive Subgradient Methods for Online Learning and Stochastic Optimization.  _Journal of Machine Learning Research_ , 2011.
  * Duchi & Singer (2009) Duchi, J. C. and Singer, Y.  Efficient Learning using Forward-Backward Splitting.  In _Advances in Neural Information Processing Systems_ , 2009.
  * Germain et al. (2015) Germain, M., Gregor, K., Murray, I., and Larochelle, H.  MADE: Masked Autoencoder for Distribution Estimation.  In _International Conference on Machine Learning_ , 2015.
  * Guo et al. (2019) Guo, F.-M., Liu, S., Mungall, F. S., Lin, X., and Wang, Y.  Reweighted Proximal Pruning for Large-Scale Language Representation.  _arXiv preprint_ , 2019.
  * Hiriart-Urruty & Lemaréchal (1993) Hiriart-Urruty, J.-B. and Lemaréchal, C.  _Convex Analysis and Minimization Algorithms II: Advanced Theory and Bundle Methods_.  Springer, 1993.
  * Ioannou et al. (2017) Ioannou, Y., Robertson, D., Cipolla, R., and Criminisi, A.  Deep Roots: Improving CNN Efficiency with Hierarchical Filter Groups.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2017.
  * Ke et al. (2019) Ke, N. R., Bilaniuk, O., Goyal, A., Bauer, S., Larochelle, H., Schölkopf, B., Mozer, M. C., Pal, C., and Bengio, Y.  Learning Neural Causal Models from Unknown Interventions.  _arXiv preprint_ , 2019.
  * Kingma & Ba (2015) Kingma, D. P. and Ba, J.  Adam: A Method for Stochastic Optimization.  _International Conference on Learning Representations_ , 2015.
  * Krizhevsky et al. (2012) Krizhevsky, A., Sutskever, I., and Hinton, G. E.  ImageNet Classification with Deep Convolutional Neural Networks.  In _Advances in Neural Information Processing Systems_ , 2012.
  * Kyrillidis et al. (2015) Kyrillidis, A., Baldassarre, L., El Halabi, M., Tran-Dinh, Q., and Cevher, V.  Structured Sparsity: Discrete and Convex approaches.  In _Compressed Sensing and its Applications_. Springer, 2015.
  * Lachapelle et al. (2020) Lachapelle, S., Brouillard, P., Deleu, T., and Lacoste-Julien, S.  Gradient-based Neural DAG Learning.  _International Conference on Learning Representations_ , 2020.
  * Lee et al. (2014) Lee, J. D., Sun, Y., and Saunders, M. A.  Proximal Newton-type methods for minimizing composite functions.  _SIAM Journal on Optimization_ , 2014.
  * Lee & Lee (2019) Lee, S. and Lee, J.  Compressed Learning of Deep Neural Networks for OpenCL-Capable Embedded Systems.  _Applied Sciences_ , 2019.
  * Li et al. (2017) Li, H., Kadav, A., Durdanovic, I., Samet, H., and Graf, H. P.  Pruning Filters for Efficient Convnets.  _International Conference on Learning Representations_ , 2017.
  * Lin et al. (2019) Lin, S., Ji, R., Li, Y., Deng, C., and Li, X.  Toward Compact ConvNets via Structure-Sparsity Regularized Filter Pruning.  _IEEE Transactions on Neural Networks and Learning Systems_ , 2019\.
  * Melchior et al. (2019) Melchior, P., Joseph, R., and Moolekamp, F.  Proximal Adam: Robust Adaptive Update Scheme for Constrained Optimization.  _arXiv preprint_ , 2019.
  * Moreau (1962) Moreau, J. J.  Fonctions convexes duales et points proximaux dans un espace hilbertien.  _Elsevier_ , 1962.
  * Murphy (2012) Murphy, K. P.  _Machine Learning: A Probabilistic Perspective_.  MIT press, 2012.
  * Parikh & Boyd (2014) Parikh, N. and Boyd, S.  Proximal Algorithms.  _Foundations and Trends in Optimization_ , 2014.
  * Rajpurkar et al. (2016) Rajpurkar, P., Zhang, J., Lopyrev, K., and Liang, P.  SQuAD: 100,000+ Questions for Machine Comprehension of Text.  _arXiv preprint_ , 2016.
  * Rockafellar (1976) Rockafellar, R. T.  Monotone operators and the proximal point algorithm.  _SIAM journal on control and optimization_ , 1976.
  * Rockafellar & Wets (2009) Rockafellar, R. T. and Wets, R. J.-B.  _Variational Analysis_.  Springer Science & Business Media, 2009.
  * Rumelhart et al. (1986) Rumelhart, D. E., Hinton, G. E., and Williams, R. J.  Learning representations by back-propagating errors.  _Nature_ , 1986.
  * Sanh et al. (2020) Sanh, V., Wolf, T., and Rush, A. M.  Movement Pruning: Adaptive Sparsity by Fine-Tuning.  _Advances in Neural Information Processing Systems_ , 2020.
  * Schmidt et al. (2011) Schmidt, M., Roux, N. L., and Bach, F.  Convergence rates of inexact proximal-gradient methods for convex optimization.  _arXiv preprint_ , 2011.
  * Simonyan & Zisserman (2015) Simonyan, K. and Zisserman, A.  Very Deep Convolutional Networks for Large-Scale Image Recognition.  _International Conference on Learning Representations_ , 2015.
  * Tieleman & Hinton (2012) Tieleman, T. and Hinton, G.  RMSprop: Divide the gradient by a running average of its recent magnitude.  _Neural networks for machine learning_ , 2012.
  * Wen et al. (2016) Wen, W., Wu, C., Wang, Y., Chen, Y., and Li, H.  Learning Structured Sparsity in Deep Neural Networks.  In _Advances in Neural Information Processing Systems_ , 2016.
  * Wilson et al. (2017) Wilson, A. C., Roelofs, R., Stern, M., Srebro, N., and Recht, B.  The Marginal Value of Adaptive Gradient Methods in Machine Learning.  In _Advances in Neural Information Processing Systems_ , 2017.
  * Xu et al. (2019) Xu, Y., Jin, R., and Yang, T.  Non-asymptotic Analysis of Stochastic Methods for Non-Smooth Non-Convex Regularized Problems.  In _Advances in Neural Information Processing Systems_ , 2019.
  * Yang et al. (2020) Yang, Y., Yuan, Y., Chatzimichailidis, A., van Sloun, R. J., Lei, L., and Chatzinotas, S.  ProxSGD: Training Structured Neural Networks under Regularization and Constraints.  In _International Conference on Learning Representations_ , 2020.
  * Yuan & Lin (2006) Yuan, M. and Lin, Y.  Model selection and estimation in regression with grouped variables.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 2006.
  * Yun et al. (2020) Yun, J., Lozano, A. C., and Yang, E.  A General Family of Stochastic Proximal Gradient Methods for Deep Learning.  _arXiv preprint_ , 2020.
  * Zhang (2010) Zhang, C.-H.  Nearly Unbiased Variable Selection under Minimax Concave Penalty.  _The Annals of Statistics_ , 2010.
  * Zhou et al. (2016) Zhou, H., Alvarez, J. M., and Porikli, F.  Less is More: Towards Compact CNNs.  _European Conference on Computer Vision_ , 2016.
  * Zou & Hastie (2005) Zou, H. and Hastie, T.  Regularization and Variable Selection via the Elastic Net.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 2005.

##  Appendix A Proofs

In this section, we prove the main results from the paper, namely Theorems 1
and 2, and Eq. 11. In addition to these results, we also prove additional
results on the approximation error for MCP / $\ell_{2}$ (i.e. the equivalent
of Eq. 11 for MCP / $\ell_{2}$), as well as bounds on the different
$\theta^{\star}$ used in the main theorems, which can be used to narrow down
the search space for the root-finding algorithm. The proofs of these two
theorems are based on the following lemma, which gives an expression of the
weighted proximal operator of $h$ in terms of the (unweighted) proximal
operator of another function:

######  Lemma 1 (Becker et al. 2019).

Let ${\bm{A}}$ be a positive definite matrix. Then we have

$\operatorname{prox}_{\alpha
h}^{{\bm{A}}}({\bm{x}})={\bm{A}}^{-1/2}\operatorname{prox}_{\alpha
h\circ{\bm{A}}^{-1/2}}({\bm{A}}^{1/2}{\bm{x}}).$ (15)

This lemma is valid for any positive definite matrix ${\bm{A}}$; in
particular, Becker et al. (2019) makes use of it to derive proximal quasi-
Newton algorithms based on low-rank approximations of the Hessian. Here, we
will only use it in the limited case where the matrix ${\bm{D}}$ is diagonal
(and positive definite).

###  A.1 Weighted proximal operator of $\ell_{2}$

Based on Lemma 1, in order to compute the weighted proximal operator of
$h({\bm{x}})=\lambda\|{\bm{x}}\|_{2}$, it is sufficient to compute the
(unweighted) proximal operator $h\circ{\bm{A}}^{-1/2}$. The following
proposition shows how to compute this operator, up to a change of variable (to
simplify the notation), for a diagonal positive definite matrix ${\bm{D}}$.

######  Proposition 2.

Let ${\bm{D}}=\operatorname{diag}(d_{1},\ldots,d_{n})$ be a positive definite
diagonal matrix (i.e. $d_{i}>0$ for all $i$), and
$h({\bm{x}})=\lambda\|{\bm{x}}\|_{2}$ the $\ell_{2}$ penalty. The proximal
operator of $h\circ{\bm{D}}$ is given by

$\big{[}\operatorname{prox}_{\alpha
h\circ{\bm{D}}}({\bm{x}})\big{]}_{i}=\left\\{\begin{array}[]{cl}\dfrac{\theta^{\star}x_{i}}{\theta^{\star}+\alpha\lambda
d_{i}^{2}}&\textrm{if
$\|{\bm{D}}^{-1}{\bm{x}}\|_{2}>\alpha\lambda$}\\\\[10.00002pt]
0&\textrm{otherwise,}\end{array}\right.$ (16)

where $\theta^{\star}>0$ is the unique positive solution of

$G(\theta^{\star})\triangleq\sum_{i=1}^{n}\left[\frac{d_{i}x_{i}}{\theta^{\star}+\alpha\lambda
d_{i}^{2}}\right]^{2}=1.$ (17)

###### Proof.

Let us first recall the definition of the proximal operator of
$h\circ{\bm{D}}$

$\operatorname{prox}_{\alpha
h\circ{\bm{D}}}({\bm{x}})\triangleq\operatorname*{argmin}_{{\bm{z}}\in{\mathbb{R}}^{n}}\underbrace{\frac{1}{2}\|{\bm{z}}-{\bm{x}}\|_{2}^{2}+\alpha\lambda\|{\bm{D}}{\bm{z}}\|_{2}}_{\triangleq\,F({\bm{z}})}.$

The function $F$ is a convex function, whose minimizers ${\bm{z}}^{\star}$
characterize the proximal operator of interest. We start with the special case
where ${\bm{z}}^{\star}={\bm{0}}$; by definition of a global minimizer, we
have

If we introduce a new function $H({\bm{y}})\triangleq
1/2\|{\bm{D}}^{-1}{\bm{y}}\|_{2}^{2}+\alpha\lambda\|{\bm{y}}\|_{2}$, by
definition of a subgradient, the last inequality is equivalent to
${\bm{D}}^{-1}{\bm{x}}\in\partial H({\bm{0}})$, where $\partial H({\bm{0}})$
is the subdifferential of $H$ at ${\bm{0}}$:

$\partial
H({\bm{0}})=\\{{\bm{y}}\in{\mathbb{R}}^{n}\mid\|{\bm{y}}\|_{2}\leq\alpha\lambda\\}.$

This shows that ${\bm{0}}$ is a global minimizer of $F$ iff
$\|{\bm{D}}^{-1}{\bm{x}}\|_{2}\leq\alpha\lambda$; in other words,

$\operatorname{prox}_{\alpha
h\circ{\bm{D}}}({\bm{x}})={\bm{0}}\qquad\Leftrightarrow\qquad\|{\bm{D}}^{-1}{\bm{x}}\|_{2}\leq\alpha\lambda.$

Now if the minimizer ${\bm{z}}^{\star}\neq{\bm{0}}$ of $F$ is non-zero, we can
use the necessary and sufficient first-order condition of optimality for $F$
($F$ being convex and differentiable at ${\bm{z}}^{\star}$ for
${\bm{z}}^{\star}\neq{\bm{0}}$) to obtain

${\bm{z}}^{\star}-{\bm{x}}+\alpha\lambda\frac{{\bm{D}}^{2}{\bm{z}}^{\star}}{\|{\bm{D}}{\bm{z}}^{\star}\|_{2}}={\bm{0}}\qquad\Leftrightarrow\qquad{\bm{z}}^{\star}=\left[{\bm{I}}+\alpha\lambda\frac{{\bm{D}}^{2}}{\|{\bm{D}}{\bm{z}}^{\star}\|_{2}}\right]^{-1}{\bm{x}}$

Let us call $\theta^{\star}\triangleq\|{\bm{D}}{\bm{z}}^{\star}\|_{2}>0$. We
can rewrite the coordinates $z_{i}^{\star}$ of ${\bm{z}}^{\star}$ more
explicitly, as a function of $x_{i}$ and $\theta^{\star}$

$\big{[}\operatorname{prox}_{\alpha
h\circ{\bm{D}}}({\bm{x}})\big{]}_{i}\triangleq
z_{i}^{\star}=\frac{\theta^{\star}x_{i}}{\theta^{\star}+\alpha\lambda
d_{i}^{2}}.$

Although we have found an expression for ${\bm{z}}^{\star}$, the constant
$\theta^{\star}$ remains a function of ${\bm{z}}^{\star}$. This introduces
some constraints on $\theta^{\star}$, which in turn will introduce constraints
on ${\bm{x}}$ to guarantee the existence of $\theta^{\star}$. For example, we
have

$\theta^{\star
2}=\|{\bm{D}}{\bm{z}}^{\star}\|_{2}^{2}=\sum_{i=1}^{n}d_{i}^{2}z_{i}^{\star
2}=\theta^{\star
2}\sum_{i=1}^{n}\left[\frac{d_{i}x_{i}}{\theta^{\star}+\alpha\lambda
d_{i}^{2}}\right]^{2}\qquad\Leftrightarrow\qquad G(\theta^{\star})=1,$

which is the condition on $\theta^{\star}$ in Eq. 17. Moreover, the function
$G$ is convex and monotonically decreasing on the positive line
${\mathbb{R}}_{+}$. This implies that a solution $\theta^{\star}>0$ of
$G(\theta^{\star})=1$ exists (and is unique) if and only if $G(0)>1$

$G(0)=\frac{\|{\bm{D}}^{-1}{\bm{x}}\|_{2}^{2}}{(\alpha\lambda)^{2}}>1\qquad\Leftrightarrow\qquad\|{\bm{D}}^{-1}{\bm{x}}\|_{2}>\alpha\lambda.$

In other words, we can conclude that

$\big{[}\operatorname{prox}_{\alpha
h\circ{\bm{D}}}({\bm{x}})\big{]}_{i}=\frac{\theta^{\star}x_{i}}{\theta^{\star}+\alpha\lambda
d_{i}^{2}}\qquad\Leftrightarrow\qquad\|{\bm{D}}^{-1}{\bm{x}}\|_{2}>\alpha\lambda.$

∎

Finally, Theorem 1 is a corollary of the above proposition. Let’s first recall
the theorem:

######  Theorem 1 (Weighted proximal operator of $\ell_{2}$).

Let ${\bm{D}}=\operatorname{diag}(d_{1},\ldots,d_{n})$ be a positive definite
diagonal matrix (i.e. $d_{i}>0$ for all $i$), and
$h({\bm{x}})=\lambda\|{\bm{x}}\|_{2}$ the $\ell_{2}$ penalty. The weighted
proximal operator of $h$ is given by

$\big{[}\operatorname{prox}_{\alpha
h}^{{\bm{D}}}({\bm{x}})\big{]}_{i}=\left\\{\begin{array}[]{cl}\dfrac{d_{i}\theta^{\star}x_{i}}{d_{i}\theta^{\star}+\alpha\lambda}&\textrm{if
$\|{\bm{D}}{\bm{x}}\|_{2}>\alpha\lambda$}\\\\[10.00002pt]
0&\textrm{otherwise},\end{array}\right.$

where $\theta^{\star}>0$ is the unique positive solution of

$\sum_{i=1}^{n}\left[\frac{d_{i}x_{i}}{d_{i}\theta^{\star}+\alpha\lambda}\right]^{2}=1.$

###### Proof.

We get the expected result using the changes of variable
${\bm{D}}\rightarrow{\bm{D}}^{-1/2}$ and
${\bm{x}}\rightarrow{\bm{D}}^{1/2}{\bm{x}}$ in Proposition 2, together with
Lemma 1 (pre-multiplication by ${\bm{D}}^{-1/2}$). ∎

###  A.2 Approximation of the weighted proximal operator of $\ell_{2}$

In this section, we prove Eq. 11, which is recalled here:

######  Proposition 11.

Let $\varepsilon>0$, and $\tilde{\theta}>0$ such that
$|\tilde{\theta}-\theta^{\star}|\leq\varepsilon$ (recall that $\theta^{\star}$
is defined as the unique positive solution of Eq. 10). Let
$\widetilde{\operatorname{prox}}^{{\bm{D}}}_{\alpha h}({\bm{x}})$ be an
approximation of the weighted proximal operator of $\ell_{2}$ (i.e.
$h({\bm{x}})=\lambda\|{\bm{x}}\|_{2}$), replacing $\theta^{\star}$ by
$\tilde{\theta}$ in Eq. 9. Then we have for all ${\bm{x}}\in{\mathbb{R}}^{n}$:

$\big{\|}\widetilde{\operatorname{prox}}_{\alpha
h}^{{\bm{D}}}({\bm{x}})-\operatorname{prox}_{\alpha
h}^{{\bm{D}}}({\bm{x}})\big{\|}_{2}\leq\varepsilon.$ (18)

###### Proof.

Let ${\bm{D}}=\operatorname{diag}(d_{1},\ldots,d_{n})$ be a positive definite
diagonal matrix, and $d_{\min}>0$ the smallest value of the diagonal of
${\bm{D}}$. For all ${\bm{x}}\in{\mathbb{R}}^{n}$

$\displaystyle\big{\|}\widetilde{\operatorname{prox}}_{\alpha
h}^{{\bm{D}}}({\bm{x}})-\operatorname{prox}_{\alpha
h}^{{\bm{D}}}({\bm{x}})\big{\|}_{2}^{2}$
$\displaystyle\leq\sum_{i=1}^{n}\left[\frac{d_{i}\tilde{\theta}x_{i}}{d_{i}\tilde{\theta}+\alpha\lambda}-\frac{d_{i}\theta^{\star}x_{i}}{d_{i}\theta^{\star}+\alpha\lambda}\right]^{2}$
$\displaystyle=\sum_{i=1}^{n}\left[\frac{\alpha\lambda
d_{i}x_{i}(\tilde{\theta}-\theta^{\star})}{(d_{i}\tilde{\theta}+\alpha\lambda)(d_{i}\theta^{\star}+\alpha\lambda)}\right]^{2}$
$\displaystyle\leq(\tilde{\theta}-\theta^{\star})^{2}\underbrace{\vphantom{\sum_{i=1}^{n}}\left[\frac{\alpha\lambda}{d_{\min}\tilde{\theta}+\alpha\lambda}\right]^{2}}_{\leq\,1}\underbrace{\sum_{i=1}^{n}\left[\frac{d_{i}x_{i}}{d_{i}\theta^{\star}+\alpha\lambda}\right]^{2}}_{=\,1\
\textrm{(\lx@cref{creftype~refnum}{eq:weighted-prox-l2-theta}})}$
$\displaystyle\leq(\tilde{\theta}-\theta^{\star})^{2}\leq\varepsilon^{2}$

∎

The following proposition also gives bounds on $\theta^{\star}$ defined in Eq.
10, to narrow down the search space for the numerical solver (such as the
Newton-Raphson algorithm, see Section C.2 for details).

######  Proposition 3 (Bounds on $\theta^{\star}$ for the $\ell_{1}/\ell_{2}$
penalty).

Let ${\bm{D}}=\operatorname{diag}(d_{1},\ldots,d_{n})$ be a positive definite
diagonal matrix, with $d_{\min}$ and $d_{\max}$ being respectively the
smallest and largest values of the diagonal of ${\bm{D}}$. Let
$\theta^{\star}$ defined by Eq. 10. Then for all ${\bm{x}}$ such that
$\|{\bm{D}}{\bm{x}}\|_{2}>\alpha\lambda$:

$0<\frac{\|{\bm{D}}{\bm{x}}\|_{2}-\alpha\lambda}{d_{\max}}\leq\theta^{\star}\leq\frac{\|{\bm{D}}{\bm{x}}\|_{2}-\alpha\lambda}{d_{\min}}.$
(19)

###### Proof.

Using the fact that for all $i$ we have $d_{\min}\leq d_{i}\leq d_{\max}$, we
get the following inequalities:

$\left[\frac{\|{\bm{D}}{\bm{x}}\|_{2}}{d_{\max}\theta^{\star}+\alpha\lambda}\right]^{2}\leq\underbrace{\sum_{i=1}^{n}\left[\frac{d_{i}x_{i}}{d_{i}\theta^{\star}+\alpha\lambda}\right]^{2}}_{=\,1\
\textrm{(\lx@cref{creftype~refnum}{eq:weighted-
prox-l2-theta})}}\leq\left[\frac{\|{\bm{D}}{\bm{x}}\|_{2}}{d_{\min}\theta^{\star}+\alpha\lambda}\right]^{2}$

These two inequalities give us the expected bounds on $\theta^{\star}$:

$\displaystyle\frac{\|{\bm{D}}{\bm{x}}\|_{2}}{d_{\max}\theta^{\star}+\alpha\lambda}$
$\displaystyle\leq 1$ $\displaystyle\Leftrightarrow$
$\displaystyle\frac{\|{\bm{D}}{\bm{x}}\|_{2}-\alpha\lambda}{d_{\max}}$
$\displaystyle\leq\theta^{\star}$
$\displaystyle\frac{\|{\bm{D}}{\bm{x}}\|_{2}}{d_{\min}\theta^{\star}+\alpha\lambda}$
$\displaystyle\geq 1$ $\displaystyle\Leftrightarrow$
$\displaystyle\frac{\|{\bm{D}}{\bm{x}}\|_{2}-\alpha\lambda}{d_{\min}}$
$\displaystyle\geq\theta^{\star}$

∎

###  A.3 Weighted proximal operator of MCP / $\ell_{2}$

Similar to Section A.1, we start by giving the expression of the proximal
operator of $h\circ{\bm{A}}^{-1/2}$ (up to a change of variable, for
simplicity), in order to prove Theorem 2 via Lemma 1.

######  Proposition 4.

Let ${\bm{D}}=\operatorname{diag}(d_{1},\ldots,d_{n})$ be a positive definite
diagonal matrix (i.e. $d_{i}>0$ for all $i$), and
$h({\bm{x}}\,;\beta,\lambda)$ the MCP / $\ell_{2}$ penalty defined in Eq. 4.
Suppose that $\alpha$ and $\beta$ satisfy $\beta>\alpha d_{\max}^{2}$, where
$d_{\max}$ is the largest value of the diagonal of ${\bm{D}}$. The proximal
operator of $h\circ{\bm{D}}$ is given by

$\big{[}\operatorname{prox}_{\alpha
h\circ{\bm{D}}}({\bm{x}})\big{]}_{i}=\left\\{\begin{array}[]{cl}x_{i}&\textrm{if
$\|{\bm{D}}{\bm{x}}\|_{2}>\beta\lambda$}\\\\[8.00003pt]
\dfrac{\beta\theta^{\star}x_{i}}{\beta\theta^{\star}+\alpha
d_{i}^{2}(\beta\lambda-\theta^{\star})}&\textrm{\begin{tabular}[]{@{}l}if
$\|{\bm{D}}{\bm{x}}\|_{2}\leq\beta\lambda$\\\\[3.00003pt] and
$\|{\bm{D}}^{-1}{\bm{x}}\|_{2}>\alpha\lambda$\end{tabular}}\\\\[15.00002pt]
0&\textrm{\begin{tabular}[]{@{}l}if
$\|{\bm{D}}{\bm{x}}\|_{2}\leq\beta\lambda$\\\\[3.00003pt] and
$\|{\bm{D}}^{-1}{\bm{x}}\|_{2}\leq\alpha\lambda$,\end{tabular}}\end{array}\right.$
(27)

where $\theta^{\star}>0$ is the unique positive solution of

$G(\theta^{\star})\triangleq\beta^{2}\sum_{i=1}^{n}\left[\frac{d_{i}x_{i}}{\theta^{\star}(\beta-\alpha
d_{i}^{2})+\alpha\beta\lambda d_{i}^{2}}\right]^{2}=1.$ (28)

###### Proof.

Let us first recall the definition of the proximal operator of
$h\circ{\bm{D}}$

$\operatorname{prox}_{\alpha
h\circ{\bm{D}}}({\bm{x}})\triangleq\operatorname*{argmin}_{{\bm{z}}\in{\mathbb{R}}^{n}}\underbrace{\frac{1}{2}\|{\bm{z}}-{\bm{x}}\|_{2}^{2}+\alpha
h({\bm{D}}{\bm{z}}\,;\beta,\lambda)}_{\triangleq\,F({\bm{z}})}.$

Let us first consider the case where a global minimizer ${\bm{z}}^{\star}$ of
$F$ satisfies $\|{\bm{D}}{\bm{z}}^{\star}\|_{2}>\beta\lambda$. By definition
of the MCP / $\ell_{2}$ penalty in Eq. 4, this means that

$F({\bm{z}}^{\star})=\operatorname*{argmin}_{{\bm{z}}\in{\mathbb{R}}^{n}}\frac{1}{2}\|{\bm{z}}-{\bm{x}}\|_{2}^{2}+\alpha\frac{\beta\lambda^{2}}{2}.$

Since the condition $\|{\bm{D}}{\bm{z}}^{\star}\|_{2}>\beta\lambda$ implies
that ${\bm{z}}^{\star}\neq{\bm{0}}$, we can use the necessary and sufficient
first-order condition of optimality for $F$ (which is convex and
differentiable at ${\bm{z}}^{\star}\neq{\bm{0}}$) to get
${\bm{z}}^{\star}={\bm{x}}$. In other words, we have

$\operatorname{prox}_{\alpha
h\circ{\bm{D}}}({\bm{x}})={\bm{x}}\qquad\Leftrightarrow\qquad\|{\bm{D}}{\bm{x}}\|_{2}>\beta\lambda.$

Otherwise, if a global minimizer ${\bm{z}}^{\star}$ of $F$ satisfies
$\|{\bm{D}}{\bm{z}}^{\star}\|_{2}\leq\beta\lambda$, then again by definition
of $h$:

$F({\bm{z}}^{\star})=\operatorname*{argmin}_{{\bm{z}}\in{\mathbb{R}}^{n}}\frac{1}{2}\|{\bm{z}}-{\bm{x}}\|_{2}^{2}+\alpha\left(\lambda\|{\bm{D}}{\bm{z}}\|_{2}-\frac{\|{\bm{D}}{\bm{z}}\|_{2}^{2}}{2\beta}\right).$

Let us first consider the special case where ${\bm{0}}$ is a minimizer of $F$.
Similar to the proof of Proposition 2, by definition of ${\bm{0}}$ being a
global minimizer of $F$, we have

$\displaystyle\forall{\bm{z}}\in{\mathbb{R}}^{n},\ \ F({\bm{z}})\geq
F({\bm{0}})$ $\displaystyle\Leftrightarrow$
$\displaystyle\forall{\bm{z}}\in{\mathbb{R}}^{n},\ \
\frac{1}{2}\|{\bm{z}}\|_{2}^{2}+\alpha\left(\lambda\|{\bm{D}}{\bm{z}}\|_{2}-\frac{\|{\bm{D}}{\bm{z}}\|_{2}^{2}}{2\beta}\right)\geq\langle{\bm{z}},{\bm{x}}\rangle$
$\displaystyle\Leftrightarrow$
$\displaystyle\forall{\bm{y}}\in{\mathbb{R}}^{n},\ \
\frac{1}{2}\|{\bm{D}}^{-1}{\bm{y}}\|_{2}^{2}-\alpha\frac{\|{\bm{y}}\|_{2}^{2}}{2\beta}+\alpha\lambda\|{\bm{y}}\|_{2}\geq\langle{\bm{y}},{\bm{D}}^{-1}{\bm{x}}\rangle.$

The last inequality is equivalent to ${\bm{D}}^{-1}{\bm{x}}\in\partial
H({\bm{0}})$, where $\partial H({\bm{0}})$ is the subdifferential of the
function $H$ defined by:

$\displaystyle
H({\bm{y}})\triangleq\frac{1}{2}\|{\bm{D}}^{-1}{\bm{y}}\|_{2}^{2}-\alpha\frac{\|{\bm{y}}\|_{2}^{2}}{2\beta}+\alpha\lambda\|{\bm{y}}\|_{2}$
and $\displaystyle\partial
H({\bm{0}})=\\{{\bm{y}}\in{\mathbb{R}}^{n}\mid\|{\bm{y}}\|_{2}\leq\alpha\lambda\\}.$

Note that the function $H$ is convex thanks to the condition $\beta>\alpha
d_{\max}^{2}$, and therefore its subdifferential is well defined. This proves
that ${\bm{0}}$ is a global minimizer of $F$ iff
$\|{\bm{D}}^{-1}{\bm{x}}\|_{2}\leq\alpha\lambda$ (and
$\|{\bm{D}}{\bm{x}}\|_{2}\leq\beta\lambda$, given the first case considered
above). In other words,

$\operatorname{prox}_{\alpha
h\circ{\bm{D}}}({\bm{x}})={\bm{0}}\qquad\Leftrightarrow\qquad\|{\bm{D}}{\bm{x}}\|_{2}\leq\beta\lambda\
\ \textrm{and}\ \ \|{\bm{D}}^{-1}{\bm{x}}\|_{2}\leq\alpha\lambda.$

If the minimizer ${\bm{z}}^{\star}\neq{\bm{0}}$ of $F$ is non-zero and
satisfies $\|{\bm{D}}{\bm{z}}^{\star}\|_{2}\leq\beta\lambda$, then we can
again use the necessary and sufficient first-order condition of optimality for
$F$ to get

${\bm{z}}^{\star}-{\bm{x}}+\alpha\left(\frac{\lambda{\bm{D}}^{2}{\bm{z}}^{\star}}{\|{\bm{D}}{\bm{z}}^{\star}\|_{2}}-\frac{{\bm{D}}^{2}{\bm{z}}^{\star}}{\beta}\right)={\bm{0}}\qquad\Leftrightarrow\qquad{\bm{z}}^{\star}=\left[{\bm{I}}+\alpha\left(\frac{\lambda}{\|{\bm{D}}{\bm{z}}^{\star}\|_{2}}-\frac{1}{\beta}\right){\bm{D}}^{2}\right]^{-1}{\bm{x}}.$

Let us call $\theta^{\star}\triangleq\|{\bm{D}}{\bm{z}}^{\star}\|_{2}>0$. We
can rewrite the coordinates $z_{i}^{\star}$ of ${\bm{z}}^{\star}$ more
explicitly, as a function of $x_{i}$ and $\theta^{\star}$

$\big{[}\operatorname{prox}_{\alpha
h\circ{\bm{D}}}({\bm{x}})\big{]}_{i}\triangleq
z_{i}^{\star}=\frac{\beta\theta^{\star}x_{i}}{\beta\theta^{\star}+\alpha
d_{i}^{2}(\beta\lambda-\theta^{\star})}.$

Similar to the proof of Proposition 2, the constraints on $\theta^{\star}$
induce constraints on ${\bm{x}}$ to guarantee the existence of
$\theta^{\star}$. For example, we have

$\theta^{\star
2}=\|{\bm{D}}{\bm{z}}^{\star}\|_{2}^{2}=\sum_{i=1}^{n}d_{i}^{2}z_{i}^{\star
2}=\beta^{2}\theta^{\star
2}\sum_{i=1}^{n}\frac{d_{i}^{2}x_{i}^{2}}{\big{(}\theta^{\star}(\beta-\alpha
d_{i}^{2})+\alpha\beta\lambda d_{i}^{2}\big{)}^{2}}\qquad\Leftrightarrow\qquad
G(\theta^{\star})=1,$

which is the condition on $\theta^{\star}$ in Eq. 28. Moreover, since we made
the assumption that $\beta>\alpha d_{\max}^{2}$, the function $G$ is convex
and monotonically decreasing on the positive line $\mathbb{R}_{+}$. This
implies that a solution $\theta^{\star}\in(0,\beta\lambda]$ of
$G(\theta^{\star})=1$ exists (recall that we are in the case where
$\|{\bm{D}}{\bm{z}}^{\star}\|_{2}\leq\beta\lambda$), and is unique, if and
only if $G(0)>1\geq G(\beta\lambda)$. In other words, this implies

$\displaystyle
G(0)=\frac{\|{\bm{D}}^{-1}{\bm{x}}\|_{2}^{2}}{(\alpha\lambda)^{2}}$
$\displaystyle>1$ $\displaystyle\Leftrightarrow$
$\displaystyle\|{\bm{D}}^{-1}{\bm{x}}\|_{2}$ $\displaystyle>\alpha\lambda$
$\displaystyle\textrm{and}\ \
G(\beta\lambda)=\frac{\|{\bm{D}}{\bm{x}}\|_{2}^{2}}{(\beta\lambda)^{2}}$
$\displaystyle\leq 1$ $\displaystyle\Leftrightarrow$
$\displaystyle\|{\bm{D}}{\bm{x}}\|_{2}$ $\displaystyle\leq\beta\lambda.$

∎

Finally, Theorem 2 is a corollary of the above proposition. Let’s first recall
the theorem:

######  Theorem 2 (Weighted proximal operator of MCP / $\ell_{2}$).

Let ${\bm{D}}=\operatorname{diag}(d_{1},\ldots,d_{n})$ be a positive definite
diagonal matrix (i.e. $d_{i}>0$ for all $i$), and
$h({\bm{x}})=\operatorname{MCP}(\|{\bm{x}}\|_{2})$ the MCP / $\ell_{2}$
penalty. Suppose that $\alpha$ and $\beta$ satisfy $\alpha<\beta d_{\min}$,
where $d_{\min}$ is the smallest value of the diagonal of ${\bm{D}}$. The
weighted proximal operator of $h$ is given by

$\big{[}\operatorname{prox}_{\alpha
h}^{{\bm{D}}}({\bm{x}})\big{]}_{i}=\left\\{\begin{array}[]{cl}x_{i}&\textrm{if
$\|{\bm{x}}\|_{2}>\beta\lambda$}\\\\[8.00003pt]
\dfrac{d_{i}\beta\theta^{\star}x_{i}}{(d_{i}\beta-\alpha)\theta^{\star}+\alpha\beta\lambda}&\textrm{\begin{tabular}[]{@{}l}if
$\|{\bm{x}}\|_{2}\leq\beta\lambda$\\\\[3.99994pt] and
$\|{\bm{D}}{\bm{x}}\|_{2}>\alpha\lambda$\end{tabular}}\\\\[15.00002pt]
0&\textrm{\begin{tabular}[]{@{}l}if
$\|{\bm{x}}\|_{2}\leq\beta\lambda$\\\\[3.99994pt] and
$\|{\bm{D}}{\bm{x}}\|_{2}\leq\alpha\lambda$,\end{tabular}}\end{array}\right.$

where $\theta^{\star}>0$ is the unique positive solution of

$\beta^{2}\sum_{i=1}^{n}\left[\frac{d_{i}x_{i}}{(d_{i}\beta-\alpha)\theta^{\star}+\alpha\beta\lambda}\right]^{2}=1.$

###### Proof.

We get the expected result using the changes of variables
${\bm{D}}\rightarrow{\bm{D}}^{-1/2}$ and
${\bm{x}}\rightarrow{\bm{D}}^{1/2}{\bm{x}}$ in Proposition 4, together with
Lemma 1 (pre-multiplication by ${\bm{D}}^{-1/2}).$ ∎

###  A.4 Approximation of the weighted proximal operator of MCP / $\ell_{2}$

Similar to Eq. 11, the following proposition shows that the error on the
weighted proximal operator of MCP / $\ell_{2}$ induced by the numerical solver
to find $\theta^{\star}$ is controlled by the error of the numerical solver.

######  Proposition 5.

Let $\varepsilon>0$, and $\tilde{\theta}>0$ such that
$|\tilde{\theta}-\theta^{\star}|\leq\varepsilon$ (recall that $\theta^{\star}$
is defined as the unique positive solution of Eq. 13). Let
$\widetilde{\operatorname{prox}}^{{\bm{D}}}_{\alpha h}({\bm{x}})$ be an
approximation of the weighted proximal operator of MCP / $\ell_{2}$ (i.e.
$h({\bm{x}})=\operatorname{MCP}(\|{\bm{x}}\|_{2};\beta,\lambda)$), replacing
$\theta^{\star}$ by $\tilde{\theta}$ in Eq. 12.

Suppose that $\alpha$ and $\beta$ satisfy $\alpha<d_{\min}\beta$, where
$d_{\min}$ is the smallest value of the diagonal of ${\bm{D}}$ (recall that
${\bm{D}}$ is a positive definite diagonal matrix). Then we have for all
${\bm{x}}\in{\mathbb{R}}^{n}$:

$\big{\|}\widetilde{\operatorname{prox}}_{\alpha
h}^{{\bm{D}}}({\bm{x}})-\operatorname{prox}_{\alpha
h}^{{\bm{D}}}({\bm{x}})\big{\|}_{2}\leq\varepsilon.$ (29)

###### Proof.

The proof is similar to the one of Eq. 18 in Section A.2. Let
${\bm{D}}=\operatorname{diag}(d_{1},\ldots,d_{n})$ be a positive definite
diagonal matrix, and $d_{\min}>0$ the smallest value of the diagonal of
${\bm{D}}$. For all ${\bm{x}}\in{\mathbb{R}}^{n}$

$\displaystyle\big{\|}\widetilde{\operatorname{prox}}_{\alpha
h}^{{\bm{D}}}({\bm{x}})-\operatorname{prox}_{\alpha
h}^{{\bm{D}}}({\bm{x}})\big{\|}_{2}^{2}$
$\displaystyle\leq\sum_{i=1}^{n}\left[\frac{d_{i}\beta\tilde{\theta}x_{i}}{(d_{i}\beta-\alpha)\tilde{\theta}+\alpha\beta\lambda}-\frac{d_{i}\beta\theta^{\star}x_{i}}{(d_{i}\beta-\alpha)\theta^{\star}+\alpha\beta\lambda}\right]^{2}$
$\displaystyle=\sum_{i=1}^{n}\left[\frac{\alpha\beta^{2}\lambda
d_{i}x_{i}(\tilde{\theta}-\theta^{\star})}{((d_{i}\beta-\alpha)\tilde{\theta}+\alpha\beta\lambda)((d_{i}\beta-\alpha)\theta^{\star}+\alpha\beta\lambda)}\right]^{2}$
$\displaystyle\leq(\tilde{\theta}-\theta^{\star})^{2}\underbrace{\vphantom{\sum_{i=1}^{2}}\left[\frac{\alpha\beta\lambda}{(d_{\min}\beta-\alpha)\tilde{\theta}+\alpha\beta\lambda}\right]^{2}}_{\leq\,1}\underbrace{\beta^{2}\sum_{i=1}^{n}\left[\frac{d_{i}x_{i}}{(d_{i}\beta-\alpha)\theta^{\star}+\alpha\beta\lambda}\right]^{2}}_{=\,1\
\textrm{(\lx@cref{creftype~refnum}{eq:weighted-prox-mcp-l2-theta})}}$
$\displaystyle\leq(\tilde{\theta}-\theta^{\star})^{2}\leq\varepsilon^{2}$

∎

Finally, the following proposition gives bounds on $\theta^{\star}$ defined in
Eq. 13, to narrow down the search space for the numerical solver (e.g. Newton-
Raphson algorithm).

######  Proposition 6 (Bounds on $\theta^{\star}$ for the MCP / $\ell_{2}$
penalty).

Let ${\bm{D}}=\operatorname{diag}(d_{1},\ldots,d_{n})$ be a positive definite
diagonal matrix, with $d_{\min}$ and $d_{\max}$ being respectively the
smallest and largest values of the diagonal of ${\bm{D}}$. Let
$\theta^{\star}$ defined by Eq. 13, and suppose that $\alpha$ and $\beta$
satisfy $\alpha<d_{\min}\beta$. Then for all ${\bm{x}}$ such that
$\|{\bm{x}}\|_{2}\leq\beta\lambda$ and
$\|{\bm{D}}{\bm{x}}\|_{2}>\alpha\lambda$

$0<\beta\frac{\|{\bm{D}}{\bm{x}}\|_{2}-\alpha\lambda}{d_{\max}\beta-\alpha}\leq\theta^{\star}\leq\beta\frac{\|{\bm{D}}{\bm{x}}\|_{2}-\alpha\lambda}{d_{\min}\beta-\alpha}.$
(30)

###### Proof.

The proof is similar to the one of Eq. 19 in Section A.2. Using the fact that
for all $i$ we have $d_{\min}\leq d_{i}\leq d_{\max}$, we get the following
inequalities:

$\left[\frac{\beta\|{\bm{D}}{\bm{x}}\|_{2}}{(d_{\max}\beta-\alpha)\theta^{\star}+\alpha\beta\lambda}\right]^{2}\leq\underbrace{\beta^{2}\sum_{i=1}^{n}\left[\frac{d_{i}x_{i}}{(d_{i}\beta-\alpha)\theta^{\star}+\alpha\beta\lambda}\right]^{2}}_{=\,1\
\textrm{(\lx@cref{creftype~refnum}{eq:weighted-prox-
mcp-l2-theta})}}\leq\left[\frac{\beta\|{\bm{D}}{\bm{x}}\|_{2}}{(d_{\min}\beta-\alpha)\theta^{\star}+\alpha\beta\lambda}\right]^{2}$

These two inequalities give us the expected bounds on $\theta^{\star}$:

$\displaystyle\frac{\beta\|{\bm{D}}{\bm{x}}\|_{2}}{(d_{\max}\beta-\alpha)\theta^{\star}+\alpha\beta\lambda}$
$\displaystyle\leq 1$ $\displaystyle\Leftrightarrow$
$\displaystyle\beta\frac{\|{\bm{D}}{\bm{x}}\|_{2}-\alpha\lambda}{d_{\max}\beta-\alpha}$
$\displaystyle\leq\theta^{\star}$
$\displaystyle\frac{\beta\|{\bm{D}}{\bm{x}}\|_{2}}{(d_{\min}\beta-\alpha)\theta^{\star}+\alpha\beta\lambda}$
$\displaystyle\geq 1$ $\displaystyle\Leftrightarrow$
$\displaystyle\beta\frac{\|{\bm{D}}{\bm{x}}\|_{2}-\alpha\lambda}{d_{\min}\beta-\alpha}$
$\displaystyle\geq\theta^{\star}$

∎

##  Appendix B Proof of the convergence analysis

In this section, we prove the main theorem of convergence from Section 5,
extending the analysis from Xu et al. (2019); Yun et al. (2020) to inexact
proximal gradient methods, under 2 on the approximation of the weighted
proximal operator. Before getting the the proof of Theorem 3, we first recall
the definition of the Frechet subdifferential, which plays a central role in
the statement of convergence for non-differentiable and non-convex functions.

######  Definition 1 (Frechet subdifferential; Rockafellar & Wets 2009).

Let $F:{\mathbb{R}}^{N}\rightarrow\overline{{\mathbb{R}}}$ be a function
(possibly non-convex), and a point $\overline{{\bm{x}}}$ such that
$F(\overline{{\bm{x}}})$ is finite. The _Frechet subdifferential_ is the set
$\widehat{\partial}F$ of regular subgradients, satisfying:

$\widehat{\partial}F(\overline{{\bm{x}}})=\bigg{\\{}{\bm{v}}\in{\mathbb{R}}^{N}\;\Big{|}\;\liminf_{\begin{subarray}{c}{\bm{x}}\rightarrow\overline{{\bm{x}}}\\\
{\bm{x}}\neq\overline{{\bm{x}}}\end{subarray}}\frac{F({\bm{x}})-F(\overline{{\bm{x}}})-\langle{\bm{v}},{\bm{x}}-\overline{{\bm{x}}}\rangle}{{\bm{x}}-\overline{{\bm{x}}}}\geq
0\bigg{\\}}.$ (31)

Recall that
$F({\bm{x}})\triangleq\mathbb{E}_{\xi}[f({\bm{x}};\xi)]+h({\bm{x}})$. We are
now ready to prove Theorem 3:

######  Theorem 3.

Suppose that 1 and 2 are satisfied. If we run Algorithm 2 with a non-
increasing step-size $\alpha_{t}$, such that $\alpha_{0}<\delta/2L$, then the
output ${\bm{x}}_{R}$ of Algorithm 2, where $R$ is sampled uniformly in
$\\{1,\ldots,T\\}$, satisfies

$\displaystyle\mathbb{E}_{R}\big{[}\operatorname{dist}\\!\big{(}{\bm{0}},\widehat{\partial}F({\bm{x}}_{R})\big{)}^{2}\big{]}\leq\frac{C_{1}}{T}\sum_{t=0}^{T-1}\|{\bm{g}}_{t}-\nabla
f({\bm{x}}_{t})\|_{2}^{2}+\frac{C_{2}\Delta}{T}+\frac{C_{3}}{T},$

with $C_{1}$, $C_{2}$, and $C_{3}$ positive constants independent of $T$. Here
$\operatorname{dist}({\bm{z}},S)$ is the distance of a set $S$ to a point
${\bm{z}}$, defined as the minimal distance of any point in $S$ to ${\bm{z}}$.

###### Proof.

The proof follows the same steps as the proof of the non-asymptotic
convergence analysis from Yun et al. (2020) and Xu et al. (2019). We detail
the full proof here, including the steps from (Yun et al., 2020), for
completeness. Recall from 2 the notation

${\bm{x}}_{t+1}^{\star}=\operatorname{prox}_{\alpha
h}^{{\bm{D}}_{t}}\\!\big{(}{\bm{x}}_{t}-\alpha{\bm{D}}^{-1}{\bm{m}}_{t}\big{)}.$

#### Difference of consecutive iterates

Let us first consider the case of $t$ fixed, where 2 (2) is satisfied; that
is, we assume that ${\bm{x}}_{t+1}$ is an $\varepsilon_{t+1}$-approximation of
the true proximal update ${\bm{x}}_{t+1}^{\star}$, and $h$ is
$L^{\prime}$-smooth in a $\varepsilon_{t+1}$-ball around
${\bm{x}}_{t+1}^{\star}$. By definition of the weighted proximal operator as a
minimizer in Eq. 8, we have ${\bm{x}}_{t+1}^{\star}$ satisfying

$\displaystyle\langle{\bm{m}}_{t},{\bm{x}}_{t+1}^{\star}\rangle+h({\bm{x}}_{t+1}^{\star})+\frac{1}{2\alpha_{t}}\|{\bm{x}}_{t+1}^{\star}-{\bm{x}}_{t}\|^{2}_{{\bm{D}}_{t}}\leq\langle{\bm{m}}_{t},{\bm{x}}_{t}\rangle+h({\bm{x}}_{t})$
$\displaystyle\Leftrightarrow\qquad$
$\displaystyle\langle{\bm{m}}_{t},{\bm{x}}_{t+1}^{\star}-{\bm{x}}_{t}\rangle+h({\bm{x}}_{t+1}^{\star})+\frac{1}{2\alpha_{t}}\|{\bm{x}}_{t+1}^{\star}-{\bm{x}}_{t}\|^{2}_{{\bm{D}}_{t}}\leq
h({\bm{x}}_{t})$ (32)

Moreover since we assumed that $h$ is $L^{\prime}$-smooth in a
$\varepsilon_{t+1}$-ball ${\mathcal{B}}$ around ${\bm{x}}_{t+1}^{\star}$, and
${\bm{x}}_{t+1}\in{\mathcal{B}}$:

$h({\bm{x}}_{t+1})\leq h({\bm{x}}_{t+1}^{\star})+\langle\nabla
h({\bm{x}}_{t+1}^{\star}),{\bm{x}}_{t+1}-{\bm{x}}_{t+1}^{\star}\rangle+\frac{L^{\prime}}{2}\|{\bm{x}}_{t+1}-{\bm{x}}_{t+1}^{\star}\|_{2}^{2}$
(33)

Adding Eq. 32 & Eq. 33, we get an inequality on the difference of $h$
evaluated at two consecutive iterates:

$\langle{\bm{m}}_{t},{\bm{x}}_{t+1}^{\star}-{\bm{x}}_{t}\rangle+\frac{1}{2\alpha_{t}}\|{\bm{x}}_{t+1}^{\star}-{\bm{x}}_{t}\|^{2}_{{\bm{D}}_{t}}-\langle\nabla
h({\bm{x}}_{t+1}^{\star}),{\bm{x}}_{t+1}-{\bm{x}}_{t+1}^{\star}\rangle-\frac{L^{\prime}}{2}\|{\bm{x}}_{t+1}-{\bm{x}}_{t+1}^{\star}\|_{2}^{2}\leq
h({\bm{x}}_{t})-h({\bm{x}}_{t+1})$ (34)

On the other hand, using the $L$-smoothness of $f$, we can derive the
following two inequalities

$\displaystyle f({\bm{x}}_{t+1})$ $\displaystyle\leq
f({\bm{x}}_{t+1}^{\star})+\langle\nabla
f({\bm{x}}_{t+1}^{\star}),{\bm{x}}_{t+1}-{\bm{x}}_{t+1}^{\star}\rangle+\frac{L}{2}\|{\bm{x}}_{t+1}-{\bm{x}}_{t+1}^{\star}\|_{2}^{2}$
$\displaystyle f({\bm{x}}_{t+1}^{\star})$ $\displaystyle\leq
f({\bm{x}}_{t})+\langle\nabla
f({\bm{x}}_{t}),{\bm{x}}_{t+1}^{\star}-{\bm{x}}_{t}\rangle+\frac{L}{2}\|{\bm{x}}_{t+1}^{\star}-{\bm{x}}_{t}\|_{2}^{2}$

which, once added together, leads to a similar inequality as the one in Eq.
34, involving the difference of $f$ evaluated at two consecutive iterates:

$f({\bm{x}}_{t+1})-f({\bm{x}}_{t})\leq\langle\nabla
f({\bm{x}}_{t+1}^{\star}),{\bm{x}}_{t+1}-{\bm{x}}_{t+1}^{\star}\rangle+\langle\nabla
f({\bm{x}}_{t}),{\bm{x}}_{t+1}^{\star}-{\bm{x}}_{t}\rangle+\frac{L}{2}\|{\bm{x}}_{t+1}-{\bm{x}}_{t+1}^{\star}\|_{2}^{2}+\frac{L}{2}\|{\bm{x}}_{t+1}^{\star}-{\bm{x}}_{t}\|_{2}^{2}$
(35)

Finally, subtracting Eq. 35 to Eq. 34, we get the following inequality on the
difference of $F$ evaluated at two consecutive iterates

$\displaystyle\langle{\bm{m}}_{t}-\nabla
f({\bm{x}}_{t}),{\bm{x}}_{t+1}^{\star}-{\bm{x}}_{t}\rangle-\langle\nabla
F({\bm{x}}_{t+1}^{\star}),{\bm{x}}_{t+1}-{\bm{x}}_{t+1}^{\star}\rangle+\frac{1}{2\alpha_{t}}\|{\bm{x}}_{t+1}^{\star}-{\bm{x}}_{t}\|_{{\bm{D}}_{t}}^{2}$
$\displaystyle\qquad-\frac{L+L^{\prime}}{2}\|{\bm{x}}_{t+1}-{\bm{x}}_{t+1}^{\star}\|_{2}^{2}-\frac{L}{2}\|{\bm{x}}_{t+1}^{\star}-{\bm{x}}_{t}\|_{2}^{2}\leq
F({\bm{x}}_{t})-F({\bm{x}}_{t+1}).$ (36)

#### Bound on $\|{\bm{x}}_{t+1}^{\star}-{\bm{x}}_{t}\|_{2}^{2}$

We will eventually use this inequality to obtain an upper-bound on
$\|{\bm{x}}_{t+1}^{\star}-{\bm{x}}_{t}\|_{2}^{2}$. But first, we would like to
bound some of the remaining terms in this inequality, starting with the first
term of the left-hand side of Eq. 36. Since we have assumed in 2 that the
series of $\varepsilon_{t+1}^{2}$ was convergent, it guarantees that
$\varepsilon_{t+1}$ is bounded by some constant $D^{\prime}$

$\displaystyle\sum_{t=0}^{+\infty}\varepsilon_{t+1}^{2}<+\infty\qquad$
$\displaystyle\Rightarrow\qquad\forall t\geq 0,\
\|{\bm{x}}_{t+1}-{\bm{x}}_{t+1}^{\star}\|_{2}\leq\varepsilon_{t+1}\leq
D^{\prime}$ $\displaystyle\Rightarrow\qquad\forall t\geq 0,\
\|{\bm{x}}_{t+1}^{\star}-{\bm{x}}_{t}\|_{2}\leq\|{\bm{x}}_{t+1}-{\bm{x}}_{t}\|_{2}+\|{\bm{x}}_{t+1}-{\bm{x}}_{t+1}^{\star}\|_{2}\leq
D+D^{\prime}$

As an aside, note that while this provides an inequality over
$\|{\bm{x}}_{t+1}^{\star}-{\bm{x}}_{t}\|_{2}$, the bound we will eventually
derive from Eq. 36 will be more appropriate for proving convergence. Using the
definition of ${\bm{m}}_{t}=(1-\rho_{t}){\bm{g}}_{t}+\rho_{t}{\bm{m}}_{t-1}$,
together with Young’s inequality ($|\langle{\bm{a}},{\bm{b}}\rangle|\leq
c^{2}/2\|{\bm{a}}\|_{2}^{2}+1/2c^{2}\|{\bm{b}}\|_{2}^{2}$, for any constant
$c\neq 0$), we get

(37)

In particular, the last inequality (Eq. 37) uses 1 to bound
$\|{\bm{g}}_{t}\|_{2}$, and Lemma 1 of (Yun et al., 2020) to bound
$\|{\bm{m}}_{t-1}\|_{2}$, both by $G$. We then use Young’s inequality again to
bound the second term on the left-hand side of Eq. 36:

$\displaystyle\big{|}\langle\nabla
F({\bm{x}}_{t+1}^{\star}),{\bm{x}}_{t+1}-{\bm{x}}_{t+1}^{\star}\rangle\big{|}$
$\displaystyle\leq\frac{1}{2C}\|\nabla
F({\bm{x}}_{t+1}^{\star})\|_{2}^{2}+\frac{C}{2}\|{\bm{x}}_{t+1}-{\bm{x}}_{t+1}^{\star}\|_{2}^{2}$
$\displaystyle\leq\frac{1}{2C}\operatorname{dist}\\!\big{(}{\bm{0}},\widehat{\partial}F({\bm{x}}_{t+1}^{\star})\big{)}^{2}+\frac{C}{2}\|{\bm{x}}_{t+1}-{\bm{x}}_{t+1}^{\star}\|_{2}^{2},$
(38)

where $C$ is a positive constant defined as

$C=\frac{3(2+1/\gamma^{2}+3L+4L^{2})}{\delta/2\alpha_{0}-L}>0.$

$C$ is indeed positive because we assumed that $\alpha_{0}<\delta/2L$.
Finally, since we have assumed in 1 that ${\bm{D}}_{t}\succeq\delta{\bm{I}}$
(in other words, ${\bm{D}}_{t}-\delta{\bm{I}}$ is positive semi-definite), we
can lower-bound $\|{\bm{x}}_{t+1}^{\star}-{\bm{x}}_{t}\|_{{\bm{D}}_{t}}$ by
its Euclidean norm

$\frac{1}{2\alpha_{t}}\|{\bm{x}}_{t+1}^{\star}-{\bm{x}}_{t}\|_{{\bm{D}}_{t}}^{2}\geq\frac{\delta}{2\alpha_{t}}\|{\bm{x}}_{t+1}^{\star}-{\bm{x}}_{t}\|_{2}^{2}\geq\frac{\delta}{2\alpha_{0}}\|{\bm{x}}_{t+1}^{\star}-{\bm{x}}_{t}\|_{2}^{2},$

where the second inequality is due to the learning rate $\alpha_{t}$ being
non-increasing. Putting it all together:

Gathering all the terms involving
$\|{\bm{x}}_{t+1}^{\star}-{\bm{x}}_{t}\|_{2}^{2}$ on the left-hand side of the
inequality, and using 2, we finally get the following upper-bound:

$\displaystyle\left(\frac{\delta}{2\alpha_{0}}-L\right)\|{\bm{x}}_{t+1}^{\star}-{\bm{x}}_{t}\|_{2}^{2}$
$\displaystyle\leq
F({\bm{x}}_{t})-F({\bm{x}}_{t+1})+\frac{1}{2L}\|{\bm{g}}_{t}-\nabla
f({\bm{x}}_{t})\|_{2}^{2}+2\rho_{t}(D+D^{\prime})G$
$\displaystyle\qquad\qquad+\frac{1}{2C}\operatorname{dist}\\!\big{(}{\bm{0}},\widehat{\partial}F({\bm{x}}_{t+1}^{\star})\big{)}^{2}+\frac{L+L^{\prime}+C}{2}\varepsilon_{t+1}^{2}.$
(39)

Note that while we derived the above bound in the case where 2 (2) is
satisfied, this inequality is also valid when
${\bm{x}}_{t+1}={\bm{x}}_{t+1}^{\star}$; the terms coming from Eq. 38 would
vanish, making the above inequality looser, but valid for all $t$.

#### Bound on
$\operatorname{dist}\\!\big{(}{\bm{0}},\widehat{\partial}F({\bm{x}}_{t+1}^{\star})\big{)}^{2}$

Recall that by definition of the weighted proximal operator as a minimizer,
and by the first-order condition of optimality, we have
${\bm{x}}_{t+1}^{\star}$ satisfying

${\bm{0}}\in\frac{1}{\alpha_{t}}{\bm{D}}_{t}({\bm{x}}_{t+1}^{\star}-{\bm{x}}_{t})+{\bm{m}}_{t}+\widehat{\partial}h({\bm{x}}_{t+1}^{\star}).$

Adding $\nabla f({\bm{x}}_{t+1}^{\star})$, we have (Rockafellar & Wets, 2009;
Xu et al., 2019)

$\nabla
f({\bm{x}}_{t+1}^{\star})-{\bm{m}}_{t}-\frac{1}{\alpha_{t}}{\bm{D}}_{t}({\bm{x}}_{t+1}^{\star}-{\bm{x}}_{t})\in\nabla
f({\bm{x}}_{t+1}^{\star})+\widehat{\partial}h({\bm{x}}_{t+1}^{\star})=\widehat{\partial}F({\bm{x}}_{t+1}^{\star}).$
(40)

Since we know one element of $\widehat{\partial}F({\bm{x}}_{t+1}^{\star})$, we
can upper-bound its distance to ${\bm{0}}$:

(41)

Again, we would like to bound some of the remaining terms on the right-hand
side of the above inequality. We will start with the first term; using the
definition of ${\bm{m}}_{t}$, we have

$\displaystyle\|{\bm{m}}_{t}-\nabla f({\bm{x}}_{t+1}^{\star})\|_{2}^{2}$
$\displaystyle=\|(1-\rho_{t}){\bm{g}}_{t}-\nabla f({\bm{x}}_{t})+\nabla
f({\bm{x}}_{t})-\nabla
f({\bm{x}}_{t+1}^{\star})+\rho_{t}{\bm{m}}_{t-1}\|_{2}^{2}$ $\displaystyle\leq
4\Big{[}\|{\bm{g}}_{t}-\nabla f({\bm{x}}_{t})\|_{2}^{2}+\|\nabla
f({\bm{x}}_{t+1}^{\star})-\nabla
f({\bm{x}}_{t})\|_{2}^{2}+\|\rho_{t}{\bm{g}}_{t}\|_{2}^{2}+\|\rho_{t}{\bm{m}}_{t-1}\|_{2}^{2}\Big{]}$
$\displaystyle\leq 4\Big{[}\|{\bm{g}}_{t}-\nabla
f({\bm{x}}_{t})\|_{2}^{2}+L^{2}\|{\bm{x}}_{t+1}^{\star}-{\bm{x}}_{t}\|_{2}^{2}+2\rho_{t}^{2}G^{2}\Big{]},$
(42)

where we used the same bounds on $\|{\bm{g}}_{t}\|_{2}$ and
$\|{\bm{m}}_{t-1}\|_{2}$ as in Eq. 37. Similarly, we can bound the second term
on the right-hand side of Eq. 41:

(43)

Putting it all together, and using the bound on
$\|{\bm{x}}_{t+1}^{\star}-{\bm{x}}_{t}\|_{2}^{2}$ we derived earlier:

$\displaystyle\operatorname{dist}\\!\big{(}{\bm{0}},\widehat{\partial}F({\bm{x}}_{t+1}^{\star})\big{)}^{2}$
$\displaystyle\leq 3\Big{[}\Big{(}4+\frac{1}{L}\Big{)}\|{\bm{g}}_{t}-\nabla
f({\bm{x}}_{t})\|_{2}^{2}+\Big{(}2+\frac{1}{\gamma^{2}}+3L+4L^{2}\Big{)}\|{\bm{x}}_{t+1}^{\star}-{\bm{x}}_{t}\|_{2}^{2}$
$\displaystyle\qquad\qquad+4\rho_{t}(D+D^{\prime})G+8\rho_{t}^{2}G^{2}\Big{]}$
$\displaystyle\leq
C\big{(}F({\bm{x}}_{t})-F({\bm{x}}_{t+1})\big{)}+\Big{[}12+\frac{3}{L}+\frac{C}{2L}\Big{]}\|{\bm{g}}_{t}-\nabla
f({\bm{x}}_{t})\|_{2}^{2}+\frac{1}{2}\operatorname{dist}\\!\big{(}{\bm{0}},\widehat{\partial}F({\bm{x}}_{t+1}^{\star})\big{)}^{2}$
$\displaystyle\qquad\qquad+(2C+4)\rho_{t}(D+D^{\prime})G+8\rho_{t}^{2}G^{2}+\frac{C(L+L^{\prime}+C)}{2}\varepsilon_{t+1}^{2}$

Rearranging the terms together, we get the following bound on
$\operatorname{dist}\\!\big{(}{\bm{0}},\widehat{\partial}F({\bm{x}}_{t+1}^{\star})\big{)}^{2}$:

$\displaystyle\frac{1}{2}\operatorname{dist}\\!\big{(}{\bm{0}},\widehat{\partial}F({\bm{x}}_{t+1}^{\star})\big{)}^{2}$
$\displaystyle\leq
C\big{(}F({\bm{x}}_{t})-F({\bm{x}}_{t+1})\big{)}+\Big{[}12+\frac{3}{L}+\frac{C}{2L}\Big{]}\|{\bm{g}}_{t}-\nabla
f({\bm{x}}_{t})\|_{2}^{2}$
$\displaystyle\qquad\qquad+(2C+4)\rho_{t}(D+D^{\prime})G+8\rho_{t}^{2}G^{2}+\frac{C(L+L^{\prime}+C)}{2}\varepsilon_{t+1}^{2}.$
(44)

#### Bound on
$\operatorname{dist}\\!\big{(}{\bm{0}},\widehat{\partial}F({\bm{x}}_{t+1})\big{)}^{2}$

Although we have derived a bound on
$\operatorname{dist}\\!\big{(}{\bm{0}},\widehat{\partial}F({\bm{x}}_{t+1}^{\star})\big{)}^{2}$,
we are eventually interested in bounding the distance of
$\widehat{\partial}F({\bm{x}}_{t+1})$ to ${\bm{0}}$, at the possible
approximation ${\bm{x}}_{t+1}$ of ${\bm{x}}_{t+1}^{\star}$. Using the
triangular inequality

$\displaystyle\operatorname{dist}\\!\big{(}{\bm{0}},\widehat{\partial}F({\bm{x}}_{t+1})\big{)}^{2}$
$\displaystyle\leq\big{[}\operatorname{dist}\\!\big{(}{\bm{0}},\widehat{\partial}F({\bm{x}}_{t+1}^{\star})\big{)}+d_{H}\big{(}\widehat{\partial}F({\bm{x}}_{t+1}),\widehat{\partial}F({\bm{x}}_{t+1}^{\star})\big{)}\big{]}^{2}$
$\displaystyle\leq
2\operatorname{dist}\\!\big{(}{\bm{0}},\widehat{\partial}F({\bm{x}}_{t+1}^{\star})\big{)}^{2}+2d_{H}\big{(}\widehat{\partial}F({\bm{x}}_{t+1}),\widehat{\partial}F({\bm{x}}_{t+1}^{\star})\big{)}^{2},$

where $d_{H}$ is the Pompeiu–Hausdorff distance between two sets. In the case
of 2 (2), where $h$ (and therefore $F$) is differentiable in the
$\varepsilon_{t+1}$-ball ${\mathcal{B}}$ around ${\bm{x}}_{t+1}^{\star}$, both
of these sets are reduced to a singleton

Note that while we derived the inequality above with 2 (2), it is also valid
when ${\bm{x}}_{t+1}={\bm{x}}_{t+1}^{\star}$, albeit looser (since the
Pompeiu-Hausdorff distance would vanish in the latter case); hence this
inequality is valid for all $t$. Using the bound on
$\operatorname{dist}\\!\big{(}{\bm{0}},\widehat{\partial}F({\bm{x}}_{t+1}^{\star})\big{)}^{2}$,
we can finally get the following bound on
$\operatorname{dist}\\!\big{(}{\bm{0}},\widehat{\partial}F({\bm{x}}_{t+1})\big{)}^{2}$

$\displaystyle\operatorname{dist}\\!\big{(}{\bm{0}},\widehat{\partial}F({\bm{x}}_{t+1})\big{)}^{2}$
$\displaystyle\leq
4C\big{(}F({\bm{x}}_{t})-F({\bm{x}}_{t+1})\big{)}+4\Big{[}12+\frac{3}{L}+\frac{C}{2L}\Big{]}\|{\bm{g}}_{t}-\nabla
f({\bm{x}}_{t})\|_{2}^{2}$
$\displaystyle\qquad\qquad+4(2C+4)\rho_{t}(D+D^{\prime})G+32\rho_{t}^{2}G^{2}+\big{[}2C(L+L^{\prime}+C)+2(L+L^{\prime})^{2}\big{]}\varepsilon_{t+1}^{2}$
(45)

#### Convergence result

To prove our final convergence result, we can simply write the expectation
over squared distances (where $R$ is uniform) as an average of quantities we
have been capable of bounding:

$\mathbb{E}_{R}\big{[}\operatorname{dist}\\!\big{(}{\bm{0}},\widehat{\partial}F({\bm{x}}_{R})\big{)}^{2}\big{]}=\frac{1}{T}\sum_{t=0}^{T-1}\operatorname{dist}\\!\big{(}{\bm{0}},\widehat{\partial}F({\bm{x}}_{t+1})\big{)}^{2}$
(46)

Moreover, note that the right-hand side of the inequality in Eq. 45 involves
the difference of $F$ at two consecutive iterates, and therefore involves a
telescoping series once summed over. Using 1, we get

$\sum_{t=0}^{T-1}\big{(}F({\bm{x}}_{t})-F({\bm{x}}_{t+1})\big{)}=F({\bm{x}}_{0})-F({\bm{x}}_{T})\leq
F({\bm{x}}_{0})-F({\bm{x}}^{\star})\leq\Delta.$

We can similarly bound some of the remaining sums. Using the definition of
$\rho_{t}=\rho_{0}\mu^{t}$ (1), we have

$\displaystyle\sum_{t=0}^{T-1}\rho_{t}$
$\displaystyle\leq\rho_{0}\sum_{t=0}^{+\infty}\mu^{t}=\frac{\rho_{0}}{1-\mu}$
$\displaystyle\sum_{t=0}^{T-1}\rho_{t}^{2}$
$\displaystyle\leq\rho_{0}^{2}\sum_{t=0}^{+\infty}\mu^{2t}=\frac{\rho_{0}^{2}}{1-\mu^{2}},$

and using 2,

$\sum_{t=0}^{T-1}\varepsilon_{t+1}^{2}\leq\sum_{t=0}^{+\infty}\varepsilon_{t+1}^{2}=K.$

To conclude, using the bound from Eq. 45 inside Eq. 46, and the various bounds
on the sums above, we get the following bound on the expected squared
distance:

$\mathbb{E}_{R}\big{[}\operatorname{dist}\\!\big{(}{\bm{0}},\widehat{\partial}F({\bm{x}}_{R})\big{)}^{2}\big{]}\leq\frac{C_{1}}{T}\sum_{t=0}^{T-1}\|{\bm{g}}_{t}-\nabla
f({\bm{x}}_{t})\|_{2}^{2}+\frac{C_{2}\Delta}{T}+\frac{C_{3}}{T},$ (47)

where $C_{1}$, $C_{2}$, and $C_{3}$ are three constants independent of $T$,
defined as

$\displaystyle C_{1}$
$\displaystyle=4\Big{[}12+\frac{3}{L}+\frac{C}{2L}\Big{]}$ $\displaystyle
C_{2}$ $\displaystyle=4C$ $\displaystyle C_{3}$
$\displaystyle=\frac{4(2C+4)\rho_{0}(D+D^{\prime})G}{1-\mu}+\frac{32\rho_{0}^{2}G^{2}}{1-\mu^{2}}+K\big{[}2C(L+L^{\prime}+C)+2(L+L^{\prime})^{2}\big{]}.$

∎

We can also prove Corollary 1 (Xu et al., 2019), which we recall here:

######  Corollary 1 (Fixed mini-batch size).

If the assumptions of Theorem 3 are satisfied, with
$T=2(C_{2}\Delta+C_{3})/\varepsilon^{2}$ and with a fixed mini-batch size
$m_{t}$ with $m_{t}=2C_{1}\sigma^{2}/\varepsilon^{2}$, then the output
${\bm{x}}_{R}$ of Algorithm 2 satisfies

$\mathbb{E}_{R}\big{[}\operatorname{dist}\\!\big{(}{\bm{0}},\widehat{\partial}F({\bm{x}}_{R})\big{)}^{2}\big{]}\leq\varepsilon^{2},$

where $C_{1}$, $C_{2}$, and $C_{3}$ are the constants from Theorem 3. To have
$\mathbb{E}[\operatorname{dist}({\bm{0}},\widehat{\partial}F({\bm{x}}_{R}))]\leq\varepsilon$,
it is then sufficient to have $T=O(1/\varepsilon^{2})$, making the total
complexity $O(1/\varepsilon^{4})$.

###### Proof.

If we assume that 1 is satisfied, we know that ${\bm{g}}_{t}$ is estimated
using a mini-batch of $m_{t}$ samples of the form $\nabla
f({\bm{x}}_{t};\xi_{i_{t}})$, then we can further bound the remaining quantity
in Eq. 47:

Using Theorem 3, we then have

$\mathbb{E}_{R}\big{[}\operatorname{dist}\\!\big{(}{\bm{0}},\widehat{\partial}F({\bm{x}}_{R})\big{)}^{2}\big{]}\leq\frac{C_{1}\sigma^{2}}{m_{t}}+\frac{C_{2}\Delta}{T}+\frac{C_{3}}{T}.$
(48)

And finally substituting $T=2(C_{2}\Delta+C_{3})/\varepsilon^{2}$ and
$m_{t}=2C_{1}\sigma^{2}/\varepsilon^{2}$ in Eq. 48, we get the expected result

$\mathbb{E}_{R}\big{[}\operatorname{dist}\\!\big{(}{\bm{0}},\widehat{\partial}F({\bm{x}}_{R})\big{)}^{2}\big{]}\leq\varepsilon^{2}.$
(49)

∎

##  Appendix C Algorithmic details

In this section, we give some details about the different algorithms used in
this work, including the proximal gradient method ProxGen (Yun et al., 2020),
the Newton-Raphson algorithm in order to approximate the weighted proximal
operator, as well as the procedure used to prune neural networks with indirect
sparsity.

In Algorithms 1 and 2, we show a side-by-side comparison between ProxGen
(Algorithm 2; Yun et al., 2020) to minimize a composite loss function, and a
standard first-order adaptive method to minimize an objective function $f$.
The only difference is highlighted in red, where the weighted proximal
operator is applied to the gradient update.

Algorithm 1 Minimization of the loss function
$f({\bm{x}})$ with an adaptive optimizer

0: A loss function $\mathbb{E}_{\xi}[f({\bm{x}};\xi)]$

0: A sequence of learning rates $(\alpha_{t})_{t\geq 0}$, ${\bm{m}}_{0}$,
${\bm{D}}_{0}$

for $t=0,\ldots,T-1$ do

Draw a minibatch $\xi_{t}$

Update the mean estimate ${\bm{m}}_{t}$ with ${\bm{g}}_{t}=\nabla
f({\bm{x}}_{t};\xi_{t})$

Update the preconditioning matrix ${\bm{D}}_{t}$

Update the parameters:
${\bm{x}}_{t+1}\leftarrow{\bm{x}}_{t}-\alpha_{t}{\bm{D}}_{t}^{-1}{\bm{m}}_{t}$

end for

${\bm{x}}_{T}$, or ${\bm{x}}_{R}$, $R$ sampled uniformly in $\\{1,\ldots,T\\}$

Algorithm 2 Minimization of the composite loss function
$f({\bm{x}})+h({\bm{x}})$ with an adaptive optimizer (Yun et al., 2020)

0: A composite loss function $\mathbb{E}_{\xi}[f({\bm{x}};\xi)]+h({\bm{x}})$

0: A sequence of learning rates $(\alpha_{t})_{t\geq 0}$, ${\bm{m}}_{0}$,
${\bm{D}}_{0}$

for $t=0,\ldots,T-1$ do

Draw a minibatch $\xi_{t}$

Update the mean estimate ${\bm{m}}_{t}$ with ${\bm{g}}_{t}=\nabla
f({\bm{x}}_{t};\xi_{t})$

Update the preconditioning matrix ${\bm{D}}_{t}$

Update the parameters:
${\bm{x}}_{t+1}\leftarrow\operatorname{prox}_{\alpha_{t}h}^{{\bm{D}}_{t}}\\!\big{(}{\bm{x}}_{t}-\alpha_{t}{\bm{D}}_{t}^{-1}{\bm{m}}_{t}\big{)}$

end for

${\bm{x}}_{T}$, or ${\bm{x}}_{R}$, $R$ sampled uniformly in $\\{1,\ldots,T\\}$

###  C.1 Adaptive optimizers

Table 3 shows, for reference, the updates of the mean estimate ${\bm{m}}_{t}$
and the preconditioning matrix ${\bm{D}}_{t}$ for different standard adaptive
first-order methods.

Table 3: Standard adaptive optimizers used in Deep Learning, with the
corresponding updates for the mean estimate ${\bm{m}}_{t}$, and the
preconditioning matrix ${\bm{D}}_{t}$. Table adapted from (Melchior et al.,
2019).

Mean estimate Preconditioning matrix $\widehat{{\bm{m}}}_{t}$ ${\bm{m}}_{t}$
${\bm{v}}_{t}$ ${\bm{D}}_{t}$ SGD – ${\bm{g}}_{t}$ – ${\bm{I}}$ Momentum
$\mu\widehat{{\bm{m}}}_{t-1}+{\bm{g}}_{t}$ $\widehat{{\bm{m}}}_{t}$ –
${\bm{I}}$ (Rumelhart et al., 1986) AdaGrad – ${\bm{g}}_{t}$
${\bm{v}}_{t-1}+{\bm{g}}_{t}^{2}$ $\sqrt{{\bm{v}}_{t}}+\varepsilon$ (Duchi et
al., 2011) RMSprop – ${\bm{g}}_{t}$
$\beta{\bm{v}}_{t-1}+(1-\beta){\bm{g}}_{t}^{2}$
$\sqrt{{\bm{v}}_{t}+\varepsilon}$ (Tieleman & Hinton, 2012) Adam
$\beta_{1}\widehat{{\bm{m}}}_{t-1}+(1-\beta_{1}){\bm{g}}_{t}$
$\dfrac{\widehat{{\bm{m}}}_{t}}{1-\beta_{1}^{t}}$
$\beta_{2}{\bm{v}}_{t-1}+(1-\beta_{2}){\bm{g}}_{t}^{2}$
$\sqrt{\dfrac{{\bm{v}}_{t}}{1-\beta_{2}^{t}}}+\varepsilon$ (Kingma & Ba, 2015)

###  C.2 Approximation of the weighted proximal operator with Newton-Raphson

In Algorithm 3, we show how to use the Newton-Raphson algorithm as our routine
for finding $\theta^{\star}$, and therefore approximate the weighted proximal
operator of the $\ell_{2}$ norm. This procedure can be adapted to MCP /
$\ell_{2}$ in a straightforward way. Note that while this iterative procedure
has to be run at every gradient update during optimization, empirically this
induces a reasonable overhead, which can be controlled by the tolerance
$\varepsilon$; lower tolerance yields fewer iterations per gradient update.
See Section 6.4 for empirical evidence.

Algorithm 3 Approximation of the weighted proximal operator of the $\ell_{2}$
norm ($h({\bm{x}})=\lambda\|{\bm{x}}\|_{2}$) with Newton-Raphson

0: A point ${\bm{x}}$ such that $\|{\bm{D}}{\bm{x}}\|_{2}>\alpha\lambda$.

0: An error tolerance $\varepsilon$

Initialization:
$\theta\leftarrow(\|{\bm{D}}{\bm{x}}\|_{2}-\alpha\lambda)/d_{\max}$
$\triangleright$ Eq. 19

while $|G(\theta)|>\varepsilon$ do

$\theta\leftarrow\theta-\dfrac{G(\theta)}{G^{\prime}(\theta)}$

end while

$\big{[}\operatorname{prox}_{\alpha
h}^{{\bm{D}}}({\bm{x}})\big{]}_{i}=\dfrac{d_{i}\theta
x_{i}}{d_{i}\theta+\alpha\lambda}$ $\triangleright$ Theorem 1

Where $G$ is the function defined by

$G(\theta)=\sum_{i=1}^{n}\left[\frac{d_{i}x_{i}}{d_{i}\theta+\alpha\lambda}\right]^{2}-1$

###  C.3 Pruning with indirect sparsity

In Section 6.3, we argued that structured sparsity had benefit not only on the
layer where groups of variables were zeroed-out, but on neighboring layer as
well. Fig. 4 shows an illustration of this behaviour on a convolutional neural
network, inspired by Fig. 1, where structured sparsity has been applied
channel-wise.

[width=0.8]figures/indirect-sparsity-pruning

Figure 4: Pruning a convolutional neural network with indirect sparsity. This
illustration shows the groups of filters set to zero thanks to the structured
sparsity inducing penalty (green for the first layer, blue for the second).
The filters that can be further pruned in the first layer thanks to indirect
sparsity are shown in light-blue.

The algorithm to prune these filters operates from right to left. When a group
of filters are set to zero thanks to the structured sparsity inducing penalty,
the corresponding channel at the previous layer (e.g. second channel of Layer
2 in Fig. 4, shown in light-blue) does not contribute to the predictions of
the network anymore; here, the representation at Layer 3 is independent of the
channel in light-blue. Therefore the channel in light-blue can be ignored from
prior computations as well, meaning that the filters in light-blue can be set
to zero as well, without functionally affecting the neural network: these
filters were pruned _indirectly_ by the structure in the sparsity patterns in
Layer 2.

On top of this, the structured sparsity inducing penalty applied to the first
layer can also set some groups of variables to zero (here the group shown in
green); this procedure can then be run from the output layer, all the way to
the input layer. This leads to significantly fewer parameters: in Fig. 4, the
first layer only has 50% of non-zero parameters left, while only 25% were
pruned directly through the structured sparsity inducing penalty applied to
that layer (i.e. 75% group sparsity). Finally, note that while we showed the
effect of indirect sparsity on a convolutional neural network, this applies
similarly to linear layers with row-wise groups.