Acceleration in Distributed Optimization Under Similarity

Ye Tian∗ Gesualdo Scutari∗ Tianyu Cao∗ Alexander Gasnikov†

∗Purdue University † Moscow Institute of Physics and Technology

###### Abstract

We study distributed (strongly convex) optimization problems over a network of
agents, with no centralized nodes. The loss functions of the agents are
assumed to be similar, due to statistical data similarity or otherwise. In
order to reduce the number of communications to reach a solution accuracy, we
proposed a preconditioned, accelerated distributed method. An
$\varepsilon$-solution is achieved in
$\tilde{\mathcal{O}}\big{(}\sqrt{\frac{\beta/\mu}{(1-\rho)}}\log
1/\varepsilon\big{)}$ number of communications steps, where $\beta/\mu$ is the
relative condition number between the global and local loss functions, and
$\rho$ characterizes the connectivity of the network. This rate matches (up to
poly-log factors) for the first time lower complexity communication bounds of
distributed gossip-algorithms applied to the class of problems of interest.
Numerical results show significant communication savings with respect to
existing accelerated distributed schemes, especially when solving ill-
conditioned problems.

##  1 Introduction

We study distributed optimization over a network of $m$ agents, in the form

$\\!\\!\min_{x\in\mathbb{R}^{d}}\,\,u(x)\triangleq f(x)+r(x),\quad
f(x)\triangleq\frac{1}{m}\sum_{i=1}^{m}f_{i}(x),$ (P)

where $f_{i}:\mathbb{R}^{d}\to\mathbb{R}$ is the loss function of agent $i$,
known only to that agent; and $r:\mathbb{R}^{d}\to[-\infty,\infty]$ is an
extended-value function (known to all agents), which is instrumental to
enforce further conditions on the solution, such as sparsity or constraints.
The network of agents is modelled as undirected, fixed graph, with no
centralized node; we refer to such architectures as mesh networks.

An instance of (P) 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.
Specifically, denoting by
$\mathcal{D}_{i}=\\{z_{i}^{(1)},\ldots,z_{i}^{(n)}\\}$ the set of $n$ samples
owned by agent $i$, the empirical risk $f_{i}$ in (P) reads

$f_{i}(x)=\frac{1}{n}\sum_{j=1}^{n}\ell\big{(}x;z_{i}^{j}\big{)},\vspace{-0.2cm}$
(1)

where $\ell(x;z_{i}^{j})$ measures the mismatch between the parameter $x$ and
the sample $z_{i}^{j}$.

The lack of global knowledge of $f$ from the agents and of a centralized node
in the network calls for the design of distributed algorithms, whereby agents
alternate a computation procedure based on local information and communication
round(s) with neighboring nodes. Since the cost of communications is often the
bottleneck in distributed computing, if compared with local (parallel)
computations (e.g., [Bekkerman et al., 2011, Lian et al., 2017]), a lot of
research has been devoted to designing distributed algorithms that are
_communication efficient_.

Communication-saving via acceleration: Acceleration (in the sense of Nesterov)
has been extensively investigated as a procedure to reduce the communication
burden of distributed algorithms–Table 1 summarizes the communication
complexity of existing first-order methods over mesh networks (see Sec. 1.2
for a discussion of these works). For $L$-smooth and $\mu$-strongly convex
functions $f$ in (P), linear convergence rate is certified, with a constant
factor scaling with $\sqrt{\kappa}$ ($\kappa\triangleq L/\mu$ is the condition
number of $f$). This dependence is not improvable, meaning that it matches
lower communication complexity bounds for the class of $L$-smooth and
$\mu$-strongly convex function $f_{i}$’s [Scaman et al., 2017]. However, for
ill-conditioned problems–e.g., the typical setting of many ERM problems
wherein the optimal regularization parameter for test predictive performance
is very small–$\kappa$ can be extremely large; hence the aforementioned
scaling of the number of communications with $\sqrt{\kappa}$ is
unsatisfactory.

Exploiting function similarity: Further improvements can be obtained if extra
structure is postulated for the $f_{i}$’s. This is, e.g., the case of ERM
problems wherein each $f_{i}$ [see (1)] has an additional finite-sum
structure. This is an instance of the property known as function similarity
(see, e.g., [Shamir et al., 2014, Arjevani and Shamir, 2015, Zhang and Lin,
2015, Hendrikx et al., 2020b]):

$\left\|\nabla^{2}f_{i}(x)-\nabla^{2}F(x)\right\|\leq{\beta},\vspace{-0.1cm}$
(2)

for all $x$ in a proper domain of interest and all $i=1,\ldots,m$, where
$\beta>0$ measures the degree of similarity between the Hessian matrices of
local and global losses. For instance, in the aforementioned ERM setting, when
data are i.i.d. among machines, the $f_{i}$’s in (1) reflect statistical
similarities in the data residing at different nodes, resulting in
$\beta=\tilde{O}(1/\sqrt{n})$ w.h.p., where $n$ is the local sample size
($\tilde{O}$ hides log-factors and dependence on $d$).

The situation $1+\beta/\mu\ll\kappa$, happens in several scenarios. For
instance, consider some ill-conditioned functions. Another example are ERMs
with optimal regularization $\mu=\mathcal{O}(1/\sqrt{mn})$ and
$L=\mathcal{O}(1)$ (e.g., see [Zhang and Lin, 2015, Table 1] for ridge
regression), we have: $\kappa=\mathcal{O}(\sqrt{m\cdot n})$ while
$\beta/\mu=\mathcal{O}(\sqrt{m})$–the former grows with the local sample size
$n$, while the latter is independent. This motivated a surge of studies aiming
at exploiting function similarity coupled with acceleration to boost
communication efficiency (see Sec. 1.2): linear convergence is certified with
a number of communication steps scaling (asymptotically [Hendrikx et al.,
2020b]) with $1+\sqrt{\beta/\mu}$ [Zhang and Lin, 2015], which can be
significantly smaller than $\sqrt{\kappa}$, and matches lower complexity
bounds [Arjevani and Shamir, 2015] up to log-factors. These algorithms however
are centralized and cannot be implemented over mesh networks. This suggests
the following open question:

Is linear convergence with a number of communications scaling with
$1+\sqrt{\beta/\mu}$ achievable by any distributed algorithm over mesh
networks?

Table 1: Distributed Accelerated Algorithms over Mesh Networks: SSDA/MSDA
[Scaman et al., 2017], OPAPC [Kovalev et al., 2020], Accelerated Dual Ascent
[Uribe et al., 2020, Alg. 3], APM-C [Li et al., 2018], Mudag [Ye et al.,
2020a], Accelerated EXTRA [Li and Lin, 2020], DAccGD [Rogozin et al., 2020],
and DPAG [Ye et al., 2020b]. $L$ (resp. $\mu$) denotes the smoothness (resp.
strong convexity) constant of $F$ while $L_{\textnormal{m}x}$ (resp.
$\mu_{\textnormal{m}n}$) is the largest smoothness (smallest strong convexity)
constant of the $f_{i}$’s; $\rho$ is the network connectivity;
$\widetilde{\mathcal{O}}$ hides poly-log factors.

