Faster Coordinate Descent via Adaptive Importance Sampling

Dmytro Perekrestenko Volkan Cevher Martin Jaggi

ETH Zurich EPFL EPFL

###### Abstract

Coordinate descent methods employ random partial updates of decision variables
in order to solve huge-scale convex optimization problems. In this work, we
introduce new adaptive rules for the random selection of their updates. By
adaptive, we mean that our selection rules are based on the dual residual or
the primal-dual gap estimates and can change at each iteration. We
theoretically characterize the performance of our selection rules and
demonstrate improvements over the state-of-the-art, and extend our theory and
algorithms to general convex objectives. Numerical evidence with hinge-loss
support vector machines and Lasso confirm that the practice follows the
theory.

##  1 Introduction

Coordinate descent methods rely on random partial updates of decision
variables for scalability. Indeed, due to their space and computational
efficiency as well as their ease of implementation, these methods are the
state-of-the-art for a wide selection of standard machine learning and signal
processing applications (Fu, 1998; Hsieh et al., 2008; Wright, 2015).

Basic coordinate descent methods sample an active coordinate set for
optimization uniformly at random, cf., Stochastic Dual Coordinate Ascent
(SDCA) Shalev-Shwartz and Zhang (2013) and other variants Friedman et al.
(2007, 2010); Shalev-Shwartz and Tewari (2011). However, recent results
suggest that by employing an appropriately defined _non-uniform_ fixed
sampling strategy, the convergence can be improved both in the theory as well
as in practice Zhao and Zhang (2014); Necorara et al. (2012); Nesterov (2012).

In this work, we show that we can surpass the existing convergence rates by
exploiting adaptive sampling strategies that change the sampling probability
distribution during each iteration. For this purpose, we adopt the primal-dual
framework of Dünner et al. (2016). In contast, however, we also handle convex
optimization problems with general convex regularizers without assuming strong
convexity of the regularizer.

In particular, we consider an adaptive coordinate-wise duality gap based
sampling. Hence, our work can be viewed as a natural continuation of the work
of Csiba et al. (2015), where the authors introduce an adaptive version of
SDCA for the smoothed hinge-loss support vector machine (SVM). However, our
work generalizes the gap-based adaptive criterion of Osokin et al. (2016) in a
nontrivial way to a broader convex optimization template of the following
form:

$\min_{\boldsymbol{\alpha}\in\mathbb{R}^{n}}\
f(A\boldsymbol{\alpha})+\sum_{i}g_{i}(\alpha_{i}),\vspace{-1mm}$ (1)

where $A$ is the data matrix, $f$ is a smooth convex function, and each
$g_{i}$ is a general convex function.

The template problem class in (1) includes not only smoothed hinge-loss SVM,
but also Lasso, Ridge Regression, (the dual formulation of the) original
hinge-loss SVM, Logistic Regression, etc. As a result, our theoretical results
for adaptive sampling can also recover the existing results for fixed non-
uniform Zhao and Zhang (2014) and uniform Dünner et al. (2016) sampling as
special cases.

Contributions. Our contributions are as follows:

  * •

We introduce new adaptive and fixed non-uniform sampling schemes for random
coordinate descent for problems for the template (1).

  * •

To our knowledge, we provide the first convergence rate analysis of coordinate
descent methods with adaptive sampling for problems with general convex
regularizer (i.e., the class in (1)).

  * •

We derive convergence guarantees with arbitrary sampling distributions for
both strongly convex and the general convex cases, and identify new
convergence improvements.

  * •

We support the theory with numerical evidence (i.e., Lasso and hinge-loss SVM)
and illustrate significant performance improvements in practice.

Outline: Section 2 provides basic theoretical preliminaries. Section 3
describes our theoretical results and introduces new sampling schemes. Section
4 discusses the application of our theory to machine learning, and compares
the computational complexity of proposed sampling methods. Section 5 provides
numerical evidence for the new methods. Section 6 discusses our contributions
in the light of existing work.

##  2 Preliminaries

We recall some concepts from convex optimization used in the sequel. The
_convex conjugate_ of a function $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is
defined as
$f^{*}(\boldsymbol{v}):=\sup_{\boldsymbol{u}\in\mathbb{R}^{n}}\boldsymbol{v}^{\top}\boldsymbol{u}-f(\boldsymbol{u})$.

######  Definition 2.1.

A function $f$: $\mathbb{R}^{n}\rightarrow\mathbb{R}\cup\\{+\infty\\}$ has
$B$-bounded support if its domain is bounded by $B$:

$f(\boldsymbol{u})<+\infty\rightarrow\|\boldsymbol{u}\|\leq B.$

######  Lemma 2.2 (Duality between Lipschitzness and L-Bounded Support, Dünner
et al. (2016)).

Given a proper convex function $g$, it holds that $g$ has $L$-bounded support
if and only if $g^{*}$ is $L$-Lipschitz.

###  2.1 Our primal-dual setup

In this paper, we develop coordinate descent methods for the following primal-
dual optimization pair:

$\displaystyle\min_{\boldsymbol{\alpha}\in\mathbb{R}^{n}}\Big{[}\mathcal{O}_{A}(\boldsymbol{\alpha}):=f(A\boldsymbol{\alpha})+\sum_{i}g_{i}(\alpha_{i})\Big{]},~{}$
(A)
$\displaystyle\min_{\boldsymbol{w}\in\mathbb{R}^{d}}\Big{[}\mathcal{O}_{B}(\boldsymbol{w}):=f^{*}(\boldsymbol{w})+\sum_{i}g_{i}^{*}(-\boldsymbol{a}_{i}^{\top}\boldsymbol{w})\Big{]},~{}$
(B)

where we have $A=[\boldsymbol{a}_{1},\dots,\boldsymbol{a}_{n}]$.

Our primal-dual template generalizes the standard primal-dual SDCA setup in
Shalev-Shwartz and Zhang (2013). As a result, we can provide gap certificates
for the quality of the numerical solutions while being applicable to broader
set of problems.

##### Optimality conditions.

The first-order optimality conditions for the problems (A) and (B) are given
by

$\displaystyle\boldsymbol{w}\in\partial f(A\boldsymbol{\alpha}),\hskip
19.91692pt-\boldsymbol{a}_{i}^{\top}\boldsymbol{w}\in\partial
g_{i}(\alpha_{i})$ $\displaystyle\forall i\in[n]$ (2) $\displaystyle
A\boldsymbol{\alpha}\in\partial f^{*}(\boldsymbol{w}),\hskip
19.91692pt\alpha_{i}\in\partial
g_{i}^{*}(-\boldsymbol{a}_{i}^{\top}\boldsymbol{w})$ $\displaystyle\forall
i\in[n]$

For a proof, see Bauschke and Combettes (2011).

##### Duality gap.

The duality gap is the difference between primal and dual solutions:

$G(\boldsymbol{\alpha},\boldsymbol{w}):=\mathcal{O}_{A}(\boldsymbol{\alpha})-(-\mathcal{O}_{B}(\boldsymbol{w})),$
(3)

which provides a certificate on the approximation accuracy both the primal and
dual objective values.

While always non-negative, under strong duality the gap reaches zero only in
an optimal pair $(\boldsymbol{\alpha}^{\star},\boldsymbol{w}^{\star})$. When
$f$ is differentiable the optimality conditions (2) write as
$\boldsymbol{w}^{\star}=\boldsymbol{w}(\boldsymbol{\alpha}^{\star})=\nabla
f(A\boldsymbol{\alpha}^{\star})$. We will rely on the following relationship
in our applications:

$\boldsymbol{w}=\boldsymbol{w}(\boldsymbol{\alpha}):=\nabla
f(A\boldsymbol{\alpha})\ .$

Choosing this mapping allows for running existing algorithms based on
$\boldsymbol{\alpha}$ alone, without the algorithm needing to take care of
$\boldsymbol{w}$. Nevertheless, the mapping allows us to express the gap
purely in terms of the original variable $\boldsymbol{\alpha}$, at any time:
We define
$G(\boldsymbol{\alpha}):=G(\boldsymbol{\alpha},\boldsymbol{w}(\boldsymbol{\alpha}))=\mathcal{O}_{A}(\boldsymbol{\alpha})-(-\mathcal{O}_{B}(\boldsymbol{w}(\boldsymbol{\alpha})))$.

##### Coordinate-wise duality gaps.

For our problem structure of partially separable problems (A) and (B), it is
not hard to show that the duality gap can be written as a sum of coordinate-
wise gaps:

$G(\boldsymbol{\alpha})=\sum_{i}G_{i}(\alpha_{i}):=\sum_{i}\Big{(}g_{i}^{*}(-\boldsymbol{a}_{i}^{\top}\boldsymbol{w})+g_{i}(\alpha_{i})+\alpha_{i}\boldsymbol{a}_{i}^{\top}\boldsymbol{w}\Big{)}$
(4)

The relation holds since our mapping $\boldsymbol{w}=\nabla
f(A\boldsymbol{\alpha})$ invokes the Fenchel-Young inequality for $f$ with an
equality. Moreover, the Fenchel-Young inequality for $g_{i}$ implies that all
$G_{i}(\alpha_{i})$’s are non-negative.

###  2.2 Dual residuals

We base our fixed non-uniform and adaptive schemes on the concept of “dual
residual,” i.e., a measure of progress to the optimum of the dual variables
$\boldsymbol{\alpha}$. Here we assume that $\boldsymbol{w}=\nabla
f(A\boldsymbol{\alpha})$.

######  Definition 2.3 (Dual Residual. A generalization of Csiba et al.
(2015)).

Consider the primal-dual setting (A)-(B). Let each $g_{i}$ be
$\mu_{i}$-strongly convex with convexity parameter $\mu_{i}\geq 0$ $\forall
i\in[n]$. For the case $\mu_{i}=0$ we require $g_{i}$ to have a bounded
support. Then, given $\boldsymbol{\alpha}$, the $i$-th dual residue on
iteration $t$ is given by:

$\kappa_{i}^{(t)}:=\min_{u\in\partial
g_{i}^{*}(-\boldsymbol{a}_{i}^{\top}\boldsymbol{w}^{(t)})}|u-\alpha_{i}^{(t)}|.$

######  Remark 2.4.

Note that for $u$ to be well defined, i.e., the subgradient in (6) not to be
empty, we need the domain of $g^{*}$ to be the whole space. For $\mu>0$ this
is given by strong convexity of $g_{i}$, while for $\mu_{i}=0$ this follows
from the bounded support assumption on $g_{i}$.

######  Definition 2.5 (Coherent probability vector, Csiba et al. (2015)).

We say that probability vector $\boldsymbol{p}^{(t)}\in\mathbb{R}^{n}$ is
coherent with the dual residue vector $\boldsymbol{\kappa}^{(t)}$ if for all
$i\in[n]$, we have $\kappa_{i}^{(t)}\neq 0\hskip 2.84544pt\Rightarrow\hskip
2.84544ptp_{i}^{(t)}>0.$

######  Definition 2.6 ($t$-support set).

We call the set

$I_{t}:=\\{i\in[n]:\kappa_{i}^{(t)}\neq 0\\}\subseteq[n]\vspace{-2mm}$

a $t$-support set.

######  Lemma 2.7.

Suppose that for each $i$, $g_{i}^{*}$ is $L_{i}$-Lipschitz. Then, $\forall
i:$ $|\kappa_{i}|\leq 2L_{i}$.

###### Proof.

By Lemma 2.2, the $L_{i}$-Lipschitzness $g_{i}^{*}$ implies $L_{i}$-bounded
support of $g_{i}(\alpha_{i})$ and therefore $|\alpha_{i}|\leq L_{i}$. By
writing Lipschitzness as bounded subgradient, $|u_{i}|\leq L_{i}$, and
$|\kappa_{i}|=|\alpha_{i}-u_{i}|\leq|\alpha_{i}|+|u_{i}|\leq 2L_{i}$. ∎

###  2.3 Coordinate Descent

Algorithm 1 describes the Coordinate Descent (CD) method in the primal-dual
setting (A) and (B). The method has $3$ major steps: Coordinate selection,
line-search along the chosen coordinate, and primal-dual parameter updates.
While the standard CD methods chooses the coordinates at random with fixed
distibutions, we develop adaptive strategies in the sequel that change the
sampling distribution per iteration. The other steps remain essentially the
same.

Algorithm 1 Coordinate Descent

1:Let
$\boldsymbol{\alpha}^{(0)}:=\boldsymbol{0}\in\mathbb{R}^{n},\boldsymbol{w}^{(0)}:=\boldsymbol{w}(\boldsymbol{\alpha}^{(0)})$

2:for $t=0,1,...,T$ do

3: Sample $i\in[n]$ randomly according to $\boldsymbol{p}^{(t)}$

4: Find $\Delta\alpha_{i}$ minimizing
$\mathcal{O}_{A}(\boldsymbol{\alpha}^{(t)}+\boldsymbol{e}_{i}\Delta\alpha_{i})$

5:
$\boldsymbol{\alpha}^{(t+1)}:=\boldsymbol{\alpha}^{(t)}+\boldsymbol{e}_{i}\Delta\alpha_{i}$

6: $\boldsymbol{w}^{(t+1)}:=\boldsymbol{w}(\boldsymbol{\alpha}^{(t+1)})$

7:end for

##  3 Adaptive Sampling-based CD

Our goal is to find a $\varepsilon_{B}$-suboptimal parameter $\boldsymbol{w}$
or $\varepsilon_{A}$-suboptimal parameter $\boldsymbol{\alpha}$, i.e.,
$\mathcal{O}_{A}(\boldsymbol{\alpha})-\mathcal{O}_{A}(\boldsymbol{\alpha}^{\star})\leq\varepsilon_{A}$
or
$\mathcal{O}_{B}(\boldsymbol{w})-\mathcal{O}_{B}(\boldsymbol{w}^{\star})\leq\varepsilon_{B}$,
for the following pair of dual optimization problems (A) and (B).

###  3.1 Key lemma

This subsection introduces a lemma that characterizes the relationship between
any sampling distribution for the coordinates, denoted as $\boldsymbol{p}$,
and the convergence rate of the CD method. For this purpose, we build upon the
(Csiba et al., 2015, Lemma 3) to relax the strong-convexity restrictions on
$g_{i}$’s. That is, we derive a convergence result for the general convex
$g_{i}$ with coordinate-dependent strong convexity constants $\mu_{i}$. In
contrast to (Csiba et al., 2015, Lemma 3), we can have $\mu_{i}=0$ when
$g_{i}$ has bounded support.

