Communication-Efficient Distributed Optimization of Self-Concordant Empirical Loss

Introduction

Many optimization problems in data science (including statistics, machine learning, data mining, etc.) are formulated with a large amount of data as input. They are typically solved by iterative algorithms which need to access the whole dataset or at least part of it during each iteration. With the amount of data we collect and process growing at a fast pace, it happens more often that the dataset involved in an optimization problem cannot fit into the memory or storage of a single computer (machine). To solve such “big data” optimization problems, we need to use distributed algorithms that rely on inter-machine communication.

In this paper, we focus on distributed optimization problems generated through sample average approximation (SAA) of stochastic optimization problems. Consider the problem

Our goal is to minimize the overall empirical loss defined with all $mn$ samples:

For stability and generalization purposes, we often add a regularization term $(\lambda/2)\|w\|_{2}^{2}$ to make the empirical loss function strongly convex. More specifically, we modify the definition of $f_{i}(w)$ as

For example, when $\phi$ is the hinge loss, this formulation yields the support-vector machine .

Since the functions $f_{i}(w)$ can be accessed only locally, we consider distributed algorithms that alternate between a local computation procedure at each machine, and a communication round involving simple map-reduce type of operations . Compared with local computation at each machine, the cost of inter-machine communication is much higher in terms of both speed/delay and energy consumption (e.g., ); thus it is often considered as the bottleneck for distributed computing. Our goal is to develop communication-efficient distributed algorithms, which try to use a minimal number of communication rounds to reach certain precision in minimizing $f(w)$.

For a concrete discussion, we make the following assumption:

where $f^{\prime\prime}(w)$ denotes the Hessian of $f$ at $w$, and $I$ is the $d\times d$ identity matrix.

Functions that satisfy Assumption A are often called $L$-smooth and $\lambda$-strongly convex. The value $\kappa=L/\lambda\geq 1$ is called the condition number of $f$, which is a key quantity in characterizing the complexity of iterative algorithms. We focus on ill-conditioned cases where $\kappa\gg 1$.

Another popular technique for distributed optimization is to use the alternating direction method of multipliers (ADMM); see, e.g., [8, Section 8]. Under the assumption that each local function $f_{i}$ is $L$-smooth and $\lambda$-strongly convex, the ADMM approach can achieve linear convergence, and the best known complexity is $\mathcal{O}(\sqrt{\kappa}\log(1/\epsilon))$ . This turns out to be the same order as for accelerated gradient methods. In this case, ADMM can actually be considered as an accelerated primal-dual first-order method; see the discussions in [10, Section 4].

The polynomial dependence of the iteration complexity on the condition number can be unsatifactory. For machine learning applications, both the precision $\epsilon$ and the regularization parameter $\lambda$ should decrease while the overall sample size $mn$ increases, typically on the order of $\Theta(1/\sqrt{mn})$ (e.g., ). This translates into the condition number $\kappa$ being $\Theta(\sqrt{mn})$. In this case, the iteration complexity, and thus the number of communication rounds, scales as $(mn)^{1/4}$ for both accelerated gradient methods and ADMM (with careful tuning of the penalty parameter). This suggests that the number of communication rounds grows with the total sample size.

Despite the rich literature on distributed optimization (e.g., ), most algorithms involve high communication cost. In particular, their iteration complexity have similar or worse dependency on the condition number as the methods discussed above. It can be argued that the iteration complexity $\mathcal{O}(\sqrt{\kappa}\log(1/\epsilon))$ cannot be improved in general for distributed first-order methods — after all, it is optimal for centralized first-order methods under the same assumption that $f(w)$ is $L$-smooth and $\lambda$-strongly convex . Thus in order to obtain better communication efficiency, we need to look into further problem structure and/or alternative optimization methods. And we need both in this paper.

First, we note that the above discussion on iteration complexity does not exploit the fact that each function $f_{i}$ is generated by, or can be considered as, SAA of a stochastic optimization problem. Since the data $z_{i,j}$ are i.i.d. samples from a common distribution, the local empirical loss functions $f_{i}(w)=(1/n)\sum_{j=1}^{n}\phi(w,z_{i,j})$ will be similar to each other if the local sample size $n$ is large. Under this assumption, Zhang et al. studied a one-shot averaging scheme that approximates the minimizer of function $f$ by simply averaging the minimizers of $f_{i}$. For a fixed condition number, the one-shot approach is communication efficient because it achieves optimal dependence on the overall sample size $mn$ (in the sense of statistical lower bounds). But their conclusion doesn’t allow the regularization parameter $\lambda$ to decrease to zero as $n$ goes to infinity (see discussions in ).

Exploiting the stochastic nature alone seems not enough to overcome ill-conditioning in the regime of first-order methods. This motivates the development of distributed second-order methods. Recently, Shamir et al. proposed a distributed approximate Newton-type (DANE) method. Their method takes advantage of the fact that, under the stochastic assumptions of SAA, the Hessians $f^{\prime\prime}_{1},f^{\prime\prime}_{2},\dots,f^{\prime\prime}_{m}$ are similar to each other. For quadratic loss functions, DANE is shown to converge in \widetilde{\mathcal{O}}\bigl{(}(L/\lambda)^{2}n^{-1}\log(1/\epsilon)\bigr{)} iterations with high probability, where the notation $\widetilde{\mathcal{O}}(\cdot)$ hides additional logarithmic factors involving $m$ and $d$. If $\lambda\sim 1/\sqrt{mn}$ as in machine learning applications, then the iteration complexity becomes $\widetilde{\mathcal{O}}(m\log(1/\epsilon))$, which scales linearly with the number of machines $m$, not the total sample size $mn$. However, the analysis in does not guarantee that DANE has the same convergence rate on non-quadratic functions.

2 Outline of our approach

where $\Delta w_{k}$ and $\delta(w_{k})$ are the Newton step and the Newton decrement, respectively, defined as

Since $f$ is the average of $f_{1},\ldots,f_{m}$, its gradient and Hessian can be written as

In order to compute $\Delta w_{k}$ in a distributed setting, the naive approach would require all the machines to send their gradients and Hessians to a master node (say machine 1). However, the task of transmitting the Hessians (which are $d\times d$ matrices) can be prohibitive for large dimensions $d$. A better alternative is to use the conjugate gradient (CG) method to compute $\Delta w_{k}$ as the solution to a linear system $f^{\prime\prime}(w_{k})\Delta w_{k}=f^{\prime}(w_{k})$. Each iteration of the CG method requires a matrix-vector product of the form

The overall method has two levels of loops: the outer-loop of the damped Newton method, and the inner loop of the CG method for computing the inexact Newton steps. A similar approach (using a distributed truncated Newton method) was proposed in for ERM of linear predictors, and it was reported to perform very well in practice. However, the total number of CG iterations (each takes a round of communication) may still be high.

First, consider the outer loop complexity. It is well-known that Newton-type methods have asymptotic superlinear convergence. However, in classical analysis of Newton’s method (e.g., [9, Section 9.5.3]), the number of steps needed to reach the superlinear convergence zone still depends on the condition number; more specifically, it scales quadratically in $\kappa$. To solve this problem, we resort to the machinery of self-concordant functions . For self-concordant empirical losses, we show that the iteration complexity of the inexact damped Newton method has a much weaker dependence on the condition number.

Second, consider the inner loop complexity. The convergence rate of the CG method also depends on the condition number $\kappa$: it takes $\mathcal{O}(\sqrt{\kappa}\log(1/\varepsilon))$ CG iterations to compute an $\varepsilon$-precise Newton step. Thus we arrive at the dilemma that the overall complexity of the CG-powered inexact Newton method is no better than accelerated gradient methods or ADMM. To overcome this difficulty, we exploit the stochastic nature of the problem and propose to use a preconditioned CG (PCG) method for solving the Newton system. Roughly speaking, if the local Hessians $f^{\prime\prime}_{1}(w_{k}),\ldots,f^{\prime\prime}_{m}(w_{k})$ are “similar” to each other, then we can use any local Hessian $f^{\prime\prime}_{i}(w_{k})$ as a preconditioner. Without loss of generality, let $P=f^{\prime\prime}_{1}(w_{k})+\mu I$, where $\mu$ is an estimate of the spectral norm $\|f^{\prime\prime}_{1}(w_{k})-f^{\prime\prime}(w_{k})\|_{2}$. Then we use CG to solve the pre-conditioned linear system

