L1-Regularized Distributed Optimization: A Communication-Efficient Primal-Dual Framework

Introduction

In this paper, we consider standard regularized loss minimization problems, including as our main focus $L_{1}$-regularized optimization problems of the form

One promising distributed method is CoCoA$\!{}^{\bf\textbf{+}}$ , a recently proposed primal-dual framework that demonstrates competitive performance, provides a flexible communication scheme, and enables the use of off-the-shelf single-machine solvers internally. However, by solving the problem in the dual, CoCoA$\!{}^{\bf\textbf{+}}$ (like SDCA, prox-SDCA, and numerous other primal-dual methods ) is only equipped to handle strongly convex regularizers, which prevents it from being directly applied to $L_{1}$-regularized objectives. Moreover, by requiring the data to be distributed by data point rather than by feature, communication can become a prohibitive bottleneck for CoCoA$\!{}^{\bf\textbf{+}}$ as the number of features grows large, which is precisely the setting of interest for $L_{1}$ regularization.

In this work, we take a different perspective and propose a framework that can run either in the dual, or on the primal directly. From this change in perspective we derive several new primal-dual distributed optimization methods, in particular for sparsity-inducing regularizers. Our approach uses ideas from CoCoA$\!{}^{\bf\textbf{+}}$, though leveraging these ideas in this new setting requires significant theoretical and algorithmic modifications, particularly in handling non-strongly convex regularizers. The proposed primal-dual framework and associated rates are novel contributions even in the non-distributed case.

By building on the CoCoA$\!{}^{\bf\textbf{+}}$ framework, proxCoCoA$\!{}^{\bf\textbf{+}}$ comes with several benefits, including the use of arbitrary local solvers on each machine, and the analysis of and ability to solve subproblems to arbitrary accuracies. However in contrast to CoCoA$\!{}^{\bf\textbf{+}}$, we consider a much broader class of optimization problems. This results in a more general framework that: (1) specifically incorporates the case of $L_{1}$ regularization; (2) allows for the flexibility of distributing the data by either feature or data point; and (3) can be run on either the primal or dual formulation, which we show to have significant theoretical and practical implications.

We derive convergence rates for the general class of problems considered in this work, leveraging a novel approach in the analysis of primal-dual rates for non-strongly convex regularizers. The proposed technique is a significant improvement over simple smoothing techniques used in, e.g., that enforce strong convexity by adding a small $L_{2}$ term to the objective. Our results include primal-dual rates and certificates for both strongly convex and non-strongly convex regularizers and losses, and we show how earlier rates of CoCoA and CoCoA$\!{}^{\bf\textbf{+}}$ can be derived as a special case of our new rates / methods.

The proposed framework yields order-of-magnitude speedups (as much as 50$\times$ faster) as compared to other state-of-the-art methods for $L_{1}$-regularized optimization. We demonstrate these performance gains in an extensive experimental comparison on real-world distributed datasets. We additionally show significant improvements over CoCoA$\!{}^{\bf\textbf{+}}$ when considering strongly convex objectives. All algorithms for comparison are implemented in Apache Spark and run on Amazon EC2 clusters. Our code is available at: github.com/gingsmith/proxcocoa.

Setup

A great variety of methods in machine learning and signal processing are posed as the minimization of a weighted sum of two convex functions, where the first term is a convex function of a linear predictor and the second term is a regularizer:

The above setting encompasses all convex loss functions depending on linear predictors ${\bf y}_{j}^{T}{\boldsymbol{\alpha}}$, together with most common convex regularizers, including all separable functions, such as $L_{1}$- or general $L_{p}$-norms, or the elastic net given by $\frac{\eta}{2}\left\lVert{\cdot}\right\rVert_{2}+(1-\eta)\left\lVert{\cdot}\right\rVert_{1}$.

The proxCoCoA++\!{}^{\bf\textbf{+}} Algorithmic Framework

The proxCoCoA$\!{}^{\bf\textbf{+}}$ framework is given in Algorithm 1. This framework builds on the recent CoCoA$\!{}^{\bf\textbf{+}}$ framework , though with a more general objective, a modified subproblem, and where we allow the method to be applied to either the primal or dual formulation. To distribute the method, we assign each machine to work only on local coordinates of the weight vector ${\boldsymbol{\alpha}}$, and access only data that is stored locally. Machines share state through the vector ${\bf v}:=A{\boldsymbol{\alpha}}$. This vector is communicated at each round after using local solvers in parallel to find (possibly) approximate solutions to the subproblems defined in (2). Solving the primal problem (A) directly with proxCoCoA$\!{}^{\bf\textbf{+}}$ will result in distributing the data column-wise (by feature), and having the vector ${\bf v}$ be of length equal to the number of data points. This can greatly reduce communication costs as the number of features grows (see Section 6). Most importantly, the proposed setup will prepare us to handle non-strongly convex regularizers in both theory and practice, as we further explain in the following sections.

For each machine, we define a data-local subproblem of the original optimization problem (A). This simpler problem can be solved on machine $k$ and only requires accessing data which is already available locally, i.e., columns $A_{i}$ such that $i\in\mathcal{P}_{k}$. The subproblem depends only on the previous shared vector ${\bf v}:=A{\boldsymbol{\alpha}}$ and the local data:

with ${\bf w}:=\nabla f({\bf v})$. We denote the change of local variables $\alpha_{i}$ for indices $i\in\mathcal{P}_{k}$ as $\Delta{\boldsymbol{\alpha}}_{[k]}$. For a given aggregation parameter $\gamma\in(0,1]$, the subproblem relaxation parameter $\sigma^{\prime}$ will be set as $\sigma^{\prime}:=\gamma K$, but can also be improved in a data-dependent way as we discuss in Appendix E.

Our local subproblems have the appealing property of being very similar in structure to the global problem (A), with the main difference being that they are defined on a smaller (local) subset of the data. For the user of our framework, this presents a major advantage in that existing single machine-solvers can be directly re-used in our distributed framework (Algorithm 1) by employing them on the subproblems $\mathcal{G}^{\sigma^{\prime}}_{k}$. Therefore, problem-specific tuned solvers which have already been developed, along with associated speed improvements (such as multi-core implementations), can be easily leveraged in the distributed setting. We quantify the dependence on local solver performance in more detail in our convergence analysis (Section 4).

