Distributed Optimization with Arbitrary Local Solvers

Motivation

Regression and classification techniques, represented in the general class of regularized loss minimization problems , are among the most central tools in modern big data analysis, machine learning, and signal processing. For these tasks, much effort from both industry and academia has gone into the development of highly tuned and customized solvers. However, with the massive growth of available datasets, major roadblocks still persist in the distributed setting, where data no longer fits in the memory of a single computer, and computation must be split across multiple machines in a network .

On typical real-world systems, communicating data between machines is several orders of magnitude slower than reading data from main memory, e.g., when leveraging commodity hardware. Therefore when trying to translate the highly tuned existing single machine solvers to the distributed setting, great care must be taken to avoid this significant communication bottleneck .

While several distributed solvers for the problems of interest have been recently developed, they are often unable to fully leverage the competitive performance of their highly tuned and customized single machine counterparts, which have already received much more research attention. More importantly, it is unfortunate that distributed solvers cannot automatically benefit from improvements made to the single machine solvers, and therefore are forced to lag behind the most recent developments.

In this paper we make a step towards resolving these issues by proposing a general communication-efficient distributed framework that can employ arbitrary single machine local solvers and thus directly leverage their benefits and problem-specific improvements. Our framework works in rounds, where in each round the local solvers on each machine find a (possibly weak) solution to a specified subproblem of the same structure as the original master problem. On completion of each round, the partial updates between the machines are efficiently combined by leveraging the primal-dual structure of the global problem . The framework therefore completely decouples the local solvers from the distributed communication. Through this decoupling, it is possible to balance communication and computation in the distributed setting, by controlling the desired accuracy and thus computational effort spent to determine each subproblem solution. Our framework holds with this abstraction even if the user wishes to use a different local solver on each machine.

The proposed framework allows for distributed optimization with the use of arbitrary local solvers on each machine. This abstraction makes the resulting framework highly flexible, and means that it can easily leverage the benefits of well-studied, problem-specific single machine solvers. In addition to increased flexibility and ease-of-use, this can result in large performance gains, as single machine solvers for the problems of interest have typically been thoroughly tuned for optimal performance. Moreover, any performance improvements that continue to be made to these local solvers will be automatically translated by the framework into the distributed setting.

On real-world compute systems, the cost of communication versus computation typical varies by many orders of magnitude, from high-performance computing environments to very slow disk-based distributed workflow systems such as MapReduce/Hadoop. For optimization algorithms, it is thus essential to accommodate varying amounts of work performed locally per round, while still providing convergence guarantees. Our framework provides exactly such control.

In this paper we extend and improve upon the CoCoA method. Our theoretical convergence rates apply to both smooth and non-smooth losses, and for both CoCoA as well as the more general CoCoA$\!{}^{\bf\textbf{\footnotesize+}}$ framework here. The framework also exhibits favorable strong scaling properties for the class of problems considered, as the number of machines $K$ increases and the data size is kept fixed. More precisely, while the convergence rate of CoCoA degrades as $K$ is increased, the stronger theoretical convergence rate here is—in the worst case—independent of $K$. Since the number of communicated vectors is only one per round and worker, this favorable scaling might be surprising. Indeed, for existing methods, splitting data among more machines generally increases communication requirements , which can severely affect overall runtime.

We additionally strengthen the rates by showing stronger primal-dual convergence for both algorithmic frameworks, which are almost tight to their dual-only (or primal-only) counterparts. Primal-dual rates for CoCoA had not previously been analyzed in the general convex case. Our primal-dual rates allow efficient and practical certificates for the optimization quality, e.g., for stopping criteria.

Finally, we give an extensive experimental comparison highlighting the impact of using various arbitrary solvers locally on each machine on several real-world, distributed datasets. We compare the performance of CoCoA and CoCoA$\!{}^{\bf\textbf{\footnotesize+}}$ across these datasets and choices of solvers, in particular illustrating the performance on a massive 280 GB dataset. Our code is available in an open source C++ library, at: https://github.com/optml/CoCoA.

2 Outline

The rest of the paper is organized as follows. Section 2 provides context and states the problem of interest, including necessary assumptions and their consequences. In Section 3 we formulate the algorithm in detail and explain how to implement it efficiently in practice. The main theoretical results are presented in Section 4, followed by discussion of relevant related work in Section 5. Practical experiments demonstrating the strength of the proposed framework are given in Section 6. Finally, we prove the main results in the appendix, in Section A.4.

Background and Problem Formulation

To set our framework in context, we first state traditional complexity measures and convergence rates for single machine algorithms, and then demonstrate that these must be adapted to more accurately represent the performance of an algorithm in the distributed setting.

When running an iterative optimization algorithm $\mathcal{A}$ on a single machine, its performance is typically measured by the total runtime:

Here, $\mathcal{T}_{\mathcal{A}}$ stands for the time it takes to perform a single iteration of algorithm $\mathcal{A}$, and $\mathcal{I}_{\mathcal{A}}(\epsilon)$ is the number of iterations $\mathcal{A}$ needs to attain an $\epsilon$-accurate objective.While for many algorithms the cost of a single iteration will vary throughout the iterative process, this simple model will suffice for our purpose to highlight the key issues associated with extending algorithms to a distributed framework.

