Coordinate Descent with Arbitrary Sampling I: Algorithms and Complexity

Introduction

With the dawn of the big data age, there has been a growing interest in solving optimization problems of unprecedented sizes. It was soon realized that traditional approaches, which work extremely well for problems of moderate sizes and when solutions of high accuracy are required, are not efficient for modern problems of large enough size and for applications where only rough or moderate accuracy solutions are sufficient. The focus of the optimization, numerical analysis and machine learning communities, and of practitioners in the sciences and industry, shifted to first-order (gradient) algorithms .

However, once the size of problems becomes truly big, it is necessary to turn to methods which are able to output a reasonably good solutions after an amount of work roughly equivalent to reading the data describing the problem a few times. For this to be possible, methods need to be able to progress while reading only a small part of the data describing the problem, which often means that a single iteration needs to be based on less information than that contained in the gradient of the objective (loss) function. The most popular methods of this type are stochastic gradient methods , randomized coordinate descent methods and semi-stochastic gradient descent methods .

In this paper we focus on randomized coordinate descent methods. After the seminal work of Nesterov , which provided an early theoretical justification of these methods for unconstrained convex minimization, the study has been successively extended to $L1$-regularized , proximal , parallel , distributed and primal-dual variants of coordinate descent. Accelerated coordinate decent—characterized by its $O(1/k^{2})$ complexity for non-strongly convex problems—was studied in . However, these methods are of theoretical nature only due to the fact that they rely on the need to perform full-dimensional vector operation at every iteration, which destroys the main advantage of coordinate descent – its ability to reduce the problem into subproblems of smaller sizes. A theoretically and practically efficient accelerated coordinate descent methods were proposed recently by Lee and Sidford and Fercoq and Richtárik , the latter work (APPROX algorithm) combining acceleration with parallelism and proximal setup. An accelerated distributed coordinate descent algorithm is obtained by specializing APPROX to a distributed sampling. All above mentioned papers only consider unconstrained or separably constrained problems. Some progress on linearly-coupled constraints has been made by Necoara et al in . Asynchronous variants of parallel coordinate descent methods were developed by Liu, Wright et al .

Virtually all existing work in stochastic optimization deals with a uniform sampling. In the context of coordinate descent, this means that the random subset (sampling) of coordinates chosen and updated at every iteration has the property that each coordinate is chosen equally likely. The possibility to assign different selection probabilities to different coordinates—also known as importance sampling—was considered in and recently in . However, these works consider the serial case only: a single coordinate is updated in each iteration. Randomized coordinate descent methods updating a subset of coordinates following an arbitrary distribution (i.e., using an arbitrary sampling) were first investigated by Richtárik and Takáč (NSync method) for strongly convex and smooth objective functions, and subsequently by Qu, Richtárik and Zhang (QUARTZ method), for strongly convex and possibly nonsmooth functions, and in a primal-dual framework.

In this paper we give the first fully unified analysis of gradient type algorithms which contain randomized coordinate descent on one end of the spectrum and gradient (or accelerated gradient) descent on the other hand. All our complexity results match or improve on the state of the art in all cases where specialized algorithms for specific samplings already exist. Moreover, we managed to substantially simplify the analysis for the sake of making the material accessible to a wide community.

2 Problem Formulation

In this paper we consider the composite optimization problem

3 Contributions

We now summarize the main contributions of this work.

We propose ALPHA (Algorithm 1) – a randomized gradient-type method for solving the convex composite optimization problem (1). In each iteration, ALPHA picks and updates a random subset of the blocks $\{1,2,\dots,n\}$, using an arbitrary sampling. That is, we allow for the distribution of the random set-valued mapping to be arbitrary (and as explained further below, analyze the iteration complexity of the method).

To the best of our knowledge, there are only two methods in the literature with the “arbitrary sampling” property, the NSync method of Richtárik and Takáč (focusing on the simple problem of minimizing a smooth strongly convex function) and the QUARTZ method of Qu, Richtárik and Zhang (a primal-dual method; considering strongly convex but possibly nonsmooth functions appearing in machine learning). Hence, our work is complementary to this development.

