#  Linearly-Convergent FISTA Variant for Composite Optimization with Duality

Casey Garner Shuzhong Zhang School of Mathematics, University of Minnesota,
<EMAIL_ADDRESS>Department of Industrial and Systems Engineering,
<EMAIL_ADDRESS>
###### Abstract

Many large-scale optimization problems can be expressed as composite
optimization models. Accelerated first-order methods such as the fast
iterative shrinkage-thresholding algorithm (FISTA) have proven effective for
numerous large composite models. In this paper, we present a new variation of
FISTA, to be called C-FISTA, which obtains global linear convergence for a
broader class of composite models than many of the latest FISTA variants. We
demonstrate the versatility and effectiveness of C-FISTA by showing C-FISTA
outperforms current first-order solvers on both group Lasso and group logistic
regression models. Furthermore, we utilize Fenchel duality to prove C-FISTA
provides global linear convergence for a large class of convex models without
the loss of global linear convergence.

Keywords: composite optimization, accelerated first-order algorithm, Fenchel
duality, group Lasso.

##  1  Introduction

The following composite optimization model,

$\begin{array}[]{lll}&\min&H(B(x))+R(x)\\\ &\mbox{s.t.}&x\in{\cal
X}\subseteq\mathbb{R}^{n},\end{array}$ (1)

where $\mathcal{X}$ is a closed, convex subset of $\mathbb{R}^{n}$,
$H,R:\mathcal{X}\rightarrow\mathbb{R}$ are convex functions and
$B:\mathbb{R}^{m}\rightarrow\mathbb{R}^{n}$ is a smooth mapping, has received
ample attention in the literature [4, 5, 12, 13, 26, 35]. A special case of
(1) of immense importance is the additive composite optimization model,

$\begin{array}[]{lll}&\min&H(x)+R(x)\\\ &\mbox{s.t.}&x\in{\cal
X}\subseteq\mathbb{R}^{n},\end{array}$ (2)

which envelopes a plethora of models including: compressive sensing [8], Lasso
[34], and group Lasso and group logistic regression models [28, 36] as well as
several machine learning constructs such as support vector machines. Due to
their dimensionality, many large-scale optimization models have rendered
second-order methods computationally impractical; thus, efficient and
accelerated first-order algorithms have become essential for tackling numerous
problems. With the advent of Nesterov’s seminal work [25] much effort has been
expended toward the development of accelerated first-order methods. For a
comprehensive overview of this line of research see the monograph of
d’Aspremont et al. and the references therein [1]. As an example, utilizing
Nesterov acceleration, Beck and Teboulle developed the influential fast
iterative shrinkage-thresholding algorithm (FISTA) [4] which obtained the
optimal sublinear convergence rate $\mathcal{O}(1/k^{2})$ for (2) with $H$ and
$R$ convex and $R$ possibly nonnsmooth. In recent years, under the assumption
$H$ and potentially $R$ are strongly convex, many accelerated versions of
FISTA for (2) have been developed which obtain the optimal linear convergence
rate proved by Nesterov [7, 8, 15, 16, 31].

Furthermore, recent work has been done to determine conditions which lighten
the assumption of strong convexity while maintaining accelerated and at times
linear convergence [3, 12, 24]. In [12] the authors demonstrate local linear
convergence to a first-order stationary point for (1) without assuming any
strong convexity on $H$. The authors extended the work in [40] by utilizing a
specific error bound condition (Definition 3.1 in [12]) in tandem with an
assumption on the quadratic growth of the objective function.

In this paper we develop an accelerated composite version of FISTA, C-FISTA,
similar to the recent FISTA variants which extends to solving model (1) as
well as model (2). Also, in line with the recent works [12, 13], we forgo any
assumption guaranteeing strong convexity of the objective function in (1) and
prove C-FISTA obtains global linear convergence under certain error bound
conditions.

###  1.1 Main Contributions

This paper presents three main contributions. The first contribution of this
paper is the development of an accelerated version of FISTA which obtains
global linear convergence for (1) under specific error bound conditions. Our
algorithm, C-FISTA, differs from the aforementioned FISTA variants because
they were either exclusively designed for the additive composite model (2) or
only obtained local instead of global linear convergence for the more
generalized model (1). We demonstrate that C-FISTA is a generalization of the
FISTA variant GFISTA [7, 8] when $B(x)=x$ proving C-FISTA is an extension of
recent algorithmic developments. Further, we present an alternative
convergence analysis than those presented for the other accelerated FISTA
variants and exclude any reference to the strong convexity parameter of $R$ as
done in [7, 15, 16].

The second contribution of this paper is our utilization of Fenchel duality to
develop dual models which can be solved via C-FISTA, and our leveraging of
Fenchel duality theory [19, 32, 33] to enable efficient computation of the
subproblems of C-FISTA. Additionally, we outline a dual algorithmic approach
for a general convex model, which can be solved by C-FISTA which generalizes
the approaches of Han and Lou [19] and Auslender [2] among others [17, 18].
Sections 4 and 5 provide the relevant background on Fenchel duality and
provide an example on how we utilize the theory in the implementation of
C-FISTA.

The third contribution of this paper is the presentation of globally linearly
convergent algorithms for solving various Lasso and logistic regression models
with C-FISTA. In our numerical experiments we demonstrate C-FISTA outperforms
the state-of-the-art SLEP software package [22] for solving group, sparse-
group and overlapping sparse-group Lasso models and sparse-group logistic
regression models. We also compare C-FISTA with ADMM on the Lasso models,
which has proven linear convergence for various Lasso formulations [20], and
demonstrate superior convergence after tuning ADMM in our numerical
experiments.

###  1.2 Organization of the Paper

In Section 2 we present C-FISTA and prove the global linear convergence of the
algorithm under certain error bound conditions. Section 3 motivates C-FISTA by
describing the Lasso and logistic regression models which are solvable with
C-FISTA in the special case $B(x)=x$ in (1). In Section 4 we present the
necessary background on Fenchel duality and make connections to our
accelerated FISTA algorithm, and Section 5 provides an example of how Fenchel
duality informs our application of C-FISTA to solve the crucial proximal
mapping step. Section 6.3 contains the algorithms for solving the Lasso and
sparse-group logistic regression models with C-FISTA, and our numerical
experiments comparing C-FISTA to other efficient first-order algorithms. The
paper concludes in Section 7 with some final remarks and potential avenues of
future investigation.

##  2 A Generalized Composite Optimization Algorithm

In this paper, our primary focus revolves around the composite optimization
model,

$\begin{array}[]{lll}&\min&H(B(x))+R(x)\\\ &\mbox{s.t.}&x\in{\cal
X}\subseteq\mathbb{R}^{n},\end{array}$

where $\mathcal{X}$ is a closed, convex subset of $\mathbb{R}^{n}$,
$H:\mathbb{R}^{m}\rightarrow\mathbb{R}$ is strongly convex with parameter
$\mu>0$ and has Lipschitz gradient with parameter $L>0$,
$R:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is convex but potentially non-smooth,
and $B:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$ is of the form
$B(x)=\left(B_{1}(x),\cdots,B_{m}(x)\right)^{\top}$ with Jacobian, $J_{B}(x)$,
at $x\in\mathbb{R}^{n}$. For convenience in our analysis we denote
$F(x):=H(B(x))+R(x)$. A key assumption for feasible implementation is that the
following proximal mapping with respect to $R$ is efficiently computable:

$\mbox{\rm Prox}_{tR}(u):=\mbox{arg}\min_{x\in{\cal
X}}\frac{1}{2}\|x-u\|^{2}+tR(x),$

where $t>0$. We now state our main assumptions and error bound conditions for
our analysis:

  * (A1)

$B_{i}\in C^{2}(\mathbb{R}^{n})$ is convex, Lipschitz continuous with
parameter $\bar{L}_{i}\geq 0$, and has Lipschitz gradient constant $L_{i}\geq
0$ for $i=1,\ldots,m$.

  * (A2)

There exists $\tau>0$ such that $\tau\|x-y\|^{2}\leq\|B(x)-B(y)\|^{2}$ for all
$x,y\in\text{span}(\mathcal{X})$.

  * (A3)

There exists $\xi\geq 0$ such that
$\forall\;x^{\prime},x^{\prime\prime}\in\text{span}(\mathcal{X})$ with
$\bar{y}=B(x^{\prime})$,

$\bigg{|}\nabla_{\bar{y}}H(\bar{y})^{\top}\left(J_{B}(x^{\prime\prime})-J_{B}(x^{\prime})\right)(x^{\prime\prime}-x^{\prime})\bigg{|}\leq\xi\|x^{\prime\prime}-x^{\prime}\|^{2}.$

We note that under assumptions (A1) - (A3) it is possible $H(B(x))$ is not
strongly convex in $x$. Letting $\bar{x}:=B(x)$, we do see, however, that the
descent inequality and the strongly convex inequality still hold for $H$ with
respect to $\bar{x}$ and $\bar{y}$, i.e. for all
$\bar{x},\bar{y}\in\mathbb{R}^{m}$,

$\displaystyle H(\bar{y})\leq
H(\bar{x})+\nabla_{\bar{x}}H(\bar{x})^{\top}(\bar{y}-\bar{x})+\frac{L}{2}\|\bar{y}-\bar{x}\|^{2},$
$\displaystyle H(\bar{y})\geq
H(\bar{x})+\nabla_{\bar{x}}H(\bar{x})^{\top}(\bar{y}-\bar{x})+\frac{\mu}{2}\|\bar{y}-\bar{x}\|^{2},$