######  Lemma 3.1.

Let $f$ be $1/\beta$-smooth and each $g_{i}$ be $\mu_{i}$-strongly convex with
convexity parameter $\mu_{i}\geq 0$ $\forall i\in[n]$. For the case
$\mu_{i}=0$, we require $g_{i}$ to have a bounded support. Then for any
iteration $t$, any sampling distribution $\boldsymbol{p}^{(t)}$ and any
arbitrary $s_{i}\in[0,1]$ $\forall i\in[n]$, the iterates of the CD method
satisfy

$\displaystyle\mathbb{E}[\mathcal{O}_{A}(\boldsymbol{\alpha}^{(t+1)})\,|\,\boldsymbol{\alpha}^{(t)}]\leq\mathcal{O}_{A}(\boldsymbol{\alpha}^{(t)})-\sum_{i}s_{i}p_{i}^{(t)}G_{i}(\boldsymbol{\alpha}^{(t)})$
(5)
$\displaystyle-\sum_{i}p_{i}^{(t)}\Big{(}\frac{\mu_{i}(s_{i}-s_{i}^{2})}{2}-\frac{s_{i}^{2}\|\boldsymbol{a}_{i}\|^{2}}{2\beta}\Big{)}\lvert\kappa_{i}^{(t)}\rvert^{2},$

here $\kappa_{i}^{(t)}$ is $i$-th dual residual (see Def. 2.3).

The proof is provided in Appendix A.

######  Remark 3.2.

If in addition to the conditions of Lemma 3.1 we require
$\boldsymbol{p}^{(t)}$ to be coherent with $\boldsymbol{\kappa}^{(t)}$, then
for any $\theta\in[0,\min_{i\in I_{t}}p_{i}^{(t)}]$ it holds that

$\mathbb{E}[\mathcal{O}_{A}(\boldsymbol{\alpha}^{(t+1)})\,|\,\boldsymbol{\alpha}^{(t)}]\leq\mathcal{O}_{A}(\boldsymbol{\alpha}^{(t)})-\theta
G(\boldsymbol{\alpha}^{(t)})+\frac{\theta^{2}n^{2}}{2}F^{(t)},$ (6)
$F^{(t)}:=\frac{1}{n^{2}\beta\theta}\sum_{i\in
I_{t}}\Big{(}\frac{\theta(\mu_{i}\beta+\|\boldsymbol{a}_{i}\|^{2})}{p_{i}^{(t)}}-\mu_{i}\beta\Big{)}\lvert\kappa_{i}^{(t)}\rvert^{2}.$

###### Proof.

Since $s_{i}$ in Lemma 3.1 is an arbitrary number $\in[0,1]$, we take
$s_{i}=\frac{\theta}{p_{i}^{(t)}}$ for points with $i\in I_{t}$ and $s_{i}=0$
for all other points, here $\theta\in[0,\min_{i}p_{i}^{(t)}]$. Then, (5)
becomes the following, finalizing the proof:

$\displaystyle\mathbb{E}[\mathcal{O}_{A}(\boldsymbol{\alpha}^{(t+1)})\,|\,\boldsymbol{\alpha}^{(t)}]\leq\mathcal{O}_{A}(\boldsymbol{\alpha}^{(t)})-\theta\sum_{i\in
I_{t}}G_{i}(\boldsymbol{\alpha}^{(t)})$ $\displaystyle-\sum_{i\in
I_{t}}\Big{(}\frac{\mu_{i}\theta}{2}-\frac{\theta^{2}}{p_{i}^{(t)}}\frac{\mu_{i}\beta+\|\boldsymbol{a}_{i}\|^{2}}{2\beta}\Big{)}\lvert\kappa_{i}^{(t)}\rvert^{2}$
$\displaystyle=\mathcal{O}_{A}(\boldsymbol{\alpha}^{(t)})-\theta
G(\boldsymbol{\alpha}^{(t)})$ $\displaystyle-\frac{\theta}{2\beta}\sum_{i\in
I_{t}}\Big{(}\mu_{i}\beta-\frac{\theta(\mu_{i}\beta+\|\boldsymbol{a}_{i}\|^{2})}{p_{i}^{(t)}}\Big{)}\lvert\kappa_{i}^{(t)}\rvert^{2}$

∎

###  3.2 Why is the generalization important?

Consider the key lemma with $\mu_{i}=0$. Then, Remark 3.2 implies the
following:

$\mathbb{E}[\mathcal{O}_{A}(\boldsymbol{\alpha}^{(t+1)})\,|\,\boldsymbol{\alpha}^{(t)}]\leq\mathcal{O}_{A}(\boldsymbol{\alpha}^{(t)})-\theta
G(\boldsymbol{\alpha}^{(t)})+\frac{\theta^{2}n^{2}}{2}F^{(t)},$ (7)

where

$F^{(t)}:=\frac{1}{n^{2}\beta}\sum_{i\in
I_{t}}\Big{(}\frac{\lvert\kappa_{i}^{(t)}\rvert^{2}\|\boldsymbol{a}_{i}\|^{2}}{p_{i}^{(t)}}\Big{)}.$
(8)

Contrary to the strongly convex case, $F^{(t)}$ is positive. Hence, the
sampling distributions derived in Csiba et al. (2015) with strongly convex
$g_{i}$ are not optimal.

The following theorem generalizes (Zhao and Zhang, 2014, Theorem 5), and
(Dünner et al., 2016, Theorem 9) to allow adaptive sampling:

######  Theorem 3.3.

Assume $f$ is a $\frac{1}{\beta}$-smooth function. Then, if $g^{*}_{i}$ is
$L_{i}$-Lipschitz for each $i$ and $\boldsymbol{p}^{(t)}$ is coherent with
$\boldsymbol{\kappa}^{(t)}$, then the CD iterates satisfy

$\mathbb{E}[\varepsilon_{A}^{(t)}]\leq\frac{2F^{\circ}n^{2}+\frac{2\varepsilon_{A}^{(0)}}{p_{\min}}}{\frac{2}{p_{\min}}+t}.\vspace{-1mm}$
(9)

Moreover, we obtain a _duality gap_
$G(\bar{\boldsymbol{\alpha}})\leq\varepsilon$ after an overall number of
iterations $T$ whenever

$T\geq\max\Bigg{\\{}0,\frac{1}{p_{\min}}\log\Big{(}\frac{2\varepsilon_{A}^{(0)}}{n^{2}p_{\min}F^{\circ}}\Big{)}\Bigg{\\}}+\frac{5F^{\circ}n^{2}}{\varepsilon}-\frac{1}{p_{\min}}.$
(10)

Moreover, when $t\geq T_{0}$ with

$T_{0}:=\max\Bigg{\\{}0,\frac{1}{p_{\min}}\log\Big{(}\frac{2\varepsilon_{A}^{(0)}}{n^{2}p_{\min}F^{\circ}}\Big{)}\Bigg{\\}}+\frac{4F^{\circ}n^{2}}{\varepsilon}-\frac{2}{p_{\min}}$
(11)

we have the _suboptimality_ bound of
$\mathbb{E}[\mathcal{O}_{A}(\boldsymbol{\alpha}^{(t)})-\mathcal{O}_{A}(\boldsymbol{\alpha}^{\star})]\leq\varepsilon/2$.
Here $\varepsilon^{(0)}_{A}$ is the initial dual suboptimality and $F^{\circ}$
is an upper bound on $\mathbb{E}[F^{(t)}]$ taken over the random choice of the
sampled coordinate at $1,\dots,T_{0}$ algorithm iterations.

The proof is provided in Appendix A.

######  Remark 3.4.

We recover (Dünner et al., 2016, Theorem 9) as a special case of Theorem 3.3
by setting $p_{i}^{(t)}=\frac{1}{n}$. We recover (Zhao and Zhang, 2014,
Theorem 5) by setting $p_{i}^{(t)}=\frac{L_{i}}{\sum_{j}L_{j}}$.

####  3.2.1 Strategy I: Gap-wise sampling

Based on the results above, we first develop sampling strategies for the CD
method based on the decomposibly of the duality gap, i.e., sampling each
coordinate according to its duality gap.

######  Definition 3.5 (Nonuniformity measure, Osokin et al. (2016)).

The nonuniformity measure $\chi(\boldsymbol{x})$ of a vector
$\boldsymbol{x}\in\mathbb{R}^{n}$, is defined as:

$\chi(\boldsymbol{x}):=\sqrt{1+n^{2}\text{Var}[\boldsymbol{p}]},$

where $\boldsymbol{p}:=\frac{\boldsymbol{x}}{\|\boldsymbol{x}\|_{1}}$ is the
normalized probability vector.

######  Lemma 3.6.

Let $\boldsymbol{x}\in\mathbb{R}^{n}_{+}.$ Then, it holds that

$\|\boldsymbol{x}\|_{2}=\frac{\chi(\boldsymbol{x})}{\sqrt{n}}\|\boldsymbol{x}\|_{1}.\vspace{-2mm}$

###### Proof.

The proof follows from Def. 3.5 and

$\textstyle\text{Var}[\boldsymbol{p}]=\mathbb{E}[\boldsymbol{p}^{2}]-\mathbb{E}[\boldsymbol{p}]^{2}=\frac{1}{n}\|\boldsymbol{p}\|_{2}^{2}-\frac{1}{n^{2}}.\vspace{-5mm}$

∎

######  Theorem 3.7.

Let $f$ be a $\frac{1}{\beta}$-smooth function. Then, if $g^{*}_{i}$ is
$L_{i}$-Lipschitz for each $i$ and
$p_{i}^{(t)}:=\frac{G_{i}(\boldsymbol{\alpha}^{(t)})}{G(\boldsymbol{\alpha}^{(t)})}$,
then the iterations of the CD method satisfies

$\mathbb{E}[\varepsilon_{A}^{(t)}]\leq\frac{2F^{\circ}_{g}n^{2}+2n\varepsilon_{A}^{(0)}}{t+2n},$
(12)

where $F^{\circ}_{g}$ is an upper bound on
$\mathbb{E}\Big{[}F^{(t)}_{g}\Big{]},$ where the expectation is taken over the
random choice of the sampled coordinate at iterations $1,\dots,t$ of the
algorithm. Here $\overrightarrow{\boldsymbol{G}}$ and
$\overrightarrow{\boldsymbol{F}}$ are defined as follows:

$\overrightarrow{\boldsymbol{G}}:=(G_{i}(\boldsymbol{\alpha}^{(t)}))_{i=1}^{n},\hskip
14.22636pt\overrightarrow{\boldsymbol{F}}:=(\|\boldsymbol{a}_{i}\|^{2}\lvert\kappa_{i}^{(t)}\rvert^{2})_{i=1}^{n},$

and $F^{(t)}_{g}$ is defined analogously to (8):

$F^{(t)}_{g}:=\frac{\chi(\overrightarrow{\boldsymbol{F}})}{n\beta(\chi(\overrightarrow{\boldsymbol{G}}))^{3}}\sum_{i}\|\boldsymbol{a}_{i}\|^{2}\lvert\kappa_{i}^{(t)}\rvert^{2}.\vspace{-2mm}$
(13)

The proof is provided in Appendix A.

##### Gap-wise vs. Uniform Sampling:

Here we compare the rates obtained by Theorem 3.7 for gap-wise sampling and
Theorem 3.3 for uniform sampling. According to the Theorem 3.3, the rate for
any distribution can be written as follows

$\mathbb{E}[\varepsilon_{A}^{(t)}]\leq\frac{2F^{\circ}n^{2}+\frac{2\varepsilon_{A}^{(0)}}{p_{\min}}}{\frac{2}{p_{\min}}+t}=\frac{\frac{2}{\beta}\mathbb{E}\Big{[}\sum_{i}\frac{\lvert\kappa_{i}^{(t)}\rvert^{2}\|\boldsymbol{a}_{i}\|^{2}}{p_{i}^{t}}\Big{]}+\frac{2\varepsilon_{A}^{(0)}}{p_{\min}}}{\frac{2}{p_{\min}}+t}.$

For the uniform distribution ($p_{i}=1/n$), this yields

$\mathbb{E}[\varepsilon_{A}^{(t)}]\leq\frac{\frac{2n}{\beta}\mathbb{E}\Big{[}\sum_{i}\lvert\kappa_{i}^{(t)}\rvert^{2}\|\boldsymbol{a}_{i}\|^{2}\Big{]}+{2n\varepsilon_{A}^{(0)}}}{2n+t}.$
(14)

The rate of gap-wise sampling depends on non-uniformity measures
$\chi(\overrightarrow{\boldsymbol{G}})$ and
$\chi(\overrightarrow{\boldsymbol{F}})$:

$\mathbb{E}[\varepsilon_{A}^{(t)}]\leq\frac{\frac{2n}{\beta}\mathbb{E}\Big{[}\frac{\chi(\overrightarrow{\boldsymbol{F}})}{(\chi(\overrightarrow{\boldsymbol{G}}))^{3}}\sum_{i}\lvert\kappa_{i}^{(t)}\rvert^{2}\|\boldsymbol{a}_{i}\|^{2}\Big{]}+2n\varepsilon_{A}^{(0)}}{2n+t}.$

In the best case for gap-wise sampling the variance in
$(\lvert\kappa_{i}^{(t)}\rvert^{2}\|\boldsymbol{a}_{i}\|^{2})_{i=1}^{n}$ is
$0$, $\chi(\overrightarrow{\boldsymbol{F}})\approx 1$, and variance of gaps is
maximal $\chi(\overrightarrow{\boldsymbol{G}})\approx\sqrt{n}$. When this
condition holds, the convergence rate becomes the following:

$\mathbb{E}[\varepsilon_{A}^{(t)}]\leq\frac{\frac{2}{\beta\sqrt{n}}\mathbb{E}\Big{[}\sum_{i}\lvert\kappa_{i}^{(t)}\rvert^{2}\|\boldsymbol{a}_{i}\|^{2}\Big{]}+2n\varepsilon_{A}^{(0)}}{2n+t}.$

