Distributed Optimization Based on Gradient-tracking Revisited: Enhancing Convergence Rate via Surrogation

Introduction

We study distributed optimization over networks in the form:

Distributed optimization in the form (P) has found a wide range of applications in several areas, including network information processing, telecommunications, multi-agent control, and machine learning. An instance of particular interest to this work is the distributed Empirical Risk Minimization (ERM) whereby the goal is to minimize the average loss over some dataset, distributed across the nodes of the network (cf. Sec. 2.1.2). Letting $\mathcal{D}^{(i)}=\{\mathbf{z}_{1}^{(i)},\ldots,\mathbf{z}_{n}^{(i)}\}$ the dataset of $n$ examples available at node $i$’s side, the local empirical loss reads $f_{i}(\mathbf{x})=1/n\sum_{j=1}^{n}f(\mathbf{x};\mathbf{z}_{j}^{(i)})$, where $f(\mathbf{x};\mathbf{z}_{j}^{(i)})$ measures the fit between the parameter $\mathbf{x}$ and the sample $\mathbf{z}_{j}^{(i)}$. Data sets are usually large and high-dimensional, which makes routing local data to other agents (let alone to a centralized node) infeasible or highly inefficients. Given the cost of communications (especially if compared with the speed of local processing), the challenge in such a network setting is designing communication efficient distributed algorithms.

Motivated by the aforementioned applications, our focus pertains to such a design in two possible settings (one being a special case of the other) : 1) The scenario where no significant relationship can be assumed among the local functions $f_{i}$–this is what the literature of distributed optimization has extensively studied, and will be refereed to as the unrelated setting—and 2) the case where the $f_{i}$’s are related, e.g., because they reflect statistical similarity in the data residing at different nodes. For instance, in the distributed ERM problem above, when data are i.i.d. among machines, one can show that quantities such as the gradients and Hessian matrices of the local functions differ only by $\beta=\mathcal{O}(1/\sqrt{n})$, due to concentrations of measure effects –we will refer to this as $\beta$-related setting (cf. Sec. 2.1.2). If properly exploited in the algorithmic design, such similarity can speed up the optimization/learning process over general purpose optimization algorithms.

Problem (P) in the two settings above has been extensively studied in the centralized environment, including star-networks wherein there is a master node connected to all the other workers. Our interest is in the following (non-accelerated) algorithms:

1) Unrelated setting: (P) can be solved on star-networks employing the standard proximal gradient method: to reach precision $\epsilon>0$ on the objective value, one needs \mathcal{O}\big{(}\kappa_{g}\log(1/\epsilon)\big{)} iterations (which is also the number of communication rounds between the master and the workers), where $\kappa_{g}$ is the condition number of $F$.

2) $\beta$-related setting: When the agents’ functions $f_{i}$ are sufficiently similar, a linear rate proportional to $\kappa_{g}$ may be highly suboptimal. For instance, in the extreme case where all $f_{i}$’s are identical ($\beta=0$), the number of iterations/communications to an $\epsilon>0$ solution would remain the same as for $\beta=\mathcal{O}(L)$. In fact, when $1+\beta/\mu<\kappa_{g}$, faster rates can be obtained exploiting the similarity of the $f_{i}$’s. Specifically, proposed DANE: a mirror-descent type algorithm over star-networks, where each worker $i$ replaces the quadratic term in its local proximal-gradient update with the Bregman divergence of the reference function $f_{i}+\beta/2\|\bullet\|^{2}$; and the master averages the solutions of the workers. DANE is applicable to (P) with $G=0$: For quadratic losses, it achieves an $\epsilon$-solution in \mathcal{O}\big{(}(\beta/\mu)^{2}\cdot\log(1/\epsilon)\big{)} iterations/communications (it is assumed $\beta/\mu\geq 1)$ while no improvement is proved over the proximal gradient if the $f_{i}$’s are not quadratic. More recently, proposed CEASE, which achieves DANE’s rate for (P) with $G\neq 0$ and nonquadratic losses. Using recent results in , it is not difficult to check that the mirror-descent algorithm implemented at the master (thus without averaging workers’ iterates) with the Bregman divergence of $f_{1}+\beta/2\|\bullet\|^{2}$ ($f_{1}$ is the local function at the master) achieves an $\epsilon>0$ solution in \widetilde{\mathcal{O}}\big{(}\beta/\mu\cdot\log(1/\epsilon)\big{)} iterations/communications, improving thus on DANE/CEASE’s rates.

1 Major contributions

We provide the first linear convergence rate analysis of a distributed algorithm, SONATA (Successive cONvex Approximation algorithm over Time-varying digrAphs), applicable to the composite, constrained formulation (P) over (time-varying, directed) graphs. SONATA was earlier proposed in the companion paper for nonconvex problems. It combines the use of surrogate functions in the agents’ subproblems with a perturbed (push-sum) consensus mechanism that aims at locally tracking the gradient of $F$. Surrogate functions replace the more classical first order approximation of the local $f_{i}$’s, which is the omnipresent choice in current distributed algorithms, offering the potential to better suit the geometry of the problem. For instance, (approximate) Newton-type subproblems or mirror descent-type updates naturally fit our surrogate models; they are the key enabler of provably faster rates in the $\beta$-related setting. We comment SONATA’s rates below (cf. Table 3).

Unrelated setting (Table 3): When the network is sufficiently connected or it has a star-topology, SONATA reaches an $\epsilon$-solution on the objective value in \mathcal{O}\big{(}\kappa_{g}\log(1/\epsilon)\big{)} iterations/communications, which matches the rate of the centralized proximal-gradient algorithm. For arbitrary network connectivity, the same iteration complexity is achieved at the cost of $\mathcal{O}((1-\rho)^{-1/2})$ rounds of communications per iteration (employing Chebishev acceleration), where $\rho\in[0,1)$ is the second largest eigenvalue modulus of the mixing matrix. Our rates improve on those of existing distributed algorithms which show a much more pessimistic dependence on the optimization parameters and are proved under more restrictive assumptions–contrast Table 2 with Table 3. Linear rates over time-varying digraphs are reported in Table 4 (cf. Sec. 4.2).

$\beta$-related setting (Table 3): When the agents’ functions are sufficiently similar (specifically, $1+\beta/\mu<\kappa_{g}$), the use of a mirror descent-type surrogate over linearization of the $f_{i}$’s provably yields faster rates, at higher computation costs. This improves on the rate of existing distributed algorithms, which are oblivious of function similarity (cf. Table 2). Notice that this is achieved without exchanging any Hessian matrix over the network but leveraging function homogeneity via surrogation. When customized over star-topologies, SONATA’s rates improve on DANE/CEASE’s ones too.

2 Related works

Early works on distributed optimization aimed at decentralizing the (sub)gradient algorithm. The Distributed Gradient Descent (DGD) was introduced in for unconstrained instances of (P) and in for least squares, bot over undirected graphs. A refined convergence rate analysis of DGD can be found in . Subsequent variants of DGD include the projected (sub)gradient algorithm and the push-sum gradient consensus algorithm , the latter implementable over digraphs. While different, the updates of the agents’ variables in the above algorithms can be abstracted as a combination of one (or multiple) consensus step(s) (weighted average with neighbors variables) and a local (sub)gradient descent step, controlled by a step-size (in some schemes, followed by a proximal operation). A diminishing step-size is used to reach exact consensus on the solution, converging thus at a sublinear rate. With a fixed step-size $\alpha$, linear rate of the iterates is achievable, but it can only converge to a $\mathcal{O}(\alpha)$-neighborhood of the solution .

Several subsequent attempts have been proposed to cope with this speed-accuracy dilemma, leading to algorithms converging to the exact solution while employing a constant step-size. Based upon the mechanism put forth to cancel the steady state error in the individual gradient direction, existing proposals can be roughly organized in three groups, namely: i) primal-based distributed methods leveraging the idea of gradient tracking ; ii) distributed schemes using ad-hoc corrections of the local optimization direction ; and iii) primal-dual-based methods . We elaborate next on these works, focusing on schemes achieving linear rate– Table 1 organizes these schemes based upon the setting their convergence is established while Table 2 reports the explicit expression of the rates.

