Distributed optimization with arbitrary local solvers

With the growth of data and necessity for distributed optimization methods, solvers that work well on a single machine must be re-designed to leverage distributed computation. Recent work in this area has been limited by focusing heavily on developing highly specific methods for the distributed environment. These special-purpose methods are often unable to fully leverage the competitive performance of their well-tuned and customized single machine counterparts. Further, they are unable to easily integrate improvements that continue to be made to single machine methods. To this end, we present a framework for distributed optimization that both allows the flexibility of arbitrary solvers to be used on each (single) machine locally and yet maintains competitive performance against other state-of-the-art special-purpose distributed methods. We give strong primal–dual convergence rate guarantees for our framework that hold for arbitrary local solvers. We demonstrate the impact of local solver selection both theoretically and in an extensive experimental comparison. Finally, we provide thorough implementation details for our framework, highlighting areas for practical performance gains.


Motivation
Regression and classification techniques, represented in the general class of regularized loss minimization problems [67], 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 [3,7,12,18,22,28,32,34,37,46,50,58,60,62,76].
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 [26,70].
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 [26,35,70].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.

Contributions
Reusability of Existing Local Solvers.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.
Adaptivity to Communication Cost.On real-world compute systems, the cost of communication versus computation typical varies by many orders of magnitude, from highperformance 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.
Strong Theoretical Guarantees.In this paper we extend and improve upon the Co-CoA [26] method.Our theoretical convergence rates apply to both smooth and nonsmooth losses, and for both CoCoA as well as the more general CoCoA + 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 [1,58], which can severely affect overall runtime.
Primal-Dual Convergence.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.
Experimental Results.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 + 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.

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 A on a single machine, its performance is typically measured by the total runtime: TIME(A) = I A ( ) × T A . (T-A) Here, T A stands for the time it takes to perform a single iteration of algorithm A, and I A ( ) is the number of iterations A needs to attain an -accurate objective. 1n 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 [37]), where The extra term c is the time required to perform one round of communication. 2 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 T 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 [26,35,61,[70][71][72] 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 [26,70] 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 [26], 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 [35,61].
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 Θ ∈ [0, 1], where Θ = 0 describes an exact solution of the subproblem, and Θ = 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 T A (Θ) represents the time the local algorithm A needs to obtain accuracy Θ on the local subproblem.Note that the number of outer iterations I( , Θ) is independent of choice of the inner algorithm A, which will also be reflected by our convergence analysis presented in Section 4. Our convergence rates will hold for any local solver A achieving local quality Θ.For strongly convex problems, the general form will be I( , . The inverse dependence on 1 − Θ suggests that there is a limit to how much we can gain by solving local subproblems to high accuracy -Θ close to 0. There will always be order of log(1/ ) outer iterations needed.Hence, excessive local accuracy should not be necessary.On the other hand, if Θ → 1, meaning that the cost

Smooth
and quality of the local solver diminishes, then the number of rounds I( , Θ) 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 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 descent 3 .Running gradient descent for a single iteration would happen to attain a particular value of Θ.Note that we typically do not set the value Θ 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 Θ such that T A (Θ) = O(1) × 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 Θ, both theoretically and empirically, later in Sections 4 and 6.

Problem Formulation
Let the training data {x i ∈ R d , y i ∈ R} n i=1 be the set of input-output pairs, where y i can be real valued, but also from a discrete set in the case of classification problems.We will assume without loss of generality that ∀i : x i ≤ 1.Many common tasks in machine learning and signal processing can be cast as the following optimization problem: where i is some convex loss function and g is a regularizer.Note that y i is typically hidden in the formulation of functions i .Table 1 lists several common loss functions together with their convex conjugates * i [57].The dual optimization problem for formulation (1) -as a special case of Fenchel duality -can be written as follows [57,73]: where X = [x 1 , x 2 , . . ., x n ] ∈ R d×n , and * i and g * are the convex conjugate functions of i and g, respectively.The convex (Fenchel) conjugate of a function φ : R k → R is defined as the function φ * : R k → R, with φ * (u) := sup s∈R k {s T u − φ(s)}.
For simplicity throughout the paper, let us denote such that D(α) It is well known [19,48,57,63] that the first-order optimality conditions give rise to a natural mapping that relates pairs of primal and dual variables.The mapping employs the linear map given by the data X, and maps any dual variable α ∈ R n to a primal candidate vector w ∈ R d as follows: λn Xα , where we denote v(α) := 1 λn Xα .For this mapping, under the assumptions that we make in Section 2.2 below, it holds that if α is an optimal solution of (2), then w(α ) is an optimal solution of (1).In particular, strong duality holds between the primal and dual problems.If we define the duality gap function as then Gap(α ) = 0, which ensures that by solving the dual problem (2) we also solve the original primal problem of interest (1).As we will later see, there are many benefits to leveraging this primal-dual relationship, including the ability to use the duality gap as a certificate of solution quality, and, in the distributed setting, as a tool through which we can effectively distribute computation.
Notation.We assume that to solve problem (2) we have a network of K machines at our disposal.The data {x i , y i } n i=1 is residing on the K machines in a distributed way, with every machine only holding a part of the whole dataset.In the same way we split the dual variables α i , with each corresponding to an individual data point x i .The given data distribution is described using a partition P 1 , . . ., P K that corresponds to the indices of data and dual variables residing on machine k.Formally, P k ⊆ {1, 2, . . ., n} for each k, P k ∩ P l = ∅ whenever k = l, and K k=1 P k = {1, 2, . . ., n}.In order to efficiently use this structure in the text, we introduce the following notation.For any h ∈ R n , let h [k] be the vector in R n defined in such a way that (h . Analogously, we write X [k] for the matrix consisting only of the columns i ∈ P k , padded with zeros in all other columns.

Technical Assumptions
Here we first state the properties and assumptions used throughout the paper.We assume that for all i ∈ {1, . . ., n}, the function i is convex, i.e., ∀λ ∈ [0, 1] and ∀x, y ∈ R we have We also assume the function g to be 1-strongly convex, i.e., for all w, u ∈ R d it holds that g(w + u) ≥ g(w) + ∇g(w), u + 1 2 u 2 , where ∇g(w) is any subgradient4 of the function g.Here, • denotes the standard Euclidean norm.
Note that we use subgradients in the definition of strong convexity.This is due to the fact that while we will need the function g to be strongly convex in our analysis, we do not require smoothness.An example used in practice is g(w) = w 2 + λ w 1 for some λ ∈ R. Also note that in the problem formulation (1) we have a regularization parameter λ, which controls the strong convexity parameter of the entire second term.Hence, fixing the strong convexity parameter of g to 1 is not restrictive in this regard.For instance, this setting has been used previously in [13,48,57].
The following assumptions state properties of the functions i , which we use only in certain results in the paper.We always explicitly state when we require each assumption.Assumption 2.1 ((1/γ)-Smoothness) Functions i : R → R are 1/γ-smooth, if ∀i ∈ {1, . . ., n} and ∀x, h ∈ R it holds that where ∇ i (x) denotes the gradient of the function i .
Assumption 2.2 (L-Lipschitz Continuity) Functions i : R → R are L-Lipschitz continuous, if ∀i ∈ {1, . . ., n} and ∀x, h ∈ R it holds that Remark 1 As a consequence of having 1/γ-smoothness of i and 1-strong convexity of g, we have that the functions * i (•) are γ-strongly convex and g * (•) is 1-smooth [53].These are the properties we will ultimately use as we will be solving the dual problem (2).Note that 1-smoothness of g * : R d → R means that for all x, h ∈ R d , Following lemma, a consequence of 1-smoothness of g * and the definition of f , will be crucial in deriving a meaningful local subproblem for the proposed distributed framework.
Lemma 2.3 Let f be defined in (3).Then for all α, h ∈ R n we have Remark 2 Note that although the above inequality appears as a consequence of the problem structure (2) and of strong convexity of g, there are other ways to satisfy it.Hence, our dual analysis holds for all optimization problems of the form max α D(α), where D(α) = −f (α) − R(α), and where f satisfies inequality (8).However, for the duality gap analysis we naturally do require (8) that the dual problem arises from the primal problem with g being strongly convex.

The Framework
In this section we start by giving a general view of the proposed framework, explaining the most important concepts needed to make the framework efficient.In Section 3.1 we discuss the formulation of the local subproblems, and in Section 3.2 specific details and best practices for implementation.
The data distribution plays a crucial role in Algorithm 1, where in each outer iteration indexed by t, machine k runs an arbitrary local solver on a problem described only by the data that particular machine owns and other fixed constants or linear functions.
The crucial property is that the optimization algorithm on machine k changes only coordinates of the dual optimization variable α t corresponding to the partition P k to obtain an approximate solution to the local subproblem.We will formally specify this in Assumption 4.1.After each such step, updates from all machines are aggregated to form a new iterate α t+1 .The aggregation parameter ν will typically be between ν = 1/K, corresponding to averaging, and ν = 1, to adding.
for k ∈ {1, 2, . . ., K} in parallel over machines do 4: end for 6: 7: end for Here we list the core conceptual properties of Algorithm 1, which are important qualities that allow it to run efficiently.
Locality.The local subproblem G k (LO) is defined purely based on the data points residing on machine k, as well as a single shared vector in R d (representing the state of the α t variables of the other machines).Each local solver can then run independently and in parallel, i.e., there is no need for communication while solving the local subproblems.Local changes.The optimization algorithm used to solve the local subproblem G k outputs a vector h t [k] with nonzero elements only in coordinates corresponding to variables α [k] stored locally (i.e., i ∈ P k ).Efficient maintenance.Given the description of the local problem G k ( • ; α t ) at time t, the new local problem G k ( • ; α t+1 ) at time t + 1 can be formed on each machine, requiring only communication of a single vector in R d from each machine k to the master node, and vice versa, back to each machine 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 α 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 α 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 G k in the following section.

The Local Subproblems
We can define a data-local subproblem of the original dual optimization problem (2), which can be solved on machine k and only requires accessing data which is already available locally, i.e., datapoints with i ∈ P k .More formally, each machine k is assigned the following local subproblem, depending only on the previous shared primal vector w ∈ R d , and the change in the local dual variables α i with i We are now ready to define the local objective G σ k ( • ; α) as follows: The role of the parameter σ ≥ 1 is to measure the "difficulty" of the data partition, in a sense which we will discuss in detail in Section 3.3 below.
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 h [k] expresses the update proposed by machine k.In this spirit, note also that the approximation coincides with D at the reference point α, i.e.K k=1 G σ k (0; α) = D(α).We will discuss the interpretation and properties of these subproblems in more detail below in Section 3.3.

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.
Recall that the local subproblem objective G σ k ( • ; α) was defined in (LO).We will now equivalently rewrite this optimization problem, to clarify how it is expressed only using local information.To do so, we use our simplifying notation v = v(α) := 1  λn Xα for given α.As we see in the reformulation, it is precisely this vector v ∈ R d which contains all the necessary shared information between the machines.Given the vector v, the subproblem (LO) writes equivalently as Here for the reformulation of the gradient term, we have simply used the chain rule on the objective f (recall the definition f (α Practical Distributed Framework.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 α [k] ).No information about the other machines variables α or their data is necessary.The only requirement for the method to work is that between the rounds, the changes in the α [k] variables on each machine and the resulting global change in v are kept consistent, in the sense that v t = v(α t ) := 1 λn Xα t must always hold.Note that for the evaluation of ∇g * (v), the vector 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 worker machine 1 stored data ↵ [1] worker machine K stored data ↵ [1] worker machine K stored data The first two iterations of the improved framework (practical implementation).
In the following more detailed formulation of the CoCoA + framework shown in Algorithm 2 (equivalent reformulation of Algorithm 1), the crucial communication pattern of the framework finally becomes more clear: Per round, only a single vector (the update on v ∈ R d ) needs to be sent over the communication network.The reduce-all operation in line 10 means that each machine sends their vector ∆v t k ∈ R d to the network, which performs the addition operation of the K vectors to the old v t .The resulting vector v t+1 is then communicated back to all machines, so that all have the same copy of v t+1 before the beginning of the next round.
The framework as shown below in Algorithm 2 clearly maintains the consistency of α t and v t = v t (α 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.
Algorithm 2 Improved CoCoA+ Framework, Practical Implementation for k ∈ {1, 2, . . ., K} in parallel over machines do 5: end for 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 ν (controlling the level of adding versus averaging the resulting updates from each machine) and σ (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 σ is to measure the difficulty of the given data partition.For the convergence results discussed below to hold, σ 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 {P k } K k=1 , i.e., In this notation, it is easy to see that the crucial quantity defining σ min above is written as The following lemma shows that if the aggregation and subproblem parameters ν and σ satisfy (10), then the sum of the subproblems k G σ k will closely approximate the global objective function D.More precisely, this sum is a block-separable lower bound on D.
The following lemma gives a simple choice for the subproblem parameter σ , which is trivial to calculate for all values of the aggregation parameter ν ∈ R, and safe in the sense of the desired condition (10) above.Later we will show experimentally (Section 6) that the choice of this safe upper bound for σ only has a minimal effect on the overall performance of the algorithm.Lemma 3.2 For any aggregation parameter ν ∈ [0, 1], the choice of the subproblem parameter σ := νK is valid for (10), i.e., νK ≥ σ min .

Main Results
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.

Quality of Local Solutions
The notion of approximation quality provided by the local solvers is measured according to the following: Assumption 4.1 (Quality of local solution) Let Θ ∈ [0, 1) and α ∈ R n be fixed, and let h [k] be the optimal solution of a local subproblem G k ( • ; α).We assume the local optimization procedure run on every node k ∈ [K], in each iteration t produces a (possibly The assumption specifies the (relative) accuracy Θ obtained on solving the local subproblem G k .Considering the two extreme examples, setting Θ = 0 would require to find the exact maximum, while Θ = 1 states that no improvement was achieved at all by the local solver.Intuitively, we would prefer Θ to be small, but spending many computational resources to drive Θ to 0 can be excessive in practice, since 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 Θ 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 Θ ∈ [0, 1] 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).

Complexity Bounds
Now we are ready to state the main results.Theorem 4.2 covers the case when ∀i the loss function i is 1/γ smooth and Theorem 4.3 covers the case when i is L-Lipschitz continuous.For simplicity in the rates, we define the following two quantities: Furthermore, after T iterations with we have the expected duality gap Theorem 4.3 (Lipschitz continuous loss functions) Consider Algorithm 2 with Assumption 4.1.Let i (•) be L-Lipschitz continuous, and Gap > 0 be the desired duality gap (and hence an upper-bound on primal sub-optimality).Then after T iterations, where ) , we have that the expected duality gap satisfies at the averaged iterate 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 − Θ suggests that it is worth pushing the rate of local accuracy Θ 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 [26]) by more aggressive adding instead, that is ν = 1 instead of ν = 1/K.Note that the safe subproblem parameter for the averaging case (ν = 1/K) is σ := 1, while for adding (ν = 1) it is given by σ := 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 ν and σ for the two cases, as we will illustrate more clearly in the next subsection.

Discussion and Interpretations of Convergence Results
As the above theorems suggest, it is not possible to meaningfully change the aggregation parameter ν in isolation.It comes naturally coupled with a particular subproblem.
In this section, we explain a simple way to be able to have the aggregation parameter ν = 1, that is to aggressively add up the updates from each machine.The motivation for this comes from a common practical setting.When solving the SVM dual (Hinge loss: i (a) = max{0, y i − a}), the optimization problem comes with "box constraints", i.e., for all i ∈ {1, . . ., n}, we have α i ∈ [0, 1] (see Table 1).The particular values of α i being 0 or 1 have a particular interpretation in the context of original problem (1).If we used ν < 1, we would never be able reach the upper boundary of any variable α i , when starting the algorithm at 0. This example illustrates some of the downsides of averaging vs adding updates, coming from the fact that the step-size from using averaging (by being 1/K times shorter) can result in 1/K times slower convergence.
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: On the other hand, if we would just average results (as proposed in [26]), we would obtain following corollary: iterations, we have Comparing the leading terms in Equations ( 17) and (18) we see that the leading term for the ν = 1 choice is O(λγn + σ max K), which is always better than for the ν = 1/K case, when the leading term is O(λγn + σ max K).This strongly suggests that adding in Framework 2 is preferable, especially when λγn σ 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 λ.When the regularizer is O(1/n), the difference is negligible.This behaviour is also present in practice, as we will point out in Section 6.

Discussion and Related Work
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.
Single-Machine Solvers.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 [52].We refer the reader to [8,[39][40][41] 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 hyperparameters 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) [56].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 [57,68].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 [54], SAGA [16] and others [17,55] 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 [27,29,30,42,69] 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 [2].Finally, another class of extensions to SGD are stochastic quasi-Newton methods [6,11].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 [47].However, the method needs to modified substantially before it can be implemented in a distributed environment.
SGD-based Algorithms.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 [18,43,50].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.
One-Shot Communication Schemes.At the other extreme, there are distributed methods using only a single round of communication, such as [24,36,38,75,78].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 [60].Additional relevant lower bounds on the minimum number of communication rounds necessary for a given approximation quality are presented in [1,3].
Mini-Batch Methods.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) [14,15,37,46,48,50,51,57,66,70] 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 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.
Distributed Batch Solvers.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) [9] 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 w and the original optimization objective [9], which is typically hard to tune in many applications.Nevertheless, the method has been used for distributed SVM training in, e.g., [23].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 [74], 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 + 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: The purpose of our subproblems as defined in (LO) is to form a data-dependent blockseparable quadratic approximation to the smooth part of the original (dual) objective (2), while leaving the non-smooth part R intact (recall that R(α) was defined to collect the * i functions, and is separable over the coordinate blocks).Now if hypothetically each of our regularized quadratic subproblems (LO) were to be minimized exactly, the resulting steps could be interpreted as block-wise proximal Newton-type steps on each coordinate block k of the dual (2), where the Newton-subproblem is modified to also contain the proximal part R.This connection only holds for the special case of adding (ν = 1), and would correspond to a carefully adapted step-size in the block-wise Newton case.
One of the main crucial differences of our proposed CoCoA + 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 Θ 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.[60].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 [60] requires additional strong assumptions on the data partitioning, such that the local Hessian approximations are consistent between the machines.Distributed Methods Allowing Local Optimization.Developing distributed optimization methods that allow for arbitrary weak local optimizers requires carefully devising datalocal subproblems to be solved after each communication round.
By making use of the primal-dual structure in the line of work of [31,45,[70][71][72], the CoCoA and CoCoA + 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 [70], 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 [71] 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.DisDCAp can be recovered as a special case of the CoCoA + framework when using CD as a local solver, if |P k | = n/K and when using the conservative bound σ := K, see also [31,35].The convergence theory presented here therefore also covers that method, and extends it to arbitrary local solvers.
Inexact Block Coordinate Descent.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 Θ 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 [65].See also [64].

Numerical Experiments
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 + 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 ν and the subproblem parameter σ , 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.t 5.The details of these datasets are listed in Table 2.For solving subproblems, we compare numerous local solver methods, as listed in Table 3.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.

Exploration of Local Solvers within the Framework
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 Θ.For simplicity, we choose the subproblem parameter σ := νK (see Lemma 3.2) as a simple obtainable and theoretically safe value for our framework.

Comparison of Different Local Solvers
Here we compare the performance of the various local solvers listed in Table 3.We here show results for quadratic loss function i (a) = 1 2 (a−y i ) 2 with three different values of the regularization parameter, λ=10 −3 , 10 −4 , and 10 −5 , for g(.) being the default Euclidean squared norm regularizer.The dataset is RCV and we ran the framework for a maximum of T := 100 communication rounds.We set ν = 1 (adding) and choose H which gave the best performance in CPU time (see Table 4) for each solver.From Figure 2, we find that the coordinate descent (CD) local solver always outperforms the other solvers, even though it may be slower than L-BFGS at the beginning.The reason for this is that CD, as compared to the other methods, does not need to spend time evaluating the full (batch) gradient and function values.Also note that some of the solvers cannot guarantee strict decrease of the duality gap, and sometimes this fluctuation can be very dramatic.

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 [49,56] and L-BFGS [10], 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 λ, 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.

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 ν, either ν := 1 K (averaging partial solutions) or ν := 1 (adding partial solutions).As discussed in Section 4, adding the local updates (ν = 1) will lead to less iterations than taking averaging, due to choosing different σ in the subproblems.We verify this experimentally by considering several of the local solvers listed in Table 3.

The Effect of the Subproblem Parameter σ
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, λ = 10 −5 .For ν = 1 (adding the local updates), we consider several different values of σ , ranging from 1 to 8. The value σ = 8 represents the safe upper bound of νK, as given in Lemma 3.2.
Decreasing σ improves performance in terms of communication until a certain point, after which the algorithm diverges.For the rcvtest dataset, the optimal convergence occurs around σ = 5, and diverges fast for σ ≤ 3.For epsilon dataset, σ around 6 is the best choice and the algorithm will not converge to optimal solution if σ ≤ 5.However, more importantly, the "safe" upper bound of σ := νK = 8 has only slightly worse performance than the practically best (but "un-safe") value of σ .Figure 12.The effect of σ on convergence for the rcvtest and epsilon datasets distributed across 8 machines.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 ν.Looking at Figure 13, we see that when choosing ν = 1, the performance improves as the number of machines increases.However, when ν = 1 K , the algorithm slows down as K increases.The observations support our analysis in Section 4.

Performance on a Big Dataset
As shown in Figure 14, we test the algorithm on the splice-site.tdataset, whose size is about 280 GB.We show experiments for three different loss functions , namely logistic loss, hinge loss and least squares loss (see Table 1).We set λ = 10 −6 for the squared norm regularizer.The dataset is distributed across K = 4 machines and we use CD as the local solver with H = 50, 000.In all the cases, an optimal solution can be reached in about 20 minutes and again, we observe that setting the aggregation parameter ν := 1 leads to faster convergence than ν := 1 K (averaging).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.

Comparison with other distributed methods
Finally, we compare the CoCoA + framework with several competing distributed optimization algorithms.The DiSCO algorithm [77] 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 [77], in our implementation of DISCO we apply the Stochastic Average Gradient (SAG) method [54] as the solver to get the initial solutions for each local machine and solve the linear system during PCG.DiSCO-F [33], improves on the computational efficiency of original DiSCO, based on data partitioned across features rather than examples.The DANE algorithm [59] 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 + 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 + By taking any h ∈ R n for which h T Gh ≤ 1, in view of (10), we get Therefore we can conclude that νh T X T Xh ≤ νK for all h included in the definition (10) of σ 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.
Lemma A.1 Let * i be strongly6 convex with convexity parameter γ ≥ 0 with respect to the norm • , ∀i ∈ [n].Then for all iterations t of Algorithm 1 under Assumption 4.1, and any s ∈ [0, 1], it holds that where Proof.For sake of notation, we will write α instead of α t , w instead of w(α t ) and u instead of ≤ E D(α) ≤ ν D(α) where Clearly, (A14) implies that (A15) holds for t = t 0 .Now imagine that it holds for any t ≥ t 0 then we show that it also has to hold for t + 1.Indeed, using we obtain Now, we will upperbound D as follows where in the last inequality we have used the fact that geometric mean is less or equal to arithmetic mean.
If α is defined as ( 16) then we obtain that

10 :
reduce all to compute v t+1 := v t + ν K k=1 ∆v t k communication 11: end for 3.3 Compatibility of the Subproblems for Aggregating Updates

Figure 2 .
Figure 2. Performance of different local solvers.

Figure 3 .
Figure 3. Varying the number of iterations of CD as a local solver.

Figure 4 .
Figure 4. Varying the number of iterations of L-BFGS as a local solver.

Figure 5 .
Figure 5. Adding (blue solid line) vs Averaging (red dashed line) for CD as the local solver.

Figure 6 .
Figure 6.Adding (blue solid line) vs Averaging (red dashed line) for APPROX as the local solver.

Figure 7 .
Figure 7. Adding (blue solid line) vs Averaging (red dashed line) for Gradient Descent as the local solver.

Figure 8 .
Figure 8. Adding (blue solid line) vs Averaging (red dashed line) for L-BFGS as the local solver.

Figure 9 .
Figure 9. Adding (blue solid line) vs Averaging (red dashed line) for Conjugate Gradient Method as the local solver.

Figure 10 .
Figure 10.Adding (blue solid line) vs Averaging (red dashed line) for BB as the local solver.

Figure 11 .
Figure 11.Adding (blue solid line) vs Averaging (red dashed line) for FISTA as the local solver.

Figure 13 .
Figure 13.The effect of increasing the number of machines K on the time (s) to reach a solution with expected duality gap.

Figure 14 .
Figure 14.Performance of Algorithm 2 on splice-site.tdataset, with three different loss functions.

Table 1 .
Examples of commonly used loss functions.

Table 2 .
Datasets used for numerical experiments.

Table 3 .
Local solvers used in numerical experiments.

Table 4 .
Optimal H for different local solvers for RCV dataset.