Distributed Coordinate Descent Method for Learning with Big Data

Introduction

Randomized coordinate descent methods (CDMs) are increasingly popular in many learning tasks, including boosting, large scale regression and training linear support vector machines. CDMs update a single randomly chosen coordinate at a time by moving in the direction of the negative partial derivative (for smooth losses). Methods of this type, in various settings, were studied by several authors, including Hsieh et al. (2008); Shalev-Shwartz & Tewari (2009); Nesterov (2012); Richtárik & Takáč (2012c); Necoara et al. (2012); Tappenden et al. (2013b); Shalev-Shwartz & Zhang (2013b); Lu & Xiao (2013).

It is clear that in order to utilize modern shared-memory parallel computers, more coordinates should be updated at each iteration. One way to approach this is via partitioning the coordinates into blocks, and operating on a single randomly chosen block at a time, utilizing parallel linear algebra libraries. This approach was pioneered by Nesterov (2012) for smooth losses, and was extended to regularized problems in (Richtárik & Takáč, 2012c). Another popular approach involves working with a random subset of coordinates (Bradley et al., 2011). These approaches can be combined, and theory was developed for methods that update a random subset of blocks of coordinates at a time Richtárik & Takáč (2012a); Fercoq & Richtárik (2013). Further recent works on parallel coordinate descent include (Richtárik & Takáč, 2012b; Mukherjee et al., 2013; Fercoq, 2013; Tappenden et al., 2013a; Shalev-Shwartz & Zhang, 2013a).

However, none of these methods are directly scalable to problems of sizes so large that a single computer is unable to store the data describing the instance, or is unable to do so efficiently (e.g., in memory). In a big data scenario of this type, it is imperative to split the data across several nodes (computers) of a cluster, and design efficient methods for this memory-distributed setting.

Hydra. In this work we design and analyze the first distributed coordinate descent method: Hydra: HYbriD cooRdinAte descent. The method is “hybrid” in the sense that it uses parallelism at two levels: i) across a number of nodes in a cluster and ii) utilizing the parallel processing power of individual nodesWe like to think of each node of the cluster as one of the many heads of the mythological Hydra..

Assume we have $c$ nodes (computers) available, each with parallel processing power. In Hydra, we initially partition the coordinates $\{1,2,\dots,d\}$ into $c$ sets, ${\cal P}_{1},\dots,{\cal P}_{c}$, and assign each set to a single computer. For simplicity, we assume that the partition is balanced: $|{\cal P}_{k}|=|{\cal P}_{l}|$ for all $k,l$. Each computer owns the coordinates belonging to its partition for the duration of the iterative process. Also, these coordinates are stored locally. The data matrix describing the problem is partitioned in such a way that all data describing features belonging to ${\cal P}_{l}$ is stored at computer $l$. Now, at each iteration, each computer, independently from the others, chooses a random subset of $\tau$ coordinates from those they own, and computes and applies updates to these coordinates. Hence, once all computers are done, $c\tau$ coordinates will have been updated. The resulting vector, stored as $c$ vectors of size $s=d/c$ each, in a distributed way, is the new iterate. This process is repeated until convergence. It is important that the computations are done locally on each node, with minimum communication overhead. We comment on this and further details in the text.

The main insight. We show that the parallelization potential of Hydra, that is, its ability to accelerate as $\tau$ is increased, depends on two data-dependent quantities: i) the spectral norm of the data ($\sigma$) and ii) a partition-induced norm of the data ($\sigma^{\prime}$). The first quantity completely describes the behavior of the method in the $c=1$ case. If $\sigma$ is small, then utilization of more processors (i.e., increasing $\tau$) leads to nearly linear speedup. If $\sigma$ is large, speedup may be negligible, or there may be no speedup whatsoever. Hence, the size of $\sigma$ suggests whether it is worth to use more processors or not. The second quantity, $\sigma^{\prime}$, characterizes the effect of the initial partition on the algorithm, and as such is relevant in the $c>1$ case. Partitions with small $\sigma^{\prime}$ are preferable. For both of these quantities we derive easily computable and interpretable estimates ($\omega$ for $\sigma$ and $\omega^{\prime}$ for $\sigma^{\prime}$), which may be used by practitioners to gauge, a-priori, whether their problem of interest is likely to be a good fit for Hydra or not. We show that for strongly convex losses, Hydra outputs an $\epsilon$-accurate solution with probability at least $1-\rho$ after $\tfrac{d\beta}{c\tau\mu}\log(\tfrac{1}{\epsilon\rho})$ iterations (we ignore some small details here), where a single iteration corresponds to changing of $\tau$ coordinates by each of the $c$ nodes; $\beta$ is a stepsize parameter and $\mu$ is a strong convexity constant.