i) Gradient-tracking-based methods: In these schemes, each agent updates its own variables along a direction that tracks the global gradient $\nabla F$. This idea was proposed independently in the NEXT algorithm for Problem (P) and in AUG-DGM for strongly convex, smooth, unconstrained optimization. The work introduced SONATA, extending NEXT over (time-varying) digraphs. A convergence rate analysis of was later developed in , with considering also (time-varying) digraphs. Other algorithms based on the idea of gradient tracking and implementable over digraphs are ADD-OPT and . Subsequent schemes, , the Push-Pull , and the $\mathcal{AB}$ algorithms, relaxed previous conditions on the mixing matrices used in the consensus and gradient tracking steps over digraphs, which neither need to be row- nor column-stochastic. All the schemes above but NEXT and SONATA are applicable only to smooth, unconstrained instances of (P), with each $f_{i}$ strongly convex. This latter assumption is restrictive in some applications, such as distributed machine learning, where not all $f_{i}$ are strongly convex but $F$ is so.

ii) Ad-hoc gradient correction-based methods: These methods developed specific corrections of the plain DGD direction. Specifically, EXTRA and its variant over digraphs, EXTRA-PUSH , introduce two different weight matrices for any two consecutive iterations as well as leverage history of gradient information. They are applicable only to it smooth, unconstrained problems; when each $f_{i}$ is strongly convex, they generate iterates that converge linearly to the minimizer of $F$. To deal with an additive convex nonsmooth term in the objective, proposed PG-EXTRA, which is thus applicable to (P) over undirected graphs, possibly with different local nonsmooth functions. However, linear convergence is not certified. A different approach is to use a linearly increasing number of consensus steps rather than correcting directly the gradient direction; this has been studied in for unconstrained minimization of smooth, strongly convex $f_{i}$’s over undirected graphs.

iii) Primal-dual methods: A common theme of these schemes is employing a prima-dual reformulation of the original multiagent problem whereby dual variables associated to a properly defined (augmented) Lagrangian function serve the purpose of correcting the plain DGD local direction. Examples of such algorithms include: i) distributed ADMM methods and their inexact implementations ; ii) distributed Augmented Lagrangian-based methods with randomized primal variable updates ; and iii) a distributed dual ascent method employing tracking of the average of the primal variable . All these schemes are applicable only to smooth, unconstrained optimization over undirected graphs, with handling time-varying graphs. The extension of these methods to digraphs seems not straightforward, because it is not clear how to enforce consensus via constraints over directed networks.

To summarize, the above literature review shows that currently there exists no distributed algorithm for the general formulation (P) that provably converges at linear rate to the exact solution, in the presence of a nonsmooth function $G$ or constraints (cf. Table 1); let alone mentioning digraphs. Furthermore, when it comes to the dependence of the rate on the optimization parameters, Table 3 shows that, even restricting to unconstrained, smooth minimization, SONATA’s rates improve on existing ones–in particular, SONATA provably obtains fast convergence if the agents’ objective functions (e.g., data) are sufficiently similar.

3 Paper organization

Sec. 2 introduces the main assumptions on the optimization problem and network, along with some motivating examples from machine learning. The SONATA algorithm over undirected graphs is studied in Sec. 3; in particular, linear convergence is proved in Sec. 3.3, while a detailed discussion on the rate expression and its scalability properties is provided in Sec. 3.4. The case of time-varying, possibly directed, graphs is considered in Sec. 4. Finally, some numerical results supporting our theoretical findings are reported in Sec. 5. The study of SONATA when $F$ is nonconvex can be found in the technical report .

Problem & Network Setting

This section summarizes the assumptions on the optimization problem and network setting. We also introduce a general learning problem over networks, which will be used as case study throughout the paper.

Our algorithmic design and convergence results pertain to two problem settings, namely: i) the one where the local functions $f_{i}$ are generic and unrelated (cf. Sec. 2.1.1), and ii) the case where they are related (cf. Sec. 2.1.2). These two settings are formally introduced below.

Consider the following standard assumption.

The set $\emptyset\neq\mathcal{K}\subseteq^{d}$ is closed and convex;

Each $f_{i}:\mathcal{O}\to$ is twice differentiable on the open set $\mathcal{O}\supseteq\mathcal{K}$ and convex;

$G:\mathcal{K}\to$ is convex possibly nonsmooth.

for some $\mu_{i}\geq 0$ and $0<L_{i}<\infty$. Unlike existing works (cf. Table 1), we do not require each $f_{i}$ to be strongly convex but just $F$ (cf. A3). Also, twice differentiability of $f_{i}$ is not really necessary, but assumed here to simplify our derivations.

Under Assumption A, we define the global conditional number associated to (P):

Related quantities determining the (linear) convergence rate of existing distributed algorithms are (cf. Table 2):

Example 1: Consider the following instance of Problem (P):

which all grow indefinitely as $\texttt{b}/\texttt{a}$ or $m$ increase. $\hfill\square$

In the setting above, our goal is to design linearly convergent distributed algorithms whose iterations complexity is proportional to $\kappa_{g}$, instead of the larger quantities in (3).

1.2 The β𝛽\beta-related setting

This setting considers explicitly the case where the functions $f_{i}$ are similar, in the sense defined below .

The local functions $f_{i}$’s (satisfying Assumption A) are called $\beta$-related if $\left\|\nabla^{2}F(\mathbf{x})-\nabla^{2}f_{i}(\mathbf{x})\right\|_{2}\leq\beta$, for all $\mathbf{x}\in\mathcal{K}$ and some $\beta\geq 0$.

The more similar the $f_{i}$’s, the smaller $\beta$. For arbitrary $f_{i}$’s, $\beta$ is of the order of

The interesting case is when $1+\beta/\mu<<\kappa_{g}$; a specific example is discussed next.

Consider a stochastic learning setting whereby the ultimate goal is to minimize some population objective

To solve (6), the $m$ agents have access only to a finite number, say $N=nm$, of i.i.d. samples from the distribution $\mathcal{P}$, evenly and randomly distributed over the network. Using the notation introduced in Sec. 1, the ERM problem reads:

where $f_{i}$ is regularized empirical loss of agent $i$, $\lambda$-strongly convex. Clearly (7) is an instance of (P), satisfying Assumption A.

For the ERM problems (7) we derive next the associated $\beta/\mu$ and contrasts with $\kappa_{g}$. $\widehat{F}$ is $\lambda$-strongly convex; therefore, we can set $\mu=\lambda$. The optimal choice of $\lambda$ is the one minimizing the statistical error resulting in using $\widehat{\mathbf{x}}$ as proxy for $\mathbf{x}^{\star}$. We have [36, Th. 7], with high probability, F(\widehat{\mathbf{x}})-F({\mathbf{x}}^{\star})\leq\frac{\lambda}{2}\|\boldsymbol{\theta}^{\star}\|^{2}+\mathcal{O}\big{(}\frac{G_{f}^{2}}{\lambda\,N}\big{)}\leq\mathcal{O}\big{(}\lambda\,B^{2}+\frac{G_{f}^{2}}{\lambda\,N}\big{)}, where $G_{f}$ is the Lipschitz constant of $f(\bullet;\mathbf{z})$ on $\mathcal{H}\bigcap\mathcal{B}_{B}$, for all $\mathbf{z}\in\mathcal{Z}$. The optimal choice of $\lambda$ and resulting minimum error rate are then

An estimate of $\beta$ can be obtained exploring the statistical similarity of the local empirical losses $f_{i}$ in (7). Under the additional assumption that $\nabla^{2}f(\bullet;\mathbf{z})$ is $M$-Lipchitz on $\mathcal{H}$, for all $\mathbf{z}\in\mathcal{Z}$, a minor modification of [58, Lemma 6] applied to (6)-(7), yields: with high probability,

where $\widetilde{\mathcal{O}}$ hides the log-factor dependence. Note that when $f(\bullet;\mathbf{z})$ is quadratic (i.e., $M=0$), $\beta$ scales favorably with the dimension $d$.

Based on (8)-(9), an estimate of $\beta/\mu$ and $\kappa_{g}$ for (7) reads:

Note that $\kappa_{g}$ increases with the local sample size $n$ while $\beta/\mu$ does not (neglecting log-factors). It turns out that algorithms converging at a rate depending on $\kappa_{g}$ exhibit a speed-accuracy dilemma: small statistical errors in (8) (larger $n$) are achieved at the cost of more iterations (larger $\kappa_{g}$). In this setting, it is thus desirable to design distributed algorithms whose rate depends on $\beta/\mu$ rather than $\kappa_{g}$.

2 Network setting

We will consider separately two network settings: i) the case where the underlying communication graph is fixed and undirected; and ii) the more general setting of time-varying directed graphs.

When the network of the agent is modeled as a fixed, undirected graph, we write $\mathcal{G}\triangleq(\mathcal{V},\mathcal{E})$, where $\mathcal{V}\triangleq\{1,\ldots,m\}$ denotes the vertex set–the set of agents–while $\mathcal{E}\triangleq\{(i,j)\,|\,i,j\in\mathcal{V}\}$ represents the set of edges–the communication links; $(i,j)\in\mathcal{E}$ iff there exists a communication link between agent $i$ and $j$. We make the following standard assumption on the graph connectivity.