Complexity analysis.

We study the iteration complexity of ALPHA. That is, for an arbitrary (but “proper”) sampling, we provide bounds on the number of iterations needed to approximately solve the problem, in expectation. Our general bounds are formulated in Section 6: Theorem 6.1 covers the non-accelerated variant with $O(1/k)$ rate and Theorem 6.2 covers the accelerated variant with $O(1/k^{2})$ rate, where $k$ is the iteration counter. To the best of our knowledge, these are the first complexity results for a randomized coordinate descent methods utilizing an arbitrary sampling for problem (1).

Expected separable overapproximation.

Besides the dependence of the complexity bound on the iteration counter $k$, it is important to study its dependence on the sampling $\hat{S}$ and the objective function. Our results make this dependence explicit: they hold under the assumption that $f$ admits an expected separable overapproximation (ESO) with respect to the sampling $\hat{S}$ (Assumption 2.1). This is an inequality involving $f$ and $\hat{S}$ which determines certain important parameters $v_{1},\dots,v_{n}$ which are needed to run the method (they determine the stepsizes) and which also appear in the complexity bounds. In some cases it is possible to design a sampling which optimizes the complexity bound.

In the case of a serial sampling, which is by far the most common type of sampling studied in conjunction with randomized coordinate descent methods ($\hat{S}$ is serial if $|\hat{S}|=1$ with probability 1), the parameter $v_{i}$ can simply be set to the Lipschitz constant of the block-derivative of $f$ corresponding to block $i$. In particular, if $n=1$, then $v_{1}$ is the Lipschitz constant of the gradient of $f$ . The situation is more complicated in the case of a parallel sampling ($\hat{S}$ is parallel if it is not serial; that is, if we allow for multiple blocks to be updated at every iteration) – and this why there is a need for the ESO inequality. Intuitively speaking, the parameters $v_{1},\dots,v_{n}$ capture certain smoothness properties of the gradient of $f$ in a random subspace spanned by the blocks selected by the sampling $\hat{S}$.

The ESO concept is of key importance in the design and analysis of randomized coordinate descent methods . We provide a systematic study of ESO inequalities in a companion paper .

Simple complexity analysis in the smooth case.

In order to make the exposition more accessible, we first focus on ALPHA applied to problem (1) with $\psi\equiv 0$ (we call this the “smooth case”). In this simpler setting, it is possible to provide a simplified complexity analysis – we do this in Section 3; see Theorem 3.1 (non-accelerated variant with $O(1/k)$ rate) and Theorem 3.2 (accelerated variant with $O(1/k^{2})$ rate). For convenience, ALPHA specialized to the smooth case is formulated as Algorithm 2. Our analysis in this case is different from the one we give in Section 6, where we analyze the method in the general proximal setup.

Flexibility.

ALPHA is a remarkably flexibleWe have named the method ALPHA because of this flexibility: “ALPHA” as a single source from which one obtains diversity. algorithm, encoding a number of classical, recent and new algorithms in special cases, depending on the choice of the parameters of the method: sampling $\hat{S}$ and “stepsize sequence” $\{\theta_{k}\}$. We devote Section 4 to highlighting several of the many algorithms ALPHA reduces to in special cases, focusing on the smooth case for simplicity (special cases in the general proximal setting are discussed in the appendix). In particular, if $\hat{S}=\{1,2,\dots,n\}$ with probability 1, ALPHA reduces to a deterministic method: gradient descent (GD; Algorithm 3) or accelerated gradient descent (AGD; Algorithm 4), depending on the choice of the sequence $\{\theta_{k}\}$. For a non-deterministic sampling, we obtain parallel coordinate descent (PCD; Algorithm 5) and accelerated parallel coordinate descent (APCD; Algorithm 6) with arbitrary sampling – which is new. If a uniform sampling is used, PCD reduces to the PCDM algorithm . If a distributed Distributed sampling is a structured uniform sampling first introduced in and further studied in and in the companion paper . In a distributed sampling, the blocks $\{1,2,\dots,n\}$ are first partitioned into $c$ sets of equal cardinality (it is useful to think of $c$ to be equal to the number of compute nodes in a distributed computing environment). The sampling is constructed by letting each node choose a subset of a fixed size (say $\tau$) of the blocks it owns, uniformly at random and independently from others, and then taking the union of these random sets. This union is a random subset of the set of blocks; and is called the $(c,\tau)$-distributed sampling. sampling is used instead, PCD reduces to Hydra . Similarly, if a uniform sampling is used, APCD reduces to APPROX (in fact, our version of APPROX is a bit more flexible with respect to choice of $\theta_{0}$, which leads to a better complexity result). APCD specialized to a distributed sampling reduces to Hydra2 .