The above definition of the local objective functions $\mathcal{G}^{\sigma^{\prime}}_{k}$ are such that they closely approximate the global objective in (A) as the “local” variable $\Delta{\boldsymbol{\alpha}}_{[k]}$ varies, which we will see in the analysis (Lemma 8 in the appendix). In fact, if the subproblem were solved exactly, this could be interpreted as a data-dependent, block-separable proximal step, applied to the $f$ part of the objective (A) as follows:

where $L=\sum_{i\in[n]}g_{i}({\boldsymbol{\alpha}}_{i}+\Delta{\boldsymbol{\alpha}}_{i})\,.$

However, note that in contrast to traditional proximal methods, our algorithm does not assume that the prox subproblems be solved to high accuracy, as we instead allow the use of local solvers of any approximation quality $\Theta$. This notion is made precise with the following assumption.

We assume that there exists $\Theta\in[0,1)$ such that $\forall k\in[K]$, the local solver at any outer iteration $t$ produces a (possibly) randomized approximate solution $\Delta{\boldsymbol{\alpha}}_{[k]}$, which satisfies

In practice, the time spent solving the local subproblems in parallel should be chosen comparable to the required time of a communication round, for best overall efficiency on a given system. We study this trade-off both in theory (Section 4) and experiments (Section 6).

1 Primal-Dual Context

Exploiting primal-dual structure is not a requirement to optimize (A); indeed, we have shown above how to solve this optimization problem directly. However, noting the relationship between primal and dual objectives has many benefits, including computation of the duality gap, which allows us to have a certificate of approximation quality. It is also useful as an analysis tool and helps relate this work to the prior work of . To leverage this structure, starting from our original formulation (A) with objective function $\mathcal{D}({\boldsymbol{\alpha}}):=f(A{\boldsymbol{\alpha}})+\sum_{i=1}^{n}g_{i}(\alpha_{i})$, the dual problem is given by

acts as a certificate of approximation quality, as the distance to the true optimum $\mathcal{P}({\bf w}^{\star})$ is always bounded above by the duality gap. A globally defined and finite duality gap $G({\boldsymbol{\alpha}})$ for any problem (A) can be obtained by bounding the region of interest for the iterates ${\boldsymbol{\alpha}}$. This “Lipschitzing” trick will make the conjugates $g^{*}_{i}$ globally defined and Lipschitz , as we prove in Section 4.

Previous work of CoCoA$\!{}^{\bf\textbf{+}}$ mapped machine learning tasks to $\mathcal{P}({\bf w})$ (B), and then solved this problem in the dual. While this can still be accomplished with the machinery of proxCoCoA$\!{}^{\bf\textbf{+}}$ (see Section F), here our main focus is to instead solve the original objective $\mathcal{D}({\boldsymbol{\alpha}})$ (A) directly. This can have a large practical impact for the described applications in the distributed setting, as it implies that we can distribute the data by feature rather than by data point. Further, we will communicate a vector equal in size to the number of data points, as opposed to the number of features. When the number of features is high (as is common in sparsity-inducing models) this can significantly reduce communication and improve overall performance, as we demonstrate in Section 6. Further, it allows us to directly leverage state-of-the-art coordinate-wise primal methods, such as glmnet and extensions . From a theoretical perspective, solving $\mathcal{D}({\boldsymbol{\alpha}})$ will allow us to consider non-strongly convex regularizers, which were not covered in CoCoA$\!{}^{\bf\textbf{+}}$, as we discuss in Section 4.

Convergence Analysis

In this section we provide convergence rates for the proposed framework, and introduce an important theoretical technique in analyzing non-strongly convex terms in the primal-dual setting. For simplicity of presentation, we assume in the analysis that the data partition is balanced; i.e., $n_{k}=n/K$ for all $k$. Furthermore, we assume that the columns of A satisfy $\|{\bf x}_{i}\|\leq 1$ for all $i\in[n]$. We present results for the case where $\gamma:=1$ in Algorithm 1, and where the subproblems (2) are defined using the corresponding safe bound $\sigma^{\prime}:=K$. This case delivers the fastest convergence rates in the distributed setting, which in particular don’t degrade as the number of machines $K$ grows and $n$ remains fixed.

Our first main theorem provides convergence guarantees for objectives with non-strongly convex regularizers, including models such as Lasso and sparse logistic regression. Providing primal-dual rates and globally defined primal-dual accuracy certificates requires a theoretical technique that we introduce below, in which we show how to satisfy the following notion of $L$-bounded support.

A function $h$ has $L$-bounded support if its effective domain is bounded by $L$, i.e.,

As we explain in Section F of the appendix, our assumption of $L$-bounded support for the $g_{i}$ functions can be interpreted as an assumption that their conjugates are globally $L$-Lipschitz.

Consider Algorithm 1 with $\gamma:=1$, and let $\Theta$ be the quality of the local solver as in Assumption 1. Let $g_{i}$ have $L$-bounded support, and $f$ be $(1/{\tau})$-smooth. Then after $T$ iterations where

we have that the expected duality gap satisfies

where $\overline{\boldsymbol{\alpha}}$ is the averaged iterate returned by Algorithm 1.

Note that the absolute value function $g_{i}=|\cdot|$ for $L_{1}$ regularization does not have $L$-bounded support, and thus violates the assumptions yielding convergence in Theorem 1. Its dual, the indicator function of the interval, is not defined globally, and thus does not always allow a finite duality gap. To address this, existing approaches typically use a simple smoothing technique as in : by adding a small amount of $L_{2}$ to the $L_{1}$-norm, it becomes strongly convex; see, e.g., . This Nesterov smoothing technique is undesirable in practice, as it changes the iterates, the convergence rate, and the tightness of the resulting duality gap. Further, the amount of smoothing can be difficult to tune and can have a large influence on the performance of the method at hand. We show examples of this issue with experiments in Section 6.

In contrast, our approach preserves all solutions of the original objective, leaves the iterate sequence unchanged, and allows for direct reusability of existing $L_{1}$ solvers. It also removes the need for additional parameter tuning. To achieve this, we modify the function $|\cdot|$ by imposing an additional weak constraint that is inactive in our region of interest. Formally, we replace $g_{i}(\cdot)=|\cdot|$ by

For large enough $B$, this problem yields the same solution as the original $L_{1}$-regularized objective. Note that this only affects convergence theory, in that it allows us to present a strong primal-dual rate (Theorem 1 for $L$=$B$). The modification of $g_{i}$ does not affect the algorithms for the original problems. Whenever a monotone optimizer is used, we will never leave the level set defined by the objective at the starting point. We provide further details on this technique in Section D.3, and illustrate how to leverage it for a variety of applications (see Section C of the appendix and also ).