In this setting, communication network is modeled as a time-varying digraph: time is slotted, and at time-frame $\nu$, the digraph reads $\mathcal{G}^{\nu}=\left(\mathcal{V},\mathcal{E}^{\nu}\right)$, where the set of edges $\mathcal{E}^{\nu}$ represents the agents’ communication links: $(i,j)\in\mathcal{E}^{\nu}$ there is a link going from agent $i$ to agent $j$. We make the following standard assumption on the “long-term” connectivity property of the graphs.

Assumption B$\,{}^{\prime}$ (On the network). The graph sequence $\{\mathcal{G}^{\nu}\}$, $\nu=0,1,\ldots$, is $B$-strongly connected, i.e., there exists a finite integer $B>0$ such that the graph with edge set $\cup_{t=\nu B}^{(\nu+1)B-1}\mathcal{E}^{t}$ is strongly connected, for all $\nu=0,1,\ldots$.

The network setting covers, as special case, star-networks, i.e., architectures with a centralized node (a.k.a. master node) connected to all the others (a.k.a. workers). This is the typical computational architecture of several federated learning systems.

The SONATA algorithm over undirected graphs

In words, each agent $i$, given the current iterates $\mathbf{x}_{i}^{\nu}$ and $\mathbf{y}_{i}^{\nu}$, first solves a strongly convex optimization problem wherein $\widetilde{F}_{i}$ is an approximation of the sum-cost $F$ at $\mathbf{x}_{i}^{\nu}$; $\widetilde{f}_{i}$ in (11a) is a strongly convex function, which plays the role of a surrogate of $f_{i}$ (cf. Assumption C below) while $\mathbf{y}_{i}^{\nu}$ acts as approximation of the gradient of $F$ at $\mathbf{x}_{i}^{\nu}$, that is, $\nabla F(\mathbf{x}_{i}^{\nu})\approx\mathbf{y}_{i}^{\nu}$ (see discussion below). Then, agent $i$ updates $\mathbf{x}_{i}^{\nu}$ along the local direction $\mathbf{d}_{i}^{\nu}$ [cf. (11b)], using the step-size $\alpha\in(0,1]$; the resulting point $\mathbf{x}_{i}^{\nu+1/2}$ is broadcast to its neighbors. The update $\mathbf{x}_{i}^{\nu+1/2}\to\mathbf{x}_{i}^{\nu+1}$ is obtained via the consensus step (11c) while the $y$-variables are updated via the perturbed consensus (11d), aiming at tracking $\nabla F(\mathbf{x}_{i}^{\nu})$.

The main assumptions underlying the convergence of SONATA are discussed next.

The surrogate functions satisfy the following conditions.

Each $\widetilde{f}_{i}:\mathcal{O}\times\mathcal{O}\to$ is $C^{2}$ and satisfies

$\nabla\widetilde{f}_{i}(\mathbf{x};\mathbf{x})=\nabla f_{i}(\mathbf{x})$, for all $\mathbf{x}\in\mathcal{K}$;

$\nabla\widetilde{f}_{i}(\bullet;\mathbf{x})$ is $\widetilde{L}_{i}$-Lipschitz continuous on $\mathcal{K}$, for all $\mathbf{x}\in\mathcal{K}$;

$\widetilde{f}_{i}(\bullet;\mathbf{x})$ is $\widetilde{\mu}_{i}$-strongly convex on $\mathcal{K}$, for all $\mathbf{x}\in\mathcal{K}$;

where $\nabla\widetilde{f}_{i}(\mathbf{x};\mathbf{z})$ is the partial gradient of $\widetilde{f}_{i}$ at $(\mathbf{x},\mathbf{z})$ with respect to the first argument.

The assumption states that $\widetilde{f}_{i}$ should be regarded as a surrogate of $f_{i}$ that preserves at each iterate $\mathbf{x}^{\nu}_{i}$ the first order properties of $f_{i}$. Conditions (i)-(iii) are certainly satisfied if one uses the classical linearization of $f_{i}$, that is,

with $\tau_{i}>0$. Using (13), one can rewrite (11a) as:

which can be interpreted as a mirror-descent update (with step-size one) for the composite minimization of $g(\mathbf{x}_{i})\triangleq f_{i}(\mathbf{x}_{i})+(\mathbf{y}_{i}^{\nu})^{\top}(\mathbf{x}_{i}-\mathbf{x}_{i}^{\nu})$, based on the Bregman distance associated with the reference function $\omega(\mathbf{x}_{i})\triangleq f_{i}(\mathbf{x}_{i})+{{\tau_{i}}}/{2}\|\mathbf{x}_{i}\|^{2}$.

We refer the reader to as good sources of examples of nonlinear surrogates satisfying Assumption C; here we only anticipate that, when the $f_{i}$’s are sufficiently similar, higher order models such as (13) yield indeed faster rates of SONATA than those achievable using linear surrogates (12). Further intuition is provided next.

For instance, (15) holds with $D_{i}=\max\{|\widetilde{\mu}_{i}-L|,|\widetilde{L}_{i}-\mu|\}$. Roughly speaking, the smaller $D_{i}$ the better $\widetilde{F}_{i}$ in (11a) approximates $F$. To see this, compare $F$ and $\widetilde{F}_{i}$ up to the second order: there exist $\theta_{1},\theta_{2}\in(0,1)$ such that

Noting that $\nabla\widetilde{f}_{i}(\mathbf{x}_{i}^{\nu};\mathbf{x}_{i}^{\nu})=\nabla f_{i}(\mathbf{x}_{i}^{\nu})$ [Assumption C(i)] and $\nabla\widetilde{F}_{i}(\mathbf{x}_{i}^{\nu};\mathbf{x}_{i}^{\nu})=\mathbf{y}_{i}^{\nu}$, and anticipating $\|\nabla F(\mathbf{x}_{i}^{\nu})-\mathbf{y}_{i}^{\nu}\|\to 0$ as $\nu\to\infty$ (see discussion below), it follows that $\widetilde{F}_{i}$ approximates $F$ asymptotically, up to the first order. A better match, is achieved when $D_{i}$ is sufficiently small. One can then expect that, if the local functions are sufficiently similar ($\beta$ is small), surrogates $\widetilde{f}_{i}$ exploiting higher order information of $f_{i}$, such as (13), may be more effective than mere linearization. Our theoretical findings confirm the above intuition–see Sec. 3.4.

In the consensus and tracking steps, the weights $w_{ij}$’s satisfy the following standard assumption.

The weight matrix $\mathbf{W}\triangleq(w_{ij})_{i,j=1}^{m}$ has a sparsity pattern compliant with $\mathcal{G}$, that is

$w_{ij}>0$, if $(i,j)\in\mathcal{E}$; and $w_{ij}=0$ otherwise;

Furthermore, $\mathbf{W}$ is doubly stochastic, that is, $\mathbf{1}^{\top}\mathbf{W}=\mathbf{1}^{\top}$ and $\mathbf{W}\mathbf{1}=\mathbf{1}$.

Several rules have been proposed in the literature compliant with Assumption D, such as the Laplacian, the Metropolis-Hasting, and the maximum-degree weights rules .

Finally, we comment the anticipated gradient tracking property of the $y$-variables, that is, $\|\nabla F(\mathbf{x}_{i}^{\nu})-\mathbf{y}_{i}^{\nu}\|\to 0$ as $\nu\to\infty$. Define the average processes

Summing (11d) over $i\in[m]$ and invoking the doubly stochasticity of $\mathbf{W}$; we have

Applying (18) inductively and using the initial condition $\mathbf{y}^{0}_{i}=\nabla f_{i}(\mathbf{x}_{i}^{0})$, $i\in[m],$ yield

That is, the average of all the $\mathbf{y}_{i}^{\nu}$’s in the network is equal to that of the $\nabla f_{i}(\mathbf{x}_{i}^{\nu})$’s, at every iteration $\nu$. Assuming that consensus on $\mathbf{x}_{i}^{\nu}$’s and $\mathbf{y}_{i}^{\nu}$’s is asymptotically achieved, that is, $\|\mathbf{x}_{i}^{\nu}-\mathbf{x}_{j}^{\nu}\|\underset{\nu\to\infty}{\longrightarrow}0$ and $\|\mathbf{y}_{i}^{\nu}-\mathbf{y}_{j}^{\nu}\|\underset{\nu\to\infty}{\longrightarrow}0$, $i\neq j$, (19) would imply the desired gradient tracking property $\|\nabla F(\mathbf{x}_{i}^{\nu})-\mathbf{y}_{i}^{\nu}\|\to 0$ as $\nu\to\infty$, for all $i\in[m]$.