where the preconditioning (multiplication by $P^{-1}$) can be computed locally at machine $1$ (the master node). The convergence rate of PCG depends on the condition number of the matrix $P^{-1}f^{\prime\prime}(w_{k})$, which is close to 1 if the spectral norm $\|f^{\prime\prime}_{1}(w_{k})-f^{\prime\prime}(w_{k})\|_{2}$ is small.

To exactly characterize the similarity between $f^{\prime\prime}_{1}(w_{k})$ and $f^{\prime\prime}(w_{k})$, we rely on stochastic analysis in the framework of SAA or ERM. We show that with high probability, $\|f^{\prime\prime}_{1}(w_{k})-f^{\prime\prime}(w_{k})\|_{2}$ decreases as $\widetilde{\mathcal{O}}(\sqrt{d/n})$ in general, and $\widetilde{\mathcal{O}}(\sqrt{1/n})$ for quadratic loss. Therefore, when $n$ is large, the preconditioning is very effective and the PCG method converges to sufficient precision within a small number of iterations. The stochastic assumption is also critical for obtaining an initial point $w_{0}$ which further brings down the overall iteration complexity.

Combining the above ideas, we propose and analyze an algorithm for Distributed Self-Concordant Optimization (DiSCO, which also stands for Distributed Second-Order method, or Distributed Stochastic Convex Optimization). We show that several popular empirical loss functions in machine learning, including ridge regression, regularized logistic regression and a (new) smoothed hinge loss, are actually self-concordant. For ERM with these loss functions, Table 1 lists the number of communication rounds required by DiSCO and several other algorithms to find an $\epsilon$-optimal solution. As the table shows, the communication cost of DiSCO weakly depends on the number of machines $m$ and on the feature dimension $d$, and is independent of the local sample size $n$ (excluding logarithmic factors). Comparing to DANE , DiSCO not only improves the communication efficiency on quadratic loss, but also handles non-quadratic classification tasks.

The rest of this paper is organized as follows. In Section 2, we review the definition of self-concordant functions, and show that several popular empirical loss functions used in machine learning are either self-concordant or can be well approximated by self-concordant functions. In Section 3, we analyze the iteration complexity of an inexact damped Newton method for minimizing self-concordant functions. In Section 4, we show how to compute the inexact Newton step using a distributed PCG method, describe the overall DiSCO algorithm, and discuss its communication complexity. In Section 5, we present our main theoretical results based on stochastic analysis, and apply them to linear regression and classification. In Section 6, we report experiment results to illustrate the advantage of DiSCO in communication efficiency, compared with other algorithms listed in Table 1. Finally, we discuss the extension of DiSCO to distributed minimization of composite loss functions in Section 7, and conclude the paper in Section 8.

Self-concordant empirical loss

The reader may refer to the books for detailed treatment of self-concordance. In particular, the following lemma [34, Corollary 4.1.2] states that any self-concordant function can be rescaled to become standard self-concordant.

If a function $f$ is self-concordant with parameter $M_{f}$, then $\frac{M_{f}^{2}}{4}f$ is standard self-concordant (with parameter $2$).

In the rest of this section, we show that several popular regularized empirical loss functions for linear regression and binary classification are either self-concordant or can be well approximated by self-concordant functions.

First we consider regularized linear regression (ridge regression) with

For binary classification, we consider the following regularized empirical loss function

For logistic regression, we minimize the objective function (7) where $\varphi$ is the logistic loss: $\varphi(t)=\log(1+e^{-t})$. We can calculate the second and the third derivatives of $\varphi(t)$:

Smoothed hinge loss

In classification tasks, it is sometimes more favorable to use the hinge loss $\varphi(t)=\max\{0,1-t\}$ than using the logistic loss. We consider a family of smoothed hinge loss functions $\varphi_{p}$ parametrized by a positive number $p\geq 3$. The function is defined by

We plot the functions $\varphi_{p}$ for $p=3,5,10,20$ on Figure 1. As the plot shows, $\varphi_{p}(t)$ is zero for $t>2$, and it is a linear function with unit slope for $t<-\frac{p-3}{p-1}$. These two linear zones are connected by three smooth non-linear segments on the interval $[-\frac{p-3}{p-1},2]$.

The smoothed hinge loss $\varphi_{p}$ satisfies the condition of Lemma 2 with $Q=p-2$ and $\alpha=\frac{1}{p-2}$. To see this, we note that the third derivative of $\varphi_{p}(t)$ is nonzero only when $t\in[-\frac{p-3}{p-1},1-\frac{p-3}{p-1}]$ and when $t\in$. On the first interval, we have

Inexact damped Newton method

The explicit account of approximation errors is essential for distributed optimization. In particular, if $f(w)=(1/m)\sum_{i=1}^{m}f_{i}(w)$ and the components $f_{i}$ locate on separate machines, then we can only perform Newton updates approximately with limited communication budget. Even in a centralized setting on a single machine, analysis of approximation errors can be important if the Newton system is solved by iterative algorithms such as the conjugate gradient method.

Before presenting the convergence analysis, we need to introduce two auxiliary functions

These two functions are very useful for characterizing the properties of self-concordant functions; see [34, Section 4.1.4] for a detailed account. Here, we simply note that $\omega(0)=\omega_{*}(0)=0$, both are strictly increasing for $t\geq 0$, and $\omega_{*}(t)\to\infty$ as $t\to 1$.

We also need to define two auxiliary vectors

The norm of the first vector, $\|\widetilde{u}_{k}\|_{2}=\sqrt{f^{\prime}(w_{k})^{T}[f^{\prime\prime}(w_{k})]^{-1}f^{\prime}(w_{k})}$, is the exact Newton decrement. The norm of the second one is $\|{\widetilde{v}_{k}}\|_{2}=\delta_{k}$, which is computed during each iteration of Algorithm 1. Note that we do not compute $\widetilde{u}_{k}$ or $\widetilde{v}_{k}$ in Algorithm 1. They are introduced solely for the purpose of convergence analysis. The following Theorem is proved in Appendix A.

For any $k\geq 0$, we have $f(w_{k+1})\leq f(w_{k})-\frac{1}{2}\omega(\|{\widetilde{u}_{k}}\|_{2})$.

If $\|{\widetilde{u}_{k}}\|_{2}\leq 1/6$, then we have $\omega(\|{\widetilde{u}_{k+1}}\|_{2})\leq\frac{1}{2}\omega(\|{\widetilde{u}_{k}}\|_{2})$.

As mentioned before, when $\epsilon_{k}=0$, the vector $v_{k}=[f^{\prime\prime}(w_{k})]^{-1}f^{\prime}(w_{k})$ becomes the exact Newton step. In this case, we have $\widetilde{v}_{k}=\widetilde{u}_{k}$, and it can be shown that $f(w_{k+1})\leq f(w_{k})-\omega(\|{\widetilde{u}_{k}}\|_{2})$ for all $k\geq 0$ and the exact damped Newton method has quadratic convergence when $\|{\widetilde{u}_{k}}\|_{2}$ is small (see [34, Section 4.1.5]). With the approximation error $\epsilon_{k}$ specified in (10), we have

Appendix A shows that when $\beta$ is sufficiently small, the above inequality leads to the conclusion in part (a). Compared with the exact damped Newton method, the guaranteed reduction of the objective value per iteration is cut by half.