Algorithm | Problem | Similarity | Gossip Matrix | Rate (# comm.)
---|---|---|---|---
SSDA/MSDA, OPAPC, Accelerated Dual Ascent |  $f_{i}$ scv, $r=0$ | ✗ |  |  $\mathcal{O}\left(\sqrt{\frac{L_{\textnormal{m}x}}{\mu_{\textnormal{m}n}}}\sqrt{\frac{1}{1-\rho}}\log(\frac{1}{\varepsilon})\right)$
APM-C |  $f_{i}$ scv, $r=0$ | ✗ | PSD |  $\mathcal{O}\left(\sqrt{\frac{L_{\textnormal{m}x}}{\mu_{\textnormal{m}n}}}\sqrt{\frac{1}{1-\rho}}\log^{2}(\frac{1}{\varepsilon})\right)$
Mudag |  $f$ scv, $r=0$ | ✗ | PSD |  $\widetilde{O}\left(\sqrt{\frac{L}{\mu}}\sqrt{\frac{1}{1-\rho}}\log(\frac{1}{\varepsilon})\right)$
Accelerated EXTRA |  $f_{i}$ scv, $r=0$ | ✗ |  |  $\widetilde{O}\left(\sqrt{\frac{L_{mx}}{\mu_{mn}}}\sqrt{\frac{1}{1-\rho}}\log(\frac{1}{\varepsilon})\right)$
DAccGD |  $f_{i}$ scv, $r=0$ | ✗ |  |  $\widetilde{O}\left(\sqrt{\frac{L}{\mu}}{\frac{1}{1-\rho}}\log^{2}(\frac{1}{\varepsilon})\right)$
DPAG |  $f$ scv, $r\neq 0$ | ✗ | PSD |  $\widetilde{O}\left(\sqrt{\frac{L}{\mu}}\sqrt{\frac{1}{1-\rho}}\log(\frac{1}{\varepsilon})\right)$
ACC-SONATA (this work) |  $f$ scv, $r\neq 0$ | ✓ |  | $\widetilde{O}\left(\sqrt{\frac{\beta}{\mu}}\sqrt{\frac{1}{1-\rho}}\log(\frac{1}{\varepsilon})\right)$

###  1.1 Contributions

We provide a positive answer to the above question.
(i) Algorithm design: We proposed Accelerated-SONATA (ACC-SONATA), an inexact
accelerated proximal method (outer loop) for (P), embedded with the
distributed algorithm known as SONATA [Sun et al., 2019] (inner loop), which
approximatively solves the proximal subproblem over mesh networks, according
to a properly designed notion of inexactness. SONATA couples local
preconditioning via surrogation of $f_{i}$ with a gradient tracking mechanism,
estimating locally the gradient of the global loss $f$. At high level, the
outer loop ensures acceleration while SONATA exploits function similarity to
boost the convergence rate of the inner loop. A direct acceleration of the
mirror method, achieving $\widetilde{\mathcal{O}}(\beta/\mu)$ over star-
networks [Lu et al., 2020], does not seem possible in general [Dragomir et
al., 2019].

(ii) New analysis: While ACC-SONATA is conceptually inspired by the proximal
acceleration framework for centralized optimization [d’Aspremont et al.,
2021], such as Catalyst [Hongzhou et al., 2015], formally it is not an
instance of any of existing accelerations. It is also different from
distributed algorithms accelerated á la Catalyst [Li and Lin, 2020, Hendrikx
et al., 2020a], which do not exploit any function similarity, use a gradient-
tracking mechanism, or deal with composite optimization ($r\neq 0$). Hence, a
new convergence analysis is needed for ACC-SONATA, which represents the
technical contribution of this work. We hinge on new potential functions (for
the inner and outer loop) that incorporate consensus errors on the local
optimization, extrapolation variables, and gradient-tracking variables. Such
potentials also shed light on the appropriate warm-restart strategy and
termination criterion of the inner-loop SONATA, which are implementable in a
distributed setting (this is not the case if using criteria in [Hongzhou et
al., 2018], in particular when $r\neq 0$). We remark that the proposed
analysis, although developed for ACC-SONATA, is fairly general and potentially
applicable to a variety of other distributed algorithms, replacing SONATA in
the inner loop. This will be the subject of future investigation.

(iii) Guarantees: By a proper choice of local surrogations (mirror-prox-like),
ACC-SONATA provably achieves an $\varepsilon$-solution (on the objective
value) of (P) in
$\tilde{\mathcal{O}}\big{(}\sqrt{\frac{\beta/\mu}{1-\rho}}\log
1/\varepsilon\big{)}$ number of communications steps, where $\rho$
characterizes the connectivity of the network. This matches for the first time
lower communication complexity bounds, up to poly-log factors. Quite
interestingly, when the agents’ losses $f_{i}$’s are heterogeneous functions,
a different choice of surrogate (linearization of each $f_{i}$) guarantees
linear convergence with communication scaling of $\sqrt{\kappa}$–this compares
favorable with existing first-order accelerated methods (see Table 1), with
ACC-SONATA being applicable on general (P) with $r\neq 0$ and with no
restrictions on gossip matrices. Our numerical results on synthetic and real
data supports the competitive performance of our method over the state of the
art.

###  1.2 Related works

Acceleration under function similarity (2) has been extensively investigated
over master/workers architectures while it remains unexplored over mesh
networks, as inferred by the following discussion.

$\bullet$ Master/workers systems: Several papers explored the idea of
statistical preconditioning to decrease the communication complexity of
solving (P) over such architectures; examples follow. DANE [Shamir et al.,
2014] is a mirror-descent type algorithm for $\eqref{eq:problem}$ with $r=0$
whereby workers perform a local data preconditioning via a suitably chosen
Bregman divergence, and the master averages the solutions of the workers. For
quadratic losses, DANE achieves communication complexity
$\widetilde{\mathcal{O}}((\beta/\mu)^{2}\log 1/\varepsilon)$. More recently,
[Fan et al., 2019] proposed CEASE, which achieves DANE’s complexity but for
nonquadratic losses and $r\neq 0$. Applying the convergence analysis of mirror
descent in [Lu et al., 2020] to CEASE enhances its rate to
$\widetilde{\mathcal{O}}((\beta/\mu)\log 1/\varepsilon)$.

Further improvements are achievable employing acceleration. Efforts include:
DiSCO [Zhang and Lin, 2015], an inexact damped Newton method coupled with a
preconditioned conjugate gradient (to compute the Newton direction), which
achieves communication complexity
$\widetilde{\mathcal{O}}((1+\sqrt{\beta/\mu})\log 1/\varepsilon)$ for self-
concordant losses (and $r\neq 0$); AIDE [Reddi et al., 2016], which uses the
Catalyst framework [Hongzhou et al., 2015], matching the rate of DiSCO for
quadratic losses (and $r=0$); DANE-HB [Yuan and Li, 2020], a variant of DANE
equipped with Heavy Ball momentum and matching for quadratic functions the
communication complexity of DiSCO and AIDE; and SPAG[Hendrikx et al., 2020b],
a preconditioned direct accelerated method, achieving for nonquadradic losses
asymptotically the convergence rate
${\mathcal{O}}((1-1/\sqrt{\beta/\mu})^{k})$ ($k$ is the iteration index).

None of above methods are implementable over mesh networks, because they all
rely on the presence of a master node. Notice also that, although designed for
mesh networks, our proposed method, ACC-SONATA compares favorably also with
the aforementioned schemes (specifically designed for star networks) by
achieving communication complexity
$\mathcal{O}(\sqrt{\beta/\mu}\log(\beta/\mu)\log 1/\varepsilon)$ for
nonquadratic losses.

$\bullet$ Acceleration over mesh networks: Given the focus of this work, we
comment next only distributed algorithms over mesh networks employing some
form of acceleration and provably convergent–they are summarized in Table 1.
Although substantially different–some are primal [Ye et al., 2020a, Ye et al.,
2020b, Li and Lin, 2020, Kovalev et al., 2020, Rogozin et al., 2020] others
are dual or penalty-based [Scaman et al., 2017, Uribe et al., 2020, Li et al.,
2018] methods, and applicable to special instances of (P) (mainly with $r=0$)
and subject to special design constraints (e.g., positive semidefinite gossip
matrix)–they all achieve linear convergence rate, with communication
complexity scaling some with $\sqrt{\kappa_{\ell}}$
($\kappa_{\ell}=L_{\textnormal{m}x}/\mu_{\textnormal{m}n}$ is the “local”
condition number) and others with $\sqrt{\kappa}$ ($\kappa=L/\mu$ is the
condition number of $f$). Note that in general $\kappa\ll\kappa_{\ell}$; hence
the latter group is preferable to the former. By using only gradient
information of the local $f_{i}$’s, none of these methods can take advantage
of function similarity. This means that their rates still scale as
$\sqrt{\kappa}$ no matter how small $\beta$ is (even $\beta=0$, i.e., all
$f_{i}$’s identical), which is highly sub-optimal when considering, e.g., ill-
conditioned losses, and contrasts with the lower complexity bound
$\mathcal{O}(\sqrt{\beta/\mu})$ (cf. Sec. 2.1).

To our knowledge, Network-DANE [Li et al., 2020] and SONATA [Sun et al., 2019]
are the only two methods that leverage statistical similarity to enhance
convergence over mesh networks, with the latter achieving communication
complexity scaling with $\widetilde{\mathcal{O}}(\beta/\mu)$ for nonquadratic
losses and $r\neq 0$. These methods however are not accelerated, missing thus
the more favorable scaling $\sqrt{\beta/\mu}$. The proposed accelerated
method, ACC-SONATA, fills exactly this gap.

##  2 Setup and Background

We study (P) under the following assumptions.

######  Assumption 1.

Given problem (P),

  1. (i)

$r:\mathbb{R}^{d}\to\mathbb{R}$ is a proper, closed, convex function; let
$\textnormal{dom}\,r$ denote its domain;

  2. (ii)

Each $f_{i}$ is convex and twice differentiable over (an open set containing)
$\textnormal{dom}g$;

  3. (iii)

$f$ is $\mu$-strongly convex and $L$-smooth on $\textnormal{dom}g$, with
$0<\mu\leq L<\infty$. Define $\kappa=L/\mu$.

Note that (i)-(ii) implies that each $f_{i}$ is $\mu_{i}$-strongly convex and
$L_{i}$-smooth, with $0\leq\mu_{i}\leq L_{i}<\infty$ (not each $f_{i}$ need to
be strongly convex); let $L_{\max}=\max_{i\in[m]}L_{i},$, were we denoted
$[m]\triangleq\\{1,\ldots,m\\}.$

Function similarity is captured by the following.

######  Definition 2.

The functions $f_{i}$’s are $\beta$-similar:

$\max_{i\in[m]}\,\,\sup_{x\in\textnormal{dom}g}\left\lVert\nabla^{2}f_{i}(x)-\nabla^{2}f(x)\right\rVert\leq\beta.\vspace{-0.2cm}$

Note that Assumption 2 always holds with

$\beta\leq\max_{i\in[m]}\max\big{\\{}|L_{i}-\mu_{i}|,|\mu-
L_{i}|\big{\\}}.\vspace{-0.1cm}$

The case of interest is of course when $1+\beta/\mu<L/\mu$, which is the
typical situation of ill-conditioned $f$’s.

Network model: The network of agents is an undirected, graph
$\mathcal{G}=(\mathcal{E},\mathcal{V})$, where $\mathcal{V}=[m]$ denotes the
vertex set (the set of agents) while $\mathcal{E}$ is the set of edges;
$\\{i,j\\}\in\mathcal{E}$ if there is a communication link between agent $i$
and agent $j$. For the sake of notation, we assume $\\{i,j\\}\in\mathcal{E}$
for any $i\in[m].$ We make the blanket assumption that the graph $\mathcal{G}$
is connected.

The distributed algorithms of interest employ gossip communications–each node
averages the values of its neighbors’ variables. The weights of this averaging
process (collected into a matrix $W\in\mathbb{R}^{m\times m}$) satisfy the
following standard assumptions.

######  Assumption 3.

The matrix ${W}\in\mathbb{R}^{m\times m}$ belongs to the class
$W=P_{M}(\overline{{W}})$, for some $M\in\mathbb{N}_{++}$ and
$\overline{{W}}\in\mathcal{R}^{m\times m}_{+}$, where $P_{M}$ is a polynomial
with degree at most $M$ with $P_{M}(1)=1$, and $\overline{{W}}$ satisfies the
following conditions:

  * (i)

$\overline{W}$ is compliant with $\mathcal{G}$, that is its $(i,j)$-entries
$\bar{w}_{ij}$ satisfies: $\bar{w}_{ii}>0$, for all $i\in[m]$;
$\bar{w}_{ij}>0$, if $(i,j)\in\mathcal{E}$; and $\bar{w}_{ij}=0$ otherwise;

  * (ii)

${1}^{\top}\overline{W}={1}^{\top}$ and $\overline{W}{1}=1$ (doubly
stochasticity).

Define $\rho\triangleq\left\lVert W-11^{\top}/m\right\rVert<1$.

The above class of gossip matrices captures single ($M=1$) and multiple
($M>1$) rounds of communications per optimization step (notice that
$P_{M}(1)=1$ is to ensure the doubly stochasticity of $W$ when $\overline{W}$
is so). Several rules have been proposed in the literature fulfilling
Assumption 3, including the Laplacian, the Metropolis-Hasting, and the
maximum-degree weights rules as well as Chebyshev [Wien, 2011, Scaman et al.,
2017] or Jacobi [Berthier et al., 2020] polynomials-based accelerations.

###  2.1 Lower complexity bounds over mesh networks under similarity

We informally state here lower communication complexity bounds over mesh
networks for the class of problems (P) satisfying Assumptions 1 and 2, and
certain distributed gossip algorithms of interest:

$\Omega\left(\sqrt{\frac{\beta/\mu}{1-\rho}}\log\left(\frac{\mu\left\lVert
x^{\star}\right\rVert^{2}}{\varepsilon}\right)\right).$ (3)

The formal statement of this result can be found in the supplementary material
(cf. Theorem 7, Sec. A). The next section is devoted to the design of the
first distributed algorithm matching such a lower complexity bound (up to
poly-log factors). As anticipated, our scheme hinges on the SONATA algorithm
[Sun et al., 2019], which we recall next.

###  2.2 A building block: SONATA algorithm

The instance of SONATA used in this work is summarized in Algorithm 3 (assumed
to be applied to (P)). Each agent $i$ owns local copies $x_{i}$ of the shared
optimization variable $x$ along with the auxiliary variable $y_{i}$ that is a
local proxy of $\nabla f$, which is not available at the agents’ sides. In
parallel and iteratively, agents update their $x$-variables, first solving in
(S.1) a local approximation of (P) wherein $\tilde{f}_{i}(x;x_{i}^{k})$ is a
surrogate of $f_{i}$ at $x_{i}^{k}$ and the linear term $y_{i}^{k}-\nabla
f_{i}(x_{i}^{k})$ is an estimate of $\sum_{j\neq i}\nabla f_{j}(x_{i}^{k})$.
This is followed by a communication step, (S.2), instrumental to enforce
asymptotic consensus among the $x$-variables and track $\nabla f$ via the
$y$-ones.

Several surrogate functions are feasible, see [Sun et al., 2019]. Here, we
focus on the following two:

$\displaystyle\widetilde{f}_{i}(x;x_{i}^{k})=f_{i}(x)+\frac{\beta}{2}\|x-x_{i}^{k}\|^{2};$
(4)
$\displaystyle\widetilde{f}_{i}(x;x_{i}^{k})=f_{i}(x_{i}^{k})+\left\langle\nabla
f_{i}(x_{i}^{k}),x-x_{i}^{k}\right\rangle+\frac{L}{2}\|x-x_{i}^{k}\|^{2}.$ (5)

Note that the use of the linearization (5) corresponds to perform at each
agent’s side a proximal gradient step; this is the typical update of the
majority of existing distributed algorithms (as those in Table 1). Such a
choice does not permit to take advantage of function similarity, when any. In
fact, when (5) is employed and a weight matrix $W$ satisfying Assumption 3 is
used for the consensus and tracking steps, SONATA applied to (P) achieves an
$\varepsilon$-solution (in terms of objective value) in
$\widetilde{\mathcal{O}}\big{(}\kappa\,\frac{1}{1-\rho}\log
1/\varepsilon\big{)}$ number of communications. Communication saving under
function similarity is provably achievable instead using surrogate (4),
resulting in a communication complexity of
$\widetilde{\mathcal{O}}\big{(}\frac{\beta}{\mu}\cdot\frac{1}{1-\rho}\log
1/\varepsilon\big{)}.$ This motivated us to use SONATA as building block of
our accelerated method aiming at exploiting function similarity.

Algorithm 1 SONATA$(\left\\{f_{i}\right\\}_{i\in[m]},\,x^{0},\,y^{0},\,T)$

Input: $\left\\{f_{i}(x)\right\\}_{i\in[m]},\,\,r(x)$ [cf. (P)];

$x^{0}=(x_{i}^{0})_{i\in[m]}$ [init. points],

$y^{0}=(y_{i}^{0})_{i\in[m]}$ [grad.-tracking init.],

$T>0$ [# iterations];
Output:
$x^{T}=\big{(}x_{i}^{T}\big{)}_{i\in[m]},\,\,y^{T}=\big{(}y_{i}^{T}\big{)}_{i\in[m]}$;

for $k=0,1,2,\ldots,T-1$ do

(S.1) Local computations: for all $i\in[m]$,

$\displaystyle{x}_{i}^{k+1/2}=\mathop{\rm argmin}_{x\in\mathbb{R}^{d}}$
$\displaystyle\,\,\tilde{f}_{i}(x;x_{i}^{k})$ $\displaystyle+\left\langle
y_{i}^{k}-\nabla f_{i}(x_{i}^{k}),x-x_{i}^{k}\right\rangle+r(x);$

(S.2) Communications: for all $i\in[m]$,

$\displaystyle x_{i}^{k+1}=\sum_{j=1}^{m}w_{ij}x_{j}^{k+1/2},$ $\displaystyle
y_{j}^{k+1}=\sum_{j=1}^{m}w_{ij}\big{(}y_{j}^{k}+\nabla
f_{j}(x_{j}^{k+1})-\nabla f_{j}(x_{j}^{k})\big{)}.$

end for

##  3 Accelerated SONATA

We are ready to introduce our main algorithm, Algorithm 2. At high level, the
scheme can be interpreted as a successive application of SONATA for the
inexact minimization of the function $u_{k}$ in (6), wherein the $z$-variable
in the quadratic term plays the role of the extrapolation à la Nesterov, to
gain acceleration. The use of SONATA in the inner loop allows us to take
advantage of function similarity, if any, by choosing surrogates as in (4) and
a suitable value for the coefficients $\alpha\in(0,1)$ and $\delta>0$ (see
Theorem 4). Notice the warm restart of SONATA every $T$ iterations, with in
particular the gradient tracking initialization unconventionally chosen, as
recommended by our convergence analysis.

Algorithm 2 Accelerated SONATA

Input: $\beta$, $\mu$, $\delta>0$, $\alpha=\sqrt{\mu/(\mu+\delta)}$;

$x_{i}^{0}=z_{i}^{0}=z_{i}^{-1}=0$, $y_{i}^{0}=\nabla f_{i}(x_{i}^{0})$
Output: $x^{K}=(x_{i}^{K})_{i\in[m]}$
for $k=0,1,2,\ldots,K-1$ do

Set: $f_{i}^{k}(x)=f_{i}(x)+\frac{\delta}{2}\left\lVert
x-z_{i}^{k}\right\rVert^{2},$

$u_{k}(x)=\frac{1}{m}\sum_{i=1}^{m}f_{i}^{k}(x)+r(x);\vspace{-0.3cm}$ (6)

(S.1) Inner loop via SONATA:

$\displaystyle\big{(}x^{k+1},\,y^{k+1}\big{)}=$
$\displaystyle\texttt{SONATA}\Big{(}\left\\{f_{i}^{k}\right\\}_{i\in[m]},\,x^{k},\,y^{k}+\delta\,\left(z^{k-1}-z^{k}\right),\,T\Big{)};$

(S.2) Extrapolation step:

$z_{i}^{k+1}=x_{i}^{k+1}+\frac{1-\alpha}{1+\alpha}\,(x_{i}^{k+1}-x_{i}^{k}),\quad\forall
i\in[m].$

end for

Algorithm rationale: The genesis of the algorithmic design can be traced back
to the idea of acceleration of a centralized inexact proximal method (outer
loop) (see, e.g., [d’Aspremont et al., 2021]), whose proximal subproblems are
approximately solved in a distributed fashion via the SONATA algorithm (inner
loop), satisfying a suitable notion of inexactness (defined in this paper) for
proximal operations. In fact, assuming one can absorb consensus errors on the
agents’ variables $x_{i}$’s and momentum vectors $z_{i}$’s into such a
criterion of solution approximation, we can approximate (S.1) and (S.2) as

$\displaystyle
x_{i}^{k}\approx\overline{x}^{k}\triangleq\frac{1}{m}\sum_{i=1}^{m}x_{i}^{k}\quad\text{and}\quad
z_{i}^{k}\approx\bar{z}^{k}\triangleq\frac{1}{m}\sum_{i=1}^{m}z_{i}^{k},$
$\displaystyle\texttt{(S.1)${}^{\prime}$:}\quad\bar{x}^{k+1}\approx\mathop{\rm
argmin}_{x\in\mathbb{R}^{d}}u(x)+\frac{\delta}{2}\left\lVert
x-\bar{z}^{k}\right\rVert^{2},$
$\displaystyle\texttt{(S.2)${}^{\prime}$:}\quad\bar{z}^{k+1}=\bar{x}^{k+1}+\frac{1-\alpha}{1+\alpha}\,(\bar{x}^{k+1}-\bar{x}^{k}),$

where we used the fact that the minimization of $u_{k}$ in (6) and that of the
function in (S.1)′ have the same solution. The dynamics above are a resemble
of an inexact proximal acceleration [d’Aspremont et al., 2021, Hongzhou et
al., 2015, Hongzhou et al., 2018].

Challenges: Despite the above connection, existing convergence analyses of
centralized accelerated methods break down when applied to ACC-SONATA.
Specifically, (i) the notions of approximate solutions for proximal problems
as in [d’Aspremont et al., 2021, Hongzhou et al., 2018] cannot be satisfied
here, because of the aforementioned consensus errors, let alone their
practical verification in a distributed setting and in the presence of the
nonsmooth function $r$; and (ii) the potential functions used therein are not
adequate to certify linear convergence of the outer loop of ACC-SONATA at the
desired accelerated rate, they do not capture unavoidable consensus and
gradient tracking errors coming out of the inexact, distributed minimization
of $u_{k}$ in (6) via SONATA. Furthermore, the convergence proof of SONATA as
in [Sun et al., 2019] is not directly applicable to study the inner loop, due
to the unconventional restart of the gradient tracking variables. Also,
R-linear convergence of the objective-value gap and consensus/tracking errors
therein seems no longer sufficient to provably obtain acceleration of the
outer loop. Our convergence analysis addresses all these challenges–we refer
to the supplementary material for details.

###  3.1 Convergence analysis

Communication complexity of Acc-SONATA is established in Theorem 4 and Theorem
5 below, pertaining to the use in the inner algorithm SONATA of the surrogates
(4) and (5), respectively. The explicit expression of the constants hidden in
the big-O notation can be found in the supplementary material.

######  Theorem 4.

Consider problem (P) under Assumption 1, with optimal value function
$u^{\star}$ and $\beta>\mu$ w.l.o.g.. Let
$\\{x^{k}\triangleq(x_{i}^{k})_{i\in[m]}\\}$ be the sequence generated by ACC-
SONATA under Assumption 3, with

$\rho\leq\mathcal{O}\left(\left(1+\frac{\kappa-1}{\beta/\mu}\right)^{-2}\right),\vspace{-0.1cm}$
(7)

and the following tuning:

$\delta=\beta-\mu,\quad T=\mathcal{O}\left(\log{\beta}/{\mu}\right),$

and agents’ surrogate functions (4) in SONATA. Define
$\bar{x}^{k}\triangleq\frac{1}{m}\sum_{i=1}^{m}x_{i}^{k}$ and the optimality
gap

$\Delta(x^{k})\triangleq\max\left(\frac{1}{m}\sum_{i=1}^{m}u(x_{i}^{k})-u^{\star},\frac{1}{m}\sum_{i=1}^{m}\left\lVert
x_{i}^{k}-\bar{x}^{k}\right\rVert^{2}\right).$

Then, there holds

$\Delta(x^{k})=\mathcal{O}\left(\left(1-c\,\sqrt{{\mu}/{\beta}}\right)^{k}\right),$
(8)

where $c\in(0,1)$ is some universal constant. Therefore,
$\Delta(x^{K})\leq\varepsilon$, $\varepsilon>0$, in

$\mathcal{O}\left(\sqrt{{\beta}/{\mu}}\cdot
T\cdot\log{1}/{\varepsilon}\right)$ (9)

total (inner plus outer loop) communication steps.

Note that (9) states linear convergence with optimal dependence on
$\beta/\mu$, up to the log-factor $T$. This is achieved under (7), which
requires the network to be sufficiently connected. If the network is not part
of the design, (7) might not be satisfied by the topology under consideration.
Still, (7) can be enforced by running multiple communication rounds (S.2) in
the inner loop SONATA. For example, if Chebyshev polynomial are used for
$W=P_{M}(\overline{W})$, where $\overline{W}$ is a gossip matrix satisfying
Assumption 3 with $\bar{\rho}=\|\overline{W}-11^{\top}/m\|$,
$\mathcal{O}(1/\sqrt{1-\bar{\rho}})$ rounds of communications per iteration
suffices to satisfy (7), resulting in an overall communication complexity

$\mathcal{O}\left(\,\sqrt{\frac{{\beta}/{\mu}}{\
{1-\bar{\rho}}}}\,\log\bigg{(}1+\frac{\kappa-1}{\beta/\mu}\bigg{)}\,T\,\log{1}/{\varepsilon}\right).$

This matches the lower complexity bound in (3) (cf. Theorem 7 in Appendix A)
up to the log-factors.

Consistently with the behaviour of SONATA, if surrogates in the form (5) are
used by the agents, acceleration is still achieved, but with a linear rate
scaling with $\sqrt{\kappa}$ rather than $\sqrt{\beta/\mu}$. This is
formalized next.

######  Theorem 5.

Consider problem (P) under Assumption 1, with optimal value function
$u^{\star}$ and $\kappa>1$ w.l.o.g.. Let
$\\{x^{k}\triangleq(x_{i}^{k})_{i\in[m]}\\}$ be the sequence generated by ACC-
SONATA under Assumption 3, with

$\rho\leq\mathcal{O}\left(\left(2+\frac{\beta/\mu-1}{\kappa}\right)^{-2}\right),\vspace{-0.1cm}$
(10)

and the following tuning:

$\delta=L-\mu,\quad T=\mathcal{O}\left(\log\kappa\right),$

and agents’ surrogate functions (5) in SONATA. Then, there holds

$\Delta(x^{k})=\mathcal{O}\left(\left(1-c\,\frac{1}{\sqrt{\kappa}}\right)^{k}\right),$
(11)

where $c\in(0,1)$ is some universal constant. Therefore,
$\Delta(x^{K})\leq\varepsilon$, $\varepsilon>0$, in

$\mathcal{O}\left(\sqrt{\kappa}\,T\log{1}/{\varepsilon}\right)$ (12)

total (inner plus outer loop) communication steps.

Notice that the more favorable acceleration (9) with respect to (12) and thus
the consequent communication saving, comes in general at the cost of
computationally more demanding subproblems to solve at the agents’ side. This
is common to all the methods (including centralized) exploiting similarity and
discussed in Sec. 1.2, and seems an unavoidable tradeoff. These methods are in
fact designed with the goal of saving communications.

ACC-SONATA on star-networks: Although ACC- SONATA has been designed
specifically for mesh networks, it readily applies to master/workers
architectures; details can be found in the supplementary material; here we
only remark that a direct application of Theorem 4 and Theorem 5 leads to the
following communication complexity to solve (P) over master/workers
architectures

$\mathcal{O}\left(\sqrt{\frac{\beta}{\mu}}\log\bigg{(}\frac{\beta}{\mu}\bigg{)}\log\frac{1}{\varepsilon}\right)\,\,\text{and}\,\,\mathcal{O}\left(\sqrt{\kappa}\log(\kappa)\log\frac{1}{\varepsilon}\right),$

respectively. Quite interestingly, these rates compare favorably with those of
the centralized methods discussed in Sec. 1.2.

Figure 1: Comparison of distributed accelerated algorithms on ridge regression
(synthetic data). (left panel): optimality gap versus total number of
communications, for given $\beta/\mu$ and $\kappa$; (mid panel): number of
communications to reach a precision of $10^{-4}$ versus $\beta/\mu$, for fixed
$\kappa$; (right panel): same quantity versus $\kappa$, for fixed $\beta/\mu$.

Figure 2: Comparison of distributed accelerated algorithms on hinge loss
minimization (COV1 dataset). (left panel): optimality gap versus total number
of communications; (right panel): number of communications to reach a
precision of $10^{-4}$ versus (total) sample.

##  4 Numerical Results

We present some numerical results on synthetic and real data, corroborating
our complexity analysis (Theorems 4 and 5). Additional experiments on
different problem classes and data sets are reported in the supplementary
material.

1) Ridge regression: Our first experiment concerns a ridge regression problem
over a network of agents, modeled as a Erdos-Renyi graph with $m=30$ nodes and
edge probability $p=0.5$. The problem is an instance of (P) with
$f_{i}(x)=1/(2n)\|A_{i}x-b_{i}\|^{2}+\lambda\|x\|^{2}$ [agent $i$ owns data
$(A_{i},b_{i})$] and $r=0$. Data are generated as follows [Sun et al., 2019]:
Each row of ${A}_{i}$ is i.i.d, drawn from $\mathcal{N}({0},{\Sigma})$, where
${\Sigma}=\sum_{j=1}^{d}\lambda_{j}{u_{j}}{u_{j}}^{\top}$. The $\lambda_{j}$’s
are uniformly distributed in [$\mu_{0}$, $L_{0}$], with $\mu_{0}=1$ and
$L_{0}=1000$, and $\mathbf{u}_{1},\ldots,\mathbf{u}_{d}$ are obtained via the
QR decomposition of a $d\times d$ random matrix with standard Gaussian i.i.d
elements. We set ${b}_{i}={A}_{i}{x}^{\star}+{w}_{i}$, where
${x}^{\star}\sim\mathcal{N}(5\cdot{1}_{d},{I})$ is the ground truth and
${w}_{i}\sim\mathcal{N}({0},0.1\cdot{I})$ is the additive noise ($1_{d}$ is
the $d$-dimensional vector of all ones).

Algorithms: We simulated two instances of ACC- SONATA, corresponding to the
choices of the surrogates (4) and (5) in SONATA (inner-loop); we termed them
as ACC-SONATA-F and ACC-SONATA-L, respectively (F stands for “full local
function” while L for “linearization”). The solution of the agents’
subproblems solved in ACC-SONATA-F and ACC-SONATA-L is computed in closed
form. According to Theorems 4 and 5, ACC-SONATA-F is expected to outperform
first-order methods, including ACC-SONATA-L, whenever $1+\beta/\mu<\kappa$,
while ACC-SONATA-L should be competitive otherwise. The free parameters of
these two instances are tuned as suggested by the theory, with
$T=\lceil\log(\beta/\mu)\rceil$ for ACC-SONATA-F and
$T=\lceil\log(\kappa)\rceil$ for ACC-SONATA-L, where $L$ and $\mu$ are
estimated by the data (quadratic function) and so $\beta$ using Definition 2.
The weight matrix $W$ according to the Metropolis-Hasting rule. We compare our
algorithms with the following, widely tested in the literature (Table 1):
APM-C, Mudag, and ACC-EXTRA. The tuning of these schemes follows the
recommendations as in their respective papers. In particular, Mudag and APM-C
require the gossip matrix to be positive definite, so we set $(W+I)/2$, with
$W$ being the matrix used in the other algorithms.

In Fig. 2(left-panel), we fix $\beta/\mu<\kappa$, and plot the optimality gap
$\frac{1}{m}\sum_{i=1}^{m}\|x_{i}^{k}-{x}_{\textnormal{rg}}\|^{2}$ versus the
total number of communications, for each of the algorithms, where
${x}_{\textnormal{rg}}$ is the solution of the ridge regression problem
(computed in closed form). All the schemes achieve linear convergence. As
predicted, ACC-SONATA-F outperforms the other accelerate methods that do not
take advantage of function similarity. Quite interestingly, ACC-SONATA-L
compares quite favorably also with directed acceleration methods, such as
Mudag and AMP-C, while sharing similar computational costs.

To investigate the impact of $\kappa$ and $\beta/\mu$ on the convergence rate
of the algorithms, in the next experiment we consider the following two
scenarios:

(i) Changing $\beta/\mu$ with (almost) fixed $\kappa$: We generate instances
of ridge regression with decreasing $\beta$ and (almost) fixed $\kappa$,
setting $\lambda=0$ and increasing the local sample size $n$ (20 values)
within $[100,40000]$; the empirical $\kappa$ remains approximately constant
close to the nominal value $L_{0}/\mu_{0}=1000$. Fig. 2(mid-panel) captures
this scenario, we plot the number of communications to drive the optimality
gap below $10^{-4}$ versus $\beta/\mu$, and $\kappa\approx 930$.

(ii) Changing $\kappa$ with fixed $\beta/\mu$: We generate a sequence of
instances of ridge regression with varying $\kappa$ (acting on $\lambda$)
while keeping $\beta/\mu$ constant by changing the local sample size $n$ to
compensate for the variations of $\mu$ (due to $\lambda$). Fig. 2(right-panel)
plots the number of communications to drive the optimality gap below $10^{-4}$
versus $\kappa$, for $\beta/\mu\approx 232$.

The following comments are in order. First, the mid-panel confirms what
predicted by Theorem 4: the number of communications of ACC-SONATA-F scales
roughly with $\sqrt{\beta/\mu}$ while that of first-order distributed schemes
is fairly invariant with $\beta/\mu$. This is because methods using only
gradient information (including ACC-SONATA-L) cannot benefit from statistical
similarity. On the other hand, the right-panel shows that communication
complexity of accelerated first-order methods (including ACC-SONATA-L, as
stated in Theorem 5) deteriorates when $\kappa$ grows whereas that of ACC-
SONATA-F is (almost) invariant. It is interesting to remark that ACC-SONATA-F
remains competitive, outperforming the other schemes, even when
$\beta/\mu\approx\kappa$.

2) Hinge loss minimization on real data: We consider the minimization problem
$\min_{x\in\mathbb{R}^{d}}\frac{1}{N}\sum_{i=1}^{N}l(y_{i}\langle
x,z_{i}\rangle)+\frac{\lambda}{2}\left\lVert x\right\rVert^{2}$ over the same
network of agents considered in the previous example, where $\ell$ is the
smooth hinge loss as in [Shamir et al., 2014]. We experimented on the data set
COV1 (see, e.g., [Shalev-Shwartz and Zhang, 2013]). The tuning of the
algorithm is the same a described in the previous experiment, based upon
estimation of the quantities $\mu$, $L$, $L_{\textnormal{mx}}$,
$\mu_{\textnormal{mn}}$, and $\beta$ from the data–see supplementary material
for details. We notice that $\mu\approx\mu_{\textnormal{mn}}\approx\lambda$.
The solution of the agents’ subproblems in ACC-SONATA-F is estimated up to the
accuracy $10^{-10}$ by running the gradient algorithm with step-size equal to
0.03.

