Asynchronous Stochastic Coordinate Descent: Parallelism and Convergence Properties

Introduction

We consider the convex optimization problem

Formulations of the type (1) arise in many data analysis and machine learning problems, for example, the linear primal or nonlinear dual formulation of support vector machines , the LASSO approach to regularized least squares, and regularized logistic regression. Algorithms based on gradient and approximate / partial gradient information have proved effective in these settings. We mention in particular gradient projection and its accelerated variants , proximal gradient and accelerated proximal gradient methods for regularized objectives, and stochastic gradient methods . These methods are inherently serial, in that each iteration depends on the result of the previous iteration. Recently, parallel multicore versions of stochastic gradient and stochastic coordinate descent have been described for problems involving large data sets; see for example .

This paper proposes an asynchronous stochastic proximal coordinate-descent algorithm, called AsySPCD, for composite objective functions. The basic step of AsySPCD, executed repeatedly by each core of a multicore system, is as follows: Choose an index $i\in\{1,2,\dotsc,n\}$; read $x$ from shared memory and evaluate the $i$th element of $\nabla f$; subtract a short, constant, positive of this partial gradient from $(x)_{i}$; and perform a proximal operation on $(x)_{i}$ to account for the regularization term $g_{i}(\cdot)$. We use a simple model of computation that matches well to modern multicore architectures. Each core performs its updates on centrally stored vector $x$ in an asynchronous, uncoordinated fashion, without any form of locking. A consequence of this model is that the version of $x$ that is read by a core in order to evaluate its gradient is usually not the same as the version to which the update is made later, because $x$ is updated in the interim by other cores. (Generally, we denote by $\hat{x}$ the version of $x$ that is used by a core to evaluate its component of $\nabla f(\hat{x})$.) We assume, however, that indefinite delays do not occur between reading and updating: There is a bound $\tau$ such no more than $\tau$ component-wise updates to $x$ are missed by a core, between the time at which it reads the vector $\hat{x}$ and the time at which it makes its update to the chosen element of $x$. A similar model of parallel asynchronous computation was used in Hogwild! and AsySCD . However, there is a key difference in this paper: We do not assume that the evaluation vector $\hat{x}$ is a version of $x$ that actually existed in the shared memory at some point in time. Rather, we account for the fact that the components of $x$ may be updated by multiple cores while in the process of being read by another core, so that $\hat{x}$ may be a “hybrid” version that never actually existed in memory. Our new model, which we call an “inconsistent read” model, is significantly closer to the reality of asynchronous computation, and dispenses with the somewhat unsatisfying “consistent read” assumption of previous work. It also requires a quite distinct style of analysis; our proofs differ substantially from those in previous related works.

We show that, for suitable choices of steplength, our algorithm converges at a linear rate if an “optimal strong convexity” property (2) holds. It attains sublinear convergence at a “$1/k$” rate for general convex functions. Our analysis also defines a sufficient condition for near-linear speedup in the number of cores used. This condition relates the value of delay parameter $\tau$ (which corresponds closely to the number of cores / threads used in the computation) to the problem dimension $n$. A parameter that quantifies the cross-coordinate interactions in $\nabla f$ also appears in this relationship. When the Hessian of $f$ is nearly diagonal, the minimization problem (1) is almost separable, so higher degrees of parallelism are possible.

We review related work in Section 2. Section 3 specifies the proposed algorithm. Convergence results are described in Section 4, with proofs given in the appendix. Computational experience is reported in Section 5. A summary and conclusions appear in Section 6.

We use the following notation in the remainder of the paper.

$\Omega$ denotes the intersection of ${\rm dom}(f)$ and ${\rm dom}(g)$

$S$ denotes the set on which $F$ attains its optimal value, which is denoted by $F^{*}$.

$\mathcal{P}_{S}(\cdot)$ denotes Euclidean-norm projection onto $S$.

Given a matrix $A$, we use $A_{\cdot j}$ to denote its $j$th column and $A_{i\cdot}$ to denote its $i$th row.

$\|\cdot\|$ denotes the Euclidean norm $\|\cdot\|_{2}$.

$f^{*}_{j}:=f(\mathcal{P}_{S}(x_{j}))$ and $g^{*}_{j}:=g(\mathcal{P}_{S}(x_{j}))$.

