# Self-Healing First-Order Distributed Optimization with Packet Loss

Israel L. Donato Ridgley1,4 Member, IEEE
Randy A. Freeman1,3,4 Member, IEEE  and Kevin M. Lynch2,3,4 Fellow, IEEE All
authors are affiliated with Northwestern University, Evanston, IL 60208 USA
(e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>kmlynch@northwestern.edu).1Department of Electrical & Computer Engineering;
2Department of Mechanical Engineering; 3Northwestern Institute on Complex
Systems; 4Center for Robotics and BiosystemsThis material is based upon work
supported by the National Science Foundation under Grant No. CMMI-2024774.

###### Abstract

We describe SH-SVL, a parameterized family of first-order distributed
optimization algorithms that enable a network of agents to collaboratively
calculate a decision variable that minimizes the sum of cost functions at each
agent. These algorithms are _self-healing_ in that their convergence to the
correct optimizer can be guaranteed even if they are initialized randomly,
agents join or leave the network, or local cost functions change. We also
present simulation evidence that our algorithms are self-healing in the case
of dropped communication packets. Our algorithms are the first single-
Laplacian methods for distributed convex optimization to exhibit all of these
characteristics. We achieve self-healing by sacrificing internal stability, a
fundamental trade-off for single-Laplacian methods.

##  1 Introduction

The distributed optimization problem involves a network of $n$ agents that
each calculates a decision vector that minimizes a global additive objective
function of the form $f(\cdot)=\sum_{i}f_{i}(\cdot)$, where $f_{i}$ denotes
the local convex objective function known only to agent $i$. Specifically,
each agent maintains a local estimate $x_{i}$ of the global minimizer

$x_{\text{opt}}=\operatorname*{arg\,min}_{\theta}\sum_{i}f_{i}(\theta),$ (1)

which we assume is unique. The agents reach consensus $x_{i}=x_{\text{opt}}$
by computing the gradients of their local objective functions
$\nabla\mkern-1.5muf_{i}(x_{i})$ and passing messages along the links of the
communication network.

Distributed optimization problems of this form have broad application. For
example, a distributed set of servers or sensors could perform a learning task
(e.g., classification) using their local data without uploading the data to a
central server for bandwidth, resiliency, or privacy reasons [1]. Swarms of
robots can use distributed optimization to plan motions to solve the
rendezvous problem [2].

The optimization of a collective cost function in a network setting has seen
considerable interest over the last decade [3, 4, 5, 6, 7, 8, 9, 10].
Recently, several authors have adapted methods from control theory to study
distributed optimization algorithms as linear systems in feedback with
uncertainties constrained by integral quadratic constraints (IQCs) [11, 12,
3]. These works have made it possible to more easily compare the various known
algorithms across general classes of cost functions and graph topologies.
Furthermore, the framework has been extended to stochastic optimization
algorithms by using techniques from Markov Jump Linear Systems in the
centralized setting [13] as well as the setting of synchronous packet loss in
autonomous vehicle fleets [14].

The work [3] uses these techniques to describe several recent distributed
gradient-tracking optimization algorithms within a common framework, then
describes a new algorithm (SVL) within that framework that achieves a superior
worst-case convergence rate. However, all of the algorithms considered in [3]
share a common undesirable trait: to reach the correct solution, their states
must start in a particular subspace of the overall global state space and
remain in it at every time step. If for any reason the state trajectories
leave this subspace (e.g., incorrect initialization, dropped packets,
computation errors, agents leaving the network, changes to objective functions
due to continuous data collection), then the system will no longer converge to
the minimizer. Such methods cannot automatically recover from disturbances or
other faults that displace their trajectories from this subspace; in other
words, they are not self-healing. The method presented in [15] is self-
healing; however, the decaying gradient step size limits the rate of
convergence and causes disturbances to generate increasingly long transients.

In this paper, we extend our results from dynamic average consensus estimators
[16, 17] to design SH-SVL, a family of self-healing distributed optimization
algorithms whose trajectories need not evolve on a pre-defined subspace. In
practice, this means that our algorithms can be arbitrarily initialized,
agents can join or leave the network at will, and agents can change their
objective functions as necessary, such as when they collect new data. An
important consequence of the self-healing property is that our algorithm can
be modified with a low-overhead packet-loss protocol which allows the
algorithm to recover from lost or corrupted packets.

We refer to distributed optimization algorithms that communicate one or two
variables (having the same vector dimension as the decision variable $x_{i}$)
per time step as single- and double-Laplacian methods, respectively. Examples
of single-Laplacian methods are SVL and NIDS, while examples of double-
Laplacian methods are uEXTRA and DIGing [3, 5, 6, 7, 8]. SH-SVL is the first
self-healing single-Laplacian methods for convex optimization that converge to
the exact (rather than an approximate) solution with a fixed step size (see
[16, 18] for the specific case of average consensus). They achieve self-
healing by sacrificing internal stability, a fundamental trade-off for single-
Laplacian gradient tracking-methods. In particular, each agent will have an
internal state that grows linearly in time in steady state, but because such
growth is not exponential it will not cause any numerical issues unless the
algorithm runs over a long time horizon. Double-Laplacian methods can achieve
both internal stability and self-healing, but they require twice as much
communication per time step and converge no faster than single-Laplacian
methods [3, 17].

This paper uses the same algorithms as [19], but provides theoretical and
numerical proof for convergence of SH-SVL (Algorithm 2) with packet loss. The
performance of SH-SVL is explored in the setting of synchronous packet loss on
arbitrary graphs and edgewise packet loss in small directed cycles. We perform
a numerical study to investigate the performance when algorithm parameters are
optimized for the specific level of packet loss versus no packet loss. In
addition, we show that SH-SVL can recover the same performance as SVL when
parameters are optimized.

##  2 Preliminaries and Main Results

###  2.1 Notation and terminology

Let $\mathbb{1}_{n}$ be the $n$-dimensional column vector of all ones, $I_{n}$
be the identity matrix in $\mathbb{R}^{n\times n}$, and
$\Pi_{n}=\frac{1}{n}\mathbb{1}\mathbb{1}^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}$
be the projection matrix onto the vector $\mathbb{1}_{n}$. We drop the
subscript $n$ when the size is clear from context. We refer to the one-
dimensional linear subspace of $\mathbb{R}^{n}$ spanned by the vector
$\mathbb{1}_{n}$ as the _consensus direction_ or the _consensus subspace_. We
refer to the $(n-1)$-dimensional subspace of $\mathbb{R}^{n}$ associated with
the projection matrix $(I_{n}-\Pi_{n})$ as the _disagreement direction_ or
subspace. The $i$-th column of the identity matrix is denoted $e_{i}$, with a
size determined by context. The $i$-th row of the matrix $A$ is written
$A_{r_{i}}=e_{i}^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}A$.
Given $A\in\mathbb{R}^{n\times n}$, we let $D(A)\in\mathbb{R}^{n\times n}$
denote the diagonal matrix whose diagonal elements are the same as A.

Subscripts denote the agent index whereas superscripts denote the time index.
The symbol $\otimes$ represents the Kronecker product. $A^{+}$ indicates the
Moore-Penrose inverse of $A$. Symmetric quadratic forms
$x^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}\mkern-1.5muAx$
are written as
$[\star]^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}\mkern-1.5muAx$
to save space when $x$ is long. The local decision variables are
$d$-dimensional and represented as a row vector, i.e.,
$x_{i}\in\mathbb{R}^{1\times d}$, and the local gradients are a map
$\nabla\mkern-1.5muf_{i}:\mathbb{R}^{1\times d}\rightarrow\mathbb{R}^{1\times
d}$. The symbol $||\cdot||$ refers to the Euclidean norm of vectors and the
spectral norm of matrices. The condition number of the symmetric matrix $A$ is
defined by
$\operatorname{cond}(A)=\lambda_{\textnormal{max}}(A)/\lambda_{\textnormal{min}}(A)$.
The expected value operator is $\mathbb{E}[\cdot]$.

We model a network of $n$ agents participating in a distributed computation as
a weighted digraph $\mathcal{G}=(\mathcal{V},\mathcal{E})$, where
$\mathcal{V}=\\{1,...,n\\}$ is the set of $n$ nodes (or vertices) and
$\mathcal{E}$ is the set of edges such that if $(i,j)\in\mathcal{E}$ then node
$i$ can receive information from $j$. We define an ordering on $\mathcal{E}$
such that the set is in ascending order according to the receiving node index
followed by the sending node index. We make use of the weighted graph
Laplacian $\mathcal{L}\in\mathbb{R}^{n\times n}$ associated with $\mathcal{G}$
such that $-\mathcal{L}_{ij}$ is the weight on edge $(i,j)\in\mathcal{E}$,
$\mathcal{L}_{ij}=0$ when $(i,j)\not\in\mathcal{E}$ and $i\neq j$, and the
diagonal elements of $\mathcal{L}$ are $\mathcal{L}_{ii}=-\sum_{j\neq
i}\mathcal{L}_{ij}$, so that $\mathcal{L}\mathbb{1}=0$. We define
$\sigma=||I-\Pi-\mathcal{L}||$, which is a parameter related to the edge
weights and the graph connectivity.

Throughout this work we stack variables and objective functions such that

$x^{k}=\begin{bmatrix}x_{1}^{k}\\\ \vdots\\\
x_{n}^{k}\end{bmatrix}\in\mathbb{R}^{n\times
d}\;\;\text{and}\;\;\nabla\mkern-1.5muF(x^{k})=\begin{bmatrix}\nabla\mkern-1.5muf_{1}(x_{1}^{k})\\\
\vdots\\\
\nabla\mkern-1.5muf_{n}(x_{n}^{k})\end{bmatrix}\in\mathbb{R}^{n\times d}.$

If a fixed point, $x^{\star}$, exists for the signal $x^{k}$, then we denote
the error signal $\tilde{x}^{k}=x^{k}-x^{\star}$.

###  2.2 Assumptions

  1. (A1)

Given $0<\mu\leq L$, we assume that the local gradients are sector bounded on
the interval $(\mu,L)$, meaning that they satisfy the quadratic inequality

$[\star]^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}\begin{bmatrix}-2\mu
LI_{d}&(L+\mu)I_{d}\\\
(L+\mu)I_{d}&-2I_{d}\end{bmatrix}\begin{bmatrix}(x_{i}-x_{\text{opt}})^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}\\\
(\nabla\mkern-1.5muf_{i}(x_{i})-\nabla\mkern-1.5muf_{i}(x_{\text{opt}}))^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}\end{bmatrix}\geq
0$

for all $x_{i}\in\mathbb{R}^{1\times d}$, where $x_{\text{opt}}$ satisfies
$\sum_{i=1}^{n}\nabla\mkern-1.5muf_{i}(x_{\text{opt}})=0$. We define the
condition ratio as $\kappa=L/\mu$, which captures the variation in the
curvature of the objective function.

  2. (A2)

The graph $\mathcal{G}$ is strongly connected.

  3. (A3)

The graph $\mathcal{G}$ is weight balanced, meaning that
$\mathbb{1}^{\mkern-1.5mu{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}\mkern-1.5mu\mathcal{L}=0$.

  4. (A4)

The weights of $\mathcal{G}$ are such that $\sigma=||I-\Pi-\mathcal{L}||<1$.

######  Remark 1.

Assumption (A1) is known as a sector IQC (for a more detailed description see
[11]). (A1) is satisfied when the local objective functions are $\mu$-strongly
convex with $L$-Lipschitz continuous gradients, though the assumption itself
is weaker than those conditions.

######  Remark 2.

Throughout this paper we assume without loss of generality that the dimension
of the local decision and state variables is $d=1$.

######  Remark 3.

Under appropriate conditions on the communications network, the agents can
self-balance their weights in a distributed way to satisfy (A3); for example,
they can use a scalar consensus filter like push-sum (see Algorithm 12 in
[20]). Additionally, agents can normalize their edge weights to enforce that
their weighted in-degrees (and thus their weighted out-degrees) sum to less
than one in order to satisfy (A4).

######  Lemma 1.

Given assumptions (A2) – (A4), the output of the Laplacian is sector bounded
on $(1-\sigma,1+\sigma)$ and equivalently the signals $y$ and $\mathcal{L}y$
satisfy the following quadratic constraint

$[\star]^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}\begin{bmatrix}(\sigma^{2}-1)(I_{n}-\Pi_{n})&(I_{n}-\Pi_{n})\\\
(I_{n}-\Pi_{n})&-(I_{n}-\Pi_{n})\end{bmatrix}\begin{bmatrix}y\\\
\mathcal{L}y\end{bmatrix}\geq 0$ (2)

for all $y\in\mathbb{R}^{n}$.

###### Proof.

From assumptions (A2) \- (A4) and the definition of the spectral norm, we have
that

$\displaystyle\sigma$
$\displaystyle=\max_{y}\frac{||(I-\Pi-\mathcal{L})y||}{||y||}$
$\displaystyle=\max_{y}\frac{||(I-\Pi-\mathcal{L})y||}{||(I-\Pi)y||}$
$\displaystyle\sigma^{2}$
$\displaystyle=\max_{y}\frac{y^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}(I-\Pi-\mathcal{L}^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}})(I-\Pi-\mathcal{L})y}{y^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}(I-\Pi)y}$
$\displaystyle=\max_{y}\frac{y^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}(I-\Pi-\mathcal{L}^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}-\mathcal{L}+\mathcal{L}^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}\mathcal{L})y}{y^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}(I-\Pi)y}$
$\displaystyle\geq\frac{y^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}(I-\Pi-\mathcal{L}^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}-\mathcal{L}+\mathcal{L}^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}\mathcal{L})y}{y^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}(I-\Pi)y}.$

From here, rearranging terms yields (2). ∎

###  2.3 Results

In the following sections we present a parameterized family of distributed,
synchronous, discrete-time algorithms to be be run on each agent such that,
under assumptions (A1)-(A4), we achieve the following:

Accurate convergence:

in the absence of disturbances or other faults, the local estimates $x_{i}$
converge to the optimizer $x_{\text{opt}}$ with an R-linear rate.

Self-healing:

the system state trajectories need not evolve on a pre-defined subspace and
will recover from events such as arbitrary initialization, temporary node
failure, computation errors, or changes in local objectives.

Packet-loss protocol:

if agents are permitted to store a $d$-dimensional memory state for each of
their neighbors, they can implement a packet-loss protocol that allows
computations to continue in the event communication is temporarily lost. This
extends the self-healing of the network to packet loss in a way that is not
possible if the system state trajectories are required to evolve on a pre-
defined subspace.

First we reproduce Algorithm 1 from our prior work [19] along with its
relevant results for stability and convergence to aid in understanding the
convergence analysis of Algorithm 2 (SH-SVL) with packet loss. We analyze
Algorithm 2 (SH-SVL) with packet loss under two scenarios: synchronous packet
loss for general graph topologies and edgewise packet loss where the graph
must be known. We show numerical convergence rates for the later case only for
small directed cycles. Numerical convergence rates are shown for both
algorithms and the case of sub-optimal parameter choice during packet loss is
considered.

##  3 Synthesis of Self-Healing
Distributed Optimization Algorithms

###  3.1 Canonical first-order methods

As a motivation for our algorithms, we use the canonical form first described
in [21] and later used as the SVL template [3]. When the communication graph
is constant, many single-Laplacian methods such as SVL, EXTRA and Exact
Diffusion can be described in this form [4, 9, 10, 21, 3]. Algorithms
representable by the SVL template can also be expressed as a state space
system $G$ in feedback with an uncertain and nonlinear block containing the
objective function gradients $\nabla\mkern-1.5muF(\cdot)$ and the Laplacian
$\mathcal{L}$ shown in Figure 1, where

$\displaystyle G=\left[\begin{array}[]{c|c|c}A&B_{u}&B_{v}\\\ \hline\cr
C_{x}&D_{xu}&D_{xv}\\\ \hline\cr
C_{y}&D_{yu}&D_{yv}\end{array}\right]=\left[\begin{array}[]{cc|c|c}1&\beta&-\alpha&-\gamma\\\
0&1&0&-1\\\ \hline\cr 1&0&0&-\delta\\\ \hline\cr
1&0&0&0\end{array}\right]\otimes I_{n}.$ (10)

$G$${\begin{bmatrix}\nabla\mkern-1.5muF(\cdot)&0\\\
0&\mathcal{L}\end{bmatrix}}$${\begin{bmatrix}u^{k}\\\
v^{k}\end{bmatrix}}$${\begin{bmatrix}x^{k}\\\ y^{k}\end{bmatrix}}$ Figure 1:
Distributed optimization algorithms represented as a feedback interconnection
of an LTI system $G$ and an uncertain block containing the gradients and the
graph Laplacian.

The LTI system $G$ has two states $w_{1}^{k}$ and $w_{2}^{k}$, inputs $u^{k}$
and $v^{k}$, and outputs $x^{k}$ and $y^{k}$ such that

$\displaystyle u_{i}^{k}=\nabla
f_{i}(x_{i}^{k}),\;\;\;v_{i}^{k}=\sum_{j=1}^{n}\mathcal{L}_{ij}^{k}y_{j}^{k}.$
(11)

Additionally, the design parameters of SVL are $\alpha,\beta,\gamma,$ and
$\delta$; $\alpha$ controls the influence of the gradient step, $\gamma$ and
$\delta$ control the influence of the consensus step, and $\beta$ controls the
coupling between the two internal states.

We would like to alert the reader to a small notational difference between our
work and [3]: in this work, the variable $x$ is the input to the gradients and
the variable $y$ is the input to the Laplacian (as shown in Figure 1), whereas
in [3] $y$ is the input to the gradients and $z$ is the input to the
Laplacian. (We cannot use $z$ because we already use it as the frequency
variable of the $z$-transform.)

Algorithms described by Figure 1 and (10) contain two discrete-time
integrators in the LTI block $G$: one integrator is necessary to force the
steady state error to zero, and the other is responsible for the agents coming
to consensus. In the SVL template, the output of the graph Laplacian feeds
into the integrator responsible for consensus.

Algorithms representable by the SVL template, and more broadly all existing
first-order methods with a single Laplacian, require that the system
trajectories evolve on a pre-defined subspace. From our work with average
consensus estimators [16, 17], we know that this drawback arises from the
positional order of the Laplacian and integrator blocks. When the Laplacian
feeds into the integrator, the output of the Laplacian cannot drive the
integrator state away from the consensus subspace, which leads to an
observable but uncontrollable mode. If the integrator state is initialized on
the consensus subspace, or it is otherwise disturbed there, the estimate of
the optimizer will contain an uncorrectable error. Switching the order of the
Laplacian and integrator renders the integrator state controllable but causes
it to become inherently unstable because the integrator output in the
consensus direction is disconnected from the rest of the system. We exploit
this trade-off to develop self-healing distributed optimization algorithms
with only a single Laplacian.

###  3.2 Factorization and integrator location

To switch the order of the Laplacian and the integrator, we first factor an
integrator out of the $G(z)$ block of Figure 1,

$\displaystyle G(z)=\begin{bmatrix}\dfrac{-\alpha}{z-1}&-\dfrac{\delta
z^{2}+(\gamma-2\delta)z+(\beta+\delta-\gamma)}{(z-1)^{2}}\\\\[8.0pt]
\dfrac{-\alpha}{z-1}&-\dfrac{\gamma
z+(\beta-\gamma)}{(z-1)^{2}}\end{bmatrix}\otimes I_{n}$ (12)
$\displaystyle=\begin{bmatrix}\dfrac{-\alpha}{z-1}&-\dfrac{z-1+\zeta}{z-1}\\\\[8.0pt]
\dfrac{-\alpha}{z-1}&\dfrac{-\gamma z-(\beta-\gamma)}{(z-1)(\delta
z+\eta-\delta)}\end{bmatrix}\begin{bmatrix}1&0\\\\[8.0pt] 0&\dfrac{\delta
z+\eta-\delta}{z-1}\end{bmatrix}\otimes I_{n},$ (13)

where

$\displaystyle\eta=\gamma-\delta\zeta\;\;\text{and}\;\;\zeta=\begin{cases}\dfrac{\beta}{\gamma}&\text{if
}\delta=0\\\
\dfrac{\gamma-\sqrt{\gamma^{2}-4\beta\delta}}{2\delta}&\text{otherwise.}\end{cases}$
(14)

Swapping the order of the component matrices yields our new family of
algorithms (where $G_{\mkern-1.5mue}$ replaces $G$):

$\displaystyle G_{\mkern-1.5mue}(z)$
$\displaystyle=\begin{bmatrix}1&0\\\\[8.0pt] 0&\dfrac{\delta
z+\eta-\delta}{z-1}\end{bmatrix}\begin{bmatrix}\dfrac{-\alpha}{z-1}&-\dfrac{z-1+\zeta}{z-1}\\\\[8.0pt]
\dfrac{-\alpha}{z-1}&\dfrac{-\gamma z-(\beta-\gamma)}{(z-1)(\delta
z+\eta-\delta)}\end{bmatrix}\otimes I_{n}$
$\displaystyle=\begin{bmatrix}-\alpha\dfrac{1}{z-1}&-\dfrac{z-1+\zeta}{z-1}\\\\[8.0pt]
-\alpha\dfrac{\delta z+\eta-\delta}{(z-1)^{2}}&-\dfrac{\gamma
z+\beta-\gamma}{(z-1)^{2}}\end{bmatrix}\otimes I_{n}.$ (15)

$\mathcal{L}$${\dfrac{\delta z+\eta-\delta}{z-1}}I_{n}$ Figure 2: The output
of the integrator now feeds into the Laplacian, converting an uncontrollable
and observable mode in the original SVL template to a controllable and
unobservable one.

Now the output of the integrator feeds directly into the Laplacian, as
depicted in Figure 2. We assume that our parameter choices satisfy

$\displaystyle\gamma^{2}$ $\displaystyle\geq 4\beta\delta$ (16)

so that the zeros of $G_{\mkern-1.5mue}$ remain real and thus the system can
be implemented with real-valued signals. The corresponding distributed
algorithm is described in Algorithm 1, where $w_{1}$ and $w_{2}$ are the
internal states of $G_{\mkern-1.5mue}$, and the compact state space form is

$G_{\mkern-1.5mue}=\left[\begin{array}[]{cc|c|c}1&0&-\alpha&-\zeta\\\
1&1&0&-1\\\ \hline\cr 1&0&0&-1\\\
\hline\cr\delta&\eta&0&0\end{array}\right]\otimes I_{n}.$ (17)

The design parameters of Algorithm 1 (and later SH-SVL) are
$\alpha,\delta,\zeta,$ and $\eta$; $\alpha$ controls the influence of the
gradient step, $\zeta$ controls the influence of the consensus step, and
$\delta$ and $\eta$ control the mixture of internal states being transmitted
to neighbors.

######  Remark 4.

The factorization in (13) is not unique; we chose it because it leads to a
method still having only two internal states per agent. There may be other
useful factorizations.

Initialization: Each agent $i\in\\{1,...,n\\}$ chooses
$w_{1i}^{0},w_{2i}^{0}\in\mathbb{R}^{1\times d}$ arbitrarily.
$\mathcal{L}\in\mathbb{R}^{n\times n}$ is the graph Laplacian.

for  _$k=0,1,2,...$_ do

for  _$i\in\\{1,...,n\\}$_ do

Local communication

$y_{i}^{k}=\delta w_{1i}^{k}+\eta w_{2i}^{k}$

[2pt] $v_{i}^{k}=\sum_{j=1}^{n}\mathcal{L}_{ij}y_{j}^{k}$

Local gradient computation

$x_{i}^{k}=w_{1i}^{k}-v_{i}^{k}$ // Update the optimizer estimate.

[2pt] $u_{i}^{k}=\nabla\mkern-1.5muf_{i}(x_{i}^{k})$

Local state update

$w_{1i}^{k+1}=w_{1i}^{k}-\alpha u_{i}^{k}-\zeta v_{i}^{k}$

[2pt] $w_{2i}^{k+1}=w_{1i}^{k}+w_{2i}^{k}-v_{i}^{k}$

end for

end for

return _$x_{i}$_

Algorithm 1 Self-Healing Distributed Gradient Descent

######  Remark 5.

In contrast to algorithms like SVL, Algorithm 1 does not require specific
initial conditions and system trajectories are not restricted to a pre-defined
subspace. If agents change their local objective functions or drop out of the
computation, the system does not need to be reset and the system will converge
to the new minimizer. In the case of agents dropping out, the connection
topology must still be strongly connected, otherwise Algorithm 1 with a weight
balancer will converge to the minimizer for only a subset of the objective
functions and consensus across the network will not be achieved.

##  4 Stability and Convergence Rates Using IQCs

###  4.1 Projection onto the disagreement subspace

As written, our family of algorithms is internally unstable. We use the
projection matrix $(I-\Pi)$ to eliminate the instability from the global
system without affecting $x^{k}$. This procedure is a centralized calculation
that cannot be implemented in a distributed fashion, but it allows us to
analyze the convergence properties of the distributed algorithm.

Consider the steady-state values
$(w_{1}^{\star},x^{\star},u^{\star},v^{\star})$ and suppose $w_{2}^{k}$
contains a component in the $\mathbb{1}$ direction. Then that component does
not affect the aforementioned values because it is an input to the Laplacian
$\mathcal{L}$ (and lies in its nullspace); however, it grows linearly in time
due to the $w_{2}$ update. Thus the system has an internal instability that is
unobservable from the output of the bottom block in Figure 1. Since the
component of $w_{2}^{k}$ in the consensus direction is unobservable to the
variables $(w_{1}^{k},x^{k},u^{k},v^{k})$, we can throw it away without
affecting their trajectories. Using the transformation
$\hat{w}_{2}^{k}=(I-\Pi)w_{2}^{k}$, our state updates become

$\displaystyle w_{1}^{k+1}$ $\displaystyle=w_{1}^{k}-\alpha u^{k}-\zeta v^{k}$
(18) $\displaystyle\hat{w}_{2}^{k+1}$
$\displaystyle=(I-\Pi)w_{1}^{k}+(I-\Pi)\hat{w}_{2}^{k}-(I-\Pi)v^{k}$ (19)
$\displaystyle x^{k}$ $\displaystyle=w_{1}^{k}-v^{k}$ (20)
$\displaystyle\hat{y}^{k}$ $\displaystyle=\delta
w_{1}^{k}+\eta\hat{w}_{2}^{k}$ (21) $\displaystyle u^{k}$
$\displaystyle=\nabla\mkern-1.5muF(x^{k})$ (22) $\displaystyle v^{k}$
$\displaystyle=\mathcal{L}\hat{y}^{k},$ (23)

where $y^{k}$ was replaced with $\hat{y}^{k}$ in (21) and (23) to accommodate
$\hat{w}_{2}^{k}$. These updates lead to the state-space system

$G_{m}=\left[\begin{array}[]{cc|c|c}I&0&-\alpha I&-\zeta I\\\
I-\Pi&I-\Pi&0&-(I-\Pi)\\\ \hline\cr I&0&0&-I\\\ \hline\cr\delta I&\eta
I&0&0\end{array}\right].$ (24)

###  4.2 Existence and optimality of a fixed point

Now that we have eliminated the inherent instability of the global system, we
can state the following about the fixed points:

######  Theorem 1.

For the system described by $G_{m}$, there exists a unique fixed point
$(w_{1}^{\star},\hat{w}_{2}^{\star},x^{\star},\hat{y}^{\star},u^{\star},v^{\star})$,
and the fixed point has $x^{\star}$ in the consensus subspace such that
$x_{i}^{\star}=x_{\textnormal{opt}}$ for all $i\in\\{1,\dots,n\\}$, i.e., the
fixed point of the system is optimal.

###### Proof.

First, assume that the fixed point
$(w_{1}^{\star},\hat{w}_{2}^{\star},x^{\star},\hat{y}^{\star},u^{\star},v^{\star})$
exists. To prove that the variable $x^{\star}$ lies in the consensus
direction, we show that $(I-\Pi)x^{\star}=0$. From (19) and (20) we have that

$\displaystyle(I-\Pi)w_{1}^{\star}$ $\displaystyle=(I-\Pi)v^{\star}$ (25)
$\displaystyle(I-\Pi)x^{\star}$
$\displaystyle=(I-\Pi)w_{1}^{\star}-(I-\Pi)v^{\star}$ (26) $\displaystyle=0.$
(27)

Thus $x_{i}^{\star}=x_{j}^{\star}$ for all $i,j\in\\{1,\dots,n\\}$. Next we
show that $x_{i}^{\star}=x_{\text{opt}}$. From (18) then plugging in (23), we
have

$\displaystyle-\alpha u^{\star}-\zeta v^{\star}$ $\displaystyle=0$ (28)
$\displaystyle u^{\star}$
$\displaystyle=-\frac{\zeta}{\alpha}v^{\star}=-\frac{\zeta}{\alpha}\mathcal{L}\hat{y}^{\star}$
(29)
$\displaystyle\mathbb{1}^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}u^{\star}$
$\displaystyle=-\frac{\zeta}{\alpha}\mathbb{1}^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}\mathcal{L}\hat{y}^{\star}$
(30) $\displaystyle\sum_{i=1}^{n}u_{i}^{\star}$ $\displaystyle=0$ (31)
$\displaystyle\rightarrow\sum_{i=1}^{n}\nabla\mkern-1.5muf_{i}(x_{i}^{\star})$
$\displaystyle=0$ (32) $\displaystyle\rightarrow x_{i}^{\star}$
$\displaystyle=x_{\text{opt}}\;\forall\;i\in\\{1,\dots,n\\}.$ (33)

Thus the fixed point is optimal.

Next, to construct the fixed point we define

$\displaystyle x^{\star}$ $\displaystyle=\mathbb{1}x_{\text{opt}},$
$\displaystyle u^{\star}$ $\displaystyle=\nabla\mkern-1.5muf(x^{\star})$ (34)
$\displaystyle v^{\star}$ $\displaystyle=-\frac{\alpha}{\zeta}u^{\star},$
$\displaystyle w_{1}^{\star}$ $\displaystyle=x^{\star}+v^{\star}.$

Then $\hat{w}_{2}^{\star}$ is the solution to the equation

$\zeta\eta\mathcal{L}\hat{w}_{2}^{\star}=-\alpha(I-\delta\mathcal{L})u^{\star}.$
(35)

Since $\hat{w}_{2}^{k}=\mathcal{L}^{+}\mathcal{L}w_{2}^{k}$ (i.e.,
$\hat{w}_{2}^{k}$ is in the row space of $\mathcal{L}$), we write
$\hat{w}_{2}^{\star}$ in closed form as

$\hat{w}_{2}^{\star}=\frac{\alpha}{\zeta\eta}\mathcal{L}^{+}(\delta\mathcal{L}-I)u^{\star}.$
(36)

Finally, setting $\hat{y}^{\star}=\delta
w_{1}^{\star}+\eta\hat{w}_{2}^{\star}$ completes the proof. ∎

######  Remark 6.

If the graph is switching but converges in time such that the limit of the
sequence of Laplacians exists, as with a weight balancer, then a solution to
(35) still exists and an optimal fixed point can still be found. If the
sequence of Laplacians does not have a limit, then each switch of Laplacian
will introduce a transient to the system, after which the system will converge
to the new optimal fixed point. The dependence of the fixed point on the
Laplacian in (35) is directly related to the order of the integrator block and
the Laplacian block and is an inherent drawback of known self-healing methods.
The fixed points of non-self-healing methods can be independent of the graph
Laplacian, which enables them to converge in the presence of switching
Laplacians [3].

###  4.3 Convergence

Following the approaches in [11, 12, 3], we prove stability using a set of
linear matrix inequalities. First we split our modified system from (24) into
consensus and disagreement components. We define

$\displaystyle A_{m}$ $\displaystyle=A_{p}\otimes\Pi+A_{q}\otimes(I-\Pi)$ (37)
$\displaystyle B_{mu}$ $\displaystyle=B_{pu}\otimes\Pi+B_{qu}\otimes(I-\Pi)$
(38) $\displaystyle B_{mv}$
$\displaystyle=B_{pv}\otimes\Pi+B_{qv}\otimes(I-\Pi)$ (39) $\displaystyle
A_{p}$ $\displaystyle=\begin{bmatrix}1&0\\\ 0&0\end{bmatrix},$ $\displaystyle
B_{pu}$ $\displaystyle=\begin{bmatrix}-\alpha\\\ 0\end{bmatrix},$
$\displaystyle B_{pv}$ $\displaystyle=\begin{bmatrix}-\zeta\\\ 0\end{bmatrix}$
(40) $\displaystyle A_{q}$ $\displaystyle=\begin{bmatrix}1&0\\\
1&1\end{bmatrix},$ $\displaystyle B_{qu}$
$\displaystyle=\begin{bmatrix}-\alpha\\\ 0\end{bmatrix},$ $\displaystyle
B_{qv}$ $\displaystyle=\begin{bmatrix}-\zeta\\\ -1\end{bmatrix}.$

We also define the matrices

$M_{0}=\begin{bmatrix}-2\mu L&L+\mu\\\
L+\mu&-2\end{bmatrix}\;\;\text{and}\;\;M_{1}=\begin{bmatrix}\sigma^{2}-1&1\\\
1&-1\end{bmatrix}.$ (41)

Notice that $M_{0}$ is associated with the sector bound on $\nabla F$ from
(A1) and that $M_{1}$ is associated with the sector bound on $\mathcal{L}$
with inputs from Lemma 1.

We now make a statement analogous to Theorem 10 in [3].

######  Theorem 2.

If there exists $P,Q\in\mathbb{R}^{2\times 2}$,
$\lambda_{0},\lambda_{1}\in\mathbb{R}$, and $\rho\in(0,1)$, with $P,Q\succ 0$
and $\lambda_{0},\lambda_{1}\geq 0$ such that

$\displaystyle[\star]^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}\left[\begin{array}[]{cc|c}P&0&0\\\
0&-\rho^{2}P&0\\\ \hline\cr
0&0&\lambda_{0}M_{0}\end{array}\right]\left[\begin{array}[]{cc}A_{p}&B_{pu}\\\
I&0\\\ \hline\cr C_{x}&D_{xu}\\\ 0&I\end{array}\right]\preceq 0,$ (49)
$\displaystyle[\star]^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}\left[\begin{array}[]{cc|c|c}Q&0&0&0\\\
0&-\rho^{2}Q&0&0\\\ \hline\cr 0&0&\lambda_{0}M_{0}&0\\\ \hline\cr
0&0&0&\lambda_{1}M_{1}\\\
\end{array}\right]\left[\begin{array}[]{ccc}A_{q}&B_{qu}&B_{qv}\\\ I&0&0\\\
\hline\cr C_{x}&D_{xu}&D_{xv}\\\ 0&I&0\\\ \hline\cr C_{y}&D_{yu}&D_{yv}\\\
0&0&I\end{array}\right]\preceq 0,$ (60)

then the following is true for the trajectories of $G_{m}$:

$\Bigg{\|}\begin{bmatrix}w_{1}^{k}-w_{1}^{\star}\\\\[2.84544pt]
\hat{w}_{2}^{k}-\hat{w}_{2}^{\star}\end{bmatrix}\Bigg{\|}^{2}\leq\operatorname{cond}(T)\rho^{2k}\Bigg{\|}\begin{bmatrix}w_{1}^{0}-w_{1}^{\star}\\\\[2.84544pt]
\hat{w}_{2}^{0}-\hat{w}_{2}^{\star}\end{bmatrix}\Bigg{\|}^{2}$ (61)

for a fixed point
$(w_{1}^{\star},\hat{w}_{2}^{\star},x^{\star},\hat{y}^{\star},u^{\star},v^{\star})$,
where $T=P\otimes\Pi_{n}+Q\otimes(I_{n}-\Pi_{n})$. Thus the output $x^{k}$ of
Algorithm 1 converges to the optimizer with the linear rate $\rho$.

###### Proof.

Equation (61) follows directly from Theorem 4 of [11]. First, we define the
quantity
$\xi^{k}=[(w_{1}^{k})^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}\;(\hat{w}_{2}^{k})^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}]^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}$.
We then begin by taking the Kronecker product of (49) with $\Pi_{n}$ and
multiplying on both sides by
$[(\tilde{\xi}^{k})^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}\;(\tilde{u}^{k})^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}]$
and its transpose respectively, as well as taking the Kronecker product of
(60) with $I_{n}-\Pi_{n}$ and multiplying on both sides by
$[(\tilde{\xi}^{k})^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}\;(\tilde{u}^{k})^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}\;(\tilde{v}^{k})^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}]$
and its transpose respectively. Making use of the state space equations
(18)-(23), we get the equations

$\displaystyle(\tilde{\xi}^{k+1})^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}(P\otimes\Pi_{n})\tilde{\xi}^{k+1}-\rho^{2}(\tilde{\xi}^{k})^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}(P\otimes\Pi_{n})\tilde{\xi}^{k}$
(62)
$\displaystyle+\lambda_{0}[\star]^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}(M_{0}\otimes\Pi_{n})\begin{bmatrix}\tilde{x}^{k}\\\
\tilde{u}^{k}\end{bmatrix}\leq 0$
$\displaystyle(\tilde{\xi}^{k+1})^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}(Q\otimes
I_{n}-\Pi_{n})\tilde{\xi}^{k+1}-\rho^{2}(\tilde{\xi}^{k})^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}(Q\otimes
I_{n}-\Pi_{n})\tilde{\xi}^{k}$ (63)
$\displaystyle+\lambda_{0}[\star]^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}(M_{0}\otimes
I_{n}-\Pi_{n})\begin{bmatrix}\tilde{x}^{k}\\\
\tilde{u}^{k}\end{bmatrix}+\lambda_{1}[\star]^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}(M_{1}\otimes
I_{n}-\Pi_{n})\begin{bmatrix}\tilde{\hat{y}}^{k}\\\
\tilde{v}^{k}\end{bmatrix}\leq 0.$

We then add together these equations to get

$\displaystyle(\tilde{\xi}^{k+1})^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}T\tilde{\xi}^{k+1}-\rho^{2}(\tilde{\xi}^{k})^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}T\tilde{\xi}^{k}$
(64)
$\displaystyle+\lambda_{0}[\star]^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}(M_{0}\otimes
I_{n})\begin{bmatrix}\tilde{x}^{k}\\\
\tilde{u}^{k}\end{bmatrix}+\lambda_{1}[\star]^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}(M_{1}\otimes
I_{n}-\Pi_{n})\begin{bmatrix}\tilde{\hat{y}}^{k}\\\
\tilde{v}^{k}\end{bmatrix}\leq 0.$

Notice that the terms associated with $M_{0}$ and $M_{1}$ are quadratic
constraints and nonnegative, therefore we can subtract them from both sides.
Next we multiply by $\rho^{-2k}$ and sum over $k$ which causes a telescoping
sum and yields

$(\tilde{\xi}^{k})^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}T\tilde{\xi}^{k}\leq\rho^{2k}(\tilde{\xi}^{0})^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}T\tilde{\xi}^{0},$
(65)

which leads directly to the statement of (61).

Since the states of $G_{m}$ are converging at a linear rate $\rho$, the rest
of the signals in the system (including $x^{k}$) converge to an optimal fixed
point at the same rate. Additionally, the trajectories of $G_{m}$ and
$G_{\mkern-1.5mue}$ (Algorithm 1) are the same, save for $\hat{w}_{2}^{k}$ and
$\hat{y}^{k}$, so $x^{k}$ in Algorithm 1 also converges to the optimizer with
linear rate $\rho$. ∎

######  Remark 7.

In Algorithm 1, the signal $x^{k}$ converges to the optimizer with linear rate
$\rho$ as the internal signals $w_{2}^{k}$ and $y^{k}$ grow linearly. The
signals $w_{2}^{k}$ and $y^{k}$ will not become large enough to cause overflow
errors until well after the algorithm has converged to the optimizer.

######  Remark 8.

We tested the convergence rates for our algorithm with Zames-Falb IQCs in
place of sector IQCs but saw no improvement.

######  Remark 9.

The explicit dependence of LMIs (49) and (60) on $\mu$ and $L$ can be removed
to facilitate comparison of objective functions based on condition ratio
$\kappa$ alone. First, we define a new signal $g^{k}=\frac{1}{L}u^{k}$ and a
new parameter $\alpha^{\prime}=L\alpha$, such that $\alpha
u^{k}=\alpha^{\prime}g^{k}$. We generate new LMIs that only depend on $\kappa$
by multiplying LMIs (49) and (60) on both sides respectively by the matrices

$\displaystyle\frac{1}{\sqrt{\mu L}}\begin{bmatrix}1&0&0\\\ 0&1&0\\\
0&0&L\end{bmatrix}\text{ and }\frac{1}{\sqrt{\mu L}}\begin{bmatrix}1&0&0&0\\\
0&1&0&0\\\ 0&0&L&0\\\ 0&0&0&1\end{bmatrix}.$

This corresponds to using the IQC

$[\star]^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}\begin{bmatrix}-2&\kappa+1\\\
\kappa+1&-2\kappa\end{bmatrix}\begin{bmatrix}(x_{i}-x_{\text{opt}})^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}\\\
L(\nabla\mkern-1.5muf_{i}(x_{i})-\nabla\mkern-1.5muf_{i}(x_{\text{opt}}))^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}\end{bmatrix}\geq
0.$

If LMIs (49) and (60) are feasible then these new LMIs will be feasible as
well with $P^{\prime}=\frac{1}{\mu L}P$ and $Q^{\prime}=\frac{1}{\mu L}Q$.

##  5 Robustness to Packet Loss

###  5.1 Packet-loss protocol

We next give our agents some additional memory so that they can substitute
previously transmitted values when a packet is lost. Each agent
$i\in\\{1,\dots,n\\}$ maintains an edge state $r_{ij}^{k}$ for each
$j\in\mathcal{N}_{\text{in}}(i)$ (the set of neighbors who transmit to $i$).
Whenever agent $i$ receives a message from agent $j$, it updates the state
$r_{ij}$ accordingly; however, if at time $k$ no message from neighbor $j$ is
received, agent $i$ must estimate what would have likely been transmitted. One
potential strategy is to substitute in the last message received, but because
$y_{j}$ is growing linearly in quasi steady state, this naive strategy would
ruin steady-state accuracy. Instead we must account for the linear growth
present in our algorithm, which we can do by analyzing the quantity
$y_{j}^{k}-y_{j}^{k-1}$ at the quasi fixed point
$(w_{1}^{\star},x^{\star},u^{\star},v^{\star})$:

$\displaystyle y_{j}^{k}-y_{j}^{k-1}$
$\displaystyle=\delta(w_{1j}^{\star}-w_{1j}^{\star})+\eta(w_{2j}^{k}-w_{2j}^{k-1})$
(66) $\displaystyle=\eta(w_{1j}^{\star}-v_{j}^{\star})$ (67)
$\displaystyle=\eta x_{j}^{\star}\approx\eta x_{i}^{k}$ (68)

Therefore, when a packet is not received by a neighbor, agent $i$ scales its
estimate of the optimizer and adds it to its previously received (or
estimated) message. The packet-loss protocol is summarized in Algorithm 2 (SH-
SVL). By construction, the modifications included in SH-SVL will not alter the
quasi fixed points of Algorithm 1 and, in the absence of dropped packets, the
state trajectories of SH-SVL are equivalent to those of Algorithm 1. In the
following sections, we present a similar analysis to that of the base
algorithm for the simplified case of synchronous packet loss with an unknown
graph and the general case of edgewise packet loss where knowledge of the
graph topology must be included in the analysis.

######  Remark 10.

SH-SVL can be modified to include a forgetting factor $q$. If agent $i$ does
not receive a packet from neighbor $j$ in $q$ time steps, then agent $i$
assumes that the communication link has been severed and clears $r_{ij}$ from
memory.

Initialization: Each agent $i\in\\{1,...,n\\}$ chooses
$w_{1i}^{0},w_{2i}^{0}\in\mathbb{R}^{1\times d}$ arbitrarily.
$\mathcal{L}\in\mathbb{R}^{n\times n}$ is the graph Laplacian. All $r_{ij}$
are initialized the first time a message is received from a neighbor.

for  _$k=0,1,2,...$_ do

for  _$i\in\\{1,...,n\\}$_ do

Local communication

$y_{i}^{k}=\delta w_{1i}^{k}+\eta w_{2i}^{k}$

[2pt] for  _$j\in\mathcal{N}_{\text{in}}(i)\cup i$_ do

if  _Packet from $j$ received by $i$_ then

$r_{ij}^{k}=y_{j}^{k}$

else

$r_{ij}^{k}=\eta x_{i}^{k-1}+r_{ij}^{k-1}$

end if

end for

$v_{i}^{k}=\sum_{j=1}^{n}\mathcal{L}_{ij}r_{ij}^{k}$

Local gradient computation

$x_{i}^{k}=w_{1i}^{k}-v_{i}^{k}$ // Update the optimizer estimate.

[2pt] $u_{i}^{k}=\nabla\mkern-1.5muf_{i}(x_{i}^{k})$

Local state update

$w_{1i}^{k+1}=w_{1i}^{k}-\alpha u_{i}^{k}-\zeta v_{i}^{k}$

[2pt] $w_{2i}^{k+1}=w_{1i}^{k}+w_{2i}^{k}-v_{i}^{k}$

end for

end for

return $x_{i}$

Algorithm 2 Self-Healing Distributed Gradient Descent with Packet-loss
protocol (SH-SVL)

###  5.2 Synchronous Loss

The case of synchronous packet loss is applicable to noisy channels where all
communications are affected simultaneously, such as in a storm or other wide-
ranging phenomena. Since all communications are affected simultaneously, the
specific topology of the network is allowed to be unknown, as in the base
algorithm, which allows for the analysis of algorithm performance over a range
of connectivity parameters and sizes instead of specific graph topologies. At
each time step $k$, we assume packets are lost in the entire network with an
independent and identically distributed probability $p_{\ell}$. We note that
the following analysis can be adapted to the case where the packet loss
probability follows a Markov chain; we refer the reader to the body of work on
Markov jump linear systems for more information [22].

To begin the analysis, we first assume that each agent maintains a memory
state for its own output and considers its own output lost when all packets in
the network are lost. For the purposes of the centralized analysis, we assume
without loss of generality that each agent $i$ maintains a state for every
node in the network and stores it in a vector $r_{i}^{k}$. Thus if $l$ packets
are lost consecutively at time $k>l$, we have

$\displaystyle r_{i}^{k}=y^{k-l}+\eta\sum_{j=1}^{l}x_{i}^{k-j}\mathbb{1}$ (69)
$\displaystyle
v_{i}^{k}=\mathcal{L}_{i}^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}r_{i}^{k}=\mathcal{L}_{i}^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}y^{k-1}.$
(70)

Since the $x_{i}$ terms are annihilated by the Laplacian, in the analysis we
can replace all vectors $r_{i}^{k}$ with a single vector $\hat{r}^{k}$, whose
update is given by:

$\hat{r}^{k+1}=\begin{cases}\delta w_{1}^{k+1}+\eta\hat{w}_{2}^{k+1}&\text{if
Success}\\\ \hat{r}^{k}&\text{if Loss}.\end{cases}$ (71)

From Equation (71), it is straightforward to see that the fixed point from
Theorem 1 is maintained with $\hat{r}^{\star}=\delta
w_{1}^{\star}+\eta\hat{w}_{2}^{\star}$ and that the trajectories of $x^{k}$
are still unaffected by the projection of $\hat{w}_{2}^{k}$. With these
modification, the system matrix at each time is
$G_{\iota_{k}}\in\\{G_{s},G_{\ell}\\}$ where $G_{s}$ corresponds to successful
packet transmission in the network (with probability $1-p_{\ell}$) and
$G_{\ell}$ corresponds to synchronous packet loss in the network (with
probability $p_{\ell}$). They are as follows:

$\displaystyle A_{\iota_{k}}$
$\displaystyle=A_{\iota_{k}p}\otimes\Pi+A_{\iota_{k}q}\otimes(I-\Pi),$
$\displaystyle B_{\iota_{k}u}$
$\displaystyle=B_{\iota_{k}pu}\otimes\Pi+B_{\iota_{k}qu}\otimes(I-\Pi),$
$\displaystyle B_{\iota_{k}v}$
$\displaystyle=B_{\iota_{k}pv}\otimes\Pi+B_{\iota_{k}qv}\otimes(I-\Pi),$
$\displaystyle A_{sp}$ $\displaystyle=\begin{bmatrix}1&0&0\\\ 0&0&0\\\
\delta&0&0\end{bmatrix},$ $\displaystyle\,B_{spu}$
$\displaystyle=\begin{bmatrix}-\alpha\\\ 0\\\ -\alpha\delta\end{bmatrix},$
$\displaystyle\,B_{spv}$ $\displaystyle=\begin{bmatrix}-\zeta\\\ 0\\\
-\delta\zeta\end{bmatrix},$ $\displaystyle A_{sq}$
$\displaystyle=\begin{bmatrix}1&0&0\\\ 1&1&0\\\
\delta+\eta&\eta&0\end{bmatrix},$ $\displaystyle B_{squ}$
$\displaystyle=\begin{bmatrix}-\alpha\\\ 0\\\ -\alpha\delta\end{bmatrix},$
$\displaystyle B_{sqv}$ $\displaystyle=\begin{bmatrix}-\zeta\\\ -1\\\
-\delta\zeta-\eta\end{bmatrix},$

$\displaystyle A_{\ell p}$ $\displaystyle=\begin{bmatrix}1&0&0\\\ 0&0&0\\\
0&0&1\end{bmatrix},$ $\displaystyle\,B_{\ell pu}$
$\displaystyle=\begin{bmatrix}-\alpha\\\ 0\\\ 0\end{bmatrix},$
$\displaystyle\,B_{\ell pv}$ $\displaystyle=\begin{bmatrix}-\zeta\\\ 0\\\
0\end{bmatrix},$ $\displaystyle A_{\ell q}$
$\displaystyle=\begin{bmatrix}1&0&0\\\ 1&1&0\\\ 0&0&1\end{bmatrix},$
$\displaystyle B_{\ell qu}$ $\displaystyle=\begin{bmatrix}-\alpha\\\ 0\\\
0\end{bmatrix},$ $\displaystyle B_{\ell qv}$
$\displaystyle=\begin{bmatrix}-\zeta\\\ -1\\\ 0\end{bmatrix},$

$\displaystyle C_{x}$ $\displaystyle=\begin{bmatrix}I&0&0\end{bmatrix},$
$\displaystyle\quad C_{y}$ $\displaystyle=\begin{bmatrix}0&0&I\end{bmatrix},$
$\displaystyle\quad D_{xv}$ $\displaystyle=-I.$

$G_{\iota_{k}}=\left[\begin{array}[]{c|c|c}A_{\iota_{k}}&B_{\iota_{k}u}&B_{\iota_{k}v}\\\
\hline\cr C_{x}&0&D_{xv}\\\ \hline\cr C_{y}&0&0\end{array}\right].$ (72)

######  Theorem 3.

If packets are lost synchronously with an i.i.d. probability $p_{\ell}$ and
there exists $P,Q\in\mathbb{R}^{3\times 3}$,
$\lambda_{0},\lambda_{1}\in\mathbb{R}$, and $\rho\in(0,1)$, with $P,Q\succ 0$
and $\lambda_{0},\lambda_{1}\geq 0$ such that

$\displaystyle\sum_{o\in\\{s,l\\}}p_{o}[\star]^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}\left[\begin{array}[]{cc|c}P&0&0\\\
0&-\rho^{2}P&0\\\ \hline\cr
0&0&\lambda_{0}M_{0}\end{array}\right]\left[\begin{array}[]{cc}A_{op}&B_{opu}\\\
I&0\\\ \hline\cr C_{x}&0\\\ 0&I\end{array}\right]\preceq 0,$ (80)
$\displaystyle\sum_{o\in\\{s,l\\}}p_{o}[\star]^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}\left[\begin{array}[]{cc|c|c}Q&0&0&0\\\
0&-\rho^{2}Q&0&0\\\ \hline\cr 0&0&\lambda_{0}M_{0}&0\\\ \hline\cr
0&0&0&\lambda_{1}M_{1}\\\ \end{array}\right]\times$ (85)
$\displaystyle\left[\begin{array}[]{ccc}A_{oq}&B_{oqu}&B_{oqv}\\\ I&0&0\\\
\hline\cr C_{x}&0&D_{xv}\\\ 0&I&0\\\ \hline\cr C_{y}&0&0\\\
0&0&I\end{array}\right]\preceq 0,$ (92)

then the following is true for the trajectories of $G_{\iota_{k}}$:

$\mathbb{E}\Bigg{[}\Bigg{\|}\begin{bmatrix}w_{1}^{k}-w_{1}^{\star}\\\\[2.84544pt]
\hat{w}_{2}^{k}-\hat{w}_{2}^{\star}\\\\[2.84544pt]
r^{k}-r^{\star}\end{bmatrix}\Bigg{\|}^{2}\Bigg{]}\leq\operatorname{cond}(T)\rho^{2k}\Bigg{\|}\begin{bmatrix}w_{1}^{0}-w_{1}^{\star}\\\\[2.84544pt]
\hat{w}_{2}^{0}-\hat{w}_{2}^{\star}\\\\[2.84544pt]
r^{0}-r^{\star}\end{bmatrix}\Bigg{\|}^{2}$ (93)

for a fixed point
$(w_{1}^{\star},\hat{w}_{2}^{\star},r^{\star},x^{\star},\hat{y}^{\star},u^{\star},v^{\star})$,
where $T=P\otimes\Pi_{n}+Q\otimes(I_{n}-\Pi_{n})$. Thus the output $x^{k}$ of
SH-SVL converges in mean square to the optimizer with the linear rate $\rho$
and additionally converges almost surely.

###### Proof.

The proof of Theorem 3 follows directly from Theorem 2 above by redefining the
quantity
$\xi^{k}=[(w_{1}^{k})^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}\;(\hat{w}_{2}^{k})^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}\;(r^{k})^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}]^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}$
and taking expectations where appropriate. See also Theorem 1 in [13]. Almost
sure convergence can be seen by summing both sides of Equation (93) over $k$,
then applying Markov’s inequality and the Borel-Cantelli lemma. ∎

######  Remark 11.

The same system $G_{\iota_{k}}$ can be used in a periodic loss setting to
analyze the case when multiple local computation steps are taken for every
communication step. The setup follows directly from the results in [23];
however we do not include it here because the convergence rate is slower than
the setting where one computation step is performed for every communication
step.

###  5.3 Edgewise Packet Loss

In the general case of edgewise packet loss, packets are allowed to be lost
according to a Markovian process on the edges. For our analysis, we consider
the case where packets are lost with an independent and identically
distributed probability $p_{\ell}$ on each edge; however, the analysis can be
extended in a straightforward manner using any Markovian packet loss process
desired. First, we must define some additional quantities. Let

$\displaystyle E=\begin{bmatrix}\mathcal{L}_{r_{1}}&0&\dots&0\\\
0&\mathcal{L}_{r_{2}}&\dots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&0&\dots&\mathcal{L}_{r_{n}}\end{bmatrix},\;\phi^{k}=\begin{bmatrix}r_{1}^{k}\\\\[2.84544pt]
r_{2}^{k}\\\ \vdots\\\ r_{n}^{k}\end{bmatrix}.$ (94)

Then the state space equations for SH-SVL are

$\displaystyle w_{1}^{k+1}=w_{1}^{k}-\zeta E\phi^{k}-\alpha u^{k}$ (95)
$\displaystyle w_{2}^{k+1}=w_{1}^{k}+w_{2}^{k}-E\phi^{k}$ (96)
$\displaystyle\phi_{(i,j)}^{k+1}=\begin{cases}\delta w_{1j}^{k+1}+\eta
w_{2j}^{k+1}&\text{if }(i,j)\text{ success}\\\ \phi_{(i,j)}^{k}+\eta
x_{i}^{k}&\text{if }(i,j)\text{ loss}\end{cases}$ (97) $\displaystyle
x^{k}=w_{1}^{k}-E\phi^{k}$ (98) $\displaystyle u^{k}=\nabla F(x).$ (99)

To proceed with convergence analysis, we must project the unstable state
$w_{2}^{k}$ onto the disagreement subspace. The state space equations with the
projection then become

$\displaystyle w_{1}^{k+1}=w_{1}^{k}-\zeta E\hat{\phi}^{k}-\alpha u^{k}$ (100)
$\displaystyle\hat{w}_{2}^{k+1}=(I-\Pi)w_{1}^{k}+(I-\Pi)\hat{w}_{2}^{k}-(I-\Pi)E\hat{\phi}^{k}$
(101) $\displaystyle\hat{\phi}_{(i,j)}^{k+1}=\begin{cases}\delta
w_{1j}^{k+1}+\eta\hat{w}_{2j}^{k+1}&\text{if }(i,j)\text{ success}\\\
\phi_{(i,j)}^{k}+(I-\Pi)_{r_{i}}x^{k}&\text{if }(i,j)\text{ loss}\end{cases}$
(102) $\displaystyle x^{k}=w_{1}^{k}-E\hat{\phi}^{k}$ (103) $\displaystyle
u^{k}=\nabla F(x).$ (104)

Now we show that the trajectories for $x^{k}$ are identical for the systems
defined by both Equations (95) and (100), and that the fixed point for
Equation (100) is the same as in Theorem 1.

######  Theorem 4.

The trajectories for $w_{1}^{k}$ and $x^{k}$ are unchanged from the
unprojected system (95) to the projected system (100).

###### Proof.

We will prove the theorem by induction. First assume that the initial
conditions for $w_{1}^{0}$ are the same for the two systems. Furthermore,
assume that

$\displaystyle\hat{w}_{2}^{0}=(I-\Pi)w_{2}^{0}$ (105)
$\displaystyle\hat{\phi}_{(i,j)}^{0}=\phi_{(i,j)}^{0}-\frac{\eta}{n}\mathbb{1}^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}w_{2}^{0}.$
(106)

Since the projected system is inherently centralized, we can make these
modifications to the initial conditions without loss of generality. It is
clear from the initial conditions and the properties of the Laplacian that
$E\phi^{0}=E\hat{\phi}^{0}$. From the update equations we can see that $x^{0}$
and $w_{1}^{1}$ are the same for both systems and that
$\hat{w}_{2}^{1}=(I-\Pi)w_{2}^{1}$. Finally, looking at the update for
$\hat{\phi}$, we can see that

$\displaystyle\hat{\phi}^{1}_{(i,j)}$ $\displaystyle=\begin{cases}\delta
w_{1j}^{1}+\eta
w_{2j}^{1}-\frac{\eta}{n}\mathbb{1}^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}w_{2}^{1}&\text{if
}(i,j)\text{ success}\\\ (\phi_{(i,j)}^{0}+\eta
x_{i}^{0})-\frac{\eta}{n}\mathbb{1}^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}(w_{2}^{0}+x^{0})&\text{if
}(i,j)\text{ loss}\end{cases}$ $\displaystyle=\begin{cases}\delta
w_{1j}^{1}+\eta
w_{2j}^{1}-\frac{\eta}{n}\mathbb{1}^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}w_{2}^{1}&\text{if
}(i,j)\text{ success}\\\ (\phi_{(i,j)}^{0}+\eta
x_{i}^{0})-\frac{\eta}{n}\mathbb{1}^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}w_{2}^{1}&\text{if
}(i,j)\text{ loss}\end{cases}$
$\displaystyle=\phi^{1}_{(i,j)}-\frac{\eta}{n}\mathbb{1}^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}w_{2}^{1}.$

This completes the proof.

∎

Following directly from Theorem 1, the fixed point for the projected system is

$\displaystyle x^{\star}$ $\displaystyle=\mathbb{1}x_{\text{opt}},$
$\displaystyle u^{\star}$ $\displaystyle=\nabla\mkern-1.5muf(x^{\star})$ (107)
$\displaystyle w_{1}^{\star}$ $\displaystyle=x^{\star}+v^{\star},$
$\displaystyle\hat{w}_{2}^{\star}$
$\displaystyle=\frac{\alpha}{\zeta\eta}\mathcal{L}^{+}(\delta\mathcal{L}-I)u^{\star}$
$\displaystyle\hat{\phi}^{\star}_{(i,j)}$ $\displaystyle=\delta
w_{1j}^{\star}+\eta\hat{w}_{2j}^{\star}.$

In order for the system to be efficiently implemented, we must exploit the
fact that the state $\hat{\phi}_{(i,j)}^{k}$ is always a success and absorb it
into the state space equations. First, we define the modified quantities
$\bar{E}$ and $\hat{\bar{\phi}}^{k}$ as follows

$\displaystyle\bar{E}=\begin{bmatrix}\bar{\mathcal{L}}_{r_{1}}&0&\dots&0\\\
0&\bar{\mathcal{L}}_{r_{2}}&\dots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&0&\dots&\bar{\mathcal{L}}_{r_{n}}\end{bmatrix},\;\hat{\bar{\phi}}^{k}=[\hat{\phi}_{(i,j)}^{k}]_{(i,j)\in\mathcal{E}},$
(108)

where $\bar{E}\in\mathbb{R}^{n,m}$, $\hat{\bar{\phi}}^{k}\in\mathbb{R}^{m}$,
and $\bar{\mathcal{L}}_{i}$ is a row vector whose elements are
$\mathcal{L}_{ij}$ if $(i,j)\in\mathcal{E}$. Note that the trajectories of the
system have not been modified but the state space is now of size $2n+m$.

Consider the collections

$\displaystyle\omega^{k}$ $\displaystyle=\\{\tau:\tau\in\mathcal{E}\text{ and
}\tau\text{ successful at time }k\\}$ (109) $\displaystyle\Omega$
$\displaystyle=\\{\omega^{k}:k\in\mathbb{Z}_{\geq 0}\\}.$ (110)

Note that $|\Omega|=2^{m}.$ Then we have that

$\mathbb{P}(\omega^{k}=\omega)=p_{\omega}=(1-p_{\ell})^{|\omega|}p_{\ell}^{m-|\omega|}.$
(111)

Now we can define the state matrices as follows:

$\displaystyle A_{1}$ $\displaystyle=\begin{bmatrix}I-\zeta\delta
D(\mathcal{L})&-\zeta\eta D(\mathcal{L})&-\zeta E\end{bmatrix}$ (112)
$\displaystyle A_{2}$ $\displaystyle=(I-\Pi)\begin{bmatrix}(I-\delta
D(\mathcal{L}))&(I-\eta D(\mathcal{L}))&-E\end{bmatrix}$ (113) $\displaystyle
B_{1}$ $\displaystyle=-\alpha I$ (114) $\displaystyle C$
$\displaystyle=\begin{bmatrix}(I-\delta D(\mathcal{L}))&-\eta
D(\mathcal{L})&-E\end{bmatrix}.$ (115)

If $\tau=(i,j)\in\omega^{k}$ (i.e., edge $\tau$ is successful):

$\displaystyle A_{\omega_{k}\phi_{\tau}1}$
$\displaystyle=e_{j}^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}[\delta(I-\zeta\delta
D(\mathcal{L}))+\eta(I-\Pi)(I-\delta D(\mathcal{L}))]$ (116) $\displaystyle
A_{\omega_{k}\phi_{\tau}2}$
$\displaystyle=e_{j}^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}[\eta(I-\Pi)(I-\eta
D(\mathcal{L})-\delta\zeta\eta D(\mathcal{L})]$ (117) $\displaystyle
A_{\omega_{k}\phi_{\tau}\phi}$
$\displaystyle=e_{j}^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}[-(\delta\zeta
I+\eta(I-\Pi))E]$ (118) $\displaystyle B_{\omega_{k}\phi_{\tau}}$
$\displaystyle=-\alpha\delta
e_{j}^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}.$
(119)

If $\tau=(i,j)\in\mathcal{E}\setminus\omega^{k}$ (i.e., edge $\tau$ is
unsuccessful):

$\displaystyle A_{\omega_{k}\phi_{\tau}1}$
$\displaystyle=e_{i}^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}\eta(I-\Pi)(I+\delta
D(\mathcal{L}))$ (120) $\displaystyle A_{\omega_{k}\phi_{\tau}2}$
$\displaystyle=e_{i}^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}\eta^{2}(I-\Pi)D(\mathcal{L})$
(121) $\displaystyle A_{\omega_{k}\phi_{\tau}\phi}$
$\displaystyle=e_{\tau}^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}+e_{i}^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}\eta(I-\Pi)E$
(122) $\displaystyle B_{\omega_{k}\phi_{\tau}}$ $\displaystyle=0.$ (123)

Finally,

$\displaystyle A_{\omega_{k}\phi}$
$\displaystyle=\begin{bmatrix}A_{\omega_{k}\phi_{\tau}1}&A_{\omega_{k}\phi_{\tau}2}&A_{\omega_{k}\phi_{\tau}\phi}\end{bmatrix}_{\tau\in\mathcal{E}}$
(124) $\displaystyle B_{\omega_{k}\phi}$
$\displaystyle=[B_{\omega_{k}\phi_{\tau}}]_{\tau\in\mathcal{E}}$ (125)

$G_{\omega_{k}}=\left[\begin{array}[]{c|c}A_{1}&B_{1}\\\ A_{2}&0\\\
A_{\omega_{k}\phi}&B_{\omega_{k}\phi}\\\ \hline\cr C&0\end{array}\right].$
(126)

######  Theorem 5.

If packets are lost edgewise with an i.i.d. probability $p_{\ell}$ and there
exists $P\in\mathbb{R}^{2n+m\times 2n+m}$, $\lambda_{0}\in\mathbb{R}$, and
$\rho\in(0,1)$, with $P\succ 0$ and $\lambda_{0}\geq 0$ such that

$\displaystyle\sum_{\omega\in\Omega}p_{\omega}[\star]^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}\left[\begin{array}[]{cc|c}P&0&0\\\
0&-\rho^{2}P&0\\\ \hline\cr
0&0&\lambda_{0}M_{0}\end{array}\right]\left[\begin{array}[]{cc}A_{\omega}&B_{\omega}\\\
I&0\\\ \hline\cr C&0\\\ 0&I\end{array}\right]\preceq 0,$ (134)

then the following is true for the trajectories of $G_{\omega_{k}}$:

$\mathbb{E}\Bigg{[}\Bigg{\|}\begin{bmatrix}w_{1}^{k}-w_{1}^{\star}\\\\[2.84544pt]
\hat{w}_{2}^{k}-\hat{w}_{2}^{\star}\\\\[2.84544pt]
\hat{\bar{\phi}}^{k}-\hat{\bar{\phi}}^{\star}\end{bmatrix}\Bigg{\|}^{2}\Bigg{]}\leq\operatorname{cond}(P)\rho^{2k}\Bigg{\|}\begin{bmatrix}w_{1}^{0}-w_{1}^{\star}\\\\[2.84544pt]
\hat{w}_{2}^{0}-\hat{w}_{2}^{\star}\\\\[2.84544pt]
\hat{\bar{\phi}}^{0}-\hat{\bar{\phi}}^{\star}\end{bmatrix}\Bigg{\|}^{2}$ (135)

for a fixed point
$(w_{1}^{\star},\hat{w}_{2}^{\star},\hat{\bar{\phi}}^{\star},x^{\star},u^{\star})$.
Thus the output $x^{k}$ of SH-SVL converges almost surely to the optimizer
with the linear rate $\rho$.

###### Proof.

The result follows directly from Theorem 3. ∎

##  6 Numerical Results

The linear matrix inequalities provided in Theorems 2, 3, and 5 must be solved
numerically in order to certify a worst-case convergence rate $\rho$. To do
so, the associated semidefinite program was implemented in the Julia
programming language using JuMP and solved with MOSEK [24]. Given a set of
parameters, a bisection search was used to find $\rho$. The Nelder-Mead
algorithm (from the Optim library [25]) was used to find algorithm parameters
($\alpha,\delta,\zeta,\eta$) that minimized $\rho$. We found that the Nelder-
Mead algorithm needed to be initialized with parameters that certified a
converging algorithm to find useful parameters and that it could effectively
be warm started by providing near-optimal parameters. The worst-case
convergence rate of our algorithm is subject to the same lower bound,
$\rho\geq\max((\kappa-1)/(\kappa+1),\sigma)$, found in [3].

###  6.1 Convergence Rates

Figure 3 and 4 show the performance of SH-SVL in the synchronous loss setting
over a range of graph connectivity and for multiple rates of packet loss,
including the case with no packet loss where Algorithm 1 and 2 are equivalent,
for condition ratios of $\kappa=10$ and $\kappa=20$ respectively. The figures
also include the performance of SVL and the theoretical lower bound [3] on
performance for reference. Notably, in the lossless setting with parameters
optimized for a given $\kappa$ and $\sigma$, SH-SVL marginally outperforms the
SVL algorithm in the worst case. This could be due to the fact that the SVL
algorithm does not optimize for $\rho$ directly. Also of note is that when the
parameters are optimized for the specific level of packet loss, SH-SVL
converges even in the presence of extraordinarily high packet loss on graphs
with very low connectivity (the value of $\sigma$ is high).

Figure 3: Performance of SH-SVL in the synchronous packet loss setting for a
range of independent packet loss chances. SVL included for reference. The
objective function condition ratio is $\kappa=10$. Figure 4: Performance of
SH-SVL in the synchronous packet loss setting for a range of independent
packet loss chances. SVL included for reference. The objective function
condition ratio is $\kappa=20$.

However, if the rate of packet loss is higher than the parameters are
optimized for, the performance of SH-SVL can be very fragile. Figure 5 depicts
the results when the optimal parameters for the lossless case are used in the
lossy case. Interestingly, if the packet loss rate is lower than that which
the parameters are optimized for, the worst case performance remains the same.
Figure 6 depicts the results when optimal parameters for $p_{\ell}=0.5$ are
used. This would imply that system designers could decide the tradeoff between
typical operating convergence speeds and the maximum allowable packet loss in
the system.

Figure 5: Performance of SH-SVL in the synchronous packet loss setting for a
range of independent packet loss chances. Parameters are optimized for
$\kappa=10$ and given $\sigma$ value but no packet loss. Figure 6: Performance
of SH-SVL in the synchronous packet loss setting for a range of independent
packet loss chances. Parameters are optimized for $\kappa=10$, $p_{\ell}=0.5$,
and the given $\sigma$ value.

Figure 7 shows the performance of SH-SVL in both the edgewise and synchronous
loss settings over the range of packet loss probabilities and for different
sizes of directed cycles. The theoretical lower bound is included for
reference. Directed cycles were chosen as the graph topology to demonstrate
the performance of SH-SVL because they have the minimum number of edges for
the graph to be connected, which makes the size of the semidefinite program
smaller and also gives them a very low connectivity which makes them serve as
a rough lower bound on the performance for a given graph size. Only sizes up
to $n=5$ are shown due to the fact that finding optimal parameters using
Nelder-Mead for such large problems starts to take a prohibitively long time.
The 3-, 4-, and 5-node directed cycles have edge weights of $0.9$ and
$\sigma=[0.854,0.906,0.936]$ respectively. From Figures 3 and 7, it is clear
that the LMI from Theorem 6 yields a faster convergence rate than the LMIs
from Theorem 3 for the same values of $\kappa$ and $p_{\ell}$. This is
entirely due to the fact that the LMIs in Theorem 6 are less conservative
because the Laplacian is absorbed into the linear system whereas in Theorem 3,
the Laplacian is parameterized. Similar LMIs to Theorem 6 can be derived for
the synchronous loss case (results plotted as dashed lines in Figure 7); the
results suggest that the synchronous loss case is a good upper bound on
performance for the edgewise loss setting.

Figure 7: Performance of SH-SVL in the edgewise and synchronous packet loss
setting on directed cycles with the Laplacian absorbed into the linear
dynamics. Edgewise loss resultst are solid lines while synchronous loss
results are in dashed lines. $\kappa=10$.

###  6.2 Classification example

To test the performance of our algorithm under packet loss, we solved a
classification problem using the COSMO chip dataset [26]. The problem setup is
as follows: a network of $n$ agents would like to collaboratively compute a
binary classifier that identifies whether or not a computer chip will pass a
quality assurance test using data gathered independently by each agent in the
network; furthermore, they would like to do so in a distributed fashion
without sharing their datasets. In order to simulate this problem, we divide
the COSMO chip dataset into $n$ local subsets and denote agent $i$’s set of
local data indices as $S_{i}$. To improve the performance of their classifier,
the agents employ a polynomial embedding where each 2-dimensional data point
$d_{j}=[d_{j1},d_{j2}]$ is embedded in a 28-dimensional space given by

$M(d_{j})=[1,d_{j1},d_{j2},d_{j1}^{2},d_{j1}d_{j2},d_{j2}^{2},d_{j1}^{3},\dots,d_{j1}d_{j2}^{5},d_{j2}^{6}].$

The agents use the logistic loss function with $L_{2}$-regularization,
yielding local cost functions given by

$f_{i}(x_{i})=\sum_{j\in
S_{i}}\log(1+e^{-l_{j}x_{i}^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}M(d_{j})})+\frac{1}{n}||x_{i}||^{2},$
(136)

where $l_{j}\in\\{-1,1\\}$ is the label of data point $j$, and each agent
computes the label of unseen data by using the operation
$l=\operatorname{sgn}(x_{i}^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}M(d))$.
Using the cost in (136), the corresponding sector bound $(\mu,L)$ is
approximated as $\mu=\frac{2}{n}$ and

$\displaystyle
L_{i}\leq\Big{\|}\frac{2}{n}I+\frac{1}{4}M_{i}^{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}M_{i}\Big{\|},\;\;L=\max_{i}L_{i},$
(137)

where the rows of $M_{i}$ are $M(d_{j})$ for $j\in S_{i}$. The agents’
connection topology is described by an $n=7$ node directed ring lattice, shown
in Figure 8, such that $(i,j)\in\mathcal{E}$ when $j\in\\{i+1,i+3,i+5\\}\mod
n$. All edge weights in the graph are set to 1/4 and
$\sigma=||I-\Pi-\mathcal{L}||=0.562$.

$1$$2$$3$$4$$5$$6$$7$ Figure 8: The directed network topology for the
classification example. All edge weights are 1/4.

Using $(\mu,L)$ and $\sigma$, we computed the optimal parameters for our
algorithm using Nelder-Mead assuming no packet loss and computed the SVL
parameters as detailed in [3]. We then simulated both SH-SVL and SVL with and
without packet loss and took the maximum error between the distributed
algorithms and a centralized solution found using Convex.jl and MOSEK. We ran
SH-SVL using random initial conditions on the interval $[0,1]$ and the SVL
algorithm using zero initial conditions. For the packet loss run of SVL, we
held the previous message on each edge so that the fixed points would be
unaffected. At each time step, packets had a 30% chance of being lost,
independent of each other and time. The results of these simulations are shown
in Figure 9. In this scenario, SH-SVL has similar performance as SVL and SH-
SVL with lossy channels still converges to the optimum at a similar rate as
SH-SVL with lossless channels, despite the high rate of packet loss that
causes SVL to converge with high error. The worst-case convergence rate for
both SH-SVL and SVL is $\rho=0.9544$ and the actual rates in this scenario are
$0.9514$ and $0.9455$ respectively. With a packet-loss rate of $30\%$ and
parameters optimized for lossless channels, the worst-case convergence rate
for SH-SVL is $\rho=1$; however, the actual rate from Figure 9 is $0.9632$,
demonstrating the fact that the worst-case results for packet loss can be very
conservative and that the optimal parameters may not be as fragile as Figure 5
suggests.

Figure 9: Simulation of SH-SVL and SVL in lossless as well as lossy channels.
The lossy channels are modeled with an independent 30% packet loss. Error is
the maximum error.

##  7 Summary and Future Work

In this paper, we extended the analysis of our of our self-healing algorithms
to the case where packets are lost by constructing linear matrix inequalities
for both the case of synchronous packet loss and general edgewise packet loss.
We performed a numerical study which certified that SH-SVL will converge even
in the case of extremely high packet loss on networks with poor connectivity.

In the future, we plan to extend our synthesis and analysis techniques to the
more general cases of asynchronous updates, nonsmooth objectives, and
nonconvex objectives.

## References

  * [1] Pedro A Forero, Alfonso Cano and Georgios B Giannakis  “Consensus-Based Distributed Support Vector Machines.”  In _Journal of Machine Learning Research_ 11.5, 2010
  * [2] Gerasimos G. Rigatos  “Distributed gradient and particle swarm optimization for multi-robot motion planning” Copyright - Cambridge University Press; Last updated - 2015-08-15  In _Robotica_ 26.3, 2008, pp. 357–370
  * [3] A. Sundararajan, B. Van Scoy and L. Lessard  “Analysis and Design of First-Order Distributed Optimization Algorithms Over Time-Varying Graphs”  In _IEEE Transactions on Control of Network Systems_ 7.4, 2020, pp. 1597–1608  DOI: 10.1109/TCNS.2020.2988009
  * [4] Wei Shi, Qing Ling, Gang Wu and Wotao Yin  “EXTRA: An Exact First-Order Algorithm for Decentralized Consensus Optimization”  In _SIAM Journal on Optimization_ 25.2, 2015, pp. 944–966  DOI: 10.1137/14096668X
  * [5] Z. Li, W. Shi and M. Yan  “A Decentralized Proximal-Gradient Method With Network Independent Step-Sizes and Separated Convergence Rates”  In _IEEE Transactions on Signal Processing_ 67.17, 2019, pp. 4494–4506  DOI: 10.1109/TSP.2019.2926022
  * [6] D. Jakovetić  “A Unification and Generalization of Exact Distributed First-Order Methods”  In _IEEE Transactions on Signal and Information Processing over Networks_ 5.1, 2019, pp. 31–46  DOI: 10.1109/TSIPN.2018.2846183
  * [7] Angelia Nedić, Alex Olshevsky and Wei Shi  “Achieving Geometric Convergence for Distributed Optimization Over Time-Varying Graphs”  In _SIAM Journal on Optimization_ 27.4, 2017, pp. 2597–2633  DOI: 10.1137/16M1084316
  * [8] G. Qu and N. Li  “Harnessing Smoothness to Accelerate Distributed Optimization”  In _IEEE Transactions on Control of Network Systems_ 5.3, 2018, pp. 1245–1260  DOI: 10.1109/TCNS.2017.2698261
  * [9] K. Yuan, B. Ying, X. Zhao and A.. Sayed  “Exact Diffusion for Distributed Optimization and Learning—Part I: Algorithm Development”  In _IEEE Transactions on Signal Processing_ 67.3, 2019, pp. 708–723  DOI: 10.1109/TSP.2018.2875898
  * [10] K. Yuan, B. Ying, X. Zhao and A.. Sayed  “Exact Diffusion for Distributed Optimization and Learning—Part II: Convergence Analysis”  In _IEEE Transactions on Signal Processing_ 67.3, 2019, pp. 724–739  DOI: 10.1109/TSP.2018.2875883
  * [11] Laurent Lessard, Benjamin Recht and Andrew Packard  “Analysis and Design of Optimization Algorithms via Integral Quadratic Constraints”  In _SIAM J. Optim._ 26.1, 2016, pp. 57–95
  * [12] A. Sundararajan, B. Hu and L. Lessard  “Robust convergence analysis of distributed optimization algorithms”  In _2017 55th Annual Allerton Conference on Communication, Control, and Computing (Allerton)_ , 2017, pp. 1206–1212  DOI: 10.1109/ALLERTON.2017.8262874
  * [13] Bin Hu, Peter Seiler and Anders Rantzer  “A unified analysis of stochastic optimization methods using jump system theory and quadratic constraints”  In _Conference on Learning Theory_ , 2017, pp. 1157–1189  PMLR
  * [14] Pete Seiler and Raja Sengupta  “An h/sub/spl infin//approach to networked control”  In _IEEE Transactions on Automatic control_ 50.3 IEEE, 2005, pp. 356–364
  * [15] Ji Ma, Xiao Yu and Weiyao Lan  “Distributed Optimization Design of Multi-Agent Systems With Nonidentical Packet Losses”  In _2021 40th Chinese Control Conference (CCC)_ , 2021, pp. 5604–5609  IEEE
  * [16] I.. Donato Ridgley, R.. Freeman and K.. Lynch  “Private and Hot-Pluggable Distributed Averaging”  In _IEEE Control Systems Letters_ 4.4, 2020, pp. 988–993  DOI: 10.1109/LCSYS.2020.2996957
  * [17] Solmaz S. Kia et al.  “Tutorial on Dynamic Average Consensus: The Problem, its Applications, and the Algorithms”  In _IEEE Control Systems Magazine_ 39.3, 2019, pp. 40–72
  * [18] Christoforos N. Hadjicostis, Nitin H. Vaidya and Alejandro D. Domínguez-García  “Robust Distributed Average Consensus via Exchange of Running Sums”  In _IEEE Transactions on Automatic Control_ 61.6, 2016, pp. 1492–1507  DOI: 10.1109/TAC.2015.2471695
  * [19] Israel L. Donato Ridgley, Randy A. Freeman and Kevin M. Lynch  “Self-Healing First-Order Distributed Optimization”  In _2021 60th IEEE Conference on Decision and Control (CDC)_ , 2021, pp. 3850–3856  DOI: 10.1109/CDC45484.2021.9683487
  * [20] Christoforos N. Hadjicostis, Alejandro D. Domínguez-García and Themistoklis Charalambous  “Distributed Averaging and Balancing in Network Systems” 13, Foundations and Trends (R) in Systems and Control  Now Publishers, 2018
  * [21] A. Sundararajan, B. Van Scoy and L. Lessard  “A Canonical Form for First-Order Distributed Optimization Algorithms”  In _2019 American Control Conference (ACC)_ , 2019, pp. 4075–4080  DOI: 10.23919/ACC.2019.8814838
  * [22] Oswaldo Luiz Valle Costa, Marcelo Dutra Fragoso and Ricardo Paulino Marques  “Discrete-time Markov jump linear systems”  Springer Science & Business Media, 2005
  * [23] J Micah Fry, Mazen Farhood and Peter Seiler  “IQC-based robustness analysis of discrete-time linear time-varying systems”  In _International Journal of Robust and Nonlinear Control_ 27.16 Wiley Online Library, 2017, pp. 3135–3157
  * [24] MOSEK ApS  “MOSEK Optimizer API for C 8.0.0.81”, 2017  URL: http://docs.mosek.com/8.0/capi/index.html
  * [25] Patrick Kofod Mogensen and Asbjørn Nilsen Riseth  “Optim: A mathematical optimization package for Julia”  In _Journal of Open Source Software_ 3.24, 2018, pp. 615  DOI: 10.21105/joss.00615
  * [26] Michael Garstka, Mark Cannon and Paul Goulart  “COSMO: A conic operator splitting method for large convex problems”  In _European Control Conference_ , 2019  DOI: 10.23919/ECC.2019.8796161

[]Israel L. Donato RidgleyIsrael Donato Ridgley is a Ph.D. candidate in
Electrical Engineering at Northwestern University (Evanston, Illinois)
expected to graduate in Summer 2023. He is an affiliate of the Center for
Robotics and Biosystems. His research interests include distributed control
and optimization, mulit-agent systems, and privacy in multi-agent systems.

[]Randy A. FreemanRandy Freeman received a Ph.D. in Electrical Engineering
from the University of California at Santa Barbara in 1995. Since then he has
been a faculty member at Northwestern University (Evanston, Illinois), where
he is currently Professor of Electrical and Computer Engineering. He is a
member of the Northwestern Institute on Complex Systems and an affiliate of
the Center for Robotics and Biosystems. His research interests include
nonlinear systems and control, distributed control and optimization, multi-
agent systems, robust control, and optimal control.

[]Kevin M. Lynch(S’90–M’96–SM’05–F’10) received the B.S.E. degree in
electrical engineering from Princeton University, Princeton, NJ, USA, and the
Ph.D. degree in robotics from Carnegie Mellon University, Pittsburgh, PA, USA.

He is a professor of mechanical engineering at Northwestern University, where
he directs the Center for Robotics and Biosystems and is a member of the
Northwestern Institute on Complex Systems. He is a coauthor of the textbooks
_Principles of Robot Motion_ (Cambridge, MA, USA: MIT Press, 2005) and _Modern
Robotics: Mechanics, Planning, and Control_ (Cambridge, U.K.: Cambridge Univ.
Press, 2017) and the associated online courses and videos. His research
interests include robot manipulation and locomotion, self-organizing
multiagent systems, and physically collaborative human–robot systems.