Fig. 2 (left-panel) plots the optimality gap
$\frac{1}{m}\sum_{i=1}^{m}\|x_{i}^{k}-{x}_{\textnormal{op}}\|^{2}$ versus the
total number of communications, for each of the algorithms, where
${x}_{\textnormal{op}}$ is an estimate of the solution of the problem
(obtained running the centralized gradient algorithm). The results show that
the proposed methods are competitive also on real data, with ACC-SONATA-F
outperforming the others, when enough samples are present at the agents’
sides. The right-panel plots the number of communications to reach an accuracy
of $10^{-4}$ versus the total sample size (since $m$ is fixed, the local
samples size varies). This corresponds to decrease $\beta$ while keeping $\mu$
(roughly) constant. As predicted, we observe that the communication saving
experienced by ACC-SONATA-F improves with the local sample size, as the method
takes advantage of the local function structure, a feature that first-order
methods are lacking.

## References

  * [Arjevani and Shamir, 2015] Arjevani, Y. and Shamir, O. (2015).  Communication complexity of distributed convex learning and optimization.  In Advances in neural information processing systems, pages 1756–1764.
  * [Bekkerman et al., 2011] Bekkerman, R., Bilenko, M., and Langford, J. (2011).  Scaling up Machine Learning: Parallel and Distributed Approaches.  Cambridge University Press.
  * [Berthier et al., 2020] Berthier, R., Bach, F., and Gaillard, P. (2020).  Accelerated gossip in networks of given dimension using jacobi polynomial iterations.  SIAM J. on Mathematics of Data Science, 1:24–47.
  * [Chang and Lin, 2011] Chang, C.-C. and Lin, C.-J. (2011).  Libsvm: a library for support vector machines.  ACM transactions on intelligent systems and technology (TIST), 2(3):1–27.
  * [d’Aspremont et al., 2021] d’Aspremont, A., Scieur, D., and Taylor, A. B. (2021).  Acceleration methods.  ArXiv preprint, arXiv:2101.09545.
  * [Deng, 2012] Deng, L. (2012).  The mnist database of handwritten digit images for machine learning research.  IEEE Signal Processing Magazine, 29(6):141–142.
  * [Dragomir et al., 2019] Dragomir, R.-A., Taylor, A., d’Aspremont, A., and Bolt, J. (2019).  Optimal complexity and certification of bregman first-order methods.  arXiv preprint, arXiv:1911.08510.
  * [Fan et al., 2019] Fan, J., Guo, Y., and Wang, K. (2019).  Communication-efficient accurate statistical estimation.  arXiv:1906.04870.
  * [Hendrikx et al., 2020a] Hendrikx, H., Bach, F., and Massoulie, L. (2020a).  Dual-free stochastic decentralized optimization with variance reduction.  In Advances in Neural Information Processing Systems.
  * [Hendrikx et al., 2020b] Hendrikx, H., Xiao, L., Bubeck, S., Bach, F., and Massoulie, L. (2020b).  Statistically preconditioned accelerated gradient method for distributed optimization.  In International Conference on Machine Learning, pages 4203–4227. PMLR.
  * [Hongzhou et al., 2015] Hongzhou, L., Mairal, J., and Harchaoui, Z. (2015).  A universal catalyst for first-order optimization.  In Advances in Neural Information Processing Systems, pages 3384–3392.
  * [Hongzhou et al., 2018] Hongzhou, L., Mairal, J., and Harchaoui, Z. (2018).  Catalyst acceleration for first-order convex optimization: from theory to practice.  J. of Machine Learning Research, pages 1–54.
  * [Kovalev et al., 2020] Kovalev, D., Salim, A., and Richtarik, P. (2020).  Optimal and practical algorithms for smooth and strongly convex decentralized optimization.  Advances in Neural Information Processing Systems, 33.
  * [Li et al., 2020] Li, B., Cen, S., Chen, Y., and Chi, Y. (2020).  Communication-efficient distributed optimization in networks with gradient tracking and variance reduction.  J. of Machine Learning Research, (180):1–51.
  * [Li et al., 2018] Li, H., Fang, C., Yin, W., and Lin, Z. (2018).  A sharp convergence rate analysis for distributed accelerated gradient methods.  arXiv preprint arXiv:1810.01053.
  * [Li and Lin, 2020] Li, H. and Lin, Z. (2020).  Revisiting extra for smooth distributed optimization.  SIAM Journal on Optimization, 30(3):1795–1821.
  * [Lian et al., 2017] Lian, X., Zhang, C., Zhang, H., Hsieh, C.-J., Zhang, W., and Liu., J. (2017).  Can decentralized algorithms outperform centralized algorithms? A case study for decentralized parallel stochastic gradient descent.  In Advances in Neural Information Processing Systems, pages 5330–5340.
  * [Lu et al., 2020] Lu, H., Freund, R. M., and Nesterov, Y. (2020).  Relatively smooth convex optimization by first-order methods, and applications.  SIAM J. on Optimization, 28(1):333–354.
  * [Reddi et al., 2016] Reddi, S. J., Konecny, J., Richtarik, P., Poczos, B., and Smola, A. (2016).  Aide: Fast and communication efficient distributed optimization.  arXiv preprint arXiv:1608.06879.
  * [Rogozin et al., 2020] Rogozin, A., Lukoshkin, V., Gasnikov, A., Kovalev, D., and Shulgin, E. (2020).  Towards accelerated rates for distributed optimization over time-varying networks.  arXiv preprint arXiv:2009.11069.
  * [Scaman et al., 2017] Scaman, K., Bach, F., Bubeck, S., Lee, Y. T., and Massoulie, L. (2017).  Optimal algorithms for smooth and strongly convex distributed optimization in networks.  In international conference on machine learning, pages 3027–3036. PMLR.
  * [Shalev-Shwartz and Zhang, 2013] Shalev-Shwartz, S. and Zhang, T. (2013).  Stochastic dual coordinate ascent methods for regularized loss.  J. of Machine Learning Research, (14):567–599.
  * [Shamir et al., 2014] Shamir, O., Srebro, N., and Zhang, T. (2014).  Communication-efficient distributed optimization using an approximate newton-type method.  In International conference on machine learning, pages 1000–1008. PMLR.
  * [Sun et al., 2019] Sun, Y., Daneshmand, A., and Scutari, G. (2019).  Distributed optimization based on gradient-tracking revisited: Enhancing convergence rate via surrogation.  arXiv preprint arXiv:1905.02637.
  * [Uribe et al., 2020] Uribe, C. A., Lee, S., Gasnikov, A., and Nedić, A. (2020).  A dual approach for optimal algorithms in distributed optimization over networks.  In 2020 Information Theory and Applications Workshop (ITA), pages 1–37. IEEE.
  * [Wien, 2011] Wien, A. (2011).  Iterative solution of large linear systems.  Lecture Notes, TU Wien.
  * [Ye et al., 2020a] Ye, H., Luo, L., Zhou, Z., and Zhang, T. (2020a).  Multi-consensus decentralized accelerated gradient descent.  arXiv preprint arXiv:2005.00797.
  * [Ye et al., 2020b] Ye, H., Zhou, Z., Luo, L., and Zhang, T. (2020b).  Decentralized accelerated proximal gradient descent.  Advances in Neural Information Processing Systems, 2020.
  * [Yuan and Li, 2020] Yuan, X.-T. and Li, P. (2020).  On convergence of distributed approximate newton methods: Globalization, sharper bounds and beyond.  J. Mach. Learn. Res., 21:206–1.
  * [Zhang and Lin, 2015] Zhang, Y. and Lin, X. (2015).  Disco: Distributed optimization for self-concordant empirical loss.  In International conference on machine learning, pages 362–370. PMLR.