$F^{*}:=F(\mathcal{P}_{S}(x))$ denotes the optimal objective value. (Note that $F^{*}=f^{*}_{j}+g^{*}_{j}$ for any $j$.)

We use $(x)_{i}$ for the $i$th element of $x$, and $\nabla_{i}f(x)$ for the $i$th element of $\nabla f(x)$.

Similarly, for the vector function $g$, we denote

Note that the proximal operator is nonexpansive, that is, $\|\mathcal{P}_{g}(x)-\mathcal{P}_{g}(y)\|\leq\|x-y\|$.

We define the following optimal strong convexity condition for a convex function $f$ with respect to the optimal set $S$, with parameter $l>0$:

This condition is significantly weaker than the usual strong convexity condition; a strongly convex function $F(.)$ is an optimally strongly convex function, but the converse is not true in general. We provide several examples of optimally strongly convex functions that are not strongly convex:

$F(x)=f(Ax)$, where $f$ is a strongly convex function and $A$ is any matrix, possibly one with a nontrivial kernel.

$F(x)=f(Ax)+{\bf 1}_{X}(x)$ with strongly convex $f$, and arbitrary $A$, where ${\bf 1}_{X}(x)$ is an indicator function defined on a polyhedron set $X$. Note first that $y^{*}:=Ax^{*}$ is unique for any $x^{*}\in S$, from the strong convexity of $f$. The optimal solution set $S$ is defined by

The inequality (2) clearly holds for $x\notin X$, since the left-hand side is infinite in this case. For $x\in X$, we have by the famous theorem of Hoffman that there exists $c>0$ such that

Then from the strong convexity of $f(x)$, we have that there exists a positive number $l$ such that for any $x\in X$

Squared hinge loss $F(x)=\sum_{i}\max(0,a^{T}_{i}x-y_{i})^{2}$. To verify optimal strong convexity, we reformulate this problem as

Note that optimal strong convexity (2) is a weaker version of the “essential strong convexity” condition used in . A concept called “restricted strong convexity” proposed in (See Lemma 4.6) is similar in that it requires a certain quantity to increase quadratically with distance from the solution set, but different in that the objective is assumed to be differentiable. Anitescu defines a “quadratic growth condition” for (smooth) nonlinear programming in which the objective is assumed to grow at least quadratically with distance to a local solution in some feasible neighborhood of that solution. Since our setting (unconstrained, nonsmooth, convex) is quite different, we believe the use of a different term is warranted here.

Throughout this paper, we make the following assumption.

Lipschitz Constants

The coordinate Lipschitz constant $L_{\mbox{\rm\scriptsize max}}$ is defined for $x$, $i$, $t$ satisfying the same conditions as above:

We denote the ratio between these two quantities by $\Lambda$:

Besides bounding the nonlinearity of $f$ along various directions, the quantities $L_{\mbox{\rm\scriptsize res}}$ and $L_{\mbox{\rm\scriptsize max}}$ capture the interactions between the various components in the gradient $\nabla f$. In the case of twice continuously differentiable $f$, we can understand these interactions by observing the diagonal and off-diagonal terms of the Hessian $\nabla^{2}f(x)$. Let us consider upper bounds on the ratio $\Lambda$ in various situations. For simplicity, we suppose that $f$ is quadratic with positive semidefinite Hessian $Q$.

If $Q$ is sparse with at most $p$ nonzeros per row/column, we have that

so that $\Lambda\leq\sqrt{p}$ in this situation.

If $Q$ is diagonally dominant, we have for any column $i$ that

which, by taking the maximum of both sides, implies that $\Lambda\leq 2$ in this case.

Related Work

We have surveyed related work on coordinate descent and stochastic gradient methods in a recent report . Our discussion there included non-stochastic, cyclic coordinate descent methods , synchronous parallel methods that distribute the work of function and gradient evaluation , and asynchronous parallel stochastic gradient methods (including the randomized Kaczmarz algorithm) . We make some additional comments here on related topics, and include some recent references from this active research area.

