# ADMM-Tracking Gradient
for Distributed Optimization over
Asynchronous and Unreliable Networks

Guido Carnevale, Nicola Bastianello, Giuseppe Notarstefano, Ruggero Carli This
work was supported in part by the Italian Ministry of Foreign Affairs and
International Cooperation, grant number BR22GR01.Guido Carnevale and Giuseppe
Notarstefano are with the Department of Electrical, Electronic and Information
Engineering, Alma Mater Studiorum - Universita‘ di Bologna, Bologna, Italy.
Emails<EMAIL_ADDRESS>Bastianello is with the School of
Electrical Engineering and Computer Science, KTH Royal Institute of
Technology, Stockholm, Sweden. Email<EMAIL_ADDRESS>Carli is with the
Department of Information Engineering of the University of Padova, Via G.
Gradenigo 6/B, 35131 Padova, Italy. Email<EMAIL_ADDRESS>
###### Abstract

In this paper, we propose (i) a novel distributed algorithm for consensus
optimization over networks and (ii) a robust extension tailored to deal with
asynchronous agents and packet losses. The key idea is to achieve dynamic
consensus on (i) the agents’ average and (ii) the global descent direction by
iteratively solving an online auxiliary optimization problem through a
distributed implementation of the Alternating Direction Method of Multipliers
(ADMM). Such a mechanism is suitably interlaced with a local proportional
action steering each agent estimate to the solution of the original consensus
optimization problem. First, in the case of ideal networks, by using tools
from system theory, we prove the linear convergence of the scheme with
strongly convex costs. Then, by exploiting the averaging theory, we extend
such a first result to prove that the robust extension of our method preserves
linear convergence in the case of asynchronous agents and packet losses.
Further, by using the notion of Input-to-State Stability, we also guarantee
the robustness of the schemes with respect to additional, generic errors
affecting the agents’ updates. Finally, some numerical simulations confirm our
theoretical findings and show that the proposed methods outperform the
existing state-of-the-art distributed methods for consensus optimization.

##  I Introduction

Multi-agent systems have become ubiquitous in many engineering domains,
ranging from robotics to power grids, from traffic networks to the Internet of
Things. Recent technological advances have indeed led to the widespread
adoption of devices equipped with computational and communication resources.
These devices (or agents) then form interconnected systems capable of
cooperatively performing a variety of tasks, e.g. learning, resource
allocation, and exploration. Many of these tasks can be formulated as the
consensus (or cost-coupled) optimization problem

$\displaystyle\min_{x\in\mathbb{R}^{n}}\sum_{i=1}^{N}f_{i}(x),$

with $f_{i}$ the local costs of the $N$ agents. This motivated a wide interest
in algorithms to enable multi-agent systems to collaborate towards a solution
of consensus problems [1, 2, 3].

Different approaches to designing distributed algorithms for consensus
optimization have been explored in the literature, with the main ones being
decentralized (sub)gradient, Gradient Tracking (GT), and alternating direction
method of multipliers (ADMM) [1]. The class of decentralized (sub)gradient
methods has the drawback that exact convergence is only achieved using
diminishing stepsizes [2]. Gradient Tracking algorithms were then introduced
to guarantee exact convergence with fixed step-sizes, see [4, 5, 6], and have
proved their versatility in different setups, see e.g. [7, 8, 9, 10, 11, 12,
13]. Besides these gradient-based methods, ADMM, popularized by [14], has
received wide attention and has proven especially suited to distributed setups
[15, 16, 17, 18].

However, the successful deployment of distributed algorithms designed
according to these approaches is hindered by the practical challenges that
multi-agent systems raise. First of all, the agents cooperating may have
highly heterogeneous computational resources, implying that they may perform
local updates at different rates. One option is to enforce synchronization of
the agents so that they update in lockstep. But besides the technical
challenge of synchronization, it would result in faster agents sitting idle
before they can initiate a new update [19]. For these reasons, we are
interested in asynchronous distributed algorithms. A second challenge is that
of unreliable communications between the agents, which may be subject to
packet losses or quantization errors, especially when wireless transmission is
employed [20]. It is therefore important to characterize the resilience of
distributed algorithms to imperfect communications. Finally, in applications
such as learning, the large size of the data set stored by an agent implies
that evaluating the local gradient may be too computationally expensive. For
this reason, stochastic gradients are often employed to approximate the local
gradient over a subset of the data, which introduces additive errors [21, 22,
23].

Let us now turn back to the distributed algorithms discussed above; among
them, ADMM has arguably proven to be the most robust. In particular, it
guarantees exact, almost sure convergence in the presence of asynchrony and
packet losses [16, 18, 24, 25], and can reach an approximate solution in the
presence of additive communication errors [26, 27, 28, 29]. The drawback of
ADMM is that the agents need to solve a local minimization problem to update
their state, which in general is more computationally expensive than the
gradient evaluations of GT. Gradient Tracking methods also exhibit convergence
in the presence of asynchrony and packet losses [7, 30, 8]. However, as
highlighted in [31, 32, 33, 34, 35], they are not robust to additive errors
e.g. due to unreliable communications. Indeed, GT algorithms are designed to
reconstruct the global gradient $\nabla\sum_{i=1}^{N}f_{i}$ using local
gradients $\nabla f_{i}$, which they accomplish employing a dynamic average
consensus block (see [36, 37]). The average consensus, though, introduces a
marginally stable dynamics [34], which causes the lack of robustness to
errors.

The objective of this paper then is to propose a novel gradient tracking
algorithm that is provably robust to asynchrony and unreliable communications.
Our starting point is the insight that the dynamic average consensus block in
GT algorithms can be replaced by a different, robust consensus protocol, as
done with the ratio (or push-sum) consensus in [38, 30]. Specifically, our
algorithm uses an ADMM-based dynamic consensus protocol, which was shown to be
robust to both asynchrony and unreliable communications [29], differently from
ratio consensus. The protocol is derived by formulating dynamic consensus as
an online, quadratic problem and applying the robust distributed ADMM to it
[29]. Since the consensus problem is quadratic, the agents’ updates have a
closed form and no minimization needs to be solved, avoiding the larger
computational footprint of ADMM for general problems. Our contribution then
brings together the best of both Gradient Tracking, with its lightweight
computations, and ADMM, with its robustness.

In the following, we summarize our contributions. We propose a novel
distributed gradient tracking method that employs an ADMM-based dynamic
consensus protocol. The algorithm then inherits the robustness properties of
ADMM to both asynchrony and unreliable communications. We analyze the
convergence of the proposed algorithm – for synchronous and reliable networks
– by interpreting it as a singularly perturbed system given by the
interconnection between (i) a slow subsystem given by an approximation of the
gradient descent and a consensus step, and (ii) a fast one related to the
states of the ADMM-based consensus protocol. With these tools, we are able to
show linear convergence to the optimal solution of the consensus optimization
problem. Then, this result allows us to analyze a robust version of the
algorithm tailored for networks with asynchrony and packet losses. In this
stochastic setup, we apply averaging theory to the singularly perturbed system
that represents our algorithm, and by studying the associated averaged system,
we are able to show that the linear convergence of the scheme is preserved.
Moreover, we are further able to prove that the proposed schemes exhibit
Input-to-State Stability (ISS) with respect to additive errors, due, e.g., to
unreliable communications or the use of stochastic gradients. Finally, we
showcase the performance of the algorithms in numerical simulations for
quadratic and logistic regression problems over random networks. We compare
its performance with a standard and a robust version of the GT, which lack the
same robustness properties. Preliminary results related to this paper appeared
in [39]. However, the preliminary work [39] does not consider the scenario
over asynchronous and unreliable networks and does not analyze the algorithm
robustness. Moreover, in [39], the convergence proof in the case of ideal
networks is omitted.

The paper is organized as follows. Section II introduces the problem setup. In
Section III, we design the novel distributed algorithm and provide its
convergence properties. Section IV is devoted to analyzing the algorithm. In
Section V, we present a robust version of the algorithm to address the
scenario with imperfect networks, while Section VI is devoted to analyzing
this algorithm. Section VII characterizes the robustness of the proposed
algorithms. Finally, Section VIII compares our methods and Gradient Tracking.

#### Notation

The identity matrix in $\mathbb{R}^{m\times m}$ is $I_{m}$. , while $0_{m}$ is
the all-zero matrix in $\mathbb{R}^{m\times m}$. The vector of $N$ ones is
denoted by $1_{N}$, while $\mathbf{1}_{N,n}:=1_{N}\otimes I_{n}$ with
$\otimes$ being the Kronecker product. For a finite set $S$, we denote by
$|S|$ its cardinality. The symbol $\textsc{col}(v_{i})_{i\in\\{1,\ldots,N\\}}$
denotes the vertical concatenation of the column vectors $v_{1},\dots,v_{N}$.
We denote as
$\text{blkdiag}(M_{1},\dots,M_{N})\in\mathbb{R}^{\sum_{i=1}^{N}n_{i}}$ the
block diagonal matrix whose $i$-th block is given by
$M_{i}\in\mathbb{R}^{n_{i}\times n_{i}}$.

##  II Problem Description and Preliminaries

###  II-A Problem Setup

We consider a network of $N$ agents that aim to solve

$\displaystyle\min_{x\in\mathbb{R}^{n}}\sum_{i=1}^{N}f_{i}(x),$ (1)

where $x\in\mathbb{R}^{n}$ is the (common) decision variable, while
$f_{i}:\mathbb{R}^{n}\to\mathbb{R}$ is the objective function of agent
$i\in\\{1,\ldots,N\\}$. In the following, we will also use the function
$f:\mathbb{R}^{n}\to\mathbb{R}$ defined as $f(x):=\sum_{i=1}^{N}f_{i}(x)$. Our
goal is to design an algorithm to solve (1) in a distributed way, namely with
update laws implementable over a network of agents using only (i) local
information and (ii) neighboring communication. Indeed, we consider a network
of agents communicating according to an undirected graph
$\mathcal{G}=(\\{1,\ldots,N\\},\mathcal{E})$, with
$\mathcal{E}\subset\\{1,\ldots,N\\}\times\\{1,\ldots,N\\}$ such that $i$ and
$j$ can exchange information only if $(i,j)\in\mathcal{E}$. The set of
neighbors of $i$ is
$\mathcal{N}_{i}:=\\{j\in\\{1,\ldots,N\\}\mid(i,j)\in\mathcal{E}\\}$, while
its degree is $d_{i}:=|\mathcal{N}_{i}|$ and
$\mathbf{d}:=\sum_{i=1}^{N}d_{i}$. Notice that it holds
$\mathbf{d}=2|\mathcal{E}|$.

The following assumptions formalize the considered setup.

######  Assumption II.1.

[Graph] The graph $\mathcal{G}$ is connected. $\square$

######  Assumption II.2.

[Objective functions] The objective function $f$ is $\mathbf{c}$-strongly
convex, while the gradients $\nabla f_{i}(\cdot)$ are $L$-Lipschitz continuous
for all $i\in\\{1,\ldots,N\\}$. $\square$

We notice that Assumption II.2 ensures that problem (1) has a unique minimizer
and we denote it as $x^{\star}\in\mathbb{R}^{n}$.

###  II-B Fundamentals of singularly perturbed systems

In the next, we will propose two novel distributed algorithms and, to assess
their convergence property, we will resort to a system-theoretic perspective
based on _singular perturbation_ theory. To this end, we provide below a
generic convergence result extending [40, Theorem II.5] for the time-varying
case. Indeed, we will also consider the case in which the agents are
asynchronous and their communication is subject to packet losses and, thus,
the result [40, Theorem II.5] does not suit our purposes since it only
considers time-invariant dynamics.

######  Theorem II.3 (Global exponential stability for time-varying singularly
perturbed systems).

Consider the system

$\displaystyle x^{{{t}+1}}$ $\displaystyle=x^{{t}}+\gamma
f(x^{{t}},z^{{t}},{t})$ (2a) $\displaystyle z^{{{t}+1}}$
$\displaystyle=g(z^{{t}},x^{{t}},{t}),$ (2b)

with $x^{{t}}\in\mathcal{D}\subseteq\mathbb{R}^{n}$,
$z^{{t}}\in\mathbb{R}^{m}$,
$f:\mathcal{D}\times\mathbb{R}^{m}\times\mathbb{N}\rightarrow\mathbb{R}^{n}$,
$g:\mathbb{R}^{m}\times\mathbb{R}^{n}\times\mathbb{N}\rightarrow\mathbb{R}^{m}$,
$\gamma>0$. Assume that $f$ and $g$ are Lipschitz continuous uniformly in
${t}$ with respect to $x^{{t}}$ and $z^{{t}}$ with Lipschitz constants
$L_{f},L_{g}>0$, respectively. Assume that there exists
$z^{\text{eq}}:\mathbb{R}^{n}\to\mathbb{R}^{m}$ such that for any
$x\in\mathbb{R}^{n}$

$\displaystyle 0$ $\displaystyle=\gamma f(0,z^{\text{eq}}(0),{t})$
$\displaystyle z^{\text{eq}}(x)$ $\displaystyle=g(z^{\text{eq}}(x),x,{t}),$

for all ${t}\in\mathbb{N}$ and with $z^{\text{eq}}$ being Lipschitz continuous
with Lipschitz constant $L_{z^{\text{eq}}}>0$. Let

$x^{{{t}+1}}=x^{{t}}+\gamma f(x^{{t}},z^{\text{eq}}(x^{{t}}),{t})$ (3)

be the _reduced system_ and

$\tilde{z}^{{t}+1}=g(\tilde{z}^{{t}}+z^{\text{eq}}(x),x,{t})-z^{\text{eq}}(x)$
(4)

be the _boundary layer system_ with $\tilde{z}^{{t}}\in\mathbb{R}^{m}$.

Assume that there exist a continuous function
$U:\mathbb{R}^{m}\times\mathbb{N}\to\mathbb{R}$ and $\bar{\gamma}_{1}>0$ such
that, for any $\gamma\in(0,\bar{\gamma}_{1})$, there exist
$b_{1},b_{2},b_{3},b_{4}>0$ such that

$\displaystyle b_{1}\left\|\tilde{z}\right\|^{2}\leq U(\tilde{z},{t})\leq
b_{2}\left\|\tilde{z}\right\|^{2}$ (5a) $\displaystyle
U(g(\tilde{z}+z^{\text{eq}}(x),x,{t})-z^{\text{eq}}(x),{{t}+1})-U(\tilde{z},{t})$
$\displaystyle\leq-b_{3}\left\|\tilde{z}\right\|^{2}$ (5b)
$\displaystyle|U(\tilde{z}_{1},{t})-U(\tilde{z}_{2},{t})|\leq
b_{4}\left\|\tilde{z}_{1}-\tilde{z}_{2}\right\|\left\|\tilde{z}_{1}\right\|$
$\displaystyle\hskip
14.22636pt+b_{4}\left\|\tilde{z}_{1}-\tilde{z}_{2}\right\|\left\|\tilde{z}_{2}\right\|,$
(5c)

for any $\tilde{z},\tilde{z}_{1},\tilde{z}_{2}\in\mathbb{R}^{m}$,
$x\in\mathbb{R}^{n}$, ${t}\in N$. Further, assume there exist a continuous
function $W:\mathcal{D}\times\mathbb{N}\to\mathbb{R}$ and $\bar{\gamma}_{2}>0$
such that, for any $\gamma\in(0,\bar{\gamma}_{2})$, there exist
$c_{1},c_{2},c_{3},c_{4}>0$ such that

$\displaystyle c_{1}\left\|x\right\|^{2}\leq W(x,{t})\leq
c_{2}\left\|x\right\|^{2}$ (6a) $\displaystyle W(x+\gamma
f(x,z^{\text{eq}}(x),{t}),{t})-W(x,{t})\leq-c_{3}\gamma\left\|x\right\|^{2}$
(6b) $\displaystyle|W(x_{1},{t})-W(x_{2},{t})|\leq
c_{4}\left\|x_{1}-x_{2}\right\|\left\|x_{1}\right\|$ $\displaystyle\hskip
14.22636pt+c_{4}\left\|x_{1}-x_{2}\right\|\left\|x_{2}\right\|,$ (6c)

for any $x,x_{1},x_{2},x_{3}\in\mathcal{D}$, ${t}\in\mathbb{N}$.

Then, there exist
$\bar{\gamma}\in(0,\min\\{\bar{\gamma}_{1},\bar{\gamma}_{2}\\})$,
$\kappa_{1}>0$, and $\kappa_{2}>0$ such that, for all
$\gamma\in(0,\bar{\gamma})$, it holds

$\displaystyle\left\|\begin{bmatrix}x^{{t}}-x^{\star}\\\
z^{{t}}-z^{\text{eq}}(x^{{t}})\end{bmatrix}\right\|\leq\kappa_{1}\left\|\begin{bmatrix}x^{0}\\\
z^{0}-z^{\text{eq}}(x^{0})\end{bmatrix}\right\|e^{-\kappa_{2}t},$

for any $(x^{0},z^{0})\in\mathcal{D}\times\mathbb{R}^{m}$.

The proof of Theorem II.3 is given in Appendix -B.

##  III ADMM-Tracking Gradient: Algorithm Design and Convergence Properties

Let $x_{i}^{{t}}\in\mathbb{R}^{n}$ be the estimate of the solution to problem
(1) maintained by agent $i$ at iteration ${t}\in\mathbb{N}$. We follow a
control-oriented design for the update of $x_{i}^{{t}}$. Then, let
$u_{i}^{t}\in\mathbb{R}^{n}$ be the $i$-th control input and consider the
single integrator dynamics

$\displaystyle x_{i}^{{t}+1}=x_{i}^{{t}}+u_{i}^{t}.$ (7)

The control law determining $u_{i}^{t}$ should have the twofold purpose of
removing (i) the consensus error with respect to the other agents’ estimates,
and (ii) the optimality error related to problem (1). Hence, one may design
$u_{i}^{t}$ as a proportional action with respect to the above errors, namely

$\displaystyle
u_{i}^{t}=\left(\frac{1}{N}\sum_{j=1}^{N}x_{j}^{{t}}-x_{i}^{{t}}\right)-\frac{\delta}{N}\sum_{j=1}^{N}\nabla
f_{j}(x_{j}^{{t}}),$ (8)

where $\delta>0$ is the stepsize. However, in a distributed setting, agent $i$
cannot access $\frac{1}{N}\sum_{j=1}^{N}x_{j}^{{t}}$ and
$\frac{1}{N}\sum_{j=1}^{N}\nabla f_{j}(x_{j}^{{t}})$. Therefore, we modify the
control law (8) by employing two auxiliary variables
$y_{i}^{{t}},s_{i}^{{t}}\in\mathbb{R}^{n}$ aimed at reconstructing
$\frac{1}{N}\sum_{j=1}^{N}x_{j}^{{t}}$ and $\frac{1}{N}\sum_{j=1}^{N}\nabla
f_{j}(x_{j}^{{t}})$, respectively, obtaining

$\displaystyle u_{i}^{t}=(y_{i}^{{t}}-x_{i}^{{t}})-s_{i}^{{t}}.$ (9)