On a single machine, most of the state-of-the-art first-order optimization methods can achieve quick convergence in practice in terms of (T-A) by performing a large amount of relatively fast iterations. In the distributed setting, however, time to communicate between two machines can be several orders of magnitude slower than even a single iteration of such an algorithm. As a result, the overall time needed to perform a single iteration can increase significantly.

Distributed timing can therefore be better illustrated using the following practical distributed efficiency model (see also ), where

The extra term $c$ is the time required to perform one round of communication.For simplicity, we assume here a fixed network architecture, and compare only to classes of algorithms that communicate a single vector in each iteration, rendering $c$ to be a constant, assuming we have a fixed number of machines. Most first-order algorithms would fall into this class.As a result, an algorithm that performs well in the setting of (T-A) does not necessarily perform well in the distributed setting (T-B), especially when implemented in a straightforward or naïve way. In particular, if $c\gg\mathcal{T}_{\mathcal{A}}$, we could intuitively expect less potential for improvement, as most of the time in the method will be spent on communication, not on actual computational effort to solve the problem. In this setting, novel optimization procedures are needed that carefully consider the amount of communication and the distribution of data across multiple machines.

One approach to this challenge is to design novel optimization algorithms from scratch, designed to be efficient in the distributed setting. This approach has one obvious practical drawback: There have been numerous highly efficient solvers developed, fine-tuned to particular problems of interest, as long as the problem fits onto a single machine. These solvers are ideal if run on a single machine, but with the growth of data and necessity of data distribution, they must be re-designed to work in modern data regimes.

Recent work has attempted to address this issue by designing algorithms that reduce the communication bottleneck by allowing infrequent communication, while utilizing already existing algorithms as local sub-procedures. The presented work here builds on the promising approach of in this direction. See Section 5 for a detailed discussion of the related literature.

The core idea is that one can formulate a local subproblem for each individual machine, and run an arbitrary local solver dependent only on local data for a number of iterations — obtaining a partial local update. After each worker returns its partial update, a global update is formed by their aggregation.

The big advantage of this is that companies and practitioners do not have to implement new algorithms that would be suitable for the distributed setting. We provide a way for them to utilize their existing algorithms that work on a single machine, and provide a novel communication protocol on top of this.

In the original work on CoCoA , authors provide convergence analysis only for the case when the overall update is formed as an average of the partial updates, and note that in practice it is possible to improve performance by making a longer step in the same direction. The main contribution of this work is a more general convergence analysis of various settings, which enables us to do better than averaging. In one case, we can even sum the partial updates to obtain the overall update, which yields the best result, both in theory and practice. We will see that this can result in significant performance gains, see also .

In the analysis, we will allow local solvers of arbitrarily weak accuracy, each working on its own subproblem which is defined in a completely data-local way for each machine. The relative accuracy obtained by each local solver will be denoted by $\Theta\in$, where $\Theta=0$ describes an exact solution of the subproblem, and $\Theta=1$ means that the local subproblem objective has not improved at all, for this call of the local solver. This paradigm results in a substantial change in how we analyze efficiency in the distributed setting. The formula practitioners are interested in minimizing thus changes to become:

Here, the function $\mathcal{T}_{\mathcal{A}}(\Theta)$ represents the time the local algorithm $\mathcal{A}$ needs to obtain accuracy $\Theta$ on the local subproblem. Note that the number of outer iterations $\mathcal{I}(\epsilon,\Theta)$ is independent of choice of the inner algorithm $\mathcal{A}$, which will also be reflected by our convergence analysis presented in Section 4. Our convergence rates will hold for any local solver $\mathcal{A}$ achieving local quality $\Theta$. For strongly convex problems, the general form will be $\mathcal{I}(\epsilon,\Theta)=\frac{\mathcal{O}(\log(1/\epsilon))}{1-\Theta}\,.$ The inverse dependence on $1-\Theta$ suggests that there is a limit to how much we can gain by solving local subproblems to high accuracy — $\Theta$ close to . There will always be order of $\log(1/\epsilon)$ outer iterations needed. Hence, excessive local accuracy should not be necessary. On the other hand, if $\Theta\rightarrow 1$, meaning that the cost and quality of the local solver diminishes, then the number of rounds $\mathcal{I}(\epsilon,\Theta)$ will explode, which is to be expected.

To illustrate the strength of the paradigm (T-C) compared to (T-B), suppose we run gradient descent as local solver $\mathcal{A}$ for just a single iteration. It turns out that within our framework, choosing this local solver would lead to a method which is equivalent to running naively distributed gradient descentNote that this is not obvious at this point. They are identical, subject to choice of local subproblems as specified in Section 3.1.. Running gradient descent for a single iteration would happen to attain a particular value of $\Theta$. Note that we typically do not set the value $\Theta$ explicitly. It is implicitly chosen by the number of iterations or stopping criterion specified by the user for the local solver. There is no reason to think that the value attained by single iteration of gradient descent would be optimal. For instance, it may be the case that running gradient descent for, say, $200$ iterations, instead of just one, would give substantially better result in practice, due to better communication efficiency. These kinds of considerations are discussed at length in Section 6.