Part (b) of Theorem 1 suggests a linear rate of convergence when $\|{\widetilde{u}_{k}}\|_{2}$ is small. This is slower than the quadratic convergence rate of the exact damped Newton method, due to the allowed approximation errors in computing the Newton step. Actually, superlinear convergence can be established if we set the tolerances $\epsilon_{k}$ to be small enough; see Appendix B for detailed analysis. However, when $v_{k}$ is computed through a distributed iterative algorithm (like the distributed PCG algorithm in Section 4.2), a smaller $\epsilon_{k}$ would require more local computational effort and more rounds of inter-machine communication. The choice in equation (10) is a reasonable trade-off in practice.

Using Theorem 1, we can derive the iteration complexity of Algorithm 1 for obtaining an arbitrary accuracy. We present this result as a corollary.

Here $\lceil t\rceil$ denotes the smallest nonnegative integer that is larger or equal to $t$.

Since $\omega(t)$ is strictly increasing for $t\geq 0$, part (a) of Theorem 1 implies that if $\|\widetilde{u}_{k}\|_{2}>1/6$, one step of Algorithm 1 decreases the value of $f(w)$ by at least a constant $\frac{1}{2}\omega(1/6)$. So within at most $K_{1}=\lceil\frac{f(w_{0})-f(w_{\star})}{\frac{1}{2}\omega(1/6)}\rceil$ iterations, we are guaranteed that $\|{\widetilde{u}_{k}}\|_{2}\leq 1/6$.

According to [34, Theorem 4.1.13], if $\|{\widetilde{u}_{k}}\|_{2}<1$, then we have

Moreover, it is easy to check that $\omega_{*}(t)\leq 2\,\omega(t)$ for $0\leq t\leq 1/6$. Therefore, using part (b) of Theorem 1, we conclude that when $k\geq K_{1}$,

Bounding the right-hand side of the above inequality by $\epsilon$, we have $f(w_{k})-f(w_{\star})\leq\epsilon$ whenever $k\geq K_{1}+\left\lceil\log_{2}\left(\frac{2\,\omega(1/6)}{\epsilon}\right)\right\rceil=K$, which is the desired result.

We note that when $\|{\widetilde{u}_{k}}\|_{2}\leq 1/6$ (as long as $k\geq K_{1}$), we have $f(w_{k})-f(w_{\star})\leq 2\,\omega(1/6)$. Thus for $\epsilon>2\,\omega(1/6)$, it suffices to have $k\geq K_{1}$. ∎

We discuss two stopping criteria for Algorithm 1. The first one is based on the strong convexity of $f$, which leads to the inequality (e.g., [34, Theorem 2.1.10])

Therefore, we can use the stopping criterion $\|f^{\prime}(w_{k})\|_{2}\leq\sqrt{2\lambda\epsilon}$, which implies $f(w_{k})-f(w_{\star})\leq\epsilon$. However, this choice can be too conservative in practice (see discussions in [9, Section 9.1.2]).

Another choice for the stopping criterion is based on self-concordance. Using the fact that $\omega_{*}(t)\leq t^{2}$ for $0\leq t\leq 0.68$ (see [9, Section 9.6.3]), we have

provided $\|{\widetilde{u}_{k}}\|_{2}\leq 0.68$. Since we do not compute $\|{\widetilde{u}_{k}}\|_{2}$ (the exact Newton decrement) directly in Algorithm 1, we can use $\delta_{k}$ as an approximation. Using the inequality (11), and noticing that $\|{\widetilde{v}_{k}}\|_{2}=\delta_{k}$, we conclude that

implies $f(w_{k})-f(w_{\star})\leq\epsilon$ when $\epsilon\leq 0.68^{2}$. Since $\delta_{k}$ is computed at each iteration of Algorithm 1, this can serve as a good stopping criterion.

2 Scaling for non-standard self-concordant functions

Using the sequence $\{\epsilon_{k}\}$ defined in (10), the condition for computing $v_{k}$ in Step 1 is

In other words, the precision requirement in Step 1 is scaling invariant.

Step 2 of Algorithm 1 can be rewritten as

The DiSCO algorithm

In this section, we adapt the inexact damped Newton method (Algorithm 1) to a distributed system, in order to minimize

where each function $f_{i}$ can only be evaluated locally at machine $i$ (see background in Section 1). This involves two questions: (1) how to set the initial point $w_{0}$ and (2) how to compute the inexact Newton step $v_{k}$ in a distributed manner. After answering these two questions, we will present the overall DiSCO algorithm and analyze its communication complexity.

In accordance with the averaging structure in (18), we choose the initial point based on averaging. More specifically, we let

where each $\widehat{w}_{i}$ is the solution to a local optimization problem at machine $i$:

Here we comment on the computational cost of solving (20) locally at each machine. Suppose each $f_{i}(w)$ has the form in (3), then the local optimization problems in (20) become

The finite average structure of the above objective function can be effectively exploited by the stochastic average gradient (SAG) method or its new variant SAGA . Each step of these methods processes only one component function $\phi(w,z_{i,j})$, picked uniformly at random. Suppose $f_{i}(w)$ is $L$-smooth, then SAG returns an $\epsilon$-optimal solution with \mathcal{O}\bigl{(}(n+\frac{L+\rho}{\lambda+\rho})\log(1/\epsilon)\bigr{)} steps of stochastic updates. For ERM of linear predictors, we can also use the stochastic dual coordinate ascent (SDCA) method , which has the same complexity. We also mention some recent progress in accelerated stochastic coordinate gradient methods , which can be more efficient both in theory and practice.

2 Distributed computing of the inexact Newton step

In each iteration of Algorithm 1, we need to compute an inexact Newton step $v_{k}$ such that $\|f^{\prime\prime}(w_{k})v_{k}-f^{\prime}(w_{k})\|_{2}\leq\epsilon_{k}$. This boils down to solving the Newton system $f^{\prime\prime}(w_{k})v_{k}=f^{\prime}(w_{k})$ approximately. When the objective $f$ has the averaging form (18), its Hessian and gradient are given in (6). In the setting of distributed optimization, we propose to use a preconditioned conjugate gradient (PCG) method to solve the Newton system.

To simplify notation, we use $H$ to represent $f^{\prime\prime}(w_{k})$ and use $H_{i}$ to represent $f^{\prime\prime}_{i}(w_{k})$. Without loss of generality, we define a preconditioning matrix using the local Hessian at the first machine (the master node):

where $\mu>0$ is a small regularization parameter. Algorithm 2 describes our distributed PCG method for solving the preconditioned linear system

In particular, the master machine carries out the main steps of the classical PCG algorithm (e.g., [20, Section 10.3]), and all machines (including the master) compute the local gradients and Hessians and perform matrix-vector multiplications. Communication between the master and other machines are used to form the overall gradient $f^{\prime}(w_{k})$ and the matrix-vector products

We note that the overall Hessian $H=f^{\prime\prime}(w_{k})$ is never formed and the master machine only stores and updates the vectors $Hu^{(t)}$ and $Hv^{(t)}$.

As explained in Section 1.2, the motivation for preconditioning is that when $H_{1}$ is sufficiently close to $H$, the condition number of $P^{-1}H$ might be close to $1$, which is much smaller than that of $H$ itself. As a result, the PCG method may converge much faster than CG without preconditioning. The following lemma characterizes the extreme eigenvalues of $P^{-1}H$ based on the closeness between $H_{1}$ and $H$.

Suppose Assumption A holds. If $\|{H_{1}-H}\|_{2}\leq\mu$, then we have

Since both $P$ and $H$ are symmetric and positive definite, all eigenvalues of $P^{-1}H$ are positive real numbers (e.g., [21, Section 7.6]). The eigenvalues of $P^{-1}H$ are identical to that of $P^{-1/2}HP^{-1/2}$. Thus, it suffices to prove inequalities (22) and (23) for the matrix $P^{-1/2}HP^{-1/2}$. To prove inequality (22), we need to show that $H\preceq P=H_{1}+\mu I$. This is equivalent to $H-H_{1}\preceq\mu I$, which is a direct consequence of the assumption $\|H_{1}-H\|_{2}\leq\mu I$.