Appendix

This document serves as supplementary material of the paper entitled
“Acceleration in Distributed Optimization under Similarity”. It contains
additional theoretical and numerical results, along with all the proofs of the
theorems presented in the main paper. Specifically

  * Sec. A contains the lower complexity bounds for the class of problems and oracle algorithms of interest;

  * Sec. B presents additional numerical results for different classes of problems on synthetic and real data;

  * Sec. C provides the proof of Theorem 4 and Theorem 5;

  * Sec. D customizes ACC-SONATA to star-topologies (master/workers architectures).

##  Appendix A Lower complexity bounds over mesh networks under similarity

Given (P) over a mesh network $\mathcal{G}$, we consider the following general
class of distributed algorithms, which generalize the oracle model [Arjevani
and Shamir, 2015] for centralized schemes (star architectures) and smooth
($r=0$) instances of (P).

######  Definition 6 (Distributed oracle).

Each agent $i$ has its own local memory $\mathcal{M}_{i}=\\{0\\}$, updated as
follows:

$\bullet$ Local computation: Between communication rounds, each agent $i$
iteratively computes and adds to $\mathcal{M}_{i}$ some finite number of
points $x$, each satisfying

$\displaystyle\tau_{1}\,x+\tau_{2}\nabla
f_{i}(x)+\tau_{3}\,g\,\in\textnormal{span}\bigg{\\{}x^{\prime},\,\nabla
f_{i}(x^{\prime}),$
$\displaystyle\quad(\nabla^{2}f_{i}(x^{\prime})+D)x^{\prime\prime},(\nabla^{2}f_{i}(x^{\prime})+D)^{-1}x^{\prime\prime}\bigg{\\}},$

for given $x^{\prime},x^{\prime\prime}\in\mathcal{M}_{i}$, $g\in\partial
r(x)$, $\tau_{1},\tau_{2},\tau_{3}\geq 0$ such that
$\tau_{1}+\tau_{2}+\tau_{3}>0$; and $D$ is a diagonal matrix such that the
inverse matrices above exist;

$\bullet$ Communications: After every communication round, each agent $i$
updates $\mathcal{M}_{i}$ according to

$\mathcal{M}_{i}=\textnormal{span}\left\\{\underset{(i,j)\in\mathcal{E}}{\cup}M_{j}\right\\};\vspace{-0.3cm}$

$\bullet$ Output: $x_{i}\in\mathcal{M}_{i}$, for all $i\in[m]$.

The above procedure models a fairly general class of distributed algorithms
over graphs. Computations at each node are based on linear operations
involving current or past iterates, local gradients, and vector products with
local Hessians and their inverses; exact optimization of local subproblem
involving such quantities or proximal solutions are also incorporated. During
communications, the agents can share with their neighbors any of the vectors
they have computed up until that time.

The following result provides lower complexity bounds for solving Problem (P)
by any distributed algorithm in $\mathcal{A}$–this is an extension of
[Arjevani and Shamir, 2015, Th. 1] to the decentralized setting.

######  Theorem 7.

For any $\mu\in[0,1),\,\beta\in(0,1),$ and $\rho\in[0,1)$, there exist (i) an
instance of (P) (with sufficiently large $d$ and solution $x^{\star}$)
satisfying Assumption 1, with $r=0$, $f_{i}$’s being $\beta$-similar and
$f(x)=\frac{1}{m}\sum_{i=1}^{m}f_{i}(x)$ being $1$-smooth and $\mu$-strongly
convex; and (ii) a gossip matrix $W$ satisfying Assumption 3 over the graph
$\mathcal{G}$ with parameter $\rho$ such that for any distributed algorithm in
$\mathcal{A}$ using the matrix $W$ in the communications, the number of
communication rounds $N$ required for obtaining the solution $x_{i}$’s such
that $u({x}_{i})-u(x^{\star})\leq\epsilon$, for all $i\in[m]$, is

$N=\Omega\left(\sqrt{\frac{\beta/\mu}{1-\rho}}\log\left(\frac{\mu\left\lVert
x^{\star}\right\rVert^{2}}{\varepsilon}\right)\right).$

###### Proof.

The case of fully connected networks ($\rho=0$) has been already studied in
[Arjevani and Shamir, 2015, Th. 1]. Therefore, here we assumed $\rho>0$. Our
proof is inspired by [Scaman et al., 2017, Th. 2], with some key differences
due to the different setting of our problem.

We first introduce the local cost functions of agents numbered as
$1,2,\ldots,m.$ With $\zeta=\frac{1}{32}$, define the following two subsets of
agents:

$\mathcal{A}_{l}=\left\\{i\big{|}1\leq i\leq\lceil\zeta
m\rceil\right\\}\quad\text{and}\quad\mathcal{A}_{r}=\left\\{i\big{|}\lfloor(1-\zeta)m\rfloor+1\leq
i\leq m\right\\}.$

We then define for each agent $i$ the cost function
$f_{i}:\,\ell^{2}\to\mathbb{R}$ as

$f_{i}(x)=\begin{cases}\frac{\beta(1-\mu)}{8}\frac{m}{\lceil\zeta
m\rceil}x^{\top}A_{1}x-\frac{\beta(1-\mu)}{4}\frac{m}{\lceil\zeta
m\rceil}e_{1}^{\top}x+\frac{\mu}{2}\left\lVert x\right\rVert^{2},&1\leq
i\leq\lceil\zeta m\rceil\\\ \frac{\mu}{2}\left\lVert
x\right\rVert^{2},&\lceil\zeta m\rceil+1\leq i\leq\lfloor(1-\zeta)m\rfloor\\\
\frac{\beta(1-\mu)}{8}\frac{m}{\lceil\zeta
m\rceil}x^{\top}A_{2}x+\frac{\mu}{2}\left\lVert
x\right\rVert^{2},&\lfloor(1-\zeta)m\rfloor+1\leq i\leq m\end{cases}$

with

$A_{1}\triangleq\begin{bmatrix}1&0&0&0&0&\cdots\\\ 0&1&-1&0&0&\cdots\\\
0&-1&1&0&0&\cdots\\\ 0&0&0&1&-1&\cdots\\\ 0&0&0&-1&1&\cdots\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\ddots\end{bmatrix},\qquad
A_{2}\triangleq\begin{bmatrix}1&-1&0&0&0&\cdots\\\ -1&1&0&0&0&\cdots\\\
0&0&1&-1&0&\cdots\\\ 0&0&-1&1&0&\cdots\\\ 0&0&0&0&1&\cdots\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\ddots\end{bmatrix}.$

We define the distance between the two subsets $\mathcal{A}_{l}$ and
$\mathcal{A}_{r}$ as $d_{c}.$ It follows that, to have at least one non-zero
element in the $k$th entry of the local copies of agents in both of the above
two subsets, one must perform at least $k$ local computation steps and
$(k-1)d_{c}$ communication steps. The number of communication rounds required
to obtain $f(\hat{x})-f(x^{\star})\leq\epsilon$ is thus
$\Omega\left(d_{c}\sqrt{\frac{\beta}{\mu}}\log\left(\frac{\mu\left\lVert
x^{\star}\right\rVert^{2}}{\epsilon}\right)\right).$

We then describe the communication graph for the given $\rho.$ For $m\geq 2$,
we define $\rho_{m}=\frac{\rho}{2+\rho}+\frac{2}{2+\rho}\cos{\frac{\pi}{m}}.$
Since $\rho_{2}=\frac{\rho}{2+\rho}<\rho$ and $\lim_{k\to\infty}\rho_{k}=1,$
we know that there exist $m$ such that $\rho_{m}<\rho\leq\rho_{m+1}.$ We
discuss the cases of $m\geq 3$ and $m=2$ separately:

i)  $m\geq 3.$ We begin defining a Laplacian matrix $L_{m,a}$ for a line graph
composed of $m$ nodes; specifically:

$L_{m,a}=L_{m,a}^{\top},\quad L_{m,a}\,1=0,\quad\text{and}\quad
L_{m,a}(i,i+1)=a\,\mathbbm{1}\\{i=1\\}-1,$

with $a\in[0,1)$, where $\mathbbm{1}\\{\bullet\\}$ is the indicator function.
We then define

$W_{m,a}=I-\frac{1}{2+\rho}L_{m,a}.$

It is not difficult to check that $W_{m,a}$ satisfies Assumption 3, with
$\left\lVert W_{m,0}-1\,1^{\top}/m\right\rVert=\rho_{m}.$ Furthermore, since
$\left\lVert W_{m,1}-1\,1^{\top}/m\right\rVert=1,$ by continuity, we know that
there exists an $a\in(0,1)$ such that $\left\lVert
W_{m,a}-1\,1^{\top}/m\right\rVert=\rho.$ In addition, we have

$1-\rho\geq
1-\rho_{m+1}=\frac{2}{2+\rho}\left(1-\cos{\frac{\pi}{m+1}}\right)\geq\frac{2}{3}\left(1-\cos{\frac{\pi}{m+1}}\right)\overset{(*)}{\geq}\frac{8}{3}\frac{1}{(m+1)^{2}}.$

Note that $(*)$ is due to that $\cos{\frac{\pi}{n}}\leq 1-\frac{4}{n^{2}}$ for
$n\geq 4.$ Equivalently,

$m\geq\sqrt{\frac{8}{3(1-\rho)}}-1.$

The distance between the two subsets is thus

$d_{c}=\lfloor(1-\zeta)m\rfloor+1-\lceil\zeta
m\rceil\geq\frac{15}{16}m-1\geq\frac{15}{16}\left(\sqrt{\frac{8}{3(1-\rho)}}-1\right)-1\overset{(*)}{\geq}\frac{4}{25}\sqrt{\frac{1}{1-\rho}}.$

Note that $(*)$ is due to that $\rho>\rho_{3}>\frac{1}{2}.$

ii) $m=2.$ In this case, we have $\rho\leq\rho_{3}=\frac{1+\rho}{2+\rho}$ and
equivalently, $\rho\in(0,\frac{\sqrt{5}-1}{2}].$ We define the Laplacian
matrix for a complete graph of $3$ nodes as: $L_{a}=L_{a}^{\top}$,
$L_{a}\,1=0$, and for $i<j$, $L_{a}(i,j)=a\,\mathbbm{1}\\{i=1,\,j=2\\}-1.$ We
then set $W_{a}=I-\frac{1}{3}L_{a}.$ Due to $\left\lVert
W_{1}-1\,1^{\top}/3\right\rVert=\frac{2}{3}>\frac{\sqrt{5}-1}{2}$ and
$\left\lVert W_{0}-1\,1^{\top}/3\right\rVert=0,$ there exists an $a\in(0,1)$
such that $\left\lVert W_{a}-1\,1^{\top}/3\right\rVert=\rho.$ In this case, we
have

$d_{c}=1\geq\sqrt{\frac{3-\sqrt{5}}{2}}\frac{1}{\sqrt{1-\rho}}.$

Therefore, for any $\rho\in(0,1)$, there exists a communication matrix $W$
satisfying Assumption 3, with $\left\lVert W-1\,1^{\top}/m\right\rVert=\rho$,
such that the number of communication rounds required to obtain
$f(\hat{x})-f(x^{\star})\leq\epsilon$ is

$\Omega\left(\sqrt{\frac{\beta}{\mu\,(1-\rho)}}\log\left(\frac{\mu\left\lVert
x^{\star}\right\rVert^{2}}{\epsilon}\right)\right).$

∎

##  Appendix B Additional Numerical Results

This section complements Sec. 4 of the main paper, providing additional
numerical results and details on the experiments presented therein.
Specifically, we consider the following problems and data sets:

  * •

Sec. B.2 studies hinge-loss minimization on MNIST [Deng, 2012] and HIGGS
[Chang and Lin, 2011, from LIBSVM] datasets;

  * •

Sec. B.3 considers logistic regression problems on the SUSY dataset [Chang and
Lin, 2011, from LIBSVM] datasets[cite].

###  B.1 Setting and algorithm tuning

The setting of the experiments is the same as the one described in Sec. 4,
expect of course for the instance of (P).

The tuning of the simulated algorithms follows the direction as in the
associated papers. To set the free parameters, an estimate of the smoothness
constants $L_{i}$’s and strong-convexity constants $\mu_{i}$’s is needed. We
use the following bounds for the hinge and logistic losses. In both cases, the
local losses have the following structure:

$\displaystyle
f_{i}(x)=\frac{1}{n}\sum_{j=1}^{n}\ell\left(b_{i}^{j}\cdot\left\langle
x,a_{i}^{j}\right\rangle\right)+\frac{\lambda}{2}\left\lVert
x\right\rVert^{2},$

where $a_{i}^{j}\in\mathbb{R}^{d}$ are the feature vectors and
$b_{i}^{j}\in\\{-1,1\\}$ are the associated labels. The Hessian matrix of
$f_{i}$ is