In general, one would intuitively expect that the optimal choice would be to have $\Theta$ such that $\mathcal{T}_{\mathcal{A}}(\Theta)=\mathcal{O}(1)\times c$. In practice, however, the best strategy for any given local solver is to try several values for the number of local iterations to estimate the optimal choice. We will further discuss the importance of $\Theta$, both theoretically and empirically, later in Sections 4 and 6.

Here we list the core conceptual properties of Algorithm 1, which are important qualities that allow it to run efficiently.

The optimization algorithm used to solve the local subproblem $\mathcal{G}_{k}$ outputs a vector $\mathbf{h}_{[k]}^{t}$ with nonzero elements only in coordinates corresponding to variables ${\boldsymbol{\alpha}}_{[k]}$ stored locally (i.e., $i\in\mathcal{P}_{k}$).

Let us now comment on these properties in more detail. Locality is important for making the method versatile, and is the way we escape the restricted setting described by (T-B) that allows us much greater flexibility in designing the overall optimization scheme. Local changes result from the fact that along with data, we distribute also coordinates of the dual variable ${\boldsymbol{\alpha}}$ in the same way, and thus only make updates to the coordinates stored locally. As we will see, efficient maintenance of the subproblems can be obtained. For this, a communication efficient encoding of the current shared state ${\boldsymbol{\alpha}}$ is necessary. To this goal, we will in Section 3.2 show that communication of a single $d$-dimensional vector is enough to formulate the subproblems (LO) in each round, by carefully exploiting their partly separable structure.

Note that Algorithm 1 is the “analysis friendly” formulation of our algorithm framework, and it is not yet fully illustrative for implementation purposes. In Section 3.2 we will precisely formulate the actual communication scheme, and illustrate how the above properties can be achieved.

Before that, we will formulate the precise subproblem $\mathcal{G}_{k}$ in the following section.

We are now ready to define the local objective $\mathcal{G}^{\sigma^{\prime}}_{k}(\,\cdot\,;{\boldsymbol{\alpha}})$ as follows:

The interpretation of the above defined subproblems is that they will form a quadratic approximation of the smooth part of the true objective $D$, which becomes separable over the machines. The approximation keeps the non-smooth $R$ part intact. The variable $\mathbf{h}_{[k]}$ expresses the update proposed by machine $k$. In this spirit, note also that the approximation coincides with $D$ at the reference point ${\boldsymbol{\alpha}}$, i.e. $\sum_{k=1}^{K}\mathcal{G}^{\sigma^{\prime}}_{k}({\bf 0};{\boldsymbol{\alpha}})=D({\boldsymbol{\alpha}})$. We will discuss the interpretation and properties of these subproblems in more detail below in Section 3.3.

2 Practical Communication Efficient Implementation

We will now discuss how Algorithm 1 can efficiently be implemented in a distributed environment. Most importantly, it remains to clarify how the “local” subproblems can actually be formulated and solved by using only local information from the corresponding machine, and to make precise what information needs to be communicated in each round.

Here for the reformulation of the gradient term, we have simply used the chain rule on the objective $f$ (recall the definition $f({\boldsymbol{\alpha}}):=\lambda g^{*}(\mathbf{v})$), giving ∇f(α)_[k] = 1n X_[k]^T ∇g^* ( v) .

In summary, we have seen that each machine can formulate the local subproblem given purely local information (the local data $X_{[k]}$ as well as the local dual variables ${\boldsymbol{\alpha}}_{[k]}$). No information about the other machines variables ${\boldsymbol{\alpha}}$ or their data is necessary.

The only requirement for the method to work is that between the rounds, the changes in the ${\boldsymbol{\alpha}}_{[k]}$ variables on each machine and the resulting global change in $\mathbf{v}$ are kept consistent, in the sense that $\mathbf{v}^{t}=\mathbf{v}({\boldsymbol{\alpha}}^{t}):=\tfrac{1}{\lambda n}X{\boldsymbol{\alpha}}^{t}$ must always hold. Note that for the evaluation of $\nabla g^{*}(\mathbf{v})$, the vector $\mathbf{v}$ is all that is needed. In practice, $g$ as well as its conjugate $g^{*}$ are simple vector valued regularization functions, the most prominent example being $g(\mathbf{v})=g^{*}(\mathbf{v})=\frac{1}{2}\|\mathbf{v}\|^{2}$.

The framework as shown below in Algorithm 2 clearly maintains the consistency of ${\boldsymbol{\alpha}}^{t}$ and $\mathbf{v}^{t}=\mathbf{v}^{t}({\boldsymbol{\alpha}}^{t})$ after each round, no matter which local solver is used to approximately solve (LO’). A diagram illustrating the communication and computation involved in the first two full iterations of Algorithm 2 is given in Figure 1.

3 Compatibility of the Subproblems for Aggregating Updates