1 A special instance: SONATA on star-networks

Although the main focus of the paper is the study of SONATA over meshed-networks, it is worth discussing here its special instance over star networks. Specifically, consider a star (unidirected) graph with $m$ nodes, where one of them (the master node) connects with all the others (workers). The workers still own only one function $f_{i}$ of the sum-cost $F$. Two common approaches developed in the literature to solve (P) in this setting are: (i) based upon receiving the gradients $\nabla f_{i}$ from the workers, the master solves (P) and broadcasts the updated vector variables to the workers; (ii) based upon receiving the full gradient $\nabla F$ and the current iterate from the master, all the workers solve locally an instance of (P) and send their outcomes to the master that averages them out, producing then the new iterate. Here we follow the latter approach; the algorithm is described in Algorithm 2, which corresponds to SONATA (up to a proper initialization), with weight matrix $\mathbf{W}=\left[\mathbf{1},\,\mathbf{0}_{m,m-1}\right]\left[\mathbf{1}/m,\,\mathbf{0}_{m,m-1}\right]^{\top}$.

SONATA-star, employing linear surrogates [cf. (12)] and $\alpha=1$, reduces to the proximal gradient algorithm. When the surrogates (13) are used (and still $\alpha=1$), SONATA-star coincides with the DANE algorithm if $G=0$ and to the CEASE (with averaging) algorithm if $G\neq 0$. Nevertheless, our convergence rates improve on those of DANE and CEASE–see Sec. 3.4.1.

2 Intermediate definitions

We conclude this section introducing some quantities that will be used in the rest of the paper. We define the optimality gap as

where $\mathbf{x}^{\star}$ is the unique solution of Problem (P).

We stack the local variables and gradients in the column vectors

The average of each of the vectors above is defined as $\bar{\mathbf{x}}^{\nu}\triangleq(1/m)\cdot\sum_{i=1}^{m}\mathbf{x}_{i}^{\nu}$. The consensus disagreements on $\mathbf{x}_{i}^{\nu}$’s and $\mathbf{y}^{\nu}_{i}$’s are

respectively, while the gradient tracking error is defined as

Finally, given the weight matrix $\mathbf{W}$, we define

Under Assumptions B and D, it is well known that (see, e.g., )

where $\sigma(\bullet)$ denotes the largest singular value of its argument.

3 Linear convergence rate

Our proof of linear rate of SONATA passes through the following steps. Step 1: We begin showing that the optimality gap $p^{\nu}$ converges linearly up to an error of the order of $\mathcal{O}(\|\mathbf{x}_{\bot}^{\nu}\|^{2}+\|\mathbf{y}_{\bot}^{\nu}\|^{2})$, see Proposition 3.4. Step 2 proves that $\|\mathbf{x}_{\bot}^{\nu}\|$ and $\|\mathbf{y}_{\bot}^{\nu}\|$ are also linearly convergent up to an error $\mathcal{O}(\|\mathbf{d}^{\nu}\|)$, see Proposition 3.5. In Step 3 we close the loop establishing $\|\mathbf{d}^{\nu}\|=\mathcal{O}(\sqrt{p^{\nu}}+\|\mathbf{y}_{\bot}^{\nu}\|)$, see Proposition 3.6. Finally, in Step 4, we properly chain together the above inequalities (cf. Proposition 3.8), so that linear rate is proved for the sequences $\{p^{\nu}\}$, $\{\|\mathbf{x}_{\bot}^{\nu}\|^{2}\}$, $\{\|\mathbf{y}_{\bot}^{\nu}\|^{2}\}$, and $\{\|\mathbf{d}^{\nu}\|^{2}\}$–see Theorems 3.9 and 3.10. We will tacitly assume that Assumptions A, B, C, and D are satisfied.

Invoking the convexity of $U$ and the doubly stochasticity of $\mathbf{W}$, we can bound $p^{\nu+1}$ as

We can now bound $U(\mathbf{x}_{j}^{\nu+\frac{1}{2}})$, regarding the local optimization (11a)-(11b) as a perturbed descent on the objective, whose perturbation is due to the tracking error $\boldsymbol{\delta}^{\nu}$. In fact, Lemma 3.1 below shows that, for sufficiently small $\alpha$, the local update (11b) will decrease the objective value $U$ up to some error, related to $\boldsymbol{\delta}_{i}^{\nu}$.

Let $\{\mathbf{x}_{i}^{\nu}\}$ be the sequence generated by SONATA; there holds:

where $\mathbf{H}\triangleq\int_{0}^{1}(1-\theta)\nabla^{2}F(\theta\mathbf{x}_{i}^{\nu+\frac{1}{2}}+(1-\theta)\mathbf{x}_{i}^{\nu})d\theta$.

Invoking the optimality of $\widehat{\mathbf{x}}_{i}^{\nu}$ and defining $\widetilde{\mathbf{H}}_{i}\triangleq\int_{0}^{1}\nabla^{2}\widetilde{f}_{i}(\theta\,\widehat{\mathbf{x}}_{i}^{\nu}+(1-\theta)\,\mathbf{x}_{i}^{\nu};\mathbf{x}_{i}^{\nu})d\theta$, we have

where the equality follows from $\nabla\widetilde{f}_{i}(\mathbf{x}_{i}^{\nu};\mathbf{x}_{i}^{\nu})=\nabla f_{i}(\mathbf{x}_{i}^{\nu})$ and the integral form of the mean value theorem. Substituting (30) in (29) and using the convexity of $G$ yield

It remains to bound $\alpha\mathbf{H}-\widetilde{\mathbf{H}}_{i}$. We proceed as follows:

We can now substitute (28) into (27) and get

where in (a) we used Young’s inequality, with $\epsilon_{opt}>0$ satisfying

Next we lower bound $\|\mathbf{d}^{\nu}\|^{2}$ in terms of the optimality gap.

The following lower bound holds for $\|\mathbf{d}^{\nu}\|^{2}$:

where $D_{\operatorname{mx}}$ is defined in (24).

Invoking the optimality condition of $\widehat{\mathbf{x}}_{i}^{\nu}$, yields

Using the $\mu$-strong convexity of $F$, we can write

Rearranging the terms and summing over $i\in[m]$, yields

Using (27) in conjunction with $U(\mathbf{x}_{i}^{\nu+\frac{1}{2}})\leq\alpha U(\widehat{\mathbf{x}}_{i}^{\nu})+(1-\alpha)U(\mathbf{x}_{i}^{\nu})$ leads to

Combining (37) with (38) provides the desired result (35). ∎

As last step, we upper bound $\|\boldsymbol{\delta}^{\nu}\|^{2}$ in (33) in terms of the consensus errors $\|\mathbf{x}_{\bot}^{\nu}\|^{2}$ and $\|\mathbf{y}_{\bot}^{\nu}\|^{2}$.

The following upper bound holds for the tracking error $\|\boldsymbol{\delta}^{\nu}\|^{2}$:

where $L_{\operatorname{mx}}$ is defined in (4).

We are ready to prove the linear convergence of the optimality gap up to consensus errors. The result is summarized in Proposition 3.4 below. The proof follows readily multiplying (33) and (35) by $\widetilde{\mu}_{\operatorname{mn}}-\frac{L}{2}\alpha-\frac{1}{2}\epsilon_{opt}$ and ${6(L^{2}+\widetilde{L}_{\operatorname{mx}}^{2})}/{\mu}$, respectively, adding them together to cancel out $\|\mathbf{d}^{\nu}\|$, and using (39) to bound $\|\boldsymbol{\delta}^{\nu}\|^{2}$.

The optimality gap $p^{\nu}$ [cf. (20)] satisfies

where $\sigma(\alpha)\in(0,1)$ and $\eta(\alpha)>0$ are defined as

We upper bound $\|\mathbf{x}_{\bot}^{\nu}\|$ and $\|\mathbf{y}_{\bot}^{\nu}\|$ in terms of $\|\mathbf{d}^{\nu}\|$. We begin rewriting the SONATA algorithm (11a)-(11d) in vector-matrix form; using (21) and (25),we have

Noting that $\mathbf{x}_{\bot}^{\nu}=(\mathbf{I}-\mathbf{J})\mathbf{x}^{\nu}$ [similarly, $\mathbf{y}_{\bot}^{\nu}=(\mathbf{I}-\mathbf{J})\mathbf{y}^{\nu}$] and $(\mathbf{I}-\mathbf{J})\widehat{\mathbf{W}}=\widehat{\mathbf{W}}-\mathbf{J}$ (due to the doubly stochasticity of $\mathbf{W}$), it follows from (43) that