For the case of strongly convex $g_{i}$, including elastic net-regularized objectives, we obtain the following faster geometric convergence rate.

Consider Algorithm 1 with $\gamma:=1$, and let $\Theta$ be the quality of the local solver as in Assumption 1. Let $g_{i}$ be $\mu$-strongly convex $\forall i\in[n]$, and $f$ be $(1/{\tau})$-smooth. Then we have that $T$ iterations are sufficient for suboptimality $\epsilon_{\mathcal{D}}$, with

We provide proofs of both Theorem 1 and Theorem 2 in the appendix (Section E).

Related Work

For strongly convex regularizers, current state-of-the-art for empirical loss minimization is randomized coordinate ascent on the dual (SDCA) and its accelerated variants, e.g., . In contrast to primal stochastic gradient descent (SGD) methods, the SDCA family is often preferred as it is free of learning-rate parameters and has faster (geometric) convergence guarantees. Interestingly, a similar trend in coordinate solvers has been observed in recent Lasso literature, but with the roles of primal and dual reversed. For those problems, coordinate descent methods on the primal have become state-of-the-art, as in glmnet and extensions ; see, e.g., the overview in . However, primal-dual convergence rates for unmodified coordinate algorithms have to our knowledge been obtained only for strongly convex regularizers to date .

Coordinate descent on $L_{1}$-regularized problems (A) with $g(\cdot)=\lambda\|\cdot\|_{1}$ can be interpreted as the iterative minimization of a quadratic approximation of the smooth part of the objective (as in a one-dimensional Newton step), followed by a shrinkage step resulting from the $L_{1}$ part. In the single-coordinate update case, this is at the core of glmnet , and widely used in, e.g., solvers based on the primal formulation of $L_{1}$-regularized objectives . When changing more than one coordinate at a time, again employing a quadratic upper bound on the smooth part, this results in a two-loop method as in glmnet for the special case of logistic regression. This idea is crucial for the distributed setting.

Parallel coordinate descent for $L_{1}$-regularized objectives (with and without using mini-batches) was proposed in (Shotgun) and generalized in , and is among the best performing solvers in the parallel setting. Our framework reduces to Shotgun as a special case when the internal solver is a single coordinate update on the subproblem (2), $\gamma=1$, and for a suitable $\sigma^{\prime}$. However, Shotgun is not covered by our convergence theory, since it uses a potentially un-safe upper bound $\beta$ instead of $\sigma^{\prime}$, which isn’t guaranteed to satisfy the condition (21). Other parallel coordinate descent methods on the $L_{1}$-objective have recently been analyzed in , but not in the communication-efficient or distributed setting.

The methods most closely related to our approach are distributed variants of glmnet as in . Inspired by glmnet and , the work of introduced the idea of a block-diagonal Hessian upper approximation in the distributed $L_{1}$ context. The later work of specialized this approach to sparse logistic regression.

If hypothetically each of our quadratic subproblems $\mathcal{G}^{\sigma^{\prime}}_{k}(\Delta{\boldsymbol{\alpha}}_{[k]})$ as defined in (2) were to be minimized exactly, the resulting steps could be interpreted as block-wise Newton-type steps on each coordinate block $k$, where the Newton-subproblem is modified to also contain the $L_{1}$-regularizer . While allows a fixed accuracy for these subproblems—but not arbitrary approximation quality $\Theta$ as in our framework—the work of assumes that the quadratic subproblems are solved exactly. Therefore, these methods are not able to freely trade off communication and computation. Also, they do not allow the re-use of arbitrary local solvers. On the theoretical side, the rate results provided by are not explicit convergence rates but only asymptotic, as the quadratic upper bounds are not explicitly controlled for safety as with our $\sigma^{\prime}$.

ADMM , proximal gradient descent, and quasi-Newton methods such as L-BFGS and are also often used in distributed environments because of their relatively low communication requirements. However, they require at least a full (distributed) batch gradient computation at each round, and therefore do not allow the gradual trade-off between communication and computation provided by proxCoCoA$\!{}^{\bf\textbf{+}}$. The works of and have obtained encouraging results for distributed systems employing coordinate descent variants on $L_{1}$-problems. The latter approach distributes both columns and rows of the data matrix and can be extended to Lasso. However it only provides asymptotic improvement per step, and no convergence rate. We include experimental comparisons with ADMM, prox-GD, and orthant-wise limited memory quasi-Newton (OWL-QN) , an L-BFGS variant that can handle $L_{1}$ regularization , but which has no convergence rate.

Finally, we note that while the provided convergence rates for proxCoCoA$\!{}^{\bf\textbf{+}}$ mirror the convergence class of classical batch gradient methods in terms of the number of outer rounds, existing batch proximal gradient methods come with a weaker theory, as they do not allow general inexactness $\Theta$ for the local subproblem (2). In contrast, our shown convergence rates incorporate this approximation directly, and, moreover, hold for arbitrary local solvers of much cheaper cost than batch methods (where in each round, every machine has to process exactly a full pass through the local data). This makes proxCoCoA$\!{}^{\bf\textbf{+}}$ more flexible in the distributed setting, as it can adapt to varied communication costs on real systems. We will see in the following section that this flexibility results in significant performance gains over the competing methods.

Experimental Results

In this section we compare proxCoCoA$\!{}^{\bf\textbf{+}}$ to numerous state-of-the-art methods for large-scale $L_{1}$-regularized optimization, including:

descent with an $L_{1}$-prox • Prox-GD: full proximal gradient descent • OWL-QN: orthant-wise limited quasi-Newton

of multipliers • Mb-CD: mini-batch parallel coordinate descent, incl. Shotgun

The first three methods are optimized and implemented in Apache Spark’s MLlib (v1.5.0) . We employ coordinate descent as a local solver for proxCoCoA$\!{}^{\bf\textbf{+}}$, and apply proxCoCoA$\!{}^{\bf\textbf{+}}$ directly to the primal formulation of Lasso and elastic net, thereby mapping the problem to (A) and solving this objective directly. A comparison with Shotgun is provided as an extreme case to highlight the detrimental effects of frequent communication in the distributed environment.