In the worst case scenario, when variance is maximal in
$(\lvert\kappa_{i}^{(t)}\rvert^{2}\|\boldsymbol{a}_{i}\|^{2})_{i=1}^{n}$,
$\chi(\overrightarrow{\boldsymbol{F}})\approx\sqrt{n}$, the rate of gap-wise
sampling is better than of uniform only when the gaps are non-uniform enough
i.e., $\chi(\overrightarrow{\boldsymbol{G}})\geq n^{\frac{1}{6}}$.

###  3.3 Strategy II: Adaptive & Uniform

Instead of minimizing (10), we here find an optimal sampling distribution as
to minimize our bound:

$T\geq\frac{5F^{\circ}n^{2}}{\varepsilon}+\frac{5\varepsilon_{A}^{(0)}}{\varepsilon
p_{\min}}.$ (15)

The number of iterations $T$ is directly proportional to $F^{\circ}$ and
$1/p_{\min}$. Therefore, the optimal distribution $\boldsymbol{p}$ should
minimize $F^{\circ}$ and $1/p_{\min}$ at the same time.

We denote the distribution minimizing $1/p_{\min}$ as supportSet uniform,
which is the following rule:

$p_{i}^{(t)}:=\begin{cases}\frac{1}{m_{t}},&\text{if }\kappa_{i}^{(t)}\neq
0\\\ 0,&\text{otherwise}.\end{cases}$ (16)

Above, $m_{t}$ is a cardinality of the support set on iteration $t$. The
distribution minimizing $F^{\circ}$, called adaptive:

$p_{i}^{(t)}:=\frac{\lvert\kappa_{i}^{(t)}\rvert\|\boldsymbol{a}_{i}\|}{\sum_{j}|\kappa_{j}^{(t)}|\|\boldsymbol{a}_{j}\|}.$
(17)

The mix of the two aforementioned distributions balances two terms and gives a
good suboptimal $T$ in (15). We define mixed distribution as:

$p_{i}^{(t)}:=\begin{cases}\frac{\sigma}{m_{t}}+(1-\sigma)\frac{\lvert\kappa_{i}^{(t)}\rvert\|\boldsymbol{a}_{i}\|}{\sum_{j}|\kappa_{j}^{(t)}|\|\boldsymbol{a}_{j}\|},&\text{if
}\kappa_{i}^{(t)}\neq 0\\\ 0,&\text{otherwise}\end{cases}\vspace{-1mm}$ (18)

where $\sigma\in[0,1]$. This distribution gives us the following bounds on
$F^{\circ}$ and $1/p_{\min}$:

$F^{\circ}_{\text{mix}}\leq\frac{F^{\circ}_{\text{ada}}}{1-\sigma}\hskip
28.45274pt\frac{1}{p_{\min}}\leq\frac{m}{\sigma}.\vspace{-1mm}$

and bound on the number of iterations:

$T\geq\frac{5F^{\circ}_{\text{ada}}n^{2}}{\varepsilon(1-\sigma)}+\frac{5\varepsilon_{A}^{(0)}m}{\varepsilon\sigma}.\vspace{-1mm}$
(19)

Above $m:=\max_{t}m_{t}$. Since the process of finding
$F^{\circ}_{\text{ada}}$ is rather problematic, a good $\sigma$ can be found
by replacing $F^{\circ}_{\text{ada}}$ with its upper bound and minimizing (18)
w.r.t. $\sigma$. Another option is to use a “safe” choice of $\sigma=0.5$, and
provide a balance between two distributions. This strategy benefits from
convergence guarantees in case of unknown $F^{\circ}_{\text{ada}}$. In the
applications section we use the latter option and call this sampling variant
ada-uniform sampling.

###  3.4 Variations along the theme

Based on the discussion above, we summarize our new variants of sampling
schemes for Algorithm 1:

  * •

uniform \- sample uniformly at random.

  * •

supportSet uniform \- sample uniformly at random inside the support set,
defined in (16). The distribution is recomputed every iteration.

  * •

adaptive \- sample adaptively based on dual residual, defined in (17). The
distribution is recomputed every iteration.

  * •

ada-uniform \- sample based on a mixture between supportSet uniform and
adaptive, defined in (18). The distribution is recomputed every iteration.

  * •

importance \- sample with a fixed non-uniform variant of adaptive obtained by
bounding $\kappa_{i}^{(t)}$ with $2L_{i}$ (Lemma 2.7):
$p_{i}:=\frac{L_{i}\|\boldsymbol{a}_{i}\|}{\sum_{j}L_{j}\|\boldsymbol{a}_{j}\|}$.
The distribution is computed only once. When the data is normalized, this
sampling variant coincides with ”importance sampling” of Zhao and Zhang
(2014).

  * •

ada-gap \- sample randomly based on coordinate-wise duality gaps, defined in
Section 3.2.1. The distribution is recomputed every iteration.

  * •

gap-per-epoch \- Use ada-gap but with updates per-epoch. The gap-based
distribution is only recomputed at the beginning of each epoch and stays fixed
during each epoch.

Full descriptions of the variants are in Appendix B.

##  4 Applications

###  4.1 Lasso and Sparse Logistic Regression

The Lasso and L1-regularized Logistic Regression are quintessential problems
with a general convex regularizer. Given a data matrix
$A=[\boldsymbol{a}_{1},\dots,\boldsymbol{a}_{n}]$ and a vector
$\boldsymbol{y}\in\mathbb{R}^{d}$, the Lasso is stated as:

$\min_{\boldsymbol{\alpha}\in\mathbb{R}^{n}}\|A\boldsymbol{\alpha}-\boldsymbol{y}\|^{2}_{2}+\lambda\|\boldsymbol{\alpha}\|_{1}.$
(20)

Both problems can easily be reformulated in our primal-dual setting (A)-(B),
by choosing $g_{i}(\alpha_{i}):=\lambda|\alpha_{i}|$. We have
$f(A\boldsymbol{\alpha}):=\|A\boldsymbol{\alpha}-\boldsymbol{y}\|^{2}_{2}$ for
Lasso, and $f(A\boldsymbol{\alpha})$ is the logistic loss for classification
respectively. To our knowledge, there is no importance sampling or adaptive
sampling techniques for CD in this setting.

##### Lipschitzing trick.

In order to have duality gap convergence guarantees (Theorem 3.3) we need
$g_{i}^{*}$ to be Lipschitz continuous, which however is not the case for
$g_{i}=|.|$ being the absolute value function. We modify the function $g_{i}$
without affecting the iterate sequence of CD using the “Lipschitzing trick”
from Dünner et al. (2016), as follows.

According to Lemma 2.2, a proper convex function $g_{i}$ has bounded support
if and only if $g_{i}^{*}$ is Lipschitz continuous. We modify
$g_{i}(\alpha_{i})=\lambda|\alpha_{i}|$ by restricting its support to the
interval with radius
$B:=\frac{1}{\lambda}(f(A\boldsymbol{\alpha}^{(0)})+\lambda\|\boldsymbol{\alpha}^{(0)}\|_{1})$.
Since Algorithm 1 is monotone, we can choose $B$ big enough to guarantee that
$\boldsymbol{\alpha}^{(t)}$ will stay inside the ball during optimization,
i.e. that the algorithm’s iterate sequence will not be affected by $B$. By
modifying $g_{i}$ to bounded support of size $B$, we guarantee $g_{i}^{*}$ to
be $B$-Lipschitz continuous.

$\bar{g}_{i}(\alpha_{i}):=\begin{cases}\lambda|\alpha_{i}|,&\text{if
}|\alpha_{i}|\leq B\\\ +\infty,&\text{otherwise}\end{cases}$

The conjugate of $\bar{g}_{i}$ is:

$\bar{g}_{i}^{*}(u_{i})=\max_{\alpha_{i}:|\alpha_{i}|\leq
B}u_{i}\alpha_{i}-\lambda|\alpha_{i}|=B\big{[}|u_{i}|-\lambda\big{]}_{+}.$

##### Duality gap.

Using the gap decomposition (4) we obtain coordinate-wise duality gaps for
modified Lasso and sparse logistic regression, which now depends on the chosen
parameter $B$:

$\displaystyle G(\boldsymbol{\alpha})$
$\displaystyle=\sum_{i}\Big{(}g_{i}^{*}(-\boldsymbol{a}_{i}^{\top}\boldsymbol{w})+g_{i}(\alpha_{i})+\alpha_{i}\boldsymbol{a}_{i}^{\top}\boldsymbol{w}\Big{)}$
(21)
$\displaystyle=\sum_{i}\Big{(}B\big{[}|\boldsymbol{a}_{i}^{\top}\boldsymbol{w}|-\lambda\big{]}_{+}+\lambda|\alpha_{i}|+\alpha_{i}\boldsymbol{a}_{i}^{\top}\boldsymbol{w}\Big{)}.$

###  4.2 Hinge-Loss SVM

Our framework directly covers the original hinge-loss SVM formulation. The
importance sampling technique Zhao and Zhang (2014) are not applicable to the
original hinge-loss, but relies on a smoothed version of the hinge-loss,
changing the problem.

When $\varphi_{i}(.)$ is the hinge-loss, defined as
$\varphi_{i}(s):=[1-sy_{i}]_{+}$, our framework is directly applicable by
mapping the SVM dual problem to our template (A), that is

$\min_{\boldsymbol{\alpha}\in\mathbb{R}^{n}}\mathcal{O}_{A}(\boldsymbol{\alpha}):=\frac{1}{n}\sum_{i=1}^{n}\varphi_{i}^{*}(-\alpha_{i})+\frac{\lambda}{2}\Big{\|}\frac{1}{\lambda
n}\sum_{i=1}^{n}\alpha_{i}\boldsymbol{a}_{i}\Big{\|}_{2}^{2}.$ (22)

The conjugate of the hinge-loss is
$\varphi_{i}^{*}(\alpha_{i})=\alpha_{i}y_{i}$, with $\alpha_{i}y_{i}\in[0,1]$.
In other words,
$g_{i}^{*}(-\boldsymbol{a}_{i}^{\top}\boldsymbol{w})=\frac{1}{n}\varphi_{i}(\boldsymbol{a}_{i}^{\top}\boldsymbol{w})$
in our notation, and
$f^{*}(\boldsymbol{w})=\frac{\lambda}{2}\|\boldsymbol{w}\|_{2}^{2}$.

##### Duality gap.

Section 2.1 shows that the duality gap decomposes into a sum of coordinate-
wise gaps.

###  4.3 Computational costs

We discuss the computational costs of the proposed variants under the
different sampling schemes. Table 1 states the costs in detail, where _nnz_ is
the number of non-zero entries in the data matrix $A$. In the table, one
_epoch_ means $n$ consecutive coordinate updates, where $n$ is the number
features in the Lasso, and is the number of datapoints in the SVM.

Table 1: A summary of computational costs Algorithm | Cost per Epoch
---|---
uniform | $\mathcal{O}(\text{nnz})$
importance | $\mathcal{O}(\text{nnz}+n\log(n))$
gap-per-epoch | $\mathcal{O}(\text{nnz}+n\log(n))$
supportSet-uniform | $\mathcal{O}(n\cdot\text{nnz})$
adaptive | $\mathcal{O}(n\cdot\text{nnz})$
ada-uniform | $\mathcal{O}(n\cdot\text{nnz})$
ada-gap | $\mathcal{O}(n\cdot\text{nnz})$

##### Sampling and probability update

In each iteration, we sample a coordinate from a non-uniform probability
distribution. While the straightforward approach requires $\Theta(n)$ per
sample, it is not hard to see that this can be improved to $\Theta(\log(n))$
when using a tree data structure to maintain the probability vector Nesterov
(2013); Shalev-Shwartz and Wexler (2016). The tree structure can be built in
$\mathcal{O}(n\log(n))$.

##### Variable update and distribution generation

Computing all dual residuals $\kappa_{i}$ or all coordinate-wise duality gaps
$G_{i}$ is as expensive as an epoch of the classic CD method, i.e., we need to
do $\Theta(\text{nnz})$ operations (one matrix-vector multiplication). In
contrast, updating one coordinate $\alpha_{i}$ is cheap, being
$\Theta(\text{nnz}/n)$.

Lasso, rcv1*

(a) Gap, fixed distr.

(b) Subopt., fixed distr.

(c) Gap, adaptive

(d) Suboptimality, adaptive

Lasso, mushrooms

(e) Gap, fixed distr.

(f) Subopt., fixed distr.

(g) Gap, adaptive

(h) Suboptimality, adaptive

SVM, ionosphere

(i) Gap, fixed distr.

(j) Subopt., fixed distr.

(k) Gap, adaptive

(l) Suboptimality, adaptive

Figure 1: Lasso (first two rows) and SVM (bottom row). Comparison of different
fixed and adaptive variants of CD, reporting duality gap and suboptimality
measures vs. epochs - rcv1*, mushrooms and ionosphere datasets.

##### Total cost per epoch

In a naive implementation, the most expensive sampling schemes are adaptive,
supportSet-uniform, ada-uniform and ada-gap. Those completely recompute the
sampling distribution after each iteration, giving a total per-epoch
complexity of $\mathcal{O}(n\cdot\text{nnz})$. In contrast, the fixed non-
uniform sampling scheme importance requires to build the sampling distribution
only once, or once per epoch for gap-per-epoch (both giving
$\mathcal{O}(\text{nnz})$ operations). The complexity of $n$ samplings using
the tree structure is $\mathcal{O}(n\log(n))$, the complexity of a variable
update is $\mathcal{O}(\text{nnz})$. Overall, the asymptotic complexity
therefore is $\mathcal{O}(n\log(n)+\text{nnz})$ per epoch, compared to
$\mathcal{O}(\text{nnz})$ for simple uniform sampling.

(a) Gap

(b) Suboptimality

Figure 2: Lasso on the mushrooms dataset. Performance in terms of duality gap
and suboptimality, plotted against the total number of vector operations.

##  5 Experimental results

We provide numerical evidence for our CD sampling strategies on two key
machine learning problems: The Lasso and the hinge-loss SVM. All our
algorithms and theory are also directly applicable to sparse logistic
regression and others, but we omit experiments due to space limitations.

##### Datasets.

The experiments are performed on three standard datasets listed in Table 2,
available111www.csie.ntu.edu.tw/ cjlin/libsvmtools/datasets/ from the UCI
repository Asuncion and Newman (2007). Note that rcv1* is a randomly
subsampled222We randomly picked 10000 datapoints and 1000 features, and then
removed zero rows and columns. version the rcv1 dataset. Experiments on the
full rcv1 dataset are provided in Appendix C.