We note that, if $y_{i}^{{t}}=\frac{1}{N}\sum_{j=1}^{N}x_{j}^{{t}}$ and
$s_{i}^{{t}}=\frac{1}{N}\sum_{j=1}^{N}\nabla f_{j}(x_{j}^{{t}})$, then the
desired update (8) is recovered. Then, inspired by [29], for each iteration
${t}\geq 0$, we turn these two consensus-oriented goals into the online
optimization problem

$\displaystyle\begin{split}\min_{\begin{subarray}{c}(y_{1},\dots,y_{N})\in\mathbb{R}^{Nn}\\\
(s_{1},\dots,s_{N})\in\mathbb{R}^{Nn}\end{subarray}}&\sum_{i=1}^{N}\ell_{i}^{t}(y_{i},s_{i})\\\
\text{s.t.:}\>\>&\begin{bmatrix}y_{i}\\\
s_{i}\end{bmatrix}=\begin{bmatrix}y_{j}\\\
s_{j}\end{bmatrix}\>\forall(i,j)\in\mathcal{E},\end{split}$ (10)

where, for all $i\in\\{1,\ldots,N\\}$,
$\ell_{i}^{t}:\mathbb{R}^{n}\times\mathbb{R}^{n}\to\mathbb{R}$ reads as

$\displaystyle\ell_{i}^{t}(y_{i},s_{i})=\frac{1}{2}\left\|y_{i}-x_{i}^{{t}}\right\|^{2}+\frac{1}{2}\left\|s_{i}-\nabla
f_{i}(x_{i}^{{t}})\right\|^{2}.$

Indeed, if the graph $\mathcal{G}$ is connected, then the (unique) optimal
solution of problem (10), say it
$(y_{\star}^{t},s_{\star}^{t})\in\mathbb{R}^{2Nn}$, reads as
$(y_{\star}^{t},s_{\star}^{t})=(\mathbf{1}_{N,n}\frac{1}{N}\sum_{j=1}^{N}x_{j}^{{t}},\mathbf{1}_{N,n}\frac{1}{N}\sum_{j=1}^{N}\nabla
f_{j}(x_{j}^{{t}}))$ [29]. From this observation, we design the updates of
$y_{i}^{{t}}$ and $s_{i}^{{t}}$ by resorting to the distributed ADMM proposed
in [18] (see Appendix -A for a description of this method in a generic
framework). Hence, each agent $i$ maintains an additional variable
$z_{ij}^{{t}}\in\mathbb{R}^{2n}$ for each neighbor $j\in\mathcal{N}_{i}$ and
implements

$\displaystyle\begin{bmatrix}y_{i}^{{t}}\\\ s_{i}^{{t}}\end{bmatrix}$
$\displaystyle=\operatorname*{arg\,min}_{\begin{subarray}{c}y_{i}\in\mathbb{R}^{n}\\\
s_{i}\in\mathbb{R}^{n}\end{subarray}}\Bigg{\\{}\ell_{i}^{t}(y_{i},s_{i})-\begin{bmatrix}y_{i}^{\top}&s_{i}^{\top}\end{bmatrix}\sum_{j\in\mathcal{N}_{i}}z_{ij}^{{t}}$
$\displaystyle\hskip 56.9055pt+\frac{\rho
d_{i}}{2}\left(\left\|y_{i}\right\|^{2}+\left\|s_{i}\right\|^{2}\right)\bigg{\\}}$
$\displaystyle z_{ij}^{{t}+1}$
$\displaystyle=(1-\alpha)z_{ij}^{{t}}+\alpha\left(-z_{ji}^{{t}}+2\rho\begin{bmatrix}y_{j}^{{t}}\\\
s_{j}^{{t}}\end{bmatrix}\right),$

with $\rho>0$ and $\alpha\in(0,1)$, Being $\ell_{i}^{t}$ quadratic, the above
updates are equivalent to the closed form

$\displaystyle\begin{bmatrix}y_{i}^{{t}}\\\ s_{i}^{{t}}\end{bmatrix}$
$\displaystyle=\frac{1}{1+\rho d_{i}}\left(\begin{bmatrix}x_{i}^{{t}}\\\
\nabla
f_{i}(x_{i}^{{t}})\end{bmatrix}+\sum_{j\in\mathcal{N}_{i}}z_{ij}^{{t}}\right)$
(11a) $\displaystyle z_{ij}^{{t}+1}$
$\displaystyle=(1-\alpha)z_{ij}^{{t}}+\alpha m_{ji}^{t},$ (11b)

in which we also introduced $m_{ji}^{t}\in\mathbb{R}^{2n}$ to denote the
message from agent $j$ needed by agent $i$ to perform (11b), namely

$\displaystyle m_{ji}^{t}:=-z_{ji}^{{t}}+2\rho\begin{bmatrix}y_{j}^{{t}}\\\
s_{j}^{{t}}\end{bmatrix},$ (12)

for all $i\in\\{1,\ldots,N\\}$ and $j\in\mathcal{N}_{i}$. Moreover, by
resorting to the concepts of singular perturbations, we modify (7) adding a
tuning gain $\gamma>0$ as follows

$\displaystyle x_{i}^{{t}+1}=x_{i}^{{t}}+\gamma u_{i}^{t}.$ (13)

The parameter $\gamma$ allows us to arbitrarily reduce the speed of variations
of $x_{i}^{{t}}$ and, as we will detail in the analysis, allows for a proper
interconnection with (11). The whole distributed protocol arising by combining
(9), (13), and (11) is reported in Algorithm 1 and we name it ADMM-Tracking
Gradient.

Initialization: $x_{i}^{0}\in\mathbb{R}^{n},z_{i}^{0}\in\mathbb{R}^{2nd_{i}}$.

for ${t}=0,1,\dots$ do

$\begin{bmatrix}y_{i}^{{t}}\\\ s_{i}^{{t}}\end{bmatrix}=\frac{1}{1+\rho
d_{i}}\left(\begin{bmatrix}x_{i}^{{t}}\\\ \nabla
f_{i}(x_{i}^{{t}})\end{bmatrix}+\sum_{j\in\mathcal{N}_{i}}z_{ij}^{{t}}\right)$

$x_{i}^{{t}+1}=x_{i}^{{t}}+\gamma(y_{i}^{{t}}-x_{i}^{{t}})-\gamma\delta
s_{i}^{{t}}$

for $j\in\mathcal{N}_{i}$ do

$m_{ij}^{t}=-z_{ij}^{{t}}+2\rho\begin{bmatrix}y_{i}^{{t}}\\\
s_{i}^{{t}}\end{bmatrix}$

transmit $m_{ij}^{t}$ to $j$ and receive $m_{ji}^{t}$ to $j$

$z_{ij}^{{t}+1}=(1-\alpha)z_{ij}^{{t}}+\alpha m_{ji}^{t}$

Algorithm 1 ADMM-Tracking Gradient (Agent $i$)

We remark that Algorithm 1 can be implemented in a fully-distributed fashion
since it only requires neighboring communication and local variables. The next
theorem states the convergence features of Algorithm 1.

######  Theorem III.1.

Consider Algorithm 1 and let Assumptions II.1 and II.2 hold. Then, there exist
$\bar{\gamma},\bar{\delta},c_{1},c_{2}>0$ such that, for any $\rho>0$,
$\alpha\in(0,1)$, $\gamma\in(0,\bar{\gamma})$, $\delta\in(0,\bar{\gamma})$,
$(x_{i}^{0},z_{i}^{0})\in\mathbb{R}^{n}\times\mathbb{R}^{2nd_{i}}$, for all
$i\in\\{1,\ldots,N\\}$, it holds

$\displaystyle\left\|x_{i}^{{t}}-x^{\star}\right\|\leq c_{1}\exp(-c_{2}{t}).$
$\square$

Theorem III.1 is proved in Section IV-D. In detail, with tools from system
theory, we will show global exponential stability of
$(\mathbf{1}_{N,n}x^{\star},z^{\star})$ for the aggregate form of ADMM-
Tracking Gradient, where the quantity $z^{\star}\in\mathbb{R}^{2n\mathbf{d}}$
will be specified later.

##  IV Convergence Analysis of ADMM-Tracking Gradient

Here, we provide the analysis of Algorithm 1. First, in Section IV-A, we
interpret the aggregate form of ADMM-Tracking Gradient as a _singularly
perturbed_ system. Then, in Section IV-B and IV-C, respectively, we separately
study the identified subsystems by relying on Lyapunov theory. Finally, in
Section IV-D, we use the results collected in the previous steps to prove
Theorem III.1.

###  IV-A ADMM-Tracking Gradient as a Singularly Perturbed System

We start by providing a more compact form for the local update of ADMM-
Tracking Gradient. To this end, let $z_{i}^{{t}}$ be the vector stacking all
the variables $z_{ij}^{{t}}$ of the agent $i$, i.e.,
$z_{i}^{{t}}:=\textsc{col}(z_{ij}^{{t}})_{j\in\mathcal{N}_{i}}\in\mathbb{R}^{2nd_{i}}$,
while let
$z_{\mathcal{N}_{i}}^{{t}}:=\textsc{col}(z_{ji}^{{t}})_{j\in\mathcal{N}_{i}}\in\mathbb{R}^{2nd_{i}}$,
i.e., the vector stacking all the variables $z_{ji}^{{t}}$ with
$j\in\mathcal{N}_{i}$. Let
$h_{i}^{x}:\mathbb{R}^{n}\times\mathbb{R}^{nd_{i}}\to\mathbb{R}^{n}$,
$h_{i}^{\nabla}:\mathbb{R}^{n}\times\mathbb{R}^{nd_{i}}\to\mathbb{R}^{n}$, and
$h_{\mathcal{N}_{i}}:\mathbb{R}^{nd_{i}}\times\mathbb{R}^{nd_{i}}\to\mathbb{R}^{2nd_{i}}$
be the functions described by

$\displaystyle h_{i}^{x}(x_{i},z_{i})$
$\displaystyle=\frac{x_{i}+\begin{bmatrix}I_{n}&0_{n}\end{bmatrix}\sum_{j\in\mathcal{N}_{i}}z_{ij}}{1+\rho
d_{i}}$ $\displaystyle h_{i}^{\nabla}(x_{i},z_{i})$
$\displaystyle=\frac{\nabla
f_{i}(x_{i})+\begin{bmatrix}0_{n}&I_{n}\end{bmatrix}\sum_{j\in\mathcal{N}_{i}}z_{ij}}{1+\rho
d_{i}}$ (14) $\displaystyle
h_{\mathcal{N}_{i}}(x_{\mathcal{N}_{i}},z_{\mathcal{N}_{i}})$
$\displaystyle=\textsc{col}\\!\left(\\!\frac{\textsc{col}(x_{j}^{{t}},\nabla
f_{j}(x_{j}^{{t}}))+\sum_{j\in\mathcal{N}_{i}}z_{ji}^{{t}}}{1+\rho
d_{j}}\\!\right)_{\\!\\!j\\!\in\\!\mathcal{N}_{i}\\!}\\!\\!.$

Then, Algorithm 1 can be compactly described as

$\displaystyle x_{i}^{{t}+1}$
$\displaystyle=x_{i}^{{t}}+\gamma\left(h_{i}^{x}(x_{i}^{{t}},z_{i}^{{t}})-x_{i}^{{t}}\right)-\gamma\delta
h_{i}^{\nabla}(x_{i}^{{t}},z_{i}^{{t}})$ (15a) $\displaystyle z_{i}^{{t}+1}$
$\displaystyle=(1-\alpha)z_{i}^{{t}}+\alpha(-z_{\mathcal{N}_{i}}^{{t}}+2\rho
h_{\mathcal{N}_{i}}(x_{\mathcal{N}_{i}}^{t},z_{\mathcal{N}_{i}}^{{t}})).$
(15b)

Now, let us provide the aggregate formulation of (15). To this end, let us
introduce the permutation matrix $P\in\mathbb{R}^{2n\mathbf{d}\times
2n\mathbf{d}}$ that swaps the $ij$-th element with the $ji$-th one, the
matrices $A_{x}\in\mathbb{R}^{2n\mathbf{d}\times Nn}$,
$A_{z}\in\mathbb{R}^{n\mathbf{d}\times Nn}$,
$A\in\mathbb{R}^{2n\mathbf{d}\times Nn}$, $H\in\mathbb{R}^{Nn\times Nn}$, and
$\mathcal{H}\in\mathbb{R}^{2Nn\times 2Nn}$ defined as

$\displaystyle A_{x}$ $\displaystyle:=\begin{bmatrix}\mathbf{1}_{d_{1},n}\\\
0_{d_{1}n}\\\ &\mathmakebox[1.8em]{\raise 4.0pt\hbox{.}\mkern 6.0mu\raise
1.0pt\hbox{.}\mkern 6.0mu\raise-2.0pt\hbox{.}}\\\ &&\mathbf{1}_{d_{N},n}\\\
&&0_{d_{N}n}\end{bmatrix}\\!\\!,\hskip
2.27626ptA_{z}:=\begin{bmatrix}0_{d_{1}n}\\\ \mathbf{1}_{d_{1},n}\\\
&\mathmakebox[1.8em]{\raise 4.0pt\hbox{.}\mkern 6.0mu\raise
1.0pt\hbox{.}\mkern 6.0mu\raise-2.0pt\hbox{.}}\\\ &&0_{d_{N}n}\\\
&&\mathbf{1}_{d_{N},n}\end{bmatrix}$ $\displaystyle A$
$\displaystyle:=\begin{bmatrix}\mathbf{1}_{d_{1},2n}\\\
&\mathmakebox[1.8em]{\raise 4.0pt\hbox{.}\mkern 6.0mu\raise
1.0pt\hbox{.}\mkern 6.0mu\raise-2.0pt\hbox{.}}\\\
&&\mathbf{1}_{d_{N},2n}\end{bmatrix}$ $\displaystyle H$
$\displaystyle:=\begin{bmatrix}\frac{1}{1+\rho d_{1}}I_{n}\\\
&\mathmakebox[1.8em]{\raise 4.0pt\hbox{.}\mkern 6.0mu\raise
1.0pt\hbox{.}\mkern 6.0mu\raise-2.0pt\hbox{.}}\\\ &&\frac{1}{1+\rho
d_{N}}I_{n}\end{bmatrix}$ $\displaystyle\mathcal{H}$
$\displaystyle:=\begin{bmatrix}\frac{1}{1+\rho d_{1}}I_{2n}\\\
&\mathmakebox[1.8em]{\raise 4.0pt\hbox{.}\mkern 6.0mu\raise
1.0pt\hbox{.}\mkern 6.0mu\raise-2.0pt\hbox{.}}\\\ &&\frac{1}{1+\rho
d_{N}}I_{2n}\end{bmatrix}.$

Then, we define the stacking vectors
$x^{{t}}:=\textsc{col}(x^{{t}}_{1},\dots,x^{{t}}_{N})\in\mathbb{R}^{Nn}$ and
$z^{{t}}:=\textsc{col}(z^{{t}}_{1},\dots,z^{{t}}_{N})\in\mathbb{R}^{2n\mathbf{d}}$.
The aggregate formulation of (15) reads as

$\displaystyle x^{{{t}+1}}$
$\displaystyle=x^{{t}}+\gamma\left(H\left(x^{{t}}+A_{x}^{\top}z^{{t}}\right)-x^{{t}}\right)$
$\displaystyle\hskip 14.22636pt-\gamma\delta
H\left(G(x^{{t}})+A_{\nabla}^{\top}z^{{t}}\right)$ (16a) $\displaystyle
z^{{{t}+1}}$ $\displaystyle=z^{{t}}-\alpha(I+P-2\rho
PA\mathcal{H}A^{\top})z^{{t}}$ $\displaystyle\hskip 14.22636pt+2\alpha\rho
PA\mathcal{H}v(x^{{t}}),$ (16b)

where we introduced the operators $G:\mathbb{R}^{Nn}\to\mathbb{R}^{Nn}$ and
$v:\mathbb{R}^{Nn}\to\mathbb{R}^{2Nn}$ that, given any
$x:=\textsc{col}(x_{1},\dots,x_{N})\in\mathbb{R}^{Nn}$ with
$x_{i}\in\mathbb{R}^{n}$ for all $i\in\\{1,\ldots,N\\}$, are defined as
$G(x):=\textsc{col}(\nabla f_{1}(x_{1}),\dots,\nabla f_{N}(x_{N}))$ and
$v(x):=\textsc{col}(x_{1},\nabla f_{1}(x_{1}),\dots,x_{N},\nabla
f_{N}(x_{N}))$. Fig. 1 reports a block diagram graphically describing (16).

Figure 1: Block diagram representing (16).

From the discussion in [18], we deduce that the matrix $(I+P-2\rho
PA\mathcal{H}A^{\top})$ has some eigenvalues in $0$ that are all semi-simple.
Next, we introduce a decomposition to express $z^{{t}}$ according to a basis
of the subspace corresponding to the kernel of $(I+P-2\rho
PA\mathcal{H}A^{\top})$. To this end, let $b\in\mathbb{N}$ be the dimension of
the subspace $\mathcal{S}$ spanned by the eigenvectors of $(I+P-2\rho
PA\mathcal{H}A^{\top})$ associated to $0$, $B\in\mathbb{R}^{2n\mathbf{d}\times
b}$ be the matrix whose columns represent an orthonormal basis of
$\mathcal{S}$, and $M\in\mathbb{R}^{2n\mathbf{d}\times n_{d}}$ be the matrix
such that $M^{\top}M=I$ and $B^{\top}M=0$, with $n_{d}:=2n\mathbf{d}-b$. The
next lemma highlights some useful properties of the matrices $B$ and $M$.

######  Lemma IV.1.

Consider the matrices $B$ and $M$. Then, it holds

$\displaystyle A_{x}^{\top}B$ $\displaystyle=0,\quad A_{\nabla}^{\top}B=0$
(17a) $\displaystyle B^{\top}PA$ $\displaystyle=0$ (17b) $\displaystyle
B^{\top}(I\\!+\\!P\\!-\\!2\rho PA\mathcal{H}A^{\top})$ $\displaystyle=0.$
(17c)

###### Proof.

By invoking [18, Lemma 2], we claim that

$\displaystyle\ker(I+P-2\rho PA\mathcal{H}A^{\top})\subset\ker(A^{\top}).$
(18)

Further, by construction of $A_{x}$ and $A_{\nabla}$, the results
$\ker(A)\subseteq\ker(A_{x})$ and $\ker(A)\subseteq\ker(A_{\nabla})$ hold
true. Thus, the proof of (17a) follows. In order to prove (17b), we note that
(18) implies $(I+P)B=0$. Hence, since $P=P^{\top}$, it holds

$\displaystyle B^{\top}(I+P)=0.$ (19)

Further, the result (18) also implies

$\displaystyle B^{\top}A=0.$ (20)

Moreover, by construction, it holds

$\displaystyle\ker(A^{\top}P)\subseteq\ker(A^{\top}).$ (21)

The result (17b) follows by combining (20), (21) and the fact that
$P=P^{\top}$. Finally, as for (17c), it follows by combining (17b) and (19). ∎

Now, let us introduce the novel coordinates $\bar{z}\in\mathbb{R}^{b}$ and
$z_{\perp}\in\mathbb{R}^{n_{d}}$ defined as

$\displaystyle\bar{z}:=B^{\top}z,\quad z_{\perp}:=M^{\top}z.$ (22)