Stochastic coordinate descent can be viewed as a special case of stochastic gradient, so analysis of the latter approach can be applied, to obtain for example a sublinear $1/k$ rate of convergence in expectation for strongly convex functions; see, for example . However, stochastic coordinate descent is “special” in that it is possible to guarantee improvement in the objective at every step. Nesterov studied the convergence rate for a stochastic block coordinate descent method for unconstrained and separably constrained convex smooth optimization, proving linear convergence for the strongly convex case and a sublinear $1/k$ rate for the convex case. Richtárik and Takáč and Lu and Xiao extended this work to composite minimization, in which the objective is the sum of a smooth convex function and a separable nonsmooth convex function, and obtained similar (slightly stronger) convergence results. Stochastic coordinate descent is extended by Necoara and Patrascu to convex optimization with a single linear constraint, randomly updating two coordinates at a time to maintain feasibility.

In the class of synchronous parallel methods for coordinate descent, Richtárik and Takáč studied a synchronized parallel block (or minibatch) coordinate descent algorithm for composite optimization problems of the form (1), with a block separable regularizer $g$. At each iteration, processors update the randomly selected coordinates concurrently and synchronously. Speedup depends on the sparsity of the data matrix that defines the loss functions. A similar synchronous parallel method was studied in and ; the latter focuses on the case of $g(x)=\|x\|_{1}$. Scherrer et al. make greedy choices of multiple blocks of variables to update in parallel. Another greedy way of selecting coordinates was considered by Peng et al. , who also describe a parallel implementation of FISTA, an accelerated first-order algorithm due to Beck and Teboulle . Fercoq and Richtárik consider a variant of (1) in which $f$ is allowed to be nonsmooth. They apply Nesterov’s smoothing scheme to obtain a smoothed version and update multiple blocks of coordinates using block coordinate descent in parallel. Sublinear convergence rate is established for both strongly convex and weakly convex cases. Fercoq and Richtárik proposed a variant of Nesterov’s accelerated scheme to accelerate the synchronous parallel block coordinate algorithm of , proving an improved sublinear convergence rate for weakly convex problems. This variant avoids the disadvantage of the original Nesterov acceleration scheme , which requires $O(n)$ complexity per iteration, even on sparse data. Facchinei, Sagratella, and Scutari consider a general framework for synchronous block coordinate descent methods with separable regularizers, in which the block subproblems may be solved inexactly. However, the block to be updated at each step is not chosen randomly; it must contain a component that is furthest from optimality, in some sense.

We turn now to asynchronous parallel methods. Bertsekas and Tsitsiklis described an asynchronous method for fixed-point problems $x=q(x)$ over a separable convex closed feasible region. (The optimization problem (1) can be formulated in this way by defining $q(x):=\mathcal{P}_{\alpha g}[(I-\alpha\nabla f)(x)]$ for a fixed $\alpha>0$.) They use an inconsistent-read model of asynchronous computation, and establish linear convergence provided that components are not neglected indefinitely and that the iteration $x=q(x)$ is a maximum-norm contraction. The latter condition is quite strong. In the case of $g$ null and $f$ convex quadratic in (1) for instance, it requires the Hessian to satisfy a diagonal dominance condition — a stronger condition than strong convexity. By comparison, AsySCD guarantees linear convergence under an “essential strong convexity” condition, though it assumes a consistent-read model of asynchronous computation. Elsner et al. considered the same fixed point problem and architecture as , and describe a similar scheme. Their scheme appears to require locking of the shared-memory data structure for $x$ to ensure consistent reading and writing. Frommer and Szyld give a comprehensive survey of asynchronous methods for solving fixed-point problems.

Liu et al. followed the asynchronous consistent-read model of Hogwild! to develop an asynchronous stochastic coordinate descent (AsySCD) algorithm and proved sublinear ($1/k$) convergence on general convex functions and a linear convergence rate on functions that satisfy an “essential strong convexity” property. Sridhar et al. developed an efficient LP solver by relaxing an LP problem into a bound-constrained QP problem, which is then solved by AsySCD.

Liu et al. developed an asynchronous parallel variant of the randomized Kaczmarz algorithm for solving a general consistent linear system $Ax=b$, proving a linear convergence rate. Avron et al. proposed an asynchronous solver for the system $Qx=c$ where $Q$ is a symmetric positive definite matrix, proving a linear convergence rate. This method is essentially an asynchronous stochastic coordinate descent method applied to the strongly convex quadratic optimization problem $\min_{x}\,{1\over 2}x^{T}Qx-c^{T}x$. The paper considers both inconsistent- and consistent-read cases are considered, with slightly different convergence results.