Similarly, the second inequality (23) is equivalent to $H\succeq\frac{\lambda}{\lambda+2\mu}(H_{1}+\mu I)$, which is the same as $\frac{2\mu}{\lambda}H-\mu I\succeq H_{1}-H$. Since $H\succeq\lambda I$ (by Assumption A), we have $\frac{2\mu}{\lambda}H-\mu I\succeq\mu I$. The additional assumption $\|H_{1}-H\|_{2}\leq\mu I$ implies $\mu I\succeq H_{1}-H$, which complete the proof. ∎

By Assumption A, the condition number of the Hessian matrix is $\kappa(H)=L/\lambda$, which can be very large if $\lambda$ is small. Lemma 3 establishes that the condition number of the preconditioned linear system is

provided that $\|H_{1}-H\|_{2}\leq\mu$. When $\mu$ is small (comparable with $\lambda$), the condition number $\kappa(P^{-1}H)$ is close to one and can be much smaller than $\kappa(H)$. Based on classical convergence analysis of the CG method (e.g., ), the following lemma shows that Algorithm 2 terminates in $\mathcal{O}(\sqrt{1+\mu/\lambda})$ iterations. See Appendix C for the proof.

Suppose Assumption A holds and assume that $\|{H_{1}-H}\|_{2}\leq\mu$. Let

Then Algorithm 2 terminates in $T_{\mu}$ iterations and the output $v_{k}$ satisfies $\|{Hv_{k}-f^{\prime}(w_{k})}\|_{2}\leq\epsilon_{k}$.

When the tolerance $\epsilon_{k}$ is chosen as in (10), the iteration bound $T_{\mu}$ is independent of $f^{\prime}(w_{k})$, i.e.,

Under Assumption A, we always have $\|{H_{1}-H}\|_{2}\leq L$. If we choose $\mu=L$, then Lemma 4 implies that Algorithm 2 terminates in $\widetilde{\mathcal{O}}(\sqrt{L/\lambda})$ iterations. where the notation $\widetilde{\mathcal{O}}(\cdot)$ hides logarithmic factors. In practice, however, the matrix norm $\|{H_{1}-H}\|_{2}$ is usually much smaller than $L$ due to the stochastic nature of $f_{i}$. Thus, we can choose $\mu$ to be a tight upper bound on $\|{H_{1}-H}\|_{2}$, and expect the algorithm terminating in $\widetilde{\mathcal{O}}(\sqrt{\mu/\lambda})$ iterations. In Section 5, we show that if the local empirical losses $f_{i}$ are constructed with $n$ i.i.d. samples from the same distribution, then $\|{H_{1}-H}\|_{2}\sim 1/\sqrt{n}$ with high probability. As a consequence, the iteration complexity of Algorithm 2 is upper bounded by $\widetilde{\mathcal{O}}(1+\lambda^{-1/2}n^{-1/4})$.

We wrap up this section by discussing the computation and communication complexities of Algorithm 2. The bulk of computation is at the master machine, especially computing the vector $s^{(t)}=P^{-1}r^{(t)}$ in Step 3, which is equivalent to minimize the quadratic function $(1/2)s^{T}Ps-s^{T}r^{(t)}$. Using $P=f^{\prime\prime}_{1}(w_{k})+\mu I$ and the form of $f_{1}(w)$ in (3), this is equivalent to

This problem has the same structure as (21), and an $\epsilon$-optimal solution can be obtained with \mathcal{O}\bigl{(}(n+\frac{L+\mu}{\lambda+\mu})\log(1/\epsilon)\bigr{)} stochastic-gradient type of steps (see discussions at the end of Section 4.1).

As for the communication complexity, we need one round of communication at the beginning of Algorithm 2 to compute $f^{\prime}(w_{k})$. Then, each iteration takes one round of communication to compute $Hu^{(t)}$ and $Hv^{(t)}$. Thus, the total rounds of communication is bounded by $T_{\mu}+1$.

3 Communication efficiency of DiSCO

Putting everything together, we present the DiSCO algorithm in Algorithm 3. Here we study its communication efficiency. Recall that by one round of communication, the master machine broadcasts a message of $\mathcal{O}(d)$ bits to all machines, and every machine processes the aggregated message and sends a message of $\mathcal{O}(d)$ bits back to the master. The following proposition gives an upper bound on the number of communication rounds taken by the DiSCO algorithm.

Assume that $f$ is a standard self-concordant function and it satisfies Assumption A. Suppose the input parameter $\mu$ in Algorithm 3 is an upper bound on $\|{f^{\prime\prime}_{1}(w_{k})-f^{\prime\prime}(w_{k})}\|_{2}$ for all $k\geq 0$. Then for any $\epsilon>0$, in order to find a solution $\widehat{w}$ satisfying $f(\widehat{w})-f(w_{\star})<\epsilon$, the total number of communication rounds $T$ is bounded by

Ignoring logarithmic terms and universal constants, the rounds of communication $T$ is bounded by

First we notice that the number of communication rounds in each call of Algorithm 2 is no more than $1+T_{\mu}$, where $T_{\mu}$ is given in (25), and the extra 1 accounts for the communication round to form $f^{\prime}(w_{k})$. Corollary 1 states that in order to guarantee $f(w_{k})-f(w_{\star})\leq\epsilon$, the total number of calls of Algorithm 2 in DiSCO is bounded by $K$ given in (12). Thus the total number of communication rounds is bounded by $1+K(1+T_{\mu})$, where the extra one count is for computing the initial point $w_{0}$ defined in (19). ∎

It can be hard to give a good a priori estimate of $\mu$ that satisfies the condition in Theorem 2. In practice, we can adjust the value of $\mu$ adaptively while running the algorithm. Inspired by a line search procedure studied in , we propose an adaptive DiSCO method, described in Algorithm 4. The following proposition bounds the rounds of communication required by this algorithm.

Let $n_{k}$ be the number of calls to Algorithm 2 during the $k$th iteration of Algorithm 4. We have

The total number of calls to Algorithm 2 is

From the above proof, we see that the average number of calls to Algorithm 2 at each iteration is $2+\frac{1}{K}\log_{2}\left(\frac{\mu_{K}}{\mu_{0}}\right)$, roughly twice as the non-adaptive Algorithm 3. Ignoring logarithmic terms and universal constants, the number of communication round $T$ used by Algorithm 4 is bounded by

In general, we can update $\mu_{k}$ in Algorithm 4 as follows:

4 A simple variant without PCG iterations

We consider a simple variant of DiSCO where the approximate Newton step $v_{k}$ is computed without using the PCG method described in Algorithm 2. Instead, we simply set

which is equivalent to setting $v_{k}=s^{(0)}$ in the initialization phase of Algorithm 2, or forcing it to always exit during the first PCG iteration. (The latter choice gives the same search direction but with a slightly different scaling.) In this variant, each iteration of the inexact damped Newton method requires two communication rounds: one to form $f^{\prime}(w_{k})$ and another to compute the stepsize parameter $\delta_{k}=(v_{k}^{T}f^{\prime\prime}(w_{k})v_{k})^{1/2}$.

A distributed algorithm that is similar to this variant of DiSCO is proposed in . It does not compute $\delta_{k}$; instead it uses line search to determine the step size, which also requires extra round(s) of communication. It is shown in that this method works well in experiments, requiring less number of iterations to converge than ADMM. However, according to their theoretical analysis, its iteration complexity still depends badly on the condition number.

Here we examine the theoretical conditions under which this variant of DiSCO enjoys a low iteration complexity. Recall the two auxiliary vectors defined in Section 3:

The norm of their difference can be bounded as

From Lemma 3, we know that when $\|{H_{1}-H}\|_{2}\leq\mu$, the eigenvalues of $P^{-1}H$ are located within the interval $[\frac{\lambda}{\lambda+2\mu},1]$. Therefore, we have

