Fast Distributed Coordinate Descent for Non-Strongly Convex Losses

Introduction

In this paper we are concerned with solving regularized convex optimization problems in huge dimensions in cases when the loss being minimized is not strongly convex. That is, we consider the problem

For examples of regularizers $R$ and loss functions $f$ satisfying the above assumptions, relevant to machine learning applications, we refer the reader to .

It is increasingly the case in modern applications that the data describing the problem (encoded in $\mathbf{A}$ and $R$ in the above model) is so large that it does not fit into the RAM of a single computer. In such a case, unless the application at hand can tolerate slow performance due to frequent HDD reads/writes, it is often necessary to distribute the data among the nodes of a cluster and solve the problem in a distributed manner.

While efficient distributed methods exist in cases when the regularized loss $L$ is strongly convex (e.g., Hydra ), here we assume that $L$ is not strongly convex. Problems of this type arise frequently: for instance, in the SVM dual, $f$ is typically a non-strongly convex quadratic, $d$ is the number of examples, $n$ is the number of features, $\mathbf{A}$ encodes the data, and $R$ is a -$\infty$ indicator function encoding box constraints (e.g., see ). If $d>n$, then $L$ will not be strongly convex.

In this paper we propose “Hydra2” (Hydra “squared”; Algorithm 1) – the first distributed stochastic coordinate descent (CD) method with fast $O(1/k^{2})$ convergence rate for our problem. The method can be seen both as a specialization of the APPROX algorithm to the distributed setting, or as an accelerated version of the Hydra algorithm (Hydra converges as $O(1/k)$ for our problem). The core of the paper forms the development of new stepsizes, and new efficiently computable bounds on the stepsizes proposed for distributed CD methods in . We show that Hydra2 is able to solve a big data problem with $d$ equal to 50 billion.

The Algorithm

Assume we have $c$ nodes/computers available. In Hydra2 (Algorithm 1), the coordinates $i\in[d]:=\{1,2,\dots,d\}$ are first partitioned into $c$ sets $\{{\cal P}_{l},\;l=1,\dots,c\}$, each of size $|P_{l}|=s:=d/c$. The columns of $\mathbf{A}$ are partitioned accordingly, with those belonging to partition ${\cal P}_{l}$ stored on node $l$. During one iteration, all computers $l\in\{1,\dots,c\}$, in parallel, pick a subset $\hat{S}_{l}$ of $\tau$ coordinates from those they own, i.e., from ${\cal P}_{l}$, uniformly at random, where $1\leqslant\tau\leqslant s$ is a parameter of the method (Step 6​). From now on we denote by $\hat{S}$ the union of all these sets, $\hat{S}:=\cup_{l}\hat{S}_{l}$, and refer to it by the name distributed sampling.

The main work is done in Step 8​ where the update scalars $t_{k}^{i}$ are computed. This step is a proximal forward-backward operation similar to what is to be done in FISTA but we only compute the proximal operator for the coordinates in $\hat{S}$.

The partial derivatives are computed at $(\theta_{k}^{2}u_{k}+z_{k})$. For the algorithm to be efficient, the computation should be performed without computing the sum $\theta_{k}^{2}u_{k}+z_{k}$. As shown in , this can be done for functions $f$ of the form

where for all $j$, $\phi_{j}$ is a one-dimensional differentiable (convex) function. Indeed, let $D_{i}$ be the set of such $j$ for which $\mathbf{A}_{ji}\neq 0$. Assume we store and update $r_{u_{k}}=\mathbf{A}u_{k}$ and $r_{z_{k}}=\mathbf{A}z_{k}$, then

The average cost for this computation is $\frac{\tau}{s}\sum_{i=1}^{d}{\cal O}(|D_{i}|)$, i.e., each processor needs to access the elements of only one column of the matrix $\mathbf{A}$.

Steps 8​ and 9​ depend on a deterministic scalar sequence $\theta_{k}$, which is being updated in Step 12​ as in . Note that by taking $\theta_{k}=\theta_{0}$ for all $k$, $u_{k}$ remains equal to 0, and Hydra2 reduces to Hydra .

The output of the algorithm is $x_{k+1}=(\theta_{k}^{2}u_{k+1}+z_{k+1})$. We only need to compute this vector sum at the end of the execution and when we want to track $L(x_{k})$. Note that one should not evaluate $x_{k}$ and $L(x_{k})$ at each iteration since these computations have a non-negligible cost.

Convergence rate

where $\mathbf{D}$ is a diagonal matrix with diag. elements $\mathbf{D}_{ii}>0$ and $\hat{S}$ is the distributed sampling described above.