In this subsection, we shed more light on the local subproblems on each machine, as defined in (LO) above, and their interpretation. More formally, we will show how the aggregation parameters $\nu$ (controlling the level of adding versus averaging the resulting updates from each machine) and $\sigma^{\prime}$ (the subproblem parameter) interplay together, to in each round achieve a valid approximation to the global objective function $D$.

The role of the subproblem parameter $\sigma^{\prime}$ is to measure the difficulty of the given data partition. For the convergence results discussed below to hold, $\sigma^{\prime}$ must be chosen not smaller than

Here, $G$ is the block diagonal submatrix of the data covariance matrix $X^{T}X$, corresponding to the partition $\{\mathcal{P}_{k}\}_{k=1}^{K}$, i.e.,

In this notation, it is easy to see that the crucial quantity defining $\sigma^{\prime}_{min}$ above is written as $\mathbf{h}^{T}G\mathbf{h}=\textstyle\sum_{k=1}^{K}\|X_{[k]}\mathbf{h}_{[k]}\|^{2}$.

The following lemma shows that if the aggregation and subproblem parameters $\nu$ and $\sigma^{\prime}$ satisfy (10), then the sum of the subproblems $\sum_{k}\mathcal{G}^{\sigma^{\prime}}_{k}$ will closely approximate the global objective function $D$. More precisely, this sum is a block-separable lower bound on $D$.

For any aggregation parameter $\nu\in$, the choice of the subproblem parameter $\sigma^{\prime}:=\nu K$ is valid for (10), i.e., $\nu K\geq\sigma^{\prime}_{min}.$

In this section we state the main theoretical results of this paper. Before doing so, we elaborate on one of the most important aspects of the algorithmic framework: the quality of approximate local solutions.

The notion of approximation quality provided by the local solvers is measured according to the following:

The assumption specifies the (relative) accuracy $\Theta$ obtained on solving the local subproblem $\mathcal{G}_{k}$. Considering the two extreme examples, setting $\Theta=0$ would require to find the exact maximum, while $\Theta=1$ states that no improvement was achieved at all by the local solver. Intuitively, we would prefer $\Theta$ to be small, but spending many computational resources to drive $\Theta$ to can be excessive in practice, since $\mathcal{G}_{k}$ is actually not the problem we are interested in solving (2), but is the problem to be solved per communication round. The best choice in practice will therefore be to choose $\Theta$ such that the local solver runs for a time comparable to the time it takes for a single communication round. This freedom of choice of $\Theta\in$ is a crucial property of our proposed framework, allowing it to adapt to the full range of communication speeds on real world systems, ranging from supercomputers on one extreme to very slow communication rounds like MapReduce systems on the other extreme.

In Section 6 we study impact of different values of this parameter to the overall performance on solving (2).

2 Complexity Bounds

we have the expected duality gap E[ P( w(αT)) - D(αT) ]≤ϵGap.

we have that the expected duality gap satisfies

The most important observation regarding the above result is that we do not impose any assumption on the choice of the local solver, apart from sufficient decrease condition on the local objective in Assumption 4.1.

Let us now comment on the leading terms of the complexity results. The inverse dependence on $1-\Theta$ suggests that it is worth pushing the rate of local accuracy $\Theta$ down to zero. However, when thinking about overall complexity, we have to bear in mind that achieving high accuracy on the local subproblems might be too expensive. The optimal choice would depend on the time we estimate a round of communication would take. In general, if communication is slow, it would be worth spending more time on solving local subproblems, but not so much if communication is relatively fast. We discussed this tradeoff in Section 2.

We achieve a significant speedup by replacing the slow averaging aggregation (as in ) by more aggressive adding instead, that is $\nu=1$ instead of $\nu=1/K$. Note that the safe subproblem parameter for the averaging case ($\nu=1/K$) is $\sigma^{\prime}:=1$, while for adding ($\nu=1$) it is given by $\sigma^{\prime}:=K$, both proven in Lemma 3.2. The resulting speedup from more aggressive adding is strongly reflected in the resulting convergence rate as shown above, when plugging in the actual parameter values $\nu$ and $\sigma^{\prime}$ for the two cases, as we will illustrate more clearly in the next subsection.

3 Discussion and Interpretations of Convergence Results

As the above theorems suggest, it is not possible to meaningfully change the aggregation parameter $\nu$ in isolation. It comes naturally coupled with a particular subproblem.

For the case of aggressive adding, the convergence theorems for local objective (LO) results derived in Theorems 4.2 and 4.3 are as follows:

Let the assumptions of Theorem 4.2 be satisfied. If we run Algorithm 1 with $\nu=1,\sigma^{\prime}=K$ for

On the other hand, if we would just average results (as proposed in ), we would obtain following corollary:

Let the assumptions of Theorem 4.2 be satisfied. If we run Algorithm 1 with $\nu=1/K,\sigma^{\prime}=1$ for

Comparing the leading terms in Equations (17) and (18) we see that the leading term for the $\nu=1$ choice is $\mathcal{O}(\lambda\gamma n+\sigma_{\max}K)$, which is always better than for the $\nu=1/K$ case, when the leading term is $\mathcal{O}(\lambda\gamma n+\sigma_{\max}K)$. This strongly suggests that adding in Framework 2 is preferable, especially when $\lambda\gamma n\gg\sigma_{\max}$.