Algorithm

In our algorithm AsySPCD, multiple processors have access to a shared data structure for the vector $x$, and each processor is able to compute a randomly chosen element of the gradient vector $\nabla f(x)$. Each processor repeatedly runs the following proximal coordinate descent process. (Choice of the steplength parameter $\gamma$ is discussed further in the next section.)

Choose an index $i\in\{1,2,\dotsc,n\}$ at random, read $x$ into the local storage location $\hat{x}$, and evaluate $\nabla_{i}f(\hat{x})$;

Update component $i$ of the shared $x$ by taking a step of length $\gamma/L_{\mbox{\rm\scriptsize max}}$ in the direction $-\nabla_{i}f(\hat{x})$, follows by a proximal operation defined as follows:Our analysis assumes that no other process modifies $x_{i}$ while this proximal operation is being computed. As we explain in Section 5, our practical implementation actually assigns each coordinate $x_{i}$ to a single core, and allows only that core to update $x_{i}$, so this issue does not arise. An alternative implementation, pointed out by a referee, would be to use a “compare-and-swap” atomic instruction to implement the update. This operation would perform the update only if $x_{i}$ was not changed while the update was being computed.

Notice that each step changes just a single element of $x$, that is, the $i$th element. Unlike standard proximal coordinate descent, the value $\hat{x}$ at which the coordinate gradient is calculated usually differs from the value of $x$ to which the update is applied, because while the processor is evaluating its gradient, other processors may repeatedly update the value of $x$ stored in memory. As mentioned above, we use an “inconsistent read” model of asynchronous computation here, in contrast to the “consistent read” models of AsySCD and Hogwild! . Figure 1 shows how inconsistent reading can occur, as a result of updating of components of $x$ while it is being read. Consistent reading can be guaranteed by means of a software lock, but such a mechanism degrades parallel performance significantly. In fact, the implementations of Hogwild! and AsySCD described in the papers do not use any software lock, and in this respect the computations in those papers are not quite compatible with their analysis.

The “global” view of algorithm AsySPCD is shown in Algorithm 1. To obtain this version from the “local” version, we introduce a counter $j$ to track the total number of updates applied to $x$, so that $x_{j}$ is the state of $x$ in memory after update $j$ is performed. We use $i(j)$ to denote the component that is updated at iteration $j$, and $\hat{x}_{j}$ for value of $x$ that is used in the calculation of the gradient element $\nabla f_{i(j)}$. The components of $\hat{x}_{j}$ may have different ages. Some components may be current at iteration $j$, others may not reflect recent updates made by other processors. We assume however that there is an upper bound of $\tau$ on the age of each component, measured in terms of updates. $K(j)$ defines an iterate set such that

One can see that $d\leq j-1$, $\forall d\in K(j)$. Here we assume $\tau$ to be the upper bound on the age of all elements in $K(j)$, for all $j$, so that $\tau\geq j-\min\{d~{}|~{}d\in K(j)\}$. We assume further that $K(j)$ is ordered from oldest to newest index (that is, smallest to largest). Note that $K(j)$ is empty if $x_{j}=\hat{x}_{j}$, that is, if the step is simply an ordinary stochastic coordinate gradient update. The value of $\tau$ corresponds closely to the number of cores involved in the computation provided that computation of the update for each component of $x$ costs roughly the same.

Main Results

This section presents results on convergence of AsySPCD. The theorem encompasses both the linear rate for optimally strongly convex $f$ and the sublinear rate for general convex $f$. The result depends strongly on the delay parameter $\tau$. The proofs are highly technical, and are relegated to Appendix A. We note the proof techniques differ significantly from those used for the consistent-read algorithms of and .

We start by describing the key idea of the algorithm, which is reflected in the way that it chooses the steplength parameter $\gamma$. Denoting $\bar{x}_{j+1}$ by

so that $x_{j+1}-x_{j}=[(\bar{x}_{j+1})_{i(j)}-(x_{j})_{i(j)}]e_{i(j)}$. Thus, we have