$\nabla^{2}f_{i}(x)=\frac{1}{n}\sum_{j=1}^{n}\ell^{\prime\prime}\left(b_{i}^{j}\left\langle
x,a_{i}^{j}\right\rangle\right)(b_{i}^{j})^{2}a_{i}^{j}(a_{i}^{j})^{\top}+\lambda
I,$

where $\ell^{\prime\prime}$ denotes the second derivative of the loss
$\ell:\mathbb{R}\to\mathbb{R}$. Under the assumption that
$\ell^{\prime\prime}\leq C_{\ell}$,

$\lambda
I\preceq\nabla^{2}f_{i}(x)\preceq\frac{1}{n}\sum_{j=1}^{n}C_{\ell}(b_{i}^{j})^{2}a_{i}^{j}(a_{i}^{j})^{T}+\lambda
I\triangleq H_{i}.$

In particular, $C_{\ell}=1$ for the smooth hinge loss and $C_{\ell}=1/4$ for
the logistic loss.

Based on $H_{i}$ above, we use $\lambda$ as estimate of $\mu_{i}$ and the
largest eigenvalue of $H_{i}$ as that of $L_{i}$. Furthermore, for the smooth
constant $L$ of the average loss $f$ we use the overestimate
$\hat{L}\triangleq\frac{1}{m}\sum_{i=1}^{m}L_{i}\geq L$; and for $\beta$ we
use $\hat{\beta}\triangleq\max_{i\in[m]}\left\lVert
H_{i}-\frac{1}{m}\sum_{i=1}^{m}H_{i}\right\rVert$.

The other tuning of the algorithms is as follows:

  * •

ACC-SONATA-F: The inner iterations $T$ of SONATA-F and SONATA-L are set to
$\left\lceil\frac{7}{5}\cdot\log\frac{\hat{L}}{\mu}\right\rceil$ and
$\left\lceil\log\frac{\hat{\beta}}{\mu}\right\rceil$, respectively;

  * •

Mudag: the number of inner loops is set to
$\left\lceil\frac{1}{5\sqrt{1-\rho}}\log\frac{L_{mx}}{\mu}\right\rceil$;

  * •

ACC-EXTRA: the number inner loops is set to
$\left\lceil\frac{1}{5(1-\rho)}\log\frac{L_{mx}}{\mu(1-\rho)}\right\rceil$;

  * •

APM-C: the number of inner loops is set to
$T_{k}=\left\lceil\frac{k\sqrt{\mu/L_{mx}}}{100(\sqrt{1-\rho})}\right\rceil$.

Whenever the subproblems of the agents do not have a closed form solution, the
gradient algorithm is employed, and terminated when an accuracy of $10^{-10}$
is reached on the norm gradient.

We are now ready to describe our experiments.

###  B.2  Hinge loss minimization

Consider the following instance of (P):

$\min_{x\in\mathbb{R}^{d}}\frac{1}{m}\sum_{i=1}^{m}\frac{1}{n}\sum_{j=1}^{n}\ell_{s}(b_{i}^{j}\cdot\left\langle
x,a_{i}^{j}\right\rangle)+\frac{\lambda}{2}\left\lVert x\right\rVert^{2},$

where $\ell_{s}$ is the smooth hinge loss, defined as:

$\ell_{s}(t)=\left\\{\begin{aligned} &0&t>1,\\\
&\frac{1}{2}(t-1)^{2}&t\in[0,1],\\\ &\frac{1}{2}-t&t<0.\end{aligned}\right.$

We consider two datasets for the above problem, namely the MNIST and the
HIGGS. We use the label $1$ for images of digit $4$ and $-1$ for the others.
Results are summarized in Fig. 3 (MNIST) and Fig. 4 (HIGGS).

Specifically, Fig. 3 (resp. Fig. 4)-left-panel plots the optimality gap
$\frac{1}{m}\sum_{i=1}^{m}\|x_{i}^{k}-{x}_{\textnormal{op}}\|^{2}$ versus the
total number of communications, achieved by the different algorithms, where
${x}_{\textnormal{op}}$ is the optimal solution of the problem (estimated
running the gradient algorithm up to a precision of $10^{-10}$ on the norm
gradient). In the mid-panel, we plot the number of communications to drive the
optimality gap below $10^{-4}$ versus the total sample size; we consider four
sizes, namely: $1.8\times 10^{4}$, $3\times 10^{4}$, $4.8\times 10^{4}$ and
$6\times 10^{4}$ for the MNIST dataset and $1.2\times 10^{5}$, $2.4^{5}$,
$4.8\times 10^{5}$ and $9\times 10^{5}$ for the HIGGS dataset. The righ-panel
shows the same result for ACC-SONATA on a re-scaled y-axes, to highlight the
decreasing number of communications with the local sample size. These results
confirm the competitive performance of the proposed schemes over the state of
the art also on real data.

###  B.3 Logistic regression

We consider here a logistic regression problem on SUSY dataset:

$\min_{x\in\mathbb{R}^{d}}\frac{1}{m}\sum_{i=1}^{m}\frac{1}{n}\sum_{j=1}^{n}\ell_{r}\left(b_{i}^{j}\cdot\left\langle
x,a_{i}^{j}\right\rangle\right)+\frac{\lambda}{2}\left\lVert
x\right\rVert^{2},$ (13)

where $\ell_{r}=\log(1+e^{-t})$.

Fig. 5 (left-panel) is a plot of optimality gap
$\frac{1}{m}\sum_{i=1}^{m}\|x_{i}^{k}-{x}_{\textnormal{op}}\|^{2}$ versus
total number of communications. In the right-panel, we plot the number of
communications to obtain optimality gap $10^{-4}$ versus the size of total
data samples, namely: $9\times 10^{3}$, $9.9\times 10^{4}$, $2.1\times 10^{5}$
and $6\times 10^{5}$. Consistently with the other results, also for this class
of (nonquadratic) problems, ACC-SONATA exhibits favorable performance.

Figure 3: Hinge loss minimization, MNIST dataset. (left panel): optimality gap
versus total number of communications; (mid panel): number of communications
to reach a precision of $10^{-4}$ versus (total) sample; (right panel): Same
as in the mid panel on a different scale of the y-axes.

Figure 4: Hinge loss minimization, HIGGS dataset. (left panel): optimality gap
versus total number of communications; (mid panel): number of communications
to reach a precision of $10^{-4}$ versus (total) sample; (right panel): Same
as in the mid panel on a different scale of the y-axes.

Figure 5: Logistic regression, SUSY dataset. (left panel): optimality gap
versus total number of communications; (right panel): number of communications
to reach a precision of $10^{-4}$ versus (total) sample.

##  Appendix C Proof of Theorem 4 and Theorem 5

###  C.1 Sketch of the proof

As discussed in Sec. 3 [see (S.1)′ and (S.2)′], ACC-SONATA can be interpreted
as an accelerated version of the inexact proximal point algorithm. The
challenge here is finding a suitable notion of inexactness for proximal
operations that captures all consensus errors (on $x,y,z$-variables) while
being implementable in the distributed setting and, at the same time, retains
(for the outer loop) the convergence rate, $\alpha$, of the exact accelerated
proximal method. Our path towards this goal consists in the following two
steps:

  * •

Step 1 (inexactness and outer-loop convergence): We introduce an inner
termination condition for the SONATA algorithm [see (18), Sec. C.2], serving
as inexact notion of approximate proximal solution. This criterion hinges on a
proper potential function controlling the decrease of the optimality gap up to
consensus/tracking errors. Roughly speaking, this quantifies the amount of
errors in the minimization of $u_{k}$ [see (6)] that can be tolerated to
preserve the convergence of the outer loop at the desired rate
$\alpha=\sqrt{\frac{\mu}{\mu+\delta}}$. Convergence of the outer loop at such
a rate is established by introducing a proper potential function [see (19),
Sec. C.2]. Such a potential function certifies linear convergence of the
optimality gap

$\Delta(x^{k})=\max\left(\frac{1}{m}\sum_{i=1}^{m}u(x_{i}^{k})-u^{\star},\frac{1}{m}\sum_{i=1}^{m}\left\lVert
x_{i}^{k}-\bar{x}^{k}\right\rVert^{2}\right)$

at rate $\alpha$. In the setting of Theorem 4 (i.e., $\delta=\beta-\mu$) and
Theorem 5 (i.e., $\delta=L-\mu$), $\alpha$ reads

$\alpha=\sqrt{\frac{\mu}{\beta}}\quad\text{and}\quad\alpha=\sqrt{\frac{1}{\kappa}},$
(14)

respectively.

  * •

Step 2 (inner-loop convergence): We introduce a refined analysis of the SONATA
algorithm based on a new potential function certifying that the inner
termination criterion defined in Step 1 is met in $T=\tilde{\mathcal{O}}(1)$
number of iterations.

Combining Step 1 and Step 2, we can then conclude that ACC-SONATA achieves an
$\varepsilon$-solution of Problem (P) in

${\mathcal{O}}\left(T\frac{1}{\alpha}\,\log\frac{1}{\varepsilon}\right)=\widetilde{\mathcal{O}}\left(\frac{1}{\alpha}\,\log\frac{1}{\varepsilon}\right)=\left\\{\begin{matrix}\widetilde{\mathcal{O}}\left(\sqrt{\dfrac{\beta}{\mu}}\,\log\dfrac{1}{\varepsilon}\right),&\text{if
}\delta=\beta-\mu\,\,(\text{Theorem\leavevmode\nobreak\
\ref{thm:cata_net}}),&\vskip 6.0pt plus 2.0pt minus 2.0pt\\\
\widetilde{\mathcal{O}}\left(\sqrt{\kappa}\,\log\dfrac{1}{\varepsilon}\right),&\text{if
}\delta=L-\mu\,\,(\text{Theorem\leavevmode\nobreak\
\ref{crl:cata_net_lin}}),&\end{matrix}\right.$

total number of communications. The formal proof of Step 1 and Step 2 is given
in Sec. C.2 and Sec. C.3, respectively.

Notation: Before detailing the two steps, we introduce some notation used
throughout the proofs. For each $i\in[m]$, we denote by $x_{i}^{k,t}$ and
$y_{i}^{k,t}$ the decision variable and tracking variable of agent $i$ after
$t=0,\ldots,T-1$ inner iterations in the $k$-th outer iteration, respectively.
Clearly, it must be

$x_{i}^{k,0}=x_{i}^{k},\quad
y^{k,0}=y_{i}^{k}+\delta\left(z_{i}^{k-1}-z_{i}^{k}\right),\quad
x_{i}^{k,T}=x_{i}^{k+1},\quad\text{and}\quad y^{k,T}=y_{i}^{k+1}.$

Furthermore, we define

$x^{k+1,\star}=\mathop{\rm argmin}_{x\in\mathbb{R}^{d}}u_{k}(x).$

For any vector $x=[x_{1}^{\top},\ldots,x_{m}^{\top}]^{\top}$, we use $\bar{x}$
to denote the average of its $d$-dimensional blocks $x_{i}$’s, that is,
$\bar{x}=\frac{1}{m}\sum_{i=1}^{m}x_{i}$.

Consensus and tracking errors associated with the iterates $x_{i}^{k,t}$ and
$y_{i}^{k,t}$ are defined as

$\left\lVert
x_{\perp}^{k,t}\right\rVert^{2}=\frac{1}{m}\sum_{i\in[m]}\left\lVert
x_{i}^{k,t}-\bar{x}^{k,t}\right\rVert^{2}\quad\text{and}\quad\left\lVert
y_{\perp}^{k,t}\right\rVert^{2}=\frac{1}{m}\sum_{i\in[m]}\left\lVert
y_{i}^{k,t}-\bar{x}^{y,t}\right\rVert^{2},$ (15)

respectively.

###  C.2 Step 1: Inexactness and outer loop convergence

Our inexact notion of approximate proximal solution of the minimization of
$u_{k}$ [see (6)] is defined in terms of the (average) optimality gap,

$g^{k,t}\triangleq\frac{1}{m}\sum_{i=1}^{m}\left(u_{k}(x_{i}^{k,t})-u_{k}^{\star}\right),\quad
k=0,1,\ldots,\,\,t=0,\ldots,T,$ (16)

with $u_{k}^{\star}=\min_{x\in\mathbb{R}^{d}}u_{k}(x)$, and the consensus and
tracking errors (15), cumulatively captured by

$e^{k,t}\triangleq c_{x}\left\lVert
x_{\perp}^{k,t}\right\rVert^{2}+c_{y}\left\lVert
y_{\perp}^{k,t}\right\rVert^{2},\quad k=0,1,\ldots,\,\,t=0,\ldots,T,$ (17)

where $c_{x},c_{y}>0$ are suitably defined universal constants. Specifically,
under the following event

$g^{k,T}+e^{k,T}\leq\epsilon^{k+1},\quad k=0,1,\ldots,$ (18)

with $\\{\epsilon^{k}\\}$ being a suitably defined geometrically-vanishing
positive sequence, we establish linear decay of the following potential
function along the outer-loop iterates of ACC-SONATA:

$\displaystyle
P^{k}\triangleq\frac{1}{m}\sum_{j=1}^{m}\left(u(x_{j}^{k})-u^{\star}\right)+\frac{1}{m}\sum_{j=1}^{m}\frac{\mu}{2}\left\lVert
x_{j}^{k-1}+\frac{1}{\alpha}(x_{j}^{k}-x_{j}^{k-1})-x^{\star}\right\rVert^{2}+e^{k-1,T},\quad
k=0,1,\ldots,$ (19)

where we set $e^{-1,T}=0$.

In Step 2 we will show that (18) can be met by running SONATA for
$T=\tilde{\mathcal{O}}(1)$ iterations.

######  Proposition 8.

Consider problem (P) under Assumption 1, with optimal value function
$u^{\star}$. Let $\\{(x_{i}^{k},z_{i}^{k})_{i\in[m]}\\}$ be the sequence
generated by ACC-SONATA under (18), with error sequence
$\epsilon^{k}=P^{0}\,(1-c\alpha)^{k}$, $k=1,\ldots$, where $c$ is any given
constant in $(0,1)$. Then, the potential function $P^{k}$ in (19) satisfies:

$\displaystyle P^{k}\leq c_{2}\,P^{0}\,(1-c\cdot\alpha)^{k},$ (20)

with

$c_{2}=\frac{(2+\sqrt{c_{1}})^{2}}{\left(\sqrt{\frac{1-c\cdot\alpha}{1-\alpha}}-1\right)^{2}(1-\alpha)}\quad\text{and}\quad
c_{1}=1+\frac{\delta}{c_{x}}\frac{\frac{3}{2}(1-c\alpha)^{2}+5-4c\,\alpha}{(1-c\alpha)^{2}}.$
(21)

Therefore,

$\max\left(\frac{1}{m}\sum_{i=1}^{m}u(x_{i}^{k})-u^{\star},\frac{1}{m}\sum_{i=1}^{m}\left\lVert
x_{i}^{k}-\bar{x}^{k}\right\rVert^{2}\right)=\mathcal{O}\left((1-c\cdot\alpha)^{k}\right).$

###### Proof.

We provide a constructive proof for the choice of the error sequence
$\epsilon^{k}$ in (18) and potential function (19), yielding (31).

Under (18), using the definition of $u_{k}$, we have: for all
$x_{j}\in\text{dom}\,r$,

$\begin{split}&\frac{1}{m}\sum_{j=1}^{m}u(x_{j}^{k+1})=\frac{1}{m}\sum_{j=1}^{m}u_{k}(x_{j}^{k+1})-\frac{\delta}{2m^{2}}\sum_{i,j}\left\lVert
x_{j}^{k+1}-z_{i}^{k}\right\rVert^{2}\leq
u_{k}^{\star}+\epsilon^{k+1}-e^{k,T}-\frac{\delta}{2m^{2}}\sum_{i,j}\left\lVert
x_{j}^{k+1}-z_{i}^{k}\right\rVert^{2}\\\
&\leq\frac{1}{m}\sum_{j=1}^{m}\left(u_{k}(x_{j})-\frac{\mu+\delta}{2}\left\lVert
x_{j}-x^{k+1,\star}\right\rVert^{2}\right)+\epsilon^{k+1}-e^{k,T}-\frac{\delta}{2m^{2}}\sum_{i,j}\left\lVert
x_{j}^{k+1}-z_{i}^{k}\right\rVert^{2}\\\
&=\frac{1}{m}\sum_{j=1}^{m}u_{k}(x_{j})-\frac{1}{m}\sum_{j=1}^{m}\frac{\mu+\delta}{2}\left\lVert
x_{j}-x_{j}^{k+1}\right\rVert^{2}+\frac{1}{m}\sum_{j=1}^{m}\frac{\mu+\delta}{2}\left\lVert
x_{j}-x_{j}^{k+1}\right\rVert^{2}-\frac{1}{m}\sum_{j=1}^{m}\frac{\mu+\delta}{2}\left\lVert
x_{j}-x^{k+1,\star}\right\rVert^{2}\\\
&\qquad+\epsilon^{k+1}-e^{k,T}-\frac{\delta}{2m^{2}}\sum_{i,j}\left\lVert
x_{j}^{k+1}-z_{i}^{k}\right\rVert^{2}\\\
&=\frac{1}{m}\,\sum_{j=1}^{m}u(x_{j})+\frac{1}{m^{2}}\sum_{i,j}\frac{\delta}{2}\left\lVert
x_{j}-z_{i}^{k}\right\rVert^{2}-\frac{1}{m}\sum_{j=1}^{m}\frac{\mu+\delta}{2}\left\lVert
x_{j}-x_{j}^{k+1}\right\rVert^{2}+\frac{1}{m}\sum_{j=1}^{m}\frac{\mu+\delta}{2}\left\lVert
x_{j}-x_{j}^{k+1}\right\rVert^{2}\\\
&\qquad-\frac{1}{m}\sum_{j=1}^{m}\frac{\mu+\delta}{2}\left\lVert
x_{j}-x^{k+1,\star}\right\rVert^{2}+\epsilon^{k+1}-e^{k,T}-\frac{\delta}{2m^{2}}\sum_{i,j}\left\lVert
x_{j}^{k+1}-z_{i}^{k}\right\rVert^{2}\\\
&=\frac{1}{m}\,\sum_{j=1}^{m}u(x_{j})-\frac{\mu}{2m}\sum_{j=1}^{m}\left\lVert
x_{j}-x_{j}^{k+1}\right\rVert^{2}+\frac{\delta}{2m}\sum_{j=1}^{m}\left(\frac{1}{m}\sum_{i=1}^{m}\left\lVert
x_{j}-z_{i}^{k}\right\rVert^{2}-\left\lVert
x_{j}-x_{j}^{k+1}\right\rVert^{2}-\frac{1}{m}\sum_{i=1}^{m}\left\lVert
x_{j}^{k+1}-z_{i}^{k}\right\rVert^{2}\right)\\\
&\qquad+\frac{1}{m}\sum_{j=1}^{m}\frac{\mu+\delta}{2}\left\lVert
x_{j}-x_{j}^{k+1}\right\rVert^{2}-\frac{1}{m}\sum_{j=1}^{m}\frac{\mu+\delta}{2}\left\lVert
x_{j}-x^{k+1,\star}\right\rVert^{2}+\epsilon^{k+1}-e^{k,T}\\\
&\leq\frac{1}{m}\,\sum_{j=1}^{m}u(x_{j})-\frac{\mu}{2m}\sum_{j=1}^{m}\left\lVert
x_{j}-x_{j}^{k+1}\right\rVert^{2}+\frac{\delta}{m}\sum_{j=1}^{m}\left\langle
x_{j}-x_{j}^{k+1},x_{j}^{k+1}-\bar{z}^{k}\right\rangle\\\
&\qquad-\frac{1}{m}\sum_{j=1}^{m}(\mu+\delta)\left\langle
x_{j}-x_{j}^{k+1},x_{j}^{k+1}-x^{k+1,\star}\right\rangle+\epsilon^{k+1}-e^{k,T}.\end{split}$

We define $x^{\star}\triangleq\mathop{\rm argmin}_{x\in\mathbb{R}^{d}}u(x)$
and $u^{\star}\triangleq u(x^{\star}).$ Setting $x_{j}=x^{\star}$, $j\in[m]$,
leads to

$\displaystyle\frac{1}{m}\sum_{j=1}^{m}u(x_{j}^{k+1})\leq
u^{\star}-\frac{\mu}{2m}\sum_{j=1}^{m}\left\lVert
x^{\star}-x_{j}^{k+1}\right\rVert^{2}+\frac{\delta}{m}\sum_{j=1}^{m}\left\langle
x^{\star}-x_{j}^{k+1},x_{j}^{k+1}-\bar{z}^{k}\right\rangle$ (22)
$\displaystyle\qquad-\frac{1}{m}\sum_{j=1}^{m}(\mu+\delta)\left\langle
x^{\star}-x_{j}^{k+1},x_{j}^{k+1}-x^{k+1,\star}\right\rangle+\epsilon^{k+1}-e^{k,T}.$

Similarly, setting $x_{j}=x_{j}^{k}$, we have

$\displaystyle\frac{1}{m}\sum_{j=1}^{m}u(x_{j}^{k+1})\leq\frac{1}{m}\sum_{j=1}^{m}u(x_{j}^{k})-\frac{\mu}{2m}\sum_{j=1}^{m}\left\lVert
x_{j}^{k}-x_{j}^{k+1}\right\rVert^{2}+\frac{\delta}{m}\sum_{j=1}^{m}\left\langle
x_{j}^{k}-x_{j}^{k+1},x_{j}^{k+1}-\bar{z}^{k}\right\rangle$ (23)
$\displaystyle\qquad-\frac{1}{m}\sum_{j=1}^{m}(\mu+\delta)\left\langle
x_{j}^{k}-x_{j}^{k+1},x_{j}^{k+1}-x^{k+1,\star}\right\rangle+\epsilon^{k+1}-e^{k,T}.$

Define the optimality gap at the beginning of the $k$-th outer iteration,
pertaining to the minimization of the original objective function $u(x)$ as:

$p^{k}\triangleq\frac{1}{m}\sum_{j=1}^{m}u(x_{j}^{k})-u^{\star}.$

Then, multiplying (22) by $\alpha$ and (23) by $(1-\alpha)$, and suming the
obtained equations, yields

$\displaystyle
p^{k+1}\leq(1-\alpha)p^{k}-\frac{\alpha\mu}{2m}\sum_{j=1}^{m}\left\lVert
x^{\star}-x_{j}^{k+1}\right\rVert^{2}-\frac{(1-\alpha)\mu}{2m}\sum_{j=1}^{m}\left\lVert
x_{j}^{k}-x_{j}^{k+1}\right\rVert^{2}$ (24)
$\displaystyle\quad+\frac{\delta}{m}\sum_{j=1}^{m}\left\langle\alpha
x^{\star}+(1-\alpha)x_{j}^{k}-x_{j}^{k+1},x_{j}^{k+1}-\bar{z}^{k}\right\rangle$
$\displaystyle\quad-\frac{\mu+\delta}{m}\sum_{j=1}^{m}\left\langle\alpha
x^{\star}+(1-\alpha)x_{j}^{k}-x_{j}^{k+1},x_{j}^{k+1}-x^{k+1,\star}\right\rangle+\epsilon^{k+1}-e^{k,T}$
$\displaystyle=(1-\alpha)p^{k}-\frac{\alpha\mu}{2m}\sum_{j=1}^{m}\left\lVert
x^{\star}-x_{j}^{k+1}\right\rVert^{2}-\frac{(1-\alpha)\mu}{2m}\sum_{j=1}^{m}\left\lVert
x_{j}^{k}-x_{j}^{k+1}\right\rVert^{2}-\frac{\delta}{m}\sum_{j=1}^{m}\left\lVert
x_{j}^{k+1}-\bar{z}^{k}\right\rVert^{2}$
$\displaystyle\quad+\frac{\delta}{m}\sum_{j=1}^{m}\left\langle\alpha
x^{\star}+(1-\alpha)x_{j}^{k}-\bar{z}^{k},x_{j}^{k+1}-\bar{z}^{k}\right\rangle$
$\displaystyle\quad-\frac{\mu+\delta}{m}\sum_{j=1}^{m}\left\langle\alpha
x^{\star}+(1-\alpha)x_{j}^{k}-x_{j}^{k+1},x_{j}^{k+1}-x^{k+1,\star}\right\rangle+\epsilon^{k+1}-e^{k,T}.$

The above is an approximate descent on $p^{k},$ up to the error term
$\epsilon^{k+1}-e^{k,T}$ and the two inner-product terms. We proceed working
first on the term $\left\langle\alpha
x^{\star}+(1-\alpha)x_{j}^{k}-\bar{z}^{k},x_{j}^{k+1}-\bar{z}^{k}\right\rangle$.
Define

$v_{j}^{k}\triangleq x_{j}^{k-1}+\frac{1}{\alpha}(x_{j}^{k}-x_{j}^{k-1}).$

We have

$\displaystyle v_{j}^{k+1}$
$\displaystyle=x_{j}^{k}+\frac{1}{\alpha}(x_{j}^{k+1}-x_{j}^{k})=\left(1-\frac{1}{\alpha}\right)x_{j}^{k}+\frac{1}{\alpha}z_{j}^{k}+\frac{1}{\alpha}(x_{j}^{k+1}-z_{j}^{k})$
$\displaystyle=\alpha
z_{j}^{k}+\frac{1+\alpha}{\alpha}(z_{j}^{k}-x_{j}^{k})-(1+\alpha)z_{j}^{k}+2x_{j}^{k}+\frac{1}{\alpha}(x_{j}^{k+1}-z_{j}^{k})$
$\displaystyle=\alpha
z_{j}^{k}+\frac{1+\alpha}{\alpha}(z_{j}^{k}-x_{j}^{k})+(1-\alpha)x_{j}^{k-1}+\frac{1}{\alpha}(x_{j}^{k+1}-z_{j}^{k})$
$\displaystyle=\alpha
z_{j}^{k}+(1-\alpha)v_{j}^{k}+\frac{1}{\alpha}(x_{j}^{k+1}-z_{j}^{k})=\alpha
x_{j}^{k+1}+(1-\alpha)v_{j}^{k}+\left(\frac{1}{\alpha}-\alpha\right)\left(x_{j}^{k+1}-z_{j}^{k}\right).$

As byproduct of the above derivation, we also have

$\alpha
z_{j}^{k}+(1-\alpha)v_{j}^{k}=\frac{1}{\alpha}\left(z_{j}^{k}-(1-\alpha)x_{j}^{k}\right).$

Therefore,

$\displaystyle\frac{\mu}{2}\left\lVert
v_{j}^{k+1}-x^{\star}\right\rVert^{2}=\frac{\mu}{2}\left\lVert\alpha
x_{j}^{k+1}+(1-\alpha)v_{j}^{k}-x^{\star}\right\rVert^{2}+\frac{\mu}{2}\left(\frac{1}{\alpha}-\alpha\right)^{2}\left\lVert
x_{j}^{k+1}-z_{j}^{k}\right\rVert^{2}$ (25)
$\displaystyle\qquad+\mu\left(\frac{1}{\alpha}-\alpha\right)\left\langle\alpha
x_{j}^{k+1}+(1-\alpha)v_{j}^{k}-x^{\star},x_{j}^{k+1}-z_{j}^{k}\right\rangle$
$\displaystyle=\frac{\mu}{2}\left\lVert\alpha
x_{j}^{k+1}+(1-\alpha)v_{j}^{k}-x^{\star}\right\rVert^{2}+\frac{\mu}{2}\left(\frac{1}{\alpha}-\alpha\right)^{2}\left\lVert
x_{j}^{k+1}-z_{j}^{k}\right\rVert^{2}+\mu\left(\frac{1}{\alpha}-\alpha\right)\alpha\left\lVert
x_{j}^{k+1}-z_{j}^{k}\right\rVert^{2}$
$\displaystyle\qquad+\mu\left(\frac{1}{\alpha}-\alpha\right)\left\langle\alpha
z_{j}^{k}+(1-\alpha)v_{j}^{k}-x^{\star},x_{j}^{k+1}-z_{j}^{k}\right\rangle$
$\displaystyle\leq\frac{(1-\alpha)\mu}{2}\left\lVert
v_{j}^{k}-x^{\star}\right\rVert^{2}+\frac{\alpha\mu}{2}\left\lVert
x_{j}^{k+1}-x^{\star}\right\rVert^{2}+\frac{\delta^{2}+2\mu\delta}{2(\mu+\delta)}\left\lVert
x_{j}^{k+1}-z_{j}^{k}\right\rVert^{2}$
$\displaystyle\qquad+\frac{\mu}{\alpha}\left(\frac{1}{\alpha}-\alpha\right)\left\langle
z_{j}^{k}-(1-\alpha)x_{j}^{k}-\alpha
x^{\star},x_{j}^{k+1}-z_{j}^{k}\right\rangle.$

Combining (24) and (25), we obtain the decay of the potential function $P^{k}$
defined in (19):

$\displaystyle P^{k+1}$
$\displaystyle=p^{k+1}+\frac{1}{m}\sum_{j=1}^{m}\frac{\mu}{2}\left\lVert
v_{j}^{k+1}-x^{\star}\right\rVert^{2}+e^{k,T}$
$\displaystyle\leq(1-\alpha)\left(p^{k}+\frac{1}{m}\sum_{j=1}^{m}\frac{\mu}{2}\left\lVert
v_{j}^{k}-x^{\star}\right\rVert^{2}\right)+\frac{\delta}{m}\sum_{j=1}^{m}\left(\left\lVert
z_{j}^{k}-\bar{z}^{k}\right\rVert^{2}+2\left\langle
x_{j}^{k+1}-\bar{x}^{k+1},\bar{z}^{k}-z_{j}^{k}\right\rangle\right)$
$\displaystyle\quad+\underbrace{\left(\mu-\frac{\mu+\delta}{2}-\frac{\mu^{2}}{2(\mu+\delta)}\right)}_{\leq
0}\frac{1}{m}\sum_{j=1}^{m}\left\lVert x_{j}^{k+1}-z_{j}^{k}\right\rVert^{2}$
$\displaystyle\quad+\frac{\delta}{m}\sum_{j=1}^{m}\left(\left\langle
z_{j}^{k}-\bar{z}^{k},x_{j}^{k+1}-\bar{x}^{k+1}\right\rangle+(1-\alpha)\left\langle
z_{j}^{k}-\bar{z}^{k},x_{j}^{k}-\bar{x}^{k}\right\rangle-\left\lVert
z_{j}^{k}-\bar{z}^{k}\right\rVert^{2}\right)$
$\displaystyle\quad-\frac{\mu+\delta}{m}\sum_{j=1}^{m}\left\langle\alpha
x^{\star}+(1-\alpha)x_{j}^{k}-x_{j}^{k+1},x_{j}^{k+1}-x^{k+1,\star}\right\rangle+\epsilon^{k+1}$
$\displaystyle{\leq}(1-\alpha)\left(p^{k}+\frac{1}{m}\sum_{j=1}^{m}\frac{\mu}{2}\left\lVert
v_{j}^{k}-x^{\star}\right\rVert^{2}\right)+\frac{\delta}{m}\sum_{j=1}^{m}\left((1-\alpha)\left\langle
z_{j}^{k}-\bar{z}^{k},x_{j}^{k}-\bar{x}^{k}\right\rangle-\left\langle
z_{j}^{k}-\bar{z}^{k},x_{j}^{k+1}-\bar{x}^{k+1}\right\rangle\right)$
$\displaystyle\quad-\frac{\mu+\delta}{m}\sum_{j=1}^{m}\left\langle\alpha
x^{\star}-\alpha
v_{j}^{k+1},x_{j}^{k+1}-x^{k+1,\star}\right\rangle+\epsilon^{k+1}$
$\displaystyle\leq(1-\alpha)\underbrace{\left(p^{k}+\frac{1}{m}\sum_{j=1}^{m}\frac{\mu}{2}\left\lVert
v_{j}^{k}-x^{\star}\right\rVert^{2}+e^{k-1,T}\right)}_{=P^{k}}+\frac{\delta\alpha}{m}\left\lVert
v_{\perp}^{k}\right\rVert\left\lVert
z_{\perp}^{k}\right\rVert+\alpha\frac{\mu+\delta}{m}\sum_{j=1}^{m}\left\lVert
x^{\star}-v_{j}^{k+1}\right\rVert\left\lVert
x_{j}^{k+1}-x^{k+1,\star}\right\rVert+\epsilon^{k+1}$
$\displaystyle\leq(1-\alpha)\,P^{k}+\frac{\delta}{1+\alpha}\,\frac{1}{m}\left(\left\lVert
x_{\perp}^{k+1}\right\rVert+(1-\alpha)\left\lVert
x_{\perp}^{k}\right\rVert\right)\left(2\left\lVert
x_{\perp}^{k}\right\rVert+(1-\alpha)\left\lVert
x_{\perp}^{k-1}\right\rVert\right)\vskip 6.0pt plus 2.0pt minus 2.0pt$
$\displaystyle\quad+\alpha(\mu+\delta)\sqrt{\frac{1}{m}\sum_{j=1}^{m}\left\lVert
x^{\star}-v_{j}^{k+1}\right\rVert^{2}}\sqrt{\frac{1}{m}\sum_{j=1}^{m}\left\lVert
x_{j}^{k+1}-x^{k+1,\star}\right\rVert^{2}}+\epsilon^{k+1}$
$\displaystyle\leq(1-\alpha)\,P^{k}+\frac{\delta}{1+\alpha}\,\frac{1}{m}\left(\frac{3}{2}\left\lVert
x_{\perp}^{k+1}\right\rVert^{2}+\left(\frac{7}{2}-3\alpha+\frac{\alpha^{2}}{2}\right)\left\lVert
x_{\perp}^{k}\right\rVert^{2}+(1-\alpha)^{2}\left\lVert
x_{\perp}^{k-1}\right\rVert^{2}\right)$
$\displaystyle\quad+\sqrt{2}\alpha\sqrt{\mu+\delta}\sqrt{\frac{1}{m}\sum_{j=1}^{m}\left\lVert
x^{\star}-v_{j}^{k+1}\right\rVert^{2}}\cdot\sqrt{g^{k,T}}+\epsilon^{k+1}.$

We proceed now to bound $g^{k,T}$ and the consensus error terms by
$\epsilon^{k+1}$. For the former we readily have $g^{k,T}\leq\epsilon^{k+1}$,
due to (18). For the latter, we choose $\\{\epsilon^{k}\\}$ as

$\epsilon^{k+1}\triangleq\epsilon^{k}\left(1-c\alpha\right),\quad
k=0,1,\ldots,$ (26)

where $c$ is any constant in $(0,1)$ and $\epsilon^{0}$ is to be determined;
and use $e^{k,T}\leq\epsilon^{k+1}$, for any $k=0,1\ldots$ (still due to
(18)). We can write

$\displaystyle
P^{k+1}\leq(1-\alpha)P^{k}+\sqrt{2\mu}\sqrt{\frac{1}{m}\sum_{j=1}^{m}\left\lVert
x^{\star}-v_{j}^{k+1}\right\rVert^{2}}\cdot\sqrt{\epsilon^{k+1}}+c_{1}\epsilon^{k+1},$

where

$c_{1}=1+\frac{\delta}{c_{x}}\frac{\frac{3}{2}(1-c\alpha)^{2}+5-4c\,\alpha}{(1-c\alpha)^{2}}.$

Define $\lambda^{k}\triangleq(1-\alpha)^{k};$ we have

$\displaystyle\frac{P^{k+1}}{\lambda^{k+1}}$
$\displaystyle\leq\frac{P^{k}}{\lambda^{k}}+\sqrt{2\mu}\sqrt{\frac{1}{m}\sum_{j=1}^{m}\left\lVert
x^{\star}-v_{j}^{k+1}\right\rVert^{2}}\cdot\frac{\sqrt{\epsilon^{k+1}}}{\lambda^{k+1}}+c_{1}\frac{\epsilon^{k+1}}{\lambda^{k+1}}$
$\displaystyle\leq
P^{0}+c_{1}\sum_{t=1}^{k+1}\frac{\epsilon^{t}}{\lambda^{t}}+\sum_{t=1}^{k+1}\frac{\sqrt{2\mu\epsilon^{t}}}{\lambda^{t}}\sqrt{\frac{1}{m}\sum_{j=1}^{m}\left\lVert
x^{\star}-v_{j}^{k+1}\right\rVert^{2}}.$ (27)

Introducing the following quantities:

$\hat{v}^{k}\triangleq\sqrt{\frac{\mu}{2\lambda^{k}}}\cdot\sqrt{\frac{1}{m}\sum_{j=1}^{m}\left\lVert
v_{j}^{k}-x^{\star}\right\rVert^{2}},\quad a^{k}\triangleq
2\sqrt{\frac{\epsilon^{k}}{\lambda^{k}}},\quad\text{and}\quad S^{k}\triangleq
P^{0}+c_{1}\sum_{t=1}^{k}\frac{\epsilon^{t}}{\lambda^{t}},$ (28)

(C.2) can be rewritten as

$\displaystyle(\hat{v}^{k})^{2}\leq S^{k}+\sum_{t=1}^{k}a^{t}\hat{v}^{t}.$
(29)

Using (29) we can invoke [Hongzhou et al., 2015, Lemma A.10] and conclude

$\displaystyle
S^{k}+\sum_{t=1}^{k}a^{t}\hat{v}^{t}\leq\left(\sqrt{S^{k}}+\sum_{t=1}^{k}a^{t}\right)^{2}.$
(30)

This, together with (C.2), yields

$\displaystyle
P^{k}\leq\lambda^{k}\left(\sqrt{S^{k}}+\sum_{t=1}^{k}a^{t}\right)^{2}\overset{\eqref{eq:def_Lyap_terms}}{\leq}\lambda^{k}\left(\sqrt{L^{0}}+(2+\sqrt{c_{1}})\sum_{t=1}^{k}\sqrt{\frac{\epsilon^{t}}{\lambda^{t}}}\right)^{2}.$

Choosing $\epsilon^{0}=P^{0}$, resulting in
$\epsilon^{k}=P^{0}\cdot(1-c\cdot\alpha)^{k}$ [cf. (26)], leads to the desired
result:

$\displaystyle P^{k}\leq c_{2}\,P^{0}\,(1-c\cdot\alpha)^{k+1},$ (31)

with

$c_{2}=\frac{(2+\sqrt{c_{1}})^{2}}{\left(\sqrt{\frac{1-c\cdot\alpha}{1-\alpha}}-1\right)^{2}(1-\alpha)}.$

∎

###  C.3 Step 2: SONATA and inner-loop convergence

Given Proposition 8, to complete the proof of Theorem 4 and Theorem 5, we need
to show that the error condition (18) is satisfied if SONATA runs for
$T=\tilde{\mathcal{O}}(1)$ iterations; the explicit expression of $T$ depends
on the specific setting considered for the algorithm, and it is different for
the one specified in Theorem 4 and in Theorem 5–the two cases are studied
separately in Sec. C.3.1 and Sec. C.3.2, respectively.

To control the number of calls of SONATA we need to bound $g^{k,0}+e^{k,0}$,
which is done in the following lemma.

######  Lemma 9.

Instate the setting of Proposition 8. The following holds:

$g^{k,0}+e^{k,0}\leq
P^{0}\left(2(1-c\cdot\alpha)^{k}+\max\left(\frac{1}{\mu+\delta},\,2\,c_{y}\right)\dfrac{72\,\delta^{2}}{\mu}\,c_{2}\,(1-c\cdot\alpha)^{k-1}\right).$
(32)

###### Proof.

Since

$g^{k,0}+e^{k,0}=\frac{1}{m}\sum_{i=1}^{m}\left(u_{k}(x_{i}^{k})-u_{k}^{\star}\right)+c_{x}\left\lVert
x_{\perp}^{k,0}\right\rVert^{2}+c_{y}\left\lVert
y_{\perp}^{k,0}\right\rVert^{2},$

we begin bounding
$\frac{1}{m}\sum_{i=1}^{m}\left(u_{k}(x_{i}^{k})-u_{k}^{\star}\right)$.
Postponing the proof to the end of this section, we show that the following
holds

$\dfrac{1}{m}\displaystyle{\sum}_{i=1}^{m}u_{k}(x_{i}^{k})-u_{k}^{\star}\leq
2\epsilon^{k}-2e^{k-1,T}+\dfrac{\delta^{2}}{\mu+\delta}\left\lVert\bar{z}^{k}-\bar{z}^{k-1}\right\rVert^{2}.$
(33)

Therefore, we can write

$\displaystyle g^{k,0}+e^{k,0}$
$\displaystyle\overset{\eqref{eq:bound_i}}{\leq}2\epsilon^{k}-2e^{k-1,T}+\frac{\delta^{2}}{\mu+\delta}\left\lVert\bar{z}^{k}-\bar{z}^{k-1}\right\rVert^{2}+c_{x}\left\lVert
x_{\perp}^{k}\right\rVert^{2}+2\,c_{y}\left\lVert
y_{\perp}^{k}\right\rVert^{2}+2\,c_{y}\,\delta^{2}\left\lVert(z^{k-1}-z^{k})_{\perp}\right\rVert^{2}$
$\displaystyle\leq
2\epsilon^{k}+\max\left(\frac{1}{\mu+\delta},\,2\,c_{y}\right)\delta^{2}\frac{1}{m}\sum_{i\in[m]}\left\lVert
z_{i}^{k}-z_{i}^{k-1}\right\rVert^{2}$ $\displaystyle\leq
2\epsilon^{k}+\max\left(\frac{1}{\mu+\delta},\,2\,c_{y}\right)\frac{72\,\delta^{2}}{\mu}P^{k-2}$
$\displaystyle=P^{0}\left(2(1-c\cdot\alpha)^{k}+\max\left(\frac{1}{\mu+\delta},\,2\,c_{y}\right)\frac{72\,\delta^{2}}{\mu}\,c_{2}\,(1-c\cdot\alpha)^{k-1}\right).$

It remains to show that (33) holds. We prove it following similar ideas as
those in [Hongzhou et al., 2015, Lemma B.1], with some key differences due to
the distributed setting and different potential function. Writing

$\displaystyle u_{k}(x)-u_{k-1}(x)-\left(u_{k}(z)-u_{k-1}(z)\right)$
$\displaystyle=\frac{\delta}{2m}\sum_{i=1}^{m}\left(\left\lVert
x-z_{i}^{k}\right\rVert^{2}-\left\lVert
x-z_{i}^{k-1}\right\rVert^{2}-\left\lVert
z-z_{i}^{k}\right\rVert^{2}+\left\lVert z-z_{i}^{k-1}\right\rVert^{2}\right)$
$\displaystyle=\delta\left\langle\bar{z}^{k-1}-\bar{z}^{k},x-z\right\rangle,$

we have

$\displaystyle u_{k}(x_{i}^{k})-u_{k}(x^{k+1,\star})$
$\displaystyle=u_{k-1}(x_{i}^{k})-u_{k-1}(x^{k+1,\star})+\delta\left\langle\bar{z}^{k-1}-\bar{z}^{k},x_{i}^{k}-x^{k+1,\star}\right\rangle$
$\displaystyle=u_{k-1}(x_{i}^{k})-u_{k-1}(x^{k,\star})+u_{k-1}(x^{k,\star})-u_{k-1}(x^{k+1,\star})+\delta\left\langle\bar{z}^{k-1}-\bar{z}^{k},x_{i}^{k}-x^{k+1,\star}\right\rangle$
$\displaystyle\leq
u_{k-1}(x_{i}^{k})-u_{k-1}(x^{k,\star})-\frac{\mu+\delta}{2}\left\lVert
x^{k,\star}-x^{k+1,\star}\right\rVert^{2}+\delta\left\langle\bar{z}^{k-1}-\bar{z}^{k},x_{i}^{k}-x^{k+1,\star}\right\rangle.$

Therefore,

$\displaystyle\frac{1}{m}\sum_{i=1}^{m}\left(u_{k}(x_{i}^{k})-u_{k}(x^{k+1,\star})\right)$
(34) $\displaystyle\leq\epsilon^{k}-e^{k-1,T}-\frac{\mu+\delta}{2}\left\lVert
x^{k,\star}-x^{k+1,\star}\right\rVert^{2}+\delta\left\langle\bar{z}^{k-1}-\bar{z}^{k},\bar{x}^{k}-x^{k+1,\star}\right\rangle.$

At the same time, it holds

$\displaystyle\delta\left\langle\bar{z}^{k-1}-\bar{z}^{k},x^{k,\star}-x^{k+1,\star}\right\rangle\leq\frac{\mu+\delta}{2}\left\lVert
x^{k,\star}-x^{k+1,\star}\right\rVert^{2}+\frac{\delta^{2}}{2(\mu+\delta)}\left\lVert\bar{z}^{k-1}-\bar{z}^{k}\right\rVert^{2},$
(35)
$\displaystyle\delta\left\langle\bar{z}^{k-1}-\bar{z}^{k},\bar{x}^{k}-x^{k,\star}\right\rangle\leq\frac{\mu+\delta}{2}\left\lVert\bar{x}^{k}-x^{k,\star}\right\rVert^{2}+\frac{\delta^{2}}{2(\mu+\delta)}\left\lVert\bar{z}^{k-1}-\bar{z}^{k}\right\rVert^{2}.$

Combining (34) and (35), leads to

$\displaystyle\frac{1}{m}\sum_{i=1}^{m}\left(u_{k}(x_{i}^{k})-u_{k}(x^{k+1,\star})\right)\leq\epsilon^{k}-e^{k-1,T}+\frac{\delta^{2}}{\mu+\delta}\left\lVert\bar{z}^{k-1}-\bar{z}^{k}\right\rVert^{2}+\frac{\mu+\delta}{2}\left\lVert\bar{x}^{k}-x^{k,\star}\right\rVert^{2}$
$\displaystyle\leq\epsilon^{k}-e^{k-1,T}+\frac{\delta^{2}}{\mu+\delta}\left\lVert\bar{z}^{k-1}-\bar{z}^{k}\right\rVert^{2}+\frac{1}{m}\sum_{i=1}^{m}\left(u_{k-1}(x_{i}^{k})-u_{k-1}(x^{k,\star})\right)$
$\displaystyle\leq
2\epsilon^{k}-2e^{k-1,T}+\frac{\delta^{2}}{\mu+\delta}\left\lVert\bar{z}^{k-1}-\bar{z}^{k}\right\rVert^{2}.$

∎

####  C.3.1 Proof of Theorem 4

We begin proving convergence of SONATA in the setting of Theorem 4–we refer to
such an instance of SONATA as SONATA-F. Note that convergence established in
[Sun et al., 2019] is not directly applicable here. First, there is a mismatch
between the gradient tracking initialization therein and the one used in our
setting. Second, $R$-linear convergence of the optimality gap as in [Sun et
al., 2019] is not enough to certify that the termination criterion (18) is
satisfied after a finite number of iterations $T=\tilde{\mathcal{O}}(1)$. Our
refined convergence analysis of SONATA-F is stated in Lemma 10 below.

Notice that, with $\delta=\beta-\mu,$ we have that: (i) the objective function
$u_{k}(x)$ is $\beta$-similar and $\mu_{\textnormal{u}_{k}}$-strongly convex,
with $\mu_{\textnormal{u}_{k}}=\beta$; and (ii) every $f^{k}_{i}(x)$ is
$L_{\text{mx}}=L+2\beta-\mu$ smooth. By using the surrogates (4), SONATA-F
takes advantage of similarity and achieves linear convergence rate, scaling
with $\mathcal{O}({\beta}/{\mu_{\textnormal{u}_{k}}})=\mathcal{O}(1)$.

######  Lemma 10.

Consider the minimization of $u_{k}(x)$ wherein $\delta=\beta-\mu$, running
SONATA-F (initialized as in ACC-SONATA). With
$c_{x}=\frac{8(L+2\beta-\mu)^{2}}{\beta}$, $c_{y}=\frac{4}{\beta}$ and the
network connectivity $\rho$ satisfying

$\displaystyle\rho\leq\frac{1}{4\sqrt{1785}}\frac{\beta(2\beta-\mu)}{(L+2\beta-\mu)(L+4\beta-\mu)},$
(36)

SONATA-F achieves Q-linear convergence, that is,

$\displaystyle
g^{k,t+1}+e^{k,t+1}\leq\frac{33}{34}\left(g^{k,t}+e^{k,t}\right).$ (37)

###### Proof.

The proof follows builds on some intermediate results in [Sun et al., 2019];
when recalled here, we use the same notation as defined therein.

Consider [Sun et al., 2019, Proposition 3.4], and set therein
$\epsilon_{opt}=\frac{1}{2}(2\beta-\mu)$; we get

$\sigma(1)\leq
1-\frac{1}{17}\quad\text{and}\quad\eta(1)\leq\frac{18}{17\beta}.$

Therefore,

$\displaystyle g^{k,t+1}\leq\frac{16}{17}\,g^{k,t}+\frac{9}{17}e^{k,t}.$ (38)

According to [Sun et al., 2019, Proposition 3.5], we have

$\displaystyle\left\lVert x_{\perp}^{k,t+1}\right\rVert^{2}\leq
2\rho^{2}\left\lVert
x_{\perp}^{k,t}\right\rVert^{2}+2\rho^{2}\frac{1}{m}\left\lVert
d^{k,t}\right\rVert^{2},$ $\displaystyle\left\lVert
y_{\perp}^{k,t+1}\right\rVert^{2}\leq 3\rho^{2}\left\lVert
y_{\perp}^{k,t}\right\rVert^{2}+12L_{\text{mx}}^{2}\rho^{2}\left\lVert
x_{\perp}^{k,t}\right\rVert^{2}+3L_{\text{mx}}^{2}\rho^{2}\frac{1}{m}\left\lVert
d^{k,t}\right\rVert^{2};$

which leads to

$\displaystyle 4L_{\text{mx}}^{2}\left\lVert
x_{\perp}^{k,t+1}\right\rVert^{2}+2\left\lVert
y_{\perp}^{k,t+1}\right\rVert^{2}\leq\rho^{2}\,\left(32L_{\text{mx}}^{2}\left\lVert
x_{\perp}^{k,t}\right\rVert^{2}+6\left\lVert
y_{\perp}^{k,t}\right\rVert^{2}+14L_{\text{mx}}^{2}\frac{1}{m}\left\lVert
d^{k,t}\right\rVert^{2}\right).$ (39)

[Sun et al., 2019, Proposition 3.6] becomes

$\displaystyle\frac{1}{m}\left\lVert
d^{k,t}\right\rVert^{2}\leq\frac{6}{\beta}\frac{(4\beta-\mu)^{2}+4L_{\text{mx}}^{2}}{(2\beta-\mu)^{2}}\,g^{k,t}+\frac{3}{(2\beta-\mu)^{2}}\left\lVert
y_{\perp}^{k,t}\right\rVert^{2}.$ (40)

Combining (38), (39) and (40) leads to

$\displaystyle
g^{k,t+1}+e^{k,t+1}\leq\left(\frac{16}{17}+\rho^{2}\frac{168L_{\text{mx}}^{2}}{\beta^{2}}\frac{(4\beta-\mu)^{2}+4L_{\text{mx}}^{2}}{(2\beta-\mu)^{2}}\right)g^{k,t}$
$\displaystyle\qquad\qquad\qquad\quad+\left(\frac{9}{17}+\rho^{2}\max\left(8,\,\,3+21\frac{L_{\text{mx}}^{2}}{(2\beta-\mu)^{2}}\right)\right)e^{k,t}.$

It is easy to check that the contraction of $g^{k,t}+e^{k,t}$ holds when
$\rho$ satisfies

$\displaystyle\rho\leq\min\left(\sqrt{\frac{15}{272}},\,\,\sqrt{\frac{15}{34\left(3+\frac{21L_{\text{mx}}^{2}}{(2\beta-\mu)^{2}}\right)}},\,\,\sqrt{\frac{\beta^{2}}{5712L_{\text{mx}}^{2}}\frac{(2\beta-\mu)^{2}}{(4\beta-\mu)^{2}+4L_{\text{mx}}^{2}}}\right).$

A sufficient condition for that is

$\rho\leq\frac{1}{4\sqrt{1785}}\frac{\beta(2\beta-\mu)}{(L+2\beta-\mu)(L+4\beta-\mu)}.$

∎

Invoking Lemma 9, the number of inner iterations needed for (18) to hold, for
any $k=1,2,\ldots,$ can be bounded as

$\displaystyle T$ $\displaystyle\leq\left\lceil
34\,\log\frac{g^{k,0}+e^{k,0}}{\epsilon^{k+1}}\right\rceil\overset{\eqref{termination_zero_bound}}{\leq}\left\lceil
34\,\log\frac{P^{0}\left(2(1-c\cdot\alpha)^{k}+\frac{576(\beta-\mu)^{2}}{\mu\beta}c_{2}\,(1-c\cdot\alpha)^{k-1}\right)}{P^{0}\,(1-c\cdot\alpha)^{k+1}}\right\rceil$
(41) $\displaystyle=\left\lceil
34\,\log\frac{2(1-c\cdot\alpha)+\frac{576(\beta-\mu)^{2}}{\mu\beta}c_{2}}{(1-c\cdot\alpha)^{2}}\right\rceil$
$\displaystyle=\left\lceil
34\,\log\left(\frac{2}{1-c\cdot\sqrt{\mu/\beta}}+\frac{576(\beta-\mu)^{2}}{\mu\beta(1-c\cdot\sqrt{\mu/\beta})^{2}}\frac{\left(2+\sqrt{1+\frac{(\beta-\mu)\beta}{8(L+2\beta-\mu)^{2}}\frac{\frac{3}{2}(1-c\sqrt{\mu/\beta})^{2}+5-4c\,\sqrt{\mu/\beta}}{(1-c\sqrt{\mu/\beta})^{2}}}\right)^{2}}{\left(\sqrt{\frac{1-c\cdot\sqrt{\mu/\beta}}{1-\sqrt{\mu/\beta}}}-1\right)^{2}(1-\sqrt{\mu/\beta})}\right)\right\rceil$
$\displaystyle=\mathcal{O}\left(\log\frac{\beta}{\mu}\right).$

For $k=0$, due to $g^{0,0}+e^{0,0}\leq P^{0},$ we have $T=\left\lceil
34\,\log\frac{1}{1-c\cdot\alpha}\right\rceil$, which is smaller than the one
in (41).

This, together with Proposition 8 completes the proof of Theorem
4.$\hfill\square$

####  C.3.2 Proof of Theorem 5

We study now convergence of SONATA in the setting of Theorem 5–we refer to
such an instance as SONATA-L.

We begin noticing that, with $\delta=L-\mu,$ we have the following properties
for $u_{k}(x)$ and $f^{k}_{i}(x)$’s: (i) $u_{k}(x)$ is $(2L-\mu)$-smooth,
$\beta$-similar and $L$-strongly convex; and (ii) every $f^{k}_{i}(x)$ is
$L_{\text{mx}}=2L+\beta-\mu$ smooth. Hence, when using the linearization
surrogates (5), SONATA-L achieves linear rate, scaling as
$\mathcal{O}((2L-\mu)/L)=\mathcal{O}(1)$. We establish such a result below,
following the same path as in the proof of Lemma 10.

######  Lemma 11.

Consider the minimization of $u_{k}(x)$ wherein $\delta=L-\mu$, running
SONATA-L (initialized as in ACC-SONATA). With
$c_{x}=\frac{56(2L+\beta-\mu)^{2}}{L}$, $c_{y}=\frac{28}{L}$ and the network
connectivity $\rho$ satisfying

$\displaystyle\rho\leq\frac{1}{70\sqrt{15}}\frac{L^{2}}{\left(2L-\mu+\beta\right)^{2}},$
(42)

SONATA-L achieves Q-linear convergence, that is,

$\displaystyle
g^{k,t+1}+e^{k,t+1}\leq\frac{9}{10}\left(g^{k,t}+e^{k,t}\right).$ (43)

###### Proof.

We set $\epsilon_{opt}$ in [Sun et al., 2019, Proposition 3.4 ] as
$\epsilon_{opt}=\frac{L}{2}$; one gets therein $\sigma(1)\leq\frac{4}{5}$ and
$\eta(1)\leq\frac{7}{L}$. Therefore,

$\displaystyle g^{k,t+1}\leq\frac{4}{5}g^{k,t}+\frac{1}{2}e^{k,t}.$ (44)

In addition, [Sun et al., 2019, Proposition 3.6] becomes

$\displaystyle\frac{1}{m}\left\lVert
d^{k,t}\right\rVert^{2}\leq\frac{6}{L}\left(\frac{9}{4}+\frac{4L_{\text{mx}}^{2}}{(2L-\mu)^{2}}\right)g^{k,t}+\frac{3}{(2L-\mu)^{2}}\left\lVert
y_{\perp}^{k,t}\right\rVert^{2}.$ (45)

Combining (39), (44) and (45), leads to

$\displaystyle
g^{k,t+1}+e^{k,t+1}\leq\left(\frac{4}{5}+1176\,\rho^{2}\frac{L_{\text{mx}}^{2}}{L^{2}}\left(\frac{9}{4}+\frac{4L_{\text{mx}}^{2}}{(2L-\mu)^{2}}\right)\right)g^{k,t}+\left(\frac{1}{2}+\rho^{2}\max\left(8,\,\,3+21\frac{L_{\text{mx}}^{2}}{(2L-\mu)^{2}}\right)\right)e^{k,t}$

A contraction on $g^{k,t}+e^{k,t}$ is ensured choosing

$\displaystyle\rho\leq\min\left(\frac{1}{2\sqrt{5}},\,\,\sqrt{\frac{2}{5\left(3+21\frac{L_{\text{mx}}^{2}}{(2L-\mu)^{2}}\right)}},\,\,\frac{1}{28}\sqrt{\frac{1}{15\frac{L_{\text{mx}}^{2}}{L^{2}}\left(\frac{9}{4}+\frac{4L_{\text{mx}}^{2}}{(2L-\mu)^{2}}\right)}}\right).$

The above is satisfied if we set

$\rho\leq\frac{1}{70\sqrt{15}}\frac{L^{2}}{\left(2L-\mu+\beta\right)^{2}}.$

∎

We can conclude the proof of Theorem 4, using Lemma 9 to determine the number
of inner iterations needed for (18) to hold: for any $k=1,2,\ldots,$ we have

$\displaystyle T$ $\displaystyle\leq\left\lceil
10\,\log\frac{g^{k,0}+e^{k,0}}{\epsilon^{k+1}}\right\rceil\leq\left\lceil
10\log\frac{2(1-c\cdot\alpha)+\frac{4032(L-\mu)^{2}}{L\mu}c_{2}}{(1-c\cdot\alpha)^{2}}\right\rceil=\mathcal{O}\left(\log\kappa\right).$
(46)

For $k=0$, due to $g^{0,0}+e^{0,0}\leq P^{0},$ we have $T=\left\lceil
10\,\log\frac{1}{1-c\cdot\alpha}\right\rceil$, which is smaller than the value
of $T$ in (46) for $k\geq 1$.

##  Appendix D ACC-SSONATA over star-networks

In this section we customize ACC-SSONATA to the special setting of
master/workers architectures. Specifically, consider Problem (P) over 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$ in (P). The application of ACC-SSONATA to
such a setting boils down to customize the inner algorithm SONATA–we use the
instance in [Sun et al., 2019, Algorithm 3], which is reported in Algorithm 3
below (applied to (P)) for convenience. Note that the consensus step is now
replaced by the exact average performed by the master node (see (S.4)), based
upon reception of the local optimization variables from the workers. Also,
there is no need of the gradient-tracking mechanism, as the master node can
directly broadcast to the workers the aggregate gradient $\nabla f(x^{k})$.
Notice that SONATA-star can be interpreted as a special instance of SONATA
described in Algorithm 3 (up to a proper initialization) if the weight matrix
$W$ therein is set to $W=[1,0_{m,m-1}][1/m,0_{m,m-1}]^{\top}$, whose
associated $\rho=\|W-11^{\top}/m\|=0$.

Equipped with SONATA-star, ACC-SONATA-star reduces to Algorithm 4 below.
Convergence as discussed in the main text is a consequence of Theorem 4 and
Theorem 5. Note that the condition of $\rho$ is trivially satisfied as
$\rho=0$.

Algorithm 3 SONATA-
star$(\left\\{f_{i}\right\\}_{i\in[m]},\,(x_{i}^{0})_{i\in[m]},T)$

Input: $\left\\{f_{i}(x)\right\\}_{i\in[m]},\,\,r(x)$ [cf. (P)];

$(x_{i}^{0})_{i\in[m]}$ [init. points],

$T>0$ [# iterations];
Output: $x^{T}$;

for $k=0,1,2,\ldots,T-1$ do

(S.1): Each worker $i$ evaluates $\nabla f_{i}(x^{k})$ and sends it to the
master node;

(S.2): The master broadcasts $\nabla f(x^{k})=1/m\sum_{i=1}^{m}\nabla
f_{i}(x^{k})$ to the workers;

(S.3): Each worker $i$ computes

$x_{i}^{k+1/2}\triangleq{}\mathop{\rm
argmin}_{x_{i}\in\mathbb{R}^{d}}{\tilde{f}_{i}(x_{i};x^{k})+\big{(}\nabla
f(x^{k})-\nabla f_{i}(x^{k})\big{)}^{\top}(x_{i}-x^{k})}+r(x_{i}),$

and sends $x_{i}^{k+1/2}$ to the master;

(S.4): The master computes the average

$x^{k+1}=\frac{1}{m}\sum_{i=1}^{m}x_{i}^{k+1/2},\vspace{-0.2cm}$

and sends it back to the workers.

end for

Algorithm 4 Accelerated SONATA-star

Input: $\beta$, $\mu$, $\delta>0$, $\alpha=\sqrt{\mu/(\mu+\delta)}$;

$x^{0}=z^{0}=0$;
Output: $x^{K}$
for $k=0,1,2,\ldots,K-1$ do

Set: $f_{i}^{k}(x)=f_{i}(x)+\frac{\delta}{2}\left\lVert
x-z^{k}\right\rVert^{2},$

(S.1) Inner loop via SONATA-star:

$\displaystyle x^{k+1}=\texttt{SONATA-
star}\Big{(}\left\\{f_{i}^{k}\right\\}_{i\in[m]},\,x^{k},\,T\Big{)};$

(S.2) Extrapolation step:

$z^{k+1}=x^{k+1}+\frac{1-\alpha}{1+\alpha}\,(x^{k+1}-x^{k}).$

end for