We test the performance of each method in large-scale experiments fitting Lasso and elastic net regression models to the datasets shown in Table 1. All code is written in Apache Spark and experiments are run on public cloud Amazon EC2 m3.xlarge machines with one core per machine. For Mb-CD, Shotgun, and proxCoCoA$\!{}^{\bf\textbf{+}}$ in the primal, datasets are distributed by feature, whereas for Mb-SGD, Prox-GD, OWL-QN, ADMM, and CoCoA$\!{}^{\bf\textbf{+}}$ they are distributed by datapoint.

We carefully tune each competing method for best performance. ADMM requires the most tuning, both in selecting the penalty parameter $\rho$ and in solving the subproblems. Solving the subproblems to completion for ADMM is prohibitively slow, and we thus use iterations of conjugate gradient and improve performance by allowing early stopping. We also use a varying penalty parameter $\rho$ — practices described in [5, Sec. 4.3, 3.4.1]. For Mb-SGD, we tune the step size and mini-batch size parameters. For Mb-CD, we scale the updates at each round by $\frac{\beta}{b}$ for mini-batch size $b$ and $\beta\in[1,b]$, and tune both parameters $b$ and $\beta$. Further implementation details for all methods are given in the appendix (Section G).

In analyzing the performance of each algorithm (Figure 1), we measure the improvement to the primal objective given in (A) $({\mathcal{D}}({\boldsymbol{\alpha}}))$ in terms of wall-clock time in seconds. We see that both Mb-SGD and Mb-CD are slow to converge, and come with the additional burden of having to tune extra parameters (though Mb-CD makes clear improvements over Mb-SGD). As expected, naively distributing Shotgun (single coordinate updates per machine) does not perform well, as it is tailored to shared-memory systems and requires communicating too frequently. OWL-QN performs the best of all compared methods, but is still much slower to converge than proxCoCoA$\!{}^{\bf\textbf{+}}$, converging, e.g., 50$\times$ more slowly for the webspam dataset. The optimal performance of proxCoCoA$\!{}^{\bf\textbf{+}}$ is particularly evident in datasets with large numbers of features (e.g., url, kddb, webspam), which are exactly the datasets of interest for $L_{1}$ regularization.

Results are shown for regularization parameters $\lambda$ such that the resulting weight vector ${\boldsymbol{\alpha}}$ is sparse. However, our results are robust to varying values of $\lambda$ as well as to various problem settings, as we illustrate in Figure 2.

We note that in contrast to the compared methods, proxCoCoA$\!{}^{\bf\textbf{+}}$ comes with the benefit of having only a single parameter to tune: the subproblem approximation quality, $\Theta$, which can be controlled via the number of local subproblem iterations, $H$. We further explore the effect of this parameter in Figure 3, and provide a general guideline for choosing it in practice (see Remark 1). In particular, we see that while increasing $H$ always results in better performance in terms of rounds, smaller or larger values of $H$ may result in better performance in terms of wall-clock time, depending on the cost of communication and computation. The flexibility to tune $H$ is one of the reasons for proxCoCoA$\!{}^{\bf\textbf{+}}$’s significant performance gains.

Finally, we point out several important ways in which proxCoCoA$\!{}^{\bf\textbf{+}}$ improves upon the CoCoA$\!{}^{\bf\textbf{+}}$ framework . First, CoCoA$\!{}^{\bf\textbf{+}}$ cannot be included in the set of experiments in Figure 1 because it cannot be directly applied to the Lasso objective (CoCoA$\!{}^{\bf\textbf{+}}$ only allows for strongly convex regularizersCoCoA$\!{}^{\bf\textbf{+}}$ in is in fact limited to the case where the regularizer is equal to the $L_{2}$ norm $\frac{1}{2}\|\cdot\|_{2}^{2}$, though the extension to strongly convex regularizers is covered as a special case in our analysis.). Second, as shown in Figure 4, the performance of CoCoA$\!{}^{\bf\textbf{+}}$ degrades drastically when considering datasets with large numbers of features, such as the webspam dataset. One reason for this is that CoCoA$\!{}^{\bf\textbf{+}}$ distributes data by data point, which necessitates communicating a vector of length equal to the feature size. When the feature size is large, this can become expensive. The results shown hold despite the fact that we have tuned $H$ (the number of local solver iterations) separately for both proxCoCoA$\!{}^{\bf\textbf{+}}$ and CoCoA$\!{}^{\bf\textbf{+}}$.

Beyond communication, we also see that CoCoA$\!{}^{\bf\textbf{+}}$ is slower to converge as the regularizer becomes less strongly convex (Figure 4a). Indeed, even when the number of features is relatively low such as for the epsilon dataset, we see that the performance of CoCoA$\!{}^{\bf\textbf{+}}$ degrades significantly as the regularizer approaches pure $L_{1}$. In Figure 4, we illustrate this by implementing the Nesterov smoothing technique used in, e.g., — adding a small amount of strong convexity $\delta\|{\boldsymbol{\alpha}}\|_{2}^{2}$ to the objective for Lasso regression. We show results for decreasing levels of $\delta$. As $\delta$ decreases, the final sparsity of the problem starts to match that of running pure $L_{1}$ (Figure 4c), but the performance also degrades (Figure 4b). We note again that through the modification presented in Section 4, we can deliver strong rates without having to make these fundamental alterations to the problem of interest.

Acknowledgments

We thank Michael P. Friedlander and Martin Takáč for fruitful discussions.

References

Appendix

And analogously if the same holds for all subgradients, in the case of a general closed convex function $f$.

Appendix B Convex Conjugates

Below we list several useful properties of conjugates (see, e.g., [6, Section 3.3.2]):

Double conjugate: $(f^{*})^{*}=f$ if $f$ is closed and convex.

Value Scaling: (for $\alpha>0$) $f({\bf v})=\alpha g({\bf v})\qquad\Rightarrow\qquad f^{*}({\bf w})=\alpha g^{*}({\bf w}/\alpha)\,.$

Argument Scaling: (for $\alpha\neq 0$) $f({\bf v})=g(\alpha{\bf v})\qquad\Rightarrow\qquad f^{*}({\bf w})=g^{*}({\bf w}/\alpha)\,.$

Conjugate of a separable sum: $f({\bf v})=\sum_{i}\phi_{i}(v_{i})\qquad\Rightarrow\qquad f^{*}({\bf w})=\sum_{i}\phi_{i}^{*}(w_{i})\,.$

Given a proper convex function $f$, it holds that $f$ is $L$-Lipschitz if and only if $f^{*}$ has $L$-bounded support.