where $\nabla_{\bar{x}}H(\bar{x}):=\nabla H(y)\big{|}_{y=B(x)}$ and
$\nabla_{x}H(\bar{x})=\nabla_{x}H(B(x))=J_{B}(x)^{\top}\nabla_{\bar{x}}H(\bar{x})$.

######  Remark 2.1.

From assumption (A1) we obtain a few important inequalities. First, utilizing
the fact each $B_{i}$ is Lipschitz continuous, it follows from the limit
definition of the directional derivative that $\|\nabla
B_{i}(x)\|^{2}\leq\bar{L}_{i}^{2}$ for all $x\in\mathbb{R}^{n}$. Therefore,
letting $x\in\mathbb{R}^{n}$ be fixed we see that,

$\displaystyle\|J_{B}(x)\|$
$\displaystyle=\sup_{v\in\mathbb{R}^{n}}\left\\{\|J_{B}(x)v\|\;\bigg{|}\;\|v\|=1\right\\}$
$\displaystyle=\sup_{v\in\mathbb{R}^{n}}\left\\{\sqrt{\sum_{i=1}^{m}\left(\nabla
B_{i}(x)^{\top}v\right)^{2}}\;\bigg{|}\;\|v\|=1\right\\}$
$\displaystyle\leq\sqrt{\bar{L}_{1}^{2}+\ldots+\bar{L}_{m}^{2}}.$

So, defining $r:=\bar{L}_{1}^{2}+\cdots+\bar{L}_{m}^{2}$ we have,

$\|J_{B}(x)\|^{2}\leq r,\;\;\forall x\in\mathbb{R}^{n}.$ (3)

######  Remark 2.2.

As a result of (3) we obtain a second important inequality for our analysis.
Letting $x,y\in\mathbb{R}^{n}$, by the differentiability of $B$,

$B(y)-B(x)=\int_{0}^{1}J_{B}\left(x+t(y-x)\right)(y-x)dt,$

which by the triangle inequality and (3) implies,

$\|B(y)-B(x)\|\leq\int_{0}^{1}\|J_{B}\left(x+t(y-x)\right)\|\cdot\|y-x\|dt\leq\sqrt{r}\|y-x\|,$

and so,

$\|\bar{y}-\bar{x}\|^{2}=\|B(y)-B(x)\|^{2}\leq
r\|y-x\|^{2},\;\forall\;x,y\in\mathbb{R}^{n}.$ (4)

In order to maintain the well-posedness of (1) it behooves us to maintain the
convexity of $H(B(x))$. This along with assumption (A3) leads us to make a
distinction between two separate cases.

Case 1. $B(x)$ is affine linear, i.e. $B(x)=Ax+b$ for $A\in\mathbb{R}^{m\times
n}$ and $b\in\mathbb{R}^{m}$.

In this case, the convexity of $H(B(x))$ follows readily without the need for
any additional assumptions on $H$. Furthermore, (A1) and (A3) are satisfied
trivially with $\xi=0$ since $J_{B}(x)=A$ for all $x$, and (A2) will hold
provided $\text{null}(A)\cap\text{span}(\mathcal{X})=\\{0\\}$. In this case,
$r$ in (4) is the largest eigenvalue of $A^{\top}A$.

Case 2. $B(x)$ is not affine linear.

To maintain the convexity of the objective in (1) in this situation, we make
the added assumption $H$ is component-wise non-decreasing so that $H\circ B$
is convex.

We now present our accelerated composite FISTA algorithm, C-FISTA, for solving
(1).

Algorithm 1: C-FISTA

Step 0

Choose any $(x^{0},z^{0})\in{\cal X}\times{\cal X}$. Let $k:=0$.

Step 1

Let

$y^{k}:=\frac{1}{1+\theta}\,x^{k}+\frac{\theta}{1+\theta}\,z^{k},$ (5)

and

$\displaystyle x^{k+1}$ $\displaystyle:=$
$\displaystyle\mbox{arg}\min_{x\in{\cal
X}}H(\bar{y}^{k})+\nabla_{y}H(\bar{y}^{k})^{\top}(x-y^{k})+\frac{rL}{2}\|x-y^{k}\|^{2}+R(x)$
$\displaystyle=$ $\displaystyle\mbox{\rm
Prox}_{R/(rL)}\left(y^{k}-\frac{1}{rL}\nabla_{y}H(\bar{y}^{k})\right),$

and

$z^{k+1}:=(1-\theta)z^{k}+\theta y^{k}+\alpha(x^{k+1}-y^{k}).$ (7)

Step 2

Let $k:=k+1$ and return to Step 1 until convergence.

######  Theorem 1.

Assuming assumptions (A1) - (A3) hold with $\tau\mu-\xi>0$ and $H\circ B$
convex, then C-FISTA has an accelerated global linear rate of convergence for
(1); that is, for each iteration, we have,

$F(x^{k+1})-F(x^{*})+\left(\frac{\tau\mu-\xi}{2}\right)\|z^{k+1}-x^{*}\|^{2}\leq\\\
\left(1-\frac{\sqrt{(\tau\mu-\xi)(rL-\xi)}}{rL}\right)\left(F(x^{k})-F(x^{*})+\left(\frac{\tau\mu-\xi}{2}\right)\|z^{k}-x^{*}\|^{2}\right).$
(8)

In particular, if we perform C-FISTA for $k$ iterations, then it holds that,

$F(x^{k})-F(x^{*})\leq\left(1-\frac{\sqrt{(\tau\mu-\xi)(rL-\xi)}}{rL}\right)^{k}\left(F(x^{0})-F(x^{*})+\left(\frac{\tau\mu-\xi}{2}\right)\|z^{0}-x^{*}\|^{2}\right).$

###### Proof.

By the Lipschtiz inequality for $H$,

$\displaystyle H(\bar{x}^{k+1})$ $\displaystyle\leq$ $\displaystyle
H(\bar{y}^{k})+\nabla_{\bar{y}}H(\bar{y}^{k})^{\top}\left(\bar{x}^{k+1}-\bar{y}^{k}\right)+\frac{L}{2}\|\bar{x}^{k+1}-\bar{y}^{k}\|^{2},$
$\displaystyle\leq$ $\displaystyle
H(\bar{y}^{k})+\nabla_{\bar{y}}H(\bar{y}^{k})^{\top}\left(B(x^{k+1})-B(y^{k})\right)+\frac{rL}{2}\|x^{k+1}-y^{k}\|^{2},$
$\displaystyle=$ $\displaystyle
H(\bar{y}^{k})+\nabla_{\bar{y}}H(\bar{y}^{k})^{\top}\int_{0}^{1}J_{B}\left(y^{k}+t\left(x^{k+1}-y^{k}\right)\right)\left(x^{k+1}-y^{k}\right)dt+\frac{rL}{2}\|x^{k+1}-y^{k}\|^{2},$
$\displaystyle=$ $\displaystyle
H(\bar{y}^{k})+\nabla_{\bar{y}}H(\bar{y}^{k})^{\top}\int_{0}^{1}\left[J_{B}\left(y^{k}+t\left(x^{k+1}-y^{k}\right)\right)-J_{B}\left(y^{k}\right)\right]\left(x^{k+1}-y^{k}\right)dt$
$\displaystyle+\nabla_{\bar{y}}H(\bar{y}^{k})^{\top}J_{B}\left(y^{k}\right)\left(x^{k+1}-y^{k}\right)+\frac{rL}{2}\|x^{k+1}-y^{k}\|^{2},$

where the second line follows from (4) and the last lines are a result of the
Newton-Leibniz formula. Applying assumption (A3) we obtain,

$\displaystyle H(\bar{x}^{k+1})$ $\displaystyle\leq$ $\displaystyle
H(\bar{y}^{k})+\int_{0}^{1}\xi
t\|x^{k+1}-y^{k}\|^{2}dt+\nabla_{y}H(\bar{y}^{k})^{\top}\left(x^{k+1}-y^{k}\right)+\frac{rL}{2}\|x^{k+1}-y^{k}\|^{2},$
$\displaystyle=$ $\displaystyle
H(\bar{y}^{k})+\nabla_{y}H(\bar{y}^{k})^{\top}\left(x^{k+1}-y^{k}\right)+\frac{rL+\xi}{2}\|x^{k+1}-y^{k}\|^{2}.$

Utilizing the strong convexity of $H$, for all $x\in{\cal X}$ we have,

