# Best Subset Selection with Efficient Primal-Dual Algorithm

Shaogang Ren, Guanhua Fang, Ping Li

Cognitive Computing Lab
Baidu Research
10900 NE 8th St. Bellevue, WA 98004, USA

{renshaogang, fanggh2018<EMAIL_ADDRESS>
###### Abstract

111The initial version of the paper was submitted in February 2020.

Best subset selection is considered the ‘gold standard’ for many sparse
learning problems. A variety of optimization techniques have been proposed to
attack this non-convex and NP-hard problem. In this paper, we investigate the
dual forms of a family of $\ell_{0}$-regularized problems. An efficient
primal-dual method has been developed based on the primal and dual problem
structures. By leveraging the dual range estimation along with the incremental
strategy, our algorithm potentially reduces redundant computation and improves
the solutions of best subset selection. Theoretical analysis and experiments
on synthetic and real-world datasets validate the efficiency and statistical
properties of the proposed solutions.

##  1 Introduction

Sparse learning is a standard approach to alleviate model over-fitting issues
when the feature dimension is larger than the number of training samples. With
a training set $\\{x_{i},y_{i}\\}_{i=1}^{n}$ where $x_{i}\in\mathbb{R}^{p}$ is
the sample feature and $y_{i}$ is the corresponding label, this paper focuses
on the following generalized best subset selection problem,

$\displaystyle\min_{\beta\in\mathbb{R}^{p}}F(\beta)=f(\beta)+\lambda_{0}||\beta||_{0},$
(1) where $\displaystyle
f(\beta)=\sum_{i=1}^{n}l(\beta^{\top}x_{i},y_{i})+\lambda_{1}||\beta||_{1}+\lambda_{2}||\beta||_{2}^{2}.$

Here $l(\cdot)$ is a convex function, $\beta$ is the model parameter, and
$\lambda_{0}$, $\lambda_{1}$ and $\lambda_{2}$ are hyper-parameters/tuning
parameters. It is well-known that an $\ell_{0}$ solver
($\lambda_{0}>0,\lambda_{1}=0$) has superior statistical properties when the
signal-to-noise ratio (SNR) is high, but it may suffer over-fitting issues
when SNR is low (Gamarnik and Zadik, 2017; Mazumder et al., 2022). The
continuous-shrinkage solvers e.g., ridge/Lasso
($\lambda_{0}=0,\lambda_{1}>0$), can perform better in this case compared with
$\ell_{0}$ solver (Mazumder et al., 2022; Hastie et al., 2017). Combinations
of there hyper-parameters may adjust the model to work well in different noise
levels. Mazumder et al. (2022); Hazimeh and Mazumder (2020) used ridge/Lasso
to improve $\ell_{0}$ solutions and achieve better or comparable solutions
with less nonzeros.

Solving the vanilla $\ell_{0}$ regularized problem is known to be NP-hard
(Natarajan, 1995). By leveraging the significant computational advances in
mixed integer optimization (MIO), Bertsimas et al. (2016) performed near
optimal solutions to a special case of problem (1), for $\lambda_{1}=0$ and
$\lambda_{2}=0$. This method scaled up solutions to cases where feature sizes
are much larger than what were considered possible in the community (Furnival
and Wilson, 1974; Hazimeh and Mazumder, 2020). Their approach can achieve
approximate optimality via dual bounds but with the cost of longer computation
times. Bertsimas and Van Parys (2020) showed that cutting plane methods for
subset selection can work well with mild sample correlations and a succinctly
large $n$.

Different from the soft regularized ridge/Lasso problem given by (1),
Iterative Hard Thresholding (IHT) (Blumensath and Davies, 2009; Foucart, 2011;
Yuan et al., 2014; Shen and Li, 2017; Yuan and Li, 2020) has often been used
to solve $k$-sparse problems (2),

$\displaystyle\min_{||\beta||_{0}\leq
k}\sum_{i=1}^{n}l(\beta^{\top}x_{i},y_{i})+\lambda_{2}||\beta||_{2}^{2}\ .$
(2)

In Blumensath and Davies (2009); Foucart (2011), the authors demonstrated that
IHT can be applied to compute the compressed sensing problem. IHT-based
approaches have been studied by many researchers in the context of sparse
learning problems (Yuan et al., 2014; Jain et al., 2014, 2016; Yuan and Li,
2020). IHT methods require a specific value of the features number ($k$) to
start the algorithm. Many other approaches have also been developed to tackle
the $\ell_{0}$ regularized problems (Mazumder and Radchenko, 2017; Mazumder et
al., 2022; Soussen et al., 2015; Bian and Chen, 2020; Dedieu et al., 2021;
Yang and Yu, 2019; Dong et al., 2015; Hazimeh et al., 2020; Zhu et al., 2020).

Apart from $\ell_{0}$ solvers, extremely efficient and optimized
$\ell_{1}$-regularization (Lasso) solvers can solve an entire regularization
path (with a hundred values of the tuning parameter) in usually less than a
second (Friedman et al., 2010). Screening and coordinate incremental
techniques (Fercoq et al., 2015; Ndiaye et al., 2015; Massias et al., 2018;
Ren et al., 2020) can further scale the solutions to large datasets. Compared
to popular efficient solvers for Lasso, it seems that the high computation
cost for using $\ell_{0}$ regularized models (Bertsimas et al., 2016) might
discourage practitioners from adopting global optimization-based solvers of
(1) to daily analysis applications (Hastie et al., 2017; Mazumder et al.,
2022; Hazimeh and Mazumder, 2020). However, it is known (Loh and Wainwright,
2017; Hazimeh and Mazumder, 2020) that there is a significant gap in the
statistical quality of solutions that can be achieved via Lasso (and its
variants) and near-optimal solutions to non-convex subset-selection type
procedures. The choice of algorithm can significantly affect the quality of
solutions obtained. On many instances, algorithms that do a better job in
optimizing the non-convex subset-selection criterion (1) result in superior-
quality statistical estimators (for example, in terms of support recovery
(Hazimeh and Mazumder, 2020)).

Several recent studies attempt to further improve the efficiency of $\ell_{0}$
solvers. Along the line of dual methods, Liu et al. (2017); Yuan et al. (2020)
recently studied the duality of $k$-sparse problem, and they proved the strong
duality of $k$-sparse problem. With the derived dual form by Liu et al.
(2017); Yuan et al. (2020), a dual space based hard-thresholding method has
been proposed by the authors. In addition, a screening method has been
proposed by Atamtürk and Gomez (2020). Following coordinate descent (CD)
methods (Breheny and Huang, 2011; Mazumder et al., 2011; Friedman et al.,
2010; Nesterov, 2012) for linear regression problems, Hazimeh and Mazumder
(2020) proposed an efficient CD based method to scale up the solutions of
problem (1). Their method can be improved with the proposed switch techniques
that aim to escape from local solutions. Additionally, the combination of
$\ell_{0}$ and $\ell_{q}$ ($q=1$ or 2) are considered in many existing
literature. See Liu and Wu (2007); Soubies et al. (2017) for theoretical
analyses and details.

Following the studies in Pilanci et al. (2015); Liu et al. (2017); Yuan et al.
(2020), we investigate the dual form of the generalized sparse problem (1).
Under mild conditions, a strong duality theory has been established for
problem (1). A primal-dual algorithm is proposed to further improve the
efficiency and quality of solutions by leveraging the exploration in the dual
space along with coordinate screening and active incremental techniques
(Fercoq et al., 2015; Ndiaye et al., 2015, 2017; Atamtürk and Gomez, 2020;
Massias et al., 2018; Ren et al., 2020). In summary, our contributions on
theoretical side are three-folds. We first establish the duality theorem for
the generalized non-convex sparse learning problem (1). We second demonstrate
that the derived duality allows us to adopt the screening and coordinate
incremental strategies (Fercoq et al., 2015; Ndiaye et al., 2015, 2017;
Massias et al., 2018; Ren et al., 2020) in $\ell_{1}$ solvers to boost the
efficiency of the proposed $\ell_{0}$ algorithm. Thirdly, we provide
theoretical study of the proposed algorithms. Our theoretical analysis shows
that the generalized sparse problem (1) can be solved within polynomial
computation complexity. Experiments on both synthetic and real-world datasets
show the advantages of our method.

The rest of paper is organized as follows. In Section 2, we formulate the dual
form of the generalized sparse learning problem and also introduce the duality
theory. In Section 3, we propose the new primal-dual algorithm improved with
coordinate incremental techniques. Section 4 presents our algorithm analysis.
Experimental results are provided in Section 5. A discussion is given in
Section 6, and the concluding remark is presented in Section 7.

Notation. Symbol $\beta\in\mathbb{R}^{p}$ is used for the primal variable and
$\alpha\in\mathbb{R}^{n}$ is for the dual variable. We use $\|\beta\|$,
$\|\beta\|_{0}$ and $\|\beta\|_{1}$ to denote the $\ell_{2}$, $\ell_{0}$ and
$\ell_{1}$ norm of $\beta$, respectively. Functions $P(\beta)$ and $D(\alpha)$
represent the primal objective and the dual objective correspondingly. For
matrix $X$, $\sigma_{max}(X)$/$\sigma_{min}(X)$ denotes its largest/smallest
singular value. $\mathrm{supp}(\beta)$ is the support set of vector $\beta$,
i.e. $\mathrm{supp}(\beta)=\\{j|\beta_{j}\neq 0\\}$. $S^{c}$ represents the
complement of set $S$.

##  2 Duality Theory of Sparse Learning

This section extends the duality studies in Pilanci et al. (2015); Liu et al.
(2017); Yuan et al. (2020) to the generalized sparse learning problem (1). Let
$X=[x_{1},x_{2},...,x_{n}]^{\top}$ be the feature matrix,
$y=[y_{1},y_{2},...,y_{n}]^{\top}$ be the response vector and $n$ is the
number of samples. Let
$l^{*}_{i}(\alpha_{i})=\max_{u\in\mathcal{F}}\\{\alpha_{i}u-l_{i}(u)\\}$ be
the Fenchel conjugate (Fenchel, 1949) of convex loss function $l_{i}(u)$ and
$\mathcal{F}\subseteq\mathbb{R}$ be the feasible set of $\alpha_{i}$ regarding
$l^{*}_{i}()$. According to the expression
$l_{i}(u)=\max_{\alpha_{i}\in\mathcal{F}}\\{\alpha_{i}u-l_{i}^{*}(\alpha_{i})\\}$,
the primal problem can be reformulated into

$\displaystyle\min_{\beta}\sum_{i=1}^{n}\max_{\alpha_{i}\in\mathcal{F}}$
$\displaystyle\big{(}\alpha_{i}\beta^{\top}x_{i}-l_{i}^{*}(\alpha_{i})\big{)}+\lambda_{0}||\beta||_{0}+\lambda_{1}||\beta||_{1}+\lambda_{2}||\beta||^{2}.$
(3)

We use $L(\beta,\alpha)$ to represent the following objective

$\displaystyle~{}L(\beta,\alpha)=$
$\displaystyle\sum_{i=1}^{n}\big{(}\alpha_{i}\beta^{\top}x_{i}-l_{i}^{*}(\alpha_{i})\big{)}+\lambda_{0}||\beta||_{0}+\lambda_{1}||\beta||_{1}+\lambda_{2}||\beta||^{2}.$
(4)

Similar to the studies in Liu et al. (2017); Yuan et al. (2020), the RIP
(restricted strong condition number) bound conditions are not explicitly
required here. Without specifying $k$ in (2), our duality theory is close to
the standard duality paradigm, and thus generic primal-dual methods may be
used to further improve the solvers.

###  2.1 Strong Duality

We establish the duality theory that guarantees the original non-convex in (1)
can be solved in a dual space. Following Yuan et al. (2020); Liu et al.
(2017), we define the saddle point for the Lagrangian (4) of the generalized
sparse learning (1).

######  Definition 2.1.

(Saddle Point). A pair
$(\bar{\beta},\bar{\alpha})\in\mathbb{R}^{p}\times\mathcal{F}^{n}$ is said to
be a saddle point for $L$ (4) if the following holds

$\displaystyle L(\bar{\beta},\alpha)\leq L(\bar{\beta},\bar{\alpha})\leq
L(\beta,\bar{\alpha}).$ (5)

Different from the sparse saddle point in Yuan et al. (2020); Liu et al.
(2017) that requires $k$-sparse regarding the primal variable, the saddle
point defined here can be taken as a generalized saddle point. Given
$\alpha\in\mathcal{F}^{n}$, we further define
$\eta(\alpha):=-\frac{1}{2\lambda_{2}}\sum_{i=1}^{n}\alpha_{i}x_{i}=-\frac{1}{2\lambda_{2}}X^{\top}\alpha$,
$\eta_{0}:=\frac{2\sqrt{\lambda_{0}\lambda_{2}}+\lambda_{1}}{2\lambda_{2}}$,
and

$\displaystyle\beta_{j}(\alpha):=\mathfrak{B}(\eta_{j}(\alpha)):=\begin{cases}\mathrm{sign}\big{(}\eta_{j}(\alpha)\big{)}\big{(}|\eta_{j}(\alpha)|-\frac{\lambda_{1}}{2\lambda_{2}}\big{)}&\text{if}\
\ |\eta_{j}(\alpha)|>\eta_{0}\\\
\big{\\{}0,\mathrm{sign}\big{(}\eta_{j}(\alpha)\big{)}\big{(}|\eta_{j}(\alpha)|-\frac{\lambda_{1}}{2\lambda_{2}}\big{)}\big{\\}}&\text{if}\
\ |\eta_{j}(\alpha)|=\eta_{0}\\\ 0&\text{if}\ \
|\eta_{j}(\alpha)|<\eta_{0}\end{cases}.$ (6)

Moreover, we define

$\displaystyle\Psi(\eta_{j}(\alpha);\lambda_{0},\lambda_{1},\lambda_{2}):=\begin{cases}-\lambda_{2}(|\eta_{j}(\alpha)|-\frac{\lambda_{1}}{2\lambda_{2}})^{2}+\lambda_{0}&\text{if}\
\ |\eta_{j}(\alpha)|>\eta_{0}\\\
\big{\\{}0,-\lambda_{2}(|\eta_{j}(\alpha)|-\frac{\lambda_{1}}{2\lambda_{2}})^{2}+\lambda_{0}\big{\\}}\quad&\text{if}\
\ |\eta_{j}(\alpha)|=\eta_{0}\\\ 0&\text{if}\ \
|\eta_{j}(\alpha)|<\eta_{0}\end{cases}.$ (7)

The following theorem establishes the duality theory regarding the generalized
sparse problem (1).

######  Theorem 2.1.

Let $\bar{\beta}\in\mathbb{R}^{p}$ be a primal vector and
$\bar{\alpha}\in\mathcal{F}^{n}$ regarding L, then

  1. 1.

$(\bar{\alpha},\bar{\beta})$ is a saddle point of $L$ if and only if the
following conditions hold:

    1. (a)

$\bar{\beta}$ solves the primal problem;

    2. (b)

$\bar{\alpha}\in[\partial l_{1}(\bar{\beta}^{\top}x_{1}),\partial
l_{2}(\bar{\beta}^{\top}x_{2}),...,\partial
l_{n}(\bar{\beta}^{\top}x_{n})]^{\top}$;

    3. (c)

$\bar{\beta}_{j}=\mathfrak{B}(\eta_{j}(\bar{\alpha}))$.

  2. 2.

The mini-max relationship

$\displaystyle\max_{\alpha\in\mathcal{F}^{n}}\min_{\beta}L(\beta,\alpha)=\min_{\beta}\max_{\alpha\in\mathcal{F}^{n}}L(\beta,\alpha).$
(8)

holds if and only if there exists a saddle point $(\bar{\beta},\bar{\alpha})$
for L.

  3. 3.

The corresponding dual problem of (1) is written as

$\displaystyle\max_{\alpha\in\mathcal{F}^{n}}D(\alpha)=\max_{\alpha\in\mathcal{F}^{n}}-\sum_{i=1}^{n}l_{i}^{*}(\alpha_{i})+\sum_{j=1}^{p}\Psi(\eta_{j}(\alpha);\lambda_{0},\lambda_{1},\lambda_{2}),$
(9)

where $l^{*}$ is the conjugate function of $l$. The primal dual link is
written as $\beta_{j}(\alpha)=\mathfrak{B}(\eta_{j}(\alpha))$.

  4. 4.

(Strong duality) $\bar{\alpha}$ solves the dual problem in (22), i.e.,
$D(\bar{\alpha})\geq D(\alpha),\alpha\in\mathcal{F}^{n}$, and
$P(\bar{\beta})=D(\bar{\alpha})$ if and only if the pair
$(\bar{\beta},\bar{\alpha})$ satisfies the three conditions given by
(a)$\sim$(c).

Here $\eta_{0}$ is the threshold which controls the sparsity of the solution.
Larger $\lambda_{0}$ and $\lambda_{1}$ lead to sparser estimator. The mini-max
result in Theorem 2.1-1 gives the sufficient and necessary conditions to
guarantee the existence of a saddle point for the Lagrangian. Theorem 2.1-2
can be used to establish the duality theory, and it is the min-max side of the
problem, and it provides sufficient and necessary conditions under which one
can safely exchange a min-max for a max-min regarding $L$ (4).

######  Remark 2.1.

Applying Theorem 2.1, we have the following mini-max relationship

$\displaystyle\max_{\alpha\in\mathcal{F}^{n}}\min_{\beta}L(\beta,\alpha)=\min_{\beta}\max_{\alpha\in\mathcal{F}^{n}}L(\beta,\alpha).$
(10)

s holds if and only if there exists a primal vector
$\bar{\beta}\in\mathbb{R}^{p}$ and a dual vector
$\bar{\alpha}\in\mathcal{F}^{n}$ such that conditions (a) $\sim$ (c) in
Theorem2.1-1 are satisfied. Moreover, by calculations, it can be checked that
(10) holds automatically for $l(\cdot)$ being the square loss function.

We use $P(\beta)=F(\beta)$ to represent the primal objective, and $D(\alpha)$
for the dual objective given in (22). Theorem 2.1-3 indicates that the dual
objective function $D(\alpha)$ is concave and the following remark explicitly
gives the expression of its super-differential.

######  Remark 2.2.

The super-differential of the dual form (22) at $\alpha$ is given by $\nabla
D(\alpha)=X\beta(\alpha)-l^{*^{\prime}}(\alpha)=[\beta(\alpha)^{\top}x_{1}-{l^{*}_{1}}^{{}^{\prime}}(\alpha_{1}),...,\beta(\alpha)^{\top}x_{n}-{l^{*}_{n}}^{{}^{\prime}}(\alpha_{n})]^{\top}$.

The super-gradient can be alternatively derived through the partial derivative
of the Lagrangian L (4) regarding $\alpha$. The sparse strong duality theory
in Theorem 2.1-4 gives the sufficient and necessary conditions under which the
optimal values of the primal and dual problems coincide. We define the primal-
dual gap as

$\displaystyle\xi(\beta,\alpha)=P(\beta)-D(\alpha).$ (11)

According to Theorem 2.1-4, the primal-dual gap reaches zero at the primal-
dual pair $(\beta,\alpha)$ if and only if the conditions (a) $\sim$ (c) in
Theorem 2.1-1 hold. The duality theory developed in this section suggests a
natural way for finding the global minimum of the sparsity-constrained
minimization problem in (1) via primal-dual optimization methods. Let
$(l^{*^{\prime}})^{-1}$ be the inverse of $l^{*^{\prime}}$, we have the
following remark with $0\in\nabla D(\bar{\alpha})$ at $\bar{\alpha}$.

######  Remark 2.3.

If $(\bar{\beta},\bar{\alpha})$ satisfies the conditions in Theorem 2.1-1, we
have $\bar{\alpha}\in(l^{*^{\prime}})^{-1}(X\bar{\beta})\cap\mathcal{F}^{n}$.

Strong duality holds when both the primal and dual variables reach the optimal
values. Before attaining the optimal values, the duality gap value can be
bounded by the current dual variable estimations. The closer the current
estimation $\alpha$ and the optimal value $\bar{\alpha}$ are, the smaller
duality gap will be. Under special cases that the support of $\bar{\beta}$ is
recovered, then the objective function becomes a convex function since
$\|\beta^{t}\|_{0}$ remains a constant. Then strong duality holds naturally.
In practice, the condition (a) of Theorem 2.1-1 is hard to be satisfied
because of the non-convexity of the primal problem. However, as long as
$\beta$ reaches its optimal value, all the conditions in Theorem 2.1-1 can be
met because the dual problem is concave.

Different from Yuan et al. (2020); Liu et al. (2017), we study a generalized
sparse problem. The methodology developed here can be easily extended to plain
$\ell_{1}$ or $\ell_{0}$ problems (with the $\ell_{2}$ term), group sparse
structures or fused sparse structures, or even more complex and mixed sparse
structures that we cannot or do not need to specify the active feature number
$k$ value as in (2).

###  2.2 Properties of Generalized Sparse Learning

In this paper, we study the duality of generalized sparse learning problem.
Based on the strong duality of $\ell_{0}$ problem (1), screening methods
(Fercoq et al., 2015; Ndiaye et al., 2015, 2017) and coordinate increasing
techniques (Massias et al., 2018; Ren et al., 2020) can be implemented to gain
an improvement in algorithm efficiency. Following the Gap screening (Fercoq et
al., 2015; Ndiaye et al., 2015) for Lasso problem, we have the following
theorem regarding the duality gap.

######  Theorem 2.2.

Assume that the primal loss functions $\\{l_{i}(\cdot)\\}_{i=1}^{n}$ are
$1/\mu$-strongly smooth. The range of the dual variable is bounded via the
duality gap value, i.e.,
$\forall\alpha\in\mathcal{F}^{n},\beta\in\mathbb{R}^{p}$,
$\\{B(\alpha;r):||\alpha-\bar{\alpha}||_{2}\leq
r,r=\sqrt{\frac{2(P(\beta)-D(\alpha))}{\gamma}}\\}$. Here $\gamma$ is a
positive constant and $\gamma\geq\mu$.

Let $x_{\cdot i}$ be the $i$th column of X, according to the definition of
$\mathfrak{B}$ in (6), the activity of feature $j$ is determined by the
magnitude of $\eta_{j}(\bar{\alpha})$, i.e.,
$|\eta_{j}(\bar{\alpha})|=\frac{1}{2\lambda_{2}}|x^{\top}_{\cdot
j}\bar{\alpha}|$. With the ball region estimation for $\bar{\alpha}$ in
Theorem 2.2, we can estimate the activity of a feature with the value of
current $\alpha$. Let $r=\sqrt{\frac{2(P(\beta)-D(\alpha))}{\gamma}}$ be the
radius of the estimated ball range for $\bar{\alpha}$ using current $\alpha$
and $\beta$ solutions. Then $\big{|}|x^{\top}_{\cdot j}\alpha|-||x_{\cdot
j}||_{2}r\big{|}\leq|x^{\top}_{\cdot j}\bar{\alpha}|\leq|x^{\top}_{\cdot
j}\alpha|+||x_{\cdot j}||_{2}r$ and we get
$|\eta_{i}(\bar{\alpha})|\leq\frac{1}{2\lambda_{2}}\big{(}|x^{\top}_{\cdot
j}\alpha|+||x_{\cdot j}||_{2}r\big{)}<\eta_{0}\implies x_{\cdot j}\ \text{is
an inactive feature}.$ It implies

$\displaystyle\text{Upper Bound}:|x^{\top}_{\cdot j}\alpha|+||x_{\cdot
j}||_{2}r<2\lambda_{2}\eta_{0}=2\sqrt{\lambda_{0}\lambda_{2}}+\lambda_{1}\implies
j\notin\mathrm{supp}(\bar{\beta}),$ $\displaystyle\text{Lower
Bound}:\big{|}|x^{\top}_{\cdot j}\alpha|-||x_{\cdot
j}||_{2}r\big{|}>2\lambda_{2}\eta_{0}=2\sqrt{\lambda_{0}\lambda_{2}}+\lambda_{1}\implies
j\in\mathrm{supp}(\bar{\beta}).$

According to the derived dual objective (22) and equations (6)-(7), a
feature’s activity is determined by its product with the optimal dual variable
$\bar{\alpha}$, e.g., $\eta_{j}(\bar{\alpha})$ for feature $j$. The dual range
estimation $\bar{\alpha}$ ($B(\bar{\alpha};r)$) allows us to perform feature
screen in order to improve algorithm efficiency by following the approach for
Lasso (Fercoq et al., 2015; Ndiaye et al., 2015, 2017). As the support set $S$
is unknown, we just set $\gamma=\mu$ to ensure the safety of the feature
screening. Here safety means the screening operation does not remove any
feature belonging to $S=\mathrm{supp}(\bar{\beta})$. The framework proposed in
this paper lays a broader bridge between screening methods and the solutions
of $\ell_{0}$ regularized problems.

##  3 Algorithm

With the mini-max relationship in Theorem 2.1 regarding Lagrangian form
$L(\beta,\alpha)$ (4) , we first develop a primal-dual algorithm to update
both $\alpha$ and $\beta$. The dual objective $D(\alpha)$ is a non-smooth
function as the term $\Psi()$ regarding $\alpha$ is non-smooth due to the
truncation operation. We focus on the following simplified dual form

$\displaystyle\max_{\alpha\in\mathcal{F}^{n}}-\sum_{i=1}^{n}l_{i}^{*}(\alpha_{i})+\sum_{j=1}^{p}\Psi(\eta_{j}(\alpha);\lambda_{0},\lambda_{1},\lambda_{2}),\
\text{and}$ (12)
$\displaystyle\Psi(\eta_{j}(\alpha);\lambda_{0},\lambda_{1},\lambda_{2})=\begin{cases}-\lambda_{2}(|\eta_{j}(\alpha)|-\frac{\lambda_{1}}{2\lambda_{2}})^{2}+\lambda_{0}&\text{if}\
\ |\eta_{j}(\alpha)|\geq\eta_{0}\\\ 0&\text{if}\ \
|\eta_{j}(\alpha)|<\eta_{0}\end{cases}.$ (13)

The primal dual link is

$\displaystyle\beta_{j}(\alpha)=\mathfrak{B}(\eta_{j}(\alpha))=\begin{cases}\mathrm{sign}\big{(}\eta_{j}(\alpha)\big{)}\big{(}|\eta_{j}(\alpha)|-\frac{\lambda_{1}}{2\lambda_{2}}\big{)}&\text{if}\
\ |\eta_{j}(\alpha)|\geq\eta_{0}\\\ 0&\text{if}\ \
|\eta_{j}(\alpha)|<\eta_{0}\end{cases}.$ (14)

The super-gradient regarding the dual variable
$\nabla_{\alpha}D(\alpha)=X\beta(\alpha)-{l^{*}}^{{}^{\prime}}(\alpha)$ can be
taken as the partial derivative of $L(\beta,\alpha)$. We give the dual
problems of two objective functions in the supplements, and we will focus on
linear regression to present the proposed algorithms.

###  3.1 Primal-dual Updating for Linear Regression

We use linear regression as an example to illustrate the proposed primal-dual
inner solver of $\ell_{0}$ regularized problems. For least square problem, the
primal form is

$\displaystyle\min_{\beta\in\mathbb{R}^{p}}P(\beta)=f(\beta)+\lambda_{0}||\beta||_{0},\
\text{and}\
f(\beta)=\frac{1}{2}||y-X\beta||^{2}_{2}+\lambda_{1}||\beta||_{1}+\lambda_{2}||\beta||_{2}^{2}.$

For least square problem, $l_{i}(u;y_{i})=\frac{1}{2}(y_{i}-u)^{2}$, then
$l_{i}^{*}(\alpha_{i})=\frac{1}{2}\big{(}(y_{i}+\alpha_{i})^{2}-y_{i}^{2}\big{)}=\frac{1}{2}\alpha_{i}^{2}+y_{i}\alpha_{i}$.
Here $u=x_{i}^{\top}\beta$. Thus the dual problem is

$\displaystyle\max_{\alpha}D(\alpha)=$
$\displaystyle-\frac{1}{2}\alpha^{\top}\alpha-y^{\top}\alpha+\sum_{j=1}^{p}\Psi(-\frac{1}{2\lambda_{2}}\sum_{i=1}^{n}\bar{\alpha}_{i}x_{i};\lambda_{0},\lambda_{1},\lambda_{2}).$
(15)

The corresponding super-gradient can be easily computed, i.e.,
$g_{\alpha}=\nabla_{\alpha}D(\alpha)=X\beta-Y-\alpha.$ After the super
gradient ascent for the dual variables, we apply the primal-dual link function
to get the variable in the primal space.

The dual objective $D(\alpha)$ is non-smooth. The super-gradient
$g_{\alpha}=\nabla_{\alpha}D(\alpha)$ can be improved with a more accurate
primal variable estimation for the Lagrangian form (4) with the mini-max
relationship in Theorem 2.1. We use coordinate descent (CD) (Hazimeh and
Mazumder, 2020) to improve the estimation of primal variable as

$\displaystyle\beta_{j}=T(\mathbf{\beta};\lambda_{0},\lambda_{1},\lambda_{2})=\begin{cases}\mathrm{sign}(\tilde{\mathbf{\beta}}_{j})\frac{|\tilde{\mathbf{\beta}}_{j}|-\lambda_{1}}{1+2\lambda_{2}}\quad&\text{if}\
\frac{|\tilde{\mathbf{\beta}}_{j}|-\lambda_{1}}{1+2\lambda_{2}}\geq\sqrt{\frac{2\lambda_{0}}{1+2\lambda_{2}}}\\\
0\quad&\text{if}\
\frac{|\tilde{\mathbf{\beta}}_{j}|-\lambda_{1}}{1+2\lambda_{2}}<\sqrt{\frac{2\lambda_{0}}{1+2\lambda_{2}}}\end{cases}.$
(16)

Here $\tilde{\mathbf{\beta}}_{j}=(y-X\mathbf{\beta})^{\top}x_{\cdot
j}+\mathbf{\beta}_{j}x_{\cdot j}^{\top}x_{\cdot j}$, and $x_{\cdot j}$ is the
$j^{th}$ column of $X$, and it is also named the $j^{th}$ feature. The
operation $T()$ always decreases the primal objective, i.e. with
$\beta_{b}=T(\beta_{a})$, we always have $P(\beta_{b})\leq P(\beta_{a})$, and
hence a smaller duality gap.

Input: data $\\{X,y\\}$; $\lambda_{0}$, $\lambda_{1}$, $\lambda_{2}$; step
size $\omega_{t}$ at step $t$; initial $\mathbf{\alpha}^{0}$,
$\mathbf{\beta}^{0}$

Result: $\mathbf{\alpha}^{t}$, $\mathbf{\beta}^{t}$

1

2 $t\leftarrow 0$;

3 while  _$\overline{DGap}$ decreasing_ do

4 //Super-gradient

5
$g_{\mathbf{\alpha}}^{t}=[{\beta^{t-1}}^{\top}x_{1}-{l^{*}_{1}}^{{}^{\prime}}(\alpha^{t-1}_{1}),...,{\beta^{t-1}}^{\top}x_{n}-{l^{*}_{n}}^{{}^{\prime}}(\alpha^{t-1}_{n})]^{\top}$;

6 //Dual ascent with feasible projection

7
$\alpha^{t}=\mathcal{P}_{\mathcal{F}}(\alpha^{t-1}+\omega_{t}g_{\mathbf{\alpha}}^{t})$
;

8 $\eta(\alpha^{t})=-\frac{1}{2\lambda_{2}}\sum_{i=1}^{n}\alpha^{t}_{i}x_{i}$;

9 //Primal-dual relation

10
$\mathbf{\beta}^{t}=\mathfrak{B}(\eta_{j}(\alpha^{t});\lambda_{0},\lambda_{1},\lambda_{2})$;

11 //Primal coordinate descent

12 $\mathbf{\beta}^{t}=T(\beta^{t};\lambda_{0},\lambda_{1},\lambda_{2})$;

13 Compute duality gap $\overline{DGap}$ with $\beta_{t}$ and $\alpha_{t}$ ;

14 $t\leftarrow t+1$;

Algorithm 1 Inner solver with primal-dual updating

The proposed primal-dual updating procedure is given by Algorithm 1. The
primal coordinate descent $T()$ improves the solution from primal-dual
relation $\mathfrak{B}$. $\omega_{t}$ is the step size at $t$, and should be
decreasing with $t$. We use Algorithm 1 as the backbone solver in our primal-
dual algorithm, and $\overline{DGap}$ is the sub-problem’s duality gap
achieved by the inner solver.

Moreover, according to Remark 2.3, the optimal dual
$\bar{\alpha}\in(l^{*^{\prime}})^{-1}(X\bar{\beta})\cap\mathcal{F}^{n}$. For a
solver in primal space, we can use this equation to find a point in dual space
and then compute the duality gap to evaluate current solution. For linear
regression, with $\beta$ we have $\alpha=X\beta-y$, then we can compute the
duality gap via (11).

###  3.2 Improve Efficiency with Active Incremental Strategy

For sparse models, most of the features are redundant and they incur extra
computation costs. The derived dual problem structure and the duality property
provide an approach to implement feature screening (Fercoq et al., 2015;
Ndiaye et al., 2015, 2017) and feature active incremental strategy (Massias et
al., 2018; Ren et al., 2020). According to the analysis in Section 2.2, the
activity of a feature $x_{\cdot j}$ depends on the value of
$\eta_{j}(\bar{\alpha})$, i.e.,
$|\eta_{j}(\bar{\alpha})|=\frac{1}{2\lambda_{2}}|x^{\top}_{\cdot
j}\bar{\alpha}|$. We use the current estimation range of $\bar{\alpha}$, i.e.,
$B(\bar{\alpha};r^{s})$ to approximate the value of $\eta_{j}(\bar{\alpha})$.
Here $s$ is the step number in the outer loop of the algorithm, and
$\beta^{s}$ and $\alpha^{s}$ are the primal dual solutions at step $s$. The
ball radius $r^{s}=\sqrt{\frac{2(P(\beta^{s})-D(\alpha^{s}))}{\gamma}}$ for
$\bar{\alpha}$ depends the duality gap at step $s$ according to Theorem 2.2.

Input: $\alpha^{s}$, $r^{s}$, $\mathcal{A}^{s}$, $\mathcal{R}^{s}$

Result: $\mathcal{R}^{s+1}$, $\mathcal{A}^{s+1}$

1

2 Set $h=\lceil c\log(p)\rceil$, $\mathcal{H}\leftarrow$ Select $h$ features
according to the descending order of $x_{\cdot
j}^{\top}\alpha^{s},j\in\mathcal{R}^{s}$;

3$\mathcal{A}^{s+1}\leftarrow\mathcal{A}^{s}\cup\mathcal{H}$,
$\mathcal{R}^{s+1}\leftarrow\mathcal{R}^{s}\setminus\mathcal{H}$ ;

Algorithm 2 Algorithm for Feature Adding

Input: $X$, $y$; $\lambda_{0}$, $\lambda_{1}$, $\lambda_{2}$; $\xi$

Result: $\beta^{s}$

1

2 Choose a small set of features $\mathcal{A}_{0}$ from $X$ in the descending
order of $|X^{\top}l^{\prime}(\mathbf{0})|$, and $\mathcal{R}_{0}$ represents
the rest features;

3 $s\leftarrow 0$; $\beta^{0}\leftarrow 0$; $\alpha^{0}\leftarrow 0$;
$DoAdd=True$;

4 while  _True_ do

5 //Sub-problem solver

6 $\alpha^{s},\tilde{\beta}^{s}\leftarrow$ Solve the sub-problem with feature
set $\mathcal{A}^{s}$ via Algorithm 1;

7 $\beta^{s}\leftarrow$ put $\tilde{\beta}^{s}$ in a size $p$ vector and set
entries not in $\mathcal{A}^{s}$ zero;

8 //Dual range estimation

9 Compute the duality gap $DGap$ and the ball region $B(\alpha^{s};r^{s})$
with $r^{s}=\sqrt{\frac{2(P(\beta^{s})-D(\alpha^{s}))}{\gamma}}$ ;

10 if  _$DGap <\xi$ _ then

11 Stop; //Algorithm exits

12 Feature Screening with $B(\alpha^{s};r^{s})$ and (17);

13 if  _ $DoAdd$ _ then

14 //Feature Adding

15 if  _ $\max_{j\in\mathcal{R}^{s}}|x^{\top}_{\cdot j}\alpha^{s}|+\|x_{\cdot
j}\|_{2}r^{s}<2\lambda_{2}\eta_{0}$ _ then

16 $DoAdd=False$; Continue;

17 Feature Adding with $\alpha^{s}$ and Algorithm 2, and update
$\mathcal{A}^{s}$ and $\mathcal{R}^{s}$;

18 $s\leftarrow s+1$;

Algorithm 3 Primal-dual algorithm

The proposed primal-dual algorithm for $\ell_{0}$ is given by Algorithm 3.
Algorithm 3 starts with a small active set $\mathcal{A}$, and then increase
the active set’s size after solving each sub-problem. We use $\mathcal{R}$ to
represent the set of features not used by the sub-problem solver. The feature
adding algorithm is given by Algorithm 2. Moreover, we can derive a gap-
screening algorithm (Fercoq et al., 2015; Ndiaye et al., 2015) by using the
upper bound of $|\eta_{j}(\bar{\alpha})|$’s approximation given in Section
2.2. Based on the derivation in Section 2.2, we use the following safe
principle for feature screening.

$\displaystyle\textbf{Feature Screening:}\ |x^{\top}_{\cdot
j}\alpha^{s}|+||x_{\cdot j}||_{2}r^{s}<2\lambda_{2}\eta_{0}\implies
j\notin\mathrm{supp}(\bar{\beta}),\ \text{remove}.$ (17)

Here $2\lambda_{2}\eta_{0}=2\sqrt{\lambda_{0}\lambda_{2}}+\lambda_{1}$. The
screening rule is safe because it is derived based on concavity of the dual
problem. Base on (17), we derive a stopping condition for feature adding. If
all features in $\mathcal{R}$ satisfy (17), we stop Feature Adding. This
allows us to avoid redundant computation resulted from some inactive features.

In Algorithm 3, the initialization values of $\beta^{0}$ and $\alpha^{0}$ are
set to zero. We use $DGap$ to represent the duality gap of the original
problem attained by the primal-dual algorithm. Please note that we use $s$ to
denote the step in Algorithm 2 and Algorithm 3, in order to differ from steps
($t$s) in Algorithm 1. Empirically, feature screening does show power to
improve algorithm efficiency. But the feature active incremental strategy can
significantly avoid redundant computation to achieve the target duality gap.

##  4 Algorithm Analysis

We discuss the algorithm convergence in this section. Firstly, we present
theoretical results on the convergence and support recovery of the inner
solver for sub-problems.

###  4.1 Algorithm Analysis for Inner Solver

With the conditions in Theorem 2.1-1, we can quantify the duality gap and
develop the algorithm analysis. Different from the primal updating (Hazimeh
and Mazumder, 2020) or dual updating (Liu et al., 2017) algorithms, Algorithm
1 has both primal and dual updating steps. Let p denotes the size of the input
feature set $X$ of the sub-problem in Algorithm 1, $m_{1}=\max_{j:1\leq j\leq
p}|y^{\top}x_{\cdot j}|$, $m_{2}=\max_{j:1\leq j\leq p}||X^{\top}x_{\cdot
j}||$,
$m_{3}=\max_{i,\alpha_{i}^{t}\in\mathcal{F}}|l_{i}^{*^{\prime}}(\alpha_{i}^{t})|$,
$\varrho=\sqrt{n}||\alpha^{t}||_{\infty}-\lambda_{1}$, and $\omega_{t}$ is the
decreasing step size. We have the following theorem regarding the convergence
of Algorithm 1.

######  Theorem 4.1.

Assume that $l_{i}$ is $1/\mu$-smooth, $||x_{i}||\leq\vartheta\ \forall 1\leq
i\leq n$, and $||x_{\cdot j}||=1\ \forall 1\leq j\leq p$. By choosing
$w_{t}=\frac{1}{t\gamma}$, then the sequence generated by Algorithm 1
satisfies the following estimation error inequality:

$\displaystyle||\alpha^{t}-\bar{\alpha}||^{2}\leq
c_{1}\bigg{(}\frac{1}{t}+\frac{\ln t}{t}\bigg{)}.$

Here $c_{1}=\frac{c_{0}^{2}}{\mu^{2}}$,
$c_{0}=\frac{\sqrt{np}\vartheta}{2\lambda_{2}(1+2\lambda_{2})}(2\lambda_{2}m_{1}+\varrho+\sqrt{n}m_{2}\varrho-2\lambda_{1}\lambda_{2})+\sqrt{n}m_{3}$,
$gamma$ is same as in Theorem 2.2.

We can prove the convergence of the primal variable using the results in
Theorem 4.1. Let $S=\mathrm{supp}(\bar{\beta})$,
$N=\\{j|\eta_{j}(\bar{\alpha})=\eta_{0(j)}\\}$, and
$\bar{\delta}=2\lambda_{2}\min\bigg{\\{}\min_{j:j\in
S}\frac{|\eta_{j}(\bar{\alpha})|-\eta_{0}}{||x_{\cdot j}||},\min_{j:j\in
S^{c}\setminus N}\frac{\eta_{0}-|\eta_{j}(\bar{\alpha})|}{||x_{\cdot
j}||}\bigg{\\}}$. The following theorem gives the complexity for support
recovery and duality gap convergence.

######  Theorem 4.2.

Assume that $l_{i}$ is $1/\mu$-smooth, $||x_{i}||\leq\vartheta\ \forall 1\leq
i\leq n$, and $||x_{\cdot j}||=1\ \forall 1\leq j\leq p$. Let
$t_{1}=\frac{3c_{1}}{\bar{\delta}^{2}}\log\frac{3c_{1}}{\bar{\delta}^{2}}$,
with $t>t_{1}$, we have
$\mathrm{supp}(\beta(\alpha))=\mathrm{supp}(\bar{\beta})$ and
$||\beta(\alpha)-\bar{\beta}||\leq\frac{\sigma_{max}(X_{S})}{2\lambda_{2}}||\alpha-\bar{\alpha}||$.
Moreover, let
$t_{2}=\frac{3c_{1}c_{2}}{\epsilon}\log\frac{3c_{1}c_{2}}{\epsilon}$,
$c_{2}=c_{0}\bigg{(}1+\frac{\sigma_{max}(X_{S})}{2\mu\lambda_{2}}\bigg{)}$,
for any $\epsilon>0$ with $t>\mathrm{max}\\{t_{1},t_{2}\\}$, we have
$P(\beta^{t})-D(\alpha^{t})\leq\epsilon$.

Additional analysis on algorithms can be found in the appendix. Due to the
large magnitudes of $c_{0}$ and $c_{2}$ and small value of $\bar{\delta}$, it
could be time consuming for the solver to achieve very small duality gaps.

###  4.2 Outer Loop Analysis

The outer loop in Algorithm 3 involves both feature screening and feature
adding operations relying on the dual variable estimation for the original
problem, i.e., $||\alpha^{s}-\bar{\alpha}||_{2}\leq
r^{s},r^{s}=\sqrt{\frac{2(P(\beta^{s})-D(\alpha^{s}))}{\gamma}}$ defined in
Theorem 2.2. Here $\gamma=\mu+\frac{\sigma_{min}(X_{S})}{2\lambda_{2}}$. As
discussed in Section 2.2, we set $\gamma=\mu$ to ensure the safety of feature
screening.

######  Remark 4.1.

The screening operation (17) is safe, and it does not remove any features in
$\mathcal{A}^{s}\cap\mathrm{supp}(\bar{\beta})$ at step $s$. With additional
features added by the feature adding operation (Algorithm 2), the primal
objective $P(\beta^{s})$ always decreases after the solution of the sub-
problem regarding feature set $\mathcal{A}^{s}$ using the inner solver
(Algorithm 1).

The screening operation usually keeps the primal objective value
$P(\beta^{s})$ intact. With Remark 4.1, $\mathcal{A}^{s}$ and $P(\beta^{s})$
converge after some steps, and also Algorithm 3 converges with $DGap$ smaller
than a given threshold $\xi$. In fact, the active incremental strategy could
significantly reduce redundant operations introduced by inactive features
especially when the problem is with high sparse level (Massias et al., 2018;
Ren et al., 2020). The solution sparse level (the size of
$S=\mathrm{supp}(\bar{\beta})$) impacts the algorithm complexity.

##  5 Experiments

Experiments focus on linear regression. Our proposed algorithm can be extended
to other forms of loss function. Via experimental studies, we show the
effectiveness of our method by comparing with dual iterative hard thresholding
(Liu et al., 2017) and coordinate descent with spacer steps (Hazimeh and
Mazumder, 2020) algorithms, which is the state-of-art for $\ell_{0}$
regularization problem. The experimental environment is CPU: Intel Xeon
<EMAIL_ADDRESS>GHz; OS: Windows Server 2012 R2.

###  5.1 Simulation Study

In this study, we simulate the datasets under the linear regression setting,
i.e., $\mathbf{y}=X\beta+\epsilon$. The data matrix is generated according to
a multi-variate Gaussian $X_{n\times p}\sim\textrm{MVN}(0,\Sigma)$, and
$\Sigma=(\sigma_{i,j})$. Exponential correlation (Hazimeh and Mazumder, 2020)
is utilized to control feature relationship, i.e., $\sigma_{ij}=\rho^{|i-j|}$
with $\rho=0.4$. The noise $\epsilon$ is Gaussian white noise with
$SNR=\\{2,5,20\\}$ . For the true parameter $\beta$, $3\%$ entries ($0.03p$)
are randomly set to the values in $[-1.0,1.0]$, and the rest ($0.97p$) are set
to zero. We generate the datasets with $p=3\,000$ and $n$ varying in
$\\{200,300,400,500,600\\}$. Each setting is replicated for 50 times.

Figure 1: Running time, duality gap, nonzero number, PSSR, and estimation
error for four algorithms on simulated data with $SNR\in\\{20,5,2\\}$. x-axis
is the number of training samples, from 200 to 600.

We compare our primal-dual algorithm against dual iterative hard thresholding
(Liu et al., 2017) and coordinate descent with spacer steps (Hazimeh and
Mazumder, 2020) algorithms. We use ‘Dual-IHT’, and ‘CDSS’ to represent the two
algorithms, respectively. We use ‘PrimDual’ to represent the proposed primal
dual algorithm, and ‘DualAst’ to represent dual ascent method which is
Algorithm 1 without the primal updating steps. All algorithms are implemented
in Matlab and run in the same environment. ‘Dual-IHT’ (Liu et al., 2017) is a
primal-dual method with dual ascending using hard threshold to keep $k$
largest value of $|\beta|$ in the primal space. CDSS (Hazimeh and Mazumder,
2020) is a coordinate descent method operates in the primal space enhanced
with PSI (Partial Swap Inescapable) as the stopping criteria. We use the
duality gap ($DGap$) threshold $\epsilon=1.0E-6$ as the stopping condition for
Dual-IHT, DualAst and the proposed PrimDual algorithm. The algorithms may
require extremely long time to reach a small duality gap threshold. We also
use the duality change, i.e., $\zeta=|DGap^{t-2}-DGap^{t}|$ as a stopping
condition for the three algorithms, and we set $\zeta=1.0E-6$ in the
experiments. Moreover, we use the same learning rate $\omega=0.000\,5$ for the
three algorithms.

Two indices are adopted for evaluating the performance. The first one is the
percentage of successful support recovery (PSSR). The second one is parameter
estimation error $||\beta-\bar{\beta}||/||\bar{\beta}||$. Here $\bar{\beta}$
is the ground truth used in simulation. Figure 1 gives the performance of
these four algorithms on datasets with different SNR values. To achieve
meaningful comparison, we choose $\lambda_{0},\lambda_{1}$ and $\lambda_{2}$
to recover support number close to the ground truth value. For the dataset
with $SNR=20$, we use $\lambda_{0}=0.03,\lambda_{1}=0.02,\lambda_{2}=1.0$, and
we set $\lambda_{0}=0.1,\lambda_{1}=0.2,\lambda_{2}=1.0$ for dataset with
$SNR=2,5$. From the plots, we can see that the proposed primal-dual algorithm
can achieve similar PSSR and estimation error values (expect for Dual-IHT
since its sparsity is pre-determined), but use much less time. It shows that
the proposed primal-dual algorithm and incremental strategy significantly
reduce the redundant operations resulted from inactive features.

###  5.2  Experimental on Real-world Datasets

In this section, we present additional results to compare the four algorithms
on two datasets, News20 and E2006. In these experiments, we use learning rate
$\omega=0.000\,5$ for Dual-IHT, DualAst, and our PrimDual method. The stopping
conditions for the three algorithms are $\epsilon=1.0E-6$ and $\zeta=1.0E-6$.
We set $c=4.0$ for Algorithm 2 in our experiments. It is difficult to fairly
compare D-IHT with the other three methods as it uses different objective with
hard constraints, and we have to specify the hyper-parameter $k$ in advance.
The values of $k$ are set heuristically in the experiments.

####  5.2.1 News20 Dataset

After pre-processing, the commonly used News20 dataset contains 20 classes,
$15\,935$ samples, and $62\,061$ features in the training set. The 20 labels
in News20 dataset are transformed to response values ranging $[-10,10]$ in the
experiments. We randomly pick up $p=2\,000$ features and to form five datasets
with sample number ranging in $\\{200,300,400,500,600\\}$. We use learning
rate $\omega=0.0005$ for Dual-IHT, DualAst, and our PrimDual method. The
stopping conditions are $\epsilon=1.0E-6$ and $\zeta=1.0E-6$. The hyper-
parameters are set with $\lambda_{0}=0.1,\lambda_{1}=0.15,\lambda_{1}=1.0$.
The left column of Figure 2 shows the results of different methods on News20
dataset with $p=2\,000$. Each setting is replicated for 20 times. We can see
that under approximately the same primal objective and duality gap values, the
proposed primal-dual method uses less computation time compared against other
methods when $N$ becomes larger. Though DualAst consumes similar computation
cost as our method, it cannot achieve small duality gap values on all cases.
We notice that CDSS takes longest time in this case, and it could be due to
that the PSI stopping condition is hard to satisfy on some real-world
datasets.

Figure 2: The running time, primal objective $P(\hat{\beta})$, duality gap,
and nonzero number for different methods on News20 dataset.

The right column of Figure 2 gives other results on News20 datasets with a
larger feature number ($p=4\,000$). Hyper-parameters are set as
$\lambda_{0}=0.10,\lambda_{1}=0.15,\lambda_{2}=1.00$. We set $k=300$ for D-IHT
method. Each case is replicated for 10 times and the average result is
reported. From the plots, we can see that with longer running time, CDSS
achieves smallest duality gap values. However, the proposed method uses much
less computation cost to achieve similar solutions.

####  5.2.2 E2006

E2006 regression dataset has 150 360 features and 16 087 samples for
training222https://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/regression.html.
We also randomly select $p=4\,000$ features and $n=\\{200,300,400,500,600\\}$
samples in our experiments. Hyper-parameters are set to
$\lambda_{0}=0.10,\lambda_{1}=0.10,\lambda_{2}=1.00$, and $k=300$ for D-IHT.
Each case is replicated for 10 times and the average result is reported.
Figure 3 gives results on E2006 datasets. Again, the proposed primal-dual
algorithm is much more efficient in solving the problems to achieve similar
duality gaps compared against other methods. The results validate the duality
theory, problem properties, and the proposed incremental strategy.

Figure 3: The running time, primal objective $P(\hat{\beta})$, duality gap,
and nonzero number for different methods on E2006 datasets with $p=4\,000$.

##  6 Discussion

In this section, we first give several additional remarks on the method, and
then we compare our approach with some existing methods.

###  6.1 Additional Remarks

  * Technical Contribution: This paper investigates the dual form and strong duality of the generalized sparse learning problem (1) by following Pilanci et al. (2015); Liu et al. (2017); Yuan et al. (2020). The generalized form could overcome over-fitting issues of $\ell_{0}$ regularized problems when the data SNRs are low (Gamarnik and Zadik, 2017; Mazumder et al., 2022). A primal-dual framework has been developed to further scale up the solutions of $\ell_{0}$ regularized problems based on the derived dual form. Moreover, the proposed framework considers active coordinate incremental and screening strategies (Fercoq et al., 2015; Ndiaye et al., 2015, 2017; Atamtürk and Gomez, 2020; Massias et al., 2018; Ren et al., 2020) by leveraging the duality structure properties of problem (1). The quality of solutions can be evaluated by the duality gap (11) with the current dual solution $\alpha$ calculated through Remark 2.3.

  * Saddle Point: Different from the sparse saddle point defined in Liu et al. (2017); Yuan et al. (2020) that requires $k$-sparse regarding the primal variable, the saddle point in this paper can be taken as a standard saddle point. Without the specified $k$, our duality theory is more close to the standard duality paradigm, and hence some generic primal-dual methods can be employed to further improve the solver. The methodology developed here can be easily extended to plain $\ell_{1}$ or $\ell_{0}$ problems (with the $\ell_{2}$ term), group sparse structures, fused sparse structures, or even more complex and mixed sparse structures that we cannot or do not need to specify the $k$ values.

  * Strong Duality: Strong duality holds when both the primal and dual variables reach the optimal values. A closer distance between the current estimation $\alpha$ and the optimal value $\bar{\alpha}$ gives a smaller duality gap. Once the support of $\bar{\beta}$ is recovered, the objective function becomes a convex function since $\|\beta^{s}\|_{0}$ remains a constant. According to our theoretical analysis, the saddle point of problem (1) could be attained within polynomial computation complexity with a decreasing step size. When $l(\cdot)$ takes some special convex function, strong duality holds. However, for general $l(\cdot)$, it is hard to verify the strong duality.

###  6.2 Comparisons

We provide more details on the differences between our work and related works
in this section to highlight our technical contributions.

The proposed method is significantly different from Dual-IHT method (Liu et
al., 2017; Yuan et al., 2020). Firstly, the proposed primal-dual algorithm
focuses on a different problem (1) compared with the Dual-IHT’s objective (2).
Apart from using soft regularization rather than hard constraint, (1) also
includes the $\ell_{1}$ penalty that could be helpful in the cases with low
SNR values. Secondly, our primal-dual method perform updating in both primal
and dual spaces to approach the saddle point, and it can potentially attain
solutions with smaller duality gaps. Finally, most important of all, our
objective (1) and the derived dual form (22) allow us to employ screening and
coordinate incremental strategies (Fercoq et al., 2015; Ndiaye et al., 2015;
Massias et al., 2018; Ren et al., 2020; Atamtürk and Gomez, 2020) to boost the
efficiency of the algorithm.

There are several obvious differences between our primal-dual method and the
coordinate descent with spacer steps (CDSS) method (Hazimeh and Mazumder,
2020). Different from our primal-dual method, CDSS utilizes coordinate descent
in the primal space for parameter updating. In CDSS (Hazimeh and Mazumder,
2020), they also rely on partial swap inescapable (PSI-$\kappa$) to improve
the solution with $\kappa=1$. PSI with $\kappa\geq 2$ will introduce much more
extra computation that is usually not affordable. Our primal-dual method
employs coordinate incremental strategy to save computation cost, and
experimental results indicate that the proposed method can achieve similar
solution quality as CDSS but with less computation time.

In Atamtürk and Gomez (2020), the authors propose a screening method for
$\ell_{0}$ regularized problems. However, their objective does not include the
$\ell_{1}$ norm. Our propose method focuses on a more generalized problem that
is potentially more powerful on datasets with low SNRs. Moreover, besides the
safe screening rule proposed in Section 2.2, the proposed coordinate
incremental strategy introduced in Section 3.2 is empirically effective on
different datasets. The screening methods (Fercoq et al., 2015; Ndiaye et al.,
2015, 2017) and coordinate incremental strategies (Ren et al., 2020; Massias
et al., 2018) used for $\ell_{1}$ regularized problems can be taken as special
cases of the proposed method with $\lambda_{0}=0$.

##  7 Conclusion

In this paper we studied the dual forms of a broad family of $\ell_{0}$
regularized problems. Based on the derived dual form, a primal-dual algorithm
accelerated with active coordinate selection has been developed. Our
theoretical result show that the reformed best subset selection problem can be
solved within polynomial complexity. The developed theory and the proposed
framework can be integrated with many screening strategies. Experimental
results show that the proposed primal-dual method can reduce redundant
operations introduced by the inactive features and hence save computation
costs. The proposed framework sheds light on primal-dual algorithms that can
potentially further scale up the solutions of $\ell_{0}$ regularized non-
convex sparse problems.

## References

  * Atamtürk and Gomez [2020] Alper Atamtürk and Andres Gomez.  Safe screening rules for l0-regression from perspective relaxations.  In _Proceedings of the 37th International Conference on Machine Learning (ICML)_ , pages 421–430, Virtual Event, 2020.
  * Bertsimas and Van Parys [2020] Dimitris Bertsimas and Bart Van Parys.  Sparse high-dimensional regression: Exact scalable algorithms and phase transitions.  _The Annals of Statistics_ , 48(1):300–323, 2020\.
  * Bertsimas et al. [2016] Dimitris Bertsimas, Angela King, and Rahul Mazumder.  Best subset selection via a modern optimization lens.  _The Annals of Statistics_ , 44(2):813–852, 2016\.
  * Bian and Chen [2020] Wei Bian and Xiaojun Chen.  A smoothing proximal gradient algorithm for nonsmooth convex regression with cardinality penalty.  _SIAM J. Numer. Anal._ , 58(1):858–883, 2020\.
  * Blumensath and Davies [2009] Thomas Blumensath and Mike E Davies.  Iterative hard thresholding for compressed sensing.  _Applied and Computational Harmonic Analysis_ , 27(3):265–274, 2009.
  * Breheny and Huang [2011] Patrick Breheny and Jian Huang.  Coordinate descent algorithms for nonconvex penalized regression, with applications to biological feature selection.  _The Annals of Applied Statistics_ , 5(1):232, 2011.
  * Dedieu et al. [2021] Antoine Dedieu, Hussein Hazimeh, and Rahul Mazumder.  Learning sparse classifiers: Continuous and mixed integer optimization perspectives.  _J. Mach. Learn. Res._ , 22:135:1–135:47, 2021.
  * Dong et al. [2015] Hongbo Dong, Kun Chen, and Jeff Linderoth.  Regularization vs. relaxation: A conic optimization perspective of statistical variable selection.  _arXiv preprint arXiv:1510.06083_ , 2015.
  * Fenchel [1949] Werner Fenchel.  On conjugate convex functions.  _Canadian Journal of Mathematics_ , 1(1):73–77, 1949.
  * Fercoq et al. [2015] Olivier Fercoq, Alexandre Gramfort, and Joseph Salmon.  Mind the duality gap: safer rules for the lasso.  In _Proceedings of the 32nd International Conference on Machine Learning (ICML)_ , pages 333–342, Lille, France, 2015.
  * Foucart [2011] Simon Foucart.  Hard thresholding pursuit: An algorithm for compressive sensing.  _SIAM J. Numer. Anal._ , 49(6):2543–2563, 2011\.
  * Friedman et al. [2010] Jerome Friedman, Trevor Hastie, and Rob Tibshirani.  Regularization paths for generalized linear models via coordinate descent.  _Journal of Statistical Software_ , 33(1):1, 2010\.
  * Furnival and Wilson [1974] George M. Furnival and Robert W. Wilson.  Regressions by leaps and bounds.  _Technometrics_ , 16(4):499–511, 1974.
  * Gamarnik and Zadik [2017] David Gamarnik and Ilias Zadik.  High dimensional regression with binary coefficients. estimating squared error and a phase transtition.  In _Proceedings of the 30th Conference on Learning Theory (COLT)_ , pages 948–953, Amsterdam, The Netherlands, 2017.
  * Hastie et al. [2017] Trevor Hastie, Robert Tibshirani, and Ryan J Tibshirani.  Extended comparisons of best subset selection, forward stepwise selection, and the lasso.  _arXiv:1707.08692_ , 2017.
  * Hazimeh and Mazumder [2020] Hussein Hazimeh and Rahul Mazumder.  Fast best subset selection: Coordinate descent and local combinatorial optimization algorithms.  _Oper. Res._ , 68(5):1517–1537, 2020.
  * Hazimeh et al. [2020] Hussein Hazimeh, Rahul Mazumder, and Ali Saab.  Sparse regression at scale: Branch-and-bound rooted in first-order optimization.  _arXiv preprint arXiv:2004.06152_ , 2020.
  * Jain et al. [2014] Prateek Jain, Ambuj Tewari, and Purushottam Kar.  On iterative hard thresholding methods for high-dimensional m-estimation.  In _Advances in Neural Information Processing Systems (NIPS)_ , pages 685–693, Montreal, Canada, 2014.
  * Jain et al. [2016] Prateek Jain, Nikhil Rao, and Inderjit S. Dhillon.  Structured sparse regression via greedy hard thresholding.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , pages 1516–1524, Barcelona, Spain, 2016.
  * Liu et al. [2017] Bo Liu, Xiao-Tong Yuan, Lezi Wang, Qingshan Liu, and Dimitris N. Metaxas.  Dual iterative hard thresholding: From non-convex sparse minimization to non-smooth concave maximization.  In _Proceedings of the 34th International Conference on Machine Learning (ICML)_ , pages 2179–2187, Sydney, Australia, 2017.
  * Liu and Wu [2007] Yufeng Liu and Yichao Wu.  Variable selection via a combination of the $l_{0}$ and $l_{1}$ penalties.  _Journal of Computational and Graphical Statistics_ , 16(4):782–798, 2007.
  * Loh and Wainwright [2017] Po-Ling Loh and Martin J Wainwright.  Support recovery without incoherence: A case for nonconvex regularization.  _The Annals of Statistics_ , 45(6):2455–2482, 2017.
  * Massias et al. [2018] Mathurin Massias, Joseph Salmon, and Alexandre Gramfort.  Celer: a fast solver for the lasso with dual extrapolation.  In _Proceedings of the 35th International Conference on Machine Learning (ICML)_ , pages 3321–3330, Stockholmsmässan, Stockholm, Sweden, 2018.
  * Mazumder and Radchenko [2017] Rahul Mazumder and Peter Radchenko.  The discrete dantzig selector: Estimating sparse linear models via mixed integer linear optimization.  _IEEE Trans. Inf. Theory_ , 63(5):3053–3075, 2017.
  * Mazumder et al. [2011] Rahul Mazumder, Jerome H Friedman, and Trevor Hastie.  SparseNet: Coordinate descent with nonconvex penalties.  _Journal of the American Statistical Association_ , 106(495):1125–1138, 2011.
  * Mazumder et al. [2022] Rahul Mazumder, Peter Radchenko, and Antoine Dedieu.  Subset selection with shrinkage: Sparse linear modeling when the SNR is low.  _Oper. Res._ , 2022.
  * Natarajan [1995] Balas K. Natarajan.  Sparse approximate solutions to linear systems.  _SIAM J. Comput._ , 24(2):227–234, 1995.
  * Ndiaye et al. [2015] Eugène Ndiaye, Olivier Fercoq, Alexandre Gramfort, and Joseph Salmon.  GAP safe screening rules for sparse multi-task and multi-class models.  In _Advances in Neural Information Processing Systems (NIPS)_ , pages 811–819, Montreal, Canada, 2015.
  * Ndiaye et al. [2017] Eugène Ndiaye, Olivier Fercoq, Alexandre Gramfort, and Joseph Salmon.  Gap safe screening rules for sparsity enforcing penalties.  _J. Mach. Learn. Res._ , 18:128:1–128:33, 2017.
  * Nesterov [2012] Yurii E. Nesterov.  Efficiency of coordinate descent methods on huge-scale optimization problems.  _SIAM J. Optim._ , 22(2):341–362, 2012.
  * Parikh and Boyd [2014] Neal Parikh and Stephen P. Boyd.  Proximal algorithms.  _Found. Trends Optim._ , 1(3):127–239, 2014.
  * Pilanci et al. [2015] Mert Pilanci, Martin J. Wainwright, and Laurent El Ghaoui.  Sparse learning via boolean relaxations.  _Math. Program._ , 151(1):63–87, 2015.
  * Ren et al. [2020] Shaogang Ren, Weijie Zhao, and Ping Li.  Thunder: a fast coordinate selection solver for sparse learning.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , virtual, 2020.
  * Shen and Li [2017] Jie Shen and Ping Li.  A tight bound of hard thresholding.  _J. Mach. Learn. Res._ , 18:208:1–208:42, 2017.
  * Soubies et al. [2017] Emmanuel Soubies, Laure Blanc-Féraud, and Gilles Aubert.  A unified view of exact continuous penalties for $\ell_{2}-\ell_{0}$ minimization.  _SIAM J. Optim._ , 27(3):2034–2060, 2017.
  * Soussen et al. [2015] Charles Soussen, Jérôme Idier, Junbo Duan, and David Brie.  Homotopy based algorithms for $\ell_{0}$-regularized least-squares.  _IEEE Trans. Signal Process._ , 63(13):3301–3316, 2015.
  * Yang and Yu [2019] Yingzhen Yang and Jiahui Yu.  Fast proximal gradient descent for A class of non-convex and non-smooth sparse learning problems.  In _Proceedings of the Thirty-Fifth Conference on Uncertainty in Artificial Intelligence (UAI)_ , pages 1253–1262, Tel Aviv, Israel, 2019.
  * Yuan and Li [2020] Xiao-Tong Yuan and Ping Li.  Nearly non-expansive bounds for mahalanobis hard thresholding.  In _Proceedings of Conference on Learning Theory (COLT)_ , pages 3787–3813, Virtual Event [Graz, Austria], 2020.
  * Yuan et al. [2014] Xiao-Tong Yuan, Ping Li, and Tong Zhang.  Gradient hard thresholding pursuit for sparsity-constrained optimization.  In _Proceedings of the 31th International Conference on Machine Learning (ICML)_ , pages 127–135, Beijing, China, 2014.
  * Yuan et al. [2020] Xiao-Tong Yuan, Bo Liu, Lezi Wang, Qingshan Liu, and Dimitris N Metaxas.  Dual iterative hard thresholding.  _J. Mach. Learn. Res._ , 21:152–1, 2020.
  * Zhu et al. [2020] Junxian Zhu, Canhong Wen, Jin Zhu, Heping Zhang, and Xueqin Wang.  A polynomial algorithm for best-subset selection problem.  _Proceedings of the National Academy of Sciences_ , 117(52):33117–33123, 2020.

In this Appendix, we provide theoretical proofs supporting the results in the
main context. Section A deliveries the proofs of Theorem 2.1 and Theorem 2.2
regarding the duality; Section B gives the algorithm analysis; Section C
presents the dual forms of two loss functions.

##  Appendix A Proofs of Theorem 2.1 and Theorem 2.2

In this section, we provide the proofs to extend the duality theory [Pilanci
et al., 2015, Liu et al., 2017, Yuan et al., 2020] to the generalized sparse
learning problem (1). The derivation of duality presented here is
significantly different from the duality of hard thresholding [Liu et al.,
2017, Yuan et al., 2020] because that the generalized problem (1) uses soft-
regularization terms rather than hard constrains, and it also includes the
combination of three regularization norms, i.e., $\ell_{0}$-, $\ell_{1}$-, and
$\ell_{2}$-norms.

######  Lemma A.1.

For a given $\alpha\in\mathcal{F}^{n}$, let
$\beta(\alpha)=\operatorname*{argmin}_{\beta}L(\beta,\alpha)$. We have

$\displaystyle\min_{\beta}L(\beta,\alpha)=-\sum_{i=1}^{n}l_{i}^{*}(\alpha_{i})+\sum_{j=1}^{p}\Psi(\eta_{j}(\alpha);\lambda_{0},\lambda_{1},\lambda_{2}),$

where

$\displaystyle\Psi(\eta_{j}(\alpha);\lambda_{0},\lambda_{1},\lambda_{2})=\begin{cases}-\lambda_{2}(|\eta_{j}(\alpha)|-\frac{\lambda_{1}}{2\lambda_{2}})^{2}+\lambda_{0}&\text{if}\
\ |\eta_{j}(\alpha)|>\eta_{0}\\\
\big{\\{}0,-\lambda_{2}(|\eta_{j}(\alpha)|-\frac{\lambda_{1}}{2\lambda_{2}})^{2}+\lambda_{0}\big{\\}}\quad&\text{if}\
\ |\eta_{j}(\alpha)|=\eta_{0}\\\ 0&\text{if}\ \
|\eta_{j}(\alpha)|<\eta_{0}\end{cases}.$

To be specific, $l^{*}$ is the conjugate function of $l$,
$\eta(\alpha):=-\frac{1}{2\lambda_{2}}\sum_{i=1}^{n}\alpha_{i}x_{i}$ and
$\eta_{0}:=\frac{2\sqrt{\lambda_{0}\lambda_{2}}+\lambda_{1}}{2\lambda_{2}}$.
The link function between $\alpha$ and $\beta$ is

$\displaystyle~{}~{}~{}~{}\beta_{j}(\alpha)=\mathfrak{B}(\eta_{j}(\alpha))=\begin{cases}\mathrm{sign}\big{(}\eta_{j}(\alpha)\big{)}\big{(}|\eta_{j}(\alpha)|-\frac{\lambda_{1}}{2\lambda_{2}}\big{)}&\text{if}\
\ |\eta_{j}(\alpha)|>\eta_{0}\\\
\big{\\{}0,\mathrm{sign}\big{(}\eta_{j}(\alpha)\big{)}\big{(}|\eta_{j}(\alpha)|-\frac{\lambda_{1}}{2\lambda_{2}}\big{)}\big{\\}}&\text{if}\
\ |\eta_{j}(\alpha)|=\eta_{0}\\\ 0&\text{if}\ \
|\eta_{j}(\alpha)|<\eta_{0}\end{cases}.$

###### Proof.

With $l^{*}$ as the conjugate of $l$, the primal problem can be rewritten as

$\displaystyle\min_{\beta\in\mathbb{R}^{p}}L(\beta,\alpha)$ $\displaystyle=$
$\displaystyle\min_{\beta\in\mathbb{R}^{p}}\sum_{i=1}^{n}\big{(}\alpha_{i}\beta^{\top}x_{i}-l_{i}^{*}(\alpha_{i})\big{)}+\lambda_{0}||\beta||_{0}+\lambda_{1}||\beta||_{1}+\lambda_{2}||\beta||^{2}_{2}$
$\displaystyle=$
$\displaystyle\min_{\beta\in\mathbb{R}^{p}}\sum_{i=1}^{n}\big{(}\alpha_{i}\beta^{\top}x_{i}-l_{i}^{*}(\alpha_{i})\big{)}+\lambda_{0}||\beta||_{0}+\lambda_{1}||\beta||_{1}+\lambda_{2}||\beta||^{2}_{2}$
$\displaystyle=$
$\displaystyle-\sum_{i=1}^{n}l_{i}^{*}(\alpha_{i})+\min_{\beta\in\mathbb{R}^{p}}\sum_{i=1}^{n}\big{(}\alpha_{i}\beta^{\top}x_{i}\big{)}+\lambda_{0}||\beta||_{0}+\lambda_{1}||\beta||_{1}+\lambda_{2}||\beta||^{2}_{2}$
$\displaystyle=$
$\displaystyle-\sum_{i=1}^{n}l_{i}^{*}(\alpha_{i})+\min_{\beta\in\mathbb{R}^{p}}\sum_{j=1}^{p}\big{(}\sum_{i=1}^{n}\alpha_{i}x_{ij}\big{)}\beta_{j}+\lambda_{0}||\beta||_{0}+\lambda_{1}||\beta||_{1}+\lambda_{2}||\beta||^{2}_{2}$
$\displaystyle=$
$\displaystyle-\sum_{i=1}^{n}l_{i}^{*}(\alpha_{i})+\sum_{j=1}^{p}\Phi(\sum_{i=1}^{n}\alpha_{i}x_{ij};\lambda_{0},\lambda_{1},\lambda_{2}).$

Let $\tau_{j}=\alpha^{\top}x_{\cdot j}$, then

$\displaystyle\Phi(\tau_{j};\lambda_{0},\lambda_{1},\lambda_{2})=$
$\displaystyle\min_{\beta_{j}}\tau_{j}\beta_{j}+\lambda_{0}\mathbf{I}\\{\beta_{j}\neq
0\\}+\lambda_{1}|\beta_{j}|+\lambda_{2}\beta_{j}^{2}$ $\displaystyle=$
$\displaystyle\min_{\beta_{j}}\lambda_{2}\bigg{(}\beta_{j}^{2}+\frac{\tau_{j}}{\lambda_{2}}\beta_{j}\bigg{)}^{2}+\lambda_{0}\mathbf{I}\\{\beta_{j}\neq
0\\}+\lambda_{1}|\beta_{j}|$ $\displaystyle=$
$\displaystyle\min_{\beta_{j}}\lambda_{2}\bigg{(}\beta_{j}+\frac{\tau_{j}}{2\lambda_{2}}\bigg{)}^{2}+\lambda_{0}\mathbf{I}\\{\beta_{j}\neq
0\\}+\lambda_{1}|\beta_{j}|-\frac{\tau_{j}^{2}}{4\lambda_{2}}.$

Let

$\displaystyle
h(u)=\lambda_{2}\bigg{(}u+\frac{\tau_{j}}{2\lambda_{2}}\bigg{)}^{2}+\lambda_{0}\mathbf{I}\\{u\neq
0\\}+\lambda_{1}|u|-\frac{\tau_{j}^{2}}{4\lambda_{2}}.$

If $u\neq 0$, we have the soft-thresholding

$\displaystyle\hat{u}=-\frac{1}{2\lambda_{2}}\mathrm{sign}(\tau_{j})\max(|\tau_{j}|-\lambda_{1},0).$

With $|\tau_{j}|>\lambda_{1}$,

$\displaystyle h(\hat{u})$
$\displaystyle=\lambda_{2}\bigg{(}-\frac{\tau_{j}-\lambda_{1}\mathrm{sign}(\tau_{j})}{2\lambda_{2}}+\frac{\tau_{j}}{2\lambda_{2}}\bigg{)}^{2}+\lambda_{0}+\frac{\lambda_{1}}{2\lambda_{2}}(|\tau_{j}|-\lambda_{1})-\frac{\tau_{j}^{2}}{4\lambda^{2}_{2}}$
$\displaystyle=\lambda_{0}+\frac{\lambda_{1}|\tau_{j}|}{2\lambda_{2}}-\frac{\lambda_{1}^{2}}{4\lambda_{2}}-\frac{\tau_{j}^{2}}{4\lambda_{2}}$
$\displaystyle=\lambda_{0}-\frac{1}{4\lambda_{2}}\big{(}\lambda_{1}^{2}+\tau_{j}^{2}-2\lambda_{1}|\tau_{j}|\big{)}$
$\displaystyle=\lambda_{0}-\frac{1}{4\lambda_{2}}\big{(}|\tau_{j}|-\lambda_{1}\big{)}^{2}.$

With $h(\hat{u})<h(0)=0$, we get

$\displaystyle\lambda_{0}<\frac{1}{4\lambda_{2}}\big{(}|\tau_{j}|-\lambda_{1}\big{)}^{2}\implies|\tau_{j}|>2\sqrt{\lambda_{0}\lambda_{2}}+\lambda_{1}.$

Hence, $\hat{u}$ is the minimizer when
$|\tau_{j}|>2\sqrt{\lambda_{0}\lambda_{2}}+\lambda_{1}$. We have

$\displaystyle\Phi(\tau_{j};\lambda_{0},\lambda_{1},\lambda_{2})=-\frac{1}{4\lambda_{2}}\big{(}|\tau_{j}|-\lambda_{1}\big{)}^{2}+\lambda_{0}.$

When $2\sqrt{\lambda_{0}\lambda_{2}}+\lambda_{1}=|\tau_{j}|$, both $\hat{u}$
and 0 are minimizer. Then

$\displaystyle\Phi(\tau_{j};\lambda_{0},\lambda_{1},\lambda_{2})=\\{0,-\frac{1}{4\lambda_{2}}(|\tau_{j}|-\lambda_{1})^{2}+\lambda_{0}\\}.$

If $2\sqrt{\lambda_{0}\lambda_{2}}+\lambda_{1}>|\tau_{j}|$, 0 is the
minimizer, and $\Psi(\tau_{j},\lambda_{0},\lambda_{1},\lambda_{2})=0$.

The optimal primal can be written as

$\displaystyle\beta^{*}_{j}$
$\displaystyle=\mathfrak{A}(\tau_{j};\lambda_{0},\lambda_{1},\lambda_{2})=\begin{cases}\mathrm{sign}(\tau_{j})\frac{\lambda_{1}-|\tau_{j}|}{2\lambda_{2}}&\text{if}\
\ |\tau_{j}|>2\sqrt{\lambda_{0}\lambda_{2}}+\lambda_{1}\\\
\\{0,\mathrm{sign}(\tau_{j})\frac{\lambda_{1}-|\tau_{j}|}{2\lambda_{2}}\\}\quad&\text{if}\
\ |\tau_{j}|=2\sqrt{\lambda_{0}\lambda_{2}}+\lambda_{1}\\\ 0&\text{if}\ \
|\tau_{j}|<2\sqrt{\lambda_{0}\lambda_{2}}+\lambda_{1}\end{cases}$
$\displaystyle=\mathfrak{B}(\eta_{j}(\alpha)).$

Here
$\eta_{0}=\frac{2\sqrt{\lambda_{0}\lambda_{2}}+\lambda_{1}}{2\lambda_{2}}$.
With the optimal $\beta(\alpha)$, $L(\beta,\alpha)$can be written as

$\displaystyle\Phi(\tau_{j};\lambda_{0},\lambda_{1},\lambda_{2})$
$\displaystyle=\begin{cases}-\frac{1}{4\lambda_{2}}(|\tau_{j}|-\lambda_{1})^{2}+\lambda_{0}&\text{if}\
\ |\tau_{j}|>2\sqrt{\lambda_{0}\lambda_{2}}+\lambda_{1}\\\
\\{0,-\frac{1}{4\lambda_{2}}(|\tau_{j}|-\lambda_{1})^{2}+\lambda_{0}\\}\quad&\text{if}\
\ |\tau_{j}|=2\sqrt{\lambda_{0}\lambda_{2}}+\lambda_{1}\\\ 0&\text{if}\ \
|\tau_{j}|<2\sqrt{\lambda_{0}\lambda_{2}}+\lambda_{1}\end{cases}.$

Alternatively,

$\displaystyle\Psi(\eta_{j}(\alpha);\lambda_{0},\lambda_{1},\lambda_{2})$
$\displaystyle=\Phi(\tau_{j};\lambda_{0},\lambda_{1},\lambda_{2})$
$\displaystyle=\begin{cases}-\lambda_{2}(|\eta_{j}(\alpha)|-\frac{\lambda_{1}}{2\lambda_{2}})^{2}+\lambda_{0}&\text{if}\
\ |\eta_{j}(\alpha)|>\eta_{0}\\\
\big{\\{}0,-\lambda_{2}(|\eta_{j}(\alpha)|-\frac{\lambda_{1}}{2\lambda_{2}})^{2}+\lambda_{0}\big{\\}}\quad&\text{if}\
\ |\eta_{j}(\alpha)|=\eta_{0}\\\ 0&\text{if}\ \
|\eta_{j}(\alpha)|<\eta_{0}\end{cases}.$

Here $\eta(\alpha)=-\frac{1}{2\lambda_{2}}\sum_{i=1}^{n}\alpha_{i}x_{i}$, and
$\eta_{0}=\frac{2\sqrt{\lambda_{0}\lambda_{2}}+\lambda_{1}}{2\lambda_{2}}$. It
concludes the lemma. ∎

######  Lemma A.2.

(Saddle Point). Let $\bar{\beta}$ be a primal vector and $\bar{\alpha}$ a dual
vector. Then $(\bar{\alpha},\bar{\beta})$ is a saddle point of $L$ if and only
if the following conditions hold:

a) $\bar{\beta}$ solves the primal problem;