Using (44)-(45), Proposition 3.5 below establishes linear convergence of the consensus errors $\mathbf{x}_{\bot}^{\nu}$ and $\mathbf{y}_{\bot}^{\nu}$, up to a perturbation.

with $\rho$ and $L_{\operatorname{mx}}$ defined in (26) and (4), respectively.

We prove next (46b); (46a) follows readily from (44). Using (43a), (45), and the Lipschitz continuity of $\nabla f_{i}$ [cf. (1)], we can bound $\|\mathbf{y}_{\bot}^{\nu+1}\|$ as

where in the last inequality we used $\|\mathbf{W}\|\leq 1$. ∎

The following upper bound holds for $\|\mathbf{d}^{\nu}\|$:

where $L_{\operatorname{mx}}$ and $\widetilde{L}_{\operatorname{mx}}$, $\widetilde{\mu}_{\operatorname{mn}}$, $D_{\operatorname{mx}}$ are defined in (4) and (24), respectively.

By optimality of $\widehat{\mathbf{x}}_{i}^{\nu}$ and $\mathbf{x}^{\star}$ we have

Summing the two inequalities above yields

Rearranging terms and using the reverse triangle inequality we obtain the following bound for $\|\mathbf{d}_{i}^{\nu}\|$:

3.4 Step 4: Proof of the linear rate (chaining the inequalities)

We are now ready to prove linear rate of the SONATA algorithm. We build on the following intermediate result, introduced in .

Given the sequence $\{s^{\nu}\}$, define the transformations

for $z\!\in\!(0,1)$. If $S(z)$ is bounded, then $|s^{\nu}|=\mathcal{O}(z^{\nu})$.

We show next how to chain the inequalities (40), (46) and (47) so that Lemma 3.7 can be applied to the sequences $\{p^{\nu}\}$, $\{\|\mathbf{x}_{\bot}^{\nu}\|^{2}\}$, $\{\|\mathbf{y}_{\bot}^{\nu}\|^{2}\}$ and $\{\|\mathbf{d}^{\nu}\|^{2}\}$, establishing thus their linear convergence.

​​ Let $P^{K}(z)$, $X_{\bot}^{K}(z)$, $Y_{\bot}^{K}(z)$ and $D^{K}(z)$ denote the transformation (49) applied to the sequences $\{p^{\nu}\}$,​ $\{\|\mathbf{x}_{\bot}^{\nu}\|^{2}\}$, $\{\|\mathbf{y}_{\bot}^{\nu}\|^{2}\}$ and $\{\|\mathbf{d}^{\nu}\|^{2}\}$, respectively. Given the constants $\sigma(\alpha)$ and $\eta(\alpha)$ (defined in Proposition 3.4) and the free parameters $\epsilon_{x},\epsilon_{y}>0$ (to be determined), the following hold

Squaring (46) and using Young’s inequality yield

for arbitrary $\epsilon_{x},\epsilon_{y}>0$. The proof is completed by taking the maximum of both sides of (40), (47), and (53) over $\nu=0,\ldots,K$ and using $\max_{\nu=0,\ldots,K}|s^{\nu+1}|z^{-\nu}\geq z\cdot\max_{\nu=0,\ldots,K}|s^{\nu}|\,z^{-\nu}-z\cdot|s^{0}|$, for any sequence $\{s^{\nu}\}$ and $z\in(0,1)$. ∎

Chaining the inequalities in Proposition 3.8 in the way shown in Fig. 1, we can bound $D^{K}(z)$ as (see Appendix A for the proof)

where $\mathcal{P}(\alpha,z)$ is defined as

and $\mathcal{R}(\alpha,z)$ is a remainder, which is bounded under (51).

Therefore, as long as $\mathcal{P}(\alpha,z)<1$, (54) implies

where $B$ is a constant independent of $K$. Therefore, $D(z)\leq B$ and thus $\{\|\mathbf{d}^{\nu}\|^{2}\}$ converges R-linearly to zero at rate at least $z$ (cf. Lemma 3.7). Applying the same argument to the other inequalities in Proposition 3.8, one can conclude that also the sequences $\{p^{\nu}\}$, $\{\|\mathbf{x}_{\bot}^{\nu}\|^{2}\}$ and $\{\|\mathbf{y}_{\bot}^{\nu}\|\}$ converge R-linearly to zero.

The last step consists to showing that there exist a sufficiently small step-size $\alpha\in(0,1]$ and $z\in(0,1)$ satisfying (51), such that $\mathcal{P}(\alpha,z)<1$. This is proved in the Theorem 3.9 below.

The proof is organized in following two steps: Step 1) We first consider the “marginal” stable case by letting $z=1$, and show that there exists $\bar{\alpha}>0$ so that $\mathcal{P}(\alpha,1)<1$, for all $\alpha\in(0,\bar{\alpha})$; Step 2) Then, invoking the continuity of $\mathcal{P}(\alpha,z)$, we argue that, for any $\alpha\in(0,\bar{\alpha})$, one can find $\bar{z}(\alpha)<1$ such that \mathcal{P}\big{(}\alpha,\bar{z}(\alpha)\big{)}<1. This implies the boundedness of D^{K}\big{(}\bar{z}(\alpha)\big{)}, and thus \|\mathbf{d}^{\nu}\|^{2}=\mathcal{O}\big{(}\bar{z}(\alpha)^{\nu}\big{)} (cf. Lemma 3.7).

$\bullet$ Step 1: We begin optimizing the free parameters $\epsilon_{x}$, $\epsilon_{y}$, and $\epsilon_{opt}$. Since the goal is to find the largest $\bar{\alpha}$ so that $\mathcal{P}(\alpha,1)<1$, for all $\alpha\in(0,\bar{\alpha})$, the optimal choice of $\epsilon_{x}$, $\epsilon_{y}$, and $\epsilon_{opt}$ is the one that minimizes $\mathcal{P}(\alpha,1)$, that is,

We then set $\epsilon_{x}=\epsilon_{y}=\epsilon^{\star}$, and proceed to optimize $\epsilon_{opt}$, which appears in $\eta(\alpha)$ and $\sigma(\alpha)$. Recalling the definition of $\eta(\alpha)$ and $\sigma(\alpha)$ (cf. Proposition 3.4) and the constraint (34), the problem boils down to minimize

Let $\mathcal{P}^{\star}(\alpha,z)$ denote the value of $\mathcal{P}(\alpha,z)$ corresponding to the optimal choice of the above parameters. The expression of $\mathcal{P}^{\star}(\alpha,1)$ reads

The next theorem provides an explicit expression of the convergence rate in Theorem 3.9 in terms of the step-size $\alpha$; the constants $J$, $A_{\frac{1}{2}}$, and $\alpha^{*}$ therein are defined in (103), (101) with $\theta=1/2$, and (105), respectively.

In the setting of Theorem 3.9, suppose that the step-size $\alpha$ satisfies $\alpha\in(0,\alpha_{\operatorname{mx}})$, with $\alpha_{\operatorname{mx}}\triangleq\min\{(1-\rho)^{2}/A_{\frac{1}{2}},\widetilde{\mu}_{\operatorname{mn}}/(\widetilde{\mu}_{\operatorname{mn}}-D_{\operatorname{mn}}),1\}.$ Then, $U(\mathbf{x}_{i}^{\nu})-U^{\star}=\mathcal{O}(z^{\nu})$, for all $i\in[m]$, where

4 Discussion

Theorem 3.10 provides a unified set of convergence conditions for different choices of surrogates and network topologies. To shed light on the expression of the rate and its dependence on the key optimization and network parameters, we customize here Theorem 3.10 to specific network topologies and surrogate functions. We begin considering star-networks (cf. Sec. 3.4.1) and then move to general graph topologies with no master node (cf. Sec. 3.4.2). We will customize the rate achieved by SONATA employing the following two surrogate functions $\widetilde{f}_{i}$, representing the two extreme choices in the spectrum of admissible surrogates:

Convergence of SONATA-Star (Algorithm 2) is established in Corollary 3.11 below.

In particular, when the surrogates (62) and (63) are employed along with $\alpha=1$, the rate above reduces to the following expressions:

Linearization (62): $z\leq 1-\kappa_{g}^{-1}$. Therefore, $U(\mathbf{x}^{\nu})-U^{\star}\leq\epsilon$ in at most \mathcal{O}\Big{(}\kappa_{g}\log({1}/{\epsilon})\Big{)} iterations (communications);