This inequality has the same form as (11), which is responsible to obtain the desired low complexity result if $\frac{2\mu}{\lambda+\mu}$ is sufficiently small. Indeed, if $\frac{2\mu}{\lambda+2\mu}\leq\beta=\frac{1}{20}$ as specified in (10), the same convergence rate and complexity result stated in Theorem 1 and Corollary 1 apply. Since each iteration of the damped Newton method involves only two communication rounds (to compute $f^{\prime}(w_{k})$ and $\delta_{k}$ respectively), we have the following corollary.

Assume that $f$ is a standard self-concordant function and it satisfies Assumption A. In the DiSCO algorithm, we compute the inexact Newton step using (28). Suppose $\frac{2\mu}{\lambda+2\mu}\leq\frac{1}{20}$ and $\|{f^{\prime\prime}_{1}(w_{k})-f^{\prime\prime}(w_{k})}\|_{2}\leq\mu$ for all $k\geq 0$. Then for any $\epsilon>0$, in order to find a solution $\widehat{w}$ satisfying $f(\widehat{w})-f(w_{\star})<\epsilon$, the total number of communication rounds $T$ is bounded by

In Corollary 2, the requirement on $\mu$, which upper bounds $\|{f^{\prime\prime}_{1}(w_{k})-f^{\prime\prime}(w_{k})}\|_{2}$ for all $k\geq 0$, is quite strong. In particular, it requires $\mu$ to be a small fraction of $\lambda$ in order to satisfy $\frac{2\mu}{\lambda+2\mu}\leq\frac{1}{20}$. As we will see from the stochastic analysis in the next section, the spectral bound $\mu$ decreases on the order of $1/\sqrt{n}$. Therefore, in the standard setting where the regularization parameter $\lambda\sim 1/\sqrt{mn}$, the condition in Corollary 2 cannot be satisfied, and the convergence of this simple variant may be slow. In contrast, DiSCO with PCG iterations is much more tolerant of a relatively large $\mu$, and can achieve superlinear convergence with a smaller $\epsilon_{k}$.

Stochastic analysis

From Theorems 2 and 3 of the previous section, we see that the communication complexity of the DiSCO algorithm mainly depends on two quantities: the initial objective gap $f(w_{0})-f(w_{\star})$ and the upper bound $\mu$ on the spectral norms $\|f^{\prime\prime}_{1}(w_{k})-f^{\prime\prime}(w_{k})\|_{2}$ for all $k\geq 0$. As we discussed in Section 3.2, the initial gap $f(w_{0})-f(w_{\star})$ may grow with the number of samples due to the scaling used to make the objective function standard self-concordant. On the other hand, the upper bound $\mu$ may decrease as the number of samples increases based on the intuition that the local Hessians and the global Hessian become similar to each other. In this section, we show how to exploit the stochastic origin of the problem (SAA or ERM, as explained in Section 1) to mitigate the effect of objective scaling and quantify the notion of similarity between local and global Hessians. These lead to improved complexity results.

We focus on the setting of distributed optimization of regularized empirical loss. That is, our goal is to minimize $f(w)=(1/m)\sum_{i=1}^{m}f_{i}(w)$, where

We assume that $z_{i,j}$ are i.i.d. samples from a common distribution. Our theoretical analysis relies on refined assumptions on the smoothness of the loss function $\phi$. In particular, we assume that for any $z$ in the sampling space $\mathcal{Z}$, the function $\phi(\cdot,z)$ has bounded first derivative in a compact set, and its second derivatives are bounded and Lipschitz continuous. We formalize these statements in the following assumption.

There are finite constants $(V_{0},G,L,M)$, such that for any $z\in\mathcal{Z}$:

$\|{\phi^{\prime}(w,z)}\|_{2}\leq G$ for any $\|{w}\|_{2}\leq\sqrt{2V_{0}/\lambda}$;

For the regularized empirical loss in (30), condition $(iii)$ in the above assumption implies $\lambda I\preceq f^{\prime\prime}_{i}(w)\preceq LI$ for $i=1,\ldots,m$, which in turn implies Assumption A.

Recall that the initial point $w_{0}$ is obtained as the average of the solutions to $m$ regularized local optimization problems; see equations (19) and (20). The following lemma shows that expected value of the initial gap $f(w_{0})-f(w_{\star})$ decreases with order $1/\sqrt{n}$ as the local sample size $n$ increases. The proof uses the notion and techniques of uniform stability for analyzing the generalization performance of ERM . See Appendix D for the proof.

Here the expectation is taken with respect to the randomness in generating the i.i.d. data.

Suppose Assumption B holds. For any $r>0$ and any $i\in\{1,\ldots,m\}$, we have with probability at least $1-\delta$,

If $\phi(w,z_{i,j})$ are quadratic functions in $w$, then we have $M=0$ in Assumption B. In this case, Lemma 6 implies $\|{f^{\prime\prime}_{i}(w)-f^{\prime\prime}(w)}\|_{2}\sim\sqrt{1/n}$. For general non-quadratic loss, Lemma 6 implies $\|{f^{\prime\prime}_{i}(w)-f^{\prime\prime}(w)}\|_{2}\sim\sqrt{d/n}$. We use this lemma to obtain an upper bound on the spectral norm of the Hessian distances $\|{f^{\prime\prime}_{1}(w_{k})-f^{\prime\prime}(w_{k})}\|_{2}$, where the vectors $w_{k}$ are generated by Algorithm 1.

Suppose Assumption B holds and the sequence $\{w_{k}\}_{k\geq 0}$ is generated by Algorithm 1. Let r=\Big{(}\frac{2V_{0}}{\lambda}+\frac{2G}{\lambda}\sqrt{\frac{2V_{0}}{\lambda}}\Big{)}^{1/2}. Then with probability at least $1-\delta$, we have for all $k\geq 0$,

Substituting $\|{w_{0}}\|_{2}\leq\sqrt{2V_{0}/\lambda}$ (see Lemma 5) into the above inequality yields

Thus, we have $\|{w_{k}}\|_{2}\leq r$ for all $k\geq 0$. Applying Lemma 6 establishes the corollary. ∎

Here we remark that the dependence on $d$ of the upper bound in (34) comes from Lemma 6, where the bound needs to hold for all point in a $d$-dimensional ball with radius $r$. However, for the analysis of the DiSCO algorithm, we only need the matrix concentration bound to hold for a finite number of vectors $w_{0},w_{1},\dots,w_{K}$, instead of for all vectors satisfying $\|{w}\|_{2}\leq r$. Thus we conjecture that the bound in (34), especially its dependence on the dimension $d$, is too conservative and very likely can be tightened.

We are now ready to present the main results of our stochastic analysis. The following theorem provides an upper bound on the expected number of communication rounds required by the DiSCO algorithm to find an $\epsilon$-optimal solution. Here the expectation is taken with respect to the randomness in generating the i.i.d. data set $\{z_{i,j}\}$.

Then for any $\epsilon>0$, the total number of communication rounds $T$ required to reach $f(\widehat{w})-f(w_{\star})\leq\epsilon$ is bounded by

where $C_{1},C_{2},C_{3}$ are $\widetilde{\mathcal{O}}(1)$ or logarithmic terms:

In particular, ignoring numerical constants and logarithmic terms, we have

Suppose Algorithm 3 terminates in $K$ iterations, and let $t_{k}$ be the number of conjugate gradient steps in each call of Algorithm 2, for $k=0,1,\ldots,K-1$. For any given $\mu>0$, we define $T_{\mu}$ as in (25). Let $\mathcal{A}$ denotes the event that $t_{k}\leq T_{\mu}$ for all $k\in\{0,\ldots,K-1\}$. Let $\bar{\mathcal{A}}$ be the complement of $\mathcal{A}$, i.e., the event that $t_{k}>T_{\mu}$ for some $k\in\{0,\ldots,K-1\}$. In addition, let the probabilities of the events $\mathcal{A}$ and $\bar{\mathcal{A}}$ be $1-\delta$ and $\delta$ respectively. By the law of total expectation, we have