Given a closed convex function $f$, it holds that $f$ is $\mu$ strongly convex w.r.t. the norm $\|\cdot\|$ if and only if $f^{*}$ is $(1/{\mu})$-smooth w.r.t. the dual norm $\|\cdot\|_{*}$.

Appendix C Applications

$L_{1}$ regularization is obtained in the objective (A) by letting $g_{i}(\cdot):=\lambda|\cdot|$. Primal-dual convergence can be obtained by using the modification introduced in Section 4, which will guarantee $L$-bounded support. Formally, we replace $g_{i}(\cdot)=|\cdot|$ by

For large enough $B$, this problem yields the same solution as the original $L_{1}$-objective. We provide a detailed proof and description of this technique in Section D.3. Note that this only affects convergence theory, in that it allows us to present a strong primal-dual rate (Theorem 1 for $L$=$B$).

C.2 Elastic Net and General Strongly Convex Regularizers

Another application we can consider is elastic net regularization, $\frac{\eta}{2}\left\lVert{{\boldsymbol{\alpha}}}\right\rVert_{2}^{2}+(1-\eta)\left\lVert{{\boldsymbol{\alpha}}}\right\rVert_{1}$, for fixed parameter $\eta\in(0,1]$, which is obtained by setting g_{i}(\alpha):=\lambda\big{[}\frac{\eta}{2}\alpha^{2}+(1-\eta)|\alpha|\big{]} in (A). For the special case $\eta=0$, we obtain the $L_{1}$-norm. For elastic-net-regularized problems of the form (A), Theorem 2 gives a global linear (geometric) convergence rate, since $g_{i}$ is $\eta$-strongly convex. This holds as long as the data-fit function is smooth (see Section C.4), and directly yields a primal-dual algorithm and corresponding rate.

For the $L_{1}$-regularizer in the primal setting, the local subproblem (2) becomes a simple quadratic problem on the local data, with regularization applied only to local variables ${\boldsymbol{\alpha}}_{[k]}$. Therefore, existing fast $L_{1}$-solvers for the single-machine case, such as glmnet variants or blitz can be directly applied to each local subproblem $\mathcal{G}^{\sigma^{\prime}}_{k}(\,\cdot\,;{\bf v},{\boldsymbol{\alpha}}_{[k]})$ within Algorithm 1. The sparsity induced on the subproblem solutions of each machine naturally translates into the sparsity of the global solution, since the local variables ${\boldsymbol{\alpha}}_{[k]}$ will be concatenated.

In terms of the approximation quality parameter $\Theta$ for the local problems (Assumption 1), we can apply existing recent convergence results from the single machine case. For example, for randomized coordinate descent (as part of glmnet), [16, Theorem 1] gives a $O(1/t)$ approximation quality for any separable regularizer, including $L_{1}$ and elastic net; see also .

C.4 Smooth Data-Fit Functions

To illustrate the role of $f$ as a smooth data-fit function in this section—contrasting with its role as a regularizer in traditional CoCoA$\!{}^{\bf\textbf{+}}$ as we discuss in Section F—we consider the following examples.

Observing that the gradient of $f$ is $\nabla f({\bf v})={\bf v}-{\bf b}$, the dual-to-primal mapping is given by: ${\bf w}({\boldsymbol{\alpha}})$ $:=$ $\nabla f({\bf v}({\boldsymbol{\alpha}}))$ $=$ $A{\boldsymbol{\alpha}}-{\bf b}$, which is well known as the residual vector in least-squares regression.

Appendix D Proofs of Primal-Dual Relationship

In the following subsections we provide derivations of the primal-dual relationship of the general objectives (A) and (B), and then show how to derive this primal-dual setup for various applications.

The relation of our original formulation (A) to its dual formulation (B) is standard in convex analysis, and is a special case of the concept of Fenchel Duality. Using the combination with the linear map $A$ as in our case, the relationship is called Fenchel-Rockafellar Duality, see e.g. [4, Theorem 4.4.2] or [2, Proposition 15.18]. For completeness, we illustrate this correspondence with a self-contained derivation of the duality.

The dual problem of (A) follows by taking the infimum with respect to both ${\boldsymbol{\alpha}}$ and ${\bf v}$:

We change signs and turn the maximization of the dual problem (17) into a minimization and thus we arrive at the dual formulation $\eqref{eq:dualP}$ as claimed:

D.2 Conjugates and Smoothness of f𝑓f-Functions of Interest

(see also (15) above) is the conjugate of $f^{*}$, where

with the box constraint $-w_{j}b_{j}\in$.

Furthermore, $f^{*}({\bf w})$ is $1$-strongly convex over its domain if the labels satisfy $b_{j}\in$.

By separability of $f^{*}$, the conjugate of $f^{*}({\bf v})=\sum_{j}\phi^{*}_{j}(v_{j})$ is $f({\bf w})=\sum_{j}\phi_{j}(w_{j})$. For the losses, the conjugate pairs are $\phi_{j}(u)=\log(1+\exp(-b_{j}u))$, and $\phi^{*}_{j}(w_{j})=-w_{j}b_{j}\log(-w_{j}b_{j})+(1+w_{j}b_{j})\log(1+w_{j}b_{j})$ with $-w_{j}b_{j}\in$, see e.g. [26, Page 577].

For the strong convexity, we show 1-strong smoothness of the conjugate $f({\bf v}):=\sum_{j=1}^{d}\log{(1+\exp{(-b_{j}v_{j})})}=\sum_{j=1}^{d}h(b_{j}v_{j})$, which is an equivalent property, see Lemma 4. Using the second derivative $h^{\prime\prime}(a)=\frac{e^{-a}}{(1+e^{-a})^{2}}\leq 1$ of the function $h(a)=\log(1+e^{-a})$, we have that \nabla^{2}f({\bf v})=\mathbf{diag}\big{(}(h^{\prime\prime}(b_{j}v_{j})b_{j}^{2})_{j}\big{)}=\mathbf{diag}\big{(}(\frac{e^{-b_{j}v_{j}}}{(1+e^{-b_{j}v_{j}})^{2}}b_{j}^{2})_{j}\big{)}\leq 1, so $f({\bf v})$ is 1-smooth w.r.t. the Euclidean norm. ∎

D.3 Conjugates of Common Regularizers

For $\eta\in(0,1]$, the elastic net function $g_{i}(\alpha):=\frac{\eta}{2}\alpha^{2}+(1-\eta)|\alpha|$ is the convex conjugate of