b) $\bar{\alpha}\in[\partial l_{1}(\bar{\beta}^{\top}x_{1}),\partial
l_{2}(\bar{\beta}^{\top}x_{2}),...,\partial
l_{n}(\bar{\beta}^{\top}x_{n})]^{\top}$;

c) $\bar{\beta}_{j}=\mathfrak{B}(\eta_{j}(\bar{\alpha}))$.

###### Proof.

$\Leftarrow$: If the pair $(\bar{\beta},\bar{\alpha})$ is a saddle point for
$L$, then from the definition of conjugate convexity and inequality in the
definition of saddle point we have

$\displaystyle
P(\bar{\beta})=\max_{\alpha\in\mathcal{F}^{n}}L(\bar{\beta},\alpha)\leq
L(\bar{\beta},\bar{\alpha})\leq\min_{\beta\in\mathbb{R}^{p}}L(\beta,\bar{\alpha}).$

On the other hand, we know that for any $\beta\in\mathbb{R}^{p}$ and
$\alpha\in\mathcal{F}^{n}$

$\displaystyle
L(\beta,\alpha)\leq\max_{\alpha^{\prime}\in\mathcal{F}^{n}}L(\beta,\alpha^{\prime})=P(\beta).$

By combining the preceding two inequalities we obtain