Then, by invoking Lemma IV.1, we use the coordinates (22) to equivalently
rewrite system (16) as

$\displaystyle x^{{{t}+1}}$
$\displaystyle=x^{{t}}+\gamma\left(H\left(x^{{t}}+A_{x}^{\top}Mz_{\perp}^{{t}}\right)-x^{{t}}\right)$
$\displaystyle\hskip 14.22636pt-\gamma\delta
H\left(G(x^{{t}})+A_{\nabla}^{\top}Mz_{\perp}^{{t}}\right)$ (23a)
$\displaystyle\bar{z}^{{t}+1}$ $\displaystyle=\bar{z}^{{t}}$ (23b)
$\displaystyle z_{\perp}^{{t}+1}$ $\displaystyle=z_{\perp}^{{t}}-\alpha
M^{\top}(I+P-2\rho PA\mathcal{H}A^{\top})Mz_{\perp}^{{t}}$
$\displaystyle\hskip 14.22636pt+2\alpha\rho M^{\top}PA\mathcal{H}v(x^{{t}}).$
(23c)

Notably, the variable $\bar{z}^{{t}}$ does not affect the other updates of
(23) (and it holds $\bar{z}^{{t}}\equiv\bar{z}^{0}$ for all ${t}\geq 0$).
Thus, by ignoring the variable $\bar{z}^{{t}}$, and introducing the error
coordinate $\tilde{x}^{{t}}:=x^{{t}}-\mathbf{1}_{N,n}x^{\star}$, system (23)
is equivalent to

$\displaystyle\tilde{x}^{{t}+1}$
$\displaystyle=\tilde{x}^{{t}}+\gamma\left(H\left(x^{{t}}+A_{x}^{\top}Mz_{\perp}^{{t}}\right)-x^{{t}}\right)$
$\displaystyle\hskip 14.22636pt-\gamma\delta
H\left(G(x^{{t}})+A_{\nabla}^{\top}Mz_{\perp}^{{t}}\right)$ (24a)
$\displaystyle z_{\perp}^{{t}+1}$ $\displaystyle=z_{\perp}^{{t}}\\!-\\!\alpha
Fz_{\perp}^{{t}}\\!+\\!2\alpha\rho M^{\top}PA\mathcal{H}v(x^{{t}}),$ (24b)

where, with a slight abuse of notation, we used
$x^{{t}}=\tilde{x}^{{t}}+\mathbf{1}_{N,n}x^{\star}$ and introduced the matrix
$F\in\mathbb{R}^{n_{d}\times n_{d}}$ defined as

$\displaystyle F:=M^{\top}(I+P-2\rho PA\mathcal{H}A^{\top})M.$ (25)

Next, we provide a crucial property about the matrix $I-\alpha F$.

######  Lemma IV.2.

Consider the matrix $F$ as defined in (25). Then, the matrix $I-\alpha F$ is
Schur for any $\alpha\in(0,1)$.

###### Proof.

By using (17c) and (25), it holds

$\displaystyle\begin{bmatrix}B^{\top}\\\
M^{\top}\end{bmatrix}(I-\alpha(I+P-2\rho
PA\mathcal{H}A^{\top}))\begin{bmatrix}B&M\end{bmatrix}$
$\displaystyle=\begin{bmatrix}I&0\\\ 0&I-\alpha F\end{bmatrix}.$

The result [29, Lemma 2] ensures that all the eigenvalues of
$I-\alpha(I+P-2\rho PA\mathcal{H}A^{\top})$ are inside the circle on the
complex plane with center $1-\alpha$ and radius $\alpha$. Since we isolated
the eigenvalues equal to $1$, we deduce that the matrix $I-\alpha F$ is Schur
for any $\alpha\in(0,1)$. ∎

We interpret (24) as a singularly perturbed system in the form of (2), i.e.,
the interconnection between the slow subsystem (24a) and the fast one (24b).
Indeed, the fast scheme (24b) has an equilibrium parametrized in the slow
state through the function
$z_{\perp}^{\text{eq}}:\mathbb{R}^{Nn}\to\mathbb{R}^{n_{d}}$ defined as

$\displaystyle z_{\perp}^{\text{eq}}(\tilde{x}^{{t}})$ $\displaystyle:=2\rho
F^{-1}M^{\top}PA\mathcal{H}v(\tilde{x}^{{t}}+\mathbf{1}_{N,n}x^{\star}).$ (26)

Moreover, given any $\tilde{x}\in\mathbb{R}^{Nn}$ and the corresponding
$x=\tilde{x}+\mathbf{1}_{N,n}x^{\star}$, it is possible to show that

$\displaystyle HA_{x}^{\top}Mz_{\perp}^{\text{eq}}(\tilde{x})$
$\displaystyle=\frac{\mathbf{1}_{N,n}\mathbf{1}_{N,n}^{\top}}{N}x-Hx$ (27a)
$\displaystyle HA_{\nabla}^{\top}Mz_{\perp}^{\text{eq}}(\tilde{x})$
$\displaystyle=\frac{\mathbf{1}_{N,n}\mathbf{1}_{N,n}^{\top}}{N}G(x)-HG(x).$
(27b)

###  IV-B Boundary Layer System

In this section, we will study the so-called boundary layer system associated
to (24). Therefore, we consider an arbitrarily fixed
$\tilde{x}^{{t}}\equiv\tilde{x}\in\mathbb{R}^{Nn}$ for all ${t}\in\mathbb{N}$
and accordingly rewrite (24b) using the error coordinate
$\tilde{z}_{\perp}^{t}:=z_{\perp}^{{t}}-z_{\perp}^{\text{eq}}(\tilde{x})$.
Hence, the definition of $z_{\perp}^{\text{eq}}$ (cf. (26)) leads to

$\displaystyle\tilde{z}_{\perp}^{{{t}+1}}=\tilde{z}_{\perp}^{t}-\alpha
F\tilde{z}_{\perp}^{t}.$ (28)

The next lemma provides a function $U$ satisfying (5) for (28).

######  Lemma IV.3.

Consider (28). Then, there exists $U:\mathbb{R}^{n_{d}}\to\mathbb{R}$ such
that the conditions (5) are satisfied.

###### Proof.

We recall that the matrix $I-\alpha F$ is Schur for any $\alpha\in(0,1)$ (cf.
Lemma IV.2). Then, we arbitrarily choose
$Q_{\tilde{z}_{\perp}}\in\mathbb{R}^{n_{d}\times n_{d}}$,
$Q_{\tilde{z}_{\perp}}=Q_{\tilde{z}_{\perp}}^{\top}>0$ and, in view of the
Schur property of $I-\alpha F$, we claim that there exists
$S_{\tilde{z}_{\perp}}\in\mathbb{R}^{n_{d}\times n_{d}}$,
$S_{\tilde{z}_{\perp}}=S_{\tilde{z}_{\perp}}^{\top}>0$ such that

$\displaystyle(I-\alpha F)^{\top}S_{\tilde{z}_{\perp}}(I-\alpha
F)-S_{\tilde{z}_{\perp}}=-Q_{\tilde{z}_{\perp}}.$ (29)

We then choose the candidate Lyapunov function
$U(\tilde{z}_{\perp})=\tilde{z}_{\perp}^{\top}S_{\tilde{z}_{\perp}}\tilde{z}_{\perp}$.
Conditions (5a) and (5c) are satisfied since the chosen $U$ is a quadratic
positive definite function, while (5b) follows by applying (29). ∎

###  IV-C Reduced System

Now, we study the so-called reduced system associated to (24), i.e., the
system obtained by considering
$z_{\perp}^{{t}}=z_{\perp}^{\text{eq}}(\tilde{x}^{{t}})$ into (24a) for all
${t}\geq 0$. By using (27), such a system reads as

$\displaystyle\tilde{x}^{{t}+1}$
$\displaystyle=\tilde{x}^{{t}}-\gamma\left(I-\frac{\mathbf{1}_{N,n}\mathbf{1}_{N,n}^{\top}}{N}\right)\tilde{x}^{{t}}$
$\displaystyle\hskip
14.22636pt-\gamma\delta\frac{\mathbf{1}_{N,n}\mathbf{1}_{N,n}^{\top}}{N}G(\tilde{x}^{{t}}+\mathbf{1}_{N,n}x^{\star}).$
(30)

The next lemma provides a Lyapunov function for system (IV-C) satisfying the
conditions required in (6).

######  Lemma IV.4.

Consider (IV-C). Then, there exists $W:\mathbb{R}^{Nn}\to\mathbb{R}$ such
that, for any $\gamma\in(0,1)$ and $\delta\in(0,\min\\{2\mathbf{c}/(\gamma
NL^{2}),2\\})$, $W$ satisfies (6).

###### Proof.

We start by introducing an additional change of variables to isolate (i) the
optimality error and (ii) the consensus error related to the vector
$\tilde{x}^{{t}}$. To this end, let $R\in\mathbb{R}^{Nn\times(N-1)n}$ be the
matrix such that its columns span the space orthogonal to the one of
$\mathbf{1}_{N,n}$, i.e., such that $R^{\top}\mathbf{1}_{N,n}=0$ and
$R^{\top}R=I_{(N-1)n}$. Then, let us introduce the new coordinates
$\mu\in\mathbb{R}^{n}$ and $x_{\perp}\in\mathbb{R}^{(N-1)n}$ be defined as

$\displaystyle\mu:=\frac{\mathbf{1}_{N,n}^{\top}}{N}\tilde{x},\quad
x_{\perp}:=R^{\top}\tilde{x}.$ (31)

By using the coordinates (31),
$\mathbf{1}_{N,n}\in\ker(I_{Nn}-\tfrac{\mathbf{1}_{N,n}\mathbf{1}_{N,n}^{\top}}{N})$,
and $R^{\top}\tfrac{\mathbf{1}_{N,n}\mathbf{1}_{N,n}^{\top}}{N}=0$, we rewrite
(IV-C) as the cascade system

$\displaystyle\mu^{{{t}+1}}$
$\displaystyle=\mu^{t}-\gamma\frac{\mathbf{1}_{N,n}^{\top}}{N}G(\mathbf{1}_{N,n}\mu^{t}+Rx_{\perp}^{{t}}+\mathbf{1}_{N,n}x^{\star})$
(32a) $\displaystyle x_{\perp}^{{t}+1}$
$\displaystyle=(1-\gamma)x_{\perp}^{{t}}.$ (32b)

For the sake of compactness, let us introduce
$\tilde{G}:\mathbb{R}^{n}\times\mathbb{R}^{(N-1)n}\to\mathbb{R}^{n}$ defined
as

$\displaystyle\tilde{G}\left(\mu,x_{\perp}\right)$
$\displaystyle:=-\frac{\mathbf{1}_{N,n}^{\top}}{N}G(\mathbf{1}_{N,n}\mu+Rx_{\perp}+\mathbf{1}_{N,n}x^{\star})$
$\displaystyle\hskip
14.22636pt+\frac{\mathbf{1}_{N,n}^{\top}}{N}G(\mathbf{1}_{N,n}\mu+\mathbf{1}_{N,n}x^{\star}).$
(33)

By using this notation and adding and subtracting $\frac{\gamma}{N}\nabla
f(\mu^{t}+x^{\star})$ into (32a), we compactly rewrite system (32) as

$\displaystyle\mu^{{{t}+1}}$
$\displaystyle=\mu^{t}-\frac{\gamma\delta}{N}\nabla
f(\mu^{t}+x^{\star})+\gamma\delta\tilde{G}\left(\mu^{t},x_{\perp}^{{t}}\right)$
(34a) $\displaystyle x_{\perp}^{{t}+1}$
$\displaystyle=(1-\gamma)x_{\perp}^{{t}}.$ (34b)

Once this formulation is available, we consider the candidate Lyapunov
function $W:\mathbb{R}^{Nn}\times\mathbb{R}^{(N-1)n}\to\mathbb{R}$ defined as

$\displaystyle W(\tilde{x})$
$\displaystyle:=\frac{1}{2}\tilde{x}^{\top}\frac{\mathbf{1}_{N,n}\mathbf{1}_{N,n}^{\top}}{N^{2}}\tilde{x}+\frac{m}{2}\tilde{x}^{\top}RR^{\top}\tilde{x}$
$\displaystyle=\underbrace{\frac{1}{2}\left\|\mu\right\|^{2}}_{:=W_{\mu}(\mu)}+\underbrace{\frac{m}{2}\left\|x_{\perp}\right\|^{2}}_{:=W_{x_{\perp}}(x_{\perp})},$
(35)

where $m>0$ will be fixed later. Being $W$ quadratic and positive definite,
the conditions (6a) and (6c) are satisfied. To check condition (6b), we write
$\Delta W_{\mu}(\mu^{t}):=W_{\mu}(\mu^{{{t}+1}})-W_{\mu}(\mu^{t})$ along
(34a), thus obtaining

$\displaystyle\Delta W_{\mu}(\mu^{t})$
$\displaystyle=-\frac{\gamma\delta}{N}(\mu^{t})^{\top}\nabla
f(\mu^{t}+x^{\star})+\gamma\delta(\mu^{t})^{\top}\tilde{G}(\mu^{t},x_{\perp}^{{t}})$
$\displaystyle\hskip 14.22636pt+\gamma^{2}\delta^{2}\nabla
f(\mu^{t}+x^{\star})^{\top}\tilde{G}(\mu^{t},x_{\perp}^{{t}})$
$\displaystyle\hskip 14.22636pt+\frac{\gamma^{2}\delta^{2}}{2}\left\|\nabla
f(\mu^{t}+x^{\star})\right\|^{2}+\frac{\gamma^{2}\delta^{2}}{2}\left\|\tilde{G}(\mu^{t},x_{\perp}^{{t}})\right\|^{2}\\!\\!.$
(36)

Since $x^{\star}$ is the (unique) solution of problem (1), it holds $\nabla
f(x^{\star})=0$ which allows us to write

$\displaystyle-(\mu^{t})^{\top}\nabla f(\mu^{t}+x^{\star})$
$\displaystyle=-(\mu^{t})^{\top}\left(\nabla f(\mu^{t}+x^{\star})-\nabla
f(x^{\star})\right)$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{\leq}}-\mathbf{c}\left\|\mu^{t}\right\|^{2},$
(37)

where in $(a)$ we use the $\mathbf{c}$-strong convexity of $f$ (cf. Assumption
II.2). Then, since $\nabla f(x^{\star})=0$ and the gradients $\nabla f_{i}$
are Lipschitz continuous (cf. Assumption II.2), we get

$\displaystyle\left\|\nabla f(\mu^{t}+x^{\star})\right\|$
$\displaystyle=\left\|\nabla f(\mu^{t}+x^{\star})-\nabla f(x^{\star})\right\|$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{\leq}}L\left\|\mu^{t}\right\|.$
(38)

Then, by using the Lipschitz continuity of the gradients $\nabla f_{i}$ (cf.
Assumption II.2) and $\left\|R\right\|=1$, we write the bound

$\displaystyle\left\|\tilde{G}(\mu^{t},x_{\perp}^{{t}})\right\|\leq\tfrac{L}{\sqrt{N}}\left\|x_{\perp}^{{t}}\right\|.$
(39)

Hence, by using (37), (38), and (39), we bound (36) as

$\displaystyle\Delta W_{\mu}(\mu^{t})$
$\displaystyle\leq-\frac{\gamma\delta}{N}\mathbf{c}\left\|\mu^{t}\right\|^{2}+\gamma\delta
k_{1}L\left\|\mu^{t}\right\|\left\|x_{\perp}^{{t}}\right\|$
$\displaystyle\hskip
14.22636pt+\gamma^{2}\delta^{2}k_{1}^{2}\left\|\mu^{t}\right\|\left\|x_{\perp}^{{t}}\right\|+\gamma^{2}\delta^{2}\frac{L^{2}}{2}\left\|\mu^{t}\right\|^{2}$
$\displaystyle\hskip
14.22636pt+\gamma^{2}\delta^{2}\frac{k_{1}^{2}}{2}\left\|x_{\perp}^{{t}}\right\|^{2},$
(40)

where the notation has been shortened by introducing $k_{1}:=L/\sqrt{N}$.
Along the trajectories (34b), we bound $\Delta
W_{x_{\perp}}(x_{\perp}^{{t}}):=W_{x_{\perp}}(x_{\perp}^{{t}+1})-W_{x_{\perp}}(x_{\perp}^{{t}})$
as

$\displaystyle\Delta W_{x_{\perp}}(x_{\perp}^{{t}})$
$\displaystyle\leq-m\gamma\left(1-\frac{\gamma}{2}\right)\left\|x_{\perp}^{{t}}\right\|^{2},$
(41)

Then, by defining $\Delta W(\tilde{x}^{{t}}):=\Delta W_{\mu}(\mu^{t})+\Delta
W_{x_{\perp}}(x_{\perp}^{{t}})$ and using (40) and (41), we get

$\displaystyle\Delta W(\tilde{x}^{{t}})$
$\displaystyle\leq-\gamma\begin{bmatrix}\left\|\mu^{t}\right\|\\\
\left\|x_{\perp}^{{t}}\right\|\end{bmatrix}^{\top}Q\begin{bmatrix}\left\|\mu^{t}\right\|\\\
\left\|x_{\perp}^{{t}}\right\|\end{bmatrix},$ (42)

where we introduced the matrix $Q=Q^{\top}\in\mathbb{R}^{2\times 2}$ given by

$\displaystyle
Q:=\begin{bmatrix}\delta\mathbf{c}/N-\gamma\delta^{2}L^{2}/2&-\delta(k_{1}L+\gamma\delta
k_{1}^{2})/2\\\ -\delta(k_{1}L+\gamma\delta
k_{1}^{2})/2&m(1-\gamma/2)-\gamma\delta^{2}k_{1}^{2}/2\end{bmatrix}$

By Sylvester Criterion, it holds $Q>0$ if and only if

$\displaystyle\begin{cases}\delta\mathbf{c}/N-\gamma\delta^{2}L^{2}/2>0\\\
m>\bar{m}(\gamma,\delta),\end{cases}$ (43)

where $\bar{m}(\gamma,\delta)$ is defined as

$\displaystyle\bar{m}(\gamma,\delta):=\frac{\gamma^{2}\delta
k_{1}^{2}/2(\mathbf{c}/N-\gamma\delta L^{2}/2)+\delta(k_{1}L+\delta\gamma
k_{1}^{2})^{2}/4}{\gamma(\mathbf{c}/N-\gamma\delta L^{2}/2)(1-\gamma/2)}.$

Hence, we fix $\gamma\in(0,1)$, then pick
$\delta\in(0,\min\\{\tfrac{2\mathbf{c}}{\gamma NL^{2}},2\\})$, and, finally,
choose $m>\bar{m}(\gamma,\delta)$. In this way, both the conditions (43) are
satisfied. Thus, we use $q>0$ to denote the smallest eigenvalue of $Q$ and
bound (42) as

$\displaystyle\Delta W(\tilde{x}^{{t}})$ $\displaystyle\leq-\gamma
q(\left\|\mu^{t}\right\|^{2}+\left\|x_{\perp}^{{t}}\right\|^{2})$
$\displaystyle\leq-\gamma
q\tilde{x}^{\top}\left(\frac{\mathbf{1}_{N,n}\mathbf{1}_{N,n}^{\top}}{N^{2}}+RR^{\top}\right)\tilde{x}^{{t}},$