Therefore, $U(\mathbf{x}^{\nu})-U^{\star}\leq\epsilon$ in at most

The following comments are in order. When linearization is employed, SONATA-Star matches the iteration complexity of the centralized proximal-gradient algorithm. When the $f_{i}$’s are sufficiently similar, (65)-(66) proves that faster rates can be achieved if surrogates (63) are chosen over first-order approximations: when $\beta\ll L$, (66) is significantly faster than \mathcal{O}\big{(}\kappa_{g}\log({1}/{\epsilon})\big{)}. As case study, consider Example 2 (cf. Sec. 2.1.2): plugging (10) into Corollary 3.11 shows that using the surrogates (63) yields \widetilde{\mathcal{O}}\big{(}L\,\sqrt{{d\,m}}\cdot\log(1/\epsilon)\big{)} iterations (communications); this contrasts with \widetilde{\mathcal{O}}\big{(}L\,\sqrt{d\,m\,n}\cdot\log(1/\epsilon)\big{)}, achieved by first-order methods (and SONATA-Star using linearization), which instead increases with the sample size $n$.

Since SONATA-Star contains as special cases the DANE and CEASE algorithms, we contrast here Corollary 3.11 with their convergence rates. We recall that DANE is applicable to (P) when $G=0$: For quadratic losses, it achieves an $\epsilon$-optimal objective value in \mathcal{O}\big{(}(\beta/\mu)^{2}\cdot\log(1/\epsilon)\big{)} iterations/communications (here $\beta/\mu\geq 1$). This rate is worse than (66). For nonquadratic losses, did not show any rate improvement of DANE over plain gradient algorithms, i.e., \mathcal{O}\big{(}\kappa_{g}\cdot\log(1/\epsilon)\big{)} while SONATA-star still retains \mathcal{O}\big{(}\beta/\mu\cdot\log(1/\epsilon)\big{)}. The CEASE algorithm is proved to achieve an $\epsilon$-solution on the iterates in \mathcal{O}\big{(}(\beta/\mu)^{2}\cdot\log(1/\epsilon)\big{)} iterations/communications (with $\beta/\mu\geq 1$); SONATA reaches the same error on the iterates in \mathcal{O}\big{(}\beta/\mu\cdot\log(\kappa_{g}/\epsilon)\big{)} iterations/communications, which matches the order of the mirror-decent algorithm.

In the next section we extend the study to networks with no centralized nodes, sheding lights on the role of the network in achieving the same kind of results.

4.2 The general case

The convergence rate of SONATA over general graphs is summarized in Corollary 3.12 for the linearization surrogates (62) while Corollaries 3.13 and 3.14 consider the surrogates (63) based on local $f_{i}$, with Corollary 3.13 addressing the case $\beta\leq\mu$ and Corollary 3.14 the case $\beta>\mu$. The step-size $\alpha$ is tuned to obtain favorable rate expressions.

In the setting of Theorem 3.10, let $\{\mathbf{x}^{\nu}\}$ be the sequence generated by SONATA, using the surrogates (62) and step-size $\alpha=c\cdot\alpha_{\operatorname{mx}}$, $c\in(0,1)$, with $\alpha_{\operatorname{mx}}=\min\{1,(1-\rho)^{2}/(\rho\cdot 110\kappa_{g}(1+\beta/L)^{2})\}$. The number of iterations (communications) needed for $U(\mathbf{x}_{i}^{\nu})-U^{\star}\leq\epsilon$, $i\in[m]$, is

Instate assumptions of Theorem 3.10 and suppose $\beta\leq\mu$. Consider SONATA using the surrogates (63) and step-size $\alpha=c\cdot\alpha_{\operatorname{mx}}$, $c\in(0,1)$, with $\alpha_{\operatorname{mx}}=\min\{1,(1-\rho)^{2}/(M\rho)\}$ and $M=193\left(1+\frac{\beta}{\mu}\right)^{2}\left(\kappa_{g}+\frac{\beta}{\mu}\right)^{2}$. The number of iterations (communications) needed for $U(\mathbf{x}_{i}^{\nu})-U^{\star}\leq\epsilon$, $i\in[m]$, is

Instate assumptions of Theorem 3.10 and suppose $\beta>\mu$. Consider SONATA using the surrogates (63) and step-size $\alpha=c\cdot\alpha_{\operatorname{mx}}$, $c\in(0,1)$, with $\alpha_{\operatorname{mx}}=\min\{1,(1-\rho)^{2}/(M\rho)\}$ and $M=253\left(1+\frac{L}{\beta}\right)\left(\kappa_{g}+\frac{\beta}{\mu}\right)$. The number of iterations (communications) needed for $U(\mathbf{x}_{i}^{\nu})-U^{\star}\leq\epsilon$, $i\in[m]$, is

The proof of Corollaries 3.13 and 3.14 can be found in Appendix E.

$\bullet$ Order of the rate of centralized (nonaccelerated) methods (Case I): For a fixed optimization problem, if the network is sufficiently connected ($\rho$ “small”), its impact on the rate becomes negligible (the bottleneck is the optimization), and SONATA matches the network-independent rate order achieved on star-topologies (cf. Corollary 3.11) by the proximal gradient algorithm when linearization is employed [cf. (67)] and by the mirror-descent scheme when the local $f_{i}$’s are used in the surrogates [cf. (69) and (71)].

$\bullet$ Network-dependent rates (Case II): As expected, the convergence rate deteriorates as $\rho$ increases, i.e., the network connectivity gets worse. This translates in a less favorable dependence of the complexity on $\kappa_{g}$ and $\beta/\mu$ (by a square factor) and network scalability of the order of $\rho/(1-\rho)^{2}$. When $\beta\sqrt{\rho}=\mathcal{O}(L)$ (e.g., the network is decently connected or $\beta=\mathcal{O}(L)$), the complexity becomes $\mathcal{O}\left(\kappa_{g}^{2}(1-\rho)^{-2}\log(1/\epsilon)\right)$, which compares favorably with that of existing distributed schemes, determined instead by the more pessimistic local quantities (3). The scalability of the rate with the network connectivity, $(1-\rho)^{-2}$, can be improved leveraging multiple rounds of communications or accelerated consensus protocols, as discussed below.

$\bullet$ Linearization (62) vs. local $f_{i}$ (63) surrogates: As already observed in the setting of star-networks, the use of the local losses as surrogates employs a form of preconditioning in the local agents subproblems. When the $f_{i}$’s are sufficiently similar to each other, so that $1+\beta/\mu<\kappa_{g}$, exploiting local Hessian information via (63) provably reduces the iteration/communication complexity over linear models (62)–contrast (67) with (69) and (71). Note that these faster rates are achieved without exchanging any matrices over the network, which is a key feature of SONATA. On the other hand, when the functions $f_{i}$ are heterogeneous, the local surrogates (63) are no longer informative of the average-loss $F$ and using linearization might yield better rates. Although these design recommendations are based on sufficient conditions, numerical results seem to confirm the above conclusions–see Sec. 5.

$\bullet$ Multiple communications rounds and acceleration: The discussion above shows that rates of the order of those of centralized methods can be achieved if the network is sufficiently connected (Case I). When this is not the case, one can still achieve the same iteration complexity at the cost of multiple, finite, rounds of communications per iteration. Specifically, let $\rho_{0}$ be the connectivity of the given network and suppose we run $K$ steps of communications per iteration (computation) in (43a)-(43b); this yields an effective network with improved connectivity $\rho=\rho_{0}^{K}$. One can then choose $K$ so that the ratio $\rho_{0}^{K}/(1-\rho_{0}^{K})^{2}$ satisfies the condition triggering Case I in the Corollaries 3.12–3.14, as briefly summarized next.

1) Linearization: Invoking Corollary 3.12, one can check that the order of such a $K$ is $K=\mathcal{O}(\log(\kappa_{g}(1+\beta/L)^{2})/\log(1/\rho_{0}))=\mathcal{O}(\log(\kappa_{g}(1+\beta/L)^{2})/(1-\rho_{0}))$; therefore, SONATA using the surrogates (62) reaches an $\epsilon$-solution in $\mathcal{O}\left(\kappa_{g}\log(1/\epsilon)\right)$ iterations and $\mathcal{O}\left(\kappa_{g}\cdot(1-\rho_{0})^{-1}\log(\kappa_{g}(1+\beta/L)^{2})\log(1/\epsilon)\right)$ communications. The dependence on the network connectivity $\rho_{0}$ can be further improved leveraging Chebyshev polynomials (see, e.g., ): the final communication complexity of SONATA reads

