Semi-Stochastic Coordinate Descent

Introduction

In this paper we study the problem of unconstrained minimization of a strongly convex function represented as the average of a large number of smooth convex functions:

Many computational problems in various disciplines are of this form. In machine learning, $f_{i}(x)$ represents the loss/risk of classifier $x\in\R^{d}$ on data sample $i$, $f$ represents the empirical risk (=average loss), and the goal is to find a predictor minimizing $f$. Often, an L2-regularizer of the form $\mu\|x\|^{2}$, for $\mu>0$, is added to the loss, making it strongly convex and hence easier to minimize.

We assume that the functions $f_{i}:\R^{d}\rightarrow\R$ are differentiable and convex function, with Lipschitz continuous partial derivatives. Formally, we assume that for each $i\in[n]:=\{1,2,\dots,n\}$ and $j\in[d]:=\{1,2,\dots,d\}$ there exists $L_{ij}\geq 0$ such that for all $x\in\R^{d}$ and $h\in\R$,

where $e_{j}$ is the $j^{th}$ standard basis vector in $\R^{d}$, $\nabla f(x)\in\R^{d}$ the gradient of $f$ at point $x$ and $\langle\cdot,\cdot\rangle$ is the standard inner product. This assumption was recently used in the analysis the accelerated coordinate descent method APPROX . We further assume that $f$ is $\mu$-strongly convex. That is, we assume that there exists $\mu>0$ such that for all $x,y\in\R^{d}$,

Context.

Batch methods such as gradient descent (GD) enjoy a fast (linear) convergence rate: to achieve $\epsilon$-accuracy, GD needs $\mathcal{O}(\kappa\log(1/\epsilon))$ iterations, where $\kappa$ is a condition number. The drawback of GD is that in each iteration one needs to compute the gradient of $f$, which requires a pass through the entire dataset. This is prohibitive to do many times if $n$ is very large.

Stochastic gradient descent (SGD) in each iteration computes the gradient of a single randomly chosen function $f_{i}$ only—this constitutes an unbiased (but noisy) estimate of the gradient of $f$—and makes a step in that direction . The rate of convergence of SGD is slower, $\mathcal{O}(1/\epsilon)$, but the cost of each iteration is independent of $n$. Variants with nonuniform selection probabilities were considered in , a mini-batch variant (for SVMs with hinge loss) was analyzed in .

Recently, there has been progress in designing algorithms that achieve the fast $O(\log(1/\epsilon))$ rate without the need to scan the entire dataset in each iteration. The first class of methods to have achieved this are stochastic/randomized coordinate descent methods.

When applied to (11), coordinate descent methods (CD) can, like SGD, be seen as an attempt to keep the benefits of GD (fast linear convergence) while reducing the complexity of each iteration. A CD method only computes a single partial derivative $\nabla_{j}f(x)$ at each iteration and updates a single coordinate of vector $x$ only. When chosen uniformly at random, partial derivative is also an unbiased estimate of the gradient. However, unlike the SGD estimate, its variance is low. Indeed, as one approaches the optimum, partial derivatives decrease to zero. While CD methods are able to obtain linear convergence, they typically need $O((d/\mu)\log(1/\epsilon))$ iterations when applied to (11) directly The complexity can be improved to $O(\tfrac{d\beta}{\tau\mu}\log(1/\epsilon))$ in the case when $\tau$ coordinates are updated in each iteration, where $\beta\in[1,\tau]$ is a problem-dependent constant . This has been further studied for nonsmooth problems via smoothing , for arbitrary nonuniform distributions governing the selection of coordinates and in the distributed setting . Also, efficient accelerated variants with $O(1/\sqrt{\epsilon})$ rate were developed , capable of solving problems with 50 billion variables.. CD method typically significantly outperform GD, especially on sparse problems with a very large number variables/coordinates .

An alternative to applying CD to (11) is to apply it to the dual problem. This is possible under certain additional structural assumptions on the functions $f_{i}$. This is the strategy employed by stochastic dual coordinate ascent (SDCA) , whose rate is

The condition number $\kappa$ here is different (and larger) than the condition number appearing in the rate of GD. Despite this, this is a vast improvement on the rates achieved by both GD and SGD, and the method indeed typically performs much better in practice. Accelerated and mini-batch variants of SDCA have also been proposed. We refer the reader to QUARTZ for a general analysis involving the update of a random subset of dual coordinates, following an arbitrary distribution.