Outline. In Section 2 we describe the structure of the optimization problem we consider in this paper and state assumptions. We then proceed to Section 3, in which we describe the method. In Section 4 we prove bounds on the number of iterations sufficient for Hydra to find an approximate solution with arbitrarily high probability. A discussion of various aspects of our results, as well as a comparison with existing work, can be found in Section 5. Implementation details of our distributed communication protocol are laid out in Section 6. Finally, we comment on our computational experiments with a big data (3TB matrix) L1 regularized least-squares instance in Section 7.

The problem

We study the problem of minimizing regularized loss,

where $f$ is a smooth convex loss, and $R$ is a convex (and possibly nonsmooth) regularizer.

and write $\mathbf{M}=\mathbf{A}^{T}\mathbf{A}$, where $\mathbf{A}$ is some $n$-by-$d$ matrix.

Example. These assumptions are natural satisfied in many popular problems. A typical loss function has the form

Regularizer R𝑅R.

Example. The choice $R_{i}(t)=0$ for $t\in$ and $R_{i}(t)=+\infty$, otherwise, effectively models bound constraints, which are relevant for SVM dual. Other popular choices are $R(x)=\lambda\|x\|_{1}$ (L1-regularizer) and $R(x)=\tfrac{\lambda}{2}\|x\|_{2}^{2}$ (L2-regularizer).

Distributed coordinate descent

We consider a setup with $c$ computers (nods) and first partition the $d$ coordinates (features) into $c$ sets ${\cal P}_{1},\dots,{\cal P}_{c}$ of equal cardinality, $s:=d/c$, and assign set ${\cal P}_{l}$ to node $l$. Hydra is described in Algorithm 1. Hydra’s convergence rate depends on the partition; we comment on this later in Sections 4 and 5. Here we simply assume that we work with a fixed partition. We now comment on the steps.

Step 3. At every iteration, each of the $c$ computers picks a random subset of $\tau$ features from those that it owns, uniformly at random, independently of the choice of the other computers. Let $\hat{S}_{l}$ denote the set picked by node $l$ . More formally, we require that i) $\hat{S}_{l}\subseteq{\cal P}_{l}$, ii) $\mathbf{Prob}(|\hat{S}_{l}|=\tau)=1$, where $1\leq\tau\leq s$, and that iii) all subsets of ${\cal P}_{l}$ of cardinality $\tau$ are chosen equally likely. In summary, at every iteration of the method, features belonging to the random set $\hat{S}:=\cup_{l=1}^{c}\hat{S}_{l}$ are updated. Note that $\hat{S}$ has size $c\tau$, but that, as a sampling from the set $\{1,2,\dots,d\}$, it does not choose all cardinality $c\tau$ subsets of $\{1,2,\dots,d\}$ with equal probability. Hence, the analysis of parallel coordinate descent methods of Richtárik & Takáč (2012a) does not apply. We will say that $\hat{S}$ is a $\tau$-distributed sampling with respect to the partition $\{{\cal P}_{1},\dots,{\cal P}_{c}\}$.

Step 4. Once computer $l$ has chosen its set of $\tau$ coordinates to work on in Step 3, it will in parallel compute (Step 5) and apply (Step 6) updates to them.

Step 5. This is a critical step where updates to coordinates $i\in\hat{S}_{l}$ are computed. By $f^{\prime}_{i}(x)$ we denote the $i$-th partial derivative of $f$ at $x$. Notice that the formula is very simple as it involves one dimensional optimization.

Closed-form formulas. Often, $h_{k}^{i}$ can be computed in closed form. For $R_{i}(t)=\lambda_{i}|t|$ (weighted L1 regularizer), $h_{k}^{i}$ is the point in the interval $[\tfrac{-\lambda_{i}-f^{\prime}_{i}(x_{k})}{\mathbf{M}_{ii}\beta},\tfrac{\lambda_{i}-f^{\prime}_{i}(x_{k})}{\mathbf{M}_{ii}\beta}]$ which is closest to $-x_{k}^{i}$. If $R_{i}(t)=\tfrac{\lambda_{i}}{2}t^{2}$ (weighted L2 regularizer), then $h_{k}^{i}=-\tfrac{f^{\prime}_{i}(x_{k})}{\lambda_{i}\mathbf{M}_{ii}\beta}$.