2) Local $f_{i}$ surrogates: Considering the case $\beta\geq\mu$ (Corollary 3.14), we can show that SONATA using the surrogates (63) and employing multiple rounds of communications per iteration, reaches an $\epsilon$-solution in $\mathcal{O}\left(\beta/\mu\cdot\log(1/\epsilon)\right)$ iterations and \mathcal{O}\left(\beta/\mu\cdot\log\big{(}(\kappa_{g}+\beta/\mu)(1+L/\beta)\big{)}(1-\rho_{0})^{-1}\log(1/\epsilon)\right) communications. If Chebyshev polynomials are used to accelerate the communications, the communication complexity further improves to

The SONATA algorithm over directed time-varying graphs

In this section we extend SONATA and its convergence analysis to solve Problem (P) over directed, time-varying graphs (Assumption B$\,{}^{\prime}$). Note that (11a)-(11d) is not readily applicable to this setting, as constructing a doubly stochastic weight matrix compliant with a directed graph is generally infeasible or computationally costly–see e.g. . Conditions on the weight matrices can be relaxed if the consensus/tracking schemes (11c)-(11d) are properly changed to deal with the lack of doubly stochasticity.

Here, we consider the perturbed push-sum protocols as proposed in the companion paper (but in the Adapt-Then-Combine (ATC) form). The resulting distributed algorithm, still termed SONATA, is formally described in Algorithm 3.

In the perturbed push-sum protocols (73c)-(73d), $\mathbf{C}^{\nu}\triangleq(c^{\nu}_{ij})_{i,j=1}^{m}$ satisfies the assumption below.

Moreover, $\mathbf{C}^{\nu}$ is column stochastic, i.e., $\mathbf{1}^{\top}\mathbf{C}^{\nu}=\mathbf{1}^{\top}$, for all $\nu=0,1,\ldots.$

We conclude this section stating the counterparts of the definitions introduced in Sec. 2, adjusted here to the case of directed time-varying graphs. Using the column stochasticity of $\mathbf{C}^{\nu}$ and (73d), one can see that opposed to (18), the average gradient is now preserved on the weighted average of the $\mathbf{y}_{i}$’s:

where $\overline{\nabla\mathbf{f}}^{\nu}$ is defined in (17). This suggests to decompose $\mathbf{y}^{\nu}$ into its weighted average and the consensus error, defined respectively as

Accordingly, we define the weighted average of $\mathbf{x}^{\nu}$ and the consensus error as

In addition, we also generalize the definition of the optimality gap as

Furthermore, we will use the following lower and upper bounds of $\phi_{i}^{\nu}$ [35, Prop. 1]

​ The proof of linear convergence of SONATA (Algorithm 3) follows the same path of the one developed in Sec. ​​3.3 for the case of undirected graphs. Hence, we omit similar derivations and highlight only the key differences. We will tacitly assume that Assumptions A, B$\,{}^{\prime}$, C, and E are satisfied.

The optimality gap sequence $\{p_{\boldsymbol{\phi}}^{\nu}\}$ satisfies:

where the constants $L_{\operatorname{mx}}$ and $\widetilde{\mu}_{\operatorname{mn}}$ are defined in (4) and (24), respectively; and $\sigma(\alpha)\in(0,1)$ and $\eta(\alpha)>0$ are defined in (41).

The proof follows closely that of Proposition 3.4 and thus is omitted. For completeness, we report it in the supporting materials. Here, we only notice that, instead of (27), we built on: \sum_{i=1}^{m}\phi_{i}^{\nu+1}U(\mathbf{x}_{i}^{\nu+1})\leq\sum_{i=1}^{m}\phi_{i}^{\nu}U\big{(}\mathbf{x}_{i}^{\nu+\frac{1}{2}}\big{)}, where we used $\sum_{j=1}^{m}{c^{\nu}_{ij}\phi_{j}^{\nu}}/{\phi_{i}^{\nu+1}}=1$, for all $i\in[m]$.∎

The following bounds hold for $\|\mathbf{x}_{\boldsymbol{\phi},\bot}^{\nu}\|$ and $\|\mathbf{y}_{\boldsymbol{\phi},\bot}^{\nu}\|$:

where $B$ and $\rho_{B}$ are defined in (78), and $\epsilon_{x}$ and $\epsilon_{y}$ are arbitrary positive constants (to be determined).

Using the result in [26, Lemma 5] and [35, Lemma 3, 11], we obtain

The rest of the proof follows similar steps as [52, Lemma 2], hence it is omitted. ∎

where $L_{\operatorname{mx}}$, $\widetilde{L}_{\operatorname{mx}}$, $\widetilde{\mu}_{\operatorname{mn}}$, and $D_{\operatorname{mx}}$ are defined in (4) and (24), respectively.

The proof follows similar path of that of Proposition 3.6 and thus is omitted. ∎

2 Establishing linear rate

Let $P_{\boldsymbol{\phi}}^{K}(z)$, $D^{K}(z)$, $X_{\boldsymbol{\phi},\bot}^{K}(z)$, and $Y_{\boldsymbol{\phi},\bot}^{K}(z)$ denote the transformation (49) of the sequences $\{p_{\boldsymbol{\phi}}^{\nu}\}$, $\{\|\mathbf{d}^{\nu}\|^{2}\}$, $\{\|\mathbf{x}_{\boldsymbol{\phi},\bot}^{\nu}\|^{2}\}$ and $\{\|\mathbf{y}_{\boldsymbol{\phi},\bot}^{\nu}\|^{2}$ $\}$. Given the constants $\sigma(\alpha)$ and $\eta(\alpha)$, defined in Proposition 4.1, and the free parameters $\epsilon_{x},\epsilon_{y}>0$, the following holds:

The proof of the first two inequalities (85) and (88) follows the same steps of those used to prove Proposition 3.8. Applying [43, Lemma 21] to (81a) and (81b) respectively gives (86) and (87).

Chaining the inequalities in Proposition 4.4 as done in for (50) (cf. Fig. 1), we can bound $D^{K}(z)$ as

where $\mathcal{P}(\alpha,z)$ is defined as

and $\mathcal{R}(\alpha,z)$ is a bounded remainder term.

Comparing (95) to (55) we can see that they share the same form and only differ in coefficients. Therefore, with the same argument as in the proof of Theorem 3.9 we can easily arrive at the following conclusion.

We provide the proof in the supporting material. ∎

For sake of completeness, we provide an explicit expression of the linear rates in terms of the step-size $\alpha$ in the supporting material–see Theorem III.1. Table 4 summarizes the expression of the rates achieved by SONATA using the surrogate functions (62) and (63)–a formal statement of these results along with the proofs can be found in the supporting material-see Corollaries IV.1, V.1 and V.2.

The rate estimates in Table 4 are almost identical to those obtained in Sec. 3.4.2, with the difference that the network dependence now is expressed throughout $\rho_{B}$ rather than $\rho$. Therefore, similar comments–as those stated in Sec. 3.4.2–apply to the rates in Table 4. For example, if the network is sufficiently connected ($\rho_{B}$ “small”), its impact on the rate becomes negligible and SONATA matches the network-independent rate achieved on star-topology (cf. Corollary 3.11) or centralized settings. Specifically, when linearization surrogate (62) is used, this rate coincides with the rates of centralized proximal gradient algorithm.

Numerical Results

In this section, we corroborate numerically the complexity results proved in Corollaries 3.12–3.14. As a test problem, we consider the distributed ridge regression:

where the loss function of agent $i$ is $f_{i}(\mathbf{x})=\frac{1}{2n}\|\mathbf{A}_{i}\mathbf{x}-\mathbf{b}_{i}\|^{2}+\lambda\|\mathbf{x}\|^{2}$ [agent $i$ owns data $(\mathbf{A}_{i},\mathbf{b}_{i})$]. Problem parameters are generated as follows. Each row of the measurement matrix $\mathbf{A}_{i}$ is independently and identically drawn from distribution $\mathcal{N}(\mathbf{0},\boldsymbol{\Sigma})$; and $\mathbf{b}_{i}$ is generated according to the linear model $\mathbf{b}_{i}=\mathbf{A}_{i}\mathbf{x}^{*}+\mathbf{n}_{i}$, where $\mathbf{x}^{*}$ is the ground truth, generated according to $\mathcal{N}(5\cdot\mathbf{1},\mathbf{I})$, and $\mathbf{n}_{i}\sim\mathcal{N}(\mathbf{0},0.1\cdot\mathbf{I})$ is the measurement noise. The covariance matrix $\boldsymbol{\Sigma}$ is constructed according to the eigenvalue decomposition $\boldsymbol{\Sigma}=\sum_{j=1}^{d}\lambda_{j}\mathbf{u}_{j}\mathbf{u}_{j}^{\top}$, where the eigenvalues $\{\lambda_{j}\}_{j=1}^{d}$ are uniformly distributed in $[\mu_{0},L_{0}]$. The eigenvectors, forming $\mathbf{U}=[\mathbf{u}_{1},\ldots,\mathbf{u}_{d}]$, are obtained via the QR decomposition of a random $d\times d$ matrix with standard Gaussian i.i.d. elements. The network is generated using an Erdős-Rényi model $G(m,p)$, with $m=30$ nodes and each edge independently included in the graph with probability $p=0.5$.