Recently, there has been progress in designing primal methods which match the fast rate of SDCA. Stochastic average gradient (SAG) , and more recently SAGA , move in a direction composed of old stochastic gradients. The semi-stochastic gradient descent (S2GD) and stochastic variance reduced gradient (SVRG) methods employ a different strategy: one first computes the gradient of $f$, followed by $O(\kappa)$ steps where only stochastic gradients are computed. These are used to estimate the change of the gradient, and it is this direction which combines the old gradient and the new stochastic gradient information which used in the update.

Main result.

In this work we develop a new method—semi-stochastic coordinate descent (S2CD)—for solving (11), enjoying a fast rate similar to methods such as SDCA, SAG, S2GD, SVRG, MISO, SAGA, mS2GD and QUARTZ. S2CD can be seen as a hybrid between S2GD and CD. In particular, the complexity of our method is the sum of two terms:

evaluations $\nabla f_{i}$ (that is, $\log(1/\epsilon)$ evaluations of the gradient of $f$) and

evaluations of $\nabla_{j}f_{i}$ for randomly chosen functions $f_{i}$ and randomly chosen coordinates $j$, where $\hat{\kappa}$ is a new condition number which is larger than the condition number $\kappa$ appearing in the aforementioned methods. However, note that $\hat{\kappa}$ enters the complexity only in the term involving the cost of the evaluation of a partial derivative $\nabla_{j}f_{i}$, which can be substantially smaller than the evaluation cost of $\nabla f_{i}$. Hence, our complexity result can be both better or worse than previous results, depending on whether the increase of the condition number can or can not be compensated by the lower cost of the stochastic steps based on the evaluation of partial derivatives.

Outline.

The paper is organized as follows. In Section 2 we describe the S2CD algorithm and in Section 3 we state a key lemma and our main complexity result. The proof of the lemma is provided in Section 4 and the proof of the main result in Section 5.

S2CD Algorithm

In this section we describe the Semi-Stochastic Coordinate Descent method (Algorithm 1).

The method has an outer loop (an “epoch”), indexed by counter $k$, and an inner loop, indexed by $t$. At the beginning of epoch $k$, we compute and store the gradient of $f$ at $x_{k}$. Subsequently, S2CD enters the inner loop in which a sequence of vectors $y_{k,t}$ for $t=0,1\dots,t_{k}$ is computed in a stochastic way, starting from $y_{k,0}=x_{k}$. The number $t_{k}$ of stochastic steps in the inner loop is random, following a geometric law:

In each step of the inner loop, we seek to compute $y_{k,t+1}$, given $y_{k,t}$. In order to do so, we sample coordinate $j$ with probability $p_{j}$ and subsequentlyIn S2CD, as presented, coordinates $j$ is selected first, and then function $i$ is selected, according to a distribution conditioned on the choice of $j$. However, one could equivalently sample $(i,j)$ with joint probability $p_{ij}$. We opted for the sequential sampling for clarity of presentation purposes. sample $i$ with probability $q_{ij}$, where the probabilities are given by

Note that $L_{ij}=0$ means that function $f_{i}$ does not depend on the $j^{th}$ coordinate of $x$. Hence, $\omega_{i}$ is the number of coordinates function $f_{i}$ depends on – a measure of sparsity of the dataThe quantity $\omega:=\max_{i}\omega_{i}$ (degree of partial separability of $f$) was used in the analysis of a large class of randomized parallel coordinate descent methods in . The more informative quantities $\{\omega_{i}\}$ appear in the analysis of parallel/distributed/mini-batch coordinate descent methods .. It can be shown that $f$ has a $1$-Lipschitz gradient with respect to the weighted Euclidean norm with weights $\{v_{j}\}$ ([3, Theorem 1]). Hence, we sample coordinate $j$ proportionally to this weight $v_{j}$. Note that $p_{ij}$ is the joint probability of choosing the pair $(i,j)$.

Having sampled coordinate $j$ and function index $i$, we compute two partial derivatives: $\nabla_{j}f_{i}(x_{k})$ and $\nabla_{j}f_{i}(y_{k,t})$ (we compressed the notation here by writing $\nabla_{j}f_{i}(x)$ instead of $\langle\nabla f_{i}(x),e_{j}\rangle$), and combine these with the pre-computed value $\nabla_{j}f(x_{k})$ to form an update of the form

Note that only a single coordinate of $y_{k,t}$ is updated at each iteration.

In the entire text (with the exception of the statement of Theorem 2 and a part of Section 5.3, where $\mathbf{E}$ denotes the total expectation) we will assume that all expectations are conditional on the entire history of the random variables generated up to the point when $y_{k,t}$ was computed. With this convention, it is possible to think that there are only two random variables: $j$ and $i$. By $\mathbf{E}$ we then mean the expectation with respect to both of these random variables, and by $\mathbf{E}_{i}$ we mean expectation with respect to $i$ (that is, conditional on $j$). With this convention, we can write