Choice of $\beta$. The choice of the step-size parameter $\beta$ is of paramount significance for the performance of the algorithm, as argued for different but related algorithms by Richtárik & Takáč (2012a); Takáč et al. (2013); Fercoq & Richtárik (2013). We will discuss this issue at length in Sections 4 and 5.

Implementation issues: Note that computer $l$ needs to know the partial derivatives of $f$ at $x_{k}$ for coordinates $i\in\hat{S}_{l}\subseteq{\cal P}_{l}$. However, $x_{k}$, as well as the data describing $f$, is distributed among the $c$ computers. One thus needs to devise a fast and communication efficient way of computing these derivatives. This issue will be dealt with in Section 6.

Step 6. Here all the $\tau$ updates computed in Step 5 are applied to the iterate. Note that the updates are local: computer $l$ only updates coordinates it owns, which are stored locally. Hence, this step is communication-free.

Step 7. Here we are just establishing a way of labeling iterates. That is, starting with $x_{k}$, all $c$ computers modify $c\tau$ entries of $x_{k}$ in total, in a distributed way, and the result is called $x_{k+1}$. Our method is therefore inherently synchronous. We do not allow, in our analysis, for the various computers to proceed until all computers have updated all coordinates. In practice, a carefully designed asynchronous implementation will be faster, and our experiments in Section 7 are done with such an implementation.

Convergence rate analysis

Here we define two quantities, $\sigma^{\prime}$ and $\sigma$, which, as we shall see, play an important role in the computation of the stepsize parameter $\beta$ of Algorithm 1, and through it, in understanding its rate of convergence and potential for speedup by parallelization and distribution. As we shall see, these quantities might not be easily computable. We therefore also provide each with an easily computable and interpretable upper bound, $\omega^{\prime}$ for $\sigma^{\prime}$ and $\omega$ for $\sigma$.

for each $k,l\in\{1,2,\dots,c\}$. We now define

A useful consequence of (6) is the inequality

Let $a_{rl}$ be the $r$-th row of $\mathbf{A}_{l}$, and define

where $\omega^{\prime}(r)$ is the number of matrices $\mathbf{A}_{l}$ with a nonzero in row $r$. Likewise, define

where $\omega(r)$ is the number of nonzeros in the $r$-th row of $\mathbf{A}$.

2 Choice of the stepsize parameter β𝛽\beta

We analyze Hydra with stepsize parameter $\beta\geq\beta^{*}$, where

and $s_{1}=\max(1,s-1)$. As we shall see in Theorem 5, fixing $c$ and $\tau$, the number of iterations needed by Hydra find a solution is proportional to $\beta$. Hence, we would wish to use $\beta$ which is as small as possible, but not smaller than the safe choice $\beta=\beta^{*}$, for which convergence is proved. In practice, $\beta$ can often be chosen smaller than $\beta^{*}$, leading to larger steps and faster convergence. If the quantities $\sigma$ and $\sigma^{\prime}$ are hard to compute, then one can replace them by the easily computable upper bounds $\omega$ and $\omega^{\prime}$, respectively. However, there are cases when $\sigma$ can be efficiently approximated and is much smaller than $\omega$. In some ML datasets with $\mathbf{A}\in\{0,1\}^{n\times d}$, $\sigma$ is close to the average number of nonzeros in a row of $\mathbf{A}$, which can be significantly smaller than the maximum, $\omega$. On the other hand, if $\sigma$ is difficult to compute, $\omega$ may provide a good proxy. Similar remarks apply to $\sigma^{\prime}$. In the $\tau\geq 2$ case (which covers all interesting uses of Hydra), we may ignore $\beta_{2}^{*}$ altogether, as implied by the following result.

If $\tau\geq 2$, then $\beta^{*}\leq 2\beta_{1}^{*}$.

This eliminates the need to compute $\sigma^{\prime}$, at the expense of at most doubling $\beta$, which translates into doubling the number of iterations.

3 Separable approximation