An analogously significant improvement by an order of $K$ factor follows for the case of the sub-linear convergence rate for general Lipschitz loss functions, as shown in Theorem 4.3.

Note that the differences in the convergence rate are bigger for relatively big values of the regularizer $\lambda$. When the regularizer is $\mathcal{O}(1/n)$, the difference is negligible. This behaviour is also present in practice, as we will point out in Section 6.

In this section, we review a number of methods designed to solve optimization problems of the form of interest here, which are typically referred to as regularized empirical risk minimization (ERM) in the machine learning literature. Formally described in Section 2.1, this problem class (1) underlies many prominent methods of supervised machine learning.

Stochastic Gradient Descent (SGD) is the simplest stochastic method one can use to solve the problem of structure (1), and dates back to the work of Robbins and Monro . We refer the reader to for recent theoretical and practical assessment of SGD. Generally speaking, the method is extremely easy to implement, and converges to modest accuracies very quickly, which is often satisfactory in applications in machine learning. On the other hand, difficulty in choosing hyper-parameters make the method sometimes rather cumbersome, and is impractical if higher solution accuracy is needed.

The current state of the art for empirical loss minimization with strongly convex regularizers is randomized coordinate ascent on the dual objective — Stochastic Dual Coordinate Ascent (SDCA) . In contrast to primal SGD methods, the SDCA algorithm family is often preferred as it is free of learning-rate parameters, and has faster (geometric) convergence guarantees. This algorithm and its variants are increasingly used in practice . On the other hand, primal-only methods apply to a larger problem class, not only of form (1) that enables formation of dual problem (2) as considered here.

Another class of algorithms gaining attention in recent very few years are ‘variance reduced’ modifications of the original SGD algorithm. They are applied directly to the primal problem (1), but unlike SGD, have property that variance of estimates of the gradients tend to zero as they approach optimal solution. Algorithms such as SAG , SAGA and others come at the cost of extra memory requirements — they have to store a gradient for each training example. This can be addressed efficiently in the case of generalized linear models, but prohibits its use in more complicated models such as in deep learning. On the other hand, Stochastic Variance Reduced Gradient (SVRG) and its variants are often interpreted as ‘memory-free’ methods with variance reduction. However, these methods need to compute the full gradient occasionally to drive the variance reduction, which requires a full pass through the data and is an operation one generally tries to avoid. This and several other practical issues have been recently addressed in . Finally, another class of extensions to SGD are stochastic quasi-Newton methods . Despite their clear potential, a lack of theoretical understanding and complicated implementation issues compared to those above may still limit their adoption in the wider community. A stochastic dual Newton ascent (SDNA) method was proposed and analyzed in . However, the method needs to modified substantially before it can be implemented in a distributed environment.

For the empirical loss minimization problems of interest, stochastic subgradient descent (SGD) based methods are well-established. Several distributed variants of SGD have been proposed, many of which build on the idea of a parameter server . Despite their simplicity and accessibility in terms of implementation, the downside of this approach is that the amount of required communication is equal to the amount of data read locally, since one data point is accessed per machine per round (e.g., mini-batch SGD with a batch size of 1 per worker). These variants are in practice not competitive with the more communication-efficient methods considered in this work, which allow more local updates per communication round.

At the other extreme, there are distributed methods using only a single round of communication, such as . These methods require additional assumptions on the partitioning of the data, which are usually not satisfied in practice if the data are distributed “as is”, i.e., if we do not have the opportunity to distribute the data in a specific way beforehand. Furthermore, some cannot guarantee convergence rates beyond what could be achieved if we ignored data residing on all but a single computer, as shown in . Additional relevant lower bounds on the minimum number of communication rounds necessary for a given approximation quality are presented in .

Mini-batch methods (which instead of just one data-example use updates from several examples per iteration) are more flexible and lie within these two communication vs. computation extremes. However, mini-batch versions of both SGD and coordinate descent (CD) suffer from their convergence rate degrading towards the rate of batch gradient descent as the size of the mini-batch is increased. This follows because mini-batch updates are made based on the outdated previous parameter vector $\mathbf{w}$, in contrast to methods that allow immediate local updates like CoCoA.

Another disadvantage of mini-batch methods is that the aggregation parameter is harder to tune, as it can lie anywhere in the order of mini-batch size. The optimal choice is often either unknown, or difficult to compute. In the CoCoA setting, the parameter lies in the typically much smaller range given by $K$. In this work the aggregation parameter is further simplified and can be simply set to $1$, i.e., adding updates, which is achieved by formulating a more conservative local problem as described in Section 3.1.

With traditional batch gradient solvers not being competitive for the problem class (1), improved batch methods have also received much research attention recently, in the single machine case as well as in the distributed setting. In distributed environments, often used methods are the alternating direction method of multipliers (ADMM) as well as quasi-Newton methods such as L-BFGS, which can be attractive because of their relatively low communication requirements. Namely, communication is in the order of a constant number of vectors (the batch gradient information) per full pass through the data.