When the event $\mathcal{A}$ happens, we have $T\leq 1+K(T_{\mu}+1)$ where $T_{\mu}$ is given in (25); otherwise we have $T\leq 1+K(T_{L}+1)$, where

is an upper bound on the number of PCG iterations in Algorithm 2 when the event $\bar{\mathcal{A}}$ happens (see the analysis in Appendix F). Since Algorithm 2 always return a $v_{k}$ such that $\|f^{\prime\prime}(w_{k})v_{k}-f^{\prime}(w_{k})\|_{2}\leq\epsilon_{k}$, the outer iteration count $K$ share the same bound (12), which depends on the random variable $f(w_{0})-f(w_{\star})$. However, $T_{\mu}$ and $T_{L}$ are deterministic constants. So we have

Combining with the strong convexity of $f$, we obtain

where the additional 1 counts compensate for removing one $\lceil\cdot\rceil$ operator in (12).

where $C_{0}=1+\left\lceil\log_{2}(2\omega(1/6)/\epsilon)\right\rceil$. With the choice of $\delta$ in (35) and the definition of $T_{L}$ in (36), we have

Replacing $T_{\mu}$ by its expression in (25) and applying Corollary 3, we obtain the desired result. ∎

According to Theorem 4, we need to set the two input parameters $\rho$ and $\mu$ in Algorithm 3 appropriately to obtain the desired communication efficiency. Using the adaptive DiSCO method given in Algorithm 4, we can avoid the explicit specification of $\mu=\mu_{r,\delta}$ defined in (33) and (35). This is formalized in the following theorem.

In Algorithm 4, the parameter $\mu_{k}$ is automatically tuned such that the number of PCG iterations in Algorithm 2 is no more than $T_{\mu_{k}}$. By Corollary 3, with probability at least $1-\delta$, we have

where $\mu_{r,\delta}$ is defined in (33), and $r$ and $\delta$ are given in (35). Therefore we can use the same arguments in the proof of Theorem 4 to show that

and $C_{2}$ and $C_{3}$ are the same as given in Theorem 4. Ignoring constants and logarithmic terms, we obtain the desired result. ∎

To handle the constraint $\|{w}\|_{2}\leq D$, we need to replace the inexact damped Newton method in DiSCO with an inexact proximal Newton method, and replace the distributed PCG method for solving the Newton system with a preconditioned accelerated proximal gradient method. Further details of such an extension are given in Section 7.

The expectation bounds on the rounds of communication given in Theorems 4 and 5 are obtained by combining two consequences of averaging over a large number of i.i.d. local samples. One is the expected reduction of the initial gap $f(w_{0})-f(w_{\star})$ (Lemma 5), which helps to mitigate the effect of objective scaling used to make $f$ standard self-concordant. The other is a high-probability bound that characterizes the similarity between the local and global Hessians (Corollary 3). If the empirical loss $f$ is standard self-concordant without scaling, then we can regard $f(w_{0})-f(w_{\star})$ as a constant, and only need to use Corollary 3 to obtain a high-probability bound. This is demonstrated for the case of linear regression in Section 5.1.

For applications where the loss function needs to be rescaled to be standard self-concordant, the convexity parameter $\lambda$ as well as the “constants” $(V_{0},G,L,M)$ in Assumption B also need to be rescaled. If the scaling factor grows with $n$, then we need to rely on Lemma 5 to balance the effects of scaling. As a result, we only obtain bounds on the expected number of communication rounds. These are demonstrated in Section 5.2 for binary classification with logistic regression and a smoothed hinge loss.

1 Application to linear regression

We consider linear regression with quadratic regularization (ridge regression). More specifically, we minimize the overall empirical loss function

Thus we can apply Theorems 4 and 5 to obtain an expectation bound on the number of communication rounds for DiSCO. For linear regression, however, we can obtain a stronger result.

Since $f$ is a quadratic function, it is self-concordant with parameter , and by definition also standard self-concordant (with parameter $2$). In this case, we do not need to rescale the objective function, and can regard the initial gap $f(w_{0})-f(w_{\star})$ as a constant. As a consequence, we can directly apply Theorem 2 and Corollary 3 to obtain a high probability bound on the communication complexity, which is stronger than the expectation bounds in Theorems 4 and 5. In particular, Theorem 2 states that if

then the number of communication rounds $T$ is bounded as

Since there is no scaling, the initial gap $f(w_{0})-f(w_{\star})$ can be considered as a constant. For example, we can simply pick $w_{0}=0$ and have

By Corollary 3 and the fact that $M=0$ for quadratic functions, the condition (42) holds with probability at least $1-\delta$ if we choose

Further using $L\leq\lambda+2B_{x}^{2}$, we obtain the following corollary.

Suppose we apply DiSCO (Algorithm 3) to minimize $f(w)$ defined in (41) with the input parameter $\mu$ in (43), and let $T$ be the total number of communication rounds required to find an $\epsilon$-optimal solution. With probability at least $1-\delta$, we have

We note that the same conclusion also holds for the adaptive DiSCO algorithm (Algorithm 4), where we do not need to specify the input parameter $\mu$ based on (43). For the adaptive DiSCO algorithm, the bound in (44) holds for any $\delta\in(0,1)$.

The communication complexity guaranteed by Corollary 4 is strictly better than that of distributed implementation of the accelerated gradient method and ADMM (cf. Table 1). If we choose $\lambda=\Theta(1/\sqrt{mn})$, then Corollary 4 implies

with high probability. The DANE algorithm , under the same setting, converges in $\widetilde{\mathcal{O}}(m\log(1/\epsilon))$ iterations with high probability (and each iteration requires two rounds of communication). Thus DiSCO enjoys a better communication efficiency.

2 Application to binary classification

For binary classification, we consider the following regularized empirical loss function

then we have $\phi(w,(x,y))=\eta\varphi(yw^{T}x)$ and $\lambda=\eta\gamma$. It is easy to check that Assumption B holds with

which all containing the scaling factor $\eta$. Plugging these scaled constants into Theorems 4 and 5, we have the following corollary.

For logistic regression, the number of communication rounds required by DiSCO to find an $\epsilon$-optimal solution is bounded by

In the specific case when $\gamma=\Theta(1/\sqrt{mn})$, Corollary 5 implies

If we ignore logarithmic terms, then the expected number of communication rounds is independent of the sample size $n$, and only grows slowly with the number of machines $m$.

Smoothed Hinge Loss

If we choose $p=2+\log(1/\gamma)$, then applying Theorems 4 and 5 yields the following result.

For the smoothed hinge loss $\varphi_{p}$ defined in (8) with $p=2+\log(1/\gamma)$, the total number of communication rounds required by DiSCO to find an $\epsilon$-optimal solution is bounded by

Thus, the smoothed hinge loss enjoys the same communication efficiency as the logistic loss.

Numerical experiments

In this section, we conduct numerical experiments to compare the DiSCO algorithm with several state-of-the-art distributed optimization algorithms: the ADMM algorithm (e.g., ), the accelerated full gradient method (AFG) [34, Section 2.2], the L-BFGS quasi-Newton method (e.g., [37, Section 7.2]), and the DANE algorithm .

The algorithms ADMM, AFG and L-BFGS are well known and each has a rich literature. In particular, using ADMM for empirical risk minimization in a distributed setting is straightforward; see [8, Section 8]. For AFG and L-BFGS, we use the simple distributed implementation discussed in Section 1.1: at each iteration $k$, each machine computes the local gradients $f^{\prime}_{i}(w_{k})$ and sends it to the master machine to form $f^{\prime}(w_{k})=(1/m)\sum_{i=1}^{m}f^{\prime}_{i}(w_{k})$, and the master machine executes the main steps of the algorithm to compute $w_{k+1}$. The iteration complexities of these algorithms stay the same as their classical analysis for a centralized implementation, and each iteration usually involves one or two rounds of communication.