Therefore, we can view $\bar{x}_{j+1}-x_{j}$ as capturing the expected behavior of $x_{j+1}-x_{j}$. Note that when $g(x)=0$, we have $\bar{x}_{j+1}-x_{j}=-({\gamma}/{L_{\mbox{\rm\scriptsize max}}})\nabla f(\hat{x}_{j})$, a standard negative-gradient step. The choice of steplength parameter $\gamma$ entails a tradeoff: We would like $\gamma$ to be long enough that significant progress is made at each step, but not so long that the gradient information computed at $\hat{x}_{j}$ is stale and irrelevant by the time the update is applied to $x_{j}$. We enforce this tradeoff by means of a bound on the ratio of expected squared norms on $x_{j}-\bar{x}_{j+1}$ at successive iterates; specifically,

where $\rho>1$ is a user defined parameter. The analysis becomes a delicate balancing act in the choice of $\rho$ and steplength $\gamma$ between aggression and excessive conservatism. We find, however, that these values can be chosen to ensure steady convergence for the asynchronous method at a linear rate, with rate constants that are almost consistent with a standard short-step proximal full-gradient descent, when the optimal strong convexity condition (2) is satisfied.

Our main convergence result is the following.

Suppose that Assumption 1 is satisfied. Let $\rho$ be a constant that satisfies $\rho>1+4/\sqrt{n}$, and define the quantities $\theta$, $\theta^{\prime}$, and $\psi$ as follows:

Suppose that the steplength parameter $\gamma>0$ satisfies the following two bounds:

If the optimal strong convexity property (2) holds with $l>0$, we have for $j=1,2,\dotsc$ that

while for general smooth convex function $f$, we have

The following corollary proposes an interesting particular choice for the parameters for which the convergence expressions become more comprehensible. The result requires a condition on the delay bound $\tau$ in terms of $n$ and the ratio $\Lambda$.

then the steplength $\gamma=1/2$ will satisfy the bounds (9). In addition, when the optimal strong convexity property (2) holds with $l>0$, we have for $j=1,2,\dotsc$ that

while for the case of general convex $f$, we have

We note that the linear rate (15) is broadly consistent with the linear rate for the classical steepest descent method applied to strongly convex functions, which has a rate constant of $(1-2l/L)$, where $L$ is the standard Lipschitz constant for $\nabla f$. Suppose we assume (not unreasonably) that $n$ steps of stochastic coordinate descent cost roughly the same as one step of steepest descent, and that $l\leq L_{\mbox{\rm\scriptsize max}}$. It follows from (15) that $n$ steps of stochastic coordinate descent would achieve a reduction factor of about

so a standard argument would suggest that stochastic coordinate descent would require about $6L_{\mbox{\rm\scriptsize max}}/L$ times more computation. Since $L_{\mbox{\rm\scriptsize max}}/L\in[1/n,1]$, the stochastic asynchronous approach may actually require less computation. It may also gain an advantage from the parallel asynchronous implementation. A parallel implementation of standard gradient descent would require synchronization and careful division of the work of evaluating $\nabla f$, whereas the stochastic approach can be implemented in an asynchronous fashion.

For the general convex case, (16) defines a sublinear rate, whose relationship with the rate of standard gradient descent for general convex optimization is similar to the previous paragraph.

Note that the results in Theorem 1 and Corollary 2 are consistent with the analysis for constrained AsySCD in , but this paper considers the more general case of composite optimization and the inconsistent-read model of parallel computation.

As noted in Section 1, the parameter $\tau$ corresponds closely to the number of cores that can be involved in the computation, since if all cores are working at the same rate, we would expect each other core to make one update between the times at which $x$ is read and (later) updated. If $\tau$ is small enough that (13) holds, the analysis indicates that near-linear speedup in the number of processors is achievable. A small value for the ratio $\Lambda$ (not much greater than $1$) implies a greater degree of potential parallelism. As we note at the end of Section 1, this ratio tends to closer to $1$ than to $\sqrt{n}$ in some important applications. In these situations, the bound (13) indicates that $\tau$ can vary like $n^{1/4}$ without affecting the iteration-wise convergence rate, and yielding near-linear speedup in the number of cores. This quantity is consistent with the analysis for constrained AsySCD in but weaker than the unconstrained AsySCD (which allows the maximal number of cores being $O(n^{1/2})$). A further comparison is with the asynchronous randomized Kaczmarz algorithm which allows $O(m)$ cores to be used efficiently when solving a consistent sparse linear system.

