Asynchronous Decentralized Optimization with Constraints: Achievable Speeds of
Convergence for Directed Graphs

Firooz Shahriari-Mehr Ashkan Panahi

Chalmers University of Technology Chalmers University of Technology

###### Abstract

We address a decentralized convex optimization problem, where every agent has
its unique local objective function and constraint set. Agents compute at
different speeds, and their communication may be delayed and directed. For
this setup, we propose an asynchronous double averaging and gradient
projection (ASY-DAGP) algorithm. Our algorithm handles difficult scenarios
such as message failure, by employing local buffers and utilizing the temporal
correlation in the transmitted messages. We guarantee the convergence speed of
our algorithm using performance estimation problems (PEP). In particular, we
introduce the concept of the linear quadratic (LQ) PEP. This approach
simplifies the analysis of smooth convex optimization problems, going beyond
Lyapunov function analyses and avoiding restrictive assumptions such as
strong-convexity. Numerical experiments validate the effectiveness of our
proposed algorithm.

##  1 Introduction

Consider $M$ computational agents exchanging information over a uni-
directional communication network, represented by a directed graph
$\mathcal{G}$. Their goal is to minimize the sum of local objective functions
$f^{v}(\mathbf{x})$ under an intersection of local constraints $S^{v}$, where
$f^{v},$ and $S^{v}$ are only known to node $v$ in $\mathcal{G}$. This can be
written as the following decentralized constrained optimization problem:

$\min_{\mathbf{x}\in\mathbb{R}^{m}}\;\;\frac{1}{M}\sum\limits_{v=1}^{M}f^{v}(\mathbf{x})\qquad\mbox{{s.t.}}\quad\mathbf{x}\in\bigcap\limits_{v=1}^{M}S^{v}.$
(P)

In decentralized optimization, each node executes a series of local
computations and communicate with its neighboring nodes. The processing
capacity of each node, together with the complexity of the local objective
terms, lead nodes to execute their computations at varying speeds.
Communications among agents may also encounter random delays or even message
losses, due to the inherent uncertainties in the underlying communication
networks. In this scenario, two different strategies can be considered, known
as synchronous and asynchronous.

Synchronous optimization methods require all agents to complete their
computations and communications before proceeding to the next iteration. The
speed of these algorithms is limited by the slowest node, called the
_straggler_. Additionally, these algorithms require a synchronization clock
coordinator, which is a demanding task in large-scale decentralized networks
(Hannah et al.,, 2018). Asynchronous decentralized methods eliminate the idle
time of agents by removing the synchronization points. In such methods, nodes
operate uninterruptedly, relying on the current information available from
their neighbors. Figure 1 illustrates the distinction between asynchronous and
synchronous methods.

In this paper, we consider the alternative of asynchronous algorithms as they
are computationally efficient. We address two major challenges in this setup.
First, agents must handle delayed information received from their neighbours.
Second, agents update at different rates, which results in difficulties with
convergence to the optimal solution of the underlying optimization problem
(Assran and Rabbat,, 2020).

$k=0$init.$k=1$node 1  node 2  node 3  init.node 1  node 2  node 3
CommunicatingComputation finishedComputingIdle

Figure 1: Synchronous setup (left): In each iteration, agents complete their
computations and transmit their updated variables to their neighbors. The
nodes wait until all messages are delivered to their destination. Then, they
simultaneously initiate the next iteration. Asynchronous setup (right): Each
agent completes its computation and sends the updated variables to its
neighbors. Immediately, it starts its next computation using available
information. There is a possibility of receiving multiple messages from one
node or sometimes no messages.

Algorithmic contribution: We introduce, for the first time, an asynchronous
algorithm – the asynchronous double averaging and gradient projection (ASY-
DAGP) – which effectively handles setups with local constraints and directed
graphs, ensures guaranteed convergence, and utilizes a fixed step size. This
algorithm is robust to both communication failures and latency, owing to the
inclusion of local buffers at each node. ASY-DAGP generalizes DAGP, introduced
by Shahriari-Mehr and Panahi, (2022), which is the only synchronous algorithm
addressing local constraints over directed graphs, to our knowledge.

Theoretical contribution: We establish that ASY-DAGP converges in terms of
optimality gap with $O(1/\sqrt{K})$ and its feasibility and consensus gaps
vanish with $O(1/K)$. It is well-known that this rate cannot be improved
unless stricter assumptions such as strong convexity are made. Moreover, we
identify the role of delay and asynchrony by a single factor, which we refer
to as _delay response_.

In our asynchronous, decentralized, and constrained setup, the standard
convergence analysis of optimization algorithms by Lyapunov (potential)
functions becomes extremely challenging. Merely assuming smooth convex
functions, we employ a novel proof technique, which draws its inspiration from
Shahriari-Mehr and Panahi, (2022). In particular, we show that Shahriari-Mehr
and Panahi, (2022) analysis is a systematic relaxation of the performance
estimation problems (PEP) (Drori and Teboulle,, 2014; Taylor,, 2017), which
can be used in generic convergence studies. We call this relaxation linear
quadratic PEP (LQ-PEP). In asynchronous setups, LQ-PEP has a great advantage
that it avoids prerequisite such as bounded delays, being a feature of
algorithms e.g. reliant on augmented graphs.

Paper organization: In Section 2, we present our problem setup and detail the
ASY-DAGP algorithm. Theoretical results based on LQ-PEP analysis, introduced
in Section 4, are presented in Section 3. Numerical experiments are provided
in section 5. Finally, we conclude with a literature review in section 6.

##  2 Methodology

In this section, we introduce our problem setup and present our algorithm,
ASY-DAGP.

###  2.1 Problem setup

We solve the decentralized constrained optimization problem (P) in a setup
with the following assumptions:

  1. 1.

The local objective functions $f^{v}$ are convex, differentiable, and
$L-$smooth, with $L>0$.

  2. 2.

The local constraint sets $S^{v}$ are closed and convex.

  3. 3.

The optimization problem is feasible and there exists a feasible optimal
solution $\mathbf{x}^{*}$ satisfying the sufficient optimality condition:

$\mathbf{0}\in\sum_{v=1}^{M}\left(\partial I_{S^{v}}(\mathbf{x}^{*})+\nabla
f^{v}(\mathbf{x}^{*})\right).$ (1)

  4. 4.

The communication network of the agents is given by a fixed connected directed
graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$, where the set of vertices
$\mathcal{V}$ represent agents, and
$\mathcal{E}\subseteq\mathcal{V}\times\mathcal{V}$ is the set of directed
edges. An edge $(v,u)\in\mathcal{E}$ represents a uni-directional
communication link from node $u$ to node $v$. Node $u$’s incoming and outgoing
neighbors are shown by
$\mathcal{N}^{u}_{\text{in}}=\left\\{v|(u,v)\in\mathcal{E}\right\\}$ and
$\mathcal{N}^{u}_{\text{out}}=\left\\{v|(v,u)\in\mathcal{E}\right\\}$,
respectively. The graph includes all self-loops; that is,
$u\in\mathcal{N}^{u}_{\text{in/out}}$. Accordingly, we consider two gossip
matrices $\mathbf{W}$ and $\mathbf{Q}$, with the same sparsity pattern as the
graph adjacency matrix. $\mathbf{W}$ has zero row sums, and $\mathbf{Q}$ has
zero column sums. They can be built e.g. by the input and output graph
Laplacian matrices.

  5. 5.

Communication links are subject to delays. When node $u$ transmits via edge
$(v,u)$, node $v$ receives the message after an arbitrary delay.

  6. 6.

Some messages may be lost during communication, i.e. a message transmitted
from node $u$ across the communication link $(v,u)$ may not arrive at the
intended destination, node $v$.

  7. 7.

Agents have varying computational power. Their computations occur at different
speeds. As a result, in a given time interval, each node performs a different
number of iterations.

DAGP is the only existing synchronous algorithm corresponding to the first
four assumptions, but cannot handle the rest. We aim to adapt DAGP for
asynchronous setups, addressing the last three assumptions.

###  2.2 ASY-DAGP algorithm

Similar to DAGP, ASY-DAGP attempts to solve (1) in an iterative and
distributed fashion, by splitting, i.e. re-writing it as:

$\sum_{v=1}^{M}\mathbf{g}^{v}=\mathbf{0},\quad\mathbf{g}^{v}\in\partial
I_{S^{v}}(\mathbf{x}^{*})+\nabla f^{v}(\mathbf{x}^{*}).$ (2)

the relation $\sum_{v=1}^{M}\mathbf{g}^{v}=\mathbf{0}$ is called _null
condition_. In DAGP, null condition is satisfied by the _distributed null
projection_ mechanism that uses another set of auxiliary variables
$\mathbf{h}^{v}$, always satisfying the null condition by design. The variable
$\mathbf{g}^{v}$ tracks $\mathbf{h}^{v}$, that is upon convergence
$\mathbf{g}^{v}$ becomes equal to $\mathbf{h}^{v}$, hence satisfying the null
condition. For comparison, DAGP is detailed in Algorithm 1.

In ASY-DAGP, we attempt to perform a similar procedure as in DAGP. However,
the asynchrony affects many steps, for which we introduce suitable
modifications. ASY-DAGP is presented in Algorithm 2. In this algorithm, agents
update their local variables in parallel with each other. When an agent
executes an iteration and sends the updated variables to its neighbors, it
instantly starts its next iteration. However, the duration of an iteration can
be variable. We introduce local buffers, $\mathcal{B}^{vu}$, to store all
messages node $v$ receives from $u$. After node $v$ processes its messages,
the buffer is cleared.

Although each node broadcasts the same message to all its out neighbours, they
receive and process that message at different iterations due to asynchrony.
This leads to a variable number of messages in a buffer. When there are one or
more messages in a buffer $\mathcal{B}^{vu}$, two averages $\mathbf{a}^{vu}$
and $\mathbf{b}^{vu}$ of them are calculated by (10) and (11). Otherwise, the
last calculated averages are reused. Comparing (3) and (7), as well as (6) and
(9), we observe that $(\mathbf{a}^{vu},\mathbf{b}^{vu})$ substitutes the
message pair $(\mathbf{x}^{v},\mathbf{h}^{v}-\mathbf{g}^{v})$ in DAGP. Our
idea is that due to temporal correlation in the transmitted messages, the
former tracks the latter and in effect ASY-DAGP functions similar to DAGP.

There is however a difficulty with the above strategy. Since, each node
calculates a distinct average of their buffer, in general
$\mathbf{b}^{vu}\neq\mathbf{b}^{wu}$, for some $v\neq w$. A consequence of
this phenomenon is that the variables $\mathbf{h}^{v}$ will not always satisfy
the null condition, and the distributed null projection strategy fails. To
address this, we introduce new vectors $\mathbf{p}^{v}$ and modify the update
rule of (6). In appendix 9, we show that convergence to a consensus solution
necessitates
$\mathbf{g}^{v}=\mathbf{h}^{v}=\frac{1}{\gamma-1}\sum_{u\in\mathcal{N}^{v}_{\mathrm{in}}}q_{vu}\mathbf{p}^{u}$.
Since $\mathbf{Q}$ has zero column sums, null condition is implied. As seen,
ASY-DAGP requires an additional variable $\mathbf{p}^{v}$, but an advantage
observed by the analysis of the fixed point in appendix 9 is that the
condition $\ker(\mathbf{Q})=\ker(\mathbf{W}^{T})$ in DAGP is no longer
necessary.

Algorithm 1 DAGP (Shahriari-Mehr and Panahi,, 2022).

0: step size $\mu$, scaling parameters $\alpha$ and $\rho$, and gossip
matrices $\mathbf{W}$ and $\mathbf{Q}$.

1: Initialize $k=0$. Initialize $\mathbf{x}^{v}_{0}$ randomly, and
$\mathbf{g}^{v}_{0}$ and $\mathbf{h}^{v}_{0}$ with zero vectors, $\forall
v\in\mathcal{V}$.

2: repeat

3: Update $\mathbf{z}^{v},\mathbf{x}^{v},\mathbf{g}^{v}$, and $\mathbf{h}^{v}$
variables, $\;\forall v\in\mathcal{V}$:  $\displaystyle\mathbf{z}^{v}_{k+1}$
$\displaystyle=\;$
$\displaystyle\mathbf{x}^{v}_{k}-\sum\limits_{u\in\mathcal{N}^{v}_{\text{in}}}w_{vu}\mathbf{x}^{u}_{k}-\mu\left(\nabla
f^{v}(\mathbf{x}^{v}_{k})-\mathbf{g}^{v}_{k}\right)$ (3)
$\displaystyle\mathbf{x}_{k+1}^{v}$ $\displaystyle=\;$
$\displaystyle\mathrm{P}_{S^{v}}\left(\mathbf{z}^{v}_{k+1}\right)$ (4)
$\displaystyle\mathbf{g}_{k+1}^{v}$ $\displaystyle=\;$
$\displaystyle\mathbf{g}_{k}^{v}+\rho(\nabla
f^{v}(\mathbf{x}^{v}_{k})-\mathbf{g}^{v}_{k})$
$\displaystyle\;+\frac{\rho}{\mu}\left(\mathbf{z}^{v}_{k+1}-\mathbf{x}^{v}_{k+1}\right)+\alpha\left(\mathbf{h}_{k}^{v}-\mathbf{g}_{k}^{v}\right)$
(5) $\displaystyle\mathbf{h}_{k+1}^{v}$ $\displaystyle=\;$
$\displaystyle\mathbf{h}_{k}^{v}-\sum\limits_{u\in\mathcal{N}^{v}_{\text{in}}}q_{vu}(\mathbf{h}_{k}^{u}-\mathbf{g}_{k}^{u})$
(6)

4: Send the updated tuple
$(\mathbf{x}^{v}_{k+1},\;\mathbf{h}^{v}_{k+1}-\mathbf{g}^{v}_{k+1})$ to all
out-neighbors $u\in\mathcal{N}^{v}_{\mathrm{out}}$, $\;\forall
v\in\mathcal{V}$.

5: Update iteration index: $k=k+1$.

6: until Convergence

Algorithm 2 ASY-DAGP.