Here we briefly describe the DANE (Distributed Approximate NEwton) algorithm proposed by Shamir et al. . Each iteration of DANE takes two rounds of communication to compute $w_{k+1}$ from $w_{k}$. The first round of communication is used to compute the gradient $f^{\prime}(w_{k})=(1/m)\sum_{i=1}^{m}f^{\prime}_{i}(w_{k})$. Then each machine solves the local minimization problem

and take a second round of communication to compute $w_{k+1}=(1/m)\sum_{i=1}^{m}v_{k+1,i}$. Here $\mu\geq 0$ is a regularization parameter with a similar role as in DiSCO. For minimizing the quadratic loss in (41), the iteration complexity of DANE is $\widetilde{\mathcal{O}}((L/\lambda)^{2}n^{-1}\log(1/\epsilon))$. As summarized in Table 1, if the condition number $L/\lambda$ grows as $\sqrt{mn}$, then DANE is more efficient than AFG and ADMM when $n$ is large. However, the same complexity cannot be guaranteed for minimizing non-quadratic loss functions. According to the analysis in , the convergence rate of DANE on non-quadratic functions might be as slow as the ordinary full gradient descent method.

For comparison, we solve three binary classification tasks using logistic regression. The datasets are obtained from the LIBSVM datasets and summarized in Table 2. These datasets are selected to cover different relations between the sample size $N=mn$ and the feature dimensionality $d$: $N\gg d$ (Covtype ), $N\approx d$ (RCV1 ) and $N\ll d$ (News20 ). For each dataset, our goal is to minimize the regularized empirical loss function:

Among other methods in comparison, we manually tune the penalty parameter of ADMM and the regularization parameter $\mu$ for DANE to optimize their performance. For AFG, we used an adaptive line search scheme to speed up its convergence. For L-BFGS, we adopted the memory size $30$ (number of most recent iterates and gradients stored) as a general rule of thumb suggested in ,

We want to evaluate DiSCO not only on $w_{k}$, but also in the middle of calculating $v_{k}$, to show its progress after each round of communication. To this end, we follow equation (16) to define an intermediate solution $\widehat{w}_{k}^{t}$ for each iteration $t$ of the distributed PCG method (Algorithm 2):

2 Performance evaluation

It is important to note that different algorithms take different number of communication rounds per iteration. ADMM requires one round of communication per iteration. For AFG and L-BFGS, each iteration consists of at least two rounds of communications: one for finding the descent direction, and another one or more for searching the stepsize. For DANE, there are also two rounds of communications per iteration, for computing the gradient and for aggregating the local solutions. For DiSCO, each iteration in the inner loop takes one round of communication, and there is an additional round of communication at the beginning of each inner loop. Since we are interested in the communication efficiency of the algorithms, we plot their progress in reducing the objective value with respect to the number of communication rounds taken.

We plot the performance of ADMM, AFG, L-BFGS, DANE and DiSCO in Figure 2. According to the plots, DiSCO converges substantially faster than ADMM and AFG. It is also notably faster than L-BFGS and DANE. In particular, the convergence speed (and the communication efficiency) of DiSCO is more robust to the number of machines in the distributed system. For $m=4$, the performance of DiSCO is somewhat comparable to that of DANE. As $m$ grows to $16$ and $64$, the convergence of DANE becomes significantly slower, while the performance of DiSCO only degrades slightly. This coincides with the theoretical analysis: the iteration complexity of DANE is proportional to $m$, but the iteration complexity of DiSCO is proportional to $m^{1/4}$.

Since both DANE and DiSCO take a regularization parameter $\mu$, we study their sensitivity to the choices of this parameter. Figure 3 shows the performance of DANE and DiSCO with the value of $\mu$ varying from $10^{-5}$ to $128\times 10^{-5}$. We observe that the curves of DiSCO are relatively smooth and stable. In contrast, the curves of DANE exhibit sharp valley at particular values of $\mu$. This suggests that DiSCO is more robust to the non-optimal choice of parameters.

Extension to distributed composite minimization

Thus far, we have studied the problem of minimizing empirical loss functions that are standard self-concordant. In this section, we sketch how to extend the DiSCO algorithm to solve distributed composite minimization problems. By composite minimization, we consider the minimization of

We modify Algorithm 1 and Algorithm 2 to minimize the composite function $F(w)$. To modify Algorithm 1, we update $w_{k+1}$ using an inexact version of the proximal-Newton method (e.g., ). More specifically, the two steps in each iteration of Algorithm 1 are replaced with:

Find a vector $v_{k}$ that is an approximate solution of

Update $w_{k+1}=w_{k}-\frac{v_{k}}{1+\sqrt{v^{T}f^{\prime\prime}(w_{k})v_{k}}}$.

Note that for $\Psi(w)\equiv 0$, the above proximal-Newton method reduces to Algorithm 1. Since $v_{k}$ only needs to be an inexact solution to problem (48), we need a measure to quantify the approximation error. For this purpose, we define the following gradient mapping

If $v_{k}$ is an exact minimizer of (48), then we have $\|{g(v_{k})}\|_{2}=0$. In the distributed setting, we only need to find a vector $v_{k}$ such that $\|{g(v_{k})}\|_{2}\leq\epsilon_{k}$.

It remains to devise an distributed algorithm to compute an inexact minimizer $v_{k}$. Since the objective function in (48) is not quadratic, we can no longer employ the distributed PCG method in Algorithm 2. Instead, we propose a preconditioned accelerated proximal gradient method. In particular, we modify the algorithm on the master machine in Algorithm 2 as follows:

where $s^{(t+1)}$ is an auxiliary vector. We output $v_{k}=v^{(t+1)}$ once the condition $\|{g(v^{(t+1)})}\|_{2}\leq\epsilon_{k}$ is satisfied. Each update takes one round of communication to compute the vector $f^{\prime\prime}(w_{k})s^{(t)}$. Then, the sub-problem (49) is locally solved by the master machine. This problem has similar structure as problems (20) and (26), and can be solved in \mathcal{O}\bigl{(}(n+\frac{L+\mu}{\lambda+\mu})\log(1/\epsilon)\bigr{)} time using the methods proposed in .

If we replace the first term on the right-hand side of equation (49) by $\frac{L}{2}\|{v-v^{(t)}}\|_{2}^{2}$ and set $\mu=L$, then the above algorithm is exactly the accelerated proximal gradient algorithm , which converges in $\widetilde{\mathcal{O}}(\sqrt{L/\lambda})$ iterations. By utilizing the similarity between $f^{\prime\prime}_{1}(w_{k})$ and $f^{\prime\prime}(w_{k})$, and assuming $\|{f^{\prime\prime}_{1}(w_{k})-f^{\prime\prime}(w_{k})}\|_{2}\leq\mu$ for all $k\geq 0$, it can be shown that our algorithm in (49) converges in $\widetilde{\mathcal{O}}(1+\sqrt{\mu/\lambda})$ iterations, which is of the same order as the PCG algorithm.

In summary, to minimize the composite function $f(w)+\Psi(w)$, we replace Algorithm 1 by the inexact proximal Newton method, and replace Algorithm 2 by a distributed implementation of the above preconditioned accelerated proximal gradient method. Under the same assumptions on $f$, we can obtain similar guarantees on the communication efficiency as stated in Theorems 4 and 5.

Conclusions

We considered distributed convex optimization problems originated from SAA or ERM, which involve large amount of i.i.d. data stored on a distributed computing system. Since the cost of inter-machine communication is very high in practice, communication efficiency is a critical measure in evaluating the performance of a distributed algorithm. For algorithms based on first-order methods, including accelerated gradient methods and ADMM, the required number of communication rounds grows with the condition number of the objective function. The condition number itself often grows with the number of samples due to weaker regularization required. This causes the total number of communication rounds to grow with the overall sample size.