We first establish a useful identity for the expected value of a random quadratic form obtained by sampling the rows and columns of the underlying matrix via the distributed sampling $\hat{S}$. Note that the result is a direct generalization of Lemma 1 in (Takáč et al., 2013) to the $c>1$ case.

where $\alpha_{1}=1-\tfrac{\tau-1}{s_{1}}$, $\alpha_{2}=\tfrac{\tau-1}{s_{1}}$, $\alpha_{3}=\tfrac{\tau}{s}-\tfrac{\tau-1}{s_{1}}$.

We now use the above lemma to compute a separable quadratic upper bound on $\mathbf{E}[(h^{\hat{S}})^{T}\mathbf{M}h^{\hat{S}}]$.

For $x:=(D^{\mathbf{M}})^{1/2}h$, we have $(h^{\hat{S}})^{T}\mathbf{M}h^{\hat{S}}=(x^{\hat{S}})^{T}\mathbf{Q}x^{\hat{S}}$. Taking expectations on both sides, and applying Lemma 3, we see that $\mathbf{E}[(h^{\hat{S}})^{T}\mathbf{M}h^{\hat{S}}]$ is equal to (11) for $\mathbf{G}=\mathbf{Q}$. It remains to bound the three quadratics in (11). Since $D^{\mathbf{Q}}$ is the identity matrix, $x^{T}D^{\mathbf{Q}}x=h^{T}D^{\mathbf{M}}h$. In view of (7), the 2nd term is bounded as $x^{T}\mathbf{Q}x\leq\sigma x^{T}x=\sigma h^{T}D^{\mathbf{M}}h$. The last term, $x^{T}(\mathbf{Q}-B^{\mathbf{Q}})$, is equal to

It only remains to plug in these three bounds into (11). ∎

Inequalities of type (12) were first proposed and studied by Richtárik & Takáč (2012a)—therein called Expected Separable Overapproximation (ESO)—and were shown to be important for the convergence of parallel coordinate descent methods. However, they studied a different class of loss functions $f$ (convex smooth and partially separable) and different types of random samplings $\hat{S}$, which did not allow them to propose an efficient distributed sampling protocol leading to a distributed algorithm. An ESO inequality was recently used by Takáč et al. (2013) to design a mini-batch stochastic dual coordinate ascent method (parallelizing the original SDCA methods of Hsieh et al. (2008)) and mini-batch stochastic subgradient descent method (Pegasos of Shalev-Shwartz et al. (2011)), and give bounds on how mini-batching leads to acceleration. While it was long observed that mini-batching often accelerates Pegasos in practice, it was only shown with the help of an ESO inequality that this is so also in theory. Recently, Fercoq & Richtárik (2013) have derived ESO inequalities for smooth approximations of nonsmooth loss functions and hence showed that parallel coordinate descent methods can accelerate on their serial counterparts on a class of structured nonsmooth convex losses. As a special case, they obtain a parallel randomized coordinate descent method for minimizing the logarithm of the exponential loss. Again, the class of losses considered in that paper, and the samplings $\hat{S}$, are different from ours. None of the above methods are distributed.

4 Fast rates for distributed learning with Hydra

where $\phi^{\prime}(x)$ is a subgradient (or gradient) for $\phi$ at $x$.

We now show that Hydra decreases strongly convex $L$ with an exponential rate in $\epsilon$.

where $\beta\geq\beta^{*}$ and $\beta^{*}$ is given by (10). If $\{x_{k}\}$ are the random points generated by Hydra (Algorithm 1), then

To see this, substitute $h\leftarrow h^{\hat{S}}$ into (2), take expectations on both sides and then use Lemma 4 together with the fact that for any vector $a$, $\mathbf{E}[a^{T}h^{\hat{S}}]=\tfrac{\mathbf{E}[|\hat{S}|]}{d}=\tfrac{\tau c}{sc}=\tfrac{\tau}{s}$. The rest follows by following the steps in the proof in (Richtárik & Takáč, 2012a, Theorem 20). ∎

A similar result, albeit with the weaker rate $O(\tfrac{s\beta}{\tau\epsilon})$, can be established in the case when neither $f$ nor $R$ are strongly convex. In big data setting, where parallelism and distribution is unavoidable, it is much more relevant to study the dependence of the rate on parameters such as $\tau$ and $c$. We shall do so in the next section.

Discussion

In this section we comment on several aspects of the rate captured in (13) and compare Hydra to selected methods.