$\displaystyle
P(\bar{\beta})\leq\min_{\beta\in\mathbb{R}^{p}}L(\beta,\bar{\alpha})\leq\min_{\beta\in\mathbb{R}^{p}}P(\beta)\leq
P(\bar{\beta}).$

Therefore $P(\bar{\beta})=\min_{\beta\in\mathbb{R}^{p}}P(\beta)$, i.e.,
$\bar{\beta}$ solves the primal problem , which proves the necessary condition
(a). Moreover, the above arguments lead to

$\displaystyle
P(\bar{\beta})=\max_{\alpha\in\mathcal{F}^{n}}L(\bar{\beta},\alpha)=L(\bar{\beta},\bar{\alpha}).$

Then from the maximizing argument property of convex conjugate we have
$\bar{\alpha}_{i}\in\partial l_{i}(\bar{\beta}^{\top}x_{i})$, and it concludes
condition b). Note that

$\displaystyle
L(\beta,\bar{\alpha})=\lambda_{2}\bigg{|}\bigg{|}\beta+\frac{1}{2\lambda_{2}}\sum_{i=1}^{n}\bar{\alpha}_{i}x_{i}\bigg{|}\bigg{|}^{2}-\sum_{i=1}^{n}l_{i}^{*}(\bar{\alpha}_{i})+\lambda_{1}||\beta||_{1}+\lambda_{0}||\beta||_{0}+C.$
(18)