which ensures that also condition (6b) is satisfied since
$\frac{\mathbf{1}_{N,n}\mathbf{1}_{N,n}^{\top}}{N^{2}}+RR^{\top}>0$. ∎

###  IV-D Proof of Theorem III.1

The proof of Theorem III.1 is based on the exploitation of Theorem II.3. In
order to apply such a result we need to (i) provide a function $U$ satisfying
the conditions (5) when applied to system (28), (ii) provide a function $W$
satisfying the conditions (6) when system (IV-C) is considered, and (iii) the
Lipschitz continuity of the dynamics of (24) and
$z_{\perp}^{\text{eq}}(\cdot)$ (cf. (26)). As for points (i) and (ii), we
invoke Lemma IV.3 and Lemma IV.4, respectively, to claim that, for any
$\alpha\in(0,1)$, $\gamma\in(0,1)$, and
$\delta\in(0,\min\\{2\mathbf{c}/(\delta NL^{2}),2\\})$, these points are
satisfied. Finally, the Lipschitz continuity of the gradients of the objective
functions (cf. Assumption II.2) allows us to claim that point (iii) is
satisfied too. Therefore, we can apply Theorem II.3 which allows us to
guarantee that there exists $\bar{\gamma}>0$ such that, for any
$\gamma\in(0,\bar{\gamma})$, the point $(0,z_{\perp}^{\text{eq}}(0))$ is
globally exponentially stable for system (24). The proof follows by turning
out to the original coordinates $(x,z)$.

##  V Asynchronous and Lossy Networks

In this section, we study the case with imperfect networks. Specifically, we
consider the case in which the agents are (possibly) asynchronous and
communicate with (possible) packet losses. More formally, for all
$i\in\\{1,\ldots,N\\}$, we introduce a stochastic process
$\lambda^{t}_{i}\in\\{0,1\\}$ modeling the fact that agent $i$ is active or
not in the following sense

$\displaystyle\lambda^{t}_{i}\\!\\!$
$\displaystyle=1\\!\\!\implies\\!\\!\\!\\!\text{ agent $i$ updates and
transmits variables at ${t}$}$ $\displaystyle\lambda^{t}_{i}\\!\\!$
$\displaystyle=0\\!\\!\implies\\!\\!\\!\\!\text{ agent $i$ does not update and
transmit variables at ${t}$}.$

Further, for each pair $(i,j)\in\mathcal{E}$, let $\beta_{ij}^{t}\in\\{0,1\\}$
be the stochastic process modeling the packet losses in the sense

$\displaystyle\beta_{ij}^{t}$ $\displaystyle=1\implies\text{ $i$ receives
message $m_{ji}^{t}$ from $j$ at ${t}$}$ $\displaystyle\beta_{ij}^{t}$
$\displaystyle=0\implies\text{ $i$ does not receive message $m_{ji}^{t}$ from
$j$ at ${t}$},$

in which $m_{ij}^{t}$ has the same meaning as in (12).

######  Remark V.1.

We note that each process $\beta_{ij}^{t}$ depends on the process
$\lambda^{t}_{j}$. Indeed, the fact that agent $j$ is active at iteration
${t}$ is a necessary (but not sufficient) condition to ensure that agent $i$
successfully receives a message from agent $j$. $\square$

In light of the above remark, we also introduce the stochastic process
$\psi_{ij}^{t}:=\lambda^{t}_{i}\beta_{ij}^{t}$ for each pair
$(i,j)\in\mathcal{E}$. Now, in order to determine the class of considered
network imperfections, we provide the following definition.

######  Definition V.2.

Let $p^{t}$ be a stochastic process. Then, $p^{t}$ is said to be mean-ergodic
if it holds

$\displaystyle\lim_{T\to\infty}\frac{1}{T}\sum_{\tau=1}^{T}p^{\tau}=\mathbb{E}[p^{t}].$
$\square$

######  Assumption V.3 (Network Imperfections).

For all $i\in\\{1,\ldots,N\\}$ and $j\in\mathcal{N}_{i}$, the stochastic
processes $\lambda^{t}_{i}$ and $\psi_{ij}^{t}$ are ergodic in mean. Further,
it holds $\mathbb{E}[\psi_{ij}^{t}]>\epsilon_{\psi}$ for all
$(i,j)\in\mathcal{E}$, for some $\epsilon_{\psi}\in(0,1)$. Finally, for all
$i\in\\{1,\ldots,N\\}$, it holds
$\mathbb{E}[\lambda^{t}_{i}]>\epsilon_{\lambda}$ for some
$\epsilon_{\lambda}\in(0,1)$. $\square$

###  V-A Robust ADMM-Tracking Gradient: Algorithm Description and Convergence
Properties

To address the challenging framework described above, we propose a slightly
different version of Algorithm 1 that we call Robust ADMM-Tracking Gradient
and report in Algorithm 2.

Initialization: $x_{i}^{0}\in\mathbb{R}^{n},z_{ij}^{0}\in\mathbb{R}^{2n}\hskip
2.84544pt\forall j\in\mathcal{N}_{i}$.

for ${t}=0,1,\dots$ do

if active then

$\begin{bmatrix}y_{i}^{{t}}\\\ s_{i}^{{t}}\end{bmatrix}=\frac{1}{1+\rho
d_{i}}\left(\begin{bmatrix}x_{i}^{{t}}\\\ \nabla
f_{i}(x_{i}^{{t}})\end{bmatrix}+\sum_{j\in\mathcal{N}_{i}}z_{ij}^{{t}}\right)$

$x_{i}^{{t}+1}=x_{i}^{{t}}+\gamma\left(y_{i}^{{t}}-x_{i}^{{t}}\right)-\gamma\delta
s_{i}^{{t}}$

for $j\in\mathcal{N}_{i}$ do

$m_{ij}^{t}=-z_{ij}^{{t}}+2\rho\begin{bmatrix}y_{i}^{{t}}\\\
s_{i}^{{t}}\end{bmatrix}$

transmit $m_{ij}^{t}$ to $j$ and receive $m_{ji}^{t}$ to $j$

if $m_{ji}^{t}$ is received then

$z_{ij}^{{t}+1}=(1-\alpha)z_{ij}^{{t}}+\alpha m_{ji}^{t}$

Algorithm 2 Robust ADMM-Tracking Gradient (Agent $i$)

In Algorithm 2, we note that each variable $z_{ij}^{{t}}$ is updated only if,
at iteration ${t}$, the message from agent $j$ has been received by agent $i$.
In other words, the update (11b) is performed only if $\beta_{ij}^{t}=1$.
Further, as one may expect from the above discussion, the additional condition
to perform such an update is that agent $i$ is active at iteration ${t}$. The
next theorem assesses the convergence properties of Algorithm 2.

######  Theorem V.4.

Consider Algorithm 2 and let Assumptions II.1, II.2, and V.3 hold. Then, there
exist $\bar{\gamma}_{\text{R}},\bar{\alpha}_{\text{R}},a_{3},a_{4}>0$ such
that, for any $\rho>0$, $\alpha\in(0,\bar{\alpha}_{\text{R}})$,
$\gamma\in(0,\bar{\gamma}_{\text{R}})$, $\delta\in(0,4\mathbf{c}/(\delta
NL^{2}))$, $(x_{i}^{0},z_{i}^{0})\in\mathbb{R}^{n}\times\mathbb{R}^{2nd_{i}}$,
for all $i\in\\{1,\ldots,N\\}$, it holds

$\displaystyle\left\|x_{i}^{{t}}-x^{\star}\right\|\leq a_{3}\exp(-a_{4}{t}).$
$\square$

The proof of Theorem V.4 is provided in Section VI-D. More in detail, the
proof takes on (i) a singular perturbation interpretation to handle the
interaction between the variables $x_{i}$ and $z_{i}$, and (ii) averaging
theory to overcome the imperfections of the networks formalized in Assumption
V.3.

######  Remark V.5.

We remark that Assumption V.3 models the scenario of time-varying networks
(with or without packet losses) as a special case. Indeed, one may interpret
$\beta_{ij}^{t}$ as the probability that edge $(i,j)$ exists at iteration
${t}$. From this perspective, Assumption V.3, combined with Assumption II.1,
guarantees that the time-varying graph has some connectivity-like property in
terms of expected value. Therefore, Theorem V.4 guarantees the linear
convergence for Robust ADMM-Tracking Gradient also in the case of time-varying
networks. $\square$

##  VI Convergence analysis of Robust ADMM-Tracking Gradient

Here, we provide the analysis of ADMM-Tracking Gradient. First, in Section
VI-A, we interpret the aggregate form of Robust ADMM-Tracking Gradient as a
_singularly perturbed_ system. Then, in Section VI-B and VI-C, respectively,
we separately study the identified subsystems by relying on Lyapunov and
averaging theory. Finally, in Section VI-D, we use the results collected in
the previous steps to prove Theorem V.4.

###  VI-A Robust ADMM-Tracking Gradient as a Singularly Perturbed System

By using the functions $h_{i}^{x}$, $h_{i}^{\nabla}$, and
$h_{\mathcal{N}_{i}}$ (cf. (14)) and the stochastic processes
$\lambda^{t}_{i}$ and $\psi_{ij}^{t}$, we compactly rewrite the local update
of Robust ADMM-Tracking Gradient as

$\displaystyle x_{i}^{{t}+1}$
$\displaystyle=(1-\lambda^{t}_{i})x_{i}^{{t}}+\lambda^{t}_{i}x_{i}^{{t}}+\lambda^{t}_{i}\gamma(\left(h_{i}^{x}(x_{i}^{{t}},\begin{bmatrix}I&0\end{bmatrix}z_{i}^{{t}})\right)-x_{i}^{{t}})$
$\displaystyle\hskip 14.22636pt-\lambda^{t}_{i}\gamma\delta
h_{i}^{\nabla}(x_{i}^{{t}},\begin{bmatrix}0&I\end{bmatrix}z_{i}^{{t}})$ (44a)
$\displaystyle z_{i}^{{t}+1}$
$\displaystyle=(I-\Psi_{i}^{t})z_{i}^{{t}}+\Psi_{i}^{t}(1-\alpha)z_{i}^{{t}}$
$\displaystyle\hskip
14.22636pt+\Psi_{i}^{t}\alpha(-z_{\mathcal{N}_{i}}^{{t}}+2\rho
h_{\mathcal{N}_{i}}(x_{\mathcal{N}_{i}}^{t},z_{\mathcal{N}_{i}}^{{t}})),$
(44b)

where we introduced the matrix $\Psi_{i}^{t}\in\mathbb{R}^{2nd_{i}\times
2nd_{i}}$ defined as

$\displaystyle\Psi_{i}^{t}=\text{blkdiag}\left(\textsc{col}(\psi_{ij}^{t})_{j\in\mathcal{N}_{i}}\otimes
I_{2n}\right).$

For the sake of readability, we algebraically rearrange (44) as

$\displaystyle x_{i}^{{t}+1}$
$\displaystyle=x_{i}^{{t}}+\gamma\lambda^{t}_{i}(\left(h_{i}^{x}(x_{i}^{{t}},\begin{bmatrix}I&0\end{bmatrix}z_{i}^{{t}})\right)-x_{i}^{{t}})$
$\displaystyle\hskip
14.22636pt-\gamma\delta\lambda^{t}_{i}h_{i}^{\nabla}(x_{i}^{{t}},\begin{bmatrix}0&I\end{bmatrix}z_{i}^{{t}})$
(45a) $\displaystyle z_{i}^{{t}+1}$
$\displaystyle=z_{i}^{{t}}+\alpha\Psi_{i}^{t}(-z_{i}^{{t}}-z_{\mathcal{N}_{i}}^{{t}}+2\rho
h_{\mathcal{N}_{i}}(x_{\mathcal{N}_{i}}^{t},z_{\mathcal{N}_{i}}^{{t}})).$
(45b)

The local update (45) leads to the aggregate form

$\displaystyle x^{{{t}+1}}$
$\displaystyle=x^{{t}}+\gamma\Lambda^{t}\left(H\left(x^{{t}}+A_{x}^{\top}z^{{t}}\right)-x^{{t}}\right)$
$\displaystyle\hskip
14.22636pt-\gamma\delta\Lambda^{t}H\left(G(x^{{t}})+A_{\nabla}^{\top}z^{{t}}\right)$
(46a) $\displaystyle z^{{{t}+1}}$
$\displaystyle=z^{{t}}\\!-\\!\alpha\Psi^{t}\left((I\\!+\\!P)z^{{t}}\\!-\\!2\rho
PA\mathcal{H}(z^{{t}}\\!+\\!v(x^{{t}}))\right)\\!,$ (46b)

where $x^{{t}}$, $z^{{t}}$, $H$, $A_{x}$, $A_{\nabla}$, $P$, and $v$ have the
same meaning as in (16) and we further introduced the matrices
$\Lambda^{t}:=\text{blkdiag}(\lambda^{t}_{1}I_{n},\dots,\lambda^{t}_{N}I_{n})\in\mathbb{R}^{Nn\times
Nn}$ and
$\Psi^{t}:=\text{blkdiag}(\psi^{{t}}_{1}I_{2nd_{1}},\dots,\psi^{{t}}_{N}I_{2nd_{N}})\in\mathbb{R}^{2n\mathbf{d}}$.
Fig. 2 reports a block diagram graphically describing (46).

Figure 2: Block diagram representing (46).

By using the change of variables (22) and Lemma IV.1, we rewrite (46) as

$\displaystyle x^{{{t}+1}}$
$\displaystyle=x^{{t}}+\gamma\Lambda^{t}\left(H\left(x^{{t}}+A_{x}^{\top}M^{\top}z_{\perp}^{{t}}\right)-x^{{t}}\right)$
$\displaystyle\hskip
14.22636pt-\gamma\delta\Lambda^{t}H\left(G(x^{{t}})+A_{\nabla}^{\top}M^{\top}z_{\perp}^{{t}}\right)$
(47a) $\displaystyle\bar{z}^{{t}+1}$ $\displaystyle=\bar{z}^{{t}}-\alpha
B^{\top}\Psi^{t}\left(I+P-2\rho PA\mathcal{H}A^{\top}\right)Mz_{\perp}^{{t}}$
$\displaystyle\hskip 14.22636pt+\alpha 2\rho
B^{\top}\Psi^{t}PA\mathcal{H}v(x^{{t}})$ (47b) $\displaystyle
z_{\perp}^{{t}+1}$ $\displaystyle=z_{\perp}^{{t}}-\alpha
M^{\top}\Psi^{t}\left(I+P-2\rho PA\mathcal{H}A^{\top}\right)Mz_{\perp}^{{t}}$
$\displaystyle\hskip 14.22636pt+\alpha 2\rho
M^{\top}\Psi^{t}PA\mathcal{H}v(x^{{t}}).$ (47c)

As in Section IV, we observe that $\bar{z}$ does not affect the other states
of system (47) and, thus, we will ignore it in our analysis. Moreover, by
using $\tilde{x}^{{t}}=x^{{t}}-\mathbf{1}_{N,n}x^{\star}$ to transform
$x^{{t}}$, system (47) is equivalent to

$\displaystyle\tilde{x}^{{t}+1}$
$\displaystyle=\tilde{x}^{{t}}+\gamma\Lambda^{t}\left(H\left(x^{{t}}+A_{x}^{\top}M^{\top}z_{\perp}^{{t}}\right)-x^{{t}}\right)$
$\displaystyle\hskip
14.22636pt-\gamma\delta\Lambda^{t}H\left(G(x^{{t}})+A_{\nabla}^{\top}M^{\top}z_{\perp}^{{t}}\right)$
(48a) $\displaystyle z_{\perp}^{{t}+1}$ $\displaystyle=z_{\perp}^{{t}}-\alpha
M^{\top}\Psi^{t}\left(I+P-2\rho PA\mathcal{H}A^{\top}\right)Mz_{\perp}^{{t}}$
$\displaystyle\hskip 14.22636pt+\alpha 2\rho
M^{\top}\Psi^{t}PA\mathcal{H}v(\tilde{x}^{{t}}+\mathbf{1}_{N,n}x^{\star}),$
(48b)

where, with a slight abuse of notation, we maintained a hybrid notation with
$x^{{t}}=\mathbf{1}_{N,n}\mu^{t}+Rx_{\perp}^{{t}}+\mathbf{1}_{N,n}x^{\star}$.
System (48) is a singularly perturbed system. Specifically, (48a) plays the
role of the slow subsystem, while (48b) plays the role of the fast one. As in
system (24), the subsystem (48b) has an equilibrium
$z_{\perp}^{\text{eq}}(\tilde{x}^{{t}})$ parametrized in the slow state
$\tilde{x}^{{t}}$. Indeed, by using $x=\tilde{x}+\mathbf{1}_{N,n}x^{\star}$,
one may check this claim through the following chain of equations

$\displaystyle(I+P-2\rho
PA\mathcal{H}A^{\top})Mz_{\perp}^{\text{eq}}(\tilde{x}^{{t}})-2\rho
PA\mathcal{H}v(x)$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}2\rho\left((I+P-2\rho
PA\mathcal{H}A^{\top})(MF^{-1}M^{\top})-I\right)PA\mathcal{H}v(x)$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{=}}2\rho\left(MM^{\top}-I\right)PA\mathcal{H}v(x)$
$\displaystyle\stackrel{{\scriptstyle(c)}}{{=}}-2\rho
BB^{\top}PA\mathcal{H}v(x)=0,$ (49)

where in $(a)$ we used the definition of $z_{\perp}^{\text{eq}}$ (cf. (26)),
in $(b)$ we used the result $(I+P-2\rho
PA\mathcal{H}A^{\top})MF^{-1}M^{\top}=MM^{\top}$, in $(c)$ we used
$MM^{\top}=I-BB^{\top}$, while the last equality follows because, since
$A^{\top}B=0$, it holds $B^{\top}PA\mathcal{H}=0$.

###  VI-B Boundary Layer System: Averaging Analysis

Here, we study the boundary layer system of (48), i.e., following the same
steps of Section IV-B, we study the subsystem (48b) by considering an
arbitrarily fixed $\tilde{x}^{{t}}=\tilde{x}\in\mathbb{R}^{Nn}$ and the error
variable
$\tilde{z}_{\perp}^{t}=z_{\perp}^{{t}}-z_{\perp}^{\text{eq}}(\tilde{x})$.
Therefore, using the result in (49), we get

$\displaystyle\tilde{z}_{\perp}^{{{t}+1}}=\tilde{z}_{\perp}^{t}-\alpha
M^{\top}\Psi^{t}\left(I+P-2\rho
PA\mathcal{H}A^{\top}\right)M\tilde{z}_{\perp}^{t}.$ (50)

It is worth noting that (50) is a time-varying system in view of the presence
of $\Psi^{t}$. For this reason, we study it by resorting to averaging theory
for discrete-time systems (cf., e.g., [41] or [42, Ch. 7]). Specifically, we
will assess the convergence properties of (50) by studying its averaged
system. To this end, we need to show that the following limit exists

$\displaystyle
g_{\text{av}}(\tilde{z}_{\perp})=\lim_{T\to\infty}\\!\frac{1}{T}\\!\sum_{\tau=\bar{{t}}+1}^{\bar{{t}}+T}M^{\top}\Psi^{\tau}\left(I\\!+\\!P\\!-\\!2\rho
PA\mathcal{H}A^{\top}\right)M\tilde{z}_{\perp},$