Here we comment in detail on the influence of the various design parameters ($c$ = # computers, $s$ = # coordinates owned by each computer, and $\tau$ = # coordinates updated by each computer in each iteration), instance-dependent parameters ($\sigma,\omega,\mu_{R},\mu_{f}$), and parameters depending both on the instance and design ($\sigma^{\prime},\omega^{\prime}$), on the stepsize parameter $\beta$, and through it, on the convergence rate described in Theorem 5.

Notice that the size of $\mu_{R}>0$ mitigates the effect of a possibly large $\beta$ on the bound (13). Indeed, for large $\mu_{R}$, the factor $(\beta+\mu_{R})/(\mu_{f}+\mu_{R})$ approaches 1, and the bound (13) is dominated by the term $\tfrac{d}{c\tau}$, which means that Hydra enjoys linear speedup in $c$ and $\tau$. In the following comments we will assume that $\mu_{R}=0$, and focus on studying the dependence of the leading term $d\tfrac{\beta}{c\tau}$ on various quantities, including $\tau,c,\sigma$ and $\sigma^{\prime}$.

Search for small but safe β𝛽\beta.

As shown by Takáč et al. (2013, Section 4.1), mini-batch SDCA might diverge in the setting with $\mu_{f}=0$ and $R(x)\equiv 0$, even for a simple quadratic function with $d=2$, provided that $\beta=1$. Hence, small values of $\beta$ need to be avoided. However, in view of Theorem 5, it is good if $\beta$ is as small as possible. So, there is a need for a “safe” formula for a small $\beta$. Our formula (10), $\beta=\beta^{*}$, is serving that purpose. For a detailed introduction into the issues related to selecting a good $\beta$ for parallel coordinate descent methods, we refer the reader to the first 5 pages of (Fercoq & Richtárik, 2013).

If $c=1$, then by Lemma 9, $\sigma^{\prime}=c=1$, and hence $\beta^{*}_{2}=0$. However, for $c>1$ we may have $\beta^{*}_{2}>0$, which can hence be seen as a price we need to pay for using more nodes. The price depends on the way the data is partitioned to the nodes, as captured by $\sigma^{\prime}$. In favorable circumstances, $\sigma^{\prime}\approx 1$ even if $c>1$, leading to $\beta_{2}^{*}\approx 0$. However, in general we have the bound $\sigma^{\prime}\geq\tfrac{c\sigma}{d}$, which gets worse as $c$ increases and, in fact, $\sigma^{\prime}$ can be as large as $c$. Note also that $\xi$ is decreasing in $\tau$, and that $\xi(s,s)=0$. This means that by choosing $\tau=s$ (which effectively removes randomization from Hydra), the effect of $\beta^{*}_{2}$ is eliminated. This may not be always possible as often one needs to solve problems with $s$ vastly larger than the number of updates that can be performed on any given node in parallel. If $\tau\ll s$, the effect of $\beta^{*}_{2}$ can be controlled, to a certain extent, by choosing a partition with small $\sigma^{\prime}$. Due to the way $\sigma^{\prime}$ is defined, this may not be an easy task. However, it may be easier to find partitions that minimize $\omega^{\prime}$, which is often a good proxy for $\sigma^{\prime}$. Alternatively, we may ignore estimating $\sigma^{\prime}$ altogether by setting $\beta=2\beta_{1}^{*}$, as mentioned before, at the price of at most doubling the number of iterations.

Speedup by increasing τ𝜏\tau.

Let us fix $c$ and compare the quantities $\gamma_{\tau}:=\tfrac{\beta^{*}}{c\tau}$ for $\tau=1$ and $\tau=s$. We now show that $\gamma_{1}\geq\gamma_{s}$, which means that if all coordinates are updated at every node, as opposed to one only, then Hydra run with $\beta=\beta^{*}$ will take fewer iterations. Comparing the 1st and 3rd row of Table 2, we see that $\gamma_{1}=s+\sigma\tfrac{\sigma^{\prime}-1}{\sigma^{\prime}}$ and $\gamma_{s}=\sigma$. By Lemma 1, $\gamma_{1}-\gamma_{s}=s-\tfrac{\sigma}{\sigma^{\prime}}\geq 0$.

Price of distribution.

For illustration purposes, consider a problem with $d=10^{5}$ coordinates. In Figure 1(a) we depict the size of $\tfrac{d\beta_{1}^{*}}{c\tau}$ for $c=1$ and several choices of $\tau$, as a function of $\sigma$. We see that Hydra works better for small values of $\sigma$ and that with increasing $\sigma$, the benefit of using updating more coordinates diminishes. In Figure 1(a) we consider the same scenario, but with $c=100$ and $s=1000$, and we plot $\tfrac{d2\beta_{1}^{*}}{c\tau}$ on the $y$ axis. Note that the red dotted line in both plots corresponds to a parallel update of 1600 coordinates. In (a) all are updated on a single node, whereas in (b) we have 100 nodes, each updating 16 coordinates at a time. Likewise, the dashed blue dashed and solid black lines are also comparable in both plots. Note that the setup with $c=10$ has a slightly weaker performance, the lines are a bit lower. This is the price we pay for using $c$ nodes as opposed to a single node (obviously, we are ignoring communication cost here). However, in big data situations one simply has no other choice but to utilize more nodes.

2 Comparison with other methods

While we are not aware of any other distributed coordinate descent method, Hydra in the $c=1$ case is closely related to several existing parallel coordinate descent methods.

The Shotgun algorithm (parallel coordinate descent) of Bradley et al. (2011) is similar to Hydra for $c=1$. Some of the differences: Bradley et al. (2011) only consider $R$ equal to the $L1$ norm and their method works in dimension $2d$ instead of the native dimension $d$. Shotgun was not analyzed for strongly convex $f$, and convergence in expectation was established. Moreover, Bradley et al. (2011) analyze the step-size choice $\beta=1$, fixed independently of the number of parallel updates $\tau$, and give results that hold only in a “small $\tau$” regime. In contrast, our analysis works for any choice of $\tau$.

Hydra vs PCDM.

For $c=1$, Hydra reduces to the parallel coordinate descent method (PCDM) of Richtárik & Takáč (2012a), but with a better stepsize parameter $\beta$. We were able to achieve smaller $\beta$ (and hence better rates) because we analyze a different and more specialized class of loss functions (those satisfying (2)). In comparison, Richtárik & Takáč (2012a) look at a general class of partially separable losses. Indeed, in the $c=1$ case, our distributed sampling $\hat{S}$ reduces to the sampling considered in (Richtárik & Takáč, 2012a) ($\tau$-nice sampling). Moreover, our formula for $\beta$ (see Table 2) is essentially identical to the formula for $\beta$ provided in (Richtárik & Takáč, 2012a, Theorem 14), with the exception that we have $\sigma$ where they have $\omega$. By 9, we have $\sigma\leq\omega$, and hence our $\beta$ is smaller.

Hydra vs SPCDM.

SPCDM of Fercoq & Richtárik (2013) is PCDM applied to a smooth approximation of a nonsmooth convex loss; with a special choice of $\beta$, similar to $\beta_{1}$. As such, it extends the reach of PCDM to a large class of nonsmooth losses, obtaining $O(\tfrac{1}{\epsilon^{2}})$ rates.

Hydra vs mini-batch SDCA.

Takáč et al. (2013) studied the performance of a mini-batch stochastic dual coordinate ascent for SVM dual (“mini-batch SDCA”). This is a special case of our setup with $c=1$, convex quadratic $f$ and $R_{i}(t)=0$ for $t\in$ and $R_{i}(t)=+\infty$ otherwise. Our results can thus be seen as a generalization of the results in that paper to a larger class of loss functions $f$, more general regularizers $R$, and most importantly, to the distributed setting ($c>1$). Also, we give $O(\log\tfrac{1}{\epsilon})$ bounds under strong convexity, whereas (Takáč et al., 2013) give $O(\tfrac{1}{\epsilon})$ results without assuming strong convexity. However, Takáč et al. (2013) perform a primal-dual analysis, whereas we do not.

Distributed computation of the gradient

In this section we described some important elements of our distributed implementation.

Note that in Hydra, $x_{k}$ is stored in a distributed way. That is, the values $x_{k}^{i}$ for $i\in{\cal P}_{l}$ are stored on computer $l$. Moreover, Hydra partitions $\mathbf{A}$ columnwise as $\mathbf{A}=[\mathbf{A}_{1},\dots,\mathbf{A}_{c}]$, where $\mathbf{A}_{l}$ consists of columns $i\in{\cal P}_{l}$ of $\mathbf{A}$, and stores $\mathbf{A}_{l}$ on computer $l$. So, $\mathbf{A}$ is chopped into smaller pieces with stored in a distributed way in fast memory (if possible) across the $c$ nodes. Note that this allows the method to work with large matrices.

If we write $h_{k}^{i}=0$ if $i$ is not updated in iteration $k$, then

Note that the value $\delta g_{k,l}$ can be computed on node $l$ as all the required data is stored locally. Hence, we let each node compute $\delta g_{k,l}$, and then use a reduce all operation to add up the updates to obtain $g_{k+1}$, and pass the sum to all nodes. Knowing $g_{k+1}$, node $l$ is then able to compute $f^{\prime}_{i}(x_{k+1})$ for any $i\in{\cal P}_{l}$ as follows:

2 Advanced protocols

The basic protocol discussed above has obvious drawbacks. Here we identify them and propose modifications leading to better performance.

alternating Parallel and Serial regions (PS): The basic protocol alternates between two procedures: i) a computationally heavy one (done in parallel) with no MPI communication, and ii) MPI communication (serial). An easy fix would be to dedicate 1 thread to deal with communication and the remaining threads within the same computer for computation. We call this protocol Fully Parallel (FP). Figure 2 compares the basic (left) and FP (right) approaches.