Table 2: Datasets Dataset | $d$ | $n$ | nnz$/(nd)$ | $c_{v}=\frac{\mu(\|\boldsymbol{a}_{i}\|)}{\sigma(\|\boldsymbol{a}_{i}\|)}$
---|---|---|---|---
mushrooms | 112 | 8124 | 18.8% | 1.34
rcv1* | 809 | 7438 | 0.3% | 0.62
ionosphere | 351 | 33 | 88% | 3.07

##### Setup.

For Lasso, the regularization parameter $\lambda$ in (20) is set such that the
cardinality of the true support set is between 10% and 15 % of the total
number of features $n$. We use $\lambda=0.05$ for mushrooms, and
$\lambda=7\cdot 10^{-4}$ for rcv1*. For hinge-loss SVM, the regularization
parameter $\lambda$ is chosen such that the classification error on the test
set was comparable to training error. We use $\lambda=0.1$ for ionosphere.

##### Performance.

Figures 1 and 2 show the performance of all our studied variants of CD. We
record _suboptimality_ and _duality gap_ (see (3)) as the main measures of
algorithm performance. All reported results are averaged over 5 runs of each
algorithm.

##### Methods with fixed sampling distributions.

For the three efficient sampling schemes, our results show that CD importance
converges faster than CD uniform on both datasets for Lasso, however it is
worse than the uniform on SVM. The “mildly” adaptive strategy gap-per-epoch,
based on our coordinate-wise duality gap theory but computed only per-epoch,
significantly outperforms both of them. This is observed both in number of
epochs (Figure 1) as well as number of total vector operations (Figure 2).

##### Methods with adaptive sampling distributions

For the adaptive methods updating the probabilities after each coordinate step
in Figure 1, we show importance sampling as a baseline method (dashed line).
We see that measured per epoch, all adaptive methods outperform the fixed
sampling methods.Among all adaptive methods, the ada-gap algorithm shows
better convergence speed with both suboptimality and duality gap measures.

##### Highlights.

  * •

The experiments for Lasso show a clear advantage of non-uniform sampling over
the uniform, and superiority of the adaptive sampling over the fixed non-
uniform, which is supported by our theory.

  * •

Among the adaptive methods per iteration, the best performance for both Lasso
and SVM in terms of epochs is by ada-gap, which has proven convergence bounds
(Theorem 3.7), but also has high computational cost
($\Theta(d\cdot\text{nnz})$).

  * •

The best sampling scheme in terms of total computational cost is gap-per-
epoch, which is the epoch-wise variant of the ada-gap algorithm (based on
recomputing duality gaps once per epoch), as shown in Figure 2.

##  6 A discussion on the results

Coordinate descent methods have a rich history in the discipline of
optimization as well as many machine learning applications, cf., Wright (2015)
for a review.

For SVMs, CD related methods have been studied since their introduction, e.g.,
by Friess et al. (1998). Hsieh et al. (2008) is the first to propose CD in the
partially separable primal-dual setting for hinge-loss. Theoretical
convergence rates beyond the application of the hinge-loss can be found in the
SDCA line of work Shalev-Shwartz and Zhang (2013), which is the primal-dual
analog of the primal-only SGD algorithms. However, the main limitation of SDCA
is that it is only applicable to strongly convex regularizers, or requires
smoothing techniques Nesterov (2005) to be applicable to general regularizers
such as L1. The technique of Dünner et al. (2016) can extend to the CD
algorithms as well as the primal-dual analysis to the problem class of
interest here, using a bounded set of interest for the iterates instead of
relying on smoothing.

The convergence rate of stochastic methods (such as CD and SGD) naturally
depends on a sampling probability distribution over the coordinates or
datapoints respectively. While virtually all existing methods use sampling
uniformly at random Hsieh et al. (2008); Shalev-Shwartz and Tewari (2011);
Lacoste-Julien et al. (2013); Shalev-Shwartz and Zhang (2012), recently
Nesterov (2012); Qu and Richtárik (2016); Zhao and Zhang (2014); Allen-Zhu et
al. (2016) showed that an appropriately defined fixed non-uniform sampling
distribution, dubbed as importance sampling, can significantly improve the
convergence.

The work of Csiba et al. (2015) has taken the non-uniform sampling a step
further towards adaptive sampling. While restricted to strongly convex
regularizers, the rates provided for the AdaSDCA algorithm of Csiba et al.
(2015) \- when updating all probabilities after each step - can beat the ones
for uniform and importance sampling. A different approach for adapting the
sampling distribution is proposed in Osokin et al. (2016), where the block
coordinate Frank-Wolfe algorithm is enhanced with sampling proportional to
values of block-wise duality gaps. An adaptive variant of SGD is studied by
Papa et al. (2015), where they proposed an adaptive sampling scheme dependent
on the past iterations in a Markovian manner, without giving explicit
convergence rates. Other adaptive heuristics without proven convergence
guarantees include ACF Glasmachers and Dogan (2014) and ACiD Loshchilov et al.
(2011).

For general convex regularizers such as L1, the development of CD algorithms
includes Fu (1998); Friedman et al. (2007, 2010) and more recent extensions
also improving the theoretical convergence rates Shalev-Shwartz and Tewari
(2011); Johnson and Guestrin (2015); Zhao et al. (2014). All are restricted to
uniform sampling, and we are not aware of proven convergence rates showing
improvements of non-uniform or even adaptive sampling for unmodified L1
problems. Zhao and Zhang (2014); Allen-Zhu et al. (2016) show improved rates
for non-uniform sampling for L1 but require a smoothing modification of the
original problem, and are not covering adaptive sampling.

##### Conclusion.

In this work, we have investigated adaptive rules for adjusting the sampling
probabilities in coordinate descent. Our theoretical results provide improved
convergence rates for a more general class of algorithm schemes on one hand,
and optimization problems on the other hand, where we are able to directly
analyze CD on general convex objectives (as opposed to strongly convex
regularizers in previous works). This is particularly useful for L1 problems
and (original) hinge-loss objectives, which were not covered by previous
schemes. Our practical experiments confirm the strong performance of the
adaptive algorithms, and confirm that the behavior predicted by our theory.
Finally, we advocate the use of the computationally efficient gap-per-epoch
sampling scheme in practice. While the scheme is close to the ones supported
by our theory, an explicit primal-dual convergence analysis remains a future
research question.

##### Acknowledgements.

We thank Dominik Csiba and Peter Richtárik for helpful discussions. VC was
supported in part by the European Commission under Grant ERC Future Proof, SNF
200021-146750, and SNF CRSII2-147633.

## References

  * Allen-Zhu et al. [2016] Zeyuan Allen-Zhu, Zheng Qu, Peter Richtárik, and Yang Yuan.  Even Faster Accelerated Coordinate Descent Using Non-Uniform Sampling.  In _ICML 2016 - Proceedings of the 33th International Conference on Machine Learning_ , 2016.
  * Asuncion and Newman [2007] A Asuncion and DJ Newman.  UCI Machine Learning Repository.  _Miscellaneous_ , 2007.
  * Bauschke and Combettes [2011] Heinz H Bauschke and Patrick L Combettes.  _Convex Analysis and Monotone Operator Theory in Hilbert Spaces_.  CMS Books in Mathematics. 2011.
  * Csiba et al. [2015] Dominik Csiba, Zheng Qu, and Peter Richtárik.  Stochastic Dual Coordinate Ascent with Adaptive Probabilities.  In _ICML 2015 - Proceedings of the 32th International Conference on Machine Learning_ , 2015.
  * Dünner et al. [2016] Celestine Dünner, Simone Forte, Martin Takáč, and Martin Jaggi.  Primal-Dual Rates and Certificates.  In _ICML 2016 - Proceedings of the 33th International Conference on Machine Learning_ , 2016.
  * Friedman et al. [2007] Jerome Friedman, Trevor Hastie, Holger Höfling, and Robert Tibshirani.  Pathwise coordinate optimization.  _The Annals of Applied Statistics_ , 1(2):302–332, December 2007.
  * Friedman et al. [2010] Jerome Friedman, Trevor Hastie, and Robert Tibshirani.  Regularization Paths for Generalized Linear Models via Coordinate Descent.  _Journal of Statistical Software_ , 33(1):1–22, 2010.
  * Friess et al. [1998] Thilo-Thomas Friess, Nello Cristianini, and Colin Campbell.  The kernel-adatron algorithm: a fast and simple learning procedure for support vector machines.  _ICML 1998 - Proceedings of the 15th International Conference on Machine Learning_ , 1998.
  * Fu [1998] Wenjiang J. Fu.  Penalized regressions: The bridge versus the lasso.  _Journal of Computational and Graphical Statistics_ , 7(3):397–416, 1998.
  * Glasmachers and Dogan [2014] Tobias Glasmachers and Ürün Dogan.  Coordinate Descent with Online Adaptation of Coordinate Frequencies.  _arXiv_ , 2014.
  * Hsieh et al. [2008] Cho-Jui Hsieh, Kai-Wei Chang, Chih-Jen Lin, S Sathiya Keerthi, and S Sundararajan.  A Dual Coordinate Descent Method for Large-scale Linear SVM.  ICML - _the 25th International Conference on Machine Learning_ , 2008.
  * Johnson and Guestrin [2015] Tyler Johnson and Carlos Guestrin.  Blitz: A Principled Meta-Algorithm for Scaling Sparse Optimization.  In _ICML 2015 - Proceedings of the 32th International Conference on Machine Learning_ , 2015.
  * Lacoste-Julien et al. [2013] Simon Lacoste-Julien, Martin Jaggi, Mark Schmidt, and Patrick Pletscher.  Block-Coordinate Frank-Wolfe Optimization for Structural SVMs.  In _ICML 2013 - Proceedings of the 30th International Conference on Machine Learning_ , 2013.
  * Loshchilov et al. [2011] Ilya Loshchilov, Marc Schoenauer, and Michèle Sebag.  Adaptive Coordinate Descent.  In Natalio Krasnogor and Pier Luca Lanzi, editors, _Genetic and Evolutionary Computation Conference (GECCO 2011)_ , pages 885–992, Dublin, Ireland, 2011. ACM-SIGEVO.
  * Necorara et al. [2012] I Necorara, Yurii Nesterov, and François Glineur.  Efficiency of Randomized Coordinate Descent Methods on Optimization Problems With Linearly Coupled Constraints.  _Technical Report_ , 2012.
  * Nesterov [2005] Yurii Nesterov.  Smooth minimization of non-smooth functions.  _Mathematical Programming_ , 103(1):127–152, 2005\.
  * Nesterov [2012] Yurii Nesterov.  Efficiency of Coordinate Descent Methods on Huge-Scale Optimization Problems.  _SIAM Journal on Optimization_ , 22(2):341–362, 2012.
  * Nesterov [2013] Yurii Nesterov.  Subgradient methods for huge-scale optimization problems.  _Mathematical Programming_ , 146(1-2):275–297, 2013.
  * Osokin et al. [2016] Anton Osokin, Jean-Baptiste Alayrac, Isabella Lukasewitz, Puneet Dokania, and Simon Lacoste-Julien.  Minding the Gaps for Block Frank-Wolfe Optimization of Structured SVMs.  In _ICML 2016 - Proceedings of the 33th International Conference on Machine Learning_ , 2016.
  * Papa et al. [2015] Guillaume Papa, Pascal Bianchi, and Stéphan Clémençon.  Adaptive Sampling for Incremental Optimization Using Stochastic Gradient Descent.  In _ALT - 26th International Conference on Algorithmic Learning Theory_ , pages 317–331. 2015.
  * Qu and Richtárik [2016] Zheng Qu and Peter Richtárik.  Coordinate descent with arbitrary sampling I: algorithms and complexity.  _Optimization Methods and Software_ , 31(5):829–857, 2016.
  * Shalev-Shwartz and Tewari [2011] Shai Shalev-Shwartz and Ambuj Tewari.  Stochastic Methods for l1-regularized Loss Minimization.  _JMLR_ , 12:1865–1892, 2011.
  * Shalev-Shwartz and Wexler [2016] Shai Shalev-Shwartz and Yonatan Wexler.  Minimizing the Maximal Loss: How and Why?  _arXiv_ , 2016.
  * Shalev-Shwartz and Zhang [2012] Shai Shalev-Shwartz and Tong Zhang.  Proximal Stochastic Dual Coordinate Ascent.  _arXiv_ , 2012.
  * Shalev-Shwartz and Zhang [2013] Shai Shalev-Shwartz and Tong Zhang.  Stochastic Dual Coordinate Ascent Methods for Regularized Loss Minimization.  _JMLR_ , 14:567–599, 2013.
  * Wright [2015] Stephen J Wright.  Coordinate descent algorithms.  _Mathematical Programming_ , 151(1):3–34, 2015\.
  * Zhao and Zhang [2014] Peilin Zhao and Tong Zhang.  Stochastic Optimization with Importance Sampling.  _arXiv_ , 2014.
  * Zhao et al. [2014] Tuo Zhao, Han Liu, and Tong Zhang.  Pathwise Coordinate Optimization for Sparse Learning: Algorithm and Theory.  _arXiv_ , 2014.

##  Appendix A Proofs

######  Lemma’ 3.1.

Let $f$ be $1/\beta$-smooth and each $g_{i}$ be $\mu_{i}$-strongly convex with
convexity parameter $\mu_{i}\geq 0$ $\forall i\in[n]$. For the case
$\mu_{i}=0$, we require $g_{i}$ to have a bounded support. Then for any
iteration $t$, any sampling distribution $\boldsymbol{p}^{(t)}$ and any
arbitrary $s_{i}\in[0,1]$ $\forall i\in[n]$, the iterates of the CD method
satisfy

$\displaystyle\mathbb{E}[\mathcal{O}_{A}(\boldsymbol{\alpha}^{(t+1)})|\boldsymbol{\alpha}^{(t)}]\leq\mathcal{O}_{A}(\boldsymbol{\alpha}^{(t)})-\sum_{i}s_{i}p_{i}^{(t)}G_{i}(\boldsymbol{\alpha}^{(t)})$
(23)
$\displaystyle-\sum_{i}p_{i}^{(t)}\Big{(}\frac{\mu_{i}(s_{i}-s_{i}^{2})}{2}-\frac{s_{i}^{2}\|\boldsymbol{a}_{i}\|^{2}}{2\beta}\Big{)}|\kappa_{i}^{(t)}|^{2},$