We conclude this section with a high-probability bound. The result follows immediately from Markov’s inequality. See Theorem 3 in for a related result and complete proof.

Suppose that the conditions of Corollary 2 hold, including the choice of $\rho$. Then for $\epsilon>0$ and $\eta\in(0,1)$, we have that

provided that one of the following conditions holds. In the optimally strongly convex case (2) with $l>0$, we require

iterations, while in the general convex case, it suffices that

Experiments

This section presents some results to illustrate the effectiveness of AsySPCD, in particular, the fact that near-linear speedup can be observed on a multicore machine. We note that more comprehensive experiments can be found in and , for unconstrained and box-constrained problems. Although the analysis in assumes consistent read, it is not enforced in the implementation, so apart from the fact that we now include a prox-step to account for the regularization term, the implementations in and are quite similar to the one employed in this section.

We choose $\sigma=0.01$ with $m=6000$, $n=10000$, and $s=10$ in Figure 2 and $m=12000$, $n=20000$, and $s=20$ in Figure 3. We set $\lambda=20\sqrt{m\log(n)}\sigma$ (a value of the order of $\sqrt{m\log(n)}\sigma$ is suggested by compressed sensing theory) and the steplength $\gamma$ is set as $1$ in both figures. In both cases, we can estimate the ratio $\Lambda=L_{\mbox{\rm\scriptsize res}}/L_{\mbox{\rm\scriptsize max}}$ roughly by $1+\sqrt{n/m}\approx 2.3$, as suggested at the end of Section 1. Our final computed values of $x$ have nonzeros in the same locations as the chosen solution $x^{*}$, though the values differ, because of the noise in $b$.

The left-hand graph in each figure indicates the number of threads / cores and plots objective function value vs epoch count, where one epoch is equivalent to $n$ iterations. Note that the curves are almost overlaid, indicating that the total workload required for AsySPCD is almost independent of the number of cores used in the computation. This observation validates our result in Corollary 2, which indicates that provided $\tau$ is below a certain threshold, it does not seriously affect the rate of convergence, as a function of total computation performed. The right-hand graph in each figure shows speedup when executed on different numbers of cores. Near-linear speedup is observed in Figure 3, while there is a slight dropoff for the larger numbers of cores in Figure 2. The difference can be explain by the smaller dimension of the problem illustrated in Figure 2. Referring to our threshold value (13) that indicates dimensions above which linear speedup should be expected, we have by setting $\Lambda\approx 2.3$ (as discussed above) and $\tau=10$ (the maximum number of threads used in this experiment) that the left-hand side of (13) is approximately 3000, while the right-hand side is $100$ (for Figure 2) and approximately $141$ (for Figure 2). As expected, our analysis is quite conservative; near-linear speedup is observed even when the threshold (13) is violated significantly.

Conclusions

This paper proposes an asynchronous parallel proximal stochastic coordinate descent algorithm for minimizing composite objectives of the form (1). Sublinear convergence (at rate $1/k$) is proved for general convex functions, with stronger linear convergence results for problems that satisfy the optimal strong convexity property (2). Our analysis indicates the extent to which parallel implementations can be expected to yield near-linear speedup, in terms of a parameter that quantifies the cross-coordinate interactions in the gradient $\nabla f$ and a parameter $\tau$ that bounds the delay in updating. Our computational experience confirms that the linear speedup properties suggested by the analysis can be observed in practice.

Acknowledgments

The authors thank the editor and both referees for their valuable comments. Special thanks to Dr. Yijun Huang for her implementation of AsySPCD, which was used here to obtain computational results.

Appendix A Proofs of Main Results

This section provides the proofs for the main convergence results. We start with some preliminaries, then proceed to proofs of Theorem 1 and Corollary 2.

and formulate the update in Step 4 of Algorithm 1 in the following way:

(Note that $(x_{j+1})_{i}=(x_{j})_{i}$ for $i\neq i(j)$.) From the optimality condition for this formulation (see (41) in ), we have for all $x$ that

By rearranging this expression and substituting $\mathcal{P}_{S}(x)$ for $x$, we find that the following inequality is true for all $x$:

From the definition of $L_{\mbox{\rm\scriptsize max}}$, and using the notation (18), we have

From the definition of $\bar{x}_{j+1}$ in (5), we have