Robustness.

Since we establish a complexity result for an arbitrary sampling, one of the key contributions of this work is to show that coordinate descent methods are robust to the choice of the sampling $\hat{S}$. In many applications one is forced to sample the coordinates/blocks in a non-traditional way and up to this point the issue of whether the resulting algorithm would converge (let alone the issue of estimating its complexity) was open. For instance, in many metric learning / matrix problems one wishes to find a positive semidefinite matrix satisfying certain properties. It is often efficient to work with an algorithm which would in each iteration update all the elements in a certain row and the corresponding column of the matrix. If we think of the elements of this matrix as coordinates, then any sampling induced in this way puts more probability on the diagonal elements of than on the off-diagonal elements. An algorithm of this type was not analyzed before. The complexity of such a method would follow as a special of our general results specialized to the corresponding sampling.

Improved complexity results.

In all cases where ALPHA reduces to an existing method, our complexity bound either matches the best known bound for that method or improves upon the best known bound. For instance, while the complexity of PCDM (which coincides with PCD specialized to a uniform sampling; Algorithm 5) depends on the size of a certain level-set of $f$ (and in particular, requires the level set to be bounded), our bound does not involve this quantity (see Section A.3). Another example is APPROX (which is closely related to APCD specialized to a uniform sampling): we obtain a more compact and improved result.

Two in one.

We provide a unified complexity analysis covering the nonaccelerated and accelerated variants of ALPHA. This is achieved by establishing a certain key technical recursion (Lemma 6.4) for an arbitrary choice of the parameters $\{\theta_{k}\}$. Since the two variants of ALPHA differ in the choice of this sequence only, the analysis of both is identical up to this point. The recursion is then analyzed in two different ways, depending on the sequence $\{\theta_{k}\}$, which leads to the final complexity result.

Efficient implementation.

4 Outline of the paper

The paper is organized as follows. In Section 2 we establish notation, describe the ALPHA algorithm and comment on the key assumption: Expected Separable Overapproximation (ESO). We defer the in-depth study ESO inequalities to a companion paper . In Section 3 we give a simple complexity proof of ALPHA in the smooth case ($\psi=0$). Subsequently, in Section 4 we present four algorithms that ALPHA reduces to in special cases, and state the corresponding complexity results, which follow from our general result, in a simplified form. We do this for the benefit of the reader. In Section 5 we provide an equivalent form of writing ALPHA–one leading to an efficient implementation avoiding full dimensional operations. In Section 6 we state and prove the convergence result of ALPHA when applied to the general proximal minimization problem (1). Finally, in Section 7 we conclude and in the appendix we comment on several special cases ALPHA reduces to in the general proximal setup.

The Algorithm

That is, $v\cdot x$ is the vector obtained from $x$ by multiplying its block $i$ by $v_{i}$ for each $i\in[n]$. If all blocks are of size one ($N_{i}=1$ for all $i$), then $v\cdot x=\operatorname{Diag}(v)x$ where $\operatorname{Diag}(v)$ is the diagonal matrix with diagonal vector $v$.

2 ALPHA

In this section we describe ALPHA (Algorithm 1) – a randomized block coordinate descent method for solving (1).

We denote by $\operatorname{dom}\psi$ the domain of the proximal term $\psi$.

Let us extract the relations between the three sequences. Define

Note also that from the definition of $y_{k}$ in Algorithm 1, we have:

3 Expected separable overapproximation

Looking behind the compact notation in which the assumption is formulated, observe that the upper bound is a quadratic function of $h=(h^{1},\dots,h^{n})$, separable in the blocks $h^{i}$:

As a tool for the design and analysis of randomized coordinate descent methods, ESO was first formulated in for the complexity study of parallel coordinate descent method (PCDM). It is a powerful technical tool which provides a generic approach to establishing the convergence of randomized coordinate descent methods of many flavours . As shown in the listed papers as well as our results which follow, the convergence of Algorithm 1 can be established for arbitrary sampling $\hat{S}$ as long as the parameter vector $v$ is chosen such that $(f,\hat{S})\sim ESO(v)$. Moreover, the vector $v$ appears in the convergence result and directly influences the complexity of the method.

Since , the problem of computing efficiently a vector $v$ such that the ESO assumption 2.1 holds has been addressed in many papers for special uniform samplings relevant to practical implementation including serial sampling, $\tau$-nice sampling and distributed sampling , and also for a particular example of nonuniform parallel sampling . In this paper we focus on the complexity analysis of Algorithm 1 and refer the reader to the companion paper for a systematic study of the computation of admissible vector parameter $v$ for arbitrary sampling $\hat{S}$.

Simple complexity analysis in the smooth case

In this section we give a brief complexity analysis of ALPHA in the case when $\psi\equiv 0$; that is, when applied to the following unconstrained smooth convex minimization problem:

While our general theory, which we develop in Section 6, covers also this special case, the analysis we present here is different and simpler. When applied to problem 16, Step 8 in ALPHA has an explicit solution and the method reduces to Algorithm 2.

We now state the complexity result for ALPHA (Algorithm 2) in its nonaccelerated variant.

In particular, if we choose $\theta_{0}=\min_{i}p_{i}$, then for all $k\geqslant 1$,

The next result gives a complexity bound for ALPHA (Algorithm 2) in its accelerated variant.

In particular, if we choose $\theta_{0}=1$, then for all $k\geqslant 1$,

In the rest of this section, we provide a short proof of Theorems 3.1 and 3.2. In Section 6 we shall present complexity bounds (Theorem 6.1 and 6.2) for ALPHA (Algorithm 1) as applied to the general regularized problem (1). The proof in the general case is more involved, which is why we prefer to present the smooth case first and also provide a separate briefer proof.

We first establish two lemmas and then proceed directly to the proofs of the theorems.

Based on line 8 of Algorithm 2, we can write

or equivalently, $-\nabla f(y_{k})=\theta_{k}(v\circ p^{-1})\cdot\mathbf{B}a$. Using this notation, the update on line 10 of Algorithm 2 can be written as

Substituting these expressions to (25), we obtain the recursion:

It now only remains to apply Lemma 3.1 to the last two terms in (26), with $x\leftarrow z_{k}-y$, $w\leftarrow v\circ p^{-2}$ and $\hat{S}\leftarrow S_{k}$, and rearrange the resulting inequality. ∎

2 Proof of Theorem 3.1

Using the fact that $\theta_{k}=\theta_{0}$, for all $k$ and taking expectation in both sides of (22), we obtain the recursion

Let $\alpha_{k}=1+(k-1)\theta_{0}$. By convexity,

Finally, subtracting $f(y)$ from both sides and taking expectations, we obtain

3 Proof of Theorem 3.2

If $\theta_{0}\in(0,1]$, then the sequence $\{\theta_{k}\}_{k\geqslant 0}$ has the following properties (see ):

After dividing both sides of (22) by $\theta_{k}^{2}$, using (29) and taking expectations, we obtain:

where $\phi_{k}$ and $r_{k}$ are as in the proof of Theorem 3.1. Finally,

Note that in the last inequality we used the assumption that $C\geqslant 0$.

The many variants of ALPHA (in the smooth case)

The purpose of this section is to demonstrate that ALPHA is a very flexible method, encoding several classical as well as modern optimization methods for special choices of the parameters of the method. In order to achieve this goal, it is enough to focus on the smooth case, i.e., on Algorithm 2. Similar reasoning can be applied to the proximal case.