The above ESO assumption involves the smooth function $f$, the sampling $\hat{S}$ and the parameters $\{\mathbf{D}_{ii}\}_{i=1}^{d}$. It has been first introduced by Richtárik and Takáč for proposing a generic approach in the convergence analysis of the Parallel Coordinate Descent Methods (PCDM). Their generic approach boils down the convergence analysis of the whole class of PCDMs to the problem of finding proper parameters which make the ESO assumption hold. The same idea has been extended to the analysis of many variants of PCDM, including the Accelerated Coordinate Descent algorithm (APPROX) and the Distributed Coordinate Descent method (Hydra).

In particular, the following complexity result, under the ESO assumption 3.1, can be deduced from [4, Theorem 3] using Markov inequality:

Choose $0<\rho<1$ and $0<\epsilon<L_{0}$ and

Then under Assumption 3.1, the iterates $\{x_{k}\}_{k\geqslant 1}$ of Algorithm 1 satisfy $\mathbf{Prob}(L(x_{k})-L(x_{*})\leqslant\epsilon)\geqslant 1-\rho\enspace.$

The bound on the number of iterations $k$ in (4) suggests to choose the smallest possible $\{\mathbf{D}_{ii}\}_{i=1}^{d}$ satisfying the ESO assumption 3.1.

New Stepsizes

In this section we propose new stepsize parameters $\mathbf{D}_{ii}$ for which the ESO assumption is satisfied.

Let $\omega_{j}$ be the number of nonzeros in the $j$th row of $\mathbf{A}$ and $\omega^{\prime}_{j}$ be the number of “partitions active at row $j$”, i.e., the number of indexes $l\in\{1,\dots,c\}$ for which the set $\{i\in{\cal P}_{l}:\mathbf{A}_{ji}\neq 0\}$ is nonempty. Note that since $\mathbf{A}$ does not have an empty row or column, we know that $1\leqslant\omega_{j}\leqslant n$ and $1\leqslant\omega_{j}^{\prime}\leqslant c$. Moreover, we have the following characterization of these two quantities. For notational convenience we denote $\mathbf{M}_{j}=\mathbf{A}_{j:}^{\top}\mathbf{A}_{j:}$ for all $j\in\{1,\dots,n\}$.

Let $S^{\prime}=\{l:\mathbf{A}_{j:}^{(l)}\neq 0\}$, then $\omega_{j}^{\prime}=|S^{\prime}|$ and by Cauchy-Schwarz we have

Equality is reached when $\mathbf{A}_{j:}^{(l)}x^{(l)}=\alpha$ for some constant $\alpha$ for all $l\in S^{\prime}$ (this is feasible since the subsets $\{\mathcal{P}_{1},\dots,\mathcal{P}_{c}\}$ are disjoint). Hence, we proved (6). The characterization (5) follows from (6) by setting $c=d$. ∎

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}}$.

Below is the main result of this section.

For a convex differentiable function $f$ satisfying (2) and a distributed sampling $\hat{S}$ described in Section 2, the ESO assumption 3.1 is satisfied for

Existing Stepsizes and New Bounds

We now desribe stepsizes previously suggested in the literature. For simplicity, let $\mathbf{M}:=\mathbf{A}^{\top}\mathbf{A}$. Define:

The quantities $\sigma$ and $\sigma^{\prime}$ are identical to those defined in (although the definitions are slightly different). The following stepsize parameters have been introduced in :

where $\beta^{*}:=\beta^{*}_{1}+\beta^{*}_{2}$ and

In general, the computation of $\sigma$ and $\sigma^{\prime}$ can be done using the Power iteration method . However, the number of operations required in each iteration of the Power method is at least twice as the number of nonzero elements in $\mathbf{A}$. Hence the only computation of the parameters $\sigma$ and $\sigma^{\prime}$ would already require quite a few number of passes through the data. Instead, if one could provide some easily computable upper bound for $\sigma$ and $\sigma^{\prime}$, where by ’easily computable’ we mean computable by only one pass through the data, then we can run Algorithm 1 immediately without spending too much time on the computation of $\sigma$ and $\sigma^{\prime}$. Note that the ESO assumption will still hold when $\sigma$ and $\sigma^{\prime}$ in (13) are replaced by some upper bounds.

In , the authors established the following bound

which is independent of the partition $\{\mathcal{P}_{l}\}_{l=1,\dots,c}$. This bound holds for $\tau\geqslant 2$. Further, they showed that

Then in view of (15), (16) and Lemma 5.1, the ESO assumption 3.1 is also satisfied for the following easily computable parameters:

2 Improved bounds on existing stepsizes

In what follows, we show that both (15) and (16) can be improved so that smaller parameters $\{\mathbf{D}_{ii}\}_{i=1}^{d}$ are allowed in the algorithm.

Suppose $\tau\geqslant 2$ (note that then $s\geqslant\tau\geqslant 2$). For all $1\leqslant\eta\leqslant s$ the following holds