where $[.]_{+}$ is the positive part operator, $[s]_{+}=s$ for $s>0$, and zero otherwise. Furthermore, this $g^{*}$ is smooth, i.e. has Lipschitz continuous gradient with constant $1/\eta$.

We start by applying the definition of convex conjugate, that is:

We now distinguish two cases for the optimal: $\alpha^{\star}\geq 0$, $\alpha^{\star}<0$. For the first case we get that

Setting the derivative to we get $\alpha^{\star}=\frac{x-(1-\eta)}{\eta}$. To satisfy $\alpha^{\star}\geq 0$, we must have $x\geq 1-\eta$. Replacing with $\alpha^{\star}$ we thus get:

Similarly we can show that for $x\leq-(1-\eta)$

Finally, by the fact that $g^{*}(.)$ is convex, always positive, and $g^{*}(-(1-\eta))=g^{*}(1-\eta)=0$, it follows that $g^{*}(x)=0$ for every $x\in[-(1-\eta),1-\eta]$.

For the smoothness properties, we consider the derivative of this function $g^{*}(x)$ and see that $g^{*}(x)$ is smooth, i.e. has Lipschitz continuous gradient with constant $1/\eta$, assuming $\eta>0$. ∎

To apply the theoretical convergence result from Theorem 1 to objectives with $L_{1}$ norms, we modify the function $|\cdot|$ by imposing an additional constraint. Consider replacing $g_{i}(\cdot)=|\cdot|$ by

With this modified $L_{1}$-regularizer, the optimization problem (A) with regularization parameter $\lambda$ becomes

For large enough choice of the value $B$, this problems yields the same solution as the original objective:

As we can see, the $\bar{g}$ is nothing more than a constrained version of the absolute value to the interval $[-B,B]$. Therefore by setting $B$ to a large enough value that the interesting values of $\alpha_{i}$ will never reach, we can have continuous $\bar{g}^{*}$ and at the same time make (19) equivalent to (20).

Formally, a simple way to obtain a large enough value of $B$, so that all solutions of (20) are unaffected is the following. Note that we start the algorithm at ${\boldsymbol{\alpha}}={\bf 0}$. For every solution encountered during the execution of the algorithm, the objective values should never become worse than ${\mathcal{D}}({\bf 0})$. In other words, we restrict the ${\mathcal{D}}(\cdot)$ optimization problem to the level set given by the initial starting value. Formally, this means that for every $i$, we will always require:

(Note that $f({\boldsymbol{\alpha}})\geq 0$ holds without loss of generality). We can thus set the value of $B$ to be $\frac{f({\bf 0})}{\lambda}$.

The convex conjugate of $\bar{g}_{i}$ as defined above is

We start by applying the definition of convex conjugate:

We begin by looking at the case in which $\alpha\geq B$; in this case it’s easy to see that when $x\to+\infty$, we have:

as $\alpha-B\geq 0$. The case $\alpha\leq-B$ holds analogously. We’ll now look at the case $\alpha\in[0,B]$; in this case it is clear we must have $x^{\star}\geq 0$. It also must hold that $x^{\star}\leq 1$, since

for every $x>1$. Therefore the maximization becomes

which has maximum $\alpha$ at $x=1$. The remaining $\alpha\in[-B,0]$ case follows in similar fashion. ∎

Appendix E Convergence Proofs

In this section we provide proofs of our main convergence results. The results are motivated by , but where we have significantly generalized the problem of interest, and where we derive separate meaning by applying the problem directly to (A). We provide full details of Lemma 8 as a proof of concept, but omit details in later proofs that can be derived using the arguments in or earlier work of , and instead outline the proof strategy and highlight sections where the theory deviates.

We begin with a definition of the data-dependent aggregation parameter for proxCoCoA$\!{}^{\bf\textbf{+}}$, $\sigma^{\prime}$, which we will use in the throughout our convergence results.

In Algorithm 1, the aggregation parameter $\gamma$ controls the level of adding ($\gamma:=1$) versus averaging ($\gamma:=\tfrac{1}{K}$) of the partial solutions from all machines. For the convergence results discussed below to hold, the subproblem parameter $\sigma^{\prime}$ must be chosen not smaller than

The simple choice of $\sigma^{\prime}:=\gamma K$ is valid for (21), i.e.,

In some cases, it will be possible to give better (data-dependent) choices for $\sigma^{\prime}$, closer to the actual bound given in $\sigma^{\prime}_{min}$.

Our first lemma in the overall proof of convergence helps to relate change in local subproblems to the global objective $\mathcal{D}(\cdot)$.

In this proof we follow the line of reasoning in [17, Lemma 4] with a more general $(1/\tau)$ smoothness assumption on $f(\cdot)$. An outer iteration of proxCoCoA$\!{}^{\bf\textbf{+}}$ performs the following update:

We bound the terms $A$ and $B$ separately. First we bound A using $(1/\tau)$-smoothness of $f$:

Next we use Jensen’s inequality to bound B:

Plugging $A$ and $B$ back into (23) yields:

where the last equality is by the definition of the subproblem objective $\mathcal{G}^{\sigma^{\prime}}_{k}(.)$ as in (2). ∎

E.2 Proof of Main Convergence Result (Theorem 1)

Before proving the main convergence results, we introduce several useful quantities, including the the following lemma, which characterizes the effect of iterations of Algorithm 1 on the duality gap for any chosen local solver of approximation quality $\Theta$.

Let $g_{i}$ be stronglyNote that the case of weakly convex $g_{i}(.)$ is explicitly allowed here as well, as the Lemma holds for the case $\mu=0$. convex with convexity parameter $\mu\geq 0$ with respect to the norm $\|\cdot\|$, $\forall i\in[n]$. Then for all iterations $t$ of Algorithm 1 under Assumption 1, and any $s\in$, it holds that

The line of proof is motivated by [26, Lemma 19] and follows [17, Lemma 5], with a main addition being the extension to our generalized subproblems $\mathcal{G}^{\sigma^{\prime}}_{k}(\cdot;{\bf v},{\boldsymbol{\alpha}}_{[k]})$ along with the general mappings ${\bf w}({\boldsymbol{\alpha}}):=\nabla f({\bf v}({\boldsymbol{\alpha}}))$ with ${\bf v}({\boldsymbol{\alpha}}):=A{\boldsymbol{\alpha}}$.