0: step size $\mu$, scaling parameters $\alpha$, $\rho$, $\eta$, gossip
matrices $\mathbf{W}$ and $\mathbf{Q}$, and forgetting factor $\gamma$.

1: Initialize $\mathbf{x}^{v}$ randomly. Initialize $\mathbf{g}^{v}$,
$\mathbf{h}^{v}$, and $\mathbf{p}^{v}$ with zero vectors, $\;\forall
v\in\mathcal{V}.$

2: Initialize $\mathbf{b}^{vu}=\mathbf{0}$, $\forall v,u\in\mathcal{V}$,
$\mathbf{a}^{vv}=\mathbf{x}^{v}$, $\forall v\in\mathcal{V}$, and
$\mathbf{a}^{vu}=\mathbf{0}$, $\forall v,u\in\mathcal{V},u\neq v.$

3: Agents continuously receive information from their neighbors and store it
in their buffers $\mathcal{B}^{vu}$.

4: Agents update their variables in parallel as:

5: repeat

6: Update $\mathbf{z}^{v},\mathbf{p}^{v},$ and $\mathbf{h}^{v}$ variables,
$\displaystyle\mathbf{z}^{v}$ $\displaystyle=\;$
$\displaystyle\mathbf{x}^{v}-\sum\limits_{u\in\mathcal{N}^{v}_{\text{in}}}w_{vu}\mathbf{a}^{vu}-\mu\left(\nabla
f^{v}(\mathbf{x}^{v})-\mathbf{g}^{v}\right)$ (7) $\displaystyle\mathbf{p}^{v}$
$\displaystyle=\;$
$\displaystyle\mathbf{p}^{v}-\eta\sum\limits_{u\in\mathcal{N}^{v}_{\text{in}}}q_{vu}\mathbf{b}^{vu}+\eta(\gamma-1)\mathbf{g}^{v}$
(8) $\displaystyle\mathbf{h}^{v}$ $\displaystyle=\;$
$\displaystyle\gamma\mathbf{h}^{v}-\sum\limits_{u\in\mathcal{N}^{v}_{\text{in}}}q_{vu}\mathbf{b}^{vu}$
(9)

7: Update $\mathbf{x}^{v}$, $\mathbf{g}^{v}$ variables similar to (4)-(5).

8: Send the updated tuple $(\mathbf{x}^{v},\;\mathbf{p}^{v})$ to all out-
neighbors in $\mathcal{N}^{v}_{\mathrm{out}}$.

9: Update the estimates of $\mathbf{x}^{u}$ and $\mathbf{p}^{u}$ at node $v$
for every $u\in\mathcal{N}^{v}_{\text{in}}$. If the buffer $\mathcal{B}^{vu}$
is empty, reuse the existing value in memory; otherwise, compute the new
estimates as:
$\mathbf{a}^{vu}=\frac{1}{|\mathcal{B}^{vu}|}\sum\limits_{\mathbf{x}^{u}\in\mathcal{B}^{vu}}\mathbf{x}^{u}$
(10)
$\mathbf{b}^{vu}=\frac{1}{|\mathcal{B}^{vu}|}\sum\limits_{\mathbf{p}^{u}\in\mathcal{B}^{vu}}\mathbf{p}^{u}$
(11)

10: Clear the buffers $\mathcal{B}^{vu}$, for every
$u\in\mathcal{N}^{v}_{\text{in}}$.

11: until Convergence

##  3 Convergence guarantees

For the sake of analysis, we index the iterations of each node, individually.
Hence, iteration $k$ at different nodes happen at different real times.
However, it is harmless to pretend that each iteration occurs simultaneously
at every node, but the messages experience a corrected delay that may be
negative (non-causal). This assumption does not affect the generality of our
approach. Based on this convention, we introduce the local index sets
$T^{vu}_{k}$ containing all iteration numbers $k^{\prime}$, where a message
from $u$ is sent at $k^{\prime}$ and is processed by $v$ at time $k+1$. In
other words,
$T^{vu}_{k}=\\{k^{\prime}|k^{\prime}+\tau^{vu}_{k^{\prime}}=k+1\\}$.

Further, we need to introduce few definitions. First, we formalize a concept,
originally employed by Shahriari-Mehr and Panahi, (2022), which incorporates
frequency analysis tools into the study of optimization algorithms.

######  Definition 1 (Proper matrices).

A tuple $(\mathbf{S},\bar{\mathbf{R}},\mathbf{P})$ is _proper_ if, for any
given complex value $z$ and a positive real value $\beta$, all the statements
from Definition 5 in (Shahriari-Mehr and Panahi,, 2022) hold for the following
_forward-backward transfer matrix_ :

$\mathbf{F}_{\beta}(z)=\left[\begin{array}[]{cc}\mathbf{S}&\mathbf{I}-z^{-1}\bar{\mathbf{R}}^{T}\\\
\mathbf{I}-z\bar{\mathbf{R}}&\frac{1}{\beta}\mathbf{P}\mathbf{P}^{T}\end{array}\right]$
(12)

In particular, we denote by $z_{i}(\beta)$ the simple roots of
$\det(\mathbf{F}_{\beta}(z))$.

Next, we introduced a parameter quantifying the impact of delays and
asynchrony.

######  Definition 2 (Delay response).

For a nonempty index set $T_{k}^{vu}$, define

$\tau^{vu}_{i,k}(\beta)=\frac{1}{|T_{k}^{vu}|}\sum\limits_{k^{\prime}\in
T_{k}^{vu}}z_{i}^{|k+1-k^{\prime}|}(\beta)-1$ (13)

as the $i^{\text{th}}$ delay power spectrum of link $(v,u)$ at time $k$.
Accordingly, we define the delay response of link $(v,u)$ at time $m$ as

$\kappa^{vu}_{m}(\beta)=\sum\limits_{k\mid
T_{k}^{vu}\neq\emptyset,i}|z_{i}(\beta)|^{|k-m|}\tau^{vu}_{i,k}(\beta),$ (14)

Then, the delay response $\kappa$ is defined as the supremum of
$\kappa^{vu}_{m}(\beta)$ over every time $m$ and link $(v,u)$ and $\beta$.

Our main result is based on three matrices
$\mathbf{S},\bar{\mathbf{R}},\mathbf{P}$ which only depend on the parameters
of the algorithm and the network topology. Due to space limitation, they are
given in Appendix 10.

######  Theorem 1 (Rates of convergence).

Suppose that $(\mathbf{S},\bar{\mathbf{R}},\mathbf{P})$ is proper. There exist
a constant $C$ independent of the delays such that for $\kappa<\frac{1}{C}$
and a sufficiently large $K$, ASY-DAGP converges as follows:

$\|\bar{\mathbf{x}}^{v}_{K}-\bar{\mathbf{x}}_{K}\|_{2}^{2}\leq\frac{CC_{0}}{MK},$
(15)
$\left|\sum\limits_{v}\left(f^{v}(\bar{\mathbf{x}}^{v}_{K})-f^{v}(\mathbf{x}^{*})\right)\right|\leq\frac{CC_{0}}{\mu
K}+\sqrt{\frac{C_{2}CC_{0}}{MK}},$ (16)

where $K$ indicates that all nodes have performed $K$ iterations,
$\bar{\mathbf{x}}^{v}_{K}=\frac{1}{K}\sum_{k=0}^{K-1}\mathbf{x}^{v}_{k}$,
$\bar{\mathbf{x}}_{K}=\frac{1}{M}\sum_{v=1}^{M}\bar{\mathbf{x}}^{v}_{K}$.
$C_{0}$ is a constant dependent on the initial point, and
$C_{2}=\sqrt{\sum_{v}\|\mathbf{n}^{v}+\nabla f^{v}(\mathbf{x}^{*})\|^{2}}$,
measuring the heterogeneity of local functions.

######  Remark 1.

Note that by $K$ we mean that all nodes have performed $K$ iterations.
However, asynchrony could cause nodes to complete the $K\text{th}$ iteration
at various times. As a result, the convergence in real time hinges on the last
node completing its $K\text{th}$ iteration.

######  Remark 2.

The assumption that the tuple of matrices is proper is discussed in Remark 2
of (Shahriari-Mehr and Panahi,, 2022), where a similar transfer matrix is
studied. This study indicates that being proper matrices is not a restrictive
assumption, but dictates the range of parameters yielding convergence.

In the following section, we provide a sketch of the proof. The detailed
proof, the description of the matrices, along with their structures, are
provided in Appendix 10.

##  4 Proof outline of Theorem 1

A standard approach for studying the convergence speed of optimization
algorithms is the Lyapunov-based analysis (Polyak,, 1987), discussed in many
classical papers, e.g. (Nesterov,, 2012; Defazio et al.,, 2014; Schmidt et
al.,, 2017). This method relies on finding a suitable Lyapunov (potential)
function, which becomes difficult as the algorithm becomes complex. This is
especially the case for distributed and asynchronous optimization methods. The
DAGP analysis in Shahriari-Mehr and Panahi, (2022) is based on an alternative
approach, called aggregate lower bounding (ALB). In the following, we show the
ALB is a relaxation of the celebrated performance estimation problem (PEP),
introduced by Drori and Teboulle, (2014), which is a systematic ways to study
optimization algorithms. We refer to the relaxation underlying ALB as LQ-PEP,
based on which we provide convergence rates for ASY-DAGP.

LQ-PEP offers several advantages. It incorporates frequency analysis tools,
particularly Fourier analysis, into the study of optimization algorithms.
Additionally, LQ-PEP provides a tighter convergence bound than Lyapunov
analysis when employing the same inequalities in the analysis process.

###  4.1 Performance Estimation Problem

Let us first take a brief look at PEP. It evaluates the worst-case performance
of optimization algorithms after a fixed number $K$ iterations. An algorithm
consists of a state vector $\bm{\Psi}_{k}$ evolving over different iterations
$k=0,1,\ldots,K-1$, based on the underlying objective function $f$, and
concludes by returning a terminal state $\bar{\bm{\Psi}}_{K-1}$. Note that we
make a distinction between a state $\bm{\Psi}$ and the solution $\mathbf{X}$.
For our result in Theorem 1, a state consists of all variables of ASY-DAGP at
every node and the terminal state is the average state
$\bar{\bm{\Psi}}_{K-1}=\frac{1}{K}\sum_{k=1}^{K}\bm{\Psi}_{K}$. The
performance is quantified by a performance measure
$\Phi(\bar{\bm{\Psi}}_{K-1},\bm{\Psi}_{*})$, which is a positive-definite
function. For us, this is a combination of terms corresponding to optimality
gap, feasibility gap and consensus. PEP is defined as a search problem for the
worst performance over all objective functions $f$ from a given family
$\mathcal{F}$:

$\displaystyle\max_{f,\bm{\Psi}_{*},\\{\bm{\Psi}_{k}\\}_{k=0}^{K-1}}\quad$
$\displaystyle\Phi(\bar{\bm{\Psi}}_{K-1},\bm{\Psi}_{*})$ (17) s.t.
$\displaystyle f\ \text{be a function in family}\ \mathcal{F},$
$\displaystyle\bm{\Psi}_{*}\ \text{is the minimizer of}\ f,$
$\displaystyle\|\bm{\Psi}_{0}-\bm{\Psi}_{*}\|\leq R,$
$\displaystyle\text{algorithm generates}\
\bm{\Psi}_{1},\ldots,\bm{\Psi}_{K-1},$

where the dependence on $f$ often makes the problem infinite-dimensional.
Taylor et al., (2017) and Taylor, (2017) make an interesting observation that
the PEP in (17) can be reduced to a finite-dimensional problem (finite PEP).
The finite PEP is a complicated non-convex optimization. Therefore, its convex
relaxed versions are solved by Drori and Teboulle, (2014) and Taylor et al.,
(2017). Solving the relaxed finite PEPs naturally leads to upper bounds on the
speed of convergence, without any need for designing Lyaponov functions.

###  4.2 Linear Quadratic (LQ-)PEP

In Shahriari-Mehr and Panahi, (2022), the ALB approach is introduced, which
amounts to a novel finite-dimensional relaxation of (17) that we call LQ-PEP.
This relaxation is obtained in two steps of a) finding a quadratic upper bound
on the performance metric $\Phi(\bar{\bm{\Psi}}_{K-1},\bm{\Psi}_{*})$ and b) a
linear relaxation on the constraints in (17). In Appendix 10, we show that the
upper bound on the performance metric is given by

$\Phi(\bar{\bm{\Psi}}_{K-1},\bm{\Psi}_{*})\leq-\frac{1}{K}\sum_{k=0}^{K-1}\langle\tilde{\bm{\Psi}}_{k},\mathbf{S}\tilde{\bm{\Psi}}_{k}\rangle,$
(18)

where $\mathbf{S}$ is a fixed matrix depending on the topology of the network.
We also show that the following relation holds among the states generated by
the algorithm:

$\tilde{\bm{\bm{\Psi}}}_{k+1}-\bar{\mathbf{R}}\tilde{\bm{\bm{\Psi}}}_{k}-\mathbf{P}\tilde{\mathbf{U}}_{k}=\sum\limits_{l=0}^{K-2}\tilde{\mathbf{R}}_{k,k+1-l}\tilde{\bm{\bm{\Psi}}}_{l},$
(19)

where for simplicity we introduce
$\tilde{\bm{\bm{\Psi}}}_{k}\coloneqq\bm{\Psi}_{k}-\bm{\Psi}_{*}$ for a
suitable optimal state $\bm{\Psi}_{*}$, and $\tilde{\mathbf{U}}_{k}$ is a
sequence calculated by the states $\tilde{\bm{\bm{\Psi}}}_{k^{\prime}}$. Here,
$\mathbf{P}$ is a fixed matrix, $\bar{\mathbf{R}}$ only depends on the
parameters of the algorithm and $\tilde{\mathbf{R}}_{k,\tau}$ only depend on
the delay pattern. Without any delay, we have
$\tilde{\mathbf{R}}_{k,\tau}=\mathbf{0}$.

Then, LQ-PEP is given by the following result of the above statements:

######  Proposition 1.

Take matrices $\mathbf{S},\bar{\mathbf{R}},\mathbf{P}$ and
$\\{\tilde{\mathbf{R}}_{k,\tau}\\}$ as in (19), and a performance measure
$\Phi$ satisfying (18) and encompassing the optimality gap, the feasibility
gap and the consensus error. Then, the relation
$\Phi(\bar{\bm{\Psi}}_{K-1},\bm{\Psi}_{*})\leq\frac{C\|\tilde{\bm{\Psi}}_{0}\|^{2}}{K}$
holds if the optimal value of the following linear quadratic optimization is
zero.

$\displaystyle\min\limits_{\begin{subarray}{c}\\{\tilde{\bm{\Psi}}_{k}\\}_{k=0}^{K-1},\\\
{\\{\tilde{\mathbf{U}}_{k}\\}_{k=0}^{K-2}}\end{subarray}}\quad\frac{1}{2}\sum\limits_{k=0}^{K-1}\langle\tilde{\bm{\Psi}}_{k},\mathbf{S}\tilde{\bm{\Psi}}_{k}\rangle+\frac{C}{2}\|\tilde{\bm{\bm{\Psi}}}_{0}\|^{2}_{\mathrm{F}}$
(20) subject to
$\displaystyle\tilde{\bm{\bm{\Psi}}}_{k+1}-\bar{\mathbf{R}}\tilde{\bm{\bm{\Psi}}}_{k}-\mathbf{P}\tilde{\mathbf{U}}_{k}=\sum\limits_{l=0}^{K-2}\tilde{\mathbf{R}}_{k,k+1-l}\tilde{\bm{\bm{\Psi}}}_{l}.$
$\displaystyle k=0,1,\ldots,K-2$

The explicit construction of
$\mathbf{S},\bar{\mathbf{R}},\tilde{\mathbf{R}}_{k,l}$ and $\Phi$ is given in
Appendix 10. Note that $\mathbf{S}$ is an indefinite matrix and the condition
in (20) may only be achieved by considering the constraints. The constraints
in (20) represent a linear non-causal dynamical system, with
$\tilde{\mathbf{U}}_{k}$ as its input. If there are no delays, the dynamic
will become time invariant. Moreover, in this system, the existence of
negative delays is considered in the right hand side of the dynamics.

The proof of proposition (1) is straightforward. Suppose that the optimal
value is zero. Note that any sequence of states generated by the algorithm is
feasible and hence satisfies
$\frac{1}{2}\sum_{k=0}^{K-1}\langle\tilde{\bm{\Psi}}_{k},\mathbf{S}\tilde{\bm{\Psi}}_{k}\rangle+\frac{C}{2}\|\tilde{\bm{\bm{\Psi}}}_{0}\|^{2}_{\mathrm{F}}\geq
0$. By (18), we see that the result of proposition (1) holds.

###  4.3 Establishing Theorem 1

To obtain Theorem 1, we solve (20) and show that the conditions of LQ-PEP in
Proposition 1 hold. Here, we follow Shahriari-Mehr and Panahi, (2022) and use
the Lagrangian method of multipliers to show that this condition is equivalent
to the following system of recurrences having no non-zero solution:

$\displaystyle\tilde{\bm{\bm{\Psi}}}_{k+1}-\bar{\mathbf{R}}\tilde{\bm{\bm{\Psi}}}_{k}-\frac{1}{\beta}\mathbf{P}\mathbf{P}^{T}\bm{\Lambda}_{k}$
$\displaystyle=\sum\limits_{l=0}^{K-2}\tilde{\mathbf{R}}_{k,k+1-l}\tilde{\bm{\bm{\Psi}}}_{l}$
$\displaystyle\bm{\Lambda}_{k-1}-\bar{\mathbf{R}}^{T}\bm{\Lambda}_{k}+\mathbf{S}\tilde{\bm{\bm{\Psi}}}_{k}$
$\displaystyle=\sum\limits_{l=0}^{K-2}\tilde{\mathbf{R}}^{T}_{l,l-k+1}\bm{\Lambda}_{l},$

where $\bm{\Lambda}_{k}$s are the Lagrangian multipliers. The proof of this
fact and the rest of the proof of Theorem 1 is given in Appendix 11.

##  5 Experimental results

|
---|---
(a) Convergence over iterations. | (b) Convergence over time.

|
---|---
(c) Unconstrained setup. | (d) Setup with undirected graph.

Figure 2:  (a,b) Solving a constrained problem, and drawing a comparison to
DAGP and its throttled version. (c) Solving unconstrained logistic regression
problem, and comparing to APPG and ASY-SPA. (d) Solving a constrained problem
over undirected graphs, and comparing to ASY-PG-EXTRA.

In our experiments, we utilized directed and random networks, with edge
probability $p_{e}$. We assumed that the delays for each communication link
follow an exponential distribution with mean $\tau_{\text{comm}}^{vu}$. In the
first experiment, the computation time for each node followed a uniform
distribution within $[1,\bar{\tau}_{\text{comp}}^{v}]$. In the second one,
they were exponentially distribution with mean $\tau_{\text{comp}}^{v}$. We
used scaled Graph Laplacian matrices as gossip matrices, and the design
parameters were manually tuned for all algorithms. Each experiment was
conducted multiple times, but we present only a single instance due to the
negligible variations across separate runs. We report the objective values
computed at the average solution across all nodes. In our simulations, the
x-axis shows the global iteration counter, which increases when at least one
agent computes.

In every experiment, we simulated a distributed setup on a single machine by
sequentially performing the actions of different nodes. This approach gives
better control over parameters and facilitates the implementation of
asynchrony, delays, and message failures. As our experiments have an
illustrative nature, we postpone more realistic implementations to a future
study.

To the best of our knowledge, there is no other algorithm like ASY-DAGP that
handles a constrained optimization with asynchrony over directed graphs.
Hence, we compare ASY-DAGP to other existing algorithms in simplified setups,
which still demonstrate ASY-DAGP’s advantages. Specifically:

  * •

In the first experiment, we use a directed and constrained setup, but it is
synchronous.

  * •

In the second experiment, our setup is asynchronous and directed, but without
constraints.

  * •

The third experiment features an asynchronous setup that is constrained, but
undirected.

First experiment: In this experiment, we compare ASY-DAGP with DAGP. We
consider the following local objective functions and constraints

$\displaystyle
f^{v}(\mathbf{x})=\log\big{(}\cosh(\mathbf{a}_{v}^{T}\mathbf{x}-b_{v})\big{)},$
(21) $\displaystyle
S^{v}=\\{\mathbf{x}\;|\;\mathbf{c}_{v}^{T}\mathbf{x}-d_{v}\leq 0\\},$ (22)

where constants and coefficients are generated from normal distributions. We
initialize both algorithms randomly and assume
$M=10,\mathbf{x}\in\mathbb{R}^{5},p_{e}=0.8,\tau^{vu}_{\text{comp}}=10,\bar{\tau}^{v}_{\text{comp}}=5v,$
where $v$ indicates the node number. We set the design parameters to
$\rho=0.01,\alpha=0.1,\gamma=0.5,\eta=1.0.$ Figure 2a and 2b present the
resulting objective values. We observe both algorithms converge to the same
point, demonstrating that ASY-DAGP is capable of reaching the optimal solution
of (P). In terms of convergence over iterations, ASY-DAGP requires more
iterations. However, it achieves significantly faster wall-clock convergence.
It is worth noting that an increase in heterogeneity of the local computation
time (called throttling) - for instance, by deliberately slowing down several
nodes by a factor of two - may result in a slower convergence for the
synchronous setup. For ASY-DAGP, this situation is not harmful.

Second experiment: We compare ASY-DAGP to the APPG algorithm (Zhang and You,
2019b, ) in an asynchronous setup. We consider an unconstrained logistic
regression problem with $\ell_{2}$ regularization, similar to the one
presented in (Shahriari-Mehr and Panahi,, 2022). This problem is applied to
two digits of the MNIST dataset (LeCun et al.,, 2010). A random graph with
$p_{e}=0.6$ and $M=20$ nodes is considered. We use $N_{s}=1000$ samples from
the dataset and a regularization factor of $1/N_{s}$. The optimal solution is
computed by running the centralized gradient descent for a substantial number
of iterations. The optimality gap, i.e., the error in the objective function,
is reported in Figure 2c. In this experiment, the design parameters are set to
$\tau^{vu}_{\text{comm}}=50,\bar{\tau}^{v}_{\text{comp}}=5v,\gamma=0.5,\eta=1.0,\rho=0.1,\alpha=0.7,\mu=1.0.$
For both APPG and ASY-DAGP, we selected similar step-sizes slightly smaller
than $2/L$, provided they do not diverge. Although ASY-SPA typically requires
a diminishing step-size, we chose a small fixed one. With this choice, we may
only arrive at a neighborhood of the optimal solution.

Third experiment: In this experiment, we use the same setup as in the first
experiment, with the key difference of utilizing an undirected graph. We
compare ASY-DAGP with the ASY-PG-EXTRA algorithm in (Wu et al.,, 2017). For
both algorithms, we employ a near optimal step-size, and for the ASY-PG-EXTRA,
we set the relaxing parameter to $0.8$ for all agents. The results are shown
in Figure 2d. We plot the optimality gap with respect to the number of
completed communications. ASY-DAGP converges to the optimal solution with only
half the communications needed by ASY-PG-EXTRA to achieve the same result.

Communication failures: The goal of this experiment is to demonstrate the
effectiveness of ASY-DAGP in handling dropped messages. We consider the same
problem and parameters as in the first experiment. The messages can be lost
with a failure probability of $p$. Figure 3 investigates different values of
$p$. We observe that increasing the probability failure leads to more
iterations for convergence. We performed similar experiments on APPG, but it
failed to converge.

Figure 3: Robustness to message losses with the communication failure
probability of $p$.

The conducted experiments provide an evidence of ASY-DAGP’s efficiency under
diverse conditions. ASY-DAGP performs robustly even under extensive
communication failures, tolerating a loss of more than 75% of messages.
Additionally, ASY-DAGP surpasses its synchronous counterpart in wall-clock
convergence speed. Despite being originally designed for constrained setups
with directed communication networks, it can be successfully used for either
unconstrained problem scenarios or setups with undirected networks.

##  6 Related works

This review primarily focuses on decentralized gossip algorithms with general
directed communication networks. Given the real-world relevance, complex
dynamics, and asymmetric relationships inherent in directed graphs, they offer
a more comprehensive and realistic model for decentralized optimization than
undirected graphs. Notable papers, such as (Lian et al.,, 2018), (Peng et
al.,, 2016), (Wu et al.,, 2023), and (Even et al.,, 2021), have explored
undirected communication graphs in asynchronous setups. However, their
approach cannot be generalized to directed networks. In the following, we
review synchronous optimization algorithms, followed by their asynchronous
counterparts.

Synchronous optimization algorithms: Earlier papers consider the unconstrained
version of (P). For example, Nedić and Olshevsky, (2014, 2016) propose the
_subgradient-push_ algorithm based on the push-sum protocol (Kempe et al.,,
2003), which requires a vanishing step size. Taking a fixed step size leads to
an error proportional to the average dissimilarity between local objective
functions, a quantity related to _distribution shift_ in federated learning
literature (Reisizadeh et al.,, 2020; Fallah et al.,, 2020). In response, a
family of algorithms based on the so-called gradient tracking technique is
introduced (Xi and Khan,, 2017; Nedic et al.,, 2017). In particular, the
SONATA algorithm proposed by Scutari and Sun, (2019) combines the push-sum
protocol with a gradient tracking technique. The Push-Pull algorithm proposed
by Xin and Khan, (2018) and Pu et al., (2020) forgoes the push-sum protocol,
making it a simpler algorithm to analyse and implement. Both SONATA and Push-
Pull enjoy a linear convergence rate when the local objective functions are
smooth and strongly convex. Xi and Khan, (2016) study constraints by proposing
the DDPS algorithm which handles only identical constraints with diminishing
step sizes. The DAGP algorithm presented by Shahriari-Mehr and Panahi, (2022)
is the state-of-the-art, ensuring provable convergence for the generic case of
(P) with a fixed step size and directed graphs.

Asynchronous optimization algorithms: The study by Tsitsiklis et al., (1986)
is one of the first attempts to formulate the asynchronous setup, inspiring
many other algorithms. See the comprehensive survey by Assran et al., (2020).
Zhang and You, 2019a  and Assran and Rabbat, (2020) propose the asynchronous
version of the subgradient-push algorithm. They address the uneven update rate
issue, which leads to biased solutions. These methods require a decaying step
size similar to their synchronous versions. The asynchronous version of the
SONATA algorithm, ASY-SONATA, is proposed by Tian et al., (2020). They
introduce the perturbed sum-push protocol with a gradient tracking technique.
Finally, Zhang and You, 2019b  propose the fully asynchronous push-pull
gradient (APPG) algorithm. Both ASY-SONATA and APPG achieve a linear rate of
convergence using a constant step size: ASY-SONATA for strongly convex and
smooth problems, while APPG for smooth problems satisfying the Polyak-
Lojasiewicz condition.

Several papers address asynchrony and constraints by assuming composite
objective functions with non-smooth terms, notably (Wu et al.,, 2017) and
(Latafat and Patrinos,, 2022). These studies are based on the dual formulation
of the problem, leading to increased communication and computations due to the
incorporation of dual variables. The algorithm by Latafat and Patrinos, (2022)
is not completely asynchronous; it only permits agents to handle delayed
information. Wu et al., (2017) consider undirected graphs, and their algorithm
requires step-size adjustments, which demands the knowledge of the agents’
update rates — a stipulation not practical for real-world applications.
Furthermore, it lacks a guaranteed rate of convergence.

##  7 Conclusion

We proposed an asynchronous modification of the DAGP algorithm that
accommodates delays and missing messages, and is capable of solving smooth
optimization problems with local constraints over directed graphs. The
asynchronous updates eliminate idle time, thus leading to faster wall-clock
convergence. We presented a convergence analysis of our algorithm, based on a
relaxation of the PEP framework, which we refer to as LQ-PEP. It is closely
connected to the ALB approach in DAGP analysis. Using LQ-PEP, we summarized
the effect of the delay in a parameter that we referred to as delay response.
Lastly, our experimental results substantiate the resilience against dropped
messages, and show ASY-DAGP surpasses existing algorithms, even in their
specialized, restricted setups.

## References

  * Assran et al., (2020) Assran, M., Aytekin, A., Feyzmahdavian, H. R., Johansson, M., and Rabbat, M. G. (2020).  Advances in asynchronous parallel and distributed optimization.  Proceedings of the IEEE, 108(11):2013–2031.
  * Assran and Rabbat, (2020) Assran, M. S. and Rabbat, M. G. (2020).  Asynchronous gradient push.  IEEE Transactions on Automatic Control, 66(1):168–183.
  * Defazio et al., (2014) Defazio, A., Bach, F., and Lacoste-Julien, S. (2014).  Saga: A fast incremental gradient method with support for non-strongly convex composite objectives.  Advances in neural information processing systems, 27.
  * Drori and Teboulle, (2014) Drori, Y. and Teboulle, M. (2014).  Performance of first-order methods for smooth convex minimization: a novel approach.  Mathematical Programming, 145(1-2):451–482.
  * Even et al., (2021) Even, M., Hendrikx, H., and Massoulié, L. (2021).  Asynchronous speedup in decentralized optimization.  arXiv preprint arXiv:2106.03585.
  * Fallah et al., (2020) Fallah, A., Mokhtari, A., and Ozdaglar, A. (2020).  Personalized federated learning with theoretical guarantees: A model-agnostic meta-learning approach.  Advances in Neural Information Processing Systems, 33:3557–3568.
  * Hannah et al., (2018) Hannah, R., Feng, F., and Yin, W. (2018).  A2bcd: Asynchronous acceleration with optimal complexity.  In International Conference on Learning Representations.
  * Kempe et al., (2003) Kempe, D., Dobra, A., and Gehrke, J. (2003).  Gossip-based computation of aggregate information.  In 44th Annual IEEE Symposium on Foundations of Computer Science, 2003. Proceedings., pages 482–491. IEEE.
  * Latafat and Patrinos, (2022) Latafat, P. and Patrinos, P. (2022).  Primal-dual algorithms for multi-agent structured optimization over message-passing architectures with bounded communication delays.  Optimization Methods and Software, 37(6):2052–2079.
  * LeCun et al., (2010) LeCun, Y., Cortes, C., Burges, C., et al. (2010).  Mnist handwritten digit database.
  * Lian et al., (2018) Lian, X., Zhang, W., Zhang, C., and Liu, J. (2018).  Asynchronous decentralized parallel stochastic gradient descent.  In International Conference on Machine Learning, pages 3043–3052. PMLR.
  * Nedić and Olshevsky, (2014) Nedić, A. and Olshevsky, A. (2014).  Distributed optimization over time-varying directed graphs.  IEEE Transactions on Automatic Control, 60(3):601–615.
  * Nedić and Olshevsky, (2016) Nedić, A. and Olshevsky, A. (2016).  Stochastic gradient-push for strongly convex functions on time-varying directed graphs.  IEEE Transactions on Automatic Control, 61(12):3936–3947.
  * Nedic et al., (2017) Nedic, A., Olshevsky, A., and Shi, W. (2017).  Achieving geometric convergence for distributed optimization over time-varying graphs.  SIAM Journal on Optimization, 27(4):2597–2633.
  * Nesterov, (2012) Nesterov, Y. (2012).  Efficiency of coordinate descent methods on huge-scale optimization problems.  SIAM Journal on Optimization, 22(2):341–362.
  * Peng et al., (2016) Peng, Z., Xu, Y., Yan, M., and Yin, W. (2016).  Arock: an algorithmic framework for asynchronous parallel coordinate updates.  SIAM Journal on Scientific Computing, 38(5):A2851–A2879.
  * Polyak, (1987) Polyak, B. T. (1987).  Introduction to optimization. optimization software.  Inc., Publications Division, New York, 1:32.
  * Pu et al., (2020) Pu, S., Shi, W., Xu, J., and Nedić, A. (2020).  Push–pull gradient methods for distributed optimization in networks.  IEEE Transactions on Automatic Control, 66(1):1–16.
  * Reisizadeh et al., (2020) Reisizadeh, A., Farnia, F., Pedarsani, R., and Jadbabaie, A. (2020).  Robust federated learning: The case of affine distribution shifts.  Advances in Neural Information Processing Systems, 33:21554–21565.
  * Schmidt et al., (2017) Schmidt, M., Le Roux, N., and Bach, F. (2017).  Minimizing finite sums with the stochastic average gradient.  Mathematical Programming, 162:83–112.
  * Scutari and Sun, (2019) Scutari, G. and Sun, Y. (2019).  Distributed nonconvex constrained optimization over time-varying digraphs.  Mathematical Programming, 176:497–544.
  * Shahriari-Mehr and Panahi, (2022) Shahriari-Mehr, F. and Panahi, A. (2022).  Double averaging and gradient projection: Convergence guarantees for decentralized constrained optimization.
  * Taylor, (2017) Taylor, A. B. (2017).  Convex interpolation and performance estimation of first-order methods for convex optimization. PhD thesis, Catholic University of Louvain, Louvain-la-Neuve, Belgium.
  * Taylor et al., (2017) Taylor, A. B., Hendrickx, J. M., and Glineur, F. (2017).  Smooth strongly convex interpolation and exact worst-case performance of first-order methods.  Mathematical Programming, 161:307–345.
  * Tian et al., (2020) Tian, Y., Sun, Y., and Scutari, G. (2020).  Achieving linear convergence in distributed asynchronous multiagent optimization.  IEEE Transactions on Automatic Control, 65(12):5264–5279.
  * Tsitsiklis et al., (1986) Tsitsiklis, J., Bertsekas, D., and Athans, M. (1986).  Distributed asynchronous deterministic and stochastic gradient optimization algorithms.  IEEE transactions on automatic control, 31(9):803–812.
  * Wu et al., (2017) Wu, T., Yuan, K., Ling, Q., Yin, W., and Sayed, A. H. (2017).  Decentralized consensus optimization with asynchrony and delays.  IEEE Transactions on Signal and Information Processing over Networks, 4(2):293–307.
  * Wu et al., (2023) Wu, X., Liu, C., Magnusson, S., and Johansson, M. (2023).  Delay-agnostic asynchronous distributed optimization.  arXiv preprint arXiv:2303.18034.
  * Xi and Khan, (2016) Xi, C. and Khan, U. A. (2016).  Distributed subgradient projection algorithm over directed graphs.  IEEE Transactions on Automatic Control, 62(8):3986–3992.
  * Xi and Khan, (2017) Xi, C. and Khan, U. A. (2017).  Dextra: A fast algorithm for optimization over directed graphs.  IEEE Transactions on Automatic Control, 62(10):4980–4993.
  * Xin and Khan, (2018) Xin, R. and Khan, U. A. (2018).  A linear algorithm for optimization over directed graphs with geometric convergence.  IEEE Control Systems Letters, 2(3):315–320.
  * (32) Zhang, J. and You, K. (2019a).  Asyspa: An exact asynchronous algorithm for convex optimization over digraphs.  IEEE Transactions on Automatic Control, 65(6):2494–2509.
  * (33) Zhang, J. and You, K. (2019b).  Fully asynchronous distributed optimization with linear convergence in directed networks.  arXiv preprint arXiv:1901.08215.

##  Supplementary Material

The equation numbers and references are in agreement with the main body of the
paper.

##  8 Mathematical Notation

In this paper, we denote vectors using bold lowercase letters and matrices
with bold uppercase letters. The element at the $v^{\text{th}}$ row and
$u^{\text{th}}$ column of matrix $\mathbf{W}$ is denoted by $w_{vu}$. The
transpose of $\mathbf{W}$ is $\mathbf{W}^{T}$, and its right null space is
shown by $\ker(\mathbf{W})$. Therefore, $\mathbf{x}$ belongs to
$\ker(\mathbf{W})$ if and only if $\mathbf{W}\mathbf{x}=\mathbf{0}$. We use
$\mathbf{1}_{n}$ and $\mathbf{0}_{n}$ to respectively denote $n-$dimensional
vectors entirely consisting of ones and zeros. Furthermore, an $m\times n$
matrix containing only zero elements is shown by $\mathbf{O}_{m\times n}$. The
indices $m$ and $n$ may be omitted if there is no risk of confusion. The
Euclidean inner product of vectors is denoted by $\langle.,.\rangle$, and the
matrix inner product is denoted by
$\langle\mathbf{A},\mathbf{C}\rangle=\text{Tr}(\mathbf{A}\mathbf{C}^{T})$,
where $\text{Tr}(\ldotp)$ signifies the matrix trace operator. The Kronecker
delta function is shown by $\delta_{k,l}$. $\partial I_{S}$ indicates the
normal cone of $S$. $|T|$ indicates the cardinality of set $T$.

In this paper, we generally use subscripts to define the iteration number, and
superscripts to denote the node number. For example, $\nabla
f^{v}(\mathbf{x}^{v}_{k})$ represents the gradient of the local function of
node $v$ at its local variable at the $k^{\text{th}}$ iteration. Additionally,
to provide matrix representations of equations, we arrange vector variables as
the rows of a matrix. For instance, the matrix
$\mathbf{G}\in\mathbb{R}^{M\times m}$ contains all
$\mathbf{g}^{v}\in\mathbb{R}^{m}$ vectors, with $v\in\\{1,\ldots,M\\}$, as its
rows. For the sake of simplicity, in this paper, the statement
$\mathbf{G}\in\ker(\mathbf{W})$ implies that each column of $\mathbf{G}$ is an
element in the null space of $\mathbf{W}$.

##  9 Analysis of fixed points

The fixed point analysis is a simple way to study the convergence properties
of optimization algorithms. In DAGP, the extra assumption
$\ker(\mathbf{Q})=\ker(\mathbf{W}^{T})$ is required to show any fixed point of
DAGP is a consensus and optimal solution of problem (P). However, In ASY-DAGP,
by defining new $\mathbf{p}^{v}$ variables and designing their update
equation, defined in (8), based on the satisfaction of the optimality
condition for the fixed point, there is no need for this extra assumption on
gossip matrices anymore. For the fixed point iteration of ASY-DAGP, we have
$\mathbf{x}^{v}_{k+1}=\mathbf{x}^{v}_{k}=\mathbf{x}^{v}$. In addition, the
same relation holds true for $\mathbf{g}^{v},\mathbf{p}^{v},$ and
$\mathbf{h}^{v}$ variables. Since the algorithm is in the fixed point regime,
we have additional relation
$\mathbf{a}^{vu}_{k+1}=\mathbf{a}^{vu}_{k}=\mathbf{x}^{u}$, and
$\mathbf{b}^{vu}_{k+1}=\mathbf{b}^{vu}_{k}=\mathbf{p}^{u}.$ Let us look at the
fixed point iteration of ASY-DAGP in the following matrix form:

$\displaystyle\mathbf{Z}=\mathbf{X}-\mathbf{W}\mathbf{X}-\mu(\bm{\nabla}\mathbf{F}-\mathbf{G})$
(23) $\displaystyle\mathbf{X}=\mathbf{P}(\mathbf{Z})$ (24)
$\displaystyle\alpha(\mathbf{H}-\mathbf{G})+\rho\left(\bm{\nabla}\mathbf{F}-\mathbf{G}+\frac{1}{\mu}(\mathbf{Z}-\mathbf{X})\right)=\mathbf{O}$
(25) $\displaystyle\mathbf{Q}\mathbf{P}=(\gamma-1)\mathbf{G}$ (26)
$\displaystyle\mathbf{Q}\mathbf{P}=(\gamma-1)\mathbf{H}$ (27)

From the last two equations, $\mathbf{H}=\mathbf{G}$. Then, from (25):

$\mathbf{G}=\bm{\nabla}\mathbf{F}+\frac{1}{\mu}(\mathbf{Z}-\mathbf{X}).$ (28)

This relation proves that the fixed point is the consensus and optimal
solution of (P). Using (28) in (23), we have
$\mathbf{W}\mathbf{X}=\mathbf{O}$, which corroborate a consensus solution,
i.e. $\mathbf{x}^{v}=\mathbf{x}^{u}=\mathbf{x}$ for all $v,u\in\mathcal{V}$.
To show $\mathbf{x}$ is an optimal solution, first note that
$\mathbf{1}^{T}\mathbf{G}=\mathbf{0}^{T}$ from (26). Moreover, consider that
$\mathbf{z}^{v}-\mathbf{x}^{v}\in\partial I_{S^{v}}$, for all
$v\in\mathcal{V}.$ Then, by left multiplying (28) with $\mathbf{1}^{T}$,
considering the conic property of normal cone, we have the optimality
condition satisfied at $\mathbf{x}$. By this analysis, we showed that ASY-DAGP
does not require the limiting kernel assumption on gossip matrices.

##  10 Convergence analysis

In this section, we follow the steps of ALB, similar to (Shahriari-Mehr and
Panahi,, 2022), to derive LQ-PEP.

As mentioned in the paper, for the sake of analysis, we index the iterations
of each node, individually. Hence, iteration $k$ at different nodes happen at
different real times. However, we pretend that each iteration occurs
simultaneously at every node, but the messages experience a corrected delay
that may be negative (non-causal). This assumption does not affect the
generality of our approach. Based on this convention, we introduce the local
index sets $T^{vu}_{k}$ containing all iteration numbers $k^{\prime}$, where a
message from $u$ is sent at $k^{\prime}$ and is processed by $v$ at time
$k+1$. In other words,
$T^{vu}_{k}=\\{k^{\prime}|k^{\prime}+\tau^{vu}_{k^{\prime}}=k+1\\}$.

### Step one: combining inequalities from the assumptions

First, we define

$\displaystyle F^{v}(\mathbf{x})$ $\displaystyle\coloneqq
f^{v}(\mathbf{x})-f^{v}(\mathbf{x}^{*})-\langle\nabla
f^{v}(\mathbf{x}^{*}),\mathbf{x}-\mathbf{x}^{*}\rangle,$ $\displaystyle
T^{v}(\mathbf{x})$
$\displaystyle\coloneqq-\langle\mathbf{n}^{v},\mathbf{x}-\mathbf{x}^{*}\rangle,$

where $\mathbf{n}^{v}\in\partial I_{S^{v}}(\mathbf{x}^{*})$, i.e.
$\mathbf{n}^{v}$ is an element in the normal cone of $S^{v}$ at
$\mathbf{x}^{*}$. From the convexity assumption on $f^{v}$ and the definition
of normal cone, we have

$\displaystyle F^{v}(\mathbf{x}^{v})$ $\displaystyle\geq 0,$
$\displaystyle\qquad\forall\mathbf{x}$ (29) $\displaystyle
T^{v}(\mathbf{x}^{v})$ $\displaystyle\geq 0.$
$\displaystyle\qquad\forall\mathbf{x}\in S^{v}$ (30)

Hence, these functions may serve as performance measures for our decentralized
and constrained setup. In the following, we define
$F^{v}_{k+1}=F^{v}(\mathbf{x}^{v}_{k+1})$ and
$T^{v}_{k+1}=T^{v}(\mathbf{x}^{v}_{k+1})$ for notation simplicity.

Form convexity of $f^{v}$, we have

$F^{v}_{k}+\left\langle\nabla f^{v}(\mathbf{x}^{*})-\nabla
f^{v}(\mathbf{x}^{v}_{k}),\mathbf{x}^{v}_{k}-\mathbf{x}^{*}\right\rangle\leq
0.$ (31)

From the $L-$smoothness property of $f^{v}$, we have

$F_{k+1}^{v}-F_{k}^{v}-\langle\nabla f^{v}(\mathbf{x}_{k}^{v})-\nabla
f^{v}(\mathbf{x}^{*}),\mathbf{x}_{k+1}^{v}-\mathbf{x}^{v}_{k}\rangle\leq 0.$
(32)

From the projection step in (7), and $\mathbf{x}^{v}_{k+1}\in S^{v}$, we have

$\mu
T^{v}_{k+1}+\langle\mathbf{x}^{*}-\mathbf{x}^{v}_{k+1},\mathbf{z}_{k+1}^{v}-\mathbf{x}^{v}_{k+1}-\mu\mathbf{n}^{v}\rangle\leq
0.$ (33)

By summing the first two inequalities, multiplying the resulting expression by
$\mu$, and then adding it to (33), we obtain

$\displaystyle\mu(F^{v}_{k+1}+T^{v}_{k+1})$
$\displaystyle-\frac{L\mu}{2}\left\|\mathbf{x}^{v}_{k+1}-\mathbf{x}^{v}_{k}\right\|$
$\displaystyle+\Big{\langle}\mathbf{x}^{*}-\mathbf{x}^{v}_{k+1},\mathbf{z}^{v}_{k+1}-\mathbf{x}^{v}_{k+1}+\mu\big{(}\nabla
f^{v}(\mathbf{x}^{v}_{k})-\nabla
f^{v}(\mathbf{x}^{*})-\mathbf{n}^{v}\big{)}\Big{\rangle}\leq 0.$ (34)

### Step two: plugging algorithm’s dynamics

The objective of this step is to eliminate the gradient term from (34). By
substituting the definition of $\mathbf{z}^{v}_{k+1}$ from (7) into (34), we
obtain

$\displaystyle\mu(F_{k+1}^{v}$
$\displaystyle+T_{k+1}^{v})-\frac{L\mu}{2}\left\|\mathbf{x}_{k+1}^{v}-\mathbf{x}_{k}^{v}\right\|^{2}$
$\displaystyle+\Big{\langle}\mathbf{x}^{*}-\mathbf{x}^{v}_{k+1},\mathbf{x}_{k}^{v}-\sum\limits_{u}w_{vu}\mathbf{a}^{vu}_{k}-\mathbf{x}^{v}_{k+1}+\mu(\mathbf{g}^{v}_{k}-\nabla
f^{v}(\mathbf{x}^{*})-\mathbf{n}^{v})\Big{\rangle}\leq 0.$ (35)

Furthermore, to remove the algorithm dynamics’ dependency on
$\mathbf{z}^{v}_{k+1}$, we can plug the definition of $\mathbf{z}^{v}_{k+1}$
into the update rule or the dynamics of $\mathbf{g}^{v}_{k}$ as shown in (5).
This gives us

$\mathbf{g}_{k+1}^{v}=\mathbf{g}_{k}^{v}+\frac{\rho}{\mu}\left(\mathbf{x}_{k}^{v}-\sum\limits_{u}w_{vu}\mathbf{a}^{vu}_{k}-\mathbf{x}^{v}_{k+1}\right)+\alpha(\mathbf{h}^{v}_{k}-\mathbf{g}^{v}_{k}).$
(36)

Now, for simplicity in analysis and notation, we define
$\tilde{\mathbf{x}}_{k}^{v}\coloneqq\mathbf{x}^{v}_{k}-\mathbf{x}^{*}$,
$\tilde{\mathbf{g}}^{v}_{k}\coloneqq\mathbf{g}^{v}_{k}-(\nabla
f^{v}(\mathbf{x}^{*})+\mathbf{n}^{v})$,
$\tilde{\mathbf{h}}^{v}_{k}=\mathbf{h}^{v}_{k}-(\nabla
f^{v}(\mathbf{x}^{*})+\mathbf{n}^{v})$, and
$\tilde{\mathbf{p}}^{v}_{k}=\mathbf{p}^{v}_{k}-\frac{\gamma-1}{\sum_{u}q_{vu}}(\nabla
f^{v}(\mathbf{x}^{*})+\mathbf{n}^{v})$. This transformation adjusts the
system’s origin to the optimal solution, leading to a modification of the
equation in (35) to

$\mu\left(F_{k+1}^{v}+T_{k+1}^{v}\right)-\frac{L\mu}{2}\left\|\tilde{\mathbf{x}}_{k+1}^{v}-\tilde{\mathbf{x}}_{k}^{v}\right\|^{2}-\Big{\langle}\tilde{\mathbf{x}}^{v}_{k+1},\tilde{\mathbf{x}}_{k}^{v}-\sum\limits_{u}w_{vu}\tilde{\mathbf{a}}^{vu}_{k}-\tilde{\mathbf{x}}^{v}_{k+1}+\mu\tilde{\mathbf{g}}\Big{\rangle}\leq
0.$ (37)

Similarly, the introduction of new variables transforms the update rules for
the algorithm’s variables into

$\displaystyle\tilde{\mathbf{g}}_{k+1}^{v}$
$\displaystyle=\tilde{\mathbf{g}}_{k}^{v}+\frac{\rho}{\mu}\left(\tilde{\mathbf{x}}_{k}^{v}-\sum\limits_{u}w_{vu}\tilde{\mathbf{a}}^{vu}_{k}-\tilde{\mathbf{x}}^{v}_{k+1}\right)+\alpha(\tilde{\mathbf{h}}^{v}_{k}-\tilde{\mathbf{g}}^{v}_{k}),$
(38) $\displaystyle\tilde{\mathbf{p}}^{v}_{k+1}$
$\displaystyle=\tilde{\mathbf{p}}^{v}_{k}-\eta\sum\limits_{u}q_{vu}\tilde{\mathbf{b}}^{vu}_{k}+\eta(\gamma-1)\tilde{\mathbf{g}}^{v}_{k},$
(39) $\displaystyle\tilde{\mathbf{h}}_{k+1}^{v}$
$\displaystyle=\gamma\tilde{\mathbf{h}}_{k}^{v}-\sum\limits_{u}q_{vu}\tilde{\mathbf{b}}^{vu}_{k},$
(40)
$\tilde{\mathbf{a}}^{vu}_{k+1}=\left\\{\begin{array}[]{lc}\tilde{\mathbf{a}}^{vu}_{k}&\text{If
the buffer is empty. }\\\ \frac{1}{|T^{vu}_{k}|}\sum\limits_{k^{\prime}\in
T^{vu}_{k}}\tilde{\mathbf{x}}^{u}_{k^{\prime}}&\text{ Otherwise
}\end{array}\right.$ (41)
$\tilde{\mathbf{b}}^{vu}_{k+1}=\left\\{\begin{array}[]{lc}\tilde{\mathbf{b}}^{vu}_{k}&\text{If
the buffer is empty }\\\ \frac{1}{|T^{vu}_{k}|}\sum\limits_{k^{\prime}\in
T^{vu}_{k}}\tilde{\mathbf{p}}^{u}_{k^{\prime}}&\text{ Otherwise
}\end{array}\right.$ (42)

By defining the state vectors for ASY-DAGP as
$\mathbf{\Psi}_{k}\in\mathbb{R}^{(2M^{2}+6M)\times m}$

$\bm{\Psi}_{k}=\left[\mathbf{X}_{k+1}^{T}\;\;\mathbf{P}_{k+1}^{T}\;\;\mathbf{X}_{k}^{T}\;\;\mathbf{G}^{T}_{k}\;\;\mathbf{H}^{T}_{k}\;\;\mathbf{P}_{k}^{T}\;\;\\{\mathbf{A}_{v,k}^{T}\\}_{v=1}^{M}\;\;\\{\mathbf{B}_{v,k}^{T}\\}_{v=1}^{M}\right]^{T},$
(43)

and introducing the shifted states
$\tilde{\mathbf{\Psi}}_{k}=\mathbf{\Psi}_{k}-\mathbf{\Psi}_{*}$, where
$\mathbf{\Psi}_{*}$ represents the optimal states, we can describe the
evolution of ASY-DAGP’s variables using the subsequent non-causal linear
dynamical system

$\tilde{\bm{\bm{\Psi}}}_{k+1}-\bar{\mathbf{R}}\tilde{\bm{\bm{\Psi}}}_{k}-\mathbf{P}\tilde{\mathbf{U}}_{k}=\sum\limits_{l=0}^{K-2}\tilde{\mathbf{R}}_{k,k+1-l}\tilde{\bm{\bm{\Psi}}}_{l},\qquad
k=0,\ldots,K-2$ (44)

where $\bar{\mathbf{R}},\tilde{\mathbf{R}}_{k,k+1-l}$, and $\mathbf{P}$
matrices are defined in (52), (53), and (54), respectively. The matrix
$\tilde{\mathbf{U}}_{k}\coloneqq[\tilde{\mathbf{X}}_{k+2}^{T}\;\;\tilde{\mathbf{P}}_{k+2}^{T}]^{T}$
represents the input of the system.

Regarding the creation of $\bar{\mathbf{R}}$ and $\tilde{\mathbf{R}}_{k,0}$
matrices, we follow the approach that the ideal case, where there are no
communication delays and the agents have instant access to the variables of
other nodes, should be encapsulated within the definition of
$\bar{\mathbf{R}}$. By doing so, in scenarios without delays, the ASY-DAGP
analysis simplifies to the DAGP analysis. Therefore, we put all the delay-
independent terms in (38), (39), and (40) in the definition of
$\bar{\mathbf{R}}$. Also, in the creation of $\bar{\mathbf{R}}$ the blocks
representing update rules (41) and (42) are created based on the ideal case,
i.e. $\tilde{\mathbf{a}}^{vu}_{k+1}=\tilde{\mathbf{x}}^{u}_{k+1}$ and
$\tilde{\mathbf{b}}^{vu}_{k+1}=\tilde{\mathbf{p}}^{u}_{k+1}$. Since the
coefficient of $\tilde{\bm{\bm{\Psi}}}_{k}$ is
$\bar{\mathbf{R}}+\tilde{\mathbf{R}}_{k,0}$, and we put the ideal case in
$\bar{\mathbf{R}}$, the $\tilde{\mathbf{R}}_{k,0}$ matrix should compensate
for non-ideal cases. In $\tilde{\mathbf{R}}_{k,k+1-l}$, only the blocks
associated with $\tilde{\mathbf{A}}_{v,k}$ and $\tilde{\mathbf{B}}_{v,k}$ are
non-zero. Given the similarity of update rules in (41) and (42), we focus on
how we create the blocks representing the update rule of
$\tilde{\mathbf{A}}_{v,k}$. For clarity, we consider only one node $v$ and the
$u\text{th}$ row in $\tilde{\mathbf{A}}_{v,k}$. For simplicity in notation, we
define the _interaction element_ as the possible non-zero element connecting
$\tilde{\mathbf{a}}^{vu}_{k+1}$ and $\tilde{\mathbf{x}}^{u}_{k+1}$, which is
in the $u^{\text{th}}$ column of the matrix $\tilde{\mathbf{R}}_{k,k+1-l}$.
When computing $\tilde{\mathbf{R}}_{k,0}$, several scenarios arise:

  * •

If the buffer is empty and the ideal message is not present, interaction
element is $-1$. Additionally, the element linking $\mathbf{a}^{vu}_{k+1}$ and
$\mathbf{a}^{vu}_{k}$ is $+1$.

  * •

If several messages exist in the buffer and $k+1\in T^{vu}_{k}$, i.e. ideal
message is present in the buffer, the interaction element computes as
$-1+\frac{1}{|T^{vu}_{k}|}$.

  * •

If several messages are in the buffer, but $k+1$ is not in the indices set
$T^{vu}_{k}$, the interaction element is $-1$.

For the computation of other $\tilde{\mathbf{R}}_{k,k+1-l},\;l\neq k+1$, the
only case when an element can be non-zero is when there are multiple messages
in the buffer, excluding the ideal message. In this scenario, the element
connecting $\tilde{\mathbf{a}}^{vu}_{k+1}$ and $\tilde{\mathbf{x}}^{u}_{k+1}$
is $\frac{1}{|T^{vu}_{k}|}$ for all $l\in T^{vu}_{k}$.

### Step three: defining performance measures

As mentioned earlier, $F^{v}$ and $T^{v}$ are potential performance measures
for the decentralized constrained setup. However, in decentralized setups,
convergence to a consensus solution is of interest as well. Therefore, we add
and remove
$\frac{\zeta}{2}\sum_{u,v}\|\tilde{\mathbf{x}}_{k+1}^{u}-\tilde{\mathbf{x}}_{k+1}^{v}\|^{2}$
to (37). Then, by summing over all $v=1,\ldots,M$ and $k=0,\ldots,K-1$, we
have

$\displaystyle\sum\limits_{k=0}^{K-1}\left(\sum\limits_{v=1}^{M}\mu\left(F_{k+1}^{v}+T_{k+1}^{v}\right)+\frac{\zeta}{2}\sum\limits_{u,v=1}^{M}\|\tilde{\mathbf{x}}_{k+1}^{u}-\tilde{\mathbf{x}}_{k+1}^{v}\|^{2}\right)$
$\displaystyle-\frac{L\mu}{2}\sum\limits_{k=0}^{K-1}\sum\limits_{v=1}^{M}\left\|\tilde{\mathbf{x}}_{k+1}^{v}-\tilde{\mathbf{x}}_{k}^{v}\right\|^{2}$
$\displaystyle-\sum\limits_{k=0}^{K-1}\sum\limits_{v=1}^{M}\Big{\langle}\tilde{\mathbf{x}}^{v}_{k+1},\tilde{\mathbf{x}}_{k}^{v}-\sum\limits_{u}w_{vu}\tilde{\mathbf{a}}^{vu}_{k}-\tilde{\mathbf{x}}^{v}_{k+1}+\mu\tilde{\mathbf{g}}\Big{\rangle}$
$\displaystyle-\frac{\zeta}{2}\sum\limits_{k=0}^{K-1}\sum\limits_{u,v=1}^{M}\|\tilde{\mathbf{x}}_{k+1}^{u}-\tilde{\mathbf{x}}_{k+1}^{v}\|^{2}\leq
0.$ (45)

The first line in (10) includes all terms that measure the algorithm’s
performance over all iterations, represented by
$\sum_{k}\Phi(\bm{\Psi}_{k},\bm{\Psi}_{*})$. In the Performance Estimation
Problem (PEP), our goal is to find the worst-case performance. We define all
terms in the last three summations in (10) as $\bar{A}_{K}$, which can be
written in the following quadratic form

$\bar{A}_{K}=\sum_{k=0}^{K-1}\langle\tilde{\bm{\bm{\Psi}}}_{k},\mathbf{S}\tilde{\bm{\bm{\Psi}}}_{k}\rangle,$
(46)

where matrix $\mathbf{S}$ is defined in (55). Now, by applying the new
definitions into (10), we have

$\sum\limits_{k=0}^{K-1}\Phi(\bm{\Psi}_{k},\bm{\Psi}_{*})\leq-\bar{A}_{K}.$
(47)

Hence, $\bar{A}_{K}$ is an upper bound for the performance measure, and the
objective function in PEP can be replaced by $-\bar{A}_{K}$. We then
substitute the problem of maximizing $\Phi$ with the problem of minimizing
$\bar{A}_{K}$, which has a quadratic form. We also consider the constraint as
the linear dynamical system in (44), which contain all possible states created
by the optimization algorithm. We call the new problem LQ-PEP. Upon reviewing
the derivation of (10), we see that it uses a set of
$\mathcal{F}-$interpolation inequalities introduced by Taylor et al., (2017),
representing the set of convex and smooth functions. Therefore, LQ-PEP is just
a relaxation of PEP.

### Step four: solving LQ-PEP

To establish that $\bar{A}_{K}$ is lower bounded by
$-C\|\tilde{\bm{\bm{\Psi}}}_{0}\|_{\mathrm{F}}^{2}$, which is essentially the
solution to LQ-PEP, we need to show that zero is the solution to the following
optimization problem.

$\displaystyle\min\limits_{\\{\tilde{\bm{\Psi}}_{k}\\}_{k=0}^{K-1},\\{\tilde{\mathbf{U}}_{k}\\}_{k=0}^{K-2}}\;\;\sum\limits_{k=0}^{K-1}\langle\tilde{\bm{\Psi}}_{k},\mathbf{S}\tilde{\bm{\Psi}}_{k}\rangle+C\|\tilde{\bm{\Psi}}_{0}\|^{2}$
s.t.
$\displaystyle\tilde{\bm{\bm{\Psi}}}_{k+1}-\bar{\mathbf{R}}\tilde{\bm{\bm{\Psi}}}_{k}-\mathbf{P}\tilde{\mathbf{U}}_{k}=\sum\limits_{l=0}^{K-2}\tilde{\mathbf{R}}_{k,k+1-l}\tilde{\bm{\bm{\Psi}}}_{l}.\qquad
k=0,\ldots,K-2$ (48)

In the reference (Shahriari-Mehr and Panahi,, 2022), the dynamics of the
algorithm are depicted as a linear time-invariant (LTI) system, as there are
no communication delays. This distinction increases the complexity of the
analysis, making our study substantially different from their analysis.

As this step is rather substantial and possibly intriguing to readers on its
own, we will apply the result
$\bar{A}_{K}\geq-C\|\tilde{\bm{\bm{\Psi}}}_{0}\|_{\mathrm{F}}^{2}$ for now and
provide a full proof after the convergence rates are determined.

### Step five: providing the rates of convergence

Since $\bar{A}_{K}$ has a lower-bound $-CC_{0}$, where
$C_{0}\coloneqq\|\tilde{\bm{\bm{\Psi}}}_{0}\|_{\mathrm{F}}^{2}$, we obtain the
following rates of convergence.

Consensus: Time-averaged local solutions converge to the consensus solution
$\bar{\mathbf{x}}_{K}$ as

$\|\bar{\mathbf{x}}_{K}-\bar{\mathbf{x}}^{v}_{K}\|_{2}^{2}\leq\frac{C_{0}C}{MK}.\qquad\forall
v\in\mathcal{V}$ (49)

Feasibility gap: The consensus solution $\bar{\mathbf{x}}_{K}$ goes into each
constraint set with the following rate. Hence, the distance between
$\bar{\mathbf{x}}_{K}$ and the feasible set decays at the same rate.

$\mathrm{dist}^{2}(\bar{\mathbf{x}}_{K},S^{v})=\mathcal{O}(\frac{1}{K})\qquad\forall
v\in\mathcal{V}$ (50)

Optimality gap: The objective value converges to the optimal value as

$\left|\sum\limits_{v}f^{v}(\bar{\mathbf{x}}^{v}_{K})-\sum\limits_{v}f^{v}(\mathbf{x}^{*})\right|\leq\frac{C_{0}C}{\mu
K}+\sqrt{\frac{C_{0}CC_{2}}{MK}},$ (51)

where $C_{2}=\sqrt{\sum_{v}\|\mathbf{n}^{v}+\nabla
f^{v}(\mathbf{x}^{*})\|^{2}}$, and it only depends on the choice of the
regular optimal solution (which exists by the assumptions). ∎

Matrix definitions:

We define matrix $\bar{\mathbf{W}}\coloneqq\text{blkdiag}(\mathbf{W})$ as a
block-diagonal matrix where each block on its diagonal is given by a single
row from $\mathbf{W}.$ Similarly, we define matrix
$\bar{\mathbf{Q}}\coloneqq\text{blkdiag}(\mathbf{Q})$.

$\bar{\mathbf{R}}=\left[\begin{array}[]{cccccccc}\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}\\\
\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}\\\
\mathbf{I}&\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}\\\
-\frac{\rho}{\mu}\mathbf{I}&\mathbf{O}&\frac{\rho}{\mu}\mathbf{I}&(1-\alpha)\mathbf{I}&\alpha\mathbf{I}&\mathbf{O}&-\frac{\rho}{\mu}\bar{\mathbf{W}}&\mathbf{O}\\\
\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}&\gamma\mathbf{I}&\mathbf{O}&\mathbf{O}&\bar{\mathbf{Q}}\\\
\mathbf{O}&\mathbf{O}&\mathbf{O}&\eta(\gamma-1)\mathbf{I}&\mathbf{O}&\mathbf{I}&\mathbf{O}&-\eta\bar{\mathbf{Q}}\\\
\mathbf{1}_{M}\otimes\mathbf{I}&\mathbf{O}&\mathbf{O}_{M^{2}\times
M}&\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}\\\
\mathbf{O}&\mathbf{1}_{M}\otimes\mathbf{I}&\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}\end{array}\right]$
(52)
$\left\\{\tilde{\mathbf{R}}_{k,k+1-l}\right\\}_{i,j}=\left\\{\begin{array}[]{cl}-1&\quad
l=k+1,|T^{vu}_{k}|=0,i=(v+5)M+u,j=u\\\ +1&\quad
l=k+1,|T^{vu}_{k}|=0,i=(v+5)M+u,j=(v+5)M+u\\\ -1+\frac{1}{|T^{vu}_{k}|}&\quad
l=k+1,|T^{vu}_{k}|\neq 0,l\in T^{vu}_{k},i=(v+5)M+u,j=u\\\ -1&\quad
l=k+1,|T^{vu}_{k}|\neq 0,l\notin T^{vu}_{k},i=(v+5)M+u,j=u\\\
\frac{1}{|T^{vu}_{k}|}&\quad l\neq k+1,|T^{vu}_{k}|\neq 0,l\in
T^{vu}_{k},i=(v+5)M+u,j=u\\\ \----&\quad------\\\ -1&\quad
l=k+1,|T^{vu}_{k}|=0,i=M^{2}+(v+5)M+u,j=M+u\\\ +1&\quad
l=k+1,|T^{vu}_{k}|=0,i=M^{2}+(v+5)M+u,j=M^{2}+(v+6)M+u\\\
-1+\frac{1}{|T^{vu}_{k}|}&\quad l=k+1,|T^{vu}_{k}|\neq 0,l\in
T^{vu}_{k},i=M^{2}+(v+5)M+u,j=M+u\\\ -1&\quad l=k+1,|T^{vu}_{k}|\neq 0,l\notin
T^{vu}_{k},i=M^{2}+(v+5)M+u,j=M+u\\\ \frac{1}{|T^{vu}_{k}|}&\quad l\neq
k+1,|T^{vu}_{k}|\neq 0,l\in T^{vu}_{k},i=M^{2}+(v+5)M+u,j=M+u\\\
\----&\quad------\\\ 0&\quad\text{otherwise}\end{array}\right.$ (53)
$\mathbf{P}=\left[\begin{array}[]{cc}\mathbf{I}_{M\times
M}&\mathbf{O}_{M\times M}\\\ \mathbf{O}_{M\times M}&\mathbf{I}_{M\times M}\\\
\mathbf{O}_{(2M^{2}+4M)\times M}&\mathbf{O}_{(2M^{2}+4M)\times
M}\end{array}\right]$ (54)
$\mathbf{S}=\left[\begin{array}[]{cccccccc}\left(1-\frac{L\mu}{2}\right)\mathbf{I}-M\zeta\left(\mathbf{I}-\frac{1}{M}\mathbf{1}\mathbf{1}^{T}\right)&\mathbf{O}&-\frac{L\mu-1}{2}\mathbf{I}&-\frac{\mu}{2}\mathbf{I}&\mathbf{O}&\mathbf{O}&\frac{1}{2}\bar{\mathbf{W}}&\mathbf{O}\\\
\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}\\\
-\frac{L\mu-1}{2}\mathbf{I}&\mathbf{O}&-\frac{L\mu}{2}\mathbf{I}&\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}\\\
-\frac{\mu}{2}\mathbf{I}&\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}\\\
\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}\\\
\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}\\\
\frac{1}{2}\bar{\mathbf{W}}&\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}\\\
\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}&\mathbf{O}\end{array}\right]$
(55)

##  11 Solving LQ-PEP

In this proof, $D$ denotes a generic constant whose value might increase each
time it appears.

In order to solve the LQ-PEP and demonstrate that
$-C\|\tilde{\bm{\bm{\Psi}}}_{0}\|_{\mathrm{F}}^{2}\leq\bar{A}_{K}$, we need to
establish that zero is the optimal value of (48). Notably, $\mathbf{S}$ being
indefinite sets LQ-PEP apart from linear quadratic control. If zero is not the
optimum of that optimization, the problem becomes unbounded, pushing the
optimal value to $-\infty$. Consequently, we restrict the optimization problem
in (48) to a unit circle. This constraint guides us towards the following
optimization:

$\displaystyle\min\limits_{\\{\tilde{\bm{\Psi}}_{k}\\}_{k=0}^{K-1},\\{\tilde{\mathbf{U}}_{k}\\}_{k=0}^{K-2}}\sum\limits_{k=0}^{K-1}\frac{1}{2}\langle\tilde{\bm{\Psi}}_{k},\mathbf{S}\tilde{\bm{\Psi}}_{k}\rangle+\frac{C}{2}\|\tilde{\bm{\Psi}}_{0}\|^{2}$
s.t.
$\displaystyle\tilde{\bm{\bm{\Psi}}}_{k+1}-\bar{\mathbf{R}}\tilde{\bm{\bm{\Psi}}}_{k}-\mathbf{P}\tilde{\mathbf{U}}_{k}=\sum\limits_{l=0}^{K-2}\tilde{\mathbf{R}}_{k,k+1-l}\tilde{\bm{\bm{\Psi}}}_{l},\qquad
k=0,\ldots,K-2$
$\displaystyle\frac{1}{2}\|\tilde{\bm{\bm{\Psi}}}_{0}\|_{\mathrm{F}}^{2}+\frac{1}{2}\sum\limits_{k=0}^{K-2}\|\tilde{\mathbf{U}}_{k}\|_{\mathrm{F}}^{2}\leq\frac{1}{2}.\qquad\qquad$
(56)

If the optimal value for this new optimization is negative, the original
optimization in (48) becomes unbounded. Thus, asserting that zero is the
optimal solution for the original problem is equivalent to stating that the
optimal value of the optimization in (56) is non-negative.

To show that the optimal value of (56) is non-negative, we introduce the dual
Lagrangian multipliers $\bm{\Lambda}_{k}$ and $\beta\geq 0$. Consequently, we
construct the dual optimization problem. Upon solving this problem with
respect to both the primal and dual variables, we achieve

$\displaystyle\tilde{\bm{\bm{\Psi}}}_{k+1}-\bar{\mathbf{R}}\tilde{\bm{\bm{\Psi}}}_{k}-\mathbf{P}\tilde{\mathbf{U}}_{k}$
$\displaystyle=\sum\limits_{l=0}^{K-2}\tilde{\mathbf{R}}_{k,k+1-l}\tilde{\bm{\bm{\Psi}}}_{l}$
$\displaystyle\qquad k=0,\ldots,K-2$ (57)
$\displaystyle\bm{\Lambda}_{k-1}-\bar{\mathbf{R}}^{T}\bm{\Lambda}_{k}+\mathbf{S}\tilde{\bm{\bm{\Psi}}}_{k}$
$\displaystyle=\sum\limits_{l=0}^{K-2}\tilde{\mathbf{R}}^{T}_{l,l-k+1}\bm{\Lambda}_{l}$
$\displaystyle\qquad k=0,\ldots,K-2$ (58)
$\displaystyle\beta\tilde{\mathbf{U}}_{k}-\mathbf{P}^{T}\bm{\Lambda}_{k}$
$\displaystyle=\mathbf{O}$ $\displaystyle\qquad k=0,\ldots,K-2$ (59)

with boundary conditions

$\displaystyle\bm{\Lambda}_{-1}+(\beta+C)\tilde{\bm{\bm{\Psi}}}_{0}$
$\displaystyle=\mathbf{O}$ (60)
$\displaystyle\bm{\Lambda}_{K-2}+\mathbf{S}\tilde{\bm{\bm{\Psi}}}_{K-1}$
$\displaystyle=\mathbf{O}$ (61)

After some simplifications and substitutions, we notice that the optimal value
of this problem
$-\beta\left(\|\bm{\Psi}_{0}\|_{\mathrm{F}}^{2}+\sum_{k=0}^{K-2}\|\tilde{\mathbf{U}}_{k}\|_{\mathrm{F}}^{2}\right)$.
Since $\beta$ is positive, the optimal solution is negative unless the system
of linear recurrences in (57), (58), and (59) has no non-zero solution.

In summary, if we show the system of linear recurrences in (57), (58), and
(59) has no non-zero solution, zero will be the optimal value of the
optimization problems in (48) and (56), consequently, $\bar{A}_{K}$ has a
lower bound $-C\|\tilde{\bm{\bm{\Psi}}}_{0}\|_{\mathrm{F}}^{2}.$

To show this, first, we remove the dependence on $\tilde{\mathbf{U}}_{k}$.
Then, we have the new system of linear recurrences

$\displaystyle\tilde{\bm{\bm{\Psi}}}_{k+1}-\bar{\mathbf{R}}\tilde{\bm{\bm{\Psi}}}_{k}-\frac{1}{\beta}\mathbf{P}\mathbf{P}^{T}\bm{\Lambda}_{k}$
$\displaystyle=\sum\limits_{l=0}^{K-2}\tilde{\mathbf{R}}_{k,k+1-l}\tilde{\bm{\bm{\Psi}}}_{l}$
$\displaystyle\qquad k=0,\ldots,K-2$ (62)
$\displaystyle\bm{\Lambda}_{k-1}-\bar{\mathbf{R}}^{T}\bm{\Lambda}_{k}+\mathbf{S}\tilde{\bm{\bm{\Psi}}}_{k}$
$\displaystyle=\sum\limits_{l=0}^{K-2}\tilde{\mathbf{R}}^{T}_{l,l-k+1}\bm{\Lambda}_{l}$
$\displaystyle\qquad k=0,\ldots,K-2$ (63)

with the same boundary conditions as in (60) and (61). By defining the finite
duration $z-$transforms of $\tilde{\bm{\bm{\Psi}}}_{k}$, and
$\bm{\Lambda}_{k}$ as

$\tilde{\bm{\bm{\Psi}}}(z)=\sum\limits_{k=0}^{K-1}\tilde{\bm{\bm{\Psi}}}_{k}z^{k},\qquad\bm{\Lambda}(z)=\sum\limits_{k=0}^{K-1}\bm{\Lambda}_{k-1}z^{k},$
(64)

and taking the $z-$transform of (62) and (63), we have

$\left[\begin{array}[]{c}\tilde{\bm{\bm{\Psi}}}(z)\\\
\bm{\Lambda}(z)\end{array}\right]=\mathbf{F}_{\beta}^{-1}(z)\left\\{\left[\begin{array}[]{c}\mathbf{S}\tilde{\bm{\bm{\Psi}}}_{K-1}z^{K-1}+\bm{\Lambda}_{K-2}z^{K-1}-z^{-1}\bar{\mathbf{R}}^{T}\bm{\Lambda}_{-1}\\\
\tilde{\bm{\bm{\Psi}}}_{0}-\bar{\mathbf{R}}\tilde{\bm{\bm{\Psi}}}_{K-1}z^{K}-\frac{1}{\beta}\mathbf{P}\mathbf{P}^{T}\bm{\Lambda}_{-1}\end{array}\right]+\left[\begin{array}[]{c}\sum\limits_{k,l=0}^{K-2}\tilde{\mathbf{R}}^{T}_{l,l-k+1}\bm{\Lambda}_{l}z^{k}\\\
\sum\limits_{k,l=0}^{K-2}\tilde{\mathbf{R}}_{k,k+1-l}\tilde{\bm{\bm{\Psi}}}_{l}z^{k+1}\end{array}\right]\right\\}$
(65)

where

$\mathbf{F}_{\beta}(z)=\left[\begin{array}[]{cc}\mathbf{S}&\mathbf{I}-z^{-1}\bar{\mathbf{R}}^{T}\\\
\mathbf{I}-z\bar{\mathbf{R}}&-\frac{1}{\beta}\mathbf{P}\mathbf{P}^{T}\end{array}\right].$
(66)

From the definition of inverse $z-$transform, we have

$\displaystyle\tilde{\bm{\bm{\Psi}}}_{m}$ $\displaystyle=\frac{1}{2\pi
j}\;\oint z^{-(m+1)}\tilde{\bm{\bm{\Psi}}}(z)\mathrm{d}z,$ (67)
$\displaystyle\bm{\Lambda}_{m-1}$ $\displaystyle=\frac{1}{2\pi j}\;\oint
z^{-(m+1)}\bm{\Lambda}(z)\mathrm{d}z.$ (68)

Given the structure of $\mathbf{F}_{\beta}(z)^{-1}$ from Lemmas 1 and 2, we
compute the inverse z-transforms, based on Lemma 3, as:

$\left[\begin{array}[]{c}\tilde{\bm{\bm{\Psi}}}_{m}\\\
\bm{\Lambda}_{m-1}\end{array}\right]\;=\;\sum\limits_{p=1}^{q}\mathbf{a}_{p}\mathbf{b}_{p}^{T}\left[\begin{array}[]{c}\big{(}\mathbf{S}\tilde{\bm{\bm{\Psi}}}_{K-1}+\bm{\Lambda}_{K-2}\big{)}z_{p}^{K-m-2}+\bar{\mathbf{R}}^{T}\bm{\Lambda}_{-1}z_{p}^{m+2}\\\
\big{(}\frac{1}{\beta}\mathbf{P}\mathbf{P}^{T}\bm{\Lambda}_{-1}-\tilde{\bm{\bm{\Psi}}}_{0}\big{)}z_{p}^{m+1}-\bar{\mathbf{R}}\tilde{\bm{\bm{\Psi}}}_{K-1}z_{p}^{K-m-1}\end{array}\right]\;+\;\left[\begin{array}[]{c}\bar{\bm{\Psi}}_{m}\\\
\bar{\bm{\Lambda}}_{m-1}\end{array}\right],$ (69)

with

$\left[\begin{array}[]{c}\bar{\bm{\Psi}}_{m}\\\
\bar{\bm{\Lambda}}_{m-1}\end{array}\right]\;=\;\sum\limits_{p=1}^{q}\mathbf{a}_{p}\mathbf{b}_{p}^{T}\left[\begin{array}[]{c}\sum\limits_{l=0}^{K-2}\left(\sum\limits_{k=m+1}^{K-2}\tilde{\mathbf{R}}^{T}_{l,l-k+1}z_{p}^{|k-m-1|}-\sum\limits_{k=0}^{m}\tilde{\mathbf{R}}^{T}_{l,l-k+1}z_{p}^{|k-m-1|}\right)\bm{\Lambda}_{l}\\\
\sum\limits_{l=0}^{K-2}\left(\sum\limits_{k=m}^{K-2}\tilde{\mathbf{R}}_{k,k+1-l}z_{p}^{|k-m|}-\sum\limits_{k=0}^{m-1}\tilde{\mathbf{R}}_{k,k+1-l}z_{p}^{|k-m|}\right)\tilde{\bm{\bm{\Psi}}}_{l}\end{array}\right],$
(70)

where $\mathbf{a}_{p}\in\ker(\mathbf{F}_{\beta}(z_{p}))$ and is assumed to be
in the form of
$\mathbf{a}_{p}=\left[\mathbf{n}_{p,\tilde{\bm{\bm{\Psi}}}}^{T}\;\;\mathbf{n}_{p,\bm{\Lambda}}^{T}\right]^{T}$.
The set $\\{z_{1},\ldots,z_{q}\\}$ contains the poles of
$\mathbf{F}_{\beta}(z)^{-1}$, which are inside the unit circle. We interpret
the result in (69) as the solution to the linear system of recurrence
equations represented in (62) and (63). This solution consists of two parts.
The first part includes the solution to the homogeneous system, which is
achieved by setting the right-hand side of equations (62) and (63) to zero.
The second part, given in (70), provides a particular solution, which is
denoted by $\bar{\bm{\Psi}}_{m},\bar{\bm{\Lambda}}_{m-1}$ to simplify our
further analysis.

For the sake analysis, we also write the homogeneous solution as

$\sum\limits_{p=1}^{q}\mathbf{a}_{p}\bm{\gamma}_{p}^{T}z_{p}^{-m}\;,$

where

$\bm{\gamma}_{p}^{T}=\mathbf{b}_{p}^{T}\left[\begin{array}[]{c}\big{(}\mathbf{S}\tilde{\bm{\bm{\Psi}}}_{K-1}+\bm{\Lambda}_{K-2}\big{)}z_{p}^{K-2}+\bar{\mathbf{R}}^{T}\bm{\Lambda}_{-1}z_{p}^{2m+2}\\\
\big{(}\frac{1}{\beta}\mathbf{P}\mathbf{P}^{T}\bm{\Lambda}_{-1}-\tilde{\bm{\bm{\Psi}}}_{0}\big{)}z_{p}^{2m+1}-\bar{\mathbf{R}}\tilde{\bm{\bm{\Psi}}}_{K-1}z_{p}^{K-1}\end{array}\right].$
(71)

This general solution can also be written as

$\sum\limits_{p=1}^{q}\mathbf{a}_{p}\bm{\gamma}_{p}^{T}z_{p}^{-m}=\left[\mathbf{a}_{1}z_{1}^{-m}\;\;\ldots\;\;\mathbf{a}_{q}z_{q}^{-m}\right]\left[\begin{array}[]{c}\bm{\gamma}_{1}^{T}\\\
\vdots\\\ \bm{\gamma}_{q}^{T}\end{array}\right].$ (72)

The solution to the linear recurrence system must meet its boundary
conditions. To achieve this, we substitute the solution into equations (60)
and (61). This substitution yields

$\underbrace{\left[\begin{array}[]{cccc}\mathbf{n}_{1,\bm{\Lambda}}+(\beta+C)\mathbf{n}_{1,\tilde{\bm{\bm{\Psi}}}}&\ldots&\mathbf{n}_{q,\bm{\Lambda}}+(\beta+C)\mathbf{n}_{q,\tilde{\bm{\bm{\Psi}}}}\\\
(\mathbf{n}_{1,\bm{\Lambda}}+\mathbf{S}\mathbf{n}_{1,\tilde{\bm{\bm{\Psi}}}})z_{1}^{-(K-1)}&\ldots&(\mathbf{n}_{q,\bm{\Lambda}}+\mathbf{S}\mathbf{n}_{q,\tilde{\bm{\bm{\Psi}}}})z_{q}^{-(K-1)}\end{array}\right]}_{\mathbf{M}}\left[\begin{array}[]{c}\bm{\gamma}_{1}^{T}\\\
\vdots\\\
\bm{\gamma}_{q}^{T}\end{array}\right]=-\left[\begin{array}[]{c}\bar{\bm{\Lambda}}_{-1}+(\beta+C)\bar{\bm{\Psi}}_{0}\\\
\bar{\bm{\Lambda}}_{K-2}+\mathbf{S}\bar{\bm{\Psi}}_{K-1}\end{array}\right].$
(73)

To simplify, we define $\mathbf{M}$ as illustrated in (73). Hence, for the
solution to meet the boundary conditions, the following relation should be
satisfied

$\left[\begin{array}[]{c}\bm{\gamma}_{1}^{T}\\\ \vdots\\\
\bm{\gamma}_{q}^{T}\end{array}\right]=-\mathbf{M}^{-1}\left[\begin{array}[]{c}\bar{\bm{\Lambda}}_{-1}+(\beta+C)\bar{\bm{\Psi}}_{0}\\\
\bar{\bm{\Lambda}}_{K-2}+\mathbf{S}\bar{\bm{\Psi}}_{K-1}\end{array}\right].$
(74)

Moving forward, the last step is to impose bounds on the aforementioned
equalities to illustrate that these equations cannot be concurrently true
except when zero is the optimal solution. To accomplish this, we make the
equations simpler by looking at them column by column. Without loss of
generality, If we consider any arbitrary column of $\tilde{\bm{\bm{\Psi}}}$
and $\bm{\Lambda}$, , and demonstrate that this column equals zero, then all
the columns would be zero, completing our proof. For this reason, by
$\|\tilde{\bm{\bm{\Psi}}}\|$, we are referring to a proper vector norm
associated to one random column of $\tilde{\bm{\bm{\Psi}}}$. Using this
terminology and considering the $\ell^{\infty}$ vector norm, we define

$\displaystyle\lambda$
$\displaystyle=\;\max_{m=0,\ldots,K-1}\;\max\left\\{\|\tilde{\bm{\bm{\Psi}}}_{m}\|_{\infty},\|\bm{\Lambda}_{m-1}\|_{\infty}\right\\},$
(75) $\displaystyle\bar{\lambda}$
$\displaystyle=\;\max_{m=0,\ldots,K-1}\;\max\Big{\\{}\|\bar{\bm{\Psi}}_{m}\|_{\infty},\|\bar{\bm{\Lambda}}_{m-1}\|_{\infty}\Big{\\}}.$
(76)

Putting all equations in (69), (70), (72), and (74) together, we get

$\left[\begin{array}[]{c}\tilde{\bm{\bm{\Psi}}}_{m}\\\
\bm{\Lambda}_{m-1}\end{array}\right]\;=\;-\left[\mathbf{a}_{1}z_{1}^{-m}\;\;\ldots\;\;\mathbf{a}_{q}z_{q}^{-m}\right]\mathbf{M}^{-1}\left[\begin{array}[]{c}\bar{\bm{\Lambda}}_{-1}+(\beta+C)\bar{\bm{\Psi}}_{0}\\\
\bar{\bm{\Lambda}}_{K-2}+\mathbf{S}\bar{\bm{\Psi}}_{K-1}\end{array}\right]\;+\;\left[\begin{array}[]{c}\bar{\bm{\Psi}}_{m}\\\
\bar{\bm{\Lambda}}_{m-1}\end{array}\right].$ (77)

Taking norm, considering (75) and (76), and using the triangle inequality, we
get

$\lambda\leq\left\|\left[\mathbf{a}_{1}z_{1}^{-m}\;\;\ldots\;\;\mathbf{a}_{q}z_{q}^{-m}\right]\mathbf{M}^{-1}\left[\begin{array}[]{c}\bar{\bm{\Lambda}}_{-1}+(\beta+C)\bar{\bm{\Psi}}_{0}\\\
\bar{\bm{\Lambda}}_{K-2}+\mathbf{S}\bar{\bm{\Psi}}_{K-1}\end{array}\right]\right\|\;+\;\bar{\lambda}.$
(78)

From (73) and (74), one can show that for sufficiently large $K$, there exists
a generic constant $D$, dependant on constant $C$, such that
$\|\bm{\gamma}_{p}^{T}\|_{\infty}\leq D|z_{p}|^{K}\bar{\lambda}$, for all
$p=1,\ldots,q$. Then, we have

$\lambda\;\leq\;\bar{\lambda}+D\bar{\lambda}\,\max_{m}\,\sum\limits_{p=1}^{q}|z_{p}|^{K-m}\|\mathbf{a}_{p}\|\;\leq\;D\bar{\lambda}$
(79)

The above inequalities result in $\lambda\leq D\bar{\lambda}$. Now, we aim to
determine an upper bound for $\bar{\lambda}$. Recall the definition in (70)

$\left[\begin{array}[]{c}\bar{\bm{\Psi}}_{m}\\\
\bar{\bm{\Lambda}}_{m-1}\end{array}\right]\;=\;\sum\limits_{p=1}^{q}\mathbf{a}_{p}\mathbf{b}_{p}^{T}\left[\begin{array}[]{c}\sum\limits_{l=0}^{K-2}\left(\sum\limits_{k=m+1}^{K-2}\tilde{\mathbf{R}}^{T}_{l,l-k+1}z_{p}^{|k-m-1|}-\sum\limits_{k=0}^{m}\tilde{\mathbf{R}}^{T}_{l,l-k+1}z_{p}^{|k-m-1|}\right)\bm{\Lambda}_{l}\\\
\sum\limits_{l=0}^{K-2}\left(\sum\limits_{k=m}^{K-2}\tilde{\mathbf{R}}_{k,k+1-l}z_{p}^{|k-m|}-\sum\limits_{k=0}^{m-1}\tilde{\mathbf{R}}_{k,k+1-l}z_{p}^{|k-m|}\right)\tilde{\bm{\bm{\Psi}}}_{l}\end{array}\right],$
(80)

Taking norm, considering (75) and (76), we have

$\displaystyle\bar{\lambda}$
$\displaystyle\leq\max_{m}\sum\limits_{p}\|\mathbf{a}_{p}\|_{\infty}\left|\mathbf{b}_{p}^{T}\left[\begin{array}[]{c}\sum\limits_{l=0}^{K-2}\left(\sum\limits_{k=m+1}^{K-2}\tilde{\mathbf{R}}^{T}_{l,l-k+1}z_{p}^{|k-m-1|}-\sum\limits_{k=0}^{m}\tilde{\mathbf{R}}^{T}_{l,l-k+1}z_{p}^{|k-m-1|}\right)\bm{\Lambda}_{l}\\\
\sum\limits_{l=0}^{K-2}\left(\sum\limits_{k=m}^{K-2}\tilde{\mathbf{R}}_{k,k+1-l}z_{p}^{|k-m|}-\sum\limits_{k=0}^{m-1}\tilde{\mathbf{R}}_{k,k+1-l}z_{p}^{|k-m|}\right)\tilde{\bm{\bm{\Psi}}}_{l}\end{array}\right]\right|$
(83) $\displaystyle\leq
D\lambda\;\max_{m}\;\sum\limits_{p}\;\max_{l}\;\;\max\\{\|\mathbf{A}_{l,p,m}\|,\|\mathbf{B}_{l,p,m}\|\\}$
(84)

where

$\displaystyle\mathbf{A}_{k,p,m}$
$\displaystyle=\sum_{l=m+1}^{K-2}\tilde{\mathbf{R}}^{T}_{k,k+1-l}z_{p}^{|l-m-1|}-\sum\limits_{l=0}^{m}\tilde{\mathbf{R}}^{T}_{k,k+1-l}z_{p}^{|l-m-1|}$
(85) $\displaystyle\mathbf{B}_{l,p,m}$
$\displaystyle=b\sum_{k+m}^{K-2}\tilde{\mathbf{R}}_{k,k+1-l}z_{p}^{k-m}-\sum\limits_{k=0}^{m-1}\tilde{\mathbf{R}}_{k,k+1-l}z_{p}^{k-m}$
(86)

where $\|A\|_{1,\infty}$ is the $\ell_{1}-\ell_{\infty}$ mixed norm of matrix
$A$, given by $||A||_{1,\infty}=\max_{i}\sum_{j}|a_{ij}|$. To derive (84), we
applied the triangle inequality twice. For (83), we utilized the relationship
$|\mathbf{x}^{T}\mathbf{y}|\leq\|\mathbf{x}\|_{1}\|\mathbf{y}\|_{\infty}.$
Furthermore, to equate the powers of the elements, we incorporated $|z|^{-1}$
into the constant $D$, keeping in mind that $D$ is a generic constant and its
value can be adjusted as needed.

Now, by considering the definition of $\tilde{\mathbf{R}}_{k,k+1-l}$ in (53),
or specifically its non-zero elements, we define the $p^{\text{th}}$ delay
power spectrum of node $(v,u)$ and its delay response at time $m$,
respectively, as

$\tau^{vu}_{p,k}(\beta)=\frac{1}{|T_{k}^{vu}|}\sum\limits_{k^{\prime}\in
T_{k}^{vu}}z_{p}^{|k+1-k^{\prime}|}(\beta)-1,$ (87)
$\kappa^{vu}_{m}(\beta)=\sum\limits_{p}\sum\limits_{k\mid
T_{k}^{vu}\neq\emptyset}|z_{p}(\beta)|^{|k-m|}\tau^{vu}_{p,k}(\beta).$ (88)

Then, the delay response $\kappa$ is defined as the supremum of
$\kappa^{vu}_{m}(\beta)$ over every time $m$ and link $(v,u)$ and $\beta$.
Therefore, the inequality represented in (84) simplifies to

$\bar{\lambda}\leq D\kappa\lambda.$ (89)

From (79), we also have $\lambda\leq D\bar{\lambda}$. Combining these two
inequalities leads us to

$\lambda\leq D\kappa\lambda.$ (90)

Keep in mind that $D$ is a generic constant, and it relies on $C(\kappa)$.
Specifically, when $C$ increases, $D$ decreases, but it will not go lower than
$1/\kappa_{0}$. If we have a delay response, $\kappa$, from one delay event,
and it satisfies $\kappa<\kappa_{0}$, then there are possible values of $C$
where $D\kappa$ is less than 1. Based on (90), this makes $\lambda=0$, which
leads to $\tilde{\bm{\bm{\Psi}}}_{m}=\mathbf{O}$ and
$\bm{\Lambda}_{m-1}=\mathbf{O}$. This gives us the full proof.

## Lemmas

######  Lemma 1 (Inverse-matrix decomposition).

let $\mathbf{F}(z)$ be a square matrix of size $m$, and the set
$\\{z_{1},z_{2},\ldots,z_{P}\\}$ contains all the simple roots of
$\det(\mathbf{F}(z))=0$. Further, assume $z\mathbf{F}_{\beta}^{-1}(z)$ has a
finite limit as $z\to 0$.Then, the inverse matrix $\mathbf{F}^{-1}(z)$ can be
decomposed into

$\mathbf{F}^{-1}(z)=\sum\limits_{p=1}^{P}\frac{\mathbf{a}_{p}\mathbf{b}_{p}^{T}}{z-z_{p}},$
(91)

where

$\mathbf{a}_{p}\in\ker(\mathbf{F}(z_{p}))$
$\mathbf{b}_{p}\perp\text{image}(\mathbf{F}(z_{p})).$

###### Proof.

From Cramer’s rule for matrix inversion, considering that
$\det\mathbf{F}(z)=0$ has only simple roots, we decompose $\mathbf{F}^{-1}(z)$
into:

$\mathbf{F}^{-1}(z)=\sum\limits_{p=1}^{P}\frac{\mathbf{L}_{p}}{z-z_{p}},$ (92)

where $\mathbf{L}_{p}$ are some Known matrices. Since $\mathbf{F}(z)$ is
invertible for every element $\mathbf{b}^{\prime}$ in the range space (image)
of $\mathbf{F}(z)$, for all $z$, the following condition should be hold

$\mathbf{L}_{p}\mathbf{b}^{\prime}=\mathbf{0},$ (93)

otherwise the inverse does not exist at $z=z_{p}$.

From the fundamental theorem of linear algebra, specifically the Rank-Nullity
theorem, the rank or the dimension of the range space of $\mathbf{F}(z_{p})$
is $m-1$ since the nullity, i.e. the dimension of the null space, of
$\mathbf{F}(z_{p})$ is one. Therefore, the range space can be represented by
only one vector orthogonal to this subspace, represented by $\mathbf{b}_{p}$.

Since any vector in the range space of $\mathbf{F}(z)$, containing
$\mathbf{b}^{\prime}$, is orthogonal to $\mathbf{b}_{p}$, $\mathbf{L}_{p}$
should be in the form of $\mathbf{L}_{p}=\mathbf{a}_{p}\mathbf{b}_{p}^{T}$ to
satisfy (93). By following the same procedure for $\mathbf{F}(z)^{T}$,
considering $\det(\mathbf{F}(z)^{T})=\det(\mathbf{F}(z)),$ we observe that
$\mathbf{a}_{p}$ is an element in the null space of $\mathbf{F}(z_{p})$. ∎

######  Lemma 2 (Inversely-related roots).

let matrix $\mathbf{F}(z)$ be in the following form

$\mathbf{F}(z)=\left[\begin{array}[]{cc}\mathbf{A}&z^{-1}\mathbf{B}^{T}\\\
z\mathbf{B}&\mathbf{C}\end{array}\right],$ (94)

with symmetric matrices $\mathbf{A},\mathbf{B}$. Then, the roots of
$\det(\mathbf{F}(z))$ are _inversely-related_ , meaning that if $z_{i}$ is a
root, then its inverse $z_{i}^{-1}$ is also a root.

###### Proof.

Given the configuration of $\mathbf{F}(z)$, we have
$\mathbf{F}(z)^{T}=\mathbf{F}(z^{-1})$. Because
$\det\mathbf{F}(z)^{T}=\det\mathbf{F}(z^{-1})=\det\mathbf{F}(z)$, it follows
that if $z_{i}$ serves as a root, its reciprocal $z_{i}^{-1}$ is also a root.
∎

######  Lemma 3.

[Contour integral over the unit circle] Let $f(z)=\frac{z^{k}}{z-z_{p}}$. The
contour integral of $f(z)$ over the unit circle $\Gamma$ is given by

$\frac{1}{2\pi
j}\oint_{\Gamma}\,\frac{z^{k}}{z-z_{p}}=\begin{cases}z_{p}^{k}&|z_{p}|<1,k\geq
0\\\ -z_{p}^{k}&|z_{p}|>1,k<0\\\ 0&\text{otherwise}.\end{cases}$ (95)

###### Proof.

To compute the integral, we consider the following four cases:

  * •

$|z_{p}|>1,k\geq 0$: The function is analytic inside the contour. Based on the
Cauchy integral theorem, the integral equals $0.$

  * •

$|z_{p}|<1,k\geq 0$: $z=z_{p}$ is the only singular point of $f(z)$ inside the
contour. Based on the residue theorem, the integral equals $z_{p}^{k}$.

  * •

$|z_{p}|>1,k<0$: $z=0$ is the only singularity of order $k$ of $f(z)$. The
residue at $z=0$ is computed as

$\text{Res}[f(z),0]=\lim_{z\to
0}\frac{1}{(k-1)!}\frac{d^{k-1}}{dz^{k-1}}\left[\frac{1}{z-z_{p}}\right]=-z_{p}^{k}$

Based on the Residue Theorem, the solution of the integral is $-z_{p}^{k}$.

  * •

$|z_{p}|<1,k<0$: The function $f(z)$ has two singularities at $z=0$ and
$z=z_{p}$. Based on the Residue Theorem, considering that the residues cancel
each other out, the integral equals 0.

Combining the above cases, we obtain the solution to the contour integral as
given in (95). ∎