ADMM also comes with an additional penalty parameter balancing between the equality constraint on the primal variable vector $\mathbf{w}$ and the original optimization objective , which is typically hard to tune in many applications. Nevertheless, the method has been used for distributed SVM training in, e.g., . The known convergence rates for ADMM are weaker than the more problem-tailored methods mentioned we study here, and the choice of the penalty parameter is often unclear in practice.

Standard ADMM and quasi-Newton methods do not allow a gradual trade-off between communication and computation available here. An exception is the approach of Zhang, Lee and Shin , which is similar to our approach in spirit, albeit based on ADMM, in that they allow for the subproblems to be solved inexactly. However, this work focuses on L2-regularized problems and a few selected loss functions, and offers no complexity results.

Interestingly, our proposed CoCoA$\!{}^{\bf\textbf{\footnotesize+}}$ framework here – despite clearly aimed at cheap stochastic local solvers – does have similarities to block-wise variants of batch proximal methods as well, as explained as follows:

One of the main crucial differences of our proposed CoCoA$\!{}^{\bf\textbf{\footnotesize+}}$ framework compared to all known batch proximal methods (no matter if block-wise or not) is that the latter do require high accuracy subproblem solutions, and do not allow arbitrary solvers of weak accuracy $\Theta$ such as we do here, see also the next paragraph. Distributed Newton methods have been analyzed theoretically only when the subproblems are solved to high precision, see e.g. . This makes the local solvers very expensive and the convergence rates less general than in our framework (which allows weak local solvers). Furthermore, the analysis of requires additional strong assumptions on the data partitioning, such that the local Hessian approximations are consistent between the machines.

Developing distributed optimization methods that allow for arbitrary weak local optimizers requires carefully devising data-local subproblems to be solved after each communication round.

By making use of the primal-dual structure in the line of work of , the CoCoA and CoCoA$\!{}^{\bf\textbf{\footnotesize+}}$ frameworks proposed here are the first to allow the use of any local solver — of weak local approximation quality — in each round. Furthermore, the approach here also allows more control over the aggregation of updates between machines. The practical variant of the DisDCA Algorithm of , called DisDCA-p, also allows additive updates but is restricted to coordinate decent (CD) being the local solver, and was initially proposed without convergence guarantees. The work of has provided the first theoretical convergence analysis for an ideal case, when the distributed data parts are all orthogonal to each other — an unrealistic setting in practice. DisDCA-p can be recovered as a special case of the CoCoA$\!{}^{\bf\textbf{\footnotesize+}}$ framework when using CD as a local solver, if $|\mathcal{P}_{k}|=n/K$ and when using the conservative bound $\sigma^{\prime}:=K$, see also . The convergence theory presented here therefore also covers that method, and extends it to arbitrary local solvers.

Numerical Experiments

Our framework is related, but not identical, to running an inexact version of block coordinate ascent, applied to all block in parallel, and to the dual problem, where the level of inexactness is controlled by the parameter $\Theta$ through the use of a (possibly randomized) iterative “local” solver applied to the subproblems (local problems). For previous work on randomized block coordinate descent we refer to . See also .

In this section we explore numerous aspects of our distributed framework and demonstrate its competitive performance in practice. Section 6.1 first explores the impact of the local solver on overall performance, by comparing examples of various local solvers that can be used in the framework (the improved CoCoA$\!{}^{\bf\textbf{\footnotesize+}}$ framework as shown in Algorithms 1 and 2) as well as testing the effect of approximate solution quality. The results indicate that the choice of local solver can have a significant impact on overall performance. In Sections 6.2 and 6.3 we further explore framework parameters, looking at the impact of the aggregation parameter $\nu$ and the subproblem parameter $\sigma^{\prime}$, respectively. Finally, Section 6.5 demonstrates competitive practical performance of the overall framework on a large 280GB distributed dataset.

We conduct experiments on three datasets of moderate and large size, namely rcv1test, epsilon and splice-site.tThe datasets are available at http://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/.. The details of these datasets are listed in Table 6.

For solving subproblems, we compare numerous local solver methods, as listed in Table 6. Also, we apply Euclidean norm as regularizer $g(x)=\|x\|^{2}$ for all the experiments. All the algorithms are implemented in C++ with MPI, and experiments are run on a cluster of 4 Amazon EC2 m3.xlarge instances. Our open source code is available online at: https://github.com/optml/CoCoA.

In this section we compare the performance of our framework for various local solvers and various choices of inner iterations performed by a given local solver, resulting in different local accuracy measures $\Theta$. For simplicity, we choose the subproblem parameter $\sigma^{\prime}:=\nu K$ (see Lemma 3.2) as a simple obtainable and theoretically safe value for our framework.

1.2 Effect of the Quality of Local Solver Solutions on Overall Performance

Here we discuss how the quality of subproblem solutions affects the overall performance of Algorithm 2. In order to do so, we denote $H$ as the number of iterations the local solver is run for, within each communication round of the framework. We choose various values for $H$ on two local solvers, CD and L-BFGS , which gave the best performance in general. For CD, $H$ represents the number of local iterations performed on the subproblem. For L-BFGS, $H$ not only means the number of iterations, but also stands for the size of past information used to approximate the Hessian (i.e., the size of limited memory).