here $\kappa_{i}^{(t)}$ is $i$-th dual residual.

###### Proof.

Since in CD update
($\boldsymbol{\alpha}^{(t+1)}=\boldsymbol{\alpha}^{(t)}+\boldsymbol{e}_{i}\Delta\alpha_{i}$)
only one coordinate per iteration is changed, the one iteration improvement in
dual objective can be written as:

$\displaystyle\mathcal{O}_{A}(\boldsymbol{\alpha}^{(t)})-\mathcal{O}_{A}(\boldsymbol{\alpha}^{(t+1)})$
$\displaystyle=\underbrace{\Big{[}g_{i}(\alpha_{i}^{(t)})+f(A\boldsymbol{\alpha}^{(t)})\Big{]}}_{(A)}-\underbrace{\Big{[}g_{i}(\alpha_{i}^{(t+1)})+f(A\boldsymbol{\alpha}^{(t+1)})\Big{]}}_{(B)}.$

To bound part $(B)$ we use a suboptimal update
$\Delta\alpha_{i}=s_{i}\kappa_{i}^{(t)}$, for all $s_{i}\in[0,1]$:

$\displaystyle(B)$
$\displaystyle=g_{i}(\alpha_{i}^{(t+1)})+f(A\boldsymbol{\alpha}^{(t+1)})$
$\displaystyle\leq\min_{\Delta\alpha_{i}}\Big{[}g_{i}(\alpha_{i}^{(t)}+\Delta\alpha_{i})+f(A\boldsymbol{\alpha}^{(t)}+\boldsymbol{a}_{i}\Delta\alpha_{i})\Big{]}$
$\displaystyle\leq
g_{i}(\alpha_{i}^{(t)}+s_{i}\kappa_{i}^{(t)})+f(A\boldsymbol{\alpha}^{(t)}+\boldsymbol{a}_{i}s_{i}\kappa_{i}^{(t)}).$

Each of $g_{i}$ is $\mu_{i}$-strongly convex, therefore:

$\displaystyle
g_{i}(\alpha_{i}^{(t)}+s_{i}\kappa_{i}^{(t)})=g_{i}(\alpha_{i}^{(t)}+s_{i}(u_{i}^{(t)}-\alpha_{i}^{(t)}))$
$\displaystyle=g_{i}(s_{i}(u_{i}^{(t)})+(1-s_{i})(\alpha_{i}^{(t)}))$
$\displaystyle\leq
s_{i}g_{i}(u_{i}^{(t)})+(1-s_{i})g_{i}(\alpha_{i}^{(t)})-\frac{\mu_{i}}{2}s_{i}(1-s_{i})(\kappa_{i}^{(t)})^{2}.$

The function $f$ is $\frac{1}{\beta}$-smooth:

$\displaystyle
f(A\boldsymbol{\alpha}^{(t)}+\boldsymbol{a}_{i}s_{i}\kappa_{i}^{(t)})$
$\displaystyle\leq f(A\boldsymbol{\alpha}^{(t)})+\nabla
f(A\boldsymbol{\alpha}^{(t)})^{\top}(s_{i}\kappa_{i}^{(t)}\boldsymbol{a}_{i})+\frac{1}{2\beta}\|s_{i}\kappa_{i}^{(t)}\boldsymbol{a}_{i}\|^{2}.$

As a result:

$\displaystyle(B)$ $\displaystyle\leq
s_{i}g_{i}(u_{i}^{(t)})-s_{i}g_{i}(\alpha_{i}^{(t)})-\frac{\mu_{i}}{2}s_{i}(1-s_{i})(\kappa_{i}^{(t)})^{2}$
$\displaystyle+\underbrace{g_{i}(\alpha_{i}^{(t)})+f(A\boldsymbol{\alpha}^{(t)})}_{(A)}+\nabla
f(A\boldsymbol{\alpha}^{(t)})^{\top}(s_{i}\kappa_{i}^{(t)}\boldsymbol{a}_{i})$
$\displaystyle+\frac{1}{2\beta}\|s_{i}\kappa_{i}^{(t)}\boldsymbol{a}_{i}\|^{2}.$

With obtained results above and optimality condition
$\boldsymbol{w}(\boldsymbol{\alpha})=\nabla f(A\boldsymbol{\alpha})$, the
improvement in dual objective can be written as:

$\displaystyle\mathcal{O}_{A}(\boldsymbol{\alpha}^{(t)})-\mathcal{O}_{A}(\boldsymbol{\alpha}^{(t+1)})=(A)-(B)$
$\displaystyle\geq-
s_{i}g_{i}(u_{i}^{(t)})+s_{i}g_{i}(\alpha_{i}^{(t)})+\frac{\mu_{i}}{2}s_{i}(1-s_{i})(\kappa_{i}^{(t)})^{2}$
$\displaystyle-\boldsymbol{w}(\boldsymbol{\alpha}^{(t)})(s_{i}u_{i}^{(t)}\boldsymbol{a}_{i})+\boldsymbol{w}(\boldsymbol{\alpha}^{(t)})(s_{i}\alpha_{i}^{(t)}\boldsymbol{a}_{i})-\frac{1}{2\beta}\|s_{i}\kappa_{i}^{(t)}\boldsymbol{a}_{i}\|^{2}$
$\displaystyle=s_{i}\Big{(}-g_{i}(u_{i}^{(t)})+g_{i}(\alpha_{i}^{(t)})-\boldsymbol{w}(\boldsymbol{\alpha}^{(t)})(u_{i}^{t-1}\boldsymbol{a}_{i})$
$\displaystyle+\boldsymbol{w}(\boldsymbol{\alpha}^{(t)})(\alpha_{i}^{t-1}\boldsymbol{a}_{i})+\frac{\mu_{i}}{2}(1-s_{i})(\kappa_{i}^{(t)})^{2}-\frac{s_{i}}{2\beta}\|\boldsymbol{a}_{i}\|^{2}|\kappa_{i}^{(t)}|^{2}\Big{)}.$

Since $u_{i}^{(t)}\in\partial
g_{i}^{*}(-\boldsymbol{a}_{i}^{\top}\boldsymbol{w}(\boldsymbol{\alpha}^{(t)}))$,
the Fenchel-Young inequality becomes equality for $g_{i}(u_{i}^{(t)})$:

$g_{i}(u_{i}^{(t)})+g_{i}^{*}(-\boldsymbol{a}_{i}^{\top}\boldsymbol{w}(\boldsymbol{\alpha}^{(t)}))=-\boldsymbol{w}(\boldsymbol{\alpha}^{(t)})(u_{i}^{(t)}\boldsymbol{a}_{i})$

Using this fact, the bound on the improvement in dual objective becomes:

$\displaystyle\mathcal{O}_{A}(\boldsymbol{\alpha}^{(t)})-\mathcal{O}_{A}(\boldsymbol{\alpha}^{(t+1)})\geq
s_{i}\Big{(}g_{i}(\alpha_{i}^{(t)})+g_{i}^{*}(-\boldsymbol{a}_{i}^{\top}\boldsymbol{w}(\boldsymbol{\alpha}^{(t)}))$
$\displaystyle+\boldsymbol{w}(\boldsymbol{\alpha}^{(t)})(\alpha_{i}^{(t)}\boldsymbol{a}_{i})+\frac{\mu_{i}}{2}(1-s_{i})(\kappa_{i}^{(t)})^{2}-\frac{s_{i}}{2\beta}\|\boldsymbol{a}_{i}\|^{2}|\kappa_{i}^{(t)}|^{2}\Big{)}$

Therefore for any $s_{i}\in[0,1]$ it holds that:

$\displaystyle\mathcal{O}_{A}(\boldsymbol{\alpha}^{(t)})-\mathcal{O}_{A}(\boldsymbol{\alpha}^{(t+1)})$
(24) $\displaystyle\geq
s_{i}\Big{[}G_{i}(\boldsymbol{\alpha}^{(t)})+\frac{\mu_{i}}{2}(1-s_{i})|\kappa_{i}^{(t)}|^{2}-\frac{s_{i}}{2\beta}\|\boldsymbol{a}_{i}\|^{2}|\kappa_{i}^{(t)}|^{2}\Big{]},$

where $G_{i}$ is $i$-th coordinate-wise duality gap:

$\displaystyle
G(\boldsymbol{\alpha}^{(t)})=\sum_{i}G_{i}(\boldsymbol{\alpha}^{(t)})$
$\displaystyle
G_{i}(\boldsymbol{\alpha}^{(t)})=g_{i}^{*}(-\boldsymbol{a}_{i}^{\top}\boldsymbol{w})+g_{i}(\alpha_{i}^{(t)})+\alpha_{i}^{(t)}\boldsymbol{a}_{i}^{\top}\boldsymbol{w}.$

By taking an expectation of the both sides with respect to $i$, conditioned on
$\boldsymbol{\alpha}^{(t)}$, we obtain:

$\displaystyle\mathbb{E}[\mathcal{O}_{A}(\boldsymbol{\alpha}^{(t+1)})|\boldsymbol{\alpha}^{(t)}]\leq\mathcal{O}_{A}(\boldsymbol{\alpha}^{(t)})-\sum_{i}s_{i}p_{i}^{(t)}G_{i}(\boldsymbol{\alpha}^{(t)})$
$\displaystyle-\sum_{i}p_{i}^{(t)}\Big{(}\frac{\mu_{i}(s_{i}-s_{i}^{2})}{2}-\frac{s_{i}^{2}\|\boldsymbol{a}_{i}\|^{2}}{2\beta}\Big{)}|\kappa_{i}^{(t)}|^{2}$

and thus finalize the proof. ∎

######  Theorem’ 3.3.

Assume $f$ is $\frac{1}{\beta}$-smooth function. Then, if $g^{*}_{i}$ is
$L_{i}$-Lipschitz for each $i$ and $\boldsymbol{p}^{(t)}$ is coherent with
$\boldsymbol{\kappa}^{(t)}$, then the iterates of the CD method satisfies

$\mathbb{E}[\varepsilon_{A}^{(t)}]\leq\frac{2F^{\circ}n^{2}+\frac{2\varepsilon_{A}^{(0)}}{p_{\min}}}{\frac{2}{p_{\min}}+t}.$
(25)

Moreover, the overall number of iterations $T$ to obtain a duality gap
$G(\bar{\boldsymbol{\alpha}})\leq\varepsilon$ must satisfy the following:

$T\geq\max\Bigg{\\{}0,\frac{1}{p_{\min}}\log\Big{(}\frac{2\varepsilon_{A}^{(0)}}{n^{2}p_{\min}F^{\circ}}\Big{)}\Bigg{\\}}+\frac{5F^{\circ}n^{2}}{\varepsilon}-\frac{1}{p_{\min}}.$
(26)

Moreover, when $t\geq T_{0}$ with

$T_{0}=\max\Bigg{\\{}0,\frac{1}{p_{\min}}\log\Big{(}\frac{2\varepsilon_{A}^{(0)}}{n^{2}p_{\min}F^{\circ}}\Big{)}\Bigg{\\}}+\frac{4F^{\circ}n^{2}}{\varepsilon}-\frac{2}{p_{\min}}$
(27)

we have the suboptimality bound of
$\mathbb{E}[\mathcal{O}_{A}(\boldsymbol{\alpha}^{(t)})-\mathcal{O}_{A}(\boldsymbol{\alpha}^{\star})]\leq\varepsilon/2$,
where $\varepsilon^{(0)}_{A}$ is the initial dual suboptimality and
$F^{\circ}$ is an upper bound on $\mathbb{E}[F^{(t)}]$ taken over the random
choice of the sampled coordinate at $1,\dots,T_{0}$ algorithm iterations.

###### Proof.

According to Lemma 2.2 if $g^{*}_{i}$ is $L_{i}$-Lipschitz, then $g_{i}$ has
$L_{i}$-bounded support and the conditions of Remark 3.2 are satisfied. From
Remark 3.2 we know:

$\mathbb{E}[\mathcal{O}_{A}(\boldsymbol{\alpha}^{(t+1)})|\boldsymbol{\alpha}^{(t)}]\leq\mathcal{O}_{A}(\boldsymbol{\alpha}^{(t)})-\theta
G(\boldsymbol{\alpha}^{(t)})+\frac{\theta^{2}n^{2}}{2}F^{(t)}$ (28)

With
$\mathcal{O}_{A}(\boldsymbol{\alpha}^{(t)})-\mathcal{O}_{A}(\boldsymbol{\alpha}^{(t+1)})=\varepsilon_{A}^{(t)}-\varepsilon_{A}^{(t+1)}$
and
$\varepsilon_{A}^{(t)}=\mathcal{O}_{A}(\boldsymbol{\alpha}^{(t)})-\mathcal{O}_{A}(\boldsymbol{\alpha}^{*})\leq
G(\boldsymbol{\alpha}^{(t)}),$ this implies:

$\mathbb{E}[\varepsilon_{A}^{(t+1)}|\boldsymbol{\alpha}^{(t)}]\geq\varepsilon_{A}^{(t)}-\theta\varepsilon_{A}^{(t)}+\frac{\theta^{2}n^{2}}{2}F^{(t)}$

by taking unconditional expectation over all iterations and using definition
of $F^{\circ}$ we obtain:

$\displaystyle\mathbb{E}[\varepsilon_{A}^{(t+1)}]$
$\displaystyle\leq(1-\theta)\mathbb{E}[\varepsilon_{A}^{(t)}]+\frac{\theta^{2}n^{2}}{2}\mathbb{E}[F^{(t)}]$
$\displaystyle\leq(1-\theta)\mathbb{E}[\varepsilon_{A}^{(t)}]+\frac{\theta^{2}n^{2}}{2}F^{\circ}$

Now we will show using induction that we can bound the dual suboptimality as:

$\mathbb{E}[\varepsilon_{A}^{(t)}]\leq\frac{2F^{\circ}n^{2}}{\frac{2}{p_{\min}}+t-t_{0}},$
(29)