Both sides of the inequality are linear functions of $\eta$. It therefore suffices to verify that the inequality holds for $\eta=1$ and $\eta=s$, which can be done. ∎

The following result is an improvement on (15), which was shown in [3, Lemma 2].

If $\tau\geqslant 2$, then $\beta^{*}\leqslant(1+\frac{1}{\tau-1})\beta_{1}^{*}$.

We only need to apply Lemma 5.2 to (14) with $\eta=\sigma$, and additionally use the bound $\tfrac{\sigma^{\prime}-1}{\sigma^{\prime}}\leqslant 1$. This gives $\beta_{2}^{*}\leqslant\beta_{1}^{*}/(\tau-1)$, from which the result follows. ∎

The above lemma leads to an important insight: as long as $\tau\geqslant 2$, the effect of partitioning the data (across the nodes) on the iteration complexity of Hydra2 is negligible, and vanishes as $\tau$ increases.

Indeed, consider the complexity bound provided by stepsizes $\mathbf{D}^{2}$ given in (13), used within Theorem 3.2 and note that the effect of the partitioning on the complexity is through $\sigma^{\prime}$ only (which is not easily computable!), which appears in $\beta_{2}^{*}$ only. However, by the above lemma, $\beta_{1}^{*}\leqslant\beta^{*}\leqslant(1+\tfrac{1}{\tau-1})\beta_{1}^{*}$ for any partitioning (to $c$ sets each containing $s$ coordinates). Hence, i) each partitioning is near-optimal and ii) one does not even need to know $\sigma^{\prime}$ to run the method as it is enough to use stepsizes $\mathbf{D}^{2}$ with $\beta^{*}$ replaced by the upper bound provided by Lemma 5.3.

Note that the same reasoning can be applied to the Hydra method since the same stepsizes can be used there. Also note that for $\tau=1$ the effect of partitioning can have a dramatic effect on the complexity.

Applying the same reasoning to the new stepsizes $\alpha_{j}^{*}$ we can show that $\alpha_{j}^{*}\leqslant(1+1/(\tau-1))\alpha_{1,j}^{*}$ for all $j\in\{1,\dots,n\}$. This is not as useful as Lemma 5.3 since $\alpha_{2,j}^{*}$ does not need to be approximated ($\omega_{j}^{\prime}$ is easily computable).

We next give a tighter upper bound on $\sigma$ than (16).

The rest follows from the definition of $\sigma$. ∎

By combining Lemma 5.1, Lemma 5.3 and Lemma 5.4, we obtain the following smaller (compared to (17)) admissible and easily computable stepsize parameters:

Comparison of old and new stepsizes

So far we have seen four different admissible stepsize parameters, some old and some new. For ease of reference, let us call $\{\mathbf{D}^{1}_{ii}\}_{i=1}^{d}$ the new parameters defined in (8), $\{\mathbf{D}^{2}_{ii}\}_{i=1}^{d}$ the existing one given in (see (13)), $\{\mathbf{D}^{3}_{ii}\}_{i=1}^{d}$ the upper bound of $\{\mathbf{D}^{2}_{ii}\}_{i=1}^{d}$ used in (see (17)) and $\{\mathbf{D}^{4}_{ii}\}_{i=1}^{d}$ the new upper bound of $\{\mathbf{D}^{2}_{ii}\}_{i=1}^{d}$ defined in (20).

In the next result we compare the four stepsizes. The lemma implies that $\mathbf{D}^{1}$ and $\mathbf{D}^{2}$ are uniformly smaller (i.e., better – see Theorem 3.2) than $\mathbf{D}^{4}$ and $\mathbf{D}^{3}$. However, $\mathbf{D}^{2}$ involves quantities which are hard to compute (e.g., $\sigma$). Moreover, $\mathbf{D}^{4}$ is always smaller than $\mathbf{D}^{3}$.

Let $\tau\geqslant 2$. The following holds for all $i$:

Proof. We only need to show that $\mathbf{D}_{ii}^{1}\leqslant\mathbf{D}_{ii}^{4}$, the other relations are already proved. It follows from (19) that

Now, for all $i$, we can write $\mathbf{D}^{1}_{ii}-\mathbf{D}^{4}_{ii}=$

Note that we do not have a simple relation between $\{\mathbf{D}_{ii}^{1}\}_{i=1}^{d}$ and $\{\mathbf{D}_{ii}^{2}\}_{i=1}^{d}$. Indeed, for some coordinates $i$ the new stepsize parameters $\{\mathbf{D}_{ii}^{1}\}_{i=1}^{d}$ could be smaller than $\{\mathbf{D}_{ii}^{2}\}_{i=1}^{d}$, see Figure 1. in the next section for an illustration.

Numerical Experiments