Looking at Figures 3 and 4, we see that for both local solver and all values of $\lambda$, increasing $H$ will lead to less iterations of Algorithm 2. Of course, increasing $H$ comes at the cost of the time spent on local solvers increasing. Hence, a larger value of $H$ is not always the optimal choice with respect to total elapsed time. For example, for the rcv_test dataset, when choosing CD to solve the subproblems, choosing $H$ to be $40,000$ uses less time and provides faster convergence. When using L-BFGS, $H=10$ seems to be the best choice.

2 Averaging vs. Adding the Local Updates

In this section, we compare the performance of our algorithm using two different schemes for aggregating partial updates: adding vs. averaging. This corresponds to comparing two extremes for the parameter $\nu$, either $\nu:=\frac{1}{K}$ (averaging partial solutions) or $\nu:=1$ (adding partial solutions). As discussed in Section 4, adding the local updates ($\nu=1$) will lead to less iterations than taking averaging, due to choosing different $\sigma^{\prime}$ in the subproblems. We verify this experimentally by considering several of the local solvers listed in Table 6.

We show results for RCV dataset, and we apply quadratic loss function with three different choices for the regularization parameter, $\lambda$=$1e-03$, $1e-04$, and $1e-05$. The experiments in Figures 6–11 indicate that the “adding strategy” will always lead to faster convergence than averaging, even though the difference is minimal when we apply a large number of iterations in the local solver. All the blue solid plots (adding) outperform the red dashed plots (averaging), which indicates the advantage of choosing $\nu=1$. Another note here is that for smaller $\lambda$, we will have to spend more iterations to get the same accuracy (because the original objective function (1) is less strongly convex).

In this section we consider the effect of the choice of the subproblem parameter on convergence (Figure 12). We plot duality gap over the number of communications for RCV and epsilon datasets with quadratic loss and set $K=8$, $\lambda=10^{-5}$. For $\nu=1$ (adding the local updates), we consider several different values of $\sigma^{\prime}$, ranging from $1$ to $8$. The value $\sigma^{\prime}=8$ represents the safe upper bound of $\nu K$, as given in Lemma 3.2.

Decreasing $\sigma^{\prime}$ improves performance in terms of communication until a certain point, after which the algorithm diverges. For the rcvtest dataset, the optimal convergence occurs around $\sigma^{\prime}=5$, and diverges fast for $\sigma^{\prime}\leq 3$. For epsilon dataset, $\sigma^{\prime}$ around $6$ is the best choice and the algorithm will not converge to optimal solution if $\sigma^{\prime}\leq 5.$ However, more importantly, the “safe” upper bound of $\sigma^{\prime}:=\nu K=8$ has only slightly worse performance than the practically best (but “un-safe”) value of $\sigma^{\prime}$.

4 Scaling Property

Here we demonstrate the ability of our framework to scale with $K$ (number of machines). We compare the run time to reach a specific tolerance on duality gap ($10^{-4}$ and $10^{-2}$) for two choices of $\nu$. Looking at Figure 13, we see that when choosing $\nu=1$, the performance improves as the number of machines increases. However, when $\nu=\frac{1}{K}$, the algorithm slows down as $K$ increases. The observations support our analysis in Section 4.

5 Performance on a Big Dataset

Also, the number of communication rounds for the three different loss functions are almost the same if we set all the other parameters to be same. However, the patterns of duality gap decrease for the three loss functions are different.

6 Comparison with other distributed methods

Finally, we compare the CoCoA$\!{}^{\bf\textbf{\footnotesize+}}$ framework with several competing distributed optimization algorithms. The DiSCO algorithm is a Newton-type method, where in each iteration the updates on iterates are computed inexactly using a Preconditioned Conjugate Gradients (PCG) method. As suggested in , in our implementation of DISCO we apply the Stochastic Average Gradient (SAG) method as the solver to get the initial solutions for each local machine and solve the linear system during PCG. DiSCO-F , improves on the computational efficiency of original DiSCO, based on data partitioned across features rather than examples. The DANE algorithm is another distributed Newton-type method, where in each iteration there are two rounds of communications. Also, a subproblem is to be solved in each iteration to obtain updates. For all these algorithms, all the hyper parameters are tuned manually to optimize their performance.

The experiments are conducted on a compute cluster with $K=4$ machines. We run all algorithms using two datasets. Since not all methods are primal-based in nature, it is difficult to continue using duality gap as a measure of optimality. Instead, the norm of the gradient of the primal objective function (1) is used to compare the relative quality of the solutions obtained.

As shown in Figure 15, in terns of number of communications, CoCoA$\!{}^{\bf\textbf{\footnotesize+}}$ usually converges more rapidly than competing methods during the early iterations, but tends to get slower later on in the iterative process. This illustrates that the Newton-type methods can accelerate in the vicinity of the optimal solution, as expected. However, CoCoA$\!{}^{\bf\textbf{\footnotesize+}}$ can still beat other methods in running time. The main reason for this is the fact that the subproblems in our framework can be solved more efficiently, compared with DiSCO and DANE.