$\displaystyle H(\bar{x}^{k+1})$ $\displaystyle\leq$ $\displaystyle
H(\bar{x})-\nabla_{\bar{y}}H(\bar{y}^{k})^{\top}\left(\bar{x}-\bar{y}^{k}\right)-\frac{\mu}{2}\|\bar{x}-\bar{y}^{k}\|^{2}+\nabla_{y}H(\bar{y}^{k})^{\top}\left(x^{k+1}-y^{k}\right)+\frac{rL+\xi}{2}\|x^{k+1}-y^{k}\|^{2},$
(9) $\displaystyle\leq$ $\displaystyle
H(\bar{x})-\nabla_{\bar{y}}H(\bar{y}^{k})^{\top}\left(\bar{x}-\bar{y}^{k}\right)+\nabla_{y}H(\bar{y}^{k})^{\top}\left(x^{k+1}-y^{k}\right)+\frac{rL+\xi}{2}\|x^{k+1}-y^{k}\|^{2}-\frac{\tau\mu}{2}\|x-y^{k}\|^{2},$
$\displaystyle\leq$ $\displaystyle
H(\bar{x})+\frac{\xi}{2}\|x-y^{k}\|^{2}-\nabla_{y}H(\bar{y}^{k})^{\top}\left(x-y^{k}\right)+\nabla_{y}H(\bar{y}^{k})^{\top}\left(x^{k+1}-y^{k}\right)$
$\displaystyle+\frac{rL+\xi}{2}\|x^{k+1}-y^{k}\|^{2}-\frac{\tau\mu}{2}\|x-y^{k}\|^{2}$
$\displaystyle=$ $\displaystyle
H(\bar{x})+\nabla_{y}H(\bar{y}^{k})^{\top}\left(x^{k+1}-x\right)+\frac{rL+\xi}{2}\|x^{k+1}-y^{k}\|^{2}-\frac{\tau\mu-\xi}{2}\|x-y^{k}\|^{2},$

where the second inequality follows from (A2) and the last results from (A3)
as before. By the first order optimality conditions of (Step 1) we have,

$\left(\nabla_{y}H(\bar{y}^{k})+rL\left(x^{k+1}-y^{k}\right)+R^{\prime}(x^{k+1})\right)^{\top}\left(x-x^{k+1}\right)\geq
0,\,\,\mbox{ for all $x\in{\cal X}$},$

where $R^{\prime}(x^{k+1})$ is an element of the subgradient of $R$ at
$x^{k+1}$. The optimality conditions above along with the convexity of $R$ and
(9) give us,

$\displaystyle F(x^{k+1})$ $\displaystyle\leq$ $\displaystyle
F(x)+rL(x^{k+1}-y^{k})^{\top}(x-x^{k+1})$ (10) $\displaystyle\hskip
108.405pt+\frac{rL+\xi}{2}\|x^{k+1}-y^{k}\|^{2}-\frac{\tau\mu-\xi}{2}\|x-y^{k}\|^{2},\,\forall
x\in{\cal X}.$

Let $\theta\in(0,1)$, whose exact value is to be determined later. Take
$x=x^{k}$ in (10) and multiply by $1-\theta$ on both sides of the expression.
Then, let $x=x^{*}$ in (10) and multiply by $\theta$ on both sides. Finally,
adding up these two resultant inequalities and applying the assumption
$\tau\mu-\xi>0$ we have,

$\displaystyle F(x^{k+1})-F(x^{*})$ $\displaystyle\leq$
$\displaystyle(1-\theta)\left(F(x^{k})-F(x^{*})\right)-\theta\left(\frac{\tau\mu-\xi}{2}\right)\|x^{*}-y^{k}\|^{2}$
(11) $\displaystyle+rL(x^{k+1}-y^{k})^{\top}\left[(1-\theta)x^{k}+\theta
x^{*}-x^{k+1}\right]+\left(\frac{rL+\xi}{2}\right)\|x^{k+1}-y^{k}\|^{2}.\;\;\;\;\;\;\;\;$

Recall from (7) that,

$z^{k+1}=(1-\theta)z^{k}+\theta y^{k}+\alpha(x^{k+1}-y^{k}),$

where the value of the parameter $\alpha$ will be determined later. Thus,

$\displaystyle\|z^{k+1}-x^{*}\|^{2}$ $\displaystyle\leq$
$\displaystyle(1-\theta)\|z^{k}-x^{*}\|^{2}+\theta\|y^{k}-x^{*}\|^{2}+\alpha^{2}\|x^{k+1}-y^{k}\|^{2}$
(12) $\displaystyle+2\alpha(x^{k+1}-y^{k})^{\top}((1-\theta)z^{k}+\theta
y^{k}-x^{*}).$

Take $C>0$, whose value will be determined in a moment. Multiplying by $C$ on
both sides of (12), and adding the resulting inequality to (11) we obtain,

$\displaystyle F(x^{k+1})-F(x^{*})+C\|z^{k+1}-x^{*}\|^{2}$ $\displaystyle\leq$
$\displaystyle(1-\theta)\left(F(x^{k})-F(x^{*})+C\|z^{k}-x^{*}\|^{2}\right)$
$\displaystyle+\left(\theta
C-\frac{(\tau\mu-\xi)\theta}{2}\right)\|x^{*}-y^{k}\|^{2}+\left(\frac{rL+\xi}{2}+\alpha^{2}C\right)\|x^{k+1}-y^{k}\|^{2}$
$\displaystyle+rL(x^{k+1}-y^{k})^{\top}\left[(1-\theta)x^{k}+\theta
x^{*}-x^{k+1}+\frac{2\alpha C}{rL}\left((1-\theta)z^{k}+\theta
y^{k}-x^{*}\right)\right],$ $\displaystyle=$
$\displaystyle(1-\theta)\left(F(x^{k})-F(x^{*})+C\|z^{k}-x^{*}\|^{2}\right)+rL\|x^{k+1}-y^{k}\|^{2}$
$\displaystyle+rL(x^{k+1}-y^{k})^{\top}\left[(1-\theta)x^{k}-x^{k+1}+\theta\left((1-\theta)z^{k}+\theta
y^{k}\right)\right],$ $\displaystyle=$
$\displaystyle(1-\theta)\left(F(x^{k})-F(x^{*})+C\|z^{k}-x^{*}\|^{2}\right),$

if we choose the parameters to take the values

$C:=\frac{\tau\mu-\xi}{2},\,\,\alpha:=\sqrt{\frac{rL-\xi}{\tau\mu-\xi}},\,\,\theta:=\frac{\sqrt{(\tau\mu-\xi)(rL-\xi)}}{rL}.$

Observe $\theta\in(0,1)$. Since $\tau\leq r$, $\mu\leq L$ and $\tau\mu-\xi>0$,
it follows $0<\tau\mu-\xi\leq rL-\xi<rL$. Note in the last step above, we used
(5),

$y^{k}=\frac{1}{1+\theta}\,x^{k}+\frac{\theta}{1+\theta}\,z^{k}.$

Therefore, $(1-\theta)x^{k}+\theta\left((1-\theta)z^{k}+\theta
y^{k}\right)=y^{k}$, and so

$rL(x^{k+1}-y^{k})^{\top}\left[(1-\theta)x^{k}-x^{k+1}+\theta\left((1-\theta)z^{k}+\theta
y^{k}\right)\right]=-rL\|x^{k+1}-y^{k}\|^{2}.$

Hence, we have for all $k\geq 0$,

$F(x^{k+1})-F(x^{*})+\left(\frac{\tau\mu-\xi}{2}\right)\|z^{k+1}-x^{*}\|^{2}\leq\\\
\left(1-\frac{\sqrt{(\tau\mu-\xi)(rL-\xi)}}{rL}\right)\left(F(x^{k})-F(x^{*})+\left(\frac{\tau\mu-\xi}{2}\right)\|z^{k}-x^{*}\|^{2}\right),$
(13)

and so by induction starting from the $k$-th iteration we obtain the result. ∎

######  Remark 2.3.

If we have $B(x)=x$, then assumptions (A1) - (A3) hold trivially. Also, we
have $\xi=0$, $r=1$ and $\tau=1$ ensuring $\tau\mu-\xi>0$ which gives us the
convergence result,

$F(x^{k})-F(x^{*})\leq\left(1-\sqrt{\frac{\mu}{L}}\right)^{k}\left(F(x^{0})-F(x^{*})+\frac{\mu}{2}\|z^{0}-x^{*}\|^{2}\right).$

If $z^{0}=x^{0}$, then this is the same convergence result produced by GFISTA
[7, 8] which is an accelerated variant of FISTA for (2). Therefore, C-FISTA
generalizes GFISTA into the broader problem class given by the composite
optimization model (1).

######  Remark 2.4.

The condition $\tau\mu-\xi>0$ essentially requires that the product of the
first-order derivatives of $H$ and the second-order derivative of $B$ should
not exceed the product of their curvatures. In some cases, this condition is
easy to satisfy by variable-transformation, or scaling. For example, when
$B(x)$ is homogeneous of degree $d_{1}>0$, i.e. $B(\lambda
x)=\lambda^{d_{1}}B(x)$ for any $x$ and $\lambda>0$, and $R(x)$ is homogeneous
with degree $d_{2}>0$, i.e. $R(\lambda x)=\lambda^{d_{2}}R(x)$. Then for any
$\lambda>0$ one can scale the variable as
$x^{\prime}:=\lambda^{\frac{1}{d_{1}}}x$, and the problem is turned
equivalently into minimizing $H(B(x^{\prime}))+R(x^{\prime})=H(\lambda
B(x))+\lambda^{\frac{d_{2}}{d_{1}}}R(x)$. After the above change of variables,
the new objective has the modified parameters $\tau:=\lambda^{2}\tau$,
$\mu:=\mu$, and $\xi:=\lambda\xi$. Therefore, in this situation by
appropriately choosing $\lambda>0$ one can always satisfy the condition
$\tau\mu-\xi>0$.

##  3  Motivation: Lasso & Logistic Regression Models