where $t\geq
t_{0}=\max\Big{\\{}0,\frac{1}{p_{\min}}\log\big{(}\frac{2\varepsilon_{A}^{0}}{n^{2}p_{\min}F^{\circ}}\big{)}\Big{\\}}$.
Indeed, let’s choose $\theta=p_{\min}$, then the basis of induction at
$t=t_{0}$ is verified as:

$\displaystyle\mathbb{E}[\varepsilon_{A}^{(t)}]$
$\displaystyle\leq(1-p_{\min})^{t}\varepsilon_{A}^{(0)}+\sum_{i=0}^{t-1}(1-p_{\min})^{i}p_{\min}^{2}n^{2}\frac{F^{\circ}}{2}$
$\displaystyle\leq
e^{-tp_{\min}}\varepsilon_{A}^{(0)}+n^{2}p_{\min}\frac{F^{\circ}}{2}$
$\displaystyle\leq n^{2}p_{\min}F^{\circ}.$

Note that if in (29) instead of $F^{\circ}$ we take
$F^{\prime\circ}:=F^{\circ}+\frac{\varepsilon_{A}^{(0)}}{n^{2}p_{\min}},$ the
condition holds with $t_{0}=0$:

$\mathbb{E}[\varepsilon_{A}^{(t)}]\leq\frac{2F^{\prime\circ}n^{2}}{\frac{2}{p_{\min}}+t}=\frac{2F^{\circ}n^{2}+\frac{2\varepsilon_{A}^{(0)}}{p_{\min}}}{\frac{2}{p_{\min}}+t}.$

Now let’s prove the inductive step, for $t>t_{0}$. Suppose claim holds for
$t-1$, then

$\displaystyle\mathbb{E}[\varepsilon_{A}^{(t)}]$
$\displaystyle\leq(1-\theta)\mathbb{E}[\varepsilon_{A}^{(t-1)}]+\theta^{2}n^{2}\frac{F^{\circ}}{2}$
$\displaystyle\leq(1-\theta)\frac{2F^{\circ}n^{2}}{\frac{2}{p_{\min}}+(t-1)-t_{0}}+\theta^{2}n^{2}\frac{F^{\circ}}{2},$

choosing $\theta=\frac{2}{\frac{2}{p_{\min}}+t-1-t_{0}}\leq p_{\min}$ yeilds:

$\displaystyle\mathbb{E}[\varepsilon_{A}^{(t)}]$
$\displaystyle\leq\Bigg{(}1-\frac{2}{\frac{2}{p_{\min}}+t-1-t_{0}}\Bigg{)}\frac{2F^{\circ}n^{2}}{\frac{2}{p_{\min}}+(t-1)-t_{0}}$
$\displaystyle+\Bigg{(}\frac{2}{\frac{2}{p_{\min}}+t-1-t_{0}}\Bigg{)}^{2}\frac{F^{\circ}n^{2}}{2}$
$\displaystyle=\Bigg{(}1-\frac{2}{\frac{2}{p_{\min}}+t-1-t_{0}}\Bigg{)}\frac{2F^{\circ}n^{2}}{\frac{2}{p_{\min}}+(t-1)-t_{0}}$
$\displaystyle+\Bigg{(}\frac{1}{\frac{2}{p_{\min}}+t-1-t_{0}}\Bigg{)}\frac{2F^{\circ}n^{2}}{\frac{2}{p_{\min}}+(t-1)-t_{0}}$
$\displaystyle=\Bigg{(}1-\frac{1}{\frac{2}{p_{\min}}+t-1-t_{0}}\Bigg{)}\frac{2F^{\circ}n^{2}}{\frac{2}{p_{\min}}+(t-1)-t_{0}}$
$\displaystyle=\frac{2F^{\circ}n^{2}}{\frac{2}{p_{\min}}+(t-1)-t_{0}}\Bigg{(}\frac{\frac{2}{p_{\min}}+t-2-t_{0}}{\frac{2}{p_{\min}}+t-1-t_{0}}\Bigg{)}$
$\displaystyle\leq\frac{2F^{\circ}n^{2}}{\frac{2}{p_{\min}}+t-t_{0}}.$

This proves the bound (29) on suboptimality. To bound the duality gap we sum
the inequality (28) over the interval $t=T_{0}+1,...,T$ and obtain

$\displaystyle\mathbb{E}[\mathcal{O}_{A}(\boldsymbol{\alpha}^{(T_{0})})-\mathcal{O}_{A}(\boldsymbol{\alpha}^{(T)})]$
$\displaystyle\geq\theta\mathbb{E}\Big{[}\sum_{t=T_{0}+1}^{T}\mathcal{O}_{A}(\boldsymbol{\alpha}^{(t-1)})+\mathcal{O}_{B}(\boldsymbol{w}^{(t-1)})\Big{]}$
$\displaystyle-(T-T_{0})\frac{\theta^{2}n^{2}}{2}F^{\circ},$

by rearranging terms and choosing $\bar{\boldsymbol{w}}$ and
$\bar{\boldsymbol{\alpha}}$ to be the average vectors over
$t\in\\{T_{0},T-1\\}$ we get:

$\displaystyle\mathbb{E}[G(\bar{\boldsymbol{\alpha}})]=\mathbb{E}[\mathcal{O}_{A}(\bar{\boldsymbol{\alpha}})+\mathcal{O}_{B}(\bar{\boldsymbol{w}})]$
$\displaystyle\leq\frac{\mathbb{E}[\mathcal{O}_{A}(\boldsymbol{\alpha}^{(T_{0})})-\mathcal{O}_{A}(\boldsymbol{\alpha}^{(T)})]}{\theta(T-T_{0})}+\theta
n^{2}\frac{F^{\circ}}{2}.$

If $T\geq\frac{1}{p_{\min}}+T_{0}$ and $T_{0}\geq t_{0}$, we can set
$\theta=1/(T-T_{0})$ and combining this with (29) we get:

$\displaystyle\mathbb{E}[G(\bar{\boldsymbol{\alpha}})]$
$\displaystyle\leq\mathbb{E}[\mathcal{O}_{A}(\boldsymbol{\alpha}^{(T_{0})})-\mathcal{O}_{A}(\boldsymbol{\alpha}^{(T)})]+\frac{F^{\circ}n^{2}}{2(T-T_{0})}$
$\displaystyle\leq\mathbb{E}[\mathcal{O}_{A}(\boldsymbol{\alpha}^{(T_{0})})-\mathcal{O}_{A}(\boldsymbol{\alpha}^{*})]+\frac{F^{\circ}n^{2}}{2(T-T_{0})}$
$\displaystyle\leq\frac{2F_{T}n^{2}}{\frac{2}{p_{\min}}+t-t_{0}}+\frac{F^{\circ}n^{2}}{2(T-T_{0})}.$

A sufficient condition to bound the duality gap by $\varepsilon$ is that
$T_{0}\geq t_{0}-\frac{2}{p_{\min}}+\frac{4F^{\circ}n^{2}}{\varepsilon}$ and
$T\geq T_{0}+\frac{F^{\circ}n^{2}}{\varepsilon}$ which also implies
$\mathbb{E}[\mathcal{O}_{A}(\boldsymbol{\alpha}^{(T_{0})})-\mathcal{O}_{A}(\boldsymbol{\alpha}^{*})]\leq\varepsilon/2$.
Since we also need $T_{0}\geq t_{0}$ and $T-T_{0}\geq\frac{1}{p_{\min}}$, the
overall number of iterations should satisfy:

$T_{0}\geq\max\Big{\\{}t_{0},\frac{4F_{T}n^{2}}{\varepsilon}-\frac{2}{p_{\min}}+t_{0}\Big{\\}}$

and

$T-T_{0}\geq\max\Big{\\{}\frac{1}{p_{\min}},\frac{F^{\circ}n^{2}}{\varepsilon}\Big{\\}}.$

Using $a+b\geq\max(a,b)$ we finally can bound the total number of required
iterations to reach a duality gap of $\varepsilon$ by:

$\begin{aligned} T&\geq
T_{0}+\frac{1}{p_{\min}}+\frac{F^{\circ}n^{2}}{\varepsilon}\\\ &\geq
t_{0}+\frac{4F^{\circ}n^{2}}{\varepsilon}-\frac{1}{p_{\min}}+\frac{F^{\circ}n^{2}}{\varepsilon}\\\
&=t_{0}+\frac{5F^{\circ}n^{2}}{\varepsilon}-\frac{1}{p_{\min}}\end{aligned}.$

This concludes the proof. ∎

######  Theorem’ 3.7.

Let $f$ be a $\frac{1}{\beta}$-smooth function. Then, if $g^{*}_{i}$ is
$L_{i}$-Lipschitz for each $i$ and
$p_{i}^{(t)}:=\frac{G_{i}(\boldsymbol{\alpha}^{(t)})}{G(\boldsymbol{\alpha}^{(t)})}$,
then the iterations of the CD method satisfies

$\mathbb{E}[\varepsilon_{A}^{(t)}]\leq\frac{2F^{\circ}_{g}n^{2}+2n\varepsilon_{A}^{(0)}}{t+2n},$
(30)

where $F^{\circ}_{g}$ is an upper bound on
$\mathbb{E}\Big{[}F^{(t)}_{g}\Big{]},$ where the expectation is taken over the
random choice of the sampled coordinate at iterations $1,\dots,t$ of the
algorithm. Here $\overrightarrow{\boldsymbol{G}}$ and
$\overrightarrow{\boldsymbol{F}}$ are defined as:

$\overrightarrow{\boldsymbol{G}}:=(G_{i}(\boldsymbol{\alpha}^{(t)}))_{i=1}^{n},\hskip
14.22636pt\overrightarrow{\boldsymbol{F}}:=(\|\boldsymbol{a}_{i}\|^{2}|\kappa_{i}^{(t)}|^{2})_{i=1}^{n},$

and $F^{(t)}_{g}$ is defined analogously to (8):

$F^{(t)}_{g}:=\frac{\chi(\overrightarrow{\boldsymbol{F}})}{n\beta(\chi(\overrightarrow{\boldsymbol{G}}))^{3}}\sum_{i}\|\boldsymbol{a}_{i}\|^{2}|\kappa_{i}^{(t)}|^{2}.$
(31)

###### Proof.

We start from the result (23) of Lemma 3.1 when $\mu_{i}=0$:

$\displaystyle\mathbb{E}[\mathcal{O}_{A}(\boldsymbol{\alpha}^{(t+1)})|\boldsymbol{\alpha}^{(t)}]$
$\displaystyle\leq\mathcal{O}_{A}(\boldsymbol{\alpha}^{(t)})-\sum_{i}s_{i}p_{i}^{(t)}G_{i}(\boldsymbol{\alpha}^{(t)})$
$\displaystyle+\sum_{i}p_{i}^{(t)}\frac{s_{i}^{2}\|\boldsymbol{a}_{i}\|^{2}}{2\beta}|\kappa_{i}^{(t)}|^{2},$

by regrouping the elements and subtracting the optimal function value
$\mathcal{O}_{A}(\boldsymbol{\alpha}^{\star})$ from both sides we obtain:

$\displaystyle\mathbb{E}[\mathcal{O}_{A}(\boldsymbol{\alpha}^{(t+1)})-\mathcal{O}_{A}(\boldsymbol{\alpha}^{\star})|\boldsymbol{\alpha}^{(t)}]\leq\mathcal{O}_{A}(\boldsymbol{\alpha}^{(t)})-\mathcal{O}_{A}(\boldsymbol{\alpha}^{\star})$
$\displaystyle-\sum_{i}s_{i}p_{i}^{(t)}G_{i}(\boldsymbol{\alpha}^{(t)})+\sum_{i}\frac{p_{i}^{t}s_{i}^{2}}{2\beta}\|\boldsymbol{a}_{i}\|^{2}|\kappa_{i}^{(t)}|^{2}.$

With
$\varepsilon_{A}^{(t)}:=\mathcal{O}_{A}(\boldsymbol{\alpha}^{t})-\mathcal{O}_{A}(\boldsymbol{\alpha}^{*})$:

$\displaystyle\mathbb{E}[\varepsilon_{A}^{(t+1)}|\boldsymbol{\alpha}^{(t)}]$
$\displaystyle\leq\varepsilon_{A}^{(t)}-\sum_{i}s_{i}p_{i}^{(t)}G_{i}(\boldsymbol{\alpha}^{(t)})$
$\displaystyle+\sum_{i}\frac{p_{i}^{(t)}s_{i}^{2}}{2\beta}\|\boldsymbol{a}_{i}\|^{2}|\kappa_{i}^{(t)}|^{2}.$

We take
$p_{i}^{(t)}:=\frac{G_{i}(\boldsymbol{\alpha}^{(t)})}{G(\boldsymbol{\alpha}^{(t)})}$
and $s_{i}:=s$, then:

$\displaystyle\mathbb{E}[\varepsilon_{A}^{(t+1)}|\boldsymbol{\alpha}^{(t)}]$
$\displaystyle\leq\varepsilon_{A}^{(t)}-\frac{s}{G(\boldsymbol{\alpha}^{(t)})}\sum_{i}(G_{i}(\boldsymbol{\alpha}^{(t)}))^{2}$
$\displaystyle+\frac{s^{2}}{2\beta
G(\boldsymbol{\alpha}^{(t)})}\sum_{i}G_{i}(\boldsymbol{\alpha}^{(t)})\|\boldsymbol{a}_{i}\|^{2}|\kappa_{i}^{(t)}|^{2}.$

To simplify the following derivation we define a duality gap vector
$\overrightarrow{\boldsymbol{G}}:=(G_{i}(\boldsymbol{\alpha}^{(t)}))_{i=1}^{n}$
and residual vector
$\overrightarrow{\boldsymbol{F}}:=(\|\boldsymbol{a}_{i}\|^{2}|\kappa_{i}^{(t)}|^{2})_{i=1}^{n}$,
the inequality becomes:

$\mathbb{E}[\varepsilon_{A}^{(t+1)}|\boldsymbol{\alpha}^{(t)}]\leq\varepsilon_{A}^{(t)}-\frac{s}{G(\boldsymbol{\alpha}^{(t)})}\|\overrightarrow{\boldsymbol{G}}\|_{2}^{2}+\frac{s^{2}}{2\beta
G(\boldsymbol{\alpha}^{(t)})}\langle\overrightarrow{\boldsymbol{G}},\overrightarrow{\boldsymbol{F}}\rangle.$