To investigate the impact of $\kappa_{g}$ and $\beta$ on the convergence rate, we specifically consider the following two scenarios:

We run SONATA using surrogates (62) (linearization) and (63) (local $f_{i}$)–we term it as SONATA-L and SONATA-F, respectively. The simulations parameters of the different experiments are summarized in Table 5; and the algorithmic parameters are set according to Corollaries 3.12–3.14. The expressions are not tight in terms of the absolute constants. To show convergence rate in both Cases I and II in Corollary 3.12-3.14, we enlarged the second term in the expression of $\alpha_{\operatorname{mx}}$ by a constant factor. We measure the algorithm’s complexity using $T_{\epsilon}=\inf\left\{\nu\geq 0\,|\,\frac{1}{m}\sum_{i=1}^{m}(F(\mathbf{x}_{i}^{\nu})-F^{\star})\right.$ $\left.\leq 10^{-7}\right\}$.

In Table 6, we report the corresponding iteration complexity of SONATA for each simulation setup (s.1)-(s.6) in Table 5. Each figure is generated under one particular realization of the problem setting. Further, in order to compare the complexity of SONATA across different settings, all the simulations share the same network parameters, as well as the same data set whenever the problem parameters are the same. The results of our experiments are reported in Table 6; the curve are generated using only one random realization for visualization clarity. However, the behavior of the curves (e.g., scalability with respect to the parameters) is representative and consistent across all the random experiments we conducted.

$\bullet$ Scalability with respect to $\kappa_{g}$. Consider setting (S.I) wherein $\beta$ is fixed and $\lambda$ is changing. Figures for (s.1)-(s.3) show that when $\alpha=1$ (blue curve), the iteration complexity of SONATA-L scales linearly with respect to $\kappa_{g}$ [as predicted by Corollary 3.12], while that of SONATA-F is invariant whenever $\beta<\mu$ [as stated in Corollary 3.13]. When $\beta\geq\mu$, the iteration complexity of SONATA-F grows as $\lambda$ increases since $\beta/\mu$ decreases [cf. Corollary 3.14]. However, the increasing rate is much slower than SONATA-L, due to the fact that $(\beta/\mu)/\kappa_{g}=\beta/L\ll 1$ for large $\lambda$. When $\alpha<1$, the iteration complexity scales quadratically with respect to $\kappa_{g}$, in all settings, as predicted by our theory.

$\bullet$ Scalability with respect to $\beta$. Consider now setting (S.II), where we decrease the local sample size $n$ to increase $\beta$. In contrast to setting (S.I), Figures for (s.4) and (s.5) show that, with $\alpha=1$, the iteration complexity of SONATA-F scales linearly with $\beta/\mu$ when $\beta>\mu$, while that of SONATA-L is invariant–this is consistent with Corollaries 3.12 and 3.14. When $\alpha<1$, the iteration complexity scales quadratically with respect to $\beta/\mu$. Finally, the plot associated with (s.6) simply reveals that when $\beta<\mu$, iteration complexity of SONATA-F remains bounded, as stated in Corollary 3.13.

Appendix A Proof of (54)

Chaining the inequalities in (50) as shown in Fig. 1, we have

Notice that, under (51), $G_{P}(\alpha,z)$, $G_{X}(z)$, $G_{Y}(z)$, and $\omega_{p}$, $\omega_{x}$, $\omega_{y}$ are all bounded, which implies that the reminder $\mathcal{R}(\alpha,z)$ in (50) is bounded as well.$\hfill\square$

Appendix B Proof of Theorem 3.10

We find the smallest $z$ satisfying (51) such that $\mathcal{P}(\alpha,z)<1$, for $\alpha\in(0,\alpha_{\operatorname{mx}})$, with $\alpha_{\operatorname{mx}}\in(0,1)$ to be determined.

Let us begin considering the condition $z>\sigma(\alpha)$ in (51). To simplify the analysis, we impose instead the following stronger version

Observe that in the expression of $\mathcal{P}(\alpha,z)$, the only coefficient multiplying $\alpha^{2}$ that depends on $\alpha$ is the optimization gain $G_{P}(\alpha,z)\triangleq{\eta(\alpha)}/({z-\sigma(\alpha)}).$ Using (97), $G_{P}(\alpha,z)$ can be upper bounded as

where $A_{1,\theta}$, $A_{2,\theta}$ and $A_{3,\theta}$ are constants defined as

Condition (100) shows the rate $z$ must satisfy

Notice that, under $\epsilon_{x}=\epsilon_{y}=(\sqrt{z}-\rho)/\rho$, (101) implies $z>\rho^{2}(1+\epsilon_{x})=\rho^{2}(1+\epsilon_{y})=\rho\sqrt{z}$, which are the other two conditions on $z$ in (51). Therefore, overall, $z$ must satisfy (97) and (101). Letting $\epsilon_{opt}=\epsilon_{opt}^{\star}$ in (97), the condition simplifies to

Therefore, the overall convergence rate can be upper bounded by $\mathcal{O}(\bar{z}^{\nu})$, where

Note that as $\alpha$ increases from , the first term in the max operator above is monotonically increasing from $\rho^{2}<1$ while the second term is monotonically decreasing from $1$. Therefore, there must exist some $\alpha^{*}$ so that the two terms are equal, which is

To conclude, given the step-size satisfying $\alpha\in(0,\alpha_{\operatorname{mx}})$, the sequence $\{\|\mathbf{d}^{\nu}\|^{2}\}$ converges at rate $\mathcal{O}(z^{\nu})$, with $z$ given in (61). $\hfill\square$

Appendix C Proof of Corollary 3.11

Since $\mathbf{W}=\mathbf{J}$, we have $\boldsymbol{\delta}^{\nu}=\mathbf{0}$; then (33a) and (35) reduce to

respectively. Combining (106) and (107) and using $\alpha<2\widetilde{\mu}_{\operatorname{mn}}/(\widetilde{\mu}_{\operatorname{mn}}-D_{\operatorname{mn}})$, yield

We customize next (64) to the specific choices of the surrogate functions.

Finally, setting $\alpha=\min\{1,2\widetilde{\mu}_{\operatorname{mn}}/((\mu-\beta)_{+}+\beta)\}=1$ in the expression above, yields (65). $\square$

Appendix D Proof of Corollary 3.12

According to Theorem 3.10, the rate $z$ can be bounded as

where $J$ and $A_{\frac{1}{2}}$ are defined in (103) and (101), respectively.

Accordingly, the expressions of $J$ and $A_{\frac{1}{2}}$ read:

where in the last inequality we have used the fact that $\kappa_{g}\geq 1$.

Using the above expressions, in the sequel we upperbound $z_{1}$ and $z_{2}$.

Since $\alpha\in(0,1]$ must be chosen so that $z\in(0,1]$, we impose $\max\{z_{1},\bar{z}_{2}\}<1$, implying $\alpha\leq\min\{J^{-1},(1-\rho)^{2}/(M\rho),1\}$. Since $J^{-1}>1$ [cf. (111)], the condition on $\alpha$ reduces to $\alpha\leq\alpha_{\operatorname{mx}}\triangleq\min\{(1-\rho)^{2}/(M\rho),1\}$. Choose $\alpha=c\cdot\alpha_{\operatorname{mx}}$, for some given $c\in(0,1)$. Depending on the value of $\rho$, either $\alpha_{\operatorname{mx}}=1$ or $\alpha_{\operatorname{mx}}=(1-\rho)^{2}/(M\rho)$.

$\bullet$ Case I: $\alpha_{\operatorname{mx}}=1$. This corresponds to the case $M\rho\leq(1-\rho)^{2}$, which happens when the network is sufficiently connected ($\rho$ is small). Note that, we also have $\rho\leq 1/110$, otherwise $M\rho\geq 110\,\kappa_{g}\,\rho>1>(1-\rho)^{2}$. In this setting, $\alpha=c\cdot\alpha_{\operatorname{mx}}=c$, and