Note that in ALPHA we have the liberty to choose the sampling $\hat{S}$ and the sequence $\{\theta_{k}\}_{k\geqslant 0}$. As we have already seen, by modifying the sequence we can obtain simple (i.e., nonaccelerated) and accelerated variants of the method. By the choice of the sampling, we can force the method to be deterministic or randomized. In the latter case, there are many ways of choosing the distribution of the sampling. Here we will constrain ourselves to a basic classification between uniform samplings (samplings for which $p_{i}=p_{i^{\prime}}$ for all $i\in[n]$) and non-uniform or importance samplings. This is summarized in Table 1.

The deterministic variants of ALPHA (Algorithm 3 and 4) are obtained by choosing the sampling which always selects all blocks: $\hat{S}=[n]$ with probability 1. The ESO assumption in this special case has the form

In the randomized variants of ALPHA (Algorithm 5 and 6) we allow for the sampling to have an arbitrary distribution.

By specializing Algorithm 2 to the choice $\hat{S}=[n]$ and $\theta_{k}=1$ for all $k$, we obtain classical gradient descent (with fixed stepsize). Indeed, note that in this special case we have

Recall that the ESO assumption reduces to (31) when $\hat{S}=[n]$.

The complexity of the method is a corollary of Theorem 3.1.

For any optimal solution $x_{*}$ of (16), the output of Algorithm 3 for all $k\geqslant 1$ satisfies:

then $f(x_{k})-f(x_{*})\leqslant\epsilon$.

By letting $y=x_{k}$ in (22) we know that:

Note that for this special case $x_{k}=z_{k}$. Therefore,

where the second inequality follows from applying Theorem 3.1 to $\theta_{0}=1$ and $\hat{S}=[n]$. ∎

Corollary 4.1 is a basic result and can be found in many textbooks on convex optimization; see for example .

2 Special case 2: accelerated gradient descent

Let us still keep $\hat{S}=[n]$, but assume now the sequence $\{\theta_{k}\}_{k\geqslant 0}$ is chosen according to (19). In this case, Algorithm 2 reduces to accelerated gradient descent.

Note that only two sequences $\{x_{k},z_{k}\}_{k\geqslant 0}$ are explicitly used in Algorithm 4. This is achieved by replacing $y_{k}$ in Algorithm 2 by $(1-\theta_{k})x_{k}+\theta_{k}z_{k}$. The following result follows directly from Theorem 3.2 by letting $\theta_{0}=1$ and $p_{i}=1$ for all $i\in[n]$.

For any optimal solution $x_{*}$ of (16), the output of Algorithm 4 for all $k\geqslant 1$ satisfies:

then $f(x_{k})-f(x_{*})\leqslant\epsilon$.

Algorithm 4 is a special case of Algorithm 1 in . The complexity bound (35) was also proved in [41, Corollary 1]. See also .

3 Special case 3: Parallel coordinate descent

We now allow the method (Algorithm 2) to use an arbitrary sampling $\hat{S}$, but keep $\theta_{k}=\theta_{0}$ for all $k\geqslant 1$. This leads to Algorithm 5.

For any optimal solution $x_{*}$ of (16), the output of Algorithm 5 for all $k\geqslant 1$ satisfies:

In the special case when $\hat{S}$ is the serial uniform sampling, the three sequences $\{x_{k},y_{k},z_{k}\}_{k\geqslant 0}$ coincide, and one can show that the following bound holds:

Randomized coordinate descent with serial and importance sampling (in a form different from Algorithm 5) was considered by Nesterov . In the special case of serial uniform sampling when Algorithm 5 is the same as in , the following convergence rate was proved in :

where $\mathcal{R}(x_{0})$ is a weighted level-set distance to the set of optimal points $X_{*}$:

Our result does not require the level sets of $f$ to be bounded.

4 Special case 4: Accelerated parallel coordinate descent

To obtain the accelerated coordinate descent method, as a special case of Algorithm 2, we only need to let the sequence $\{\theta_{k}\}_{k\geqslant 0}$ satisfy (19).

The convergence result then follows directly as a corollary of Theorem 3.2.

For any optimal solution $x_{*}$ of (16), the output of Algorithm 6 for all $k\geqslant 1$ satisfies:

An accelerated coordinate descent method for unconstrained minimization in the special case of serial uniform sampling was first proposed and analyzed by Nesterov , where the following bound was proved:

Each serial sampling $\hat{S}$ is uniquely characterized by the vector of probabilities $p=(p_{1},\dots,p_{n})$ where $p_{i}$ is defined by (8). Suppose that the function $f$ has block-Lipschitz gradient with constants $L_{1},\dots,L_{n}$:

The optimal serial sampling given by (43) is not very useful without (at least some) knowledge of $x_{*}$, which is not known. However, note that the formula (43) confirms the intuition that blocks with larger $L_{i}$ and larger distance to the optimal block $\|x_{*}^{i}-x_{0}^{i}\|_{i}$ should be picked (and hence updated) more often.

Efficient implementation

Focusing on the iterates $x_{k},y_{k},z_{k}$ in Algorithm 1 only, the general algorithm can schematically be written as follows:

Consider the change of variables from $\{x_{k},y_{k},z_{k}\}$ to $\{z_{k},g_{k}\}$ where

and $\{\alpha^{k}\}_{k\geqslant 0}$ is a sequence defined by:

Note that, in all the special cases presented in Section 4 and 6, either $\theta_{k}<1$ for all $k\geqslant 1$ or $\theta_{k}=1$ for all $k\geqslant 1$. The latter case does not require the full-dimensional operation $y_{k}=(1-\theta_{k})x_{k}+\theta_{k}z_{k}$ because the three sequences equal to each other. We thus only address the case when $\theta_{k}<1$ for all $k\geqslant 1$ which implies that $\alpha_{k}\neq 0$ for all $k\geqslant 1$.

Then from $\{z_{k},g_{k}\}$ and $\{\alpha_{k}\}$ we can recover $\{x_{k},y_{k}\}$ as follows:

Moreover, $g_{k+1}$ can be computed recursively as follows:

Hence Algorithm 1 can be written in the following equivalent form.

2 Cost of a single iteration

In order to perform Step 7, it is important that we have access to $\nabla_{i}f(y_{k})=\nabla_{i}f(\alpha_{k}g_{k}+z_{k})$ without actually computing $y_{k}$. In , the authors show that this is possible for problems (1) where $f$ can be written as:

For $i\in[n]$ and $j\in[m]$ denote by $\mathbf{A}_{ji}$ the $i$th block of the $j$th row vector of the matrix $\mathbf{A}$, i.e., $\mathbf{A}_{ji}=U_{i}^{\top}\mathbf{A}^{\top}e_{j}$. For each $i\in[n]$, denote by $I_{i}$ the number of rows containing a non-zero $i$th block, i.e.,

Then for $f$ taking the form of (57) we have

where by abuse of notation $u_{k}^{j}$ and $w_{k}^{j}$ denote respectively the $j$th element of the vectors $u_{k}$ and $w_{k}$. With the knowledge of the vectors $u_{k}$ and $w_{k}$, computing $\nabla_{i}f(y_{k})$ requires $O(|I_{i}N_{i}|)$ operations. Now in order to keep record of the vectors $u_{k}$ and $w_{k}$, we use the following equality:

Since $\mathbf{A}\mathbf{U}_{i}$ is a matrix with $|I_{i}N_{i}|$ nonzero elements, the updating schemes (58) and (59) require then

operations. Note that this is also the total complexity of gradient computation $\nabla_{i}f(y_{k})$ at $k$th iteration. Denote by $\operatorname{nnz}(\mathbf{A})$ the total number of nonzero blocks of the matrix $A$, i.e.,

Let us consider the special case when $|\hat{S}|=\tau$ and $N_{i}=N/n$ for all $i\in[n]$. In this case, the expected one iteration computational complexity is:

To make it more direct to understand, let us consider the case when each block contains only one coordinate, i.e., $N=n$. Then the latter expected one iteration complexity becomes

Hence, in this case the one iteration complexity in expectation of Algorithm 7 is of order $O(\tau\bar{\omega})$ where $\bar{\omega}$ is the average number of nonzero elements of the columns of $\mathbf{A}$. Not considering the time spent on synchronization and handling read/write conflicts, the average processing time would be $O(\bar{\omega})$ if we use a parallel implementation with $\tau$ processors.