By bounding the last term using the Cauchy-Schwarz inequality
$\langle\overrightarrow{\boldsymbol{G}},\overrightarrow{\boldsymbol{F}}\rangle\leq\|\overrightarrow{\boldsymbol{G}}\|_{2}\|\overrightarrow{\boldsymbol{F}}\|_{2}$
and using Lemma 3.6 we obtain:

$\displaystyle\mathbb{E}[\varepsilon_{A}^{(t+1)}|\boldsymbol{\alpha}^{(t)}]\leq\varepsilon_{A}^{(t)}-\frac{s}{G(\boldsymbol{\alpha}^{(t)})}\|\overrightarrow{\boldsymbol{G}}\|_{2}^{2}$
$\displaystyle+\frac{s^{2}}{2\beta
G(\boldsymbol{\alpha}^{(t)})}\|\overrightarrow{\boldsymbol{G}}\|_{2}\|\overrightarrow{\boldsymbol{F}}\|_{2}$
$\displaystyle=\varepsilon_{A}^{(t)}-\frac{sG(\boldsymbol{\alpha}^{(t)})(\chi(\overrightarrow{\boldsymbol{G}}))^{2}}{n}$
$\displaystyle+\frac{s^{2}\chi(\overrightarrow{\boldsymbol{G}})\chi(\overrightarrow{\boldsymbol{F}})\sum_{i}\|\boldsymbol{a}_{i}\|^{2}|\kappa_{i}^{(t)}|^{2}}{2n\beta}$
$\displaystyle\leq\varepsilon_{A}^{(t)}-\varepsilon_{A}^{(t)}\frac{(\chi(\overrightarrow{\boldsymbol{G}}))^{2}s}{n}+\frac{s^{2}\chi(\overrightarrow{\boldsymbol{G}})\chi(\overrightarrow{\boldsymbol{F}})\sum_{i}\|\boldsymbol{a}_{i}\|^{2}|\kappa_{i}^{(t)}|^{2}}{2n\beta}$
$\displaystyle=\Big{(}1-\frac{(\chi(\overrightarrow{\boldsymbol{G}}))^{2}s}{n}\Big{)}\varepsilon_{A}^{(t)}+\frac{s^{2}\chi(\overrightarrow{\boldsymbol{G}})\chi(\overrightarrow{\boldsymbol{F}})\sum_{i}\|\boldsymbol{a}_{i}\|^{2}|\kappa_{i}^{(t)}|^{2}}{2n\beta}.$

In the third line we have used weak duality, that is
$G(\boldsymbol{\alpha}^{(t)})\geq\varepsilon_{A}^{(t)}$. Analogously to the
proof of Theorem 3.3 we now prove that the suboptimality is bounded by:

$\mathbb{E}[\varepsilon_{A}^{(t)}]\leq\frac{2F^{\circ}_{g}n^{2}+2n\varepsilon_{A}^{(0)}}{t+2n},$
(32)

where

$F^{\circ}_{g}\geq\mathbb{E}\bigg{[}\frac{\chi(\overrightarrow{\boldsymbol{F}})}{n\beta(\chi(\overrightarrow{\boldsymbol{G}}))^{3}}\sum_{i}\|\boldsymbol{a}_{i}\|^{2}|\kappa_{i}^{(t)}|^{2}\bigg{]}.$

The basis of induction at $t=0$ obviously follows from the non-negativity of
$F^{\circ}_{g}$.
Now let us prove the induction step, assume that condition (32) holds at step
$t$, then by taking
$s:=\frac{2n}{(t+2n)(\chi(\overrightarrow{\boldsymbol{G}}))^{2}}$ we get:

$\displaystyle\mathbb{E}[\varepsilon_{A}^{(t+1)}|\boldsymbol{\alpha}^{(t)}]\leq\Big{(}1-\frac{(\chi(\overrightarrow{\boldsymbol{G}}))^{2}s}{n}\Big{)}\varepsilon_{A}^{(t)}$
(33)
$\displaystyle+\frac{s^{2}\chi(\overrightarrow{\boldsymbol{G}})\chi(\overrightarrow{\boldsymbol{F}})\sum_{i}\|\boldsymbol{a}_{i}\|^{2}|\kappa_{i}^{(t)}|^{2}}{2n\beta}$
$\displaystyle\leq\Big{(}1-\frac{2}{(t+2n)}\Big{)}\frac{2F^{\circ}_{g}n^{2}+2n\varepsilon_{A}^{(0)}}{t+2n}$
$\displaystyle+\frac{2n\chi(\overrightarrow{\boldsymbol{F}})\sum_{i}\|\boldsymbol{a}_{i}\|^{2}|\kappa_{i}^{(t)}|^{2}}{\beta(t+2n)^{2}(\chi(\overrightarrow{\boldsymbol{G}}))^{3}}.$

By taking an unconditional expectation of (33) and bounding by
$\hat{C}:=F^{\circ}_{g}n+\varepsilon_{A}^{(0)}$ we obtain:

$\displaystyle\mathbb{E}[\varepsilon_{A}^{(t+1)}]$
$\displaystyle\leq\Big{(}1-\frac{2}{(t+2n)}\Big{)}\frac{2n\hat{C}}{t+2n}$
$\displaystyle+\frac{2n}{(t+2n)^{2}}\mathbb{E}\bigg{[}\frac{\chi(\overrightarrow{\boldsymbol{F}})\sum_{i}\|\boldsymbol{a}_{i}\|^{2}|\kappa_{i}^{(t)}|^{2}}{\beta(\chi(\overrightarrow{\boldsymbol{G}}))^{3}}\bigg{]}$
$\displaystyle\leq\Big{(}1-\frac{2}{(t+2n)}\Big{)}\frac{2n\hat{C}}{t+2n}+\frac{2n\hat{C}}{(t+2n)^{2}}$
$\displaystyle=\frac{2n\hat{C}}{t+2n}\Big{(}1-\frac{2}{(t+2n)}+\frac{1}{(t+2n)}\Big{)}$
$\displaystyle=\frac{2n\hat{C}}{t+2n}\frac{t+2n-1}{t+2n}$
$\displaystyle\leq\frac{2n\hat{C}}{t+2n}\frac{t+2n}{t+2n+1}$
$\displaystyle=\frac{2n\hat{C}}{t+2n+1}.$

And this concludes the proof. ∎

##  Appendix B Algorithms

###  B.1 Algorithms with fixed sampling

In this subsection we give the coordinate descent algorithms with sampling
schemes with fixed probabilities which we derived the in previous chapter. The
basic Coordinate Descent with uniform sampling for Lasso was presented in
(Shalev-Shwartz and Tewari, 2011, Algorithm 1). Coordinate descent for hinge-
loss SVM was given in Shalev-Shwartz and Zhang (2013). Here we give these
algorithms along with their enhanced fixed non-uniform sampling versions of
importance sampling and heuristic gap-per-epoch sampling, which is based on
initial coordinate-wise duality gaps at the beginning of each epoch.

##### Lasso

See Algorithms 2, 3. To describe the algorithm the ”soft-threshold” function
$s_{\tau}(w)$ is defined:

$s_{\tau}(w):=sign(w)(|w|-\tau)_{+}=sign(w)\max\Big{\\{}|w|-\tau,0\Big{\\}}$

Algorithm 2 Coordinate Descent for Lasso (uniform & importance)

1:Choose $\text{mode}\in[\text{uniform},\text{importance}]$

2:let $\boldsymbol{\alpha}^{(0)}=0,$ $\boldsymbol{w}^{(0)}=\nabla
f(A\boldsymbol{\alpha}^{(0)})$

3:switch mode do

4: case uniform

5: $p_{i}:=\frac{1}{n}~{}~{}~{}\forall i$

6: case importance

7:
$p_{i}:=\frac{L_{i}\|\boldsymbol{a}_{i}\|}{\sum_{j}L_{j}\|\boldsymbol{a}_{j}\|}~{}~{}~{}\forall
i$

8:for t = 0,1,… do

9: sample $j$ from $[n]$ according to distribution $\boldsymbol{p}$

10: let $z_{j}=\frac{\partial
f(A\boldsymbol{\alpha}^{(t)})}{\partial\alpha_{j}}$

11: $\alpha_{j}^{(t+1)}=s_{\lambda}(\alpha_{j}^{(t)}-z_{j})$

12: $\boldsymbol{w}^{(t+1)}=\nabla f(A\boldsymbol{\alpha}^{(t+1)})$

13:end for

Algorithm 3 Coordinate Descent for Lasso (gap-per-epoch)

1:let $\boldsymbol{\alpha}^{(0)}=0,$ $\boldsymbol{w}^{(0)}=\nabla
f(A\boldsymbol{\alpha}^{(0)})$

2:for t = 0,1,… do

3: if $\text{mod}(t,n)==0$ then

4: generate probabilities distribution $\boldsymbol{p}^{(t)}$:
$p_{i}=\frac{B\big{[}|\boldsymbol{a}_{i}^{\top}\boldsymbol{w}^{(0)}|-\lambda\big{]}_{+}+\lambda|\alpha_{i}^{(0)}|+\alpha_{i}^{(0)}\boldsymbol{a}_{i}^{\top}\boldsymbol{w}^{(0)}}{\sum_{j}\Big{(}B\big{[}|\boldsymbol{a}_{j}^{\top}\boldsymbol{w}^{(0)}|-\lambda\big{]}_{+}+\lambda|\alpha_{j}^{(0)}|+\alpha_{j}^{(0)}\boldsymbol{a}_{j}^{\top}\boldsymbol{w}^{(0)}\Big{)}}$

5: end if

6: sample $j$ from $[n]$ according to distribution $\boldsymbol{p}$

7: let $z_{j}=\frac{\partial
f(A\boldsymbol{\alpha}^{(t)})}{\partial\alpha_{j}}$

8: $\alpha_{j}^{(t+1)}=s_{\lambda}(\alpha_{j}^{(t)}-z_{j})$

9: $\boldsymbol{w}^{(t+1)}=\nabla f(A\boldsymbol{\alpha}^{(t+1)})$

10:end for

##### Hinge-Loss SVM

See Algorithms 4, 5

Algorithm 4 Stochastic Dual Coordinate Descent (uniform & importance)

1:Choose $\text{mode}\in[\text{uniform},\text{importance}]$

2:let $\boldsymbol{\alpha}^{(0)}=0,$ $\boldsymbol{w}^{(0)}=0$

3:switch mode do

4: case uniform

5: $p_{i}:=\frac{1}{n}~{}~{}~{}\forall i$

6: case importance

7:
$p_{i}:=\frac{L_{i}\|\boldsymbol{a}_{i}\|}{\sum_{j}L_{j}\|\boldsymbol{a}_{j}\|}~{}~{}~{}\forall
i$

8:for t = 0,1,… do

9: sample $j$ from $[n]$ according to distribution $\boldsymbol{p}$

10:
$\Delta\alpha_{j}=y_{j}\max\bigg{(}0,\min\Big{(}1,\frac{1-y_{j}\boldsymbol{a}_{j}^{\top}\boldsymbol{w}^{(t)}}{\|\boldsymbol{a}_{j}\|^{2}/(\lambda
n)}+y_{j}\alpha_{j}^{(t)}\Big{)}\bigg{)}-\alpha_{j}^{(t)}$

11:
$\boldsymbol{\alpha}^{(t+1)}=\boldsymbol{\alpha}^{(t)}+\Delta\alpha_{j}\boldsymbol{e}_{j}$

12: $\boldsymbol{w}^{(t+1)}=\boldsymbol{w}^{(t)}+(\lambda
n)^{-1}\Delta\alpha_{j}\boldsymbol{a}_{j}$

13:end for

Algorithm 5 Stochastic Dual Coordinate Descent (gap-per-epoch)

1:let $\boldsymbol{\alpha}^{(0)}=0,$ $\boldsymbol{w}^{(0)}=0$

2:for t = 0,1,… do

3: if $\text{mod}(t,n)==0$ then

4: generate probabilities distribution $\boldsymbol{p}^{(t)}$:
$p_{i}=\frac{\varphi_{i}(\boldsymbol{a}_{i}^{\top}\boldsymbol{w})-\alpha_{i}y_{i}+\alpha_{i}\boldsymbol{a}_{i}^{\top}\boldsymbol{w}}{\sum_{j=1}^{n}\Big{(}\varphi_{j}(\boldsymbol{a}_{j}^{\top}\boldsymbol{w})-\alpha_{j}y_{j}+\alpha_{j}\boldsymbol{a}_{j}^{\top}\boldsymbol{w}\Big{)}}$

5: end if

6: sample $j$ from $[n]$ according to distribution $\boldsymbol{p}$

7:
$\Delta\alpha_{j}=y_{j}\max\bigg{(}0,\min\Big{(}1,\frac{1-y_{j}\boldsymbol{a}_{j}^{\top}\boldsymbol{w}^{(t)}}{\|\boldsymbol{a}_{j}\|^{2}/(\lambda
n)}+y_{j}\alpha_{j}^{(t)}\Big{)}\bigg{)}-\alpha_{j}^{(t)}$

8:
$\boldsymbol{\alpha}^{(t+1)}=\boldsymbol{\alpha}^{(t)}+\Delta\alpha_{j}\boldsymbol{e}_{j}$

9: $\boldsymbol{w}^{(t+1)}=\boldsymbol{w}^{(t)}+(\lambda
n)^{-1}\Delta\alpha_{j}\boldsymbol{a}_{j}$

10:end for

###  B.2 Algorithms with adaptive sampling

In this subsection we consider Coordinate Descent with adaptive sampling
schemes. Here we present 4 different schemes:

  * •

supportSet-uniform discussed in Section 3.3, defined in (16).

  * •

adaptive discussed in Section 3.3, defined in (17).

  * •

ada-uniform discussed in Section 3.3, defined in (18).

  * •

ada-gap discussed in Section 3.2.1, defined in Theorem 3.7.

The algorithms with aforementioned sampling schemes are given below. For Lasso
see Algorithms 6, 7, 8 and 9. For hinge-loss SVM see Algorithms 10, 11, 12 and
13.

Algorithm 6 Coordinate Descent (supportSet-uniform)

1:let $\boldsymbol{\alpha}^{(0)}=0,$ $\boldsymbol{w}^{(0)}=\nabla
f(A\boldsymbol{\alpha}^{(0)})$