For simplicity, we write ${\boldsymbol{\alpha}}$ instead of ${\boldsymbol{\alpha}}^{(t)}$, ${\bf v}$ instead of ${\bf v}({\boldsymbol{\alpha}}^{(t)})$, ${\bf w}$ instead of ${\bf w}({\boldsymbol{\alpha}}^{(t)})$ and ${\bf u}$ instead of ${\bf u}^{(t)}$. We can estimate the expected change of the objective $\mathcal{D}({\boldsymbol{\alpha}})$ as follows. Starting from the definition of the update ${\boldsymbol{\alpha}}^{(t+1)}:={\boldsymbol{\alpha}}^{(t)}+\gamma\,\sum_{k}\Delta{\boldsymbol{\alpha}}_{[k]}$ from Algorithm 1, we apply Lemma 8, which relates the local approximation $\mathcal{G}^{\sigma^{\prime}}_{k}({\boldsymbol{\alpha}};{\bf v},{\boldsymbol{\alpha}}_{[k]})$ to the global objective $\mathcal{D}({\boldsymbol{\alpha}})$, and then bound this using the notion of quality of the local solver ($\Theta$), as in Assumption 1. This gives us:

We next upper bound the $C$ term, denoting $\Delta{\boldsymbol{\alpha}}^{\star}=\sum_{k=1}^{K}\Delta{\boldsymbol{\alpha}}^{\star}_{[k]}$. We first plug in the definition of the objective $\mathcal{D}$ in (A) and the local subproblems (2), and then substitute $s(u_{i}-\alpha_{i})$ for $\Delta{\boldsymbol{\alpha}}^{\star}_{i}$ and apply the $\mu$-strong convexity of the $g_{i}$ terms. This gives us:

From the definition of the primal and dual optimization problems (A) and (B), and definition of convex conjugates, we can write the duality gap as:

The convex conjugate maximal property from (26) implies that

The claimed improvement bound (24) then follows by plugging (31) into (27). ∎

The following Lemma provides a uniform bound on $R^{(t)}$:

If $g^{*}_{i}$ are $L$-Lipschitz continuous for all $i\in[n]$, then

[17, Lemma 6]. For general convex functions, the strong convexity parameter is $\mu=0$, and hence the definition (25) of the complexity constant $R^{(t)}$ becomes

[17, Remark 7] If all data points ${\bf x}_{i}$ are normalized such that $\|{\bf x}_{i}\|\leq 1$ $\forall i\in[n]$, then $\sigma_{k}\leq|\mathcal{P}_{k}|=n_{k}$. Furthermore, if we assume that the data partition is balanced, i.e., that $n_{k}=n/K$ for all $k$, then $\sigma\leq n^{2}/K$. This can be used to bound the constants $R^{(t)}$, above, as $R^{(t)}\leq\frac{4L^{2}n^{2}}{K}.$

Consider Algorithm 1, using a local solver of quality $\Theta$ (See Assumption 1). Let $g^{*}_{i}(\cdot)$ be $L$-Lipschitz continuous, and $\epsilon_{G}>0$ be the desired duality gap (and hence an upper-bound on suboptimality $\epsilon_{\mathcal{D}}$). Then after $T$ iterations, where

we have that the expected duality gap satisfies

This proof draws from the line of reasoning in [26, Theorem 2] and follows [17, Theorem 8] but for the more general problem setting (A). We begin by estimating the expected change of feasibility for $\mathcal{D}$. We can bound this above by using Lemma 9 and the fact that the $\mathcal{P}(\cdot)$ is always a lower bound for $-\mathcal{D}(\cdot)$, and then applying (32) to find:

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

Clearly, (38) implies that (39) holds for $t=t_{0}$. Assuming that it holds for any $t\geq t_{0}$, we show that it must also hold for $t+1$. Indeed, using

by applying the bounds (36) and (39), plugging in the definition of $s$ (40), and simplifying. We upperbound the term $D$ using the fact that geometric mean is less or equal to arithmetic mean:

If $\overline{\boldsymbol{\alpha}}$ is defined as (35), we apply the results of Lemma 9 and Lemma 10 to obtain

If $T\geq\lceil\frac{1}{\gamma(1-\Theta)}\rceil+T_{0}$ such that $T_{0}\geq t_{0}$ we have

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

The following main theorem simplifies the results of Theorem 11 and is a generalization of [17, Corollary 9] for general $f^{*}(\cdot)$ functions:

Consider Algorithm 1 with $\gamma:=1$, using a local solver of quality $\Theta$ (see Assumption 1). Let $g^{*}_{i}(\cdot)$ be $L$-Lipschitz continuous, and assume that the columns of $A$ satisfy $\|{\bf x}_{i}\|\leq 1$ $\forall i\in[n]$. Let $\epsilon_{G}>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

(where $\overline{\boldsymbol{\alpha}}$ is the averaged iterate returned by Algorithm 1).

Plug in parameters $\gamma:=1$, $\sigma^{\prime}:=\gamma K=K$ to the results of Theorem 11, and note that for balanced datasets we have $\sigma\leq\frac{n^{2}}{K}$ (see Remark 2). We can further simplify the rate by noting that $\tau=1$ for the 1-smooth losses (least squares and logistic) given as examples in this work. ∎

For pure $L_{1}$-regularized problems as discussed in Section C.1, we have that the above theorem directly delivers a primal-dual convergence with a sublinear rate. This is because in view of Lemma 7, we know that $g^{*}_{i}$ is $B$-Lipschitz for the bounded support modification introduced in Section 4.

Our second main theorem follows reasoning in and is a generalization of [17, Corollary 11]. We first introduce a lemma to simplify the proof.

Assume that $g_{i}$ are $\mu$-strongly convex $\forall i\in[n]$. We define $\sigma_{\max}=\max_{k\in[K]}\sigma_{k}$. Then after $T$ iterations of Algorithm 1, with

Given that $g_{i}(.)$ is $\mu$-strongly convex with respect to the $\|\cdot\|$ norm, we can apply (25) and the definition of $\sigma_{k}$ to find:

where $\sigma_{\max}=\max_{k\in[K]}\sigma_{k}$. If we plug the following value of $s$

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

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

The right hand side will be smaller than some $\epsilon_{\mathcal{D}}$ if

Moreover, to bound the duality gap, we have