Let $\bar{F}=\mathrm{supp}(\bar{\beta})$. Since the above analysis implies
$L(\bar{\beta},\bar{\alpha})=\min_{\beta}L(\beta,\bar{\alpha})$, with
$\eta(\bar{\alpha})=-\frac{1}{2\lambda_{2}}\sum_{i=1}^{n}\bar{\alpha}_{i}x_{i}$,
it must hold that (more details refer to the proof of Lemma A.1)
$\bar{\beta}_{j}=\mathfrak{B}(\eta_{j}(\bar{\alpha}))$. This validates the
condition (c).

$\Rightarrow$: Inversely, let us assume that $\bar{\beta}$ is a solution to
the primal problem (condition a)), and $\bar{\alpha}_{i}\in\partial
l_{i}(\bar{\beta}^{\top}x_{i})$ (condition b)). Again from the maximizing
argument property of convex conjugate we know that
$l_{i}(\bar{\beta}^{\top}x_{i})=\bar{\alpha}\bar{\beta}^{\top}x_{i}-l_{i}^{*}(\bar{\alpha}_{i})$.
This leads to

$\displaystyle L(\bar{\beta},\alpha)\leq
P(\bar{\beta})=\max_{\alpha\in\mathcal{F}^{n}}L(\bar{\beta},\alpha)=\
L(\bar{\beta},\bar{\alpha}).$ (19)