2:for t = 0,1,… do

3: calculate absolute values of dual residuals $|\kappa_{j}^{(t)}|$ for all
$j\in[n]$
$|\kappa_{j}^{(t)}|:=\Big{|}\alpha_{j}-B\cdot\text{sign}(\boldsymbol{a}_{i}^{\top}\boldsymbol{w}^{(t)})\cdot\big{[}|\boldsymbol{a}_{j}^{\top}\boldsymbol{w}^{(t)}|-\lambda\big{]}_{+}\Big{|}$

4: find $t$-support set $I_{t}=\\{i\in[n]:\kappa_{i}^{(t)}\neq
0\\}\subseteq[n]$

5: generate adaptive probabilities distribution $\boldsymbol{p}^{(t)}$, for
each $i\in[n]$:  $p_{i}^{(t)}:=\begin{cases}\frac{1}{|I_{t}|},&\text{if
}\kappa_{i}^{(t)}\neq 0\\\ 0,&\text{otherwise}\end{cases}$

6: sample $j$ from $[n]$ according to $\boldsymbol{p}^{(t)}$

7: let $z_{j}=\frac{\partial
f(A\boldsymbol{\alpha}^{(t)})}{\partial\alpha_{j}}$

8: $\alpha_{j}^{(t+1)}=s_{\lambda}(\alpha_{j}^{(t)}-z_{j})$

9: $\boldsymbol{w}^{(t+1)}:=\nabla f(A\boldsymbol{\alpha}^{(t+1)})$

10:end for

Algorithm 7 Coordinate Descent (adaptive)

1:let $\boldsymbol{\alpha}^{(0)}=0,$ $\boldsymbol{w}^{(0)}=\nabla
f(A\boldsymbol{\alpha}^{(0)})$

2:for t = 0,1,… do

3: calculate absolute values of dual residuals $|\kappa_{j}^{(t)}|$ for all
$j\in[n]$
$|\kappa_{j}^{(t)}|=\Big{|}\alpha_{j}-B\cdot\text{sign}(\boldsymbol{a}_{i}^{\top}\boldsymbol{w}^{(t)})\cdot\big{[}|\boldsymbol{a}_{j}^{\top}\boldsymbol{w}^{(t)}|-\lambda\big{]}_{+}\Big{|}$

4: generate adaptive probabilities distribution $\boldsymbol{p}^{(t)}$:
$p_{i}^{(t)}=\frac{|\kappa_{i}^{(t)}|\|\boldsymbol{a}_{i}\|}{\sum_{j}|\kappa_{j}^{(t)}|\|\boldsymbol{a}_{j}\|}$

5: sample $j$ from $[n]$ according to $\boldsymbol{p}^{(t)}$

6: let $z_{j}=\frac{\partial
f(A\boldsymbol{\alpha}^{(t)})}{\partial\alpha_{j}}$

7: $\alpha_{j}^{(t+1)}=s_{\lambda}(\alpha_{j}^{(t)}-z_{j})$

8: $\boldsymbol{w}^{(t+1)}=\nabla f(A\boldsymbol{\alpha}^{(t+1)})$

9:end for

Algorithm 8 Coordinate Descent (ada-uniform)

1:let $\boldsymbol{\alpha}^{(0)}=0,$ $\boldsymbol{w}^{(0)}=\nabla
f(A\boldsymbol{\alpha}^{(0)})$

2:for t = 0,1,… do

3: calculate absolute values of dual residuals $|\kappa_{j}^{(t)}|$ for all
$j\in[n]$
$|\kappa_{j}^{(t)}|=\Big{|}\alpha_{j}-B\cdot\text{sign}(\boldsymbol{a}_{i}^{\top}\boldsymbol{w}^{(t)})\cdot\big{[}|\boldsymbol{a}_{j}^{\top}\boldsymbol{w}^{(t)}|-\lambda\big{]}_{+}\Big{|}$

4: find $t$-support set $I_{t}=\\{i\in[n]:\kappa_{i}^{(t)}\neq
0\\}\subseteq[n]$

5: generate adaptive probabilities distribution $\boldsymbol{p}^{(t)}$:
$\begin{cases}p_{i}^{(t)}=\frac{1}{2|I_{t}|}+\frac{|\kappa_{i}^{(t)}|\|\boldsymbol{a}_{i}\|}{2\sum_{j}|\kappa_{j}^{(t)}|\|\boldsymbol{a}_{j}\|},&\text{if
}\kappa_{i}^{(t)}\neq 0\\\ p_{i}^{(t)}=0,&\text{otherwise}\end{cases}$

6: sample $j$ from $[n]$ according to $\boldsymbol{p}^{(t)}$

7: let $z_{j}=\frac{\partial
f(A\boldsymbol{\alpha}^{(t)})}{\partial\alpha_{j}}$

8: $\alpha_{j}^{(t+1)}=s_{\lambda}(\alpha_{j}^{(t)}-z_{j})$

9: $\boldsymbol{w}^{(t+1)}=\nabla f(A\boldsymbol{\alpha}^{(t+1)})$

10:end for

Algorithm 9 Coordinate Descent (ada-gap)

1:let $\boldsymbol{\alpha}^{(0)}=0,$ $\boldsymbol{w}^{(0)}=\nabla
f(A\boldsymbol{\alpha}^{(0)})$

2:for t = 0,1,… do

3: calculate feature-wise duality gaps $G_{j}^{(t)}$ for all $j\in[n]$
$G_{j}^{(t)}=B\big{[}|\boldsymbol{a}_{i}^{\top}\boldsymbol{w}^{(t)}|-\lambda\big{]}_{+}+\lambda|\alpha_{i}^{(t)}|+\alpha_{i}^{(t)}\boldsymbol{a}_{i}^{\top}\boldsymbol{w}^{(t)}$

4: generate adaptive probabilities distribution $\boldsymbol{p}^{(t)}$:
$p_{i}^{(t)}=\frac{G_{i}^{(t)}}{\sum_{j}G_{j}^{(t)}}$

5: sample $j$ from $[n]$ according to $\boldsymbol{p}^{(t)}$

6: let $z_{j}=\frac{\partial
f(A\boldsymbol{\alpha}^{(t)})}{\partial\alpha_{j}}$

7: $\alpha_{j}^{(t+1)}=s_{\lambda}(\alpha_{j}^{(t)}-z_{j})$

8: $\boldsymbol{w}^{(t+1)}=\nabla f(A\boldsymbol{\alpha}^{(t+1)})$

9:end for

Algorithm 10 Stochastic Dual Coordinate Descent (supportSet-uniform)

1:let $\boldsymbol{\alpha}^{(0)}=0,$ $\boldsymbol{w}^{(0)}=0$

2:for t = 0,1,… do

3: calculate absolute values of dual residuals $|\kappa_{j}^{(t)}|$ for all
$j\in[n]$

4: find $t$-support set $I_{t}=\\{i\in[n]:\kappa_{i}^{(t)}\neq
0\\}\subseteq[n]$

5: generate adaptive probabilities distribution $\boldsymbol{p}^{(t)}$, for
each $i\in[n]$:  $p_{i}^{(t)}:=\begin{cases}\frac{1}{|I_{t}|},&\text{if
}\kappa_{i}^{(t)}\neq 0\\\ 0,&\text{otherwise}\end{cases}$

6: sample $j$ from $[n]$ according to distribution $\boldsymbol{p}^{(t)}$

7:
$\Delta\alpha_{j}=y_{j}\max\bigg{(}0,\min\Big{(}1,\frac{1-y_{j}\boldsymbol{a}_{j}^{\top}\boldsymbol{w}^{(t)}}{\|\boldsymbol{a}_{j}\|^{2}/(\lambda
n)}+y_{j}\alpha_{j}^{(t)}\Big{)}\bigg{)}-\alpha_{j}^{(t)}$

8:
$\boldsymbol{\alpha}^{(t+1)}=\boldsymbol{\alpha}^{(t)}+\Delta\alpha_{j}\boldsymbol{e}_{j}$

9: $\boldsymbol{w}^{(t+1)}=\boldsymbol{w}^{(t)}+(\lambda
n)^{-1}\Delta\alpha_{j}\boldsymbol{a}_{j}$

10:end for

Algorithm 11 Stochastic Dual Coordinate Descent (adaptive)

1:let $\boldsymbol{\alpha}^{(0)}=0,$ $\boldsymbol{w}^{(0)}=0$

2:for t = 0,1,… do

3: calculate absolute values of dual residuals $|\kappa_{j}^{(t)}|$ for all
$j\in[n]$

4: generate adaptive probabilities distribution $\boldsymbol{p}^{(t)}$:
$p_{i}^{(t)}=\frac{|\kappa_{i}^{(t)}|\|\boldsymbol{a}_{i}\|}{\sum_{j}|\kappa_{j}^{(t)}|\|\boldsymbol{a}_{j}\|}$

5: sample $j$ from $[n]$ according to distribution $\boldsymbol{p}^{(t)}$

6:
$\Delta\alpha_{j}=y_{j}\max\bigg{(}0,\min\Big{(}1,\frac{1-y_{j}\boldsymbol{a}_{j}^{\top}\boldsymbol{w}^{(t)}}{\|\boldsymbol{a}_{j}\|^{2}/(\lambda
n)}+y_{j}\alpha_{j}^{(t)}\Big{)}\bigg{)}-\alpha_{j}^{(t)}$

7:
$\boldsymbol{\alpha}^{(t+1)}=\boldsymbol{\alpha}^{(t)}+\Delta\alpha_{j}\boldsymbol{e}_{j}$

8: $\boldsymbol{w}^{(t+1)}=\boldsymbol{w}^{(t)}+(\lambda
n)^{-1}\Delta\alpha_{j}\boldsymbol{a}_{j}$

9:end for

Algorithm 12 Stochastic Dual Coordinate Descent (ada-uniform)

1:let $\boldsymbol{\alpha}^{(0)}=0,$ $\boldsymbol{w}^{(0)}=0$

2:for t = 0,1,… do

3: calculate absolute values of dual residuals $|\kappa_{j}^{(t)}|$ for all
$j\in[n]$

4: find $t$-support set $I_{t}=\\{i\in[n]:\kappa_{i}^{(t)}\neq
0\\}\subseteq[n]$

5: generate adaptive probabilities distribution $\boldsymbol{p}^{(t)}$:
$\begin{cases}p_{i}^{(t)}=\frac{1}{2|I_{t}|}+\frac{|\kappa_{i}^{(t)}|\|\boldsymbol{a}_{i}\|}{2\sum_{j}|\kappa_{j}^{(t)}|\|\boldsymbol{a}_{j}\|},&\text{if
}\kappa_{i}^{(t)}\neq 0\\\ p_{i}^{(t)}=0,&\text{otherwise}\end{cases}$

6: sample $j$ from $[n]$ according to distribution $\boldsymbol{p}^{(t)}$

7:
$\Delta\alpha_{j}=y_{j}\max\bigg{(}0,\min\Big{(}1,\frac{1-y_{j}\boldsymbol{a}_{j}^{\top}\boldsymbol{w}^{(t)}}{\|\boldsymbol{a}_{j}\|^{2}/(\lambda
n)}+y_{j}\alpha_{j}^{(t)}\Big{)}\bigg{)}-\alpha_{j}^{(t)}$

8:
$\boldsymbol{\alpha}^{(t+1)}=\boldsymbol{\alpha}^{(t)}+\Delta\alpha_{j}\boldsymbol{e}_{j}$

9: $\boldsymbol{w}^{(t+1)}=\boldsymbol{w}^{(t)}+(\lambda
n)^{-1}\Delta\alpha_{j}\boldsymbol{a}_{j}$

10:end for

Algorithm 13 Stochastic Dual Coordinate Descent (ada-gap)

1:let $\boldsymbol{\alpha}^{(0)}=0,$ $\boldsymbol{w}^{(0)}=0$

2:for t = 0,1,… do

3: generate adaptive probabilities distribution $\boldsymbol{p}^{(t)}$:
$p_{i}=\frac{\varphi_{i}(\boldsymbol{a}_{i}^{\top}\boldsymbol{w})-\alpha_{i}y_{i}+\alpha_{i}\boldsymbol{a}_{i}^{\top}\boldsymbol{w}}{\sum_{j=1}^{n}\Big{(}\varphi_{j}(\boldsymbol{a}_{j}^{\top}\boldsymbol{w})-\alpha_{j}y_{j}+\alpha_{j}\boldsymbol{a}_{j}^{\top}\boldsymbol{w}\Big{)}}$

4: sample $j$ from $[n]$ according to distribution $\boldsymbol{p}^{(t)}$

5:
$\Delta\alpha_{j}=y_{j}\max\bigg{(}0,\min\Big{(}1,\frac{1-y_{j}\boldsymbol{a}_{j}^{\top}\boldsymbol{w}^{(t)}}{\|\boldsymbol{a}_{j}\|^{2}/(\lambda
n)}+y_{j}\alpha_{j}^{(t)}\Big{)}\bigg{)}-\alpha_{j}^{(t)}$

6:
$\boldsymbol{\alpha}^{(t+1)}=\boldsymbol{\alpha}^{(t)}+\Delta\alpha_{j}\boldsymbol{e}_{j}$

7: $\boldsymbol{w}^{(t+1)}=\boldsymbol{w}^{(t)}+(\lambda
n)^{-1}\Delta\alpha_{j}\boldsymbol{a}_{j}$

8:end for

##  Appendix C Experimental results on large dataset

A comparison of the algorithms with fixed per-epoch sampling on large dataset
rcv1 (see Table 3) is given in Figures 3 and 4. We use $\lambda=7\cdot
10^{-4}$ for Lasso and $\lambda=0.1$ for SVM.

Table 3: Datasets Dataset | $d$ | $n$ | nnz$/(nd)$ | $c_{v}=\frac{\mu(\|\boldsymbol{a}_{i}\|)}{\sigma(\|\boldsymbol{a}_{i}\|)}$
---|---|---|---|---
rcv1 | 47236 | 20242 | 0.16% | 0.57

(a) Gap

(b) Suboptimality

Figure 3: Lasso on the full rcv1 dataset. Performance in terms of duality gap
and suboptimality

(a) Gap

(b) Suboptimality

Figure 4: SVM on the full rcv1 dataset. Performance in terms of duality gap
and suboptimality