In this paper, we proposed and analyzed DiSCO, a communication-efficient distributed algorithm for minimizing self-concordant empirical loss functions, and discussed its application to linear regression and classification. DiSCO is based on an inexact damped Newton method, where the inexact Newton steps are computed by a distributed preconditioned conjugate gradient method. In a standard setting for supervised learning, its required number of communication rounds does not increase with the sample size, but only grows slowly with the number of machines in the distributed system. There are three main thrusts in our approach:

Self-concordant analysis. We showed that several popular empirical loss functions used in machine learning are either self-concordant or can be well approximated by self-concordant functions. We gave complexity analysis of the inexact damped Newton method, and characterized the conditions for both linear and superlinear convergence.

Preconditioned conjugate gradient (PCG) method. We proposed a distributed implementation of the PCG method for computing the inexact Newton step. In particular, the preconditioner based on similarity between local and global Hessians is very effective in reducing the number of communication rounds, both in theory and practice.

Stochastic analysis of communication efficiency. Our main theoretical results combine two consequences of averaging over a large number of i.i.d. samples. One is the expected reduction of the initial objective value, which counters the effect of objective scaling required to make the objective function standard self-concordant. The other is a high-probability bound that characterizes the similarity between the local and global Hessians.

Our numerical experiments on real datasets confirmed the superior communication efficiency of the DiSCO algorithm. In addition, we also proposed an extension for solving distributed optimization problems with composite empirical loss functions.

Appendix A Proof of Theorem 1

First, we recall the definitions of the two auxiliary functions

which form a pair of convex conjugate functions.

We notice that Step 2 of Algorithm 1 is equivalent to

When inequality (50) holds, Nesterov [34, Theorem 4.1.8] has shown that

Using the definition of functions $\omega$ and $\omega_{*}$, and with some algebraic operations, we obtain

By the second-order mean-value theorem, we have

Using the inequality (11), we can upper bound the second derivative $\omega^{\prime\prime}(t)$ as

Combining the two inequalities above, and using the relation $t^{2}/(1+t)\leq 2\omega(t)$ for all $t\geq 0$, we obtain

In the last inequality above, we used the fact that for any $t\geq 0$ we have $\omega((1-\beta)t)\leq(1-\beta)\omega(t)$, which is the result of convexity of $\omega(t)$ and $\omega(0)=0$; more specifically,

Substituting the above upper bound into inequality (A) yields

With inequality (52), we are ready to prove the statements of the lemma. In particular, Part (a) of the Lemma holds for any $0\leq\beta\leq 1/10$.

For part (b), we assume that $\|{\widetilde{u}_{k}}\|_{2}\leq 1/6$. According to [34, Theorem 4.1.13], when $\|{\widetilde{u}_{k}}\|_{2}<1$, it holds that for every $k\geq 0$,

Combining this sandwich inequality with inequality (52), we have

It is easy to verify that $\omega_{*}(t)-\omega(t)\leq 0.26\,\omega(t)$ for all $t\leq 1/6$, and $(4\beta+5\beta^{2})/(1-\beta)\leq 0.23$ if $\beta\leq 1/20$. Applying these two inequalities to inequality (54) completes the proof.

It should be clear that other combinations of the value of $\beta$ and bound on $\|{\widetilde{u}_{k}}\|_{2}$ are also possible. For example, for $\beta=1/10$ and $\|{\widetilde{u}_{k}}\|_{2}\leq 1/10$, we have $\omega(\|{\widetilde{u}_{k+1}}\|_{2})\leq 0.65\,\omega(\|{\widetilde{u}_{k}}\|_{2})$.

Appendix B Super-linear convergence of Algorithm 1

For any $k\geq 0$, we have $f(w_{k+1})\leq f(w_{k})-\frac{1}{2}\omega(\|{\widetilde{u}_{k}}\|_{2})$.

If $\|{\widetilde{u}_{k}}\|_{2}\leq 1/8$, then we have $\omega(\|{\widetilde{u}_{k+1}}\|_{2})\leq\sqrt{6}\,\omega^{3/2}(\|{\widetilde{u}_{k}}\|_{2})$.

Part (b) suggests superlinear convergence when $\|{\widetilde{u}_{k}}\|_{2}$ is small. This comes at the cost of a smaller approximation tolerance $\epsilon_{k}$ given in (55), compared with (10). Roughly speaking, when $\|{f^{\prime}(w_{k})}\|_{2}$ is relative large, the tolerance $\epsilon_{k}$ in (55) needs to be proportional to $\|{f^{\prime}(w_{k})}\|_{2}^{3/2}$ since $\omega(t)=\mathcal{O}(t)$. When $\|{f^{\prime}(w_{k})}\|_{2}$ is very small, the tolerance $\epsilon_{k}$ in (55) needs to be proportional to $\|{f^{\prime}(w_{k})}\|_{2}^{3}$ because $\omega(t)\sim t^{2}$ as $t\to 0$. In contrast, for linear convergence, the tolerance in (10) is proportional to $\|{f^{\prime}(w_{k})}\|_{2}$.

We start with the inequality (A), and upper bound the last two terms on its right-hand side. Since $\omega^{\prime}(t)=\frac{t}{1+t}<1$, we have

Applying these two bounds to (A), we obtain

Next we bound $\|{\widetilde{u}_{k}-\widetilde{v}_{k}}\|_{2}$ using the approximation tolerance $\epsilon_{k}$ specified in (55),

Combining the above inequality with (56), and using $r_{k}=L^{-1/2}\|{f^{\prime}(w_{k})}\|_{2}\leq\|{\widetilde{u}_{k}}\|_{2}$ with the monotonicity of $\omega(\cdot)$, we arrive at

Part (a) of the theorem follows immediately from inequality (57).

For part (b), we assume that $\|{\widetilde{u}_{k}}\|_{2}\leq 1/8$. Combining (53) with (57), we have

Let $h(t){\,:=\,}\omega_{*}(t)-\omega(t)$ and consider only $t\geq 0$. Notice that $h(0)=0$ and $h^{\prime}(t)=\frac{2t^{2}}{1-t^{2}}<\frac{128}{63}t^{2}$ for $t\leq 1/8$. Thus, we conclude that $h(t)\leq\frac{128}{189}t^{3}$ for $t\leq 1/8$. We also notice that $\omega(0)=0$ and $\omega^{\prime}(t)=\frac{t}{1+t}\geq\frac{8}{9}t$ for $t\leq 1/8$. Thus, we have $\omega(t)\geq\frac{4}{9}t^{2}$ for $t\leq 1/8$. Combining these results, we obtain

Applying this inequality to the right-hand side of (58) completes the proof. ∎

In classical analysis of inexact Newton methods , asymptotic superlinear convergence occurs with $\epsilon_{k}\sim\|{f^{\prime}(w_{k})}\|_{2}^{3/2}$ (in fact with $\epsilon\sim\|{f^{\prime}(w_{k})}\|_{2}^{s}$ for any $s>1$). This agrees with our analysis since $\omega(t)=\mathcal{O}(t)$ when $t$ is not too small. Our result can be very conservative asymptotically because $\omega(t)\sim t^{2}$ as $t\to 0$. However, using $\omega(t)$ and the associated self-concordance analysis, we are able to derive a much better global complexity result.

where $\lceil t\rceil$ denotes the smallest nonnegative integer that is larger than or equal to $t$.

By part (a) of Theorem 6, if $\omega(\|{\widetilde{u}_{k}}\|_{2})\geq 1/8$, then each iteration of Algorithm 1 decreases the function value at least by the constant $\frac{1}{2}\omega(1/8)$. So within at most $K_{1}{\,:=\,}\left\lceil\frac{f(w_{0})-f(w_{\star})}{\frac{1}{2}\omega(1/8)}\right\rceil$ iterations, we are guaranteed to have $\|{\widetilde{u}_{k}}\|_{2}\leq 1/8$.