The sufficient condition (c) guarantees that based on the expression of (18),
for any $\beta$, we have

$\displaystyle L(\bar{\beta},\bar{\alpha})\leq L(\beta,\bar{\alpha}).$ (20)

By combining the inequalities (19) and (20) we get that for any $\beta$ and
$\alpha$

$\displaystyle L(\bar{\beta},\alpha)\leq L(\bar{\beta},\bar{\alpha})\leq
L(\beta,\bar{\alpha}).$

This shows that $(\bar{\beta},\bar{\alpha})$ is a saddle point of the
Lagrangian $L$. ∎

Theorem 2.1 Let $\bar{\beta}\in\mathbb{R}^{p}$ be a primal vector and
$\bar{\alpha}\in\mathcal{F}^{n}$ regarding L, then

  1. 1.

$(\bar{\alpha},\bar{\beta})$ is a saddle point of $L$ if and only if the
following conditions hold:

    1. (a)

$\bar{\beta}$ solves the primal problem;

    2. (b)

$\bar{\alpha}\in[\partial l_{1}(\bar{\beta}^{\top}x_{1}),\partial
l_{2}(\bar{\beta}^{\top}x_{2}),...,\partial
l_{n}(\bar{\beta}^{\top}x_{n})]^{\top}$;

    3. (c)

$\bar{\beta}_{j}=\mathfrak{B}(\eta_{j}(\bar{\alpha}))$.

  2. 2.

The mini-max relationship

$\displaystyle\max_{\alpha\in\mathcal{F}^{n}}\min_{\beta}L(\beta,\alpha)=\min_{\beta}\max_{\alpha\in\mathcal{F}^{n}}L(\beta,\alpha).$
(21)

holds if and only if there exists a saddle point $(\bar{\beta},\bar{\alpha})$
for L.

  3. 3.

The corresponding dual problem of (1) is written as

$\displaystyle\max_{\alpha\in\mathcal{F}^{n}}D(\alpha)=\max_{\alpha\in\mathcal{F}^{n}}-\sum_{i=1}^{n}l_{i}^{*}(\alpha_{i})+\sum_{j=1}^{p}\Psi(\eta_{j}(\alpha);\lambda_{0},\lambda_{1},\lambda_{2}),$
(22)

where $l^{*}$ is the conjugate function of $l$. The primal dual link is
written as $\beta_{j}(\alpha)=\mathfrak{B}(\eta_{j}(\alpha))$.

  4. 4.

(Strong duality) $\bar{\alpha}$ solves the dual problem in (22), i.e.,
$D(\bar{\alpha})\geq D(\alpha),\alpha\in\mathcal{F}^{n}$, and
$P(\bar{\beta})=D(\bar{\alpha})$ if and only if the pair
$(\bar{\beta},\bar{\alpha})$ satisfies the three conditions given by
(a)$\sim$(c).

###### Proof.

According to Lemma A.2, statement-1 can be proved. We focus on statements 2-4.
The following Part (I), Part (II), and Part (III) are proofs of statement-2,
statement-3, and statement-4, respectively.

Part (I): the mimi-max relationship in statement-2.

$\Rightarrow$: Let $(\bar{\beta},\bar{\alpha})$ be a saddle point for $L$. On
one hand, note that the following holds for any $\beta^{\prime}$ and
$\alpha^{\prime}$,