Thus, $G({\boldsymbol{\alpha}}^{(t)})\leq\frac{1}{\gamma(1-\Theta)}\frac{\tau\mu+\sigma_{\max}\sigma^{\prime}}{\tau\mu}\epsilon_{\mathcal{D}}^{(t)}$. Hence if $\epsilon_{\mathcal{D}}\leq\gamma(1-\Theta)\frac{\tau\mu}{\tau\mu+\sigma_{\max}\sigma^{\prime}}\epsilon_{G}$ then $G({\boldsymbol{\alpha}}^{(t)})\leq\epsilon_{G}$. Therefore after

iterations we have obtained a duality gap less than $\epsilon_{G}$. ∎

Consider Algorithm 1 with $\gamma:=1$, using a local solver of quality $\Theta$ (See Assumption 1). Let $g_{i}(\cdot)$ be $\mu$-strongly convex $\forall i\in[n]$, and assume that the columns of $A$ satisfy $\|{\bf x}_{i}\|\leq 1$ $\forall i\in[n]$. Then we have that $T$ iterations are sufficient for suboptimality $\epsilon_{\mathcal{D}}$, with

Plug in parameters $\gamma:=1$, $\sigma^{\prime}:=\gamma K=K$ to the results of Theorem 13 and note that for balanced datasets we have $\sigma_{\max}\leq\frac{n}{K}$ (see Remark 2). We can further simplify the rate by noting that $\tau=1$ for the 1-smooth losses (least squares and logistic) given as examples in this work. ∎

For elastic net regularized problems as discussed in Section C.2, we have that the above theorem directly delivers a primal-dual convergence with a geometric rate. This is because in view of Lemma 6, we know that $g^{*}_{i}$ is $1/\eta$-smooth for any elastic net parameter $\eta\in(0,1]$.

Appendix F Recovering CoCoA++\!{}^{\bf\textbf{+}} as a Special Case

As a special case, proxCoCoA$\!{}^{\bf\textbf{+}}$ directly applies to any $L_{2}$-regularized loss-minimization problem, including those presented in . In this setting, the original machine-learning problem is mapped to what we here refer to as the “dual” problem formulation (B):

with $f^{*}(\cdot)=\tfrac{\lambda}{2}\|\cdot\|^{2}$ being the regularizer, and $g^{*}_{i}$ taking the role of loss function, acting on a linear predictor ${\bf x}_{i}^{T}{\bf w}$ (recall that ${\bf x}_{i}$ is a column of the data matrix $A$). In other words, the proxCoCoA$\!{}^{\bf\textbf{+}}$ algorithm will in this case apply to (A) as the dual of the original input problem (which will be mapped to (B)), as described in . The following remarks show that we recover the linear (geometric) convergence rates for smooth loss functions $g^{*}_{i}$, and sublinear convergence for Lipschitz losses. Note that this contrasts the discussed applications of proxCoCoA$\!{}^{\bf\textbf{+}}$ where the $g$ function has the role of the regularizer instead.

This follows since $g^{*}_{i}$ is $L$-Lipschitz if and only if $g_{i}$ has $L$-bounded support [24, Corollary 13.3.3].

This follows since $g^{*}_{i}$ is $\mu$-strongly convex if and only if $g_{i}$ is $(1/{\mu})$-smooth [14, Theorem 6].

Note that the approach of mapping the original objective to (B) does not allow general regularizers such as $L_{1}$. This is one of the reasons we have proposed swapping the roles of regularizers and losses, and running proxCoCoA$\!{}^{\bf\textbf{+}}$ on the primal of the original problem instead.

Appendix G Experiment Details

In this section we provide greater details on the experimental setup and implementations from Section 6. All experiments are run on Amazon EC2 clusters of m3.xlarge machines, with one core per machine. The code for each method is written in Apache Spark, v1.5.0. Our code is open-source and publicly available at: github.com/gingsmith/proxcocoa.

Mini-batch SGD is a standard and widely used method for parallel and distributed optimization. We use the optimized code provided in Spark’s machine learning library, MLlib, v1.5.0. We tune both the size of the mini-batch and the SGD step size using grid search. Proximal gradient descent can be seen as a specific setting of mini-batch SGD, where the mini-batch size is equal to the total number of datapoints. We thus also use the implementation in MLlib for prox-GD, and tune the step size parameter using grid search.

Mini-batch CD aims to improve mini-batch SGD by employing coordinate descent, which has encouraging theoretical and practical backings . We implement mini-batch CD in Spark and scale the updates made at each round by $\frac{\beta}{b}$ for mini-batch size $b$ and $\beta\in[1,b]$, tuning both parameters $b$ and $\beta$ via grid search.

As a special case of mini-batch CD, Shotgun is a popular method for parallel optimization. Shotgun can be seen an extreme case of mini-batch CD where the mini-batch is set to $1$ element per machine, i.e., there is a single update made by each machine per round. We see in the experiments that communicating this frequently becomes prohibitively slow in the distributed environment.

OWN-QN is a quasi-Newton method optimized in Spark’s spark.ml package. Outer iterations of OWL-QN make significant progress towards convergence, but the iterations themselves can be slow because they require processing the entire dataset. proxCoCoA$\!{}^{\bf\textbf{+}}$, the mini-batch methods, and ADMM with early stopping all improve on this by allowing the flexibility of only a subset of the dataset to be processed at each iteration. proxCoCoA$\!{}^{\bf\textbf{+}}$ and ADMM have even greater flexibility by allowing internal methods to process the dataset more than once. proxCoCoA$\!{}^{\bf\textbf{+}}$ makes this approximation quality specific, both in theoretical convergence rates and by providing general guidelines for setting the parameter.

We implement proxCoCoA$\!{}^{\bf\textbf{+}}$ with coordinate descent as a local solver. We note that since the framework and theory allow any internal solver to be used, proxCoCoA$\!{}^{\bf\textbf{+}}$ could benefit even beyond the results shown, by using existing fast $L_{1}$-solvers for the single-machine case, such as glmnet variants or blitz . The only parameter necessary to tune for proxCoCoA$\!{}^{\bf\textbf{+}}$ is the level of approximation quality, which we parameterize in the experiments using $H$, the number of local iterations of the iterative method run locally. Our theory relates local approximation quality to global convergence, and we provide a guideline for how to choose this value in practice that links the value to the systems environment at hand (Remark 1). We implement CoCoA$\!{}^{\bf\textbf{+}}$ as a special case of proxCoCoA$\!{}^{\bf\textbf{+}}$ for elastic net regularized objectives by mapping the main objective to (B) according to the steps described in Section F, and again use coordinate descent as a local solver.