uniformly in $\bar{{t}}\in\mathbb{N}$ and for any
$\tilde{z}_{\perp}\in\mathbb{R}^{n_{d}}$. In this regard, since the components
of $\Psi^{t}$ are generated by a stochastic process ergodic in mean (cf.
Assumption V.3), we get

$\displaystyle
g_{\text{av}}(\tilde{z}_{\perp})=M^{\top}E_{\Psi}\left(I+P-2\rho
PA\mathcal{H}A^{\top}\right)M\tilde{z}_{\perp},$ (51)

where we introduced the symbol $E_{\Psi}:=\mathbb{E}[\Psi^{t}]$. Then, by
using (51), the averaged system associated to (50) reads as

$\displaystyle\tilde{z}_{\perp,\text{av}}^{{t}+1}=\tilde{z}_{\perp,\text{av}}^{t}-\alpha
F_{\text{av}}\tilde{z}_{\perp,\text{av}}^{t},$ (52)

where $\tilde{z}_{\perp,\text{av}}^{t}\in\mathbb{R}^{n_{d}}$ denotes the state
of the system, while the matrix $F_{\text{av}}\in\mathbb{R}^{n_{d}\times m}$
is defined as $F_{\text{av}}:=M^{\top}E_{\Psi}\left(I+P-2\rho
PA\mathcal{H}A^{\top}\right)M$. The next lemma ensures that the origin is
globally exponentially stable for (52).

######  Lemma VI.1.

Consider (52). Then, the origin is a globally exponentially stable equilibrium
point for (52).

###### Proof.

Since $I>E_{\Psi}=E_{\Psi}^{\top}>\epsilon_{\psi}I$ with
$\epsilon_{\psi}\in(0,1)$ (cf. Assumption V.3), then

$\displaystyle\ker(E_{\Psi}\left(I\\!+\\!P\\!-\\!2\rho
PA\mathcal{H}A^{\top}\right))\\!\equiv\\!\ker(I\\!+\\!P\\!-\\!2\rho
PA\mathcal{H}A^{\top}).$

Further, we recall that the columns of $M$ span the subspace orthogonal with
respect to $\ker\left(I+P-2\rho PA\mathcal{H}A^{\top}\right)$. With these
arguments and the ones in the proof of Lemma IV.2, we guarantee that the
matrix $I-\alpha F_{\text{av}}$ is Schur for any $\alpha\in(0,1)$. The proof
follows by noting that (52) is an autonomous linear system governed by a Schur
state matrix. ∎

###  VI-C Reduced System: Averaging Analysis

In this section, we will study the reduced system associated to (48), i.e.,
the system obtained by plugging
$z_{\perp}^{{t}}=z_{\perp}^{\text{eq}}(\tilde{x}^{{t}})$ for all ${t}\geq 0$
into (48a). Then, by enforcing the definition of $z_{\perp}^{\text{eq}}$ (cf.
(26)), the results in (27),
$\mathbf{1}_{N,n}\in\ker\left(\frac{\mathbf{1}_{N,n}\mathbf{1}_{N,n}^{\top}}{N}-I\right)$,
and $\mathbf{1}_{N,n}^{\top}R=0$, such a reduced system reads as

$\displaystyle\tilde{x}^{{t}+1}$
$\displaystyle=\tilde{x}^{{t}}+\gamma\Lambda^{t}\left(I-\frac{\mathbf{1}_{N,n}\mathbf{1}_{N,n}^{\top}}{N}\right)\tilde{x}^{{t}}$
$\displaystyle\hskip
14.22636pt-\gamma\delta\Lambda^{t}\frac{\mathbf{1}_{N,n}\mathbf{1}_{N,n}^{\top}}{N}G(\tilde{x}^{{t}}+\mathbf{1}_{N,n}x^{\star}).$
(53)

To simplify the notation, we equivalently write system (VI-C) as

$\displaystyle\tilde{x}^{{t}+1}$ $\displaystyle=\tilde{x}^{{t}}+\gamma
r(\tilde{x}^{{t}},{t}),$ (54)

where we introduced $r:\mathbb{R}^{Nn}\times N\to\mathbb{R}^{Nn}$ defined as

$\displaystyle r(\tilde{x}^{{t}},{t})$
$\displaystyle=\Lambda^{t}\left(\frac{\mathbf{1}_{N,n}\mathbf{1}_{N,n}^{\top}}{N}-I\right)\tilde{x}$
$\displaystyle\hskip
14.22636pt-\delta\Lambda^{t}\frac{\mathbf{1}_{N,n}\mathbf{1}_{N,n}^{\top}}{N}G(\tilde{x}+\mathbf{1}_{N,n}x^{\star}).$
(55)

We note that $r(0,0,{t})=0$ for any ${t}\in\mathbb{N}$. As in Section VI-B,
since system (VI-C) is time-varying, we study it by resorting to averaging
theory. To this end, let $r_{\text{av}}:\mathbb{R}^{Nn}\to\mathbb{R}^{Nn}$ be

$\displaystyle
r_{\text{av}}(\tilde{x}):=\lim_{T\to\infty}\frac{1}{T}\sum_{\tau=\bar{{t}}+1}^{\bar{{t}}+T}r(\tilde{x},\tau).$

Since the components of $\Lambda^{t}$ are generated by stochastic processes
ergodic in mean (cf. Assumption V.3) and using the same considerations done to
derive (51), such a limit exists uniformly in $\bar{{t}}$ and for any
$\tilde{x}\in\mathbb{R}^{Nn}$. Specifically, it reads as

$\displaystyle r_{\text{av}}(\tilde{x})$
$\displaystyle=E_{\Lambda}\left(\frac{\mathbf{1}_{N,n}\mathbf{1}_{N,n}^{\top}}{N}-I\right)\tilde{x}$
$\displaystyle\hskip 14.22636pt-\delta
E_{\Lambda}\frac{\mathbf{1}_{N,n}\mathbf{1}_{N,n}^{\top}}{N}G(\tilde{x}+\mathbf{1}_{N,n}x^{\star}),$
(56)

where we introduced the symbol
$E_{\Lambda}:=\text{blkdiag}(E_{\lambda_{1}}I,\dots,E_{\lambda_{N}}I)$, with
$E_{\lambda_{i}}:=\mathbb{E}[\lambda^{t}_{i}]$ for all $i\in\\{1,\ldots,N\\}$.
The next lemma provides a suitable function turning out to be crucial to
conclude the proof of Theorem V.4.

######  Lemma VI.2.

Consider $r_{\text{av}}(\cdot)$ as defined in (56). Then, there exist
$\bar{\delta},\beta_{1},\beta_{2},\beta_{3},\beta_{4}>0$, and
$W_{\text{av}}:\mathbb{R}^{Nn}\to\mathbb{R}$ such that, for any
$\delta\in(0,\bar{\delta})$ and $\tilde{x}_{\text{av}}\in\mathbb{R}^{Nn}$, it
holds

$\displaystyle\beta_{1}\left\|\tilde{x}_{\text{av}}\right\|^{2}\leq
W_{\text{av}}(\tilde{x}_{\text{av}})$
$\displaystyle\leq\beta_{2}\left\|\tilde{x}_{\text{av}}\right\|^{2}$ (57a)
$\displaystyle\nabla
W_{\text{av}}(\tilde{x}_{\text{av}})^{\top}r_{\text{av}}(\tilde{x}_{\text{av}})$
$\displaystyle\leq-\beta_{3}\left\|\tilde{x}_{\text{av}}\right\|^{2}$ (57b)
$\displaystyle\left\|\nabla W_{\text{av}}(\tilde{x}_{\text{av}})\right\|$
$\displaystyle\leq\beta_{4}\left\|\tilde{x}_{\text{av}}\right\|.$ (57c)

###### Proof.

Let us rewrite (56) in a block-wise sense. To this end, for all
$i\in\\{1,\ldots,N\\}$, let $r_{i,\text{av}}:\mathbb{R}^{Nn}\to\mathbb{R}^{n}$
be defined as

$\displaystyle r_{i,\text{av}}(\tilde{x}_{\text{av}})$
$\displaystyle=E_{\lambda_{i}}\left(\frac{1}{N}\sum_{j=1}^{N}\tilde{x}_{j,\text{av}}-\tilde{x}_{i,\text{av}}\right)$
$\displaystyle\hskip 14.22636pt-\delta
E_{\lambda_{i}}\mathbf{1}_{N,n}^{\top}G(\tilde{x}^{{t}}+\mathbf{1}_{N,n}x^{\star})/N,$

where we decomposed $\tilde{x}_{\text{av}}$ according to
$\tilde{x}_{\text{av}}:=\textsc{col}(\tilde{x}_{1,\text{av}},\dots,\tilde{x}_{N,\text{av}})$
with $\tilde{x}_{i,\text{av}}\in\mathbb{R}^{n}$ for all
$i\in\\{1,\ldots,N\\}$. Then, it holds
$r_{\text{av}}(\tilde{x}):=\textsc{col}(r_{1,\text{av}}(\tilde{x}_{\text{av}}),\dots,r_{N,\text{av}}(\tilde{x}_{\text{av}}))$.
Now, let us consider the function
$W_{i,\text{av}}:\mathbb{R}^{n}\to\mathbb{R}$ defined as

$\displaystyle
W_{i,\text{av}}(\tilde{x}_{i,\text{av}})=\frac{1}{2E_{\lambda_{i}}}\left\|\tilde{x}_{i,\text{av}}\right\|^{2}.$

Then, it holds

$\displaystyle\nabla
W_{i,\text{av}}(\tilde{x}_{i,\text{av}})^{\top}r_{i,\text{av}}(\tilde{x}_{\text{av}})$
$\displaystyle=-\left\|\tilde{x}_{i,\text{av}}\right\|^{2}+\tilde{x}_{i,\text{av}}^{\top}\sum_{j=1}^{N}\tilde{x}_{j,\text{av}}/N$
$\displaystyle\hskip
14.22636pt-\delta\tilde{x}_{i,\text{av}}^{\top}\mathbf{1}_{N,n}^{\top}G(\tilde{x}_{\text{av}})/N$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}-\left\|\tilde{x}_{i,\text{av}}\right\|^{2}+\tilde{x}_{i,\text{av}}^{\top}\mu_{\text{av}}$
$\displaystyle\hskip
14.22636pt-\delta\tilde{x}_{i,\text{av}}^{\top}\mathbf{1}_{N,n}^{\top}G(\tilde{x}_{\text{av}})/N,$
(58)

where in $(a)$ we used
$\mu_{\text{av}}:=\mathbf{1}_{N,n}^{\top}\tilde{x}_{\text{av}}/N$. Let
$W_{\text{av}}(\tilde{x}_{\text{av}}):=\sum_{i=1}^{N}W_{i,\text{av}}(\tilde{x}_{i,\text{av}})$.
Being $\epsilon_{\lambda}<E_{\lambda_{i}}\leq 1$ for some
$\epsilon_{\lambda}\in(0,1)$ (cf. Assumption V.3) for all
$i\in\\{1,\ldots,N\\}$, $W_{\text{av}}$ is positive definite and, thus, the
conditions (57a) and (57c) are satisfied. In order to check (57b), we use (58)
obtaining

$\displaystyle\nabla
W_{\text{av}}(\tilde{x}_{\text{av}})^{\top}\tilde{x}_{\text{av}}$
$\displaystyle=-\left\|\tilde{x}_{\text{av}}\right\|^{2}+N\mu_{\text{av}}^{\top}\mu_{\text{av}}-\mu_{\text{av}}^{\top}\mathbf{1}_{N,n}^{\top}G(\tilde{x}_{\text{av}})$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}-\left\|\mathbf{1}_{N,n}\mu_{\text{av}}+Rx_{\perp,\text{av}}\right\|^{2}+N\mu_{\text{av}}^{\top}\mu_{\text{av}}$
$\displaystyle\hskip
14.22636pt-\delta\mu_{\text{av}}^{\top}\mathbf{1}_{N,n}^{\top}G(\tilde{x}_{\text{av}})$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{=}}-\left\|x_{\perp,\text{av}}\right\|^{2}-\delta\mu_{\text{av}}^{\top}\mathbf{1}_{N,n}^{\top}G(\tilde{x}_{\text{av}})$
$\displaystyle\stackrel{{\scriptstyle(c)}}{{=}}-\left\|x_{\perp,\text{av}}\right\|^{2}-\delta\mu_{\text{av}}^{\top}\nabla
f(\mu_{\text{av}}+x^{\star})$ $\displaystyle\hskip
14.22636pt+\delta\mu_{\text{av}}^{\top}N\tilde{G}(\mu_{\text{av}}^{t},x_{\perp,\text{av}}),$
(59)

where in $(a)$ we used $x_{\perp,\text{av}}:=R^{\top}\tilde{x}_{\text{av}}$
and decomposed $\tilde{x}_{\text{av}}$ according to
$\tilde{x}_{\text{av}}=\mathbf{1}_{N,n}\mu_{\text{av}}+Rx_{\perp,\text{av}}$,
in $(b)$ we expanded the square norm and used $R^{\top}R=I$,
$\mathbf{1}_{N,n}^{\top}\mathbf{1}_{N,n}=NI$, and
$R^{\top}\mathbf{1}_{N,n}=0$, while in $(c)$ we added and subtracted the term
$\delta\mu_{\text{av}}^{\top}\nabla f(\mu_{\text{av}}+x^{\star})$ and used the
operator $\tilde{G}(\cdot,\cdot)$ with the same meaning as in (33). By using
the strong convexity of $f$ (cf. Assumption II.2) and $\nabla f(x^{\star})=0$,
we bound (59) as

$\displaystyle
W_{\text{av}}(\tilde{x}_{\text{av}})^{\top}\tilde{x}_{\text{av}}$
$\displaystyle\leq\\!-\\!\left\|x_{\perp,\text{av}}\right\|^{2}\\!-\\!\delta\mathbf{c}\left\|\mu_{\text{av}}\right\|^{2}\\!+\\!\delta\mu_{\text{av}}^{\top}N\tilde{G}(\mu_{\text{av}}^{t},x_{\perp,\text{av}})$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{\leq}}\\!-\\!\left\|x_{\perp,\text{av}}\right\|^{2}\\!-\\!\delta\mathbf{c}\left\|\mu_{\text{av}}\right\|^{2}\\!+\\!\delta
L\sqrt{N}\left\|\mu_{\text{av}}\right\|\\!\left\|x_{\perp,\text{av}}\right\|$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{\leq}}\begin{bmatrix}\left\|\mu_{\text{av}}\right\|\\\
\left\|x_{\perp,\text{av}}\right\|\end{bmatrix}^{\top}Q_{\text{av}}(\delta)\begin{bmatrix}\left\|\mu_{\text{av}}\right\|\\\
\left\|x_{\perp,\text{av}}\right\|\end{bmatrix},$ (60)

where in $(a)$ we used the Cauchy-Schwarz inequality and the result (39) to
bound the term
$\mu_{\text{av}}^{\top}N\tilde{G}(\mu_{\text{av}}^{t},x_{\perp,\text{av}})$,
while in $(b)$ we rearranged the inequality in a matrix form by introducing
$Q_{\text{av}}(\delta)=Q_{\text{av}}(\delta)^{\top}\in\mathbb{R}^{2\times 2}$
defined as

$\displaystyle Q_{\text{av}}(\delta)=\begin{bmatrix}\delta\mathbf{c}&-\delta
L\sqrt{N}/2\\\ -\delta L\sqrt{N}/2&1\end{bmatrix}.$

By Sylvester Criterion, it holds $Q_{\text{av}}(\delta)>0$ if and only if
$\delta<\frac{4\mathbf{c}}{L^{2}N}$. Therefore, by choosing any
$\delta\in(0,\frac{4\mathbf{c}}{L^{2}N})$, the inequality (60) leads to

$\displaystyle\nabla
W_{\text{av}}(\tilde{x}_{\text{av}})^{\top}r_{\text{av}}(\tilde{x}_{\text{av}})$
$\displaystyle\leq-
q_{\text{av}}\left(\left\|\mu_{\text{av}}\right\|^{2}+\left\|x_{\perp,\text{av}}\right\|^{2}\right)$
$\displaystyle=-q_{\text{av}}\tilde{x}_{\text{av}}^{\top}\left(\frac{\mathbf{1}_{N,n}\mathbf{1}_{N,n}^{\top}}{N^{2}}+RR^{\top}\right)\tilde{x}_{\text{av}},$

where $q_{\text{av}}>0$ is the smallest eigenvalue of the positive definite
matrix $Q_{\text{av}}(\delta)$. Thus, since
$\frac{\mathbf{1}_{N,n}\mathbf{1}_{N,n}^{\top}}{N^{2}}+RR^{\top}>0$, also
(57b) is ensured and the proof is concluded. ∎

###  VI-D Proof of Theorem V.4

The proof of Theorem V.4 is based on a proper exploitation of Theorem II.3.
More in details, to apply Theorem V.4, we first need to use two results
arising in averaging theory, i.e., [41, Th. 2.2.2] and [42, Prop. 7.3].
Indeed, in order to apply Theorem II.3 we need to (i) provide a function $U$
satisfying the conditions (5) when applied to system (50), (ii) provide a
function $W$ satisfying the conditions (6) when system (VI-C) is considered,
and (iii) guarantee the Lipschitz continuity of the dynamics of (48) and
$z_{\perp}^{\text{eq}}(\cdot)$ (cf. (26)).

#### Point (i)

we recall that Lemma VI.1 guarantees that the origin is a globally
exponentially stable equilibrium point of system (52), i.e., the averaged
system associated to the boundary layer system (50) of (48). Hence, by
combining this result with the Lipschitz continuity of the gradients of the
objective functions (cf. Assumption II.2) and the fact that the components of
$\Psi^{t}$ are generated by stochastic processes ergodic in mean (cf.
Assumption V.3), we can apply [41, Th. 2.2.2] to claim the existence of
$\bar{\alpha}_{\text{R}}\in(0,1)$ such that, for any
$\alpha\in(0,\bar{\alpha}_{\text{R}})$, the origin is a globally exponentially
stable equilibrium point for (50). In turn, the function $U$ satisfying (5)
can be obtained by applying the Converse Lyapunov Theorem [42, Th. 5.8].
Indeed, we notice that both system (50) and its averaged version (52) are
independent with respect to the slow state $\tilde{x}^{{t}}$. Thus, the
uniformity with respect to $\tilde{x}^{{t}}$ of the conditions (5) is
guaranteed.

#### Point (ii)

by invoking Lemma VI.2, the fact that the components of $\Lambda^{t}$ are
generated by stochastic processes ergodic in mean (cf. Assumption V.3), and
the Lipschitz continuity of the gradients of the objective functions (cf.
Assumption II.2), we can apply [42, Prop. 7.3]. Specifically, we apply [42,
Prop. 7.3] which guarantees the existence of a function $W$ and
$\bar{\gamma}_{\text{av}}>0$ such that, for any
$\gamma\in(0,\bar{\gamma}_{\text{av}})$, $W$ satisfies the conditions (6) when
the reduced system (VI-C) is considered.

#### Point (iii)