Proximal minimization

In this section we present and prove complexity results for ALPHA (Algorithm 1) as applied to the general problem (1) involving the proximal term. We leave the discussion concerning special cases to the appendix.

In the presence of the proximal term $\psi$, the same complexity bounds as those given in Theorems 3.1 and 3.2 hold for the output of Algorithm 1, with the exception that $\theta_{0}$ is only allowed to be chosen between $(0,\min_{i}p_{i}]$. We now state the formal complexity theorems, first in the nonaccelerated and then in the accelerated case.

In the remainder of the section we will provide the complexity analysis.

Our approach is similar to that presented in , but with many modifications required because we allow for an arbitrary sampling. We begin with some technical lemmas.

2 Technical lemmas

Lemma 6.1 shows that each individual block $x^{i}_{k}$ of the variable $x_{k}$ is a convex combination of all the history blocks $z_{0}^{i},\dots,z_{k}^{i}$. Note that due to the importance sampling, the combination coefficients $\gamma_{k,0}^{i},\dots,\gamma_{k,k}^{i}$ is now block-dependent, in contrast with the block-independent coefficients proved in .

where for each $i$, the coefficients $\{\gamma_{k,l}^{i}\}_{l=0,\dots,k}$ are defined recursively by setting $\gamma_{0,0}^{i}=1$, $\gamma_{1,0}^{i}=1-\theta_{0}p_{i}^{-1}$, $\gamma_{1,1}^{i}=\theta_{0}p_{i}^{-1}$ and for $k\geqslant 1$,

Fix any $i\in[n]$. We proceed by induction on $k$. It is clear from $x_{0}=z_{0}$ that $\gamma_{0,0}^{i}=1$. Since $x_{1}=y_{0}+\theta_{0}p^{-1}\cdot(z_{1}-z_{0})$ and $y_{0}=x_{0}$, we get that $x_{1}^{i}=(1-\theta_{0}p_{i}^{-1})z_{0}^{i}+\theta_{0}p_{i}^{-1}z_{1}^{i}$ thus $\gamma_{1,0}^{i}=1-\theta_{0}p_{i}^{-1}$ and $\gamma_{1,1}^{i}=\theta_{0}p_{i}^{-1}$. Assuming that (62) holds for some $k\geqslant 1$, then

Therefore the recursive equation (63) holds. The identity (64) can then be verified by direct substitution. Next we assume that $\theta_{0}\in(0,\min_{i}p_{i}]$ and $\{\theta_{k}\}_{k\geqslant 0}$ is a decreasing positive sequence and show that the linear combination in (62) is a convex combination. Let $k\geqslant 1$. Since $\theta_{k}\leqslant 1$, we deduce from (63) that $\{\gamma_{k+1,l}^{i}\}_{l=0,\dots,k-1}$ are positive if $\{\gamma_{k,l}^{i}\}_{l=0,\dots,k-1}$ are positive. Moreover,

Then using $\theta_{k}\geqslant 0$, we conclude that $\gamma_{k+1,k}^{i}\geqslant 0$ if $\gamma_{k,k}^{i}\leqslant 1$. Besides, we have:

The next result was previously stated and used in .

For a proof, see for example [2, Lemma 3.2].

3 Recursion

We next prove an inequality similar to the one we established in the smooth case (22).

and then bound the expectation of $\hat{F}_{k+1}$ as follows:

4 Proof of Theorems 6.1 and 6.2

Using the same reasoning as that in the proof of Theorems 3.1 and 3.2, we analyze recursion of Lemma 6.4 and obtain:

If $\theta_{k}=\theta_{0}\in(0,1]$ for all $k\geqslant 0$, then

If $\theta_{k+1}=\frac{\sqrt{\theta_{k}^{4}+4\theta_{k}^{2}}-\theta_{k}^{2}}{2}$ for all $k\geqslant 0$, $\theta_{0}\in(0,1]$ and $\hat{F}_{0}\leqslant F(y)$, then