$\displaystyle\min_{\beta}L(\beta,\alpha^{\prime})\leq
L(\beta^{\prime},\alpha^{\prime})\leq\max_{\alpha\in\mathcal{F}^{n}}L(\beta^{\prime},\alpha)$

which implies

$\displaystyle\max_{\alpha\in\mathcal{F}^{n}}\min_{\beta}L(\beta,\alpha)\leq\min_{\beta}\max_{\alpha\in\mathcal{F}^{n}}L(\beta,\alpha).$
(23)

On the other hand, since $(\bar{\beta},\bar{\alpha})$ is a saddle point for
$L$, the following is true:

$\displaystyle~{}\max_{\alpha\in\mathcal{F}^{n}}\min_{\beta}L(\beta,\alpha)\leq\max_{\alpha\in\mathcal{F}^{n}}L(\bar{\beta},\alpha)$
$\displaystyle\leq\max_{\alpha\in\mathcal{F}^{n}}L(\bar{\beta},\bar{\alpha})\leq\min_{\beta}L(\beta,\bar{\alpha})\leq\max_{\alpha\in\mathcal{F}^{n}}\min_{\beta}L(\beta,\alpha).$
(24)

By combining (23) and (24) we prove the equality (21).

$\Leftarrow$: Assume that the equality in (21) holds. Let us define
$\bar{\beta}$ and $\bar{\alpha}$ such that

$\max_{\alpha\in\mathcal{F}^{n}}L(\bar{\beta},\alpha)=\min_{\beta}\max_{\alpha\in\mathcal{F}^{n}}L(\beta,\alpha)$

and

$\min_{\beta}L(\beta,\bar{\alpha})=\max_{\alpha\in\mathcal{F}^{n}}\min_{\beta}L(\beta,\alpha).$

Then we can see that for any $\alpha$,
$L(\bar{\beta},\bar{\alpha})\geq\min_{\beta}L(\beta,\bar{\alpha})=\max_{\alpha^{\prime}\in\mathcal{F}^{n}}L(\bar{\beta},\alpha^{\prime})$,
where the “=” is due to (21). In the meantime, for any $\beta$

$\displaystyle
L(\bar{\beta},\bar{\alpha})\leq\max_{\alpha\in\mathcal{F}^{n}}L(\bar{\beta},\alpha)=\min_{\beta}L(\beta^{\prime},\bar{\alpha})\leq
L(\beta,\bar{\alpha}).$

This shows that $(\bar{\beta},\bar{\alpha})$ is a saddle point for L.

Part (II): the dual form in statement-3.

According to Lemma A.1, for any $\alpha$, the $\beta$ that minimizes
$L(\beta,\alpha)$ is

$\displaystyle\beta_{j}(\alpha)=\mathfrak{B}(\eta_{j}(\alpha))=\begin{cases}\mathrm{sign}\big{(}\eta_{j}(\alpha)\big{)}\big{(}|\eta_{j}(\alpha)|-\frac{\lambda_{1}}{2\lambda_{2}}\big{)}&\text{if}\
\ |\eta_{j}(\alpha)|>\eta_{0}\\\
\big{\\{}0,\mathrm{sign}\big{(}\eta_{j}(\alpha)\big{)}\big{(}|\eta_{j}(\alpha)|-\frac{\lambda_{1}}{2\lambda_{2}}\big{)}\big{\\}}\quad&\text{if}\
\ |\eta_{j}(\alpha)|=\eta_{0}\\\ 0&\text{if}\ \
|\eta_{j}(\alpha)|<\eta_{0}\end{cases}.$ (25)

Then we have

$\displaystyle
D(\alpha)=-\sum_{i=1}^{n}l_{i}^{*}(\alpha_{i})+\sum_{j=1}^{p}\Psi(\eta_{j};\lambda_{0},\lambda_{1},\lambda_{2}),$
(26)

where

$\displaystyle\Psi(\eta_{j}(\alpha);\lambda_{0},\lambda_{1},\lambda_{2})$
$\displaystyle=\begin{cases}-\lambda_{2}(|\eta_{j}(\alpha)|-\frac{\lambda_{1}}{2\lambda_{2}})^{2}+\lambda_{0}&\text{if}\
\ |\eta_{j}(\alpha)|>\eta_{0}\\\
\big{\\{}0,-\lambda_{2}(|\eta_{j}(\alpha)|-\frac{\lambda_{1}}{2\lambda_{2}})^{2}+\lambda_{0}\big{\\}}\quad&\text{if}\
\ |\eta_{j}(\alpha)|=\eta_{0}\\\ 0&\text{if}\ \
|\eta_{j}(\alpha)|<\eta_{0}\end{cases}$
$\displaystyle=-\lambda_{2}||\beta(\alpha)||^{2}_{2}+\lambda_{0}||\beta(\alpha)||_{0}.$
(27)

Here
$\eta_{0}=\frac{2\sqrt{\lambda_{0}\lambda_{2}}+\lambda_{1}}{2\lambda_{2}}$.
Assume we have two arbitrary dual variables
$\alpha_{1},\alpha_{2}\in\mathcal{F}^{n}$ and any
$g(\alpha_{2})\in[\beta(\alpha_{2})^{\top}x_{1}-{l_{1}^{*}}^{{}^{\prime}}(\alpha_{2(1)}),...,\beta(\alpha_{2})^{\top}x_{n}-{l_{n}^{*}}^{{}^{\prime}}(\alpha_{2(n)})]$.
Here $\alpha_{2(n)}$ is the $n$th entry of $\alpha_{2}$. $L(\beta,\alpha)$ is
concave in terms of $\alpha$ given any fixed $\beta$. According to the
definition of $D(\alpha)$, we have

$\displaystyle D(\alpha_{1})=L(\beta(\alpha_{1}),\alpha_{1})\leq
L(\beta(\alpha_{2}),\alpha_{1})\leq L(\beta(\alpha_{2}),\alpha_{2})+\langle
g(\alpha_{2}),\alpha_{1}-\alpha_{2}\rangle.$

Hence $D(\alpha)$ is concave and the super gradient is as given.

Part (III): Strong duality. $\Rightarrow$: Given the conditions a)-c), we can
see that the pair ($\bar{\beta}$, $\bar{\alpha}$) forms a saddle point of $L$.
Thus based on the definitions of saddle point and dual function $D(\alpha)$,
we can show that

$\displaystyle D(\bar{\alpha})=\min_{\beta}L(\beta,\bar{\alpha})\geq
L(\bar{\beta},\bar{\alpha})\geq L(\bar{\beta},\alpha)\geq D(\alpha).$

This implies that $\bar{\alpha}$ solves the dual problem. Furthermore, Theorem
2.1-2 guarantees the following

$\displaystyle
D(\bar{\alpha})=\max_{\alpha\in\mathcal{F}^{n}}\min_{\beta}L(\beta,\alpha)=\min_{\beta}\max_{\alpha\in\mathcal{F}^{n}}L(\beta,\alpha)=P(\bar{\beta}).$

This indicates that the primal and dual optimal values are equal to each
other.

$\Leftarrow$: Assume that $\bar{\alpha}$ solves the dual problem in and
$D(\bar{\alpha})=P(\bar{\beta})$. Since $D(\bar{\alpha})\leq P(\beta)$ holds
for any $\beta$, $\bar{\beta}$ must be the sparse minimizer of $P(\beta)$. It
follows that

$\displaystyle\max_{\alpha\in\mathcal{F}^{n}}\min_{\beta}L(\beta,\alpha)=D(\bar{\alpha})=P(\bar{\beta})=\min_{\beta}\max_{\alpha\in\mathcal{F}^{n}}L(\beta,\alpha).$

From the $\Leftarrow$ argument the proof of Theorem 2.1-2 and we get that
conditions a)-c)in Theorem 2.1-1 should be satisfied for
$(\bar{\beta},\bar{\alpha})$. This completes the proof. ∎

Compared to the dual problem developed in Liu et al. [2017], Yuan et al.
[2020] regarding hard thresholding, the soft thresholding in (25) corresponds
to the combination of $\ell_{1}$ and $\ell_{0}$ penalties, and it is helpful
on datasets with low SNR values [Gamarnik and Zadik, 2017, Mazumder et al.,
2022].

Theorem 2.2 Assume that the primal loss functions
$\\{l_{i}(\cdot)\\}_{i=1}^{n}$ are $1/\mu$-strongly smooth. The range of the
dual variable is bounded via the duality gap value, i.e.,
$\forall\alpha\in\mathcal{F}^{n},\beta\in\mathbb{R}^{p}$,
$\\{B(\alpha;r):||\alpha-\bar{\alpha}||_{2}\leq
r,r=\sqrt{\frac{2(P(\beta)-D(\alpha))}{\gamma}}\\}$. Here $\gamma$ is a
positive constant and $\gamma\geq\mu$.

###### Proof.

With the assumption $l_{i}$ being $1/\mu$-smooth, its conjugate function
$l^{*}_{i}$ is $\mu$-strongly convex. With the dual problem

$\displaystyle D(\alpha)$
$\displaystyle=-\sum_{i=1}^{n}l_{i}^{*}(\alpha_{i})+\sum_{j=1}^{p}\Psi(\eta_{j}(\alpha);\lambda_{0},\lambda_{1},\lambda_{2}).$

Here

$\displaystyle\Psi(\eta_{j}(\alpha);\lambda_{0},\lambda_{1},\lambda_{2})=\begin{cases}-\lambda_{2}(|\eta_{j}(\alpha)|-\frac{\lambda_{1}}{2\lambda_{2}})^{2}+\lambda_{0}&\text{if}\
\ |\eta_{j}(\alpha)|\geq\eta_{0}\\\ 0&\text{if}\ \
|\eta_{j}(\alpha)|<\eta_{0}\end{cases}.$

With
$\eta_{0}=\frac{2\sqrt{\lambda_{0}\lambda_{2}}+\lambda_{1}}{2\lambda_{2}}$ and
$\eta_{j}(\alpha)=-\frac{1}{2\lambda_{2}}x^{\top}_{\cdot j}\alpha$, we get
$\Psi(\eta_{j}(\alpha);\lambda_{0},\lambda_{1},\lambda_{2})$ is concave for
$j,1\leq j\leq p$ regarding $\alpha$, hence $D(\alpha)$ is concave. For given
hyper-parameters $\lambda_{0},\lambda_{1},\lambda_{2}$, let
$\Psi(\alpha)=\sum_{j=1}^{p}\Psi(\eta_{j}(\alpha);\lambda_{0},\lambda_{1},\lambda_{2})$,
and
$\Psi_{j}(\alpha)=\Psi(\eta_{j}(\alpha);\lambda_{0},\lambda_{1},\lambda_{2})$.
For $1\leq j\leq p$,

$\displaystyle\Psi_{j}^{\prime}(\alpha)=$
$\displaystyle\begin{cases}\mathrm{sign}(\eta_{j}(\alpha))(|\eta_{j}(\alpha)|-\frac{\lambda_{1}}{2\lambda_{2}})x_{\cdot
j}&\text{if}\ \ |\eta_{j}(\alpha)|\geq\eta_{0}\\\ 0&\text{if}\ \
|\eta_{j}(\alpha)|<\eta_{0}\end{cases}$ $\displaystyle=$
$\displaystyle\begin{cases}\beta_{j}x_{\cdot j}&\text{if}\ \
|\eta_{j}(\alpha)|\geq\eta_{0}\\\ 0&\text{if}\ \
|\eta_{j}(\alpha)|<\eta_{0}\end{cases}.$

Hence,

$\displaystyle\Psi_{j}^{\prime\prime}(\alpha)=$
$\displaystyle\begin{cases}-\frac{1}{2\lambda_{2}}x_{\cdot j}x^{\top}_{\cdot
j}&\text{if}\ \ |\eta_{j}(\alpha)|\geq\eta_{0}\\\ 0&\text{if}\ \
|\eta_{j}(\alpha)|<\eta_{0}\end{cases}.$

Therefore,

$\displaystyle\Psi^{\prime\prime}(\alpha)=-\frac{1}{2\lambda_{2}}\sum_{j\in
S}x_{\cdot j}x^{\top}_{\cdot j}.$

Here $S$ is the support set regarding $\alpha$, i.e.
$S=\mathrm{supp}\big{(}\beta(\alpha)\big{)}=\\{j\big{|}|\eta_{j}(\alpha)|\geq\eta_{0}\\}=\\{j\big{|}|x_{\cdot
j}^{\top}\alpha|\geq 2\lambda_{2}\eta_{0}\\}$. The smallest eigenvalue of
Hessian matrix $\Psi^{\prime\prime}(\alpha)$ depends on $X_{S}$. Let
$\sigma_{min}(X_{S})$ be the smallest eigenvalue of $X_{S}$, then
$\Psi(\alpha)$ is concave, and it is also
$\nu=\frac{\sigma_{min}(X_{S})}{2\lambda_{2}}$-strongly concave at point
$\alpha$.

With $l^{*}_{i}$ $\mu$-strongly convex, $D(\alpha)$ is $\gamma$-strongly
concave with
$\gamma=\mu+\inf_{S}\frac{\sigma_{min}(X_{S})}{2\lambda_{2}}\geq\mu$. Then we
have

$\displaystyle D(\alpha_{1})\leq
D(\alpha_{2})+\langle\nabla_{\alpha}D(\alpha_{2}),\alpha_{1}-\alpha_{2}\rangle-\frac{\gamma}{2}||\alpha_{1}-\alpha_{2}||^{2}.$

Let $\alpha_{2}=\bar{\alpha}$, and $\alpha_{1}=\alpha\in\mathcal{F}^{n}$. As
$\bar{\alpha}$ maximizes $D(\alpha)$,
$\langle\nabla_{\alpha}D(\bar{\alpha}),\alpha-\bar{\alpha}\rangle\leq 0$. It
implies

$\displaystyle D(\alpha)\leq
D(\bar{\alpha})-\frac{\gamma}{2}||\alpha-\bar{\alpha}||^{2}.$

Thus we have a ball range for the dual variable

$\displaystyle\forall\alpha\in\mathcal{F}^{n},\beta\in\mathbb{R}^{p},||\alpha-\bar{\alpha}||_{2}\leq\sqrt{\frac{2(P(\beta)-D(\alpha))}{\gamma}}=:r.$

This completes the proof. ∎

##  Appendix B Algorithm Analysis

In this section, we present the complexity analysis of the inner updating
Algorithm 1.

###  B.1 Convergence of Inner Primal-dual Updating Algorithm

We will show that under certain conditions $\beta(\alpha)$ is locally smooth
around $\bar{\beta}=\beta(\bar{\alpha})$. For a given set of parameters
$\lambda=\\{\lambda_{0},\lambda_{1},\lambda_{2}\\}$, $\beta(\alpha)$
corresponds to a set of support features $\mathrm{supp}(\beta(\alpha))$. We
use $\bar{\delta}$ to represent the set in the dual feasible space that
$\mathrm{supp}(\beta(\alpha))=\mathrm{supp}(\bar{\beta})$.

######  Lemma B.1.

Let $X=[x_{1},...,x_{n}]^{T}\in\mathbb{R}^{n\times p}$ be the data matrix,
$\eta_{j}(\bar{\alpha})$ be the $j$th entry of $\eta(\bar{\alpha})$,
$S=\mathrm{supp}(\bar{\beta})$, and
$N=\\{j|\eta_{j}(\bar{\alpha})=\eta_{0(j)}\\}$. Assume that
$\\{l_{i}\\}_{i=1,...,n}$ are differentiable, and let

$\displaystyle\bar{\delta}=:2\lambda_{2}\min\bigg{\\{}\min_{j:j\in
S}\frac{|\eta_{j}(\bar{\alpha})|-\eta_{0}}{||x_{\cdot j}||},\min_{j:j\in
S^{c}\setminus N}\frac{\eta_{0}-|\eta_{j}(\bar{\alpha})|}{||x_{\cdot
j}||}\bigg{\\}},$

with $||\alpha-\bar{\alpha}||\leq\bar{\delta}$, we have
$\mathrm{supp}(\beta(\alpha))=\mathrm{supp}(\bar{\beta})$, and
$||\beta(\alpha)-\bar{\beta}||\leq\frac{\sigma_{max}(X_{S})}{2\lambda_{2}}||\alpha-\bar{\alpha}||$.

###### Proof.

For any $\alpha$, we have

$\displaystyle\eta(\alpha)=-\frac{1}{2\lambda_{2}}\sum_{i=1}^{n}\alpha_{i}x_{i}=-\frac{1}{2\lambda_{2}}X^{\top}\alpha.$
(28)

For a feature $j\in S=\mathrm{supp}(\bar{\beta})$, we have
$|\eta_{j}(\bar{\alpha})|-\eta_{0}>0$, which is

$\displaystyle|\eta_{j}(\bar{\alpha})+\eta_{j}(\alpha)-\eta_{j}(\bar{\alpha}))|>\eta_{0}.$

We try to find the space for $\alpha$s that have the same support as
$\bar{\alpha}$. We use the lower bound of above inequality,

$\displaystyle|\eta_{j}(\bar{\alpha})+\eta_{j}(\alpha)-\eta_{j}(\bar{\alpha}))|\geq|\eta_{j}(\bar{\alpha})|-|\eta_{j}(\alpha)-\eta_{j}(\bar{\alpha})|>\eta_{0},$

yields

$\displaystyle|\eta_{j}(\alpha)-\eta_{j}(\bar{\alpha})|<|\eta_{j}(\bar{\alpha})|-\eta_{0}.$

With (28),

$\displaystyle\frac{\big{|}\big{|}\alpha-\bar{\alpha}\big{|}\big{|}}{2\lambda_{2}}\big{|}\big{|}x_{\cdot
j}\big{|}\big{|}<|\eta_{j}(\bar{\alpha})|-\eta_{0}.$