Reduce All (RA): In general, reduce all operations may significantly degrade the performance of distributed algorithms. Communication taking place only between nodes close to each other in the network, e.g., nodes directly connected by a cable, is more efficient. Here we propose the Asynchronous StreamLined (ASL) communication protocol in which each node, in a given iteration, sends only 1 message (asynchronously) to a nearby computer, and also receives only one message (asynchronously) from another nearby computer. Communication hence takes place in an Asynchronous Ring. This communication protocol requires significant changes in the algorithm. Figure 3 illustrates the flow of messages at the end of the $k$-th iteration for $c=4$.

We order the nodes into a ring, denoting $l_{-}$ and $l_{+}$ the two nodes neighboring node $l$. Node $l$ only receives data from $l_{-}$, and sends data to $l_{+}$. Let us denote by $\delta G_{k,l}$ the data sent by node $l$ to $l_{+}$ at the end of iteration $k$. When $l$ starts iteration $k$, it already knows $\delta G_{k-1,l_{-}}$.Initially, we let $\delta g_{k,l}=\delta G_{k,l}=0$ for all $k\leq 0$. Hence, data which will be sent at the end of the $k$-th iteration by node $l$ is given by

ASL needs less communication per iteration. On the other hand, information is propagated more slowly to the nodes through the ring, which may adversely affect the number of iterations till convergence (note that we do not analyze Hydra with this communication protocol). Indeed, it takes $c-1$ iterations to propagate information to all nodes. Also, storage requirements have increased: at iteration $k$ we need to store the vectors $\delta g_{t,l}$ for $k-c\leq t\leq k$ on computer $l$.