Finally, the Lipschitz continuity of the gradients of the objective functions
(cf. Assumption II.2) allows us to claim that point (iii) is satisfied too.

Once these three points have been checked, we apply Theorem II.3 and, thus,
guarantee that there exists
$\bar{\gamma}_{\text{R}}\in(0,\bar{\gamma}_{\text{av}})$ such that, for any
$\gamma\in(0,\bar{\gamma}_{\text{R}})$, the point
$(0,z_{\perp}^{\text{eq}}(0))$ is globally exponentially stable for system
(48). The proof concludes by turning out to the original coordinates $(x,z)$.

##  VII Robustness against inexact computations and communications

In this section, we deal with possible errors in the algorithm updates (due
to, e.g., quantization effects, disturbances affecting the inter-agent
communication, or stochastic gradients). In this framework, for all
$i\in\\{1,\ldots,N\\}$, we introduce $e_{i,x}^{t}\in\mathbb{R}^{n}$ and
$e_{i,z}^{t}\in\mathbb{R}^{2n\mathbf{d}_{i}}$ to model the errors affecting
the update of $x_{i}$ and $z_{i}$, respectively, at time ${t}$. Then, in this
scenario and considering asynchronous updates and packet losses in the
communication (cf. Section V), we study the local updates

$\displaystyle x_{i}^{{t}+1}$
$\displaystyle=x_{i}^{{t}}+\gamma\lambda^{t}_{i}\left(h_{i}^{x}(x_{i}^{{t}},\begin{bmatrix}I&0\end{bmatrix}z_{i}^{{t}})-x_{i}^{{t}}\right)$
$\displaystyle\hskip
14.22636pt-\gamma\delta\lambda^{t}_{i}h_{i}^{\nabla}(x_{i}^{{t}},\begin{bmatrix}0&I\end{bmatrix}z_{i}^{{t}})+e_{i,x}^{t}$
(61a) $\displaystyle z_{i}^{{t}+1}$
$\displaystyle=z_{i}^{{t}}+\alpha\Psi_{i}^{t}(-z_{i}^{{t}}-z_{\mathcal{N}_{i}}^{{t}}+2\rho
h_{\mathcal{N}_{i}}(x_{\mathcal{N}_{i}}^{t},z_{\mathcal{N}_{i}}^{{t}}))+e_{i,z}^{t}.$
(61b)

Let $e\\!:=\\!\textsc{col}(e_{x},e_{z})\\!\in\\!\mathbb{R}^{(N+2\mathbf{d})n}$
be the vector staking all the local errors, where $e_{x}\in\mathbb{R}^{Nn}$
and $e_{z}\in\mathbb{R}^{2n\mathbf{d}}$ are defined as

$\displaystyle e_{x}:=\begin{bmatrix}e_{1,x}\\\ \vdots\\\
e_{N,x}\end{bmatrix},\quad e_{z}:=\begin{bmatrix}e_{1,z}\\\ \vdots\\\
e_{N,z}\end{bmatrix}.$ (62)

The next result characterizes the robustness of (61) with respect to the
errors $e_{i,x}^{t}$ and $e_{i,z}^{t}$ in terms of ISS property.

######  Theorem VII.1.

Consider (61) and let Assumptions II.1, II.2, and V.3 hold. Then, there exist
a $\mathcal{KL}$ function $\beta(\cdot)$, a $\mathcal{K}_{\infty}$ function
$\nu(\cdot)$, and some constants
$\bar{\gamma}_{\text{R}},\bar{\alpha}_{\text{R}},\bar{\delta},c^{0}>0$ such
that, for any $\rho>0$, $\alpha\in(0,\bar{\alpha}_{R})$,
$\gamma\in(0,\bar{\gamma}_{\text{R}})$, $\delta\in(0,\bar{\delta})$,
$(x_{i}^{0},z_{i}^{0})\in\mathbb{R}^{n}\times\mathbb{R}^{2nd_{i}}$, for all
$i\in\\{1,\ldots,N\\}$, it holds

$\displaystyle\left\|x_{i}^{{t}}-x^{\star}\right\|\leq\beta(c^{0},{t})+\nu(\left\|e\right\|_{\infty}),$

for any $e\in\mathcal{L}_{\infty}^{(N+\mathbf{d})n}$ and
${t}\in\mathbb{N}$.111See [43, Chapter 4] for the function classes’
definitions.

###### Proof.

By using $\mathrm{x}$ and $\mathrm{z}$ as given in (62), the matrices and the
quantities $x^{{t}}$, $z^{{t}}$, and $G(x^{{t}})$ with the same meaning as in
(16) and (46), the aggregate form of (61) reads as

$\displaystyle x^{{{t}+1}}$
$\displaystyle=x^{{t}}+\gamma\Lambda^{t}\left(H\left(x^{{t}}+A_{x}^{\top}z^{{t}}\right)-x^{{t}}\right)$
$\displaystyle\hskip
14.22636pt-\gamma\delta\Lambda^{t}H\left(G(x^{{t}})+A_{\nabla}^{\top}z^{{t}}\right)+e_{x}^{t}$
(63a) $\displaystyle z^{{{t}+1}}$
$\displaystyle=z^{{t}}\\!-\\!\alpha\Psi^{t}\left((I\\!+\\!P)z^{{t}}\\!-\\!2\rho
PA\mathcal{H}(z^{{t}}\\!+\\!v(x^{{t}}))\right)\\!+e_{z}^{t}.$ (63b)

As in Section IV and VI, we rewrite system (63) by (i) using the changes of
coordinates (22) and $\tilde{x}^{{t}}:=x^{{t}}-\mathbf{1}_{N,n}x^{\star}$, and
(ii) ignoring the evolution of $\bar{z}^{{t}}$. Thus, we get

$\displaystyle\tilde{x}^{{t}+1}$
$\displaystyle=\tilde{x}^{{t}}+\gamma\Lambda^{t}\left(H\left(x^{{t}}+A_{x}^{\top}M^{\top}z_{\perp}^{{t}}\right)-x^{{t}}\right)$
$\displaystyle\hskip
14.22636pt-\gamma\delta\Lambda^{t}H\left(G(x^{{t}})+A_{\nabla}^{\top}M^{\top}z_{\perp}^{{t}}\right)+e_{x}^{t}$
(64a) $\displaystyle z_{\perp}^{{t}+1}$ $\displaystyle=z_{\perp}^{{t}}-\alpha
M^{\top}\Psi^{t}\left(I+P-2\rho PA\mathcal{H}A^{\top}\right)Mz_{\perp}^{{t}}$
$\displaystyle\hskip 14.22636pt+\alpha 2\rho
M^{\top}\Psi^{t}PA\mathcal{H}v(x^{{t}})+M^{\top}e_{z}^{t},$ (64b)

where, with a slight abuse of notation, we maintained a hybrid notation with
$x^{{t}}=\tilde{x}^{{t}}+\mathbf{1}_{N,n}x^{\star}$. We note that system (64)
can be viewed as a perturbed version of (48). To take advantage of this
interpretation, let $\xi\in\mathbb{R}^{(N+2\mathbf{d})n}$ be defined as

$\displaystyle\xi:=\begin{bmatrix}\tilde{x}\\\
z_{\perp}-z_{\perp}^{\text{eq}}(\tilde{x})\end{bmatrix},$ (65)

which allows us to compactly rewrite system (64) as

$\displaystyle\xi^{{{t}+1}}=f_{\xi}(\xi^{t},{t})+\tilde{e}^{t},$ (66)

where $f_{\xi}:=\mathbb{R}^{(N+2\mathbf{d})n}\to\mathbb{R}^{(N+2\mathbf{d})n}$
is suitably defined to contain all the terms arising in (64), while
$\tilde{e}^{t}\in\mathbb{R}^{(N+2\mathbf{d})n}$ is defined as
$\tilde{e}^{t}:=\textsc{col}(e_{x}^{t},M^{\top}e_{z}^{t})$. In the proof of
Theorem V.4 (cf. Section VI-D), by using Theorem II.3, we have guaranteed the
existence of $\bar{\alpha}_{\text{R}},\bar{\delta_{\text{R}}}>0$ such that,
for any $\alpha\in(0,\bar{\alpha}_{\text{R}})$ and
$\gamma\in(0,\bar{\gamma}_{\text{R}})$, the origin is a globally exponentially
stable equilibrium for the nominal system described by

$\displaystyle\xi^{{{t}+1}}=f_{\xi}(\xi^{t},{t}).$ (67)

Therefore, by using the Converse Lyapunov Theorem [42, Th. 5.8], there exist
some constants $a_{5},a_{6},a_{7},a_{8}>0$ and a Lyapunov function
$V_{\xi}:\mathbb{R}^{(N+2\mathbf{d})n}\times\mathbb{N}\to\mathbb{R}$ such that

$\displaystyle a_{5}\left\|\xi\right\|^{2}\leq V_{\xi}(\xi,{t})$
$\displaystyle\leq a_{6}\left\|\xi\right\|^{2}$ (68a) $\displaystyle
V_{\xi}(f_{\xi}(\xi,{t}),{{t}+1})-V_{\xi}(\xi,{t})$ $\displaystyle\leq-
a_{7}\left\|\xi\right\|^{2}$ (68b)
$\displaystyle|V_{\xi}(\xi_{1},{t})-V_{\xi}(\xi_{2},{t})$ $\displaystyle\leq
a_{8}\left\|\xi_{1}-\xi_{2}\right\|(\left\|\xi_{1}\right\|+\left\|\xi_{2}\right\|),$
(68c)

for any $\xi,\xi_{1},\xi_{2}\in\mathbb{R}^{(N+2\mathbf{d})n}$ and
${t}\in\mathbb{N}$. Therefore, by evaluating $\Delta
V_{\xi}(\xi^{t},{t}):=V_{\xi}(\xi^{{t}+1},{{t}+1})-V_{\xi}(\xi^{t},{t})$ along
the trajectories of system (66), we get

$\displaystyle\Delta V_{\xi}(\xi^{t},{t})$
$\displaystyle=V_{\xi}(f_{\xi}(\xi^{t},{t})+\tilde{e}^{t},{{t}+1})-V_{\xi}(\xi^{t},{t})$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}V_{\xi}(f_{\xi}(\xi^{t},{t}),{{t}+1})-V_{\xi}(\xi^{t},{t})$
$\displaystyle\hskip
14.22636pt+V(f_{\xi}(\xi^{t},{t})+\tilde{e}^{t},{{t}+1})-V_{\xi}(f_{\xi}(\xi^{t},{t}),{{t}+1})$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{\leq}}-a_{7}\left\|\xi^{t}\right\|^{2}$
$\displaystyle\hskip
14.22636pt+V(f_{\xi}(\xi^{t},{t})+\tilde{e}^{t},{{t}+1})-V_{\xi}(f_{\xi}(\xi^{t},{t}),{{t}+1})$
$\displaystyle\stackrel{{\scriptstyle(c)}}{{\leq}}-a_{7}\left\|\xi^{t}\right\|^{2}$
$\displaystyle\hskip
14.22636pt+\\!a_{8}\left\|\tilde{e}^{t}\right\|\\!(\left\|f_{\xi}(\xi^{t},{t})+\tilde{e}^{t}\right\|\\!+\\!\left\|f_{\xi}(\xi^{t},{t})\right\|),$
(69)

where in $(a)$ we add and subtract the term
$V_{\xi}(f_{\xi}(\xi^{t},{t}),{{t}+1})$, in $(b)$ we apply (68b), in $(c)$ we
apply (68c). Being the gradients of $f_{i}$ Lipschitz continuous (cf.
Assumption II.2), then $f_{\xi}(\cdot,{t})$ is Lipschitz, i.e., there exists
$L_{\xi}>0$ such that

$\displaystyle\left\|f_{\xi}(\xi_{1},{t})-f_{\xi}(\xi_{2},{t})\right\|\leq
L_{\xi}\left\|\xi_{1}-\xi_{2}\right\|,$ (70)

for any $\xi_{1},\xi_{2}\in\mathbb{R}^{Nn+n_{d}}$ and ${t}\in\mathbb{N}$.
Thus, since the origin is an equilibrium for the nominal system (67), we can
use (i) $f_{\xi}(\xi^{t},{t})=0$ for any ${t}\in N$, (ii) the result (70), and
(iii) the triangle inequality to bound the right-hand side of (69) as

$\displaystyle\Delta V_{\xi}(\xi^{t},{t})$ $\displaystyle\leq-
a_{7}\left\|\xi^{t}\right\|^{2}+a_{8}\left\|\tilde{e}^{t}\right\|^{2}+2a_{8}L_{\xi}\left\|\tilde{e}^{t}\right\|\left\|\xi^{t}\right\|$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{\leq}}-(a_{7}+\nu
a_{8}L_{\xi})\left\|\xi^{t}\right\|^{2}+\left(a_{8}+\tfrac{a_{8}L_{\xi}}{\nu}\right)\left\|\tilde{e}^{t}\right\|^{2}$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{\leq}}\\!-(a_{7}\\!+\\!\nu
a_{8}L_{\xi})\\!\left\|\xi^{t}\right\|^{2}\\!+\\!\left(\\!a_{8}\\!+\\!\tfrac{a_{8}L_{\xi}}{\nu}\\!\right)\\!\\!\left\|\tilde{e}^{t}\right\|^{2},$
(71)

where in $(a)$ we use the Young inequality with parameter $\nu>0$ to bound the
term $2a_{8}L_{\xi}\left\|\tilde{e}^{t}\right\|\left\|\xi^{t}\right\|$, while
in $(b)$ we use the Cauchy-Schwarz inequality. Hence, by arbitrarily choosing
any $a_{7}^{\prime}\in(0,a_{7})$, if we set
$\nu\in(0,(a_{7}-a_{7}^{\prime})/a_{8})$, the right-hand side of (71) can be
bounded as

$\displaystyle\Delta V_{\xi}(\xi^{t},{t})$ $\displaystyle\leq-
a_{7}^{\prime}\left\|\xi^{t}\right\|^{2}+a_{9}\left\|\tilde{e}^{t}\right\|^{2},$
(72)

with $a_{9}:=a_{8}+\frac{a_{8}L_{\xi}}{\nu}$. The result (72), combined with
(68a), allows us to claim that $V_{\xi}$ is an ISS-Lyapunov function (cf. [44,
Def. 2]) for system (66) and, thus, that system (66) is ISS [44, Th. 1]. The
proof follows by (i) using the definition of ISS systems [44, Def. 1], (ii)
using the inequalities $\left\|x_{i}^{{t}}-x^{\star}\right\|\leq\xi^{t}$ and
$\left\|\tilde{e}^{t}\right\|\leq\left\|e^{t}\right\|$, and (iii) setting
$c^{0}=\xi^{0}$. ∎

######  Remark VII.2.

The above analysis can be easily adapted to claim the ISS robustness of ADMM-
Tracking Gradient in the case of synchronous agents’ updates and without
packet losses. Indeed, the key element in the proof of Theorem VII.1 is the
exploitation of the global exponential stability of
$(0,z_{\perp}^{\text{eq}}(0))$ for (48) proved in the proof of Theorem V.4
(cf. Section VI-D). Then, since in the proof of Theorem III.1 (cf. Section
IV-D) we prove that $(0,z_{\perp}^{\text{eq}}(0))$ is globally exponentially
stable for (24), i.e., for the system that equivalently describes the
evolution of ADMM-Tracking Gradient, the same steps of the proof of Theorem
VII.1 can be used to assess the ISS of ADMM-Tracking Gradient with respect to
the generic errors $e_{i,x}^{t}$ and $e_{i,z}^{t}$ considered in (61).
$\square$

##  VIII Numerical Simulations

In this section, we perform some numerical simulations to compare ADMM-
Tracking Gradient with the well-known Gradient Tracking algorithm [4, 5, 6],
which is described by

$\displaystyle x_{i}^{{t}+1}$
$\displaystyle=\sum_{j\in\mathcal{N}_{i}}w_{ij}x_{j}^{{t}}-\gamma z_{i}^{{t}}$
(73a) $\displaystyle z_{i}^{{t}+1}$
$\displaystyle=\sum_{j\in\mathcal{N}_{i}}w_{ij}z_{j}^{{t}}+\nabla
f_{i}(x^{{{t}+1}})-\nabla f_{i}(x^{{t}}),$ (73b)

for all $i\in\\{1,\ldots,N\\}$, where $\gamma>0$ is a parameter called step-
size, $x_{i}^{{t}}\in\mathbb{R}^{n}$ is the local solution estimates
$x_{i}^{{t}}\in\mathbb{R}^{n}$, $z_{i}^{{t}}\in\mathbb{R}^{Nn}$ is the so-
called tracker approximating the unavailable global gradient, while each
weight $w_{ij}\geq 0$ is the $(i,j)$-entry of a weighted adjacency matrix
$\mathcal{W}\in\mathbb{R}^{N\times N}$ matching the graph $\mathcal{G}$, i.e.,
$w_{ij}>0$ whenever $(j,i)\in\mathcal{E}$ and $w_{ij}=0$ otherwise. Here, also
self-loops are considered, i.e., $(i,i)\in\mathcal{E}$ for all
$i\in\\{1,\ldots,N\\}$ and, thus, $w_{ii}>0$ for all $i\in\\{1,\ldots,N\\}$.
First, we consider a quadratic scenario and, then, we address a logistic
regression problem in the case in which some errors affect the update
equations of the algorithms. Finally, we test Robust ADMM-Tracking Gradient
and compare it with a robust extension of (73) by considering imperfect
networks, i.e., asynchronous agents and unreliable communications, and the
same logistic regression problem mentioned above. The comparisons are done in
terms of the convergence of the optimality error norm
$\left\|x^{{t}}-\mathbf{1}_{N,n}x^{\star}\right\|$. All the simulations are
performed by considering networks of $N=10$ agents and an underlying randomly
generated Erdős-Rényi graph with connectivity parameter $0.1$. In all the
simulations, we run our schemes by randomly selecting
$x_{i}^{0}\in\mathbb{R}^{n}$ and $z_{ij}^{0}\in\mathbb{R}^{2n}$ for all
$i\in\\{1,\ldots,N\\}$ and $j\in\mathcal{N}_{i}$. As for Gradient Tracking, we
run it by setting the same $x_{i}^{0}$ set for our schemes, while, as
prescribed in [4, 5, 6], we set $z_{i}^{0}=\nabla f_{i}(x_{i}^{0})$ for all
$i\in\\{1,\ldots,N\\}$.

###  VIII-A Quadratic Scenario

In this section, we consider a quadratic setup described by

$\min_{x\in\mathbb{R}^{n}}\>\>\sum_{i=1}^{N}\left(\tfrac{1}{2}x^{\top}Q_{i}x+r_{i}^{\top}x\right),$