Hence

$\displaystyle\big{|}\big{|}\alpha-\bar{\alpha}\big{|}\big{|}<\min_{j:j\in
S}\frac{2\lambda_{2}(|\eta_{j}(\bar{\alpha})|-\eta_{0})}{||x_{\cdot j}||}.$

Similarly, for features $j$, $j\notin S$ and $j\notin N$,

$\displaystyle|\eta_{j}(\bar{\alpha})|<\eta_{0},$

yields

$\displaystyle|\eta_{j}(\alpha)-\eta_{j}(\bar{\alpha})|<\eta_{0}-|\eta_{j}(\bar{\alpha})|.$

With all $j$s, $j\in S^{c}\setminus N$

$\displaystyle\big{|}\big{|}\alpha-\bar{\alpha}\big{|}\big{|}<\min_{j:j\in
S^{c}\setminus
N}\frac{2\lambda_{2}(\eta_{0}-|\eta_{j}(\bar{\alpha})|)}{||x_{\cdot j}||}.$

Therefore, if $\alpha\in\bar{\delta}$,with

$\displaystyle||\alpha-\bar{\alpha}||\leq\bar{\delta}=2\lambda_{2}\min\bigg{\\{}\min_{j:j\in
S}\frac{|\eta_{j}(\bar{\alpha})|-\eta_{0}}{||x_{\cdot j}||},\min_{j:j\in
S^{c}\setminus N}\frac{\eta_{0}-|\eta_{j}(\bar{\alpha})|}{||x_{\cdot
j}||}\bigg{\\}}$

we have $\mathrm{supp}(\beta(\alpha))=\mathrm{supp}(\bar{\beta})$. With
$||\alpha-\bar{\alpha}||\leq\bar{\delta}$, the primal problem becomes a convex
$\ell_{1}$ regularization problem without any redundant features. The super-
gradient in Remark 2.2 becomes
$g_{\alpha}=X_{S}\beta_{S}(\alpha)-l^{*^{\prime}}(\alpha)$. With
$||g_{\alpha}||\rightarrow 0$,
$\beta_{S}(\alpha)=(X_{S}^{\top}X_{S})^{-1}X_{S}^{\top}l^{*^{\prime}}(\alpha)$.
As $S$ is fixed, with $||\alpha-\bar{\alpha}||$ a small value, we have
$\mathrm{sign}(\beta(\alpha))=\mathrm{sign}(\bar{\beta})$. It means

$\displaystyle||\beta(\alpha)-\bar{\beta}||$
$\displaystyle=||\mathfrak{B}(\eta(\alpha))-\mathfrak{B}(\eta(\bar{\alpha}))||$
$\displaystyle\leq||\eta(\alpha)-\eta(\bar{\alpha})||=\frac{1}{2\lambda_{2}}||X_{S}(\alpha-\bar{\alpha})||$
$\displaystyle\leq\frac{\sigma_{max}(X_{S})}{2\lambda_{2}}||\alpha-\bar{\alpha}||.$

It finishes the proof of the lemma. ∎

Note that the above lemma can be extended to any pair of
$\alpha_{1},\alpha_{2}\in\mathcal{F}^{n}$, and if they are close enough, they
have the same support set. Let
$\Psi(\alpha;\lambda_{0},\lambda_{1},\lambda_{2})=\sum_{j=1}^{p}\Psi(\eta_{j}(\alpha);\lambda_{0},\lambda_{1},\lambda_{2})$,
and it is easy to verify that
$\Psi(\alpha;\lambda_{0},\lambda_{1},\lambda_{2})$ is concave.

######  Lemma B.2.

Assume that the primal loss functions $\\{l_{i}()\\}_{i=1}^{n}$ are
$1/\mu$-strongly smooth. Then the following inequality holds for any
$\alpha_{1},\alpha_{2}\in\mathcal{F}^{n}$, and $g(\alpha_{2})\in\partial
D(\alpha_{2})$:

$\displaystyle D(\alpha_{1})\leq D(\alpha_{2})+\langle
g_{\alpha_{2}},\alpha_{1}-\alpha_{2}\rangle-\frac{\gamma}{2}||\alpha_{1}-\alpha_{2}||^{2}.$

Moreover, $\forall\alpha\in\mathcal{F}^{n}$, and $g_{\alpha}\in\partial
D(\alpha)$, $||\alpha-\bar{\alpha}||\leq\sqrt{\frac{2\langle
g_{\alpha},\bar{\alpha}-\alpha\rangle}{\gamma}}.$ Here, $\gamma$ is the same
in Theorem 2.2.

###### Proof.

With the assumption $l_{i}$ being $1/\mu$-smooth, its conjugate function
$l^{*}_{i}$ is $\mu$-strongly convex. With the dual problem

$\displaystyle D(\alpha)$
$\displaystyle=-\sum_{i=1}^{n}l_{i}^{*}(\alpha_{i})+\sum_{j=1}^{p}\Psi(\eta_{j}(\alpha);\lambda_{0},\lambda_{1},\lambda_{2}).$

Here

$\displaystyle\Psi(\eta_{j}(\alpha);\lambda_{0},\lambda_{1},\lambda_{2})=\begin{cases}-\lambda_{2}(|\eta_{j}|-\frac{\lambda_{1}}{2\lambda_{2}})^{2}+\lambda_{0}&\text{if}\
\ |\eta_{j}(\alpha)|\geq\eta_{0}\\\ 0&\text{if}\ \
|\eta_{j}(\alpha)|<\eta_{0}\end{cases}.$

According to the proof of Theorem 2.2,
$\Psi(\alpha)=\sum_{j=1}^{p}\Psi(\eta_{j}(\alpha);\lambda_{0},\lambda_{1},\lambda_{2})$
is $\nu$-strongly concave with
$\nu=\inf_{\alpha}\nu(\alpha)=\frac{\sigma_{min}(X_{S_{\alpha}})}{2\lambda}$,
and $S_{\alpha}=\mathrm{supp}\big{(}\beta(\alpha)\big{)}$. When
$S_{\alpha}=\emptyset$, we have $\nu(\alpha)=0$. Now let us consider two
arbitrary dual variables $\alpha_{1},\alpha_{2}\in\mathcal{F}^{n}$,

$\displaystyle\Psi(\alpha_{1})\leq\Psi(\alpha_{2})+\Psi^{\prime}(\alpha_{2})^{\top}(\alpha_{1}-\alpha_{2})-\frac{\nu}{2}||\alpha_{1}-\alpha_{2}||^{2}.$

Hence,

$\displaystyle D(\alpha_{1})=$
$\displaystyle-\sum_{i=1}^{n}l_{i}^{*}(\alpha_{1(i)})+\sum_{j=1}^{p}\Psi(\eta_{j}(\alpha_{1});\lambda_{0},\lambda_{1},\lambda_{2})$
$\displaystyle\leq$
$\displaystyle\sum_{i=1}^{n}\big{(}-l_{i}^{*}(\alpha_{2(i)})-l_{i}^{*^{\prime}}(\alpha_{2(i)})(\alpha_{1(i)}-\alpha_{2(i)})-\frac{\mu}{2}(\alpha_{1(i)}-\alpha_{2(i)})^{2}\big{)}$
$\displaystyle+\Psi(\alpha_{2})+\Psi^{\prime}(\alpha_{2})^{\top}(\alpha_{1}-\alpha_{2})-\frac{\nu}{2}||\alpha_{1}-\alpha_{2}||^{2}$
$\displaystyle\leq$ $\displaystyle D(\alpha_{2})+\langle
g_{\alpha_{2}},\alpha_{1}-\alpha_{2}\rangle-\frac{\gamma}{2}||\alpha_{1}-\alpha_{2}||^{2}.$
(29)

Here $\alpha_{1(i)}$ is the $i$th entry of $\alpha_{1}$. This proves the first
desirable inequality in the lemma. With the above inequality and using the
fact $D(\alpha)\leq D(\bar{\alpha})$ we get that

$\displaystyle D(\bar{\alpha})\leq$ $\displaystyle D(\alpha)+\langle
g_{\alpha},\bar{\alpha}-\alpha\rangle-\frac{\gamma}{2}||\alpha-\bar{\alpha}||^{2}\leq
D(\bar{\alpha})+\langle
g_{\alpha},\bar{\alpha}-\alpha\rangle-\frac{\gamma}{2}||\alpha-\bar{\alpha}||^{2},$

which leads to the second desired bound,

$\displaystyle||\alpha-\bar{\alpha}||\leq\sqrt{\frac{2\langle
g_{\alpha},\bar{\alpha}-\alpha\rangle}{\mu+\nu(\alpha)}}.$

It concludes the proof of the lemma ∎

Different from the primal updating [Hazimeh and Mazumder, 2020] or dual
updating [Liu et al., 2017] algorithms, Algorithm 1 has both primal and dual
updating steps. Let $m_{1}=\max_{j:1\leq j\leq p}|y^{\top}x_{\cdot j}|$, and
$m_{2}=\max_{j:1\leq j\leq p}||X^{\top}x_{\cdot j}||$,
$m_{3}=\max_{i,\alpha_{i}^{t}\in\mathcal{F}}|l_{i}^{*^{\prime}}(\alpha_{i}^{t})|$,
and $\varrho=\sqrt{n}||\alpha^{t}||_{\infty}-\lambda_{1}$. We have the
following theorem regarding the convergence of Algorithm 1.

Theorem 4.1 Assume that $l_{i}$ is $1/\mu$-smooth, $||x_{i}||\leq\vartheta\
\forall 1\leq i\leq n$, and $||x_{\cdot j}||=1\ \forall 1\leq j\leq p$. By
choosing $w_{t}=\frac{1}{t\gamma}$, then the sequence generated by Algorithm 1
satisfies the following estimation error inequality:

$\displaystyle||\alpha^{t}-\bar{\alpha}||^{2}\leq
c_{1}\bigg{(}\frac{1}{t}+\frac{\ln t}{t}\bigg{)}.$

Here $c_{1}=\frac{c_{0}^{2}}{\mu^{2}}$.
$c_{0}=\frac{\sqrt{np}\vartheta}{2\lambda_{2}(1+2\lambda_{2})}(2\lambda_{2}m_{1}+\varrho+\sqrt{n}m_{2}\varrho-2\lambda_{1}\lambda_{2})+\sqrt{n}m_{3}$,
$gamma$ is same as in Theorem 2.2.

###### Proof.

Let us consider $g^{t}$,
$g_{i}^{t}=x_{i}^{\top}\beta^{t}-l_{i}^{*^{\prime}}(\alpha_{i}^{t})$. After
computing the primal $\beta^{t}$ with the primal-dual relation (14), Algorithm
1 also performs primal coordinate descent starting with $\beta^{t}$ using (16)
to the improve super-gradient $g^{t}$.

Let $\breve{\beta}^{t}$ be the output of operation (14) at step $t$. From the
expression of $\beta^{t}$ (14), if $\breve{\beta}_{j}^{t}\neq 0$,
$\breve{\beta}^{t}_{j}(\alpha^{t})=\mathrm{sign}\big{(}\eta_{j}(\alpha^{t})\big{)}\big{(}|\eta_{j}(\alpha^{t})|-\frac{\lambda_{1}}{2\lambda_{2}}\big{)}$.
With $||x_{\cdot j}||=1$, $|\eta_{j}(\alpha^{t})|=\frac{|x_{\cdot
j}^{\top}\alpha^{t}|}{2\lambda_{2}}\leq\frac{\sqrt{n}\varphi}{2\lambda_{2}}$.
Here $\varphi=||\alpha^{t}||_{\infty}$. Then we have

$\displaystyle|\breve{\beta}_{j}^{t}|\leq|\eta_{j}(\alpha^{t})|-\frac{\lambda_{1}}{2\lambda_{2}}\leq\frac{\sqrt{n}\varphi-\lambda_{1}}{2\lambda_{2}}=\frac{\varrho}{2\lambda_{2}},$
(30)

with $\varrho=\sqrt{n}\varphi-\lambda_{1}$.

According to (16), with $\beta$ as the input, the non-zero output at entry $j$

$\displaystyle\grave{\beta}_{j}=T(\beta;\lambda_{0},\lambda_{1},\lambda_{2})=\mathrm{sign}(\tilde{\mathbf{\beta}}_{j})\frac{|\tilde{\mathbf{\beta}}_{j}|-\lambda_{1}}{1+2\lambda_{2}},$

with

$\displaystyle\tilde{\mathbf{\beta}}_{j}=\big{\langle}y-\sum_{i:i\neq j,i\in
S}x_{\cdot i}\mathbf{\beta}_{i},x_{\cdot
j}\big{\rangle}=(y-X\mathbf{\beta})^{\top}x_{\cdot
j}+\mathbf{\beta}_{j}x_{\cdot j}^{\top}x_{\cdot j}=y^{\top}x_{\cdot
j}-\beta^{\top}X^{\top}x_{\cdot j}+\mathbf{\beta}_{j}$

and $S=\mathrm{supp}(\beta)$. With $\beta_{j}\neq 0$ we have
$|\tilde{\mathbf{\beta}}_{j}|-\lambda_{1}\geq 0$ . Then

$\displaystyle|\grave{\beta}_{j}|$
$\displaystyle=\frac{1}{1+2\lambda_{2}}\big{(}|\tilde{\mathbf{\beta}}_{j}|-\lambda_{1}\big{)}\leq\frac{1}{1+2\lambda_{2}}\big{(}|y^{\top}x_{\cdot
j}|+|\beta^{\top}X^{\top}x_{\cdot j}|+|\mathbf{\beta}_{j}|-\lambda_{1}\big{)}$
$\displaystyle\leq\frac{1}{1+2\lambda_{2}}\big{(}|y^{\top}x_{\cdot
j}|+\sqrt{\beta^{\top}\beta x_{\cdot j}^{\top}XX^{\top}x_{\cdot
j}}+|\beta_{j}|-\lambda_{1}\big{)}.$

Let input $\beta=\breve{\beta}^{t}$, with (30) the upper bound of the output
after one round coordinate descent will be

$\displaystyle|\beta^{t}_{j}|$
$\displaystyle\leq\frac{1}{1+2\lambda_{2}}\big{(}|y^{\top}x_{\cdot
j}|+\sqrt{\breve{\beta}^{\top}\breve{\beta}x_{\cdot j}^{\top}XX^{\top}x_{\cdot
j}}+|\breve{\beta}_{j}|-\lambda_{1}\big{)}$
$\displaystyle\leq\frac{1}{1+2\lambda_{2}}\big{(}|y^{\top}x_{\cdot
j}|+\frac{1+\sqrt{n}||X^{\top}x_{\cdot
j}||}{2\lambda_{2}}\varrho-\lambda_{1}\big{)}$
$\displaystyle\leq\frac{1}{2\lambda_{2}(1+2\lambda_{2})}(2\lambda_{2}m_{1}+\varrho+\sqrt{n}m_{2}\varrho-2\lambda_{1}\lambda_{2}):=\psi.$

Here $m_{1}=\max_{j:1\leq j\leq p}|y^{\top}x_{\cdot j}|$, and
$m_{2}=\max_{j:1\leq j\leq p}||X^{\top}x_{\cdot j}||$. Let
$m_{3}=max_{i,\alpha_{i}^{t}\in\mathcal{F}}|l_{i}^{*^{\prime}}(\alpha_{i}^{t})|$.
Then

$\displaystyle|g_{i}^{t}|\leq|x_{i}^{\top}\beta^{t}|+|l_{i}^{*^{\prime}}(\alpha_{i}^{t})|\leq\sqrt{||x_{i}||^{2}||\beta^{t}||^{2}}+m_{3}\leq\sqrt{p}\psi\vartheta+m_{3}\
,\quad\forall 1\leq i\leq n.$

Hence,

$\displaystyle||g^{t}||\leq c_{0}:=\sqrt{np}\psi\vartheta+\sqrt{n}m_{3}.$ (31)

Let $h^{t}=||\alpha^{t}-\bar{\alpha}||$ and $v^{t}=\langle
g^{t},\bar{\alpha}-\alpha^{t}\rangle$. The concavity of $D$ implies $v^{t}\geq
0$. According to Lemma B.2,

$\displaystyle
h^{t}=||\alpha^{t}-\bar{\alpha}||\leq\sqrt{\frac{2v^{t}}{\gamma}}.$ (32)

Let $\omega^{t}$ be the step size of dual variable at step $t$. Then we have

$\displaystyle(h^{t})^{2}=$ $\displaystyle||\alpha^{t}-\bar{\alpha}||^{2}$
$\displaystyle=$
$\displaystyle||P_{\mathcal{F}^{n}}\big{(}\alpha^{t-1}+\omega^{t-1}g^{t-1}\big{)}-\bar{\alpha}||^{2}$
$\displaystyle\leq$
$\displaystyle||\alpha^{t-1}+\omega^{t-1}g^{t-1}-\bar{\alpha}||^{2}$
$\displaystyle=$
$\displaystyle(h^{t-1})^{2}-2\omega^{t-1}v^{t-1}+(\omega^{t-1})^{2}||g^{t-1}||^{2}$
$\displaystyle\leq$
$\displaystyle(h^{t-1})^{2}-\omega^{t-1}(\gamma)(h^{t-1})^{2}+(\omega^{t-1})^{2}c_{0}^{2}.$

The last step is due to (32) and (31). Let $\omega^{t-1}=\frac{1}{\gamma t}$.
Then we get

$\displaystyle(h^{t})^{2}\leq\big{(}1-\frac{1}{t}\big{)}(h^{t-1})^{2}+\frac{c_{0}^{2}}{\big{(}\gamma\big{)}^{2}t^{2}}.$