In this section we present preliminary computational results involving real and synthetic data sets. In particular, we first compare the relative benefits of the 4 stepsize rules, and then demonstrate that Hydra2 can significantly outperform Hydra on a big data non-strongly convex problem involving 50 billion variables.

The experiments were performed on two problems:

Dual of SVM: This problem can be formulated as finding $x\in^{d}$ that minimizes

The computation were performed on Archer (http://archer.ac.uk/), which is UK’s leading supercomputer.

In this section we compare the four different stepsize parameters and show how they influence the convergence of Hydra2 . In the experiments of this section we used the astro-ph dataset.Astro-ph is a binary classification problem which consists of abstracts of papers from physics. The dataset consists of $d=29,882$ samples and the feature space has dimension $n=99,757$. Note that for the SVM dual problem, coordinates correspond to samples whereas for the LASSO problem, coordinates correspond to features. The problem that we solve is the SVM dual problem (for more details, see ).

Figure 1 plots the values of $\mathbf{D}^{1}_{ii}$, $\mathbf{D}^{2}_{ii}$, $\mathbf{D}^{3}_{ii}$ and $\mathbf{D}^{4}_{ii}$ for $i$ in $\{1,\dots,d\}$ and $(c,\tau)=(32,10)$. The dataset we consider is normalized (the diagonal elements of $\mathbf{A}^{\top}\mathbf{A}$ are all equal to 1), and the parameters $\{\mathbf{D}^{2}_{ii}\}_{i}$ are equal to a constant for all $i$ (the same holds for $\{\mathbf{D}^{3}_{ii}\}_{i}$ and $\{\mathbf{D}^{4}_{ii}\}_{i}$). Recall that according to Theorem 3.2, a faster convergence rate is guaranteed if smaller (but theoretically admissible) stepsize parameters $\{\mathbf{D}_{ii}\}_{i}$ are used. Hence it is clear that $\{\mathbf{D}^{1}_{ii}\}_{i}$ and $\{\mathbf{D}^{2}_{ii}\}_{i}$ are better choices than $\{\mathbf{D}^{3}_{ii}\}_{i}$ and $\{\mathbf{D}^{4}_{ii}\}_{i}$. However, as we mentioned, computing the parameters $\{\mathbf{D}^{1}_{ii}\}_{i}$ would require a considerable computational effort, while computing $\{\mathbf{D}^{2}_{ii}\}_{i}$ is as easy as passing once through the dataset. For this relatively small dataset, the time is about 1 minute for computing $\{\mathbf{D}^{2}_{ii}\}_{i}$ and less than 1 second for $\{\mathbf{D}^{1}_{ii}\}_{i}$.

In order to investigate the benefit of the new stepsize parameters, we solved the SVM dual problem on the astro-ph dataset for $(c,\tau)=(32,10)$. Figure 2 shows evolution of the duality gap, obtained by using the four stepsize parameters mentioned previously. We see clearly that smaller stepsize parameters lead to faster convergence, as predicted by Theorem 3.2. Moreover, using our easily computable new stepsize parameters $\{\mathbf{D}_{ii}^{1}\}_{i}$, we achieve comparable convergence speed with respect to the existing parameters $\{\mathbf{D}_{ii}^{2}\}_{i}$.

2 Hydra2 vs Hydra

In this section, we report experimental results comparing Hydra with Hydra2 on a synthetic big data LASSO problem. We generated a sparse matrix $\mathbf{A}$ having the same block angular structure as the one used in [3, Sec 7.], with $d=50$ billion, $c=256$, $s=195,312,500$ (note that $d=cs$) and $n=5,000,000$. The average number of nonzero elements per row of $\mathbf{A}$ (i.e., $\sum_{j}\omega_{j}/n$) is $60,000$, and the maximal number of nonzero elements in a row of $\mathbf{A}$ (i.e., $\max_{j}\omega_{j}$) is $1,993,419$. The dataset size is $5$TB.

We have used 128 physical nodes of the ARCHER www.archer.ac.uk supercomputing facility, which is based on Cray XC30 compute nodes connected via Aries interconnect. On each physical node we have run two MPI processes – one process per NUMA (NonUniform Memory Access) region. Each NUMA region has 2.7 GHz, 12-core E5-2697 v2 (Ivy Bridge) series processor and each core supports 2 hardware threads (Hyperthreads). In order to minimize communication we have chosen $\tau=s/1000$ (hence each thread computed an update for 8,138 coordinates during one iteration, on average).

We show in Figure 4 the decrease of the error with respect to the number of iterations, plotted in log scale. It is clear that Hydra2 provides faster iteration convergence rate than Hydra. Moreover, we see from Figure 4 that the speedup in terms of the number of iterations is significantly large so that Hydra2 converges faster in time than Hydra even though the run time of Hydra2 per iteration is on average two times more expensive.

References