Experiments

As discussed in Section 6, the advantage of the RA protocol is the fact that Theorem 5 was proved in this setting, and hence can be used as a safe benchmark for comparison with the advanced protocols.

Table 4 compares the average time per iteration for the 3 approaches and 3 choices of $\tau$. We used $128$ nodes, each running 4 MPI processes (hence $c=512$). Each MPI process runs 8 OpenMP threads, giving 4,096 cores in total. The data matrix $\mathbf{A}$ has $n=10^{9}$ rows and $d=5\times 10^{8}$ columns, and has 3 TB, double precision. One can observe that in all cases, ASL-FP yields largest gains compared to the benchmark RA-PS protocol. Note that ASL has some overhead in each iteration, and hence in cases when computation per node is small ($\tau=10$), the speedup is only 1.62. When $\tau=10^{2}$ (in this case the durations of computation and communication were comparable), ASL-FP is 3.11 times faster than RA-PS. But the gain becomes again only moderate for $\tau=10^{3}$; this is because computation now takes much longer than communication, and hence the choice of strategy for updating the auxiliary vector $g_{k}$ is less significant. Let us remark that the use of larger $\tau$ requires larger $\beta$, and hence possibly more iterations (in the worst case).

Huge LASSO problem.

We generated a sparse matrix $\mathbf{A}$ with block angular structure, depicted in (17).

Extensions

Our results can be extended to the setting where coordinates are replaced by blocks of coordinates, as in (Nesterov, 2012), and to partially separable losses, as in (Richtárik & Takáč, 2012a).