where $Q_{i}\in\mathbb{R}^{n\times n}$ and $r_{i}\in\mathbb{R}^{n}$. Moreover,
it holds $Q_{i}=Q_{i}^{\top}>0$ for all $i\in\\{1,\ldots,N\\}$ and, thus, the
problem is strongly convex. We set $n=2$ and, for all $i\in\\{1,\ldots,N\\}$,
we randomly generate each matrix $Q_{i}$ so that all its eigenvalues belong to
the interval $[1,5]$, while the components of each vector $r_{i}$ are randomly
generated from the interval $[-10,20]$ with a uniform distribution. Since the
quadratic scenario gives rise to algorithm updates with linear form, we choose
the parameters of the algorithms as the ones minimizing the largest eigenvalue
of the matrices describing the algorithms in error coordinates. Specifically,
we choose $\gamma=0.865$, $\delta=0.865$, $\rho=0.3029$, $\alpha=0.865$ for
ADMM-Tracking Gradient, while we set $\gamma=0.106$ for Gradient Tracking.
Fig. 3 reports the simulation results and shows that ADMM-Tracking Gradient
outperforms Gradient Tracking in terms of convergence rate.

Figure 3:  Quadratic scenario: comparison between ADMM-Tracking Gradient (ATG)
and Gradient Tracking (GT).

###  VIII-B Logistic Regression Scenario with Errors

In this section, we consider a logistic regression scenario. A network of
agents aims to cooperatively train a linear classifier for a set of points in
a given feature space. Each agent $i$ is equipped with $m_{i}\in\mathbb{N}$
points $p_{i,1},\dots,p_{i,m_{i}}\in\mathbb{R}^{n-1}$ with binary labels
$l_{i,k}\in\\{-1,1\\}$ for all $k\in\\{1,\ldots,m_{i}\\}$. The problem
consists of building a linear classification model from these points solving
the minimization problem described by

$\min_{w\in\mathbb{R}^{n-1},b\in\mathbb{R}}\>\>\sum_{i=1}^{N}\sum_{k=1}^{m_{i}}\log\left(1+e^{-l_{i,k}(w^{\top}p_{i,k}+b)}\right)+\frac{C}{2}\left\|\begin{bmatrix}w\\\
b\end{bmatrix}\right\|^{2},$

where $C>0$ is the so-called regularization parameter. Notice that the
presence of the regularization makes the cost function strongly convex. We set
$n=2$ and we randomly generate both the points and labels. We empirically tune
the algorithms’ parameters by choosing $\alpha=0.9$, $\rho=0.9$, and
$\delta=1$ for ADMM-Tracking Gradient, while we set $\gamma=0.1$ for both
ADMM-Tracking Gradient and Gradient Tracking. Moreover, we consider the
presence of fixed errors affecting the update equations of the algorithms. In
detail, we consider an additive error $e_{i,x}^{t}=\epsilon 1_{n}$ for all
${t}\geq 0$ affecting the update equations (15a), (73a), and (73b), and,
analogously, an additive error $e_{i,z}^{t}=\epsilon 1_{2nd_{i}}$ for all
${t}\geq 0$ affecting (15b), where $\epsilon>0$ denotes the amplitudes of
these errors. Fig. 4 compares the algorithms’ performance for different values
of $\epsilon$. Differently from Gradient Tracking and as predicted by Theorem
VII.1, we note that ADMM-Tracking Gradient does not diverge and, hence, is
robust with respect to these errors.

Figure 4:  Logistic regression scenario with fixed errors: comparison between
ADMM-Tracking Gradient (ATG) and Gradient Tracking (GT).

###  VIII-C Logistic Regression Scenario wiht Asynchronous and Lossy Networks

In this section, we test the effectiveness of Robust ADMM-Tracking Gradient by
considering the same logistic regression problem tackled in Section VIII-B in
the case of imperfect networks with asynchronous agents and packet losses (see
Section V). For the sake of fairness, we compare Robust ADMM-Tracking Gradient
with a robust extension of (73) in which each agent $i$, at each iteration
${t}$, maintains local variables $x_{ij}^{t},z_{ij}^{t}\in\mathbb{R}^{n}$ for
each neighbor $j$ in order to compensate the (possible) packet losses. Then,
when agent $i$ is active, it updates its variables according to

$\displaystyle\begin{bmatrix}x_{ij}^{t}\\\ z_{ij}^{t}\end{bmatrix}$
$\displaystyle=\begin{cases}\textsc{col}(x_{j}^{t},z_{j}^{t})\hskip
28.45274pt\text{if }\psi_{ij}^{t}=1\\\
\textsc{col}(x_{ji}^{{t}-1},z_{j}^{{t}-1})\quad\text{if
}\psi_{ij}^{t}=0\end{cases}$ (74a) $\displaystyle x_{i}^{{t}+1}$
$\displaystyle=\sum_{j\in\mathcal{N}_{i}}w_{ij}x_{ij}^{t}-\gamma z_{i}^{{t}}$
(74b) $\displaystyle z_{i}^{{t}+1}$
$\displaystyle=\sum_{j\in\mathcal{N}_{i}}w_{ij}z_{ij}^{t}+\nabla
f_{i}(x^{{{t}+1}})-\nabla f_{i}(x^{{t}}).$ (74c)

For all ${t}\geq 0$, $i\in\\{1,\ldots,N\\}$, and $j\in\mathcal{N}_{i}$ we
generate each $\lambda^{t}_{i}$ and $\beta_{ij}^{t}$ by extracting them from
$\\{0,1\\}$ according to the binomial distribution with randomly generated
probabilities of success which, in turn, are uniformly randomly generated from
the interval $[0.3,1)$. We empirically tune the algorithms’ parameters by
choosing $\alpha=0.9$, $\rho=0.9$, and $\delta=1$ for ADMM-Tracking Gradient,
while we set $\gamma=0.1$ for both ADMM-Tracking Gradient and Gradient
Tracking. Fig. 5 shows that, differently from (74) and as predicted by Theorem
V.4, Robust ADMM-Tracking Gradient preserves its exact, linear convergence in
the case of imperfect networks.

Figure 5: Logistic regression scenario with imperfect networks: comparison
between Robust ADMM-Tracking Gradient (RATG) and (74) (RGT).

##  IX Conclusions

In this paper, we proposed a novel distributed algorithm for consensus
optimization and a robust extension to deal with imperfect networks. In
detail, we designed the algorithm by interpreting the dynamic consensus
problem as an additional optimization problem addressed through the ADMM. We
interlaced the obtained scheme with suitable local, proportional actions
giving rise to the novel distributed algorithm that we named ADMM-Tracking
Gradient. In the case of strongly convex problems, we proved the linear
convergence of the algorithm through a system theoretical analysis. Then, we
considered the case in which the agents’ updates are asynchronous and the
communication is affected by packet losses. For this scenario, we proposed a
robust version of the first scheme and proved that such a method preserves
linear convergence. Moreover, by relying on ISS concepts, we demonstrated that
our schemes are robust with respect to additional, generic errors. Finally,
numerical tests confirmed our theoretical findings.

## References

  * [1] G. Notarstefano, I. Notarnicola, and A. Camisa, “Distributed optimization for smart cyber-physical networks,” _Foundations and Trends® in Systems and Control_ , vol. 7, no. 3, pp. 253–383, 2019\.
  * [2] A. Nedić and J. Liu, “Distributed optimization for control,” _Annual Review of Control, Robotics, and Autonomous Systems_ , vol. 1, pp. 77–103, 2018\.
  * [3] T. Yang, X. Yi, J. Wu, Y. Yuan, D. Wu, Z. Meng, Y. Hong, H. Wang, Z. Lin, and K. H. Johansson, “A survey of distributed optimization,” _Annual Reviews in Control_ , vol. 47, pp. 278–305, 2019.
  * [4] W. Shi, Q. Ling, G. Wu, and W. Yin, “EXTRA: An exact first-order algorithm for decentralized consensus optimization,” _SIAM Journal on Optimization_ , vol. 25, no. 2, pp. 944–966, 2015.
  * [5] A. Nedić, A. Olshevsky, and W. Shi, “Achieving geometric convergence for distributed optimization over time-varying graphs,” _SIAM Journal on Optimization_ , vol. 27, no. 4, pp. 2597–2633, 2017.
  * [6] C. Xi, R. Xin, and U. A. Khan, “ADD-OPT: Accelerated distributed directed optimization,” _IEEE Transactions on Automatic Control_ , vol. 63, no. 5, pp. 1329–1339, 2017.
  * [7] N. Bof, R. Carli, G. Notarstefano, L. Schenato, and D. Varagnolo, “Multiagent Newton–Raphson Optimization Over Lossy Networks,” _IEEE Transactions on Automatic Control_ , vol. 64, no. 7, pp. 2983–2990, Jul. 2019\.
  * [8] G. Carnevale, I. Notarnicola, L. Marconi, and G. Notarstefano, “Triggered gradient tracking for asynchronous distributed optimization,” _Automatica_ , vol. 147, p. 110726, 2023.
  * [9] P. Di Lorenzo and G. Scutari, “NEXT: In-network nonconvex optimization,” _IEEE Transactions on Signal and Information Processing over Networks_ , vol. 2, no. 2, pp. 120–136, 2016.
  * [10] A. Daneshmand, G. Scutari, and V. Kungurtsev, “Second-order guarantees of distributed gradient algorithms,” _SIAM Journal on Optimization_ , vol. 30, no. 4, pp. 3029–3068, 2020.
  * [11] G. Carnevale and G. Notarstefano, “Nonconvex distributed optimization via lasalle and singular perturbations,” _IEEE Control Systems Letters_ , vol. 7, pp. 301–306, 2022.
  * [12] G. Carnevale, F. Farina, I. Notarnicola, and G. Notarstefano, “Gtadam: Gradient tracking with adaptive momentum for distributed online optimization,” _IEEE Transactions on Control of Network Systems_ , 2022.
  * [13] K. Yuan, W. Xu, and Q. Ling, “Can Primal Methods Outperform Primal-Dual Methods in Decentralized Dynamic Optimization?” _IEEE Transactions on Signal Processing_ , vol. 68, pp. 4466–4480, 2020.
  * [14] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers,” _Foundations and Trends® in Machine Learning_ , vol. 3, no. 1, pp. 1–122, 2010.
  * [15] J. F. C. Mota, J. M. F. Xavier, P. M. Q. Aguiar, and M. Püschel, “D-ADMM: A Communication-Efficient Distributed Algorithm For Separable Optimization,” _IEEE Transactions on Signal Processing_ , vol. 61, no. 10, pp. 2718–2723, May 2013.
  * [16] Z. Peng, Y. Xu, M. Yan, and W. Yin, “ARock: an Algorithmic Framework for Asynchronous Parallel Coordinate Updates,” _SIAM Journal on Scientific Computing_ , vol. 38, no. 5, pp. A2851–A2879, Jan. 2016.
  * [17] A. Makhdoumi and A. Ozdaglar, “Convergence Rate of Distributed ADMM Over Networks,” _IEEE Transactions on Automatic Control_ , vol. 62, no. 10, pp. 5082–5095, Oct. 2017.
  * [18] N. Bastianello, R. Carli, L. Schenato, and M. Todescato, “Asynchronous distributed optimization over lossy networks via relaxed admm: Stability and linear convergence,” _IEEE Transactions on Automatic Control_ , vol. 66, no. 6, pp. 2620–2635, 2020.
  * [19] Z. Peng, T. Wu, Y. Xu, M. Yan, and W. Yin, “Coordinate Friendly Structures, Algorithms and Applications,” _Annals of Mathematical Sciences and Applications_ , vol. 1, no. 1, pp. 57–119, 2016.
  * [20] H. Ye, L. Liang, and G. Y. Li, “Decentralized Federated Learning With Unreliable Communications,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 16, no. 3, pp. 487–500, Apr. 2022.
  * [21] H. Li, Z. Lin, and Y. Fang, “Variance Reduced EXTRA and DIGing and Their Optimal Acceleration for Strongly Convex Decentralized Optimization,” _Journal of Machine Learning Research_ , vol. 23, no. 222, pp. 1–41, 2022.
  * [22] J. Lei, P. Yi, J. Chen, and Y. Hong, “Distributed Variable Sample-Size Stochastic Optimization With Fixed Step-Sizes,” _IEEE Transactions on Automatic Control_ , vol. 67, no. 10, pp. 5630–5637, Oct. 2022\.
  * [23] J. Hu, G. Chen, H. Li, Z. Shen, and W. Zhang, “Push-LSVRG-UP: Distributed Stochastic Optimization Over Unbalanced Directed Networks With Uncoordinated Triggered Probabilities,” _IEEE Transactions on Network Science and Engineering_ , vol. 10, no. 2, pp. 934–950, Mar. 2023.
  * [24] P. Bianchi, W. Hachem, and F. Iutzeler, “A Coordinate Descent Primal-Dual Algorithm and Application to Distributed Asynchronous Optimization,” _IEEE Transactions on Automatic Control_ , vol. 61, no. 10, pp. 2947–2957, Oct. 2016.
  * [25] T.-H. Chang, M. Hong, W.-C. Liao, and X. Wang, “Asynchronous Distributed ADMM for Large-Scale Optimization- Part I: Algorithm and Convergence Analysis,” _IEEE Transactions on Signal Processing_ , vol. 64, no. 12, pp. 3118–3130, Jun. 2016.
  * [26] L. Majzoobi, F. Lahouti, and V. Shah-Mansouri, “Analysis of Distributed ADMM Algorithm for Consensus Optimization in Presence of Node Error,” _IEEE Transactions on Signal Processing_ , vol. 67, no. 7, pp. 1774–1784, Apr. 2019.
  * [27] T.-H. Chang, “A Proximal Dual Consensus ADMM Method for Multi-Agent Constrained Optimization,” _IEEE Transactions on Signal Processing_ , vol. 64, no. 14, pp. 3719–3734, Jul. 2016.
  * [28] A. Elgabli, J. Park, A. S. Bedi, C. B. Issaid, M. Bennis, and V. Aggarwal, “Q-GADMM: Quantized Group ADMM for Communication Efficient Decentralized Machine Learning,” _IEEE Transactions on Communications_ , vol. 69, no. 1, pp. 164–181, Jan. 2021.
  * [29] N. Bastianello and R. Carli, “Admm for dynamic average consensus over imperfect networks,” _IFAC-PapersOnLine_ , vol. 55, no. 13, pp. 228–233, 2022.
  * [30] Y. Tian, Y. Sun, and G. Scutari, “Achieving Linear Convergence in Distributed Asynchronous Multiagent Optimization,” _IEEE Transactions on Automatic Control_ , vol. 65, no. 12, pp. 5264–5279, Dec. 2020\.
  * [31] S. Pu, “A Robust Gradient Tracking Method for Distributed Optimization over Directed Networks,” in _2020 59th IEEE Conference on Decision and Control (CDC)_. Jeju Island, Korea (South): IEEE, Dec. 2020, pp. 2335–2341.
  * [32] M. Bin, I. Notarnicola, L. Marconi, and G. Notarstefano, “A system theoretical perspective to gradient-tracking algorithms for distributed quadratic optimization,” in _2019 IEEE 58th Conference on Decision and Control (CDC)_. IEEE, 2019, pp. 2994–2999.
  * [33] I. Notarnicola, M. Bin, L. Marconi, and G. Notarstefano, “The gradient tracking is a distributed integral action,” _IEEE Transactions on Automatic Control_ , 2023.
  * [34] M. Bin, I. Notarnicola, and T. Parisini, “Stability, linear convergence, and robustness of the wang-elia algorithm for distributed consensus optimization,” in _2022 IEEE 61st Conference on Decision and Control (CDC)_. IEEE, 2022, pp. 1610–1615.
  * [35] Y. Wang and T. Başar, “Gradient-tracking based Distributed Optimization with Guaranteed Optimality under Noisy Information Sharing,” _IEEE Transactions on Automatic Control_ , pp. 1–16, 2022.
  * [36] M. Zhu and S. Martínez, “Discrete-time dynamic average consensus,” _Automatica_ , vol. 46, no. 2, pp. 322–329, 2010.
  * [37] S. S. Kia, B. Van Scoy, J. Cortes, R. A. Freeman, K. M. Lynch, and S. Martinez, “Tutorial on dynamic average consensus: The problem, its applications, and the algorithms,” _IEEE Control Systems Magazine_ , vol. 39, no. 3, pp. 40–72, 2019.
  * [38] N. Bof, R. Carli, and L. Schenato, “Is ADMM always faster than Average Consensus?” _Automatica_ , vol. 91, pp. 311–315, May 2018.
  * [39] G. Carnevale, N. Bastianello, N. Bastianello, and G. Notarstefano, “Distributed consensus optimization via admm-tracking gradient,” in _2023 IEEE 62nd Conference on Decision and Control (CDC)_. IEEE, 2023, _to appear_.
  * [40] G. Carnevale, F. Fabiani, F. Fele, K. Margellos, and G. Notarstefano, “Tracking-based distributed equilibrium seeking for aggregative games,” _arXiv preprint arXiv:2210.14547_ , 2022.
  * [41] E.-W. Bai, L.-C. Fu, and S. S. Sastry, “Averaging analysis for discrete time and sampled data adaptive systems,” _IEEE Transactions on Circuits and Systems_ , vol. 35, no. 2, pp. 137–148, 1988.
  * [42] N. Bof, R. Carli, and L. Schenato, “Lyapunov theory for discrete time systems,” _arXiv preprint arXiv:1809.05289_ , 2018.
  * [43] H. K. Khalil, “Nonlinear systems,” _Upper Saddle River_ , 2002.
  * [44] H. Li, A. Liu, and L. Zhang, “Input-to-state stability of time-varying nonlinear discrete-time systems via indefinite difference lyapunov functions,” _ISA transactions_ , vol. 77, pp. 71–76, 2018.
  * [45] V. Chellaboina and W. M. Haddad, _Nonlinear dynamical systems and control: A Lyapunov-based approach_. Princeton University Press, 2008.

###  -A Distributed R-ADMM: Algorithm Description

This section briefly recalls the so-called Distributed R-ADMM proposed in
[18]. Let us consider a network of $N$ agents communicating according to the
graph $\mathcal{G}:=\\{\\{1,\ldots,N\\},\mathcal{E}\\}$, with
$\mathcal{E}\subset\\{1,\ldots,N\\}\times\\{1,\ldots,N\\}$, that aim to solve
in a distributed fashion the optimization problem

$\displaystyle\begin{split}\min_{(x_{1},\dots,x_{N})\in\mathbb{R}^{Nn}}&\sum_{i=1}^{N}\ell_{i}(x_{i})\\\
\text{s.t.:}\>\>&x_{i}=x_{j}\>\forall(i,j)\in\mathcal{E},\end{split}$ (75)

where $\ell_{i}:\mathbb{R}^{n}\to\mathbb{R}$ for all $i\in\\{1,\ldots,N\\}$.
When running Distributed R-ADMM, each agent $i$, at each iteration
${t}\in\mathbb{N}$, maintains a local estimate $x_{i}^{{t}}\in\mathbb{R}^{n}$
about a solution of problem (75) and an auxiliary variables $z_{ij}^{{t}}$ for
each $j\in\mathcal{N}_{i}$. These quantities are updated according to

$\displaystyle x_{i}^{{t}+1}$
$\displaystyle=\operatorname*{arg\,min}_{x_{i}\in\mathbb{R}^{n}}\left\\{\ell_{i}(x_{i})\\!-\\!x_{i}^{\top}\sum_{j\in\mathcal{N}_{i}}z_{ij}^{{t}}\\!+\\!\frac{\rho
d_{i}}{2}\left\|x_{i}\right\|^{2}\right\\}$ (76a) $\displaystyle
z_{ij}^{{t}+1}$ $\displaystyle=(1-\alpha)z_{ij}^{{t}}-\alpha
z_{ji}^{{t}}+2\alpha\rho x_{j}^{{t}+1},$ (76b)