Finally, by Lemma 6.1, for the latter two choices of $\{\theta_{k}\}$, if in addition $\theta_{0}\in(0,\min_{i}p_{i}]$, then for $k\geqslant 1$, each block of the vector $x_{k}$ is a convex combination of the corresponding blocks of the vectors $z_{0},\dots,z_{k}$. By the convexity of each function $\psi^{i}$, we get:

Conclusion

In this paper we propose a general randomized coordinate descent method which can be specialized to serial or parallel and accelerated or non-accelerated variants, with or without importance sampling. Based on the technical assumption which captures in a compact way certain smoothness properties of the function in a random subspace spanned by the sampled coordinates, we provide a unified complexity analysis which allows to derive as direct corollary the convergence results for the multiple variants of the general algorithm. We focused on the minimization of non-strongly convex function. Further study on a unified algorithm and complexity analysis for both strongly and non-strongly convex objective functions can be investigated.

References

Appendix A Special cases of ALPHA in the proximal setup

Extending the discussion presented in Section 4 which focused on the smooth case, we now present four special cases of Algorithm 1 for solving problem (1).

Specializing Algorithm 1 to $\hat{S}=[n]$ and $\theta_{k}=1$ for all $k\geqslant 1$, we obtain the classical proximal gradient descent algorithm. Note that in this case the three sequences $\{x_{k},y_{k},z_{k}\}_{k\geqslant 0}$ reduce to one sequence.

For any optimal solution $x_{*}$ of (1), the output of Algorithm 8 for all $k\geqslant 1$ satisfies:

then $F(x_{k})-F(x_{*})\leqslant\epsilon$.

The proof follows as a corollary of Theorem 6.1, with additional remark (32) which holds in this special case. The reader can refer to the proof of Corollary 4.1.

Corollary A.1 can be found in classical textbooks on convex optimization, see for example .

A.2 Special case 2: accelerated proximal gradient descent

By choosing $\hat{S}=[n]$ (with probability 1), $\theta_{0}=1$ and the sequence $\{\theta_{k}\}_{k\geqslant 0}$ according to (19), ALPHA (Algorithm 1) reduces to accelerated proximal gradient descent.

For any optimal solution $x_{*}$ of (1), the output of Algorithm 9 for all $k\geqslant 1$ satisfies:

In particular, for $0<\epsilon<\frac{1}{2}\|x_{0}-x_{*}\|^{2}_{v}$, if

then $F(x_{k})-F(x_{*})\leqslant\epsilon$.

As discussed for the unconstrained case in Section 4.2, Algorithm 9 and Corollary A.2 can be attributed to Tseng .

A.3 Special case 3: Parallel Coordinate Descent Method (PCDM)

By applying Theorem 6.1 and using the same reasoning as in the proof of Corollary 4.1 we obtain the following new convergence result for PCDM.

For any optimal solution $x_{*}$ of (1), the output of Algorithm 10 for all $k\geqslant 1$ satisfies:

In particular, for $0<\epsilon<\left(1-\tau/n\right)\left(F(x_{0})-F(x_{*})\right)+\frac{1}{2}\|x_{0}-x_{*}\|^{2}_{v}$, if

A high-probability result involving the level-set distance

was provided in for PCDM (Algorithm 10). Although not explicitly stated in the paper, it is apparent from the proof that their approach yields the following rate:

Since $\mathcal{R}_{v}(x_{0},x_{*})\geqslant\|x_{0}-x_{*}\|_{v}^{2}$, it is clear that for sufficiently large $k$, our rate (77) is better than (79) .

A.4 Special case 4: APPROX with importance sampling

Finally, let $\{\theta_{k}\}_{k\geqslant 0}$ be chosen in accordance with (19). In this case, ALPHA (Algorithm 1) reduces to an accelerated coordinate descent method with arbitrary sampling $\hat{S}$ for solving the proximal minimization problem 1.

For any optimal solution $x_{*}$ of (1), the output of Algorithm 11 for all $k\geqslant 1$ satisfies:

This is a direct corollary of Theorem 6.2 by taking $\theta_{0}=\min_{i}p_{i}$. ∎