We present a novel CoCoA$\!{}^{\bf\textbf{\footnotesize+}}$ framework for distributed optimization which enables fast and communication-efficient additive aggregation in distributed primal-dual optimization. We analyze the theoretical complexity of CoCoA$\!{}^{\bf\textbf{\footnotesize+}}$, giving strong primal-dual convergence rates with outer iterations scaling independently of the number of machines. We extended the basic theory to allow for non-smooth loss functions, arbitrary strongly convex regularizers, and primal-dual convergence results. Our experimental results show significant speedups in terms of runtime over previous methods, including the original CoCoA framework as well as other state-of-the-art methods.

Appendix A Proofs

Since $g$ is $1$-strongly convex, $g^{*}$ is $1$-smooth, and thus we can use (7) as follows

A.2 Proof of Lemma 3.1

where the first inequality follows from Jensen and the last equality also follows from the block diagonal definition of $G$ given in (11), i.e.

A.3 Proof of Lemma 3.2

Therefore we can conclude that $\nu\mathbf{h}^{T}X^{T}X\mathbf{h}\leq\nu K$ for all $\mathbf{h}$ included in the definition (10) of $\sigma^{\prime}_{\min}$, proving the claim.

A.4 Proofs of Theorems 4.2 and 4.3

Before we state the proofs of the main theorems, we will write and prove few crucial lemmas.

For sake of notation, we will write ${\boldsymbol{\alpha}}$ instead of ${\boldsymbol{\alpha}}^{t}$, $\mathbf{w}$ instead of $\mathbf{w}({\boldsymbol{\alpha}}^{t})$ and $\mathbf{u}$ instead of $\mathbf{u}^{t}$.

Now, let us estimate the expected change of the dual objective. Using the definition of the dual update ${\boldsymbol{\alpha}}^{t+1}:={\boldsymbol{\alpha}}^{t}+\nu\,\sum_{k}\mathbf{h}_{[k]}$ resulting in Algorithm 2, we have

Now, let us upper bound the $C$ term (we will denote by $\mathbf{h}^{\star}=\sum_{k=1}^{K}\mathbf{h}_{[k]}^{\star}$):

The convex conjugate maximal property implies that

Moreover, from the definition of the primal and dual optimization problems (1), (2), we can write the duality gap as

Now, the claimed improvement bound (21) follows by plugging (27) into (24). ∎

For general convex functions, the strong convexity parameter is $\gamma=0$, and hence the definition of $R^{t}$ becomes

At first let us estimate expected change of dual feasibility. By using the main Lemma A.1, we have

The choice of $s:=1$ and $t=t_{0}:=\max\{0,\lceil\frac{1}{\nu(1-\Theta)}\log(2\lambda n^{2}(D({\boldsymbol{\alpha}}^{\star})-D({\boldsymbol{\alpha}}^{0}))/(4L^{2}\sigma\sigma^{\prime}))\rceil\}$ will lead to

Clearly, (A.4.1) implies that (33) holds for $t=t_{0}$. Now imagine that it holds for any $t\geq t_{0}$ then we show that it also has to hold for $t+1$. Indeed, using

where in the last inequality we have used the fact that geometric mean is less or equal to arithmetic mean.

If $\overline{\boldsymbol{\alpha}}$ is defined as (16) then we obtain that

Now, if $T\geq\lceil\frac{1}{\nu(1-\Theta)}\rceil+T_{0}$ such that $T_{0}\geq t_{0}$ we obtain

To have right hand side of (38) smaller then $\epsilon_{Gap}$ it is sufficient to choose $T_{0}$ and $T$ such that

A.4.2 Proof of Theorem 4.2

into (41) we obtain that $\forall t:R^{t}\leq 0$. Putting the same $s$ into (21) will give us

Therefore if we denote by $\epsilon_{D}^{t}=D({\boldsymbol{\alpha}}^{\star})-D({\boldsymbol{\alpha}}^{t})$ we have that

The right hand side will be smaller than some $\epsilon_{D}$ if t ≥1ν(1-Θ) λγn+ σmaxσ’λγn log1ϵD. Moreover, to bound the duality gap, we have

Therefore ${Gap}({\boldsymbol{\alpha}}^{t})\leq\frac{1}{\nu(1-\Theta)}\frac{\lambda\gamma n+\sigma_{\max}\sigma^{\prime}}{\lambda\gamma n}\epsilon_{D}^{t}$. Hence if $\epsilon_{D}\leq\nu(1-\Theta)\frac{\lambda\gamma n}{\lambda\gamma n+\sigma_{\max}\sigma^{\prime}}\epsilon_{Gap}$ then ${Gap}({\boldsymbol{\alpha}}^{t})\leq\epsilon_{Gap}$. Therefore after t ≥1ν(1-Θ) λγn+ σmaxσ’λγn log( 1ν(1-Θ) λγn+ σmaxσ’λγn 1ϵGap ) iterations we have obtained a duality gap less than $\epsilon_{Gap}$.