Recursively applying the above inequality we get

$\displaystyle(h^{t})^{2}\leq\frac{c_{0}^{2}}{\gamma^{2}}\bigg{(}\frac{1}{t}+\frac{\ln
t}{t}\bigg{)}\leq c_{1}\bigg{(}\frac{1}{t}+\frac{\ln t}{t}\bigg{)}.$

Here
$c_{0}=\frac{\vartheta\sqrt{np}}{2\lambda_{2}(1+2\lambda_{2})}(2\lambda_{2}m_{1}+\varrho+\sqrt{n}m_{2}\varrho-2\lambda_{1}\lambda_{2})+\sqrt{n}m_{3}$.
This proves the bound in the theorem. ∎

Different from the complexity analysis in Liu et al. [2017] regarding hard
thresholding, the analysis presented here is based on the primal-dual problem
structures in (1) and (12). Moreover, the study in Liu et al. [2017] focus on
the dual updating steps regarding hard thresholding. Whereas Theorem 4.1
includes the complexity of both primal and dual updating steps given in
Algorithm 1. We further prove the convergence of primal variable and the
duality gap.

######  Lemma B.3.

For $\alpha\in\mathcal{F}^{n}$, with $\beta=\beta(\alpha)$ the primal-dual gap
can be written as

$\displaystyle
P(\beta)-D(\alpha)=\sum_{i=1}^{n}\bigg{(}l(\beta^{\top}x_{i},y_{i})+l_{i}^{*}(\alpha_{i})\bigg{)}-\alpha^{\top}X\beta\
.$ (33)

Moreover, with $\theta\in[\partial l_{1}(\beta^{\top}x_{1}),\partial
l_{2}(\beta^{\top}x_{2}),...,\partial l_{n}(\beta^{\top}x_{n})]^{\top}$, we
have

$\displaystyle P(\beta)-D(\alpha)\leq\langle g_{\alpha},\theta-\alpha\rangle.$
(34)

###### Proof.

With the definitions of $P(\beta)$, $D(\alpha)$, and (27)

$\displaystyle P(\beta)-D(\alpha)$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{n}l(\beta^{\top}x_{i},y_{i})+\lambda_{1}||\beta||_{1}+\lambda_{2}||\beta||_{2}^{2}+\lambda_{0}||\beta||_{0}-\bigg{(}-\sum_{i=1}^{n}l_{i}^{*}(\alpha_{i})+\sum_{j=1}^{p}\Psi(\eta_{j}(\alpha);\lambda_{0},\lambda_{1},\lambda_{2})\bigg{)}$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{n}l(\beta^{\top}x_{i},y_{i})+\lambda_{1}||\beta||_{1}+\lambda_{2}||\beta||_{2}^{2}+\lambda_{0}||\beta||_{0}-\bigg{(}-\sum_{i=1}^{n}l_{i}^{*}(\alpha_{i})-\lambda_{2}||\beta||^{2}_{2}+\lambda_{0}||\beta||_{0}\bigg{)}$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{n}\bigg{(}l(\beta^{\top}x_{i},y_{i})+l_{i}^{*}(\alpha_{i})\bigg{)}+\lambda_{1}||\beta||_{1}+2\lambda_{2}||\beta||_{2}^{2}\
.$ (35)

Let $S=\mathrm{supp}(\beta)$, the last two terms in above equation can
rewritten as

$\displaystyle\lambda_{1}||\beta||_{1}+2\lambda_{2}||\beta||_{2}^{2}=$
$\displaystyle
2\lambda_{2}\bigg{(}\frac{\lambda_{1}}{2\lambda_{2}}||\beta||_{1}+||\beta||_{2}^{2}\bigg{)}$
$\displaystyle=$ $\displaystyle 2\lambda_{2}\sum_{j\in
S}\bigg{(}\frac{\lambda_{1}}{2\lambda_{2}}(|\eta_{j}(\alpha)|-\frac{\lambda_{1}}{2\lambda_{2}})+(|\eta_{j}(\alpha)|-\frac{\lambda_{1}}{2\lambda_{2}})^{2}\bigg{)}$
$\displaystyle=$ $\displaystyle 2\lambda_{2}\sum_{j\in
S}\bigg{(}(|\eta_{j}(\alpha)|-\frac{\lambda_{1}}{4\lambda_{2}})^{2}-\big{(}\frac{\lambda_{1}}{4\lambda_{2}}\big{)}^{2}\bigg{)}$
$\displaystyle=$ $\displaystyle 2\lambda_{2}\sum_{j\in
S}|\eta_{j}(\alpha)|\bigg{(}|\eta_{j}(\alpha)|-\frac{\lambda_{1}}{2\lambda_{2}}\bigg{)}$
$\displaystyle\overset{}{=}$ $\displaystyle 2\lambda_{2}\sum_{j\in
S}\eta_{j}(\alpha)\mathrm{sign}(\eta_{j}(\alpha))\bigg{(}|\eta_{j}(\alpha)|-\frac{\lambda_{1}}{2\lambda_{2}}\bigg{)}$
$\displaystyle=$ $\displaystyle-\alpha^{\top}X\beta\ .$ (36)

The last step is due to the definitions of $\eta(\alpha)$ and $\beta(\alpha)$
(6). With equations (35) and (36), we prove equation (33).

For $i\in[n]$, with the maximizing argument property of convex conjugate

$\displaystyle
l_{i}(\beta^{\top}x_{i})=\beta^{\top}x_{i}l^{\prime}_{i}(\beta^{\top}x_{i})-l^{*}_{i}(l^{\prime}_{i}(\beta^{\top}x_{i})),$
$\displaystyle
l^{*}_{i}(\alpha_{i})=\alpha_{i}l^{*^{\prime}}_{i}(\alpha_{i})-l_{i}(l^{*^{\prime}}_{i}(\alpha_{i}))\
.$

Adding above two equations, and applying Fenchel-Young inequality, we get

$\displaystyle l_{i}(\beta^{\top}x_{i})+l^{*}_{i}(\alpha_{i})=$
$\displaystyle\beta^{\top}x_{i}l^{\prime}_{i}(\beta^{\top}x_{i})+\alpha_{i}l^{*^{\prime}}_{i}(\alpha_{i})-\bigg{(}l^{*}_{i}(l^{\prime}_{i}(\beta^{\top}x_{i}))+l_{i}(l^{*^{\prime}}_{i}(\alpha_{i}))\bigg{)}$
$\displaystyle\leq$
$\displaystyle\beta^{\top}x_{i}l^{\prime}_{i}(\beta^{\top}x_{i})+\alpha_{i}l^{*^{\prime}}_{i}(\alpha_{i})-l^{*^{\prime}}_{i}(\alpha_{i})l^{\prime}_{i}(\beta^{\top}x_{i}).$
(37)

Moreover, with $\theta\in[\partial l_{1}(\beta^{\top}x_{1}),\partial
l_{2}(\beta^{\top}x_{2}),...,\partial l_{n}(\beta^{\top}x_{n})]^{\top}$

$\displaystyle\langle g_{\alpha},\theta-\alpha\rangle$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{n}\big{(}\beta^{\top}x_{i}-l_{i}^{*^{\prime}}(\alpha_{i})\big{)}\big{(}l_{i}^{\prime}(\beta^{\top}x_{i})-\alpha_{i}\big{)}$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{n}\bigg{(}\beta^{\top}x_{i}l^{\prime}_{i}(\beta^{\top}x_{i})+\alpha_{i}l^{*^{\prime}}_{i}(\alpha_{i})-l^{*^{\prime}}_{i}(\alpha_{i})l^{\prime}_{i}(\beta^{\top}x_{i})-\alpha_{i}\beta^{\top}x_{i}\bigg{)}$
$\displaystyle\geq$
$\displaystyle\sum_{i=1}^{n}\bigg{(}l_{i}(\beta^{\top}x_{i})+l^{*}_{i}(\alpha_{i})-\alpha_{i}\beta^{\top}x_{i}\bigg{)}$
(38) $\displaystyle=$ $\displaystyle P(\beta)-D(\alpha).$ (39)

Step (38) is due to (37), and step (39) is by (33). It concludes the lemma. ∎

Theorem 4.2 Assume that $l_{i}$ is $1/\mu$-smooth, $||x_{i}||\leq\vartheta\
\forall 1\leq i\leq n$, and $||x_{\cdot j}||=1\ \forall 1\leq j\leq p$. Let
$t_{1}=\frac{3c_{1}}{\bar{\delta}^{2}}\log\frac{3c_{1}}{\bar{\delta}^{2}}$,
with $t>t_{1}$, we have
$\mathrm{supp}(\beta(\alpha))=\mathrm{supp}(\bar{\beta})$ and
$||\beta(\alpha)-\bar{\beta}||\leq\frac{\sigma_{max}(X_{S})}{2\lambda_{2}}||\alpha-\bar{\alpha}||$.
Moreover, let
$t_{2}=\frac{3c_{1}c_{2}}{\epsilon}\log\frac{3c_{1}c_{2}}{\epsilon}$,
$c_{2}=c_{0}\bigg{(}1+\frac{\sigma_{max}(X_{S})}{2\mu\lambda_{2}}\bigg{)}$,
for any $\epsilon>0$ with $t>\mathrm{max}\\{t_{1},t_{2}\\}$, we have
$P(\beta^{t})-D(\alpha^{t})\leq\epsilon$.

###### Proof.

With Theorem 4.1, let $||\alpha-\bar{\alpha}||^{2}\leq
c_{1}\bigg{(}\frac{1}{t}+\frac{\ln t}{t}\bigg{)}\leq\bar{\delta}^{2}$, and it
implies that
$t\geq\frac{3c_{1}}{\bar{\delta}^{2}}\log\frac{3c_{1}}{\bar{\delta}^{2}}$.
Hence, with $t\geq
t_{1}:=\frac{3c_{1}}{\bar{\delta}^{2}}\log\frac{3c_{1}}{\bar{\delta}^{2}}$,
$\mathrm{supp}(\beta(\alpha))=\mathrm{supp}(\bar{\beta})$ and
$||\beta(\alpha)-\bar{\beta}||\leq\frac{\sigma_{max}(X_{S})}{2\lambda_{2}}||\alpha-\bar{\alpha}||$
according to Lemma B.1.

Let $\theta^{t}=[\partial l_{1}((\beta^{t})^{\top}x_{1}),\partial
l_{2}((\beta^{t})^{\top}x_{2}),...,\partial
l_{n}((\beta^{t})^{\top}x_{n})]^{\top}$. According to Lemma B.3,

$\displaystyle P(\beta^{t})-D(\alpha^{t})\leq$ $\displaystyle\langle
g^{t}_{\alpha^{t}},\theta^{t}-\alpha^{t}\rangle$ $\displaystyle\leq$
$\displaystyle||g^{t}_{\alpha^{t}}||\big{(}||\theta^{t}-\bar{\alpha}||+||\bar{\alpha}-\alpha^{t}||\big{)}.$
(40)

As $\bar{\alpha}\in[\partial l_{1}(\bar{\beta}^{\top}x_{1}),\partial
l_{2}(\bar{\beta}^{\top}x_{2}),...,\partial
l_{n}(\bar{\beta}^{\top}x_{n})]^{\top}$, $l_{i}()$ is $\frac{1}{\mu}$ smooth,
and $||x_{\cdot j}||=1$, we have

$\displaystyle||\theta^{t}-\bar{\alpha}||\leq\frac{1}{\mu}||\beta^{t}-\bar{\beta}||\leq\frac{\sigma_{max}(X_{S})}{2\mu\lambda_{2}}||\alpha^{t}-\bar{\alpha}||.$

Therefore,

$\displaystyle P(\beta^{t})-D(\alpha^{t})\leq$
$\displaystyle||g^{t}_{\alpha^{t}}||\big{(}||\theta^{t}-\bar{\alpha}||+||\bar{\alpha}-\alpha^{t}||\big{)}$
$\displaystyle\leq$ $\displaystyle
c_{0}\bigg{(}1+\frac{\sigma_{max}(X_{S})}{2\mu\lambda_{2}}\bigg{)}||\alpha^{t}-\bar{\alpha}||.$
(41)

If
$t>\mathrm{max}\\{t_{1},t_{2}\\},t_{2}:=\frac{3c_{1}c_{2}^{2}}{\epsilon^{2}}\log\frac{3c_{1}c_{2}^{2}}{\epsilon^{2}}$,
and
$c_{2}:=c_{0}\bigg{(}1+\frac{\sigma_{max}(X_{S})}{2\mu\lambda_{2}}\bigg{)}$,
we have
$||\alpha^{t}-\bar{\alpha}||\leq\frac{\epsilon}{c_{0}\bigg{(}1+\frac{\sigma_{max}(X_{S})}{2\mu\lambda_{2}}\bigg{)}}$,
and with (41), $P(\beta^{t})-D(\alpha^{t})\leq\epsilon$. ∎

##  Appendix C Dual Problems of Some Loss Functions

###  C.1 Logistic Loss

The primal form of logistic regression is given by

$\displaystyle\min_{\beta\in\mathbb{R}^{p}}F(\beta)$
$\displaystyle=f(\beta)+\lambda_{0}||\beta||_{0},$ (42)
$\displaystyle~{}~{}~{}~{}f(\beta)$
$\displaystyle=\sum_{i=1}^{n}\big{(}-y_{i}x_{i}^{\top}\beta+\log(1+\exp(x_{i}^{\top}\beta))\big{)}+\lambda_{1}||\beta||_{1}+\lambda_{2}||\beta||_{2}^{2}.$

Here $l_{i}(u;y_{i})=-y_{i}u+\log(1+\exp(u))$, and
$l_{i}^{*}(\alpha_{i})=(\alpha_{i}+y_{i})\log(\alpha_{i}+y_{i})+(1-\alpha_{i}-y_{i})\log(1-\alpha_{i}-y_{i}),\alpha_{i}+y_{i}\in[0,1]$,
here $u=\beta^{\top}x_{i}$. The dual objective is

$\displaystyle D(\alpha)=$
$\displaystyle-(\alpha+y)\log(\alpha+y)-(1-\alpha-y)\log(1-\alpha-y)+\sum_{j=1}^{p}\Psi(-\frac{1}{2\lambda_{2}}\sum_{i=1}^{n}\bar{\alpha}_{i}x_{i};\lambda_{0},\lambda_{1},\lambda_{2}).$

Here $\Psi()$ is given by (13). With $\alpha_{i}+y_{i}\in[0,1]$, the dual
feasible project operator is
$\mathcal{P}_{\mathcal{F}}(\alpha_{i})=\begin{cases}\alpha_{i}\quad\ \quad-
y_{i}\leq\alpha\leq 1-y_{i}\\\ 0\quad\ \ \quad\ \alpha\leq-y_{i}\\\ 1\quad\ \
\quad\ \alpha\geq 1-y_{i}\\\ \end{cases}.$ The super gradient regarding
logistic regression is

$\displaystyle g_{\alpha}=\nabla_{\alpha}D(\alpha)=$
$\displaystyle\bigg{[}\beta^{\top}(\alpha)x_{1}-\log\big{(}\frac{\alpha_{1}+y_{1}}{1-\alpha_{1}-y_{1}}\big{)},...,\beta^{\top}\big{(}\alpha)x_{n}-\log(\frac{\alpha_{n}+y_{n}}{1-\alpha_{n}-y_{n}}\big{)}\bigg{]}^{\top}.$

There is no closed form of updating formula with coordinate descent regarding
the primal problem (42). We can apply proximal algorithm [Parikh and Boyd,
2014] to this type of primal loss functions.

###  C.2 Huber Loss

Consider a regression problem with Huber loss, i.e.,

$\displaystyle\min_{\beta\in\mathbb{R}^{p}}F(\beta)=\sum_{i=1}^{n}\left\\{l_{Huber}(y_{i}x_{i}^{T}\beta)\right\\}+\lambda_{0}\|\beta\|_{0}+\lambda_{1}\|\beta\|_{1}+\lambda_{2}\|\beta\|_{2},$
$\displaystyle
l_{Huber}(y_{i}x_{i}^{T}\beta)=\begin{cases}0,&y_{i}x_{i}^{T}\beta\geq 1\\\
1-y_{i}x_{i}^{T}\beta-\gamma/2,&y_{i}x_{i}^{T}\beta<1-\gamma\\\
\frac{1}{2\gamma}(1-y_{i}x_{i}^{T}\beta)^{2},&\textrm{otherwise}\end{cases}$

with $\gamma$ being some hyper tuning parameter. The dual function of
$l_{Huber}(\cdot)$ is

$\displaystyle
l^{\ast}_{Huber}(\alpha_{i})=\begin{cases}y_{i}\alpha_{i}+\frac{\gamma}{2}\gamma\alpha^{2},&-1\leq
y_{i}x_{i}^{T}\leq 0\\\ +\infty,&\textrm{otherwise}\end{cases}.$

Therefore, the corresponding Lagrangian function is
$L(\beta,\alpha)=\sum_{i}\\{\alpha_{i}x_{i}^{T}\beta-
l_{Huber}^{\ast}(\alpha_{i})\\}+\lambda_{0}||\beta||_{0}+\lambda_{1}||\beta||_{1}+\lambda_{2}||\beta||^{2}.$
The dual problem can be written as

$\displaystyle\max_{\alpha\in\mathcal{F}^{n}}-\sum_{i=1}^{n}l_{Huber}^{*}(\alpha_{i})+\sum_{j=1}^{p}\Psi(\eta_{j}(\alpha);\lambda_{0},\lambda_{1},\lambda_{2}).$