which means that conditioned on $j$, $G^{ij}_{kt}$ is an unbiased estimate of the $j^{th}$ partial derivative of $f$ at $y_{k,t}$. An equally easy calculation reveals that the random vector $g_{kl}^{ij}$ is an unbiased estimate of the gradient of $f$ at $y_{k,t}$:

Hence, the update step performed by S2CD is a stochastic gradient step of fixed stepsize $h$.

Before we describe our main complexity result in the next section, let us briefly comment on a few special cases of S2CD:

If $n=1$ (this can be always achieved simply by grouping all functions in the average into a single function), S2CD reduces to a stochastic CD algorithm with importance samplingA parallel CD method in which every subset of coordinates can be assigned a different probability of being chosen/updated was analyzed in ., as studied in , but written with many redundant computations. Indeed, the method in this case does not require the $x_{k}$ iterates, nor does it need to compute the gradient of $f$, and instead takes on the form:

It is possible to extend the S2CD algorithm and results to the case when coordinates are replaced by (nonoverlapping) blocks of coordinates, as in — we did not do it here for the sake of keeping the notation simple. In such a setting, we would obtain semi-stochastic block coordinate descent. In the special case with all variables forming a single block, the algorithm reduces to the S2GD method described in , but with nonuniform probabilities for the choice of $i$ — proportional to the Lipschitz constants of the gradient of the functions $f_{i}$ (this is also studied in ). As in , the complexity result then depends on the average of the Lipschitz constants.

Note that the algorithm, as presented, assumes the knowledge of $\mu$. We have done this for simplicity of exposition: the method works also if $\mu$ is replaced by some lower bound in the method, which can be set to 0 (see ). The change to the complexity results will be only minor and all our conclusions hold. Likewise, it is possible to give an $O(1/\epsilon)$ complexity result in the non-strongly convex case $f$ using standard regularization arguments (e.g., see ).

Complexity Result

In this section, we state and describe our complexity result; the proof is provided in Section 5.

An important step in our analysis is proving a good upper bound on the variance of the (unbiased) estimator $g_{kt}^{ij}=p_{j}^{-1}G_{kt}^{ij}e_{j}$ of $\nabla f(y_{k,t})$, one that we can “believe” would diminish to zero as the algorithm progresses. This is important for several reasons. First, as the method approaches the optimum, we wish $g_{kt}^{ij}$ to be progressively closer to the true gradient, which in turn will be close to zero. Indeed, if this was the case, then S2CD behaves like gradient descent with fixed stepsize $h$ close to optimum. In particular, this would indicate that using fixed stepsizes makes sense.

In light of the above discussion, the following lemma plays a key role in our analysis:

The iterates of the S2CD algorithm satisfy

The proof of this lemma can be found in Section 4.

Note that as $y_{k,t}\to x_{*}$ and $x_{k}\to x_{*}$, the bound (9) decreases to zero. This is the main feature of modern fast stochastic gradient methods: the squared norm of the stochastic gradient estimate progressively diminishes to zero, as the method progresses, in expectation. Note that the standard SGD method does not have this property: indeed, there is no reason for $\mathbf{E}_{i}\|\nabla f_{i}(x)\|^{2}$ to be small even if $x=x_{*}$.

We are now ready to state the main result of this paper.

If $0<h<1/(2\hat{L})$, then for all $k\geq 0$ we have: It is possible to replace modify the argument slightly and replace the term $\hat{L}$ appearing in the numerator by $\hat{L}-\frac{\mu}{\max_{s}p_{s}}$. However, as this does not bring any significant improvements, we decided to present the result in this simplified form.

By analyzing the above result (one can follow the steps in [6, Theorem 6]), we get the following useful corollary:

Fix the number of epochs $k\geq 1$, error tolerance $\epsilon\in(0,1)$ and let $\Delta:=\epsilon^{1/k}$ and $\hat{\kappa}:=\hat{L}/\mu$. If we run Algorithm 1 with stepsize $h$ and $m$ set as

then $\mathbf{E}[f(x_{k})-f(x_{*})]\leq\epsilon(f(x_{0})-f(x_{*}))$. In particular, for $k=\lceil\log(1/\epsilon)\rceil$ we have $\tfrac{1}{\Delta}\leq\exp(1)$, and we can pick