and note that this definition is consistent with $(\Delta_{j})_{i(j)}$ defined in (18). From (6), we have

Recalling that the indices in $K(j)$ are sorted in the increasing order from smallest (oldest) iterate to largest (newest) iterate, we use $K(j)_{t}$ to denote the $t$-th smallest entry in $K(j)$. For $T=0,1,\dotsc,|K(j)|$, we define

where the second inequality holds because $\hat{x}_{j,t}$ and $\hat{x}_{j,t+1}$ differ in only a single coordinate.

A.2 Proof of Theorem 1

We prove (10) by induction. First, note that for any vectors $a$ and $b$, we have

The second factor in the r.h.s. of (25) is bounded as follows:

where the fourth inequality uses $\|\nabla f(x_{j})-\nabla f(x_{j-1})\|\leq L_{\mbox{\rm\scriptsize res}}\|x_{j}-x_{j-1}\|$, since $x_{j}$ and $x_{j-1}$ differ in just one component.

We set $j=1$, and note that $K(0)=\emptyset$ and $K(1)\subset\{0\}$. In this case, we obtain a bound from (26)

By substituting this bound in (25) and setting $j=1$, and taking expectations, we obtain

For any positive scalars $\mu_{1}$, $\mu_{2}$, and $\alpha$, we have

By taking $j=1$ in (30), and substituting in (28), we obtain

To see the last inequality, one only needs to verify that

where the last inequality follows from the second bound for $\gamma$ in (9). We have thus shown that (10) holds for $j=1$.

To take the inductive step, we assume that (10) holds up to index $j-1$. We have for $j-1-\tau\leq d\leq j-2$ and any $\beta>0$ (using (29) again) that

Thus by setting $\beta=\rho^{(d+1-j)/2}$, we obtain

By substituting (27) into (25) and taking expectation on both sides of (25), we obtain

where the last equality follows from the definition of $\theta$ in (8). It follows that

To see the last inequality, one only needs to verify that

and the last inequality is true because of the upper bound of $\gamma$ in (9). We have thus proved (10).

Next we will show the expectation of the objective $F$ is monotonically decreasing. We have by using the definition (18) and (6) that

Consider the expectation of the last term on the right-hand side of this expression. We have

By taking expectations on both sides of (32) and substituting (33), we obtain

To see $\left({\frac{1}{\gamma}-\frac{1}{2}}\right)L_{\mbox{\rm\scriptsize max}}-{L_{\mbox{\rm\scriptsize res}}\theta\over n^{1/2}}\geq 0$ or equivalently $\left({\frac{1}{\gamma}-\frac{1}{2}}\right)-{\Lambda\theta\over n^{1/2}}\geq 0$, we note from (8) and (9) that

Therefore, we have proved the monotonicity of the expectation of the objectives, that is,

Next we prove the sublinear convergence rate for the constrained smooth convex case in (12). We have

For the expectation of $T_{1}$, defined in (35), we have

For $T_{3}$, let us look the expectation of several individual terms first

Now we take the expectation on $T_{3}$ and use the equalities above to obtain:

By substituting the upper bounds from (37), (38), and (39) into (35), we obtain

which follows from the definition (8) of $\psi$ and from the first upper bound on $\gamma$ in (9). It follows from (40) that

By substituting the definition of $S_{j+1}$ into (44), we obtain

The sublinear convergence expression (12) follows when we drop the (nonnegative) first term on the left-hand side of this expression, and rearrange.

Finally, we prove the linear convergence rate (11) for the optimally strongly convex case. All bounds proven above continue to hold, and we make use the optimal strong convexity property in (2):

By using this result together with some elementary manipulation, we obtain

By taking expectations of both sides in this expression, and comparing with (42), we obtain

where the last inequality follows from induction over $j$. We obtain (11) by substituting the definition (42) of $S_{j}$. ∎

A.3 Proof of Corollary 2

Note that for $\rho$ defined by (14), and using (13), we have

Thus from the definition of $\psi$ (8), we have that

where for the second-last inequality we used (13) to obtain

Thus, the steplength parameter choice $\gamma=1/2$ satisfies the first bound in (9). To show that the second bound in (9) holds also, we have

We can thus set $\gamma=1/2$, and by substituting this choice into (11), we obtain (15). We obtain (16) by making the same substitution into (12). ∎

References