where $\alpha\in(0,1)$ and $\rho>0$ are tuning parameters.

###  -B Proof of Theorem II.3

Let us define $\tilde{z}^{{t}}\coloneqq z^{{t}}-z^{\text{eq}}(x^{{t}})$ and
rewrite (2) as

$\displaystyle x^{{{t}+1}}$ $\displaystyle=x^{{t}}+\gamma
f(x^{{t}},\tilde{z}^{{t}}+z^{\text{eq}}(x^{{t}}),{t})$ (77a)
$\displaystyle\tilde{z}^{{t}+1}$
$\displaystyle=g(\tilde{z}^{{t}}+z^{\text{eq}}(x^{{t}}),x^{{t}},{t})-z^{\text{eq}}(x^{{t}})+\Delta
z^{\text{eq}}(x^{{{t}+1}},x^{{t}}),$ (77b)

where $\Delta
z^{\text{eq}}(x^{{{t}+1}},x^{{t}})\coloneqq-z^{\text{eq}}(x^{{{t}+1}})+z^{\text{eq}}(x^{{t}})$.
Pick $W$ as in (6). By evaluating $\Delta W(x^{{t}},{t})\coloneqq
W(x^{{{t}+1}},{{t}+1})-W(x^{{t}},{t})$ along the trajectories of (77a), we
obtain

$\displaystyle\Delta W(x^{{t}},{t})$ $\displaystyle=W(x^{{t}}+\gamma
f(x^{{t}},\tilde{z}^{{t}}+z^{\text{eq}}(x^{{t}}),{t}),{{t}+1})-W(x^{{t}},{t})$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}W(x^{{t}}+\gamma
f(x^{{t}},z^{\text{eq}}(x^{{t}}),{t}),{{t}+1})-W(x^{{t}},{t})$
$\displaystyle\hskip 14.22636pt+W(x^{{t}}+\gamma
f(x^{{t}},\tilde{z}^{{t}}+z^{\text{eq}}(x^{{t}}),{t}),{{t}+1})$
$\displaystyle\hskip 14.22636pt-W(x^{{t}}+\gamma
f(x^{{t}},z^{\text{eq}}(x^{{t}}),{t}),{{t}+1})$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{\leq}}-\gamma
c_{3}\left\|x^{{t}}\right\|^{2}+W(x^{{t}}+\gamma
f(x^{{t}},\tilde{z}^{{t}}+z^{\text{eq}}(x^{{t}}),{t}),{{t}+1})$
$\displaystyle\hskip 14.22636pt-W(x^{{t}}+\gamma
f(x^{{t}},z^{\text{eq}}(x^{{t}}),{t}),{{t}+1})$
$\displaystyle\stackrel{{\scriptstyle(c)}}{{\leq}}-\gamma
c_{3}\left\|x^{{t}}\right\|^{2}+2\gamma
c_{4}L_{f}\left\|\tilde{z}^{{t}}\right\|\left\|x^{{t}}\right\|$
$\displaystyle\hskip
14.22636pt+\gamma^{2}c_{4}L_{f}\left\|\tilde{z}^{{t}}\right\|\left\|f(x^{{t}},\tilde{z}^{{t}}+z^{\text{eq}}(x^{{t}}),{t})\right\|$
$\displaystyle\hskip
14.22636pt+\gamma^{2}c_{4}L_{f}\left\|\tilde{z}^{{t}}\right\|\left\|f(x^{{t}},z^{\text{eq}}(x^{{t}}),{t})\right\|,$
(78)

where in $(a)$ we add and subtract the term $W(x^{{t}}+\gamma
f(x^{{t}},z^{\text{eq}}(x^{{t}}),{t}),{{t}+1})$, in $(b)$ we exploit (6b) to
bound the difference of the first two terms, in $(c)$ we use (6c), the
Lipschitz continuity of $f$, and the triangle inequality. By recalling that
$f(0,z^{\text{eq}}(0),{t})=0$ we can thus write

$\displaystyle\left\|f(x^{{t}},\tilde{z}^{{t}}+z^{\text{eq}}(x^{{t}}),{t})\right\|$
$\displaystyle=\left\|f(x^{{t}},\tilde{z}^{{t}}+z^{\text{eq}}(x^{{t}}),{t})-f(0,z^{\text{eq}}(0),{t})\right\|$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{\leq}}L_{f}\left\|x^{{t}}\right\|+L_{f}\left\|\tilde{z}^{{t}}+z^{\text{eq}}(x^{{t}})-z^{\text{eq}}(0)\right\|,$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{\leq}}L_{f}(1+L_{\text{eq}})\left\|x^{{t}}\right\|+L_{f}\left\|\tilde{z}^{{t}}\right\|,$
(79)

where in $(a)$ we use the Lipschitz continuity of $f$ and $z^{\text{eq}}$, and
in $(b)$ we use the Lipschitz continuity of $z^{\text{eq}}$ together with the
triangle inequality. With similar arguments, we have

$\displaystyle\left\|f(x^{{t}},z^{\text{eq}}(x^{{t}}),{t})\right\|\leq
L_{f}(1+L_{\text{eq}})\left\|x^{{t}}\right\|.$ (80)

Using inequalities (79) and (80) we then bound (78) as

$\displaystyle\Delta W(x^{{t}},{t})$ $\displaystyle\leq-\gamma
c_{3}\left\|x^{{t}}\right\|^{2}+2\gamma
c_{4}L_{f}\left\|\tilde{z}^{{t}}\right\|\left\|x^{{t}}\right\|+\gamma^{2}c_{4}L_{f}^{2}\left\|\tilde{z}^{{t}}\right\|^{2}$
$\displaystyle\hskip
14.22636pt+2\gamma^{2}c_{4}L_{f}^{2}(1+L_{\text{eq}})\left\|\tilde{z}^{{t}}\right\|\left\|x^{{t}}\right\|$
$\displaystyle\leq-
c_{3}\left\|x^{{t}}\right\|^{2}\\!+\\!\gamma^{2}k_{3}\left\|\tilde{z}^{{t}}\right\|^{2}\\!+\\!(\gamma
k_{1}+\gamma^{2}k_{2})\left\|\tilde{z}^{{t}}\right\|\left\|x^{{t}}\right\|,$
(81)

where we introduce the constants

$\displaystyle k_{1}$ $\displaystyle\coloneqq 2c_{4}L_{f},\quad k_{2}\coloneqq
2c_{4}L_{f}^{2}(1+L_{\text{eq}}),\quad k_{3}\coloneqq c_{4}L_{f}^{2}.$

We now evaluate $\Delta U(\tilde{z}^{{t}},{t})\coloneqq
U(\tilde{z}^{{t}+1},{{t}+1})-U(\tilde{z}^{{t}},{t})$ (see (5)) along the
trajectories of (77b), obtaining

$\displaystyle\Delta U(\tilde{z}^{{t}},{t})$
$\displaystyle=U(g(\tilde{z}^{{t}}+z^{\text{eq}}(x^{{t}}),x^{{t}},{t})-z^{\text{eq}}(x^{{t}})+\Delta
z^{\text{eq}}(x^{{{t}+1}},x^{{t}}),{{t}+1})$ $\displaystyle\hskip
14.22636pt-U(\tilde{z}^{{t}},{t})$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{\leq}}U(g(\tilde{z}^{{t}}+z^{\text{eq}}(x^{{t}}),x^{{t}}),{t}-z^{\text{eq}}(x^{{t}}),{{t}+1})-U(\tilde{z}^{{t}},{t})$
$\displaystyle\hskip
14.22636pt-U(g(\tilde{z}^{{t}}+z^{\text{eq}}(x^{{t}}),x^{{t}},{t})-z^{\text{eq}}(x^{{t}}),{{t}+1})$
$\displaystyle\hskip
14.22636pt+U(g(\tilde{z}^{{t}}+z^{\text{eq}}(x^{{t}}),x^{{t}},{t})-z^{\text{eq}}(x^{{t}})+\Delta
z^{\text{eq}}(x^{{{t}+1}},x^{{t}}),{{t}+1})$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{\leq}}-b_{3}\left\|\tilde{z}^{{t}}\right\|^{2}-U(g(\tilde{z}^{{t}}+z^{\text{eq}}(x^{{t}}),x^{{t}},{t})-z^{\text{eq}}(x^{{t}}),{{t}+1})$
$\displaystyle\hskip
14.22636pt+U(g(\tilde{z}^{{t}}+z^{\text{eq}}(x^{{t}}),x^{{t}},{t})-z^{\text{eq}}(x^{{t}})+\Delta
z^{\text{eq}}(x^{{{t}+1}},x^{{t}}),{{t}+1})$
$\displaystyle\stackrel{{\scriptstyle(c)}}{{\leq}}-b_{3}\left\|\tilde{z}^{{t}}\right\|^{2}$
$\displaystyle\hskip 14.22636pt+b_{4}\left\|\Delta
z^{\text{eq}}(x^{{{t}+1}},x^{{t}})\right\|$ $\displaystyle\hskip
22.76228pt\times\left\|g(\tilde{z}^{{t}}+z^{\text{eq}}(x^{{t}}),x^{{t}},{t})-z^{\text{eq}}(x^{{t}})+\Delta
z^{\text{eq}}(x^{{{t}+1}},x^{{t}})\right\|$ $\displaystyle\hskip
14.22636pt+b_{4}\left\|\Delta
z^{\text{eq}}(x^{{{t}+1}},x^{{t}})\right\|\left\|g(\tilde{z}^{{t}}+z^{\text{eq}}(x^{{t}}),x^{{t}},{t})-z^{\text{eq}}(x^{{t}})\right\|$
$\displaystyle\stackrel{{\scriptstyle(d)}}{{\leq}}-b_{3}\left\|\tilde{z}^{{t}}\right\|^{2}+b_{4}\left\|\Delta
z^{\text{eq}}(x^{{{t}+1}},x^{{t}})\right\|^{2}$ $\displaystyle\hskip
14.22636pt+2b_{4}\left\|\Delta
z^{\text{eq}}(x^{{{t}+1}},x^{{t}})\right\|\left\|g(\tilde{z}^{{t}}+z^{\text{eq}}(x^{{t}}),x^{{t}},{t})-z^{\text{eq}}(x^{{t}})\right\|,$
(82)

where in $(a)$ we add and subtract
$U(g(\tilde{z}^{{t}}+z^{\text{eq}}(x^{{t}}),x^{{t}},{t})-z^{\text{eq}}(x^{{t}}),{{t}+1})$,
in $(b)$ we exploit (5b) to bound the first two terms, in $(c)$ we use (5c) to
bound the difference of the last two terms, and in $(d)$ we use the triangle
inequality. The definition of $\Delta z^{\text{eq}}(x^{{{t}+1}},x^{{t}})$ and
the Lipschitz continuity of $z^{\text{eq}}$ lead to

$\displaystyle\left\|\Delta z^{\text{eq}}(x^{{{t}+1}},x^{{t}})\right\|$
$\displaystyle\leq L_{\text{eq}}\left\|x^{{{t}+1}}-x^{{t}}\right\|$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{\leq}}\gamma
L_{\text{eq}}\left\|f(x^{{t}},\tilde{z}^{{t}}+z^{\text{eq}}(x^{{t}}),{t})\right\|$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{\leq}}\gamma
L_{\text{eq}}\left\|f(x^{{t}},\tilde{z}^{{t}}+z^{\text{eq}}(x^{{t}}),{t})-f(0,z^{\text{eq}}(0),{t})\right\|$
$\displaystyle\stackrel{{\scriptstyle(c)}}{{\leq}}\gamma
L_{\text{eq}}L_{f}(1+L_{\text{eq}})\left\|x^{{t}}\right\|+\gamma
L_{\text{eq}}L_{f}\left\|\tilde{z}^{{t}}\right\|,$ (83)

where in $(a)$ we use the update (77a), in $(b)$ we add the term
$f(0,z^{\text{eq}}(0),{t})$ since this is zero, and in $(c)$ we use the
triangle inequality and the Lipschitz continuity of $f$ and $z^{\text{eq}}$.
Moreover, since
$g(z^{\text{eq}}(x^{{t}}),x^{{t}},{t})=z^{\text{eq}}(x^{{t}})$, we obtain

$\displaystyle\left\|g(\tilde{z}^{{t}}+z^{\text{eq}}(x^{{t}}),x^{{t}},{t})-z^{\text{eq}}(x^{{t}})\right\|$
$\displaystyle=\left\|g(\tilde{z}^{{t}}+z^{\text{eq}}(x^{{t}}),x^{{t}},{t})-g(z^{\text{eq}}(x^{{t}}),x^{{t}},{t})\right\|$
$\displaystyle\leq L_{g}\left\|\tilde{z}^{{t}}\right\|,$ (84)

where the inequality is due to the Lipschitz continuity of $g$. Using
inequalities (83) and (84), we then bound (82) as

$\displaystyle\Delta U(\tilde{z}^{{t}},{t})$ $\displaystyle\leq-
b_{3}\left\|\tilde{z}^{{t}}\right\|^{2}+2\gamma
b_{4}L_{\text{eq}}L_{g}L_{f}(1+L_{\text{eq}})\left\|x^{{t}}\right\|\left\|\tilde{z}^{{t}}\right\|$
$\displaystyle\hskip 14.22636pt+2\gamma
b_{4}L_{\text{eq}}L_{g}L_{f}\left\|\tilde{z}^{{t}}\right\|^{2}$
$\displaystyle\hskip
14.22636pt+\gamma^{2}b_{4}L_{\text{eq}}^{2}L_{f}^{2}(1+L_{\text{eq}})^{2}\left\|x^{{t}}\right\|^{2}$
$\displaystyle\hskip
14.22636pt+2\gamma^{2}b_{4}L_{\text{eq}}^{2}L_{f}^{2}(1+L_{\text{eq}})\left\|x^{{t}}\right\|\left\|\tilde{z}^{{t}}\right\|$
$\displaystyle\hskip
14.22636pt+\gamma^{2}b_{4}L_{\text{eq}}^{2}L_{f}^{2}\left\|\tilde{z}^{{t}}\right\|^{2}$
$\displaystyle\leq(-b_{3}+\gamma
k_{6}+\gamma^{2}k_{7})\left\|\tilde{z}^{{t}}\right\|^{2}+\gamma^{2}k_{8}\left\|x^{{t}}\right\|^{2}$
$\displaystyle\hskip 14.22636pt+(\gamma
k_{4}+\gamma^{2}k_{5})\left\|x^{{t}}\right\|\left\|\tilde{z}^{{t}}\right\|,$
(85)

where we introduce the constants

$\displaystyle k_{4}$ $\displaystyle\coloneqq
2b_{4}L_{\text{eq}}L_{g}L_{f}(1+L_{\text{eq}}),\quad$ $\displaystyle
k_{5}\coloneqq 2b_{4}L_{\text{eq}}^{2}L_{f}^{2}(1+L_{\text{eq}}),$
$\displaystyle k_{6}$ $\displaystyle\coloneqq
2b_{4}L_{\text{eq}}L_{g}L_{f},\quad$ $\displaystyle k_{7}\coloneqq
b_{4}L_{\text{eq}}^{2}L_{f}^{2},$ $\displaystyle k_{8}$
$\displaystyle\coloneqq b_{4}L_{\text{eq}}^{2}L_{f}^{2}(1+L_{\text{eq}})^{2}.$

We pick the following Lyapunov candidate
$V:\mathcal{D}\times\mathbb{R}^{m}\to\mathbb{R}$:

$\displaystyle
V(x^{{t}},\tilde{z}^{{t}},{t})=W(x^{{t}},{t})+U(\tilde{z}^{{t}},{t}).$

By evaluating $\Delta V(x^{{t}},\tilde{z}^{{t}},{t})\coloneqq
V(x^{{{t}+1}},\tilde{z}^{{t}+1},{{t}+1})-V(x^{{t}},\tilde{z}^{{t}},{t})=\Delta
W(x^{{t}},{t})+\Delta U(\tilde{z}^{{t}},{t})$ along the trajectories of (77),
we can use the results (81) and (85) to write

$\displaystyle\Delta V(x^{{t}},\tilde{z}^{{t}},{t})$
$\displaystyle\leq-\begin{bmatrix}\left\|x^{{t}}\right\|\\\
\left\|\tilde{z}^{{t}}\right\|\end{bmatrix}^{\top}Q(\gamma)\begin{bmatrix}\left\|x^{{t}}\right\|\\\
\left\|\tilde{z}^{{t}}\right\|\end{bmatrix},$ (86)

where we define the matrix
$\mathcal{Q}(\gamma)=\mathcal{Q}(\gamma)^{\top}\in\mathbb{R}^{2}$ as

$\displaystyle\mathcal{Q}(\gamma)\coloneqq\begin{bmatrix}\gamma
c_{3}-\gamma^{2}k_{8}&q_{21}(\gamma)\\\ q_{21}(\gamma)&b_{3}-\gamma
k_{6}-\gamma^{2}(k_{3}+k_{7})\end{bmatrix},$

with
$q_{21}(\gamma)\coloneqq-\frac{1}{2}(\gamma(k_{1}+k_{4})+\gamma^{2}(k_{2}+k_{5}))$.
By Sylvester criterion [43], we know that $\mathcal{Q}>0$ if and only if

$\displaystyle\gamma c_{3}b_{3}>p(\gamma)$ (87)

where the polynomial $p(\gamma)$ is defined as

$\displaystyle p(\gamma)$ $\displaystyle\coloneqq
q_{21}(\gamma)^{2}+\gamma^{2}c_{3}k_{6}+\gamma^{3}c_{3}(k_{3}+k_{7})+\gamma^{2}b_{3}k_{8}$
$\displaystyle\hskip
14.22636pt-\gamma^{3}k_{6}k_{8}-\gamma^{4}k_{8}(k_{3}+k_{7}).$ (88)

We note that $p$ is a continuous function of $\gamma$ and $\lim_{\gamma\to
0}p(\gamma)/\gamma=0$. Hence, there exists some
$\bar{\gamma}\in(0,\min\\{\bar{\gamma}_{1},\bar{\gamma}_{2}\\})$ – recall that
$\bar{\gamma}_{1}$ and $\bar{\gamma}_{2}$ exist as $U$ and $W$ are taken to
satisfy (5) and (6) – so that (87) is satisfied for any
$\gamma\in(0,\bar{\gamma})$. Under such a choice of $\gamma$, and denoting by
$q_{\text{small}}>0$ the smallest eigenvalue of $\mathcal{Q}(\gamma)$, we
bound (86) as

$\displaystyle\Delta V(x^{{t}},\tilde{z}^{{t}},{t})\leq-
q_{\text{small}}\left\|\begin{bmatrix}\left\|x^{{t}}\right\|\\\
\left\|\tilde{z}^{{t}}\right\|\end{bmatrix}\right\|^{2},$

which allows us to conclude, in view of [45, Theorem 13.2], that $(0,0)$ is an
exponentially stable equilibrium point for system (77). The theorem’s
conclusion follows then by considering the definition of exponentially stable
equilibrium point and by reverting to the original coordinates
$(x^{{t}},z^{{t}})$.