If we run S2CD with the parameters set as in (13), then in each epoch the gradient of $f$ is evaluated once (this is equivalent to $n$ evaluations of $\nabla f_{i}$), and the partial derivative of some function $f_{i}$ is evaluated $2m\approx 52\hat{\kappa}=O(\hat{\kappa})$ times. If we let $C_{grad}$ be the average cost of evaluating the gradient $\nabla f_{i}$ and $C_{pd}$ be the average cost of evaluating the partial derivative $\nabla_{j}f_{i}$, then the total work of S2CD can be written as

The complexity results of methods such as S2GD/SVRG and SAG/SAGA —in a similar but not identical setup to ours (these papers assume $f_{i}$ to be $L_{i}$-smooth)—can be written in a similar form:

where $\kappa=L/\mu$ and either $L=L_{max}:=\max_{i}L_{i}$ (), or $L=L_{avg}:=\tfrac{1}{n}\sum_{i}L_{i}$ (). The difference between our result (14) and existing results (15) is in the term $\hat{\kappa}{\cal C}_{pd}$ – previous results have $\kappa{\cal C}_{grad}$ in that place. This difference constitutes a trade-off: while $\hat{\kappa}\geq\kappa$ (we comment on this below), we clearly have ${\cal C}_{pd}\leq{\cal C}_{grad}$. The comparison of the quantities $\kappa{\cal C}_{grad}$ and $\hat{\kappa}{\cal C}_{pd}$ is not straightforward and is problem dependent.

Let us now comment how do the condition numbers $\hat{\kappa}$ and $\kappa_{avg}=L_{avg}/\mu$ compare. It can be show that (see )

and, moreover, this inequality can be tight. Since $\omega_{i}\geq 1$ for all $i$, we have

Finally, $\hat{\kappa}$ can be smaller or larger than $\kappa_{max}:=L_{max}/\mu$.

Proof of Lemma 1

We will prove the following stronger inequality:

Lemma 1 follows by dropping the negative term.

We first break down the left hand side of (16) into $d$ terms each of which we will bound separately. By first taking expectation conditioned on $j$ and then taking the full expectation, we can write:

STEP 2.

We now further break each of these $d$ terms into three pieces. That is, for each $j=1,\dots,d$ we have:

STEP 3.

In this step we bound the first two terms in the right hand side of inequality (18). It will now be useful to introduce the following notation:

Let us fist examine the first term in the right-hand side of (18). Using the coordinate co-coercivity lemma (Lemma 4) with $y=x_{*}$, we obtain the inequality

Note that by (4) and (10), we have that for all $s=1,2,\dots,d$,

Continuing from (21), we can therefore further write

The same reasoning applies to the second term on the right-hand side of the inequality (18) and we have:

STEP 4.

Next we bound the third term on the right-hand side of the inequality (18). First note that since $f$ is $\mu$-strongly convex (see (3)), for all $x\in\R^{d}$ we have:

STEP 5.

We conclude by combining (17), (18), (22), (23) and (25).

Proof of the Main Result

In this section we provide the proof of our main result. In order to present the proof in an organize fashion, we first establish two technical lemmas.

It is a well known and widely used fact (see, e.g. ) that for a continuously differentiable function $\phi:\R^{d}\to\R$ and constant $L_{\phi}>0$, the following two conditions are equivalent:

The second condition is often referred to by the name co-coercivity. Note that our assumption (2) on $f_{i}$ is similar to the first inequality. In our first lemma we establish a coordinate-based co-coercivity result which applies to functions $f_{i}$ satisfying (2).

For all $x,y\in\R^{d}$ and $i=1,\dots,n$, $j=1,\dots,d$, we have:

Fix any $i,j$ and $y\in\R^{d}$. Consider the function $g_{i}:\R^{d}\rightarrow\R$ defined by:

Then since $f_{i}$ is convex, we know that $g_{i}(x)\leq 0$ for all $x$, with $g_{i}(y)=0$. Hence, $y$ minimizes $g_{i}$. We also know that for any $x\in\R^{d}$:

Since $f_{i}$ satisfies (2), so does $g_{i}$, and hence for all $x\in\R^{d}$ and $h\in\R$, we have

which together with (27) yields the result. ∎

2 Recursion

We now proceed to the final lemma, establishing a key recursion which ultimately yields the proof of the main theorem, which we present in Section 5.3.

The iterates of S2CD satisfy the following recursion:

3 Proof of Theorem 2

where the expectation now is with respect to the entire history. Then by Lemma 5 we have the following $m$ inequalities:

By summing up the above $m$ inequalities, we get:

It follows from the strong convexity assumption (3) that

Hence if $0<2h\hat{L}<1$, then we obtain:

References