References

Appendix A Proof Lemma 1

The inequality $\omega^{\prime}\leq c$ is obviously true. By considering $x$ with zeroes in all coordinates except those that belong to ${\cal P}_{k}$ (where $k$ is an arbitrary but fixed index), we see that $x^{T}\mathbf{Q}x=x^{T}B^{\mathbf{Q}}x$, and hence $\sigma^{\prime}\geq 1$.

Let $\mathbf{U}_{k}$ be a column submatrix of the $d$-by-$d$ identity matrix corresponding to columns $i\in{\cal P}_{k}$. Clearly, $\mathbf{A}_{k}=\mathbf{A}\mathbf{U}_{k}$ and $\mathbf{U}_{k}^{T}\mathbf{Q}e_{k}$ is the $k$-th diagonal block of $\mathbf{Q}$, i.e.,

This means that $\phi^{\prime}$ is block Lipschitz (with blocks corresponding to variables in ${\cal P}_{k}$), with respect to the norm $\|\cdot\|_{(k)}$, with Lipschitz constant $1$. Indeed, this is, in fact, satisfied with equality:

This is relevant because then, by Richtárik & Takáč (2012a, Theorem 7; see comment 2 following the theorem), it follows that $\phi^{\prime}$ is Lipschitz with respect to $\|\cdot\|_{w}$, where $w_{k}=1$ for all $k=1,\dots,c$, with Lipschitz constant $\omega^{\prime}$ ($\omega^{\prime}$ is the degree of partial block separability of $\phi$ with respect to the blocks ${\cal P}_{k}$). Hence,

which establishes the inequality $\sigma^{\prime}\leq\omega^{\prime}$.

We now show that $\tfrac{\sigma}{s}\leq\sigma^{\prime}$. If we let $\theta:=\max\{x^{T}B^{\mathbf{Q}}x:x^{T}x\leq 1\}$, then $x^{T}B^{\mathbf{Q}}x\leq\theta x^{T}x$ and hence $\{x\;:\;x^{T}x\leq 1\}\subseteq\{x\;:\;x^{T}B^{\mathbf{Q}}x\leq\theta\}$. This implies that

In the last step we have used the fact that $\sigma(\mathbf{Q})=\sigma\leq c=\dim(\mathbf{Q})$, proved in steps 1 and 2, applied to the setting $\mathbf{Q}\leftarrow\mathbf{Q}^{kk}$.

The chain of inequalities $1\leq\sigma\leq\omega\leq c$ is obtained as a special case of the chain $1\leq\sigma^{\prime}\leq\omega^{\prime}\leq d$ (proved above) when $c=d$ (and hence ${\cal P}_{l}=\{l\}$ for $l=1,\dots,d$). Indeed, in this case $B^{\mathbf{Q}}=D^{\mathbf{Q}}$, and so $x^{T}B^{\mathbf{Q}}x=x^{T}D^{\mathbf{Q}}x=x^{T}x$, which means that $\sigma^{\prime}=\sigma$ and $\omega^{\prime}=\omega$.

Appendix B Proof of Lemma 2

It is enough to argue that $\beta_{2}^{*}\leq\beta_{1}^{*}$. Notice that $\beta_{2}^{*}$ is increasing in $\sigma^{\prime}$. On the other hand, from Lemma 1 we know that $\sigma^{\prime}\leq c=\tfrac{d}{s}$. So, it suffices to show that

After straightforward simplification we observe that this inequality is equivalent to $(s-\tau)+(\tau-2)\sigma+\tfrac{\sigma}{d}(s+\tau)\geq 0$, which clearly holds.

Appendix C Proof of Lemma 3

In the $s=1$ case the statement is trivially true. Indeed, we must have $\tau=1$ and thus $\mathbf{Prob}(\hat{S}=\{1,2,\dots,d\})=1$, $h^{\hat{S}}=h$, and hence

This finishes the proof since $\tfrac{\tau-1}{s_{1}}=0$.

Consider now the $s>1$ case. From Lemma 3 in Richtárik & Takáč (2012a) we get

where $p_{ij}=\mathbf{Prob}(i\in\hat{S}\;\&\;j\in\hat{S})$. One can easily verify that

In particular, the first case follows from Eq (32) and the second from Eq (37) in Richtárik & Takáč (2012a). It only remains to substitute $p_{ij}$ into (22) and transform the result into the desired form.