A class of motivating composite optimization models for which C-FISTA is
applicable, in the special case $B(x)=x$, are the various Lasso formulations.
The standard Lasso model originating from Tibshirani [34],

$\begin{array}[]{lll}&\min&\frac{1}{2}\|Ax-b\|^{2}+\gamma\|x\|_{1}\\\
&\mbox{s.t.}&x\in\mathbb{R}^{n},\end{array}$

is solvable via C-FISTA whenever $A$ is of full-column rank. The results from
Section 2 further guarantee global linear convergence to the optimal solution
even when there exists a constraint on $x$. For example, the following
constrained Lasso model is solvable with Algorithm 1,

$\begin{array}[]{ll}\min&\frac{1}{2}\|Ax-b\|^{2}+\gamma\|x\|_{1}\\\
\mbox{s.t.}&l_{i}\leq x_{i}\leq u_{i},\,\,i=1,2,...,n.\end{array}$

The utility of the Lasso model is detailed extensively in the literature; see
[21, 28, 36, 37, 38]. The success of the basic model spawned the development
of several Lasso variants including: the group Lasso (GL) model [28, 37],

$\begin{array}[]{lll}&\min&\frac{1}{2}\|Ax-b\|^{2}+\gamma\sum_{j\in{\cal
G}}\|x(j)\|\\\ &\mbox{s.t.}&x\in\mathbb{R}^{n}\end{array}$

where $x(j)$ is a subvector of $x$ corresponding to the indices in $j\in{\cal
G}$ where ${\cal G}$ disjointly partitions $\\{1,2,\ldots,n\\}$, the sparse-
group Lasso (SGL) model [30, 38],

$\begin{array}[]{lll}&\min&\frac{1}{2}\|Ax-b\|^{2}+\gamma_{1}\sum_{j\in{\cal
G}}\|x(j)\|+\gamma_{2}\|x\|_{1}\\\
&\mbox{s.t.}&x\in\mathbb{R}^{n},\end{array}$

and the overlapping group and overlapping sparse-group Lasso (OSGL) models
[21, 36, 38] which are of identical form to the group Lasso and sparse-group
Lasso models respectively expect ${\cal G}$ no longer is a disjoint partition
of $\\{1,\ldots,n\\}$. As we will see in Section 5, the subproblems required
in C-FISTA for the Lasso models are solvable in closed form or by an efficient
subroutine in the case of the overlapping group/sparse-group Lasso models with
global linear convergence achieved in each variant.

Another model where C-FISTA provides global linear convergence is the sparse-
group logistic regression (SGLR) model [22, 23],

$\begin{array}[]{lll}&\min&\frac{1}{m}\sum_{i=1}^{m}\ln\left(1+\exp\left(-y_{i}\left(a_{i}^{\top}x+b\right)\right)\right)+\gamma_{1}\sum_{j\in{\cal
G}}\|x(j)\|+\gamma_{2}\|x\|_{1}\\\
&\mbox{s.t.}&x\in\mathbb{R}^{n},b\in\mathbb{R},\end{array}$

where $y\in\mathbb{R}^{n}$ is a vector containing only entries of $\pm 1$. The
sparse-group logistic regression model has applications in various machine
learning models especially in the area of classification, and the subproblems
for the logistic regression model are in the same format as the sparse-group
Lasso model detailed extensively in Section 5. Before discussing how to
compute subproblems for the Lasso models, we first discuss Fenchel duality in
Section 4 and describe how C-FISTA can be utilized via a dual approach to
solve a general constrained composite model.

##  4  A Dual Formulation and Algorithmic Approach

In this section, we outline the construction of a dual algorithm using Fenchel
duality which generates an approximate primal solution from the dual solution
at every iterate. To elucidate this approach we first provide a brief overview
of standard convex analysis results regarding conjugate functions and Fenchel
duality [32, 33]. Unless otherwise stated, we use the standard notation of
Rockafellar [32]. For generality, we will consider the following constrained
composite model:

$\displaystyle\min$ $\displaystyle h_{0}(x)+h_{1}(x)+\cdots+h_{\ell}(x)$ s.t.
$\displaystyle x\in{\cal C}_{j}\subseteq\mathbb{R}^{n},\,j=1,2,...,r,$

where the $h_{i}$’s are proper convex functions on $\mathbb{R}^{n}$ and the
${\cal C}_{j}$’s are closed convex subsets of $\mathbb{R}^{n}$. To begin, let
us introduce some basic notions in convex analysis.

Let $S\subseteq\mathbb{R}^{n}$ be a closed convex set. The support function of
$S$ is defined as,

$f_{S}(y):=\max_{x\in S}y^{\top}x,$

and the polar set of $S$ is,

$S^{\circ}:=\\{y\mid y^{\top}x\leq 1,\,\forall x\in S\\}.$

Thus, we can succinctly write, $S^{\circ}=\\{y\mid f_{S}(y)\leq 1\\}$. Let $g$
be a convex function whose domain is contained in $S$. The conjugate of $g$ is
defined as,

$g^{*}(y):=\sup_{x\in S}\left[y^{\top}x-g(x)\right],$

where the function is assumed to take value $+\infty$ anywhere outside its
domain. In deriving the Fenchel dual of (4) we need some standard results from
convex analysis.

######  Lemma 2.

Suppose $f$ is a proper convex function, then

$-f^{*}(0)=\inf_{x}\;f(x).$ (15)

######  Lemma 3.

(Theorem 16.4 [32]) Suppose $g_{1},...,g_{r}$ are proper convex functions on
$\mathbb{R}^{n}$ and the intersection of the relative interiors of the domains
of the $g_{i}$’s are non-empty, then

$\left(g_{1}+\cdots+g_{r}\right)^{*}(x)=\inf_{y_{1}+\cdots+y_{r}=x}\sum_{j=1}^{r}g_{j}^{*}(y_{j}).$
(16)

######  Lemma 4.

For a convex set $S$, suppose its indicator function is defined as,

$I_{S}(x):=\left\\{\begin{array}[]{ll}0,&\mbox{if $x\in S$};\\\
+\infty,&\mbox{if $x\not\in S$}.\end{array}\right.$

Then, the conjugate of the indicator function is simply the support function
of $S$; that is,

$I_{S}^{*}=f_{S}.$

The proofs of Lemmas 2 and 4 follow directly from the definitions of conjugate
and support functions. Thus, rewriting (4) in the equivalent form,

$\displaystyle\min$ $\displaystyle
h_{0}(x)+h_{1}(x)+\cdots+h_{\ell}(x)+I_{{\cal C}_{1}}(x)+\ldots+I_{{\cal
C}_{r}}(x)$ s.t. $\displaystyle x\in\mathbb{R}^{n},$

and utilizing Lemmas 2, 3, and 4, the Fenchel dual of (4) is:

$\displaystyle\min$ $\displaystyle h_{0}^{*}(-y_{1}-\cdots-
y_{\ell}-y_{\ell+1}-\cdots-y_{\ell+r})$
$\displaystyle+h_{1}^{*}(y_{1})+\cdots+h^{*}_{\ell}(y_{\ell})+f_{{\cal
C}_{1}}(y_{\ell+1})+\cdots+f_{{\cal C}_{r}}(y_{\ell+r})$ s.t. $\displaystyle
y_{j}\in\mathbb{R}^{n},\,j=1,2,...,\ell+r.$

Note, for the sake of presentation we have written the dual problem (4) in
such a manner that the optimal solutions of (4) and (4) differ by a negative
sign. To better understand the link between the primal and dual models, we
demonstrate how the optimality conditions of (4) and (4) are linked. Assuming
the conjugate functions of the $h_{i}$’s ($i=0,1,...,\ell$) are differentiable
and denoting $y^{*}_{j}$ for $j=1,\ldots,\ell+r$ to be optimal for (4), we
have the first-order optimality conditions,

$\left\\{\begin{array}[]{ll}-\nabla
h_{0}^{*}(-\sum_{j=1}^{\ell+r}y_{j}^{*})+\nabla
h_{i}^{*}(y^{*}_{i})=0,&i=1,...,\ell\\\ \nabla
h_{0}^{*}(-\sum_{j=1}^{\ell+r}y_{j}^{*})\in\partial f_{{\cal
C}_{j}}(y^{*}_{\ell+j}),&j=1,...,r.\end{array}\right.$ (18)

Denoting $x^{*}:=\nabla h_{0}^{*}(-\sum_{j=1}^{\ell+r}y_{j}^{*})$, the
optimality conditions imply $x^{*}=\nabla h_{i}^{*}(y^{*}_{i})$ for
$i=1,\ldots,\ell$, and so by Lemma 5 of [39], $\nabla h_{i}(x^{*})=y^{*}_{i}$,
for $i=1,...,\ell$. Similarly, because $f_{{\cal C}_{j}}^{*}=I_{{\cal
C}_{j}}$, if $x^{*}\in\partial f_{{\cal C}_{j}}(y^{*}_{\ell+j})$, then
$y^{*}_{\ell+j}\in\partial I_{{\cal C}_{j}}(x^{*})$, which means that
$y^{*}_{\ell+j}$ is a normal direction at $x^{*}$:

$(y^{*}_{\ell+j})^{\top}(x-x^{*})\leq 0,\mbox{ for all $x\in{\cal C}_{j}$}.$

Furthermore, by (18)

$\nabla h_{0}(x^{*})=-\sum_{j=1}^{\ell+r}y_{j}^{*}=-\sum_{i=1}^{\ell}\nabla
h_{i}(x^{*})-\sum_{j=1}^{r}y^{*}_{\ell+j}.$ (19)

Thus, the first-order optimality conditions for the dual problem implies (19),
which is the optimality condition for the primal problem (4). To see this,
note that

$\nabla h_{0}(x^{*})+\sum_{i=1}^{\ell}\nabla
h_{i}(x^{*})=-\sum_{j=1}^{r}y^{*}_{\ell+j},$

and so,

$\left(\nabla h_{0}(x^{*})+\sum_{i=1}^{\ell}\nabla
h_{i}(x^{*})\right)^{\top}(x-x^{*})=-\sum_{j=1}^{r}\left(y^{*}_{\ell+j}\right)^{\top}(x-x^{*})\geq
0,\,\,\forall x\in\cap_{j=1}^{r}{\cal C}_{j}.$

We know present an important relationship between the primal solution induced
by the dual iterations which outlines a dual approach to solving (4).

######  Proposition 5.

Suppose that $h_{0}$ is a strongly convex function with strong convexity
parameter $\mu$ and gradient Lipschitz constant $L$. Let $\\{y^{k}\\}$ be a
sequence converging to the dual solution $y^{*}$. For each $y^{k}$ in the dual
sequence, we recover a primal solution $x^{k}:=\nabla
h_{0}^{*}(-\sum_{j=1}^{\ell+r}y_{j}^{k})$ such that,

$\frac{1}{L}\left\|\sum_{j=1}^{\ell+r}\left(y^{k}_{j}-y^{*}_{j}\right)\right\|\leq\|x^{k}-x^{*}\|\leq\frac{1}{\mu}\left\|\sum_{j=1}^{\ell+r}\left(y^{k}_{j}-y^{*}_{j}\right)\right\|.$

###### Proof.

By the Fenchel duality relation, we know that $h_{0}^{*}$ is also a strongly
convex function with strong convexity parameter $1/L$ and gradient Lipschitz
constant $1/\mu$ (Theorem 1 [39]). Further, by the definition of $x^{k}$ and
previously stated results, $\nabla
h_{0}(x^{k})=-\sum_{j=1}^{\ell+r}y_{j}^{k}$. Therefore, $\nabla
h_{0}(x^{k})-\nabla h_{0}(x^{*})=\sum_{j=1}^{\ell+r}(y_{j}^{*}-y_{j}^{k})$,
and so the gradient Lipschitz condition on $h$ yields,

$\left\|\sum_{j=1}^{\ell+r}\left(y^{k}_{j}-y^{*}_{j}\right)\right\|=\|\nabla
h_{0}(x^{k})-\nabla h_{0}(x^{*})\|\leq L\|x^{k}-x^{*}\|.$

On the other hand, the gradient Lipschitz condition on $h^{*}$ states,

$\|x^{k}-x^{*}\|=\left\|\nabla
h_{0}^{*}\left(-\sum_{j=1}^{\ell+r}y_{j}^{k}\right)-\nabla
h_{0}^{*}\left(-\sum_{j=1}^{\ell+r}y_{j}^{*}\right)\right\|\leq\frac{1}{\mu}\left\|\sum_{j=1}^{\ell+r}\left(y^{k}_{j}-y^{*}_{j}\right)\right\|.$

∎

Therefore, the above result gives us the framework to develop a general dual
approach to solving (4). Applying any algorithm to solve the dual problem (4)
we can recover a primal solution at each iteration; furthermore, the rate of
convergence on the primal side for any algorithm will only differ by a fixed
constant, $1/\mu$, in-comparison to the rate of convergence on the dual side.

The presented dual formulation is not novel; Han and Lou [19] and Fukushima et
al. [17] have very similar derivations of the dual and many algorithms such as
Dysktra’s projection algorithm [18] utilize Fenchel duality. The clarity
provided by the previous result is that, under the assumption $h_{0}$ is
strongly convex with gradient Lipschitz, any algorithm which generates a
converging dual solution will generate a primal solution with the same rate of
convergence up to a constant factor. Connecting this result back to C-FISTA in
the special case $B(x)=x$, we can solve the dual problem (4) and obtain a
primal solution from Proposition 5. Therefore, any primal problem (4) can be
solved with C-FISTA assuming the required proximal mappings are computable and
the necessary assumptions for Theorem 1 are maintained.

Besides using Fenchel duality in a direct fashion as detailed in this section,
one can also utilize Fenchel duality theory to solve intermediate subproblems
for the primal. In the next section we describe this process using the sparse-
group Lasso model as an example.

##  5 Solving Subproblems with Fenchel Duality

One example of how duality can be leveraged to apply C-FISTA can be seen in
how we solve the subproblems for the group and sparse-group Lasso models. In
the sparse-group Lasso model, the main subproblem which needs to be solved to
apply C-FISTA is of the form,

$\begin{array}[]{lll}(S)&\min&\frac{1}{2}\|x-d\|^{2}+\gamma_{1}\|x\|+\gamma_{2}\|x\|_{1}\\\
&\mbox{s.t.}&x\in\mathbb{R}^{n},\end{array}$

with $\gamma_{1},\gamma_{2}\geq 0$. Utilizing Lemmas 2 and 3 and standard
conjugate functions [33], we see the Fenchel dual of $(S)$ is,

$\begin{array}[]{lll}(DS)&\min&\frac{1}{2}\|d+y_{1}+y_{2}\|^{2}\\\
&\mbox{s.t.}&\|y_{1}\|\leq\gamma_{1},\\\
&&\|y_{2}\|_{\infty}\leq\gamma_{2}.\end{array}$

When $y_{2}$ is fixed, consider the simple projection problem
$\min_{\|y_{1}\|\leq\gamma_{1}}\,\,\|d+y_{2}+y_{1}\|$, whose solution is
explicit:

$y_{1}^{*}=\left\\{\begin{array}[]{ll}-d-y_{2},&\mbox{ if
$\|d+y_{2}\|\leq\gamma_{1}$};\\\ -\gamma_{1}(d+y_{2})/\|d+y_{2}\|,&\mbox{ if
$\|d+y_{2}\|>\gamma_{1}$}.\end{array}\right.$

Therefore,

$\|d+y_{2}+y_{1}^{*}\|=\left\\{\begin{array}[]{ll}0,&\mbox{ if
$\|d+y_{2}\|\leq\gamma_{1}$};\\\ \|d+y_{2}\|-\gamma_{1},&\mbox{ if
$\|d+y_{2}\|>\gamma_{1}$}.\end{array}\right.$

Thus, the dual problem $(DS)$ is equivalent to
$\min_{\|y_{2}\|_{\infty}\leq\gamma_{2}}\,\,(\|d+y_{2}\|-\gamma_{1})_{+}$,
which can be solved through the associated model,

$\begin{array}[]{lll}(DS)^{\prime}&\min&\frac{1}{2}\|d+y_{2}\|^{2}\\\
&\mbox{s.t.}&\|y_{2}\|_{\infty}\leq\gamma_{2}.\end{array}$

Since the individual components of $y_{2}$ are decoupled in $(DS)^{\prime}$,
the above problem can be reduced to solving $n$ 1-dimensional problems,

$\begin{array}[]{ll}\min&(d_{i}+t)^{2}\\\
\mbox{s.t.}&|t|\leq\gamma_{2},\end{array}$

for $i=1,2,...,n$. Observe that the solutions for the above 1-dimensional
models can be found using the following thresholding operator:

${\rm Th}_{\gamma_{2}}(t):=\left\\{\begin{array}[]{rl}-\gamma_{2},&\mbox{if
$t<-\gamma_{2}$};\\\ t,&\mbox{if $-\gamma_{2}\leq t\leq\gamma_{2}$};\\\
\gamma_{2},&\mbox{if $t>\gamma_{2}$}.\\\ \end{array}\right.$

Therefore, a solution for $(DS)$ is $y_{2}^{*}={\rm Th}_{\gamma_{2}}(-d)$,
where $({\rm Th}_{\gamma_{2}})_{i}={\rm Th}_{\gamma_{2}}(-d_{i})$ for all $i$,
and,

$y_{1}^{*}=\left\\{\begin{array}[]{ll}-d-{\rm Th}_{\gamma_{2}}(-d),&\mbox{ if
$\|d+{\rm Th}_{\gamma_{2}}(-d)\|\leq\gamma_{1}$};\\\ -\gamma_{1}(d+{\rm
Th}_{\gamma_{2}}(-d))/\|d+{\rm Th}_{\gamma_{2}}(-d)\|,&\mbox{ if $\|d+{\rm
Th}_{\gamma_{2}}(-d)\|>\gamma_{1}$}.\end{array}\right.$

After solving $(DS)$ with the optimal solution $(y_{1}^{*},y_{2}^{*})$, one
recovers the optimal solution to $(S)$ as $x^{*}=d+y_{1}^{*}+y_{2}^{*}$ using
the results from Section 4. In particular,

$x^{*}=\left\\{\begin{array}[]{ll}0,&\mbox{ if $\|d+{\rm
Th}_{\gamma_{2}}(-d)\|\leq\gamma_{1}$};\\\ (\|d+{\rm
Th}_{\gamma_{2}}(-d)\|-\gamma_{1})/\|d+{\rm Th}_{\gamma_{2}}(-d)\|\left(d+{\rm
Th}_{\gamma_{2}}(-d)\right),&\mbox{ if $\|d+{\rm
Th}_{\gamma_{2}}(-d)\|>\gamma_{1}$}.\end{array}\right.$

Hence, by Fenchel duality we obtain a closed form solution to $(S)$ yielding
an exact solution to the required subproblems to apply C-FISTA to the sparse-
group Lasso model. Similar approaches can be taken to solve other potential
subproblems which arise when applying C-FISTA. In the next section, we present
the algorithms and numerical experiments for solving the models discussed in
Section 3.

##  6  Numerical Experiments

To demonstrate the practicality and efficiency of C-FISTA, we conducted
numerous experimental tests on group Lasso, sparse-group Lasso, overlapping
sparse-group Lasso and sparse-group logistic regression models. Subsection 6.1
details how to apply C-FISTA to solve the Lasso models; Subsection 6.2
presents the solution procedure for the sparse-group logistic regression
model, and Subsection 6.3 displays the results of the numerical experiments
comparing C-FISTA, the state-of-the-art software package SLEP [22], and ADMM
on the Lasso and logistic regression models.

###  6.1 Lasso Models

In this section we state how to apply C-FISTA to solve the various Lasso
models described in Section 3 beginning with the group Lasso model.

####  6.1.1 C-FISTA for Group Lasso

Applying C-FISTA to solve any model hinges on computing the proximal mapping
in (Step 1). This key subproblem in the group Lasso model is of the form,

$\begin{array}[]{ll}\min&\frac{1}{2}\|x-d\|^{2}+\gamma\|x\|\\\
\mbox{s.t.}&x\in\mathbb{R}^{n},\end{array}$ (20)

which by the first-order optimality conditions has the simple closed-form
solution,

$x^{*}=\begin{cases}0,&\text{ if }\|d\|\leq\gamma;\\\
\left(\frac{\|d\|-\gamma}{\|d\|}\right)d,&\text{ if }\|d\|>\gamma.\end{cases}$

With the solution to the group Lasso subproblems (20), we write down how to
solve the group Lasso model with C-FISTA in Algorithm 2.

Algorithm 2: C-FISTA for $(GL)$ Step 0. Choose any
$(x^{0},z^{0})\in\mathbb{R}^{n}\times\mathbb{R}^{n}$. Let $k:=0$.
Step 1.  Let,

$y^{k}:=\frac{1}{1+\theta}x^{k}+\frac{\theta}{1+\theta}z^{k}.$

Step 2. For $j=1,\ldots,J$:

$x(j)^{k+1}:=\begin{cases}0,&\|d(j)\|\leq\gamma/L;\\\
\left(\frac{\|d(j)\|-\gamma/L}{\|d(j)\|}\right)d(j),&\|d(j)\|>\gamma/L;\end{cases}$

with $d:=y^{k}-\frac{1}{L}A^{\top}(Ay^{k}-b)$ and corresponding subvectors
$d(j)$.

Step 3.

$z^{k+1}:=(1-\theta)z^{k}+\theta y^{k}+{\alpha}\left(x^{k+1}-y^{k}\right).$

Step 4. Let $k:=k+1$; return to Step 1 until convergence.

Practical application of Algorithm 2 requires the user to determine bounds for
the Lipschitz constant and strong convexity constants, $L$ and $\mu$, for
$H(x)$. For the group Lasso models $H(x):=\frac{1}{2}\|Ax-b\|^{2}$ making such
bounds readily available. By the definition of strong convexity we see
$\mu=\lambda_{\min}\left(A^{\top}A\right)$. Similarly, using the definition of
the Lipschtiz constant, we see that
$L\leq\lambda_{\max}\left(A^{\top}A\right)$ where
$\lambda_{\min}\left(A^{\top}A\right)$ and
$\lambda_{\max}\left(A^{\top}A\right)$ are the smallest and largest
eigenvalues of $A^{\top}A$ respectively. For other strongly convex functions
$H(x)$ on which C-FISTA is applicable, tight bounds for $L$ and $\mu$ might be
unavailable. In these situations conservative estimates for these bounds or a
backtracking scheme such as in [4] to estimate these parameters would be
required. In this paper we do not focus on subroutines to estimate these
bounds in our algorithms; however, the backtracking strategies applied in [4,
7, 15] could similarly be implemented in many instances of C-FISTA.

####  6.1.2 C-FISTA for Sparse-Group Lasso

The key subproblem in applying C-FISTA to solve the sparse-group Lasso model
was determined in Section 5; thus, we can write Algorithm 3 to solve the
sparse-group Lasso model with C-FISTA.

Algorithm 3: C-FISTA for $(SGL)$ Step 0. Choose any
$(x^{0},z^{0})\in\mathbb{R}^{n}\times\mathbb{R}^{n}$. Let $k:=0$.
Step 1.  Let

$y^{k}:=\frac{1}{1+\theta}x^{k}+\frac{\theta}{1+\theta}z^{k}.$

Step 2. For $j=1,\ldots,J$:

$x(j)^{k+1}:=\left\\{\begin{array}[]{ll}0,&\mbox{ if $\|d(j)+{\rm
Th}_{\gamma_{2}/L}(-d(j))\|\leq\gamma_{1}/L$};\\\ \left(\frac{\|d(j)+{\rm
Th}_{\gamma_{2}/L}(-d(j))\|-\gamma_{1}/L)}{\|d(j)+{\rm
Th}_{\gamma_{2}/L}(-d(j))\|}\right)\left(d(j)+{\rm
Th}_{\gamma_{2}/L}(-d(j))\right),&\mbox{ if $\|d(j)+{\rm
Th}_{\gamma_{2}/L}(-d(j))\|>\gamma_{1}/L$};\end{array}\right.$

with $d:=y^{k}-\frac{1}{L}\nabla H(y^{k})=y^{k}-\frac{1}{L}A^{\top}(Ay^{k}-b)$
and $d(j)$ the corresponding subvectors.

Step 3.

$z^{k+1}:=(1-\theta)z^{k}+\theta y^{k}+{\alpha}\left(x^{k+1}-y^{k}\right).$

Step 4. Let $k:=k+1$; return to Step 1 until convergence.

Since $H(x)$ is identical in all of the Lasso models, the same bounds for the
Lipschitz and strong convexity constants in Algorithm 2 apply for each of the
Lasso models. Comparing Algorithms 2 and 3 we note only Step 2 has been
updated because the $y$ and $z$-updates are independent of the objective
function. The alteration in Step 2 is solely due to the difference in the
regularization term in the group and sparse-group Lasso models which alters
the proximal mapping (Step 1). Changes in the smooth part of the objective
function, $H$, appear as changes in $d:=y^{k}-\frac{1}{L}\nabla H(y^{k})$,
recall from Remark 2.3 that $r=1$ in the special case $B(x)=x$; therefore,
altering $H$ presents few challenges to implementing C-FISTA; the key
computational hurdle lies in computing efficient proximal mappings which
depend on $R$.

####  6.1.3 C-FISTA for Overlapping Sparse-Group Lasso

The key difference between the overlapping sparse-group Lasso model and those
previously discussed is the removal of the prohibition that the groups of
variables cannot intersect. Allowing overlapping groups of variables removes
the ability to decouple the minimization in (Step 1) into subproblems over the
individual subvectors. In the sparse-group Lasso model, since the groups of
variables do not overlap, we were able to show in Section 5 there was a closed
form solution to the proximal mappings; however, with overlapping groups the
key subproblem for implementing C-FISTA,

$\min\;\frac{1}{2}\|x-d\|^{2}+\frac{\gamma_{1}}{L}\sum_{j\in{\cal
G}}\|x(j)\|+\frac{\gamma_{2}}{L}\|x\|_{1},$ (21)

is no longer decomposable and a simple closed form solution is unobtainable;
therefore, we must solve this subproblem numerically at each iteration in-
order to apply C-FISTA to the overlapping sparse-group Lasso model. In [36]
the authors apply a variation of FISTA to solve the (OSGL) model. Our approach
is distinct from their approach because we accelerate FISTA through the use of
two sequences, $\\{y^{k}\\}$ and $\\{z^{k}\\}$, while they use a single
sequence performing a Nesterov-like acceleration. In order to apply their
algorithm, the authors of [36] developed an efficient subroutine,
overlapping.c, to numerically solve (21). Their subroutine is provided in the
SLEP software package [22]. In our solving of the (OSGL) model with C-FISTA,
we apply Algorithm 3 exactly in the same fashion as for the sparse-group Lasso
model except when computing $x^{k+1}$ in Step 2 we solve (21) with
overlapping.c from the SLEP software package.

###  6.2 C-FISTA for Sparse-Group Logistic Regression

In this section we describe how to apply C-FISTA to solve the sparse-group
logistic regression model,

$\begin{array}[]{lll}&\min&\frac{1}{m}\sum_{i=1}^{m}\ln\left(1+\exp\left(-y_{i}\left(a_{i}^{\top}x+b\right)\right)\right)+\gamma_{1}\sum_{j\in{\cal
G}}\|x(j)\|+\gamma_{2}\|x\|_{1}\\\
&\mbox{s.t.}&x\in\mathbb{R}^{n},b\in\mathbb{R}.\end{array}$

The sparse-group Lasso and sparse-group logistic regression models only differ
with respect to their strongly convex term, so the procedure for solving the
sparse-group logistic regression model with C-FISTA is very similar to
Algorithm 3. Only two alterations to Algorithm 3 are required to solve the
sparse-group logistic regression model. One, we must update $d$ in Algorithm 3
to be,

$d=y^{k}-\frac{1}{L}\nabla
H(y^{k})=y^{k}+\frac{1}{mL}\sum_{i=1}^{m}\frac{y_{i}\exp(-y_{i}(a_{i}^{\top}y^{k}+b))}{1+\exp(-y_{i}(a_{i}^{\top}y^{k}+b))}\begin{pmatrix}a_{i}\\\
1\end{pmatrix},$ (22)

and second we require new estimates for the bounds of the Lipschitz and strong
convexity constants $L$ and $\mu$. Updating $d$ is straightforward, but
obtaining tight upper and lower bounds on the Lipschtiz and strong convexity
constants can become a challenging task as mentioned previously. Algorithm 4
presents C-FISTA for solving the sparse-group logistic regression model. In
our implementation of Algorithm 4, we fixed the Lipschitz and strong convexity
parameters though backtracking procedures could be utilized.

Algorithm 4: C-FISTA for $(SGLR)$ Step 0. Choose any
$(x^{0},b^{0},z^{0})\in\mathbb{R}^{n}\times\mathbb{R}\times\mathbb{R}^{n+1}$.
Let $k:=0$.
Step 1.  Let

$y^{k}:=\frac{1}{1+\theta}\begin{pmatrix}x^{k}\\\
b^{k}\end{pmatrix}+\frac{\theta}{1+\theta}z^{k}.$

Step 2. For $j\in J$:

$x(j)^{k+1}:=\left\\{\begin{array}[]{ll}0,&\mbox{ if $\|d(j)+{\rm
Th}_{\gamma_{2}/L}(-d(j))\|\leq\gamma_{1}/L$};\\\ \left(\frac{\|d(j)+{\rm
Th}_{\gamma_{2}/L}(-d(j))\|-\gamma_{1}/L)}{\|d(j)+{\rm
Th}_{\gamma_{2}/L}(-d(j))\|}\right)\left(d(j)+{\rm
Th}_{\gamma_{2}/L}(-d(j))\right),&\mbox{ if $\|d(j)+{\rm
Th}_{\gamma_{2}/L}(-d(j))\|>\gamma_{1}/L$};\end{array}\right.$

where $d$ is as defined in (22) and $d(j)$ is the subvector of $d$
corresponding to the subvector $x(j)$.
Step 3.

$b^{k+1}=y^{k+1}_{n+1}+\frac{1}{mL}\sum_{i=1}^{m}\frac{y_{i}\exp(-y_{i}(a_{i}^{\top}y^{k}+b))}{1+\exp(-y_{i}(a_{i}^{\top}y^{k}+b))}.$

Step 4.

$z^{k+1}:=(1-\theta)z^{k}+\theta y^{k}+{\alpha}\left[\begin{pmatrix}x^{k+1}\\\
b^{k+1}\end{pmatrix}-y^{k}\right].$

Step 5. Let $k:=k+1$; return to Step 1 until convergence.

Note, computing $b^{k+1}$ at each iteration is simple relative to $x^{k+1}$
because the subprobem (Step 1) is decomposable into a minimization over $x$
and $b$ respectively. Since the minimization over $b$ contains no
regularization terms, the first-order optimality conditions provide a simple
update for $b$.

###  6.3 Experimental Set-Up and Results

####  6.3.1 Experimental Set-Up for the Lasso Models

In this section we present our experimental set-up for comparing C-FISTA, SLEP
and ADMM on the Lasso models. We constructed synthetic data sets for testing
in a manner similar to the process utilized in [38]. For our tests of the
group and sparse-group Lasso models, we randomly generated a $m\times n$ data
matrix $A$ from the standard normal distribution, and we formed three
different subvector groups, $\mathcal{G}_{1},\mathcal{G}_{2},$ and
$\mathcal{G}_{3}$, where the subvectors in $\mathcal{G}_{i}$ where of sizes
$10,100,$ and 200 respectively. Thus, in our experiments we let
$A\in\mathbb{R}^{8000\times 5000}$ and had,

$\displaystyle\mathcal{G}_{1}$
$\displaystyle=\\{(1,\ldots,10),(11,\ldots,20),\ldots,(4991,\ldots,5000)\\},$
$\displaystyle\mathcal{G}_{2}$
$\displaystyle=\\{(1,\ldots,100),(101,\ldots,200),\ldots,(4901,\ldots,5000)\\},$
$\displaystyle\mathcal{G}_{3}$
$\displaystyle=\\{(1,\ldots,200),(201,\ldots,400),\ldots,(4801,\ldots,5000)\\}.$

The response vectors $b^{i}$ for $i=1,2,3$ were formed as,

$b^{i}=Ax+\delta\cdot\epsilon,$

where each component of $\epsilon$ was drawn from the standard norm
distribution, $\delta>0$ provided a scaling factor for the noise, and
$x(j)=(1,2,\ldots,10,0,\ldots,0)^{\top}$ for $j=1,2,\ldots,10$ with
$x(j)=\mathbf{0}$ for the remaining subvector groups
$j=11,\ldots,|\mathcal{G}_{i}|$. We set the positive penalty parameters
$\gamma_{1}$ and $\gamma_{2}$ to be equal and of a magnitude to ensure sparse
but non-trivial solutions.

For the overlapping sparse-group Lasso model we constructed the data matrix
$A$ in the same manner, and let the subvector groups be,

$\displaystyle\mathcal{G}_{1}$
$\displaystyle=\\{(1,\ldots,10),(6,\ldots,15),\ldots,(4986,\ldots,4995),(4991,\ldots,5000)\\},$
$\displaystyle\mathcal{G}_{2}$
$\displaystyle=\\{(1,\ldots,100),(51,\ldots,150),\ldots,(4851,\ldots,4950),(4901,\ldots,5000)\\},$
$\displaystyle\mathcal{G}_{3}$
$\displaystyle=\\{(1,\ldots,200),(101,\ldots,300),\ldots,(4701,\ldots,4900),(4801,\ldots,5000)\\}.$

Thus, each of the adjacent subvectors in $\mathcal{G}_{i}$ for $i=1,2,3$
overlapped with another subvector. The response vectors were constructed as
done for the group and sparse-group Lasso models.

####  6.3.2 Experimental Set-Up for the Logistic Regression Model

The set-up for the numerical experiements on the sparse-group logistic
regression model were similar to the set-up for the Lasso models. We generated
the data matrix $A\in\mathbb{R}^{500\times 5000}$ from the standard normal
distribution letting $a_{i}$ for $i=1,2,\ldots,500$ be the rows of $A$, and
set $y=(-\mathbf{e}^{\top},\mathbf{e}^{\top})^{\top}$ where
$\mathbf{e}\in\mathbb{R}^{250}$ is the vector of all ones. We conducted three
experiments from the generated data using the same subvector groups,
$\mathcal{G}_{1}$, $\mathcal{G}_{2}$ and $\mathcal{G}_{3}$, utilized in the
group and sparse-group Lasso tests. We set $\gamma_{1}=\gamma_{2}=0.01$ in our
tests to ensure a non-trivial solution with substantial sparsity.

###  6.4 Experimental Results

All of the numerical experiments were conducted in MATLAB R2021a. We
implemented C-FISTA for the Lasso models as detailed in Section 6.1. From the
state-of-the-art SLEP software package [22] we utilized glLeastR.m,
sgLeastR.m, overlapping_LeastR.m and sgLogisticR.m to solve the group, sparse-
group and overlapping sparse-group Lasso models and sparse-group logistic
regression model respectively. All iterative sequences in the algorithms were
initialized at zero. For comparison, all of the Lasso models were tested with
the same data matrix $A$ and response vectors $b^{i}$. Figure 1 displays the
convergence results for C-FISTA, SLEP and ADMM on the three Lasso models and
the sparse-group logistic regression model. In Figure 1, the y-axis measures
the absolute difference from the optimal objective value at each iteration of
the individual algorithms. We took the final value of the objective function
of C-FISTA as the exact optimal solution which agreed, as can be seen in
Figure 1, to within machine tolerance of the final iterates produced by ADMM
and SLEP. Note, we fine-tuned the augmented Lagrangian parameter $\rho$ in the
ADMM algorithm to achieve the optimal performance for the method.

The computational results displayed in Figure 1 clearly demonstrate the proven
global linear convergence of C-FISTA for all of the models. While ADMM and
SLEP both display linear convergence properties, C-FISTA significantly
outperformed both ADMM and SLEP in all the numerical tests. Furthermore,
C-FISTA maintained a robustness between the various Lasso models. While SLEP’s
convergence rate suffered as the subvector group sizes increased in the group
Lasso model, the convergence rate of C-FISTA did not suffer in the group or
sparse-group Lasso models as the group sizes increased. The results for the
overlapping sparse-group Lasso tests are of particular interest because SLEP
is considered the state-of-the-art first order approach for this model, and
C-FISTA was about two times quicker than SLEP in our tests. Lastly, the
numerical results from the sparse-group logistic regression tests clearly show
the vast improvement of C-FISTA over SLEP; C-FISTA demonstrated clear linear
convergence while SLEP required substantially more iterations to obtain a
solution accurate to within $\pm 1.0\times 10^{-10}$.

##  7 Conclusions

In this paper we developed an accelerated composite version of FISTA, C-FISTA,
which handles the composite optimization model (1) and is a generalization of
GFISTA [7, 8]. We proved global linear convergence for C-FISTA without having
a strongly convex objective function, and demonstrated through Fenchel duality
the breadth of convex models which are solvable via C-FISTA in Section 4.
Lastly, we demonstrated C-FISTA outperforms both ADMM and the state-of-the-art
software package SLEP on various Lasso and logistic regression models.

The following lines of directions could be of interest as future research
topics. First, the conditions underlying our global linear convergence results
may be weakened. Second, one may further study possible adaptive schemes to
implement the accelerated FISTA without requiring the exact knowledge of the
problem parameters. Finally, it is interesting to study how the new algorithm
behaves beyond the scope of convex optimization.

Figure 1: The numerical results for the Lasso models $(GL)$, $(SGL)$ and
$(OSGL)$ and the sparse-group logistic regression model $(SGLR)$ with three
different subvector groupings. The $y$-axis is given in logarithmic scale and
measures the absolute difference in the objective value and the optimal
objective value at each iteration; $\gamma=5$, $\gamma_{1}=\gamma_{2}=10$ and
$\gamma_{1}=\gamma_{2}=0.1$ for the $(GL)$, $(SGL)$ and $(OSGL)$ models
respectively while $\gamma_{1}=\gamma_{2}=0.01$ for the $(SGLR)$ model.

## References

  * [1] A. d’Aspremont, D. Scieur, A. Taylor, Acceleration Methods, January 2021, https://arxiv.org/abs/2101.09545v1
  * [2] A. Auslender, Asymptotic properties of the Fenchel’s dual functional and their applications to decomposition problems, J. Optim. Theory Appl., 73 (1992), pp. 427-450.
  * [3] A. Beck, S. Shtern, Linearly convergent away-step conditional gradient for non-strongly convex functions. Math. Program. 164, 1-27 (2017). https://doi.org/10.1007/s10107-016-1069-4
  * [4] A. Beck and M. Teboulle, A Fast Iterative Shrinkage-Threshold Algorithm Linear Inverse Problems, SIAM Journal on Imaging Sciences, 2, pages 183-202, 2009.
  * [5] J. V. Burke and M. C. Ferris. A Gauss-Newton method for convex composite optimization. Mathematical Programming, 71 (2, Ser. A): 179-194, 1995.
  * [6] E.J. Candés, X.D. Li, Y. Ma, and J. Wright, Robust Principal Component Analysis? Journal of the ACM, 58 (3), 2011.
  * [7] L. Calatroni and A. Chambolle, Backtracking strategies for accelerated descent methods with smooth composite objectives. SIAM Journal on Optimization, Society for Industrial and Applied Mathematics, 29 (3), pages 1772-1798, 2019.
  * [8] A. Chambolle and T. Pock, An introduction to continuous optimization for imaging, Acta Numerica 25, pages 161-319, 2016.
  * [9] V. Chandrasekaran, B. Recht, P.A. Parrilo, and A.S. Willsky, The convex geometry of linear inverse problems, Foundations of Computational Mathematics, 12 (6), pages 805-849, 2012.
  * [10] P. L. Combettes and V. R. Wajs, Signal recovery by proximal forward-backward splitting, Multiscale Model. Simul., 4, pages 1168-1200, 2005.
  * [11] I. Daubechies, M. Defrise, and C. D. Mol, An iterative thresholding algorithm for linear inverse problems with a sparsity constraint, Comm. Pure Appl. Math., 57, pages 1413-1457, 2004.
  * [12] D. Drusvyatskiy and A. S. Lewis, Error Bounds, Quadratic Growth, and Linear Convergence of Proximal Methods. Mathematics of Operations Research, 43, (3), pages 919-948, 2018.
  * [13] D. Drusvyatskiy and C. Paquette, Efficiency of minimizing compositions of convex functions and smooth maps, Mathematical Programming, 178, pages 503-558, 2019.
  * [14] M. A. T. Figueiredo and R. D. Nowak,An EM algorithm for wavelet-based image restoration, IEEE Trans. Image Process., 12, pages 906-916, 2003.
  * [15] M. I. Florea and S. A. Vorobyov, An Accelerated Composite Gradient Method for Large-scale Composite Objective Problems. IEEE Transactions on Signal Processing, 67(2), pages 444 - 459, 2019.
  * [16] M. I. Florea and S. A. Vorobyov,A generalized accelerated composite gradient method: Uniting Nesterov’s fast gradient method and FISTA, IEEE Transactions on Signal Processing, 68, pages 3033-3048, 2020.
  * [17] M. Fukushima, M. Haddou, H. Van Nguyen, et al. A parallel descent algorithm for convex programming. Comput Optim Applic 5, 5-37 (1996). https://doi.org/10.1007/BF00429749
  * [18] N. Gaffke and R. Mathar, A cyclic projection algorithm via duality. Metrika, 36, 29-54 (1989).
  * [19] S. P. Han and G. Lou, A Parallel Algorithm for a Class of Convex Programs, SIAM J. Control and Opt. 26, (2), 1988.
  * [20] M. Hong, Z.-Q. Luo, On the linear convergence of the alternating direction method of multipliers. Math. Program., 162, 165-199, 2017.
  * [21] L. Jacob, G. Obozinski, and J. Vert, Group lasso with overlap and graph lasso. In ICML, 2009.
  * [22] J. Liu and S. S. Ji and J. Ye SLEP: Sparse Learning and Efficient Projections, 2009. http://yelabs.net/software/SLEP/
  * [23] L. Meier, S. Geer, and P. Buhlmann. The group lasso for logistic regression. Journal of the Royal Statistical Society: Series B, 70:53-71, 2008.
  * [24] I. Necoara, Y. Nesterov & F. Glineur, Linear convergence of first order methods for non-strongly convex optimization. Mathematical Programming, 175, 69–107, 2019. https://doi.org/10.1007/s10107-018-1232-1
  * [25] Y. Nesterov, A method for solving the convex programming problem with convergence rate $\mathcal{O}(1/k^{2})$, Dokl. Akad. Nauk SSSR, 269 (1983), pp. 543-547 (in Russian).
  * [26] Y. Nesterov. Gradient methods for minimizing composite functions. Mathematical Programming, 140(1, Ser. B):125-161, 2013.
  * [27] M. Pham, A. Ninh, H. Le and Y. Liu, An Efficient Algorithm for Minimizing Multi Non-Smooth Component Functions. Journal of Computational and Graphical Statistics, 30:1, 162-170, 2021. DOI: 10.1080/10618600.2020.1804390
  * [28] Z.W. Qin, K. Scheinberg and D. Goldfarb, Efficient Block-Coordinate Descent Algorithms for the Group Lasso, Math. Prog. Comp. 5, pages 143-169, 2013.
  * [29] N. Rao, P. Shah, and S. Wright, Forward - Backward Greedy Algorithms for Atomic Norm Regularization, IEEE Transactions on Signal Processing, 63 (21), pages 5798-5811, 2015.
  * [30] N. Simon , J. Friedman , T. Hastie and R. Tibshirani, A Sparse-Group Lasso, Journal of Computational and Graphical Statistics, 22:2, 231-245, 2013, DOI: 10.1080/10618600.2012.681250.
  * [31] S. Rebegoldi and L. Calatroni, Scaled, inexact and adaptive generalized FISTA for strongly convex optimization. https://arxiv.org/abs/2101.03915, 2021.
  * [32] R. T. Rockafellar. Convex Analysis. Princeton Univ. Press, Princeton, second edition, 1970.
  * [33] K. Roos, M. Balvert, B. Gorissen, D. Hertog, A Universal and Structured Way to Derive Dual Optimization Problem Formulations, INFORMS Journal on Optimization 2(4), pages 229-255, 2020.https://doi.org/10.1287/ijoo.2019.0034
  * [34] R. Tibshirani, Regression shrinkage and selection via the Lasso. J. R. Stat. Soc. B 58, pages 267-288, 1996.
  * [35] S. J. Wright. Convergence of an inexact algorithm for composite nonsmooth optimization. IMA J. Numer. Anal., 10(3):299-321, 1990.
  * [36] L. Yuan, J. Liu and L. Ye, Efficient Methods for Overlapping Group Lasso, IEEE Trans Pattern Anal Mach Intell. 2013 Sep;35(9):2104-16. doi: 10.1109/TPAMI.2013.17. PMID: 23868773
  * [37] M. Yuan and Y. Lin, Model Selection and Estimation in Regression With Grouped Variables. Journal of the Royal Statistical Society, Series B, 68, 49-67, 2007.
  * [38] Y. J. Zhang, N. Zhang, D. F. Sun and K. C. Toh. An efficient Hessian based algorithm for solving large-scale sparse group Lasso problems, Mathematical Programming, 179(1): 223-263, 2020.
  * [39] X. Zhou, On the Fenchel Duality between Strong Convexity and Lipschitz Continuous Gradient, 2018, https://arxiv.org/abs/1803.06573v1.
  * [40] Z. Zhou, A. MC. So, A unified approach to error bounds for structured convex optimization problems. Mathematical Programming, 165, 689–728 (2017). https://doi.org/10.1007/s10107-016-1100-9