#  Discretized Distributed Optimization over Dynamic Digraphs

Mohammadreza Doostmohammadian, Wei Jiang, Muwahida Liaquat, Alireza Aghasi,
Houman Zarrabi M. Doostmohammadian is with the Faculty of Mechanical
Engineering at Semnan University, Semnan, Iran (doost@semnan.ac.ir). W. Jiang
and M. Liaquat are with the School of Electrical Engineering at Aalto
University, Espoo, Finland (firstname.lastname@aalto.fi). M. Liaquat is also
with Finnish Geospatial Research Institute (FGI). A. Aghasi is with the
Department of Electrical Engineering and Computer Science, Oregon State
University, USA (alireza.aghasi@oregonstate.edu). H. Zarrabi is with Iran
Telecom Research Centre (ITRC), Tehran, Iran (h.zarrabi@itrc.ac.ir).

###### Abstract

We consider a discrete-time model of continuous-time distributed optimization
over dynamic directed-graphs (digraphs) with applications to distributed
learning. Our optimization algorithm works over general strongly connected
dynamic networks under switching topologies, e.g., in mobile multi-agent
systems and volatile networks due to link failures. Compared to many existing
lines of work, there is no need for bi-stochastic weight designs on the links.
The existing literature mostly needs the link weights to be stochastic using
specific weight-design algorithms needed both at the initialization and at all
times when the topology of the network changes. This paper eliminates the need
for such algorithms and paves the way for distributed optimization over time-
varying digraphs. We derive the bound on the gradient-tracking step-size and
discrete time-step for convergence and prove dynamic stability using arguments
from consensus algorithms, matrix perturbation theory, and Lyapunov theory.
This work, particularly, is an improvement over existing stochastic-weight
undirected networks in case of link removal or packet drops. This is because
the existing literature may need to rerun time-consuming and computationally
complex algorithms for stochastic design, while the proposed strategy works as
long as the underlying network is weight-symmetric and balanced. The proposed
optimization framework finds applications to distributed classification and
learning.

###### Note to Practitioners

Inspired by recent advances in cloud-computing and distributed and parallel
processing along with embedded low-cost CPUs and wireless communications, this
paper considers distributed algorithms for optimization and machine learning
over wireless-connected autonomous multi-agent systems (MASs). In contrast to
the classical centralized learning methods, which are prone to single-point-
of-failure and centralized processing, in cooperative optimization the
learning is distributed among a group of data-processing agents (e.g. robots)
with communication units. This article provides an efficient algorithm to
enable MASs to collaboratively optimize a cost function, e.g., for binary
classification and distributed support vector machine (D-SVM). Sampled data
systems related to MASs and robotic networks, due to digital communications
and discretized control models, use discrete-time algorithms that need to
account for intermittent communications and dynamic networking. This requires
discretized algorithms over dynamic digraphs, practical in the presence of
packet drops (or lost communication channels). Most existing distributed
algorithms either are susceptible to change in the communication network or
impose computationally inefficient stochastic design. These algorithms need to
be rerun, for example, whenever the mobile robotic network changes due to
limited communication range. This makes most existing algorithms infeasible in
real-time applications. Our proposed method outperforms similar algorithms
over dynamic robotic networks and in the presence of link removal (packet
drops). We show this efficiency of our distributed optimization method by
simulation.

###### Index Terms:

distributed optimization, matrix perturbation theory, consensus constraint,
dynamic digraphs

##  I Introduction

Distributed optimization finds application in many emerging areas from
estimation [1] and machine learning [2] to distributed scheduling and resource
allocation [3, 4, 5, 6, 7]. In many applications, the network topology might
be time-varying [8, 2, 9] or contain uni-directional transmission links and
data-exchange [10, 2, 9, 11, 12]. This motivates us to explore optimization
strategies over general dynamic directed graphs (digraphs).

This work extends the existing works in the sense that the link weights are
not necessarily bi-stochastic or split into row and column stochastic matrices
[9, 11, 12]. This relaxes the need for weight-stochastic design algorithms
[13] for weight-symmetric networks in case of link failure or mutual data
packet drops. Recall that, any link failure or change in the network topology
fails the stochastic condition, and thus, compensation algorithms [13, 14, 15,
16, 17] are needed to update the link weights to comply with stochasticity in
the existing literature. Eliminating these weight constraints allows for the
proposed solution to be easily applied over time-varying digraph topologies
and networks subject to packet drops. To prove the convergence of our
continuous-time dynamics, we first determine the admissible range of gradient-
tracking step-size based on the perturbation analysis [18]. Proving that all
the eigenvalues have negative real values (except one zero eigenvalue for each
state dimension), the stability is proved; however, the explicit bound on the
convergence rate of the algorithm is left for future work111Our distributed
optimization algorithm involves multiple steps like local updates, consensus,
communication scheduling, auxiliary variable updates, and more. The explicit
interactions between these steps can be complicated to analyze using
perturbation theory. Therefore, due to the complexity of perturbation-based
analysis and due to the dynamic nature of the network topology it is hard to
explicitly find the bound on the convergence rate. . Next, by discretizing the
dynamics, we bound the discrete-time step-size for stability. As compared to
some other existing literature, our algorithm operates in a single time-scale
with no need for an extra inner consensus loop (as an additional time-scale)
[19, 20], and with no need for stochastic weight design [9, 11, 12]. This
follows the fact that our weight-balanced condition is more relaxed than the
stochastic condition in general. We give particular examples to show
improvement over weight-symmetric undirected unreliable networks subject to
link removals. We give illustrative examples and simulations to support this
claim.

General Notations: $\partial_{a}(z)=\frac{dz}{da}$ denotes the derivative with
respect to $a$. $\sigma(A)$ represents the eigen-spectrum of the matrix $A$.
$I_{n}$ and $\mathbf{0}_{n\times n}$ denote the identity matrix and all-zero
matrix of size $n$. $\mathbf{0}_{n}$ and $\mathbf{1}_{n}$ denote the column
vector of all ones or zeros. $\nabla F$ and $\nabla^{2}F$ denote the gradient
and second gradient (Hessian) of function $F$. ${\text{blockdiag}[A_{i}]}$
implies a block-diagonal matrix with each $A_{i}$ as its $i$th diagonal block.
$\otimes$ denotes the Kronecker product. LHP and RHP are abbreviations of
left-half-plane and right-half-plane. The operator norm of a matrix is defined
as [18]: $\lVert A\rVert=\sup_{\lVert\mathbf{x}\rVert=1}\lVert
A\mathbf{x}\rVert$.

##  II Problem Statement

###  II-A Preliminaries

We represent the information exchange between nodes by a graph
$\mathcal{G}_{q}=\\{\mathcal{V},\mathcal{E}_{q}\\}$ of order/size $n$, with
(time-invariant) set of nodes $\mathcal{V}$ and
$\mathcal{E}_{q}\subseteq\mathcal{V}\times\mathcal{V}$ as the (time-varying)
set of links ($q$ as the switching signal). A directed graph (digraph) is
_strongly connected (SC)_ if there exists a directed path between every two
nodes. The diameter $d_{g}$ of $\mathcal{G}$ is the longest shortest path
between any two nodes. Matrices ${W=\\{w_{ij}\\}}$, ${A=\\{a_{ij}\\}}$
represent two (generally different) adjacency matrices of $\mathcal{G}$. The
Laplacian matrix $\overline{W}=\\{\overline{w}_{ij}^{q}\\}$, is defined as
$\overline{w}_{ij}^{q}=w_{ij}$ for $i\neq j$ and
$\overline{w}_{ii}^{q}=-\sum_{j=1}^{n}w_{ij}$ for $i=j$ (similar definition
holds for $\overline{A}$). This implies that $\overline{W}\mathbf{1}_{n}=0$
and $\overline{A}\mathbf{1}_{n}=0$.

The work [21] equivalently defines the Laplacian as $I_{n}-\overline{D}^{-1}W$
(known as the multi-rate integrator) with the entries of the diagonal
(modified) degree matrix $\overline{D}$ as
$\overline{D}_{i,i}=\sum_{j=1}^{n}w_{ij}+c$ and $c$ as an additive fixed
constant to avoid the singularity. A network is weight-balanced if
$\overline{W}\mathbf{1}_{n}=\overline{W}^{\top}\mathbf{1}_{n}$ (and similarly
$\overline{A}\mathbf{1}_{n}=\overline{A}^{\top}\mathbf{1}_{n}$). It is column
(resp. row) stochastic if $\overline{W}^{\top}\mathbf{1}_{n}=\mathbf{1}_{n}$
(resp. $\overline{W}\mathbf{1}_{n}=\mathbf{1}_{n}$) and bi-stochastic if both
row and column stochastic. An undirected network with symmetric $W$ (and $A$)
is weight-symmetric and weight-balanced. A weight-symmetric and row (or
column) stochastic network is also bi-stochastic.

###  II-B Formulation

The distributed consensus-constraint optimization problem is as follows:

$\displaystyle\text{minimize}_{\mathbf{x}\in\mathbb{R}^{mn}}F(\mathbf{x}),\qquad
F(\mathbf{x})=\sum_{i=1}^{n}f_{i}(\mathbf{x}_{i})$ $\displaystyle\text{subject
to}~{}\mathbf{x}_{1}=\mathbf{x}_{2}=\dots=\mathbf{x}_{n}.$ (1)

Problem (1) can be reformulated as the Laplacian-constraint formulation [10],

$\displaystyle\text{minimize}_{\mathbf{x}\in\mathbb{R}^{nm}}F(\mathbf{x})=\sum_{i=1}^{n}f_{i}(\mathbf{x}_{i})~{}~{}\text{subject
to}~{}~{}\mathcal{L}\mathbf{x}=\mathbf{0}_{nm}.$ (2)

where $\mathcal{L}:=L\otimes I_{m}$ and $L$ is the Laplacian of the weight-
balanced digraph $\mathcal{G}$. Recall from consensus literature [22] that
state dynamics in the form $\dot{\mathbf{x}}=\mathcal{L}\mathbf{x}$ over
strongly-connected networks results in state consensus
$\mathbf{x}_{1}=\mathbf{x}_{2}=\dots=\mathbf{x}_{n}$. This is because, for
strongly-connected network, the equilibrium
$\mathcal{L}\mathbf{x}=\mathbf{0}_{nm}$ holds if and only if
$\mathbf{x}_{1}=\mathbf{x}_{2}=\dots=\mathbf{x}_{n}$ (i.e., consensus is
achieved). This gives the intuition why the two formulations are the same.
Therefore, for the transition from problem (1) to problem (2), the role or
condition of the network and edges is that they must satisfy strong-
connectivity. See more details in [10, Section III] Then, under certain
conditions (see [10, 20] for details), problem (1) is shown to be equivalent
with general unconstrained formulation222Let the optimal state of formlation
(3) be $\mathbf{x}^{*}$, then the optimal state of formulations (1) and (2) is
in the form $\mathbf{1}_{n}\otimes\mathbf{x}^{*}$.  as,

$\displaystyle\text{minimize}_{\mathbf{x}\in\mathbb{R}^{m}}$ $\displaystyle
F(\mathbf{x})=\sum_{i=1}^{n}f_{i}(\mathbf{x})$ (3)

In this work we solve problem (1) and the solution can be easily extended to
problems (2) and (3) (similar to [10, 20]).

######  Assumption 1.

The local cost functions $f_{i}:\mathbb{R}^{m}\mapsto\mathbb{R}$ are smooth,
and strictly convex with locally $\gamma$-Lipschitz gradient.

This assumption implies that each local cost
$f_{i}:\mathbb{R}^{m}\mapsto\mathbb{R}$ is (locally) twice differentiable and
strictly convex. Then, the optimal state of (3) is in the form
$\mathbf{x}^{*}:=\mathbf{1}_{n}\otimes\overline{\mathbf{x}}^{*}$ with
$\overline{\mathbf{x}}^{*}\in\mathbb{R}^{m}$ and
${(\mathbf{1}_{n}^{\top}\otimes
I_{m})\boldsymbol{\nabla}F(\mathbf{x}^{*})=\mathbf{0}_{m}}$ [10]. For
application in D-SVM [2], we assume that
$f_{i}(\mathbf{x}_{i})=\sum_{j=1}^{m}f_{i,j}(x_{i,j})$ with $x_{i,j}$ denoting
the $j$th element of $\mathbf{x}_{i}$.

##  III The Proposed Continuous-Time Model

###  III-A The Networked Dynamics

To solve problem (1), the following continuous-time linear dynamics is
proposed:

$\displaystyle\dot{\mathbf{x}}_{i}$
$\displaystyle=-\sum_{j=1}^{n}w^{q}_{ij}(\mathbf{x}_{i}-\mathbf{x}_{j})-\alpha\mathbf{y}_{i},$
(4) $\displaystyle\dot{\mathbf{y}}_{i}$
$\displaystyle=-\sum_{j=1}^{n}a^{q}_{ij}(\mathbf{y}_{i}-\mathbf{y}_{j})+\partial_{t}\boldsymbol{\nabla}f_{i}(\mathbf{x}_{i}),$
(5)

with $n$ as the number of nodes/agents over the network, vector
${\mathbf{x}_{i}\in\mathbb{R}^{m}}$ as the node $i$’s state, vector
${\mathbf{y}_{i}\in\mathbb{R}^{m}}$ as an auxiliary variable for gradient-
tracking (GT), ${\dot{\mathbf{x}}_{i}=\partial_{t}\mathbf{x}}_{i}$, and
${\alpha\in\mathbb{R}_{>0}}$ as the GT step-size. Note that the general time
index is denoted by $t$ and the network switching signal is denoted by index
$q$. This switching signal $q:t\mapsto Q$ is a function of time, taking
integer values in the set $Q$; in fact, $q(t)$ tells us what configuration the
graph topology and link weights have at time $t$. $q(t)$ can be periodic or
the links can randomly disappear or reappear with a fixed probability. For
example, consider $10$ strongly-connected graph topologies
$\\{\mathcal{G}_{0},\mathcal{G}_{2},\dots,\mathcal{G}_{9}\\}$; the switching
signal could be $q=\lfloor\frac{t}{10}\rfloor$ (where $\lfloor x\rfloor$
denotes greatest integer less than or equal to $x$); this switching signal
gives the index of the associated graph in the set $Q=\\{0,1,\dots,9\\}$ at
time $t$. The only constraint on these network topologies is that they must be
strongly-connected at all time-instants. Define the Hessian matrix
${H\coloneqq\text{blockdiag}[\boldsymbol{\nabla}^{2}f_{i}(\mathbf{x}_{i})]}$,
then, Eq. (4)-(5) is written in a compact form

$\displaystyle\left(\begin{array}[]{c}\dot{\mathbf{x}}\\\
\dot{\mathbf{y}}\end{array}\right)=M(t,\alpha)\left(\begin{array}[]{c}{\mathbf{x}}\\\
{\mathbf{y}}\end{array}\right),$ (10) $\displaystyle M(t,\alpha)=$
$\displaystyle\left(\begin{array}[]{cc}\overline{W}\otimes I_{m}&-\alpha
I_{mn}\\\ H(\overline{W}\otimes I_{m})&\overline{A}\otimes I_{m}-\alpha
H\end{array}\right).$ (13)

We make the following assumption on $\mathcal{G}$, $W$ and $A$.

######  Assumption 2.

The graph $\mathcal{G}$ is directed, and strongly connected at every time $t$.
The link weights are positive and $\sum_{j=1}^{n}w^{q}_{ij}<1$,
$\sum_{j=1}^{n}a^{q}_{ij}<1$. The matrices ${W}$ and ${A}$ are weight-
balanced.

The assumption that $\sum_{j=1}^{n}w^{q}_{ij}<1$, $\sum_{j=1}^{n}a^{q}_{ij}<1$
is for the sake of proof analysis on bounding $\alpha$ (see Sections III-B and
III-C). In case of having
$\max\\{\sum_{j=1}^{n}w^{q}_{ij},\sum_{j=1}^{n}a^{q}_{ij}\\}=c>1$ the $\alpha$
parameter needs to be reduced by a factor of $c$ accordingly. Also, the
matrices ${W}$ and ${A}$ can be equal (as a special case). The definition of
matrices ${W}$ and ${A}$ in the Assumption 2 are relaxed as compared to
stochastic features in many existing literature. This relaxes the stochastic
weight design in case of link failure or switching topologies.

The difference between bi-stochastic weight design and balanced weight design
is explained here. Recall that bi-stochastic link weights imply that
$\sum_{j=1}^{n}w^{q}_{ij}=\sum_{i=1}^{n}w^{q}_{ij}=1$,
$\sum_{j=1}^{n}a^{q}_{ij}=\sum_{i=1}^{n}a^{q}_{ij}=1$. On the other hand, the
weight-balanced design is more relaxed as
$\sum_{j=1}^{n}w^{q}_{ij}=\sum_{i=1}^{n}w^{q}_{ij}$,
$\sum_{j=1}^{n}a^{q}_{ij}=\sum_{i=1}^{n}a^{q}_{ij}$ (there is no need the sum
to equal to $1$). This relaxation significantly reduces the complexity of the
algorithm, particularly for time-varying network topologies. For example
consider the case that a bidirectional symmetric link between nodes $a$, $b$
is removed. Thus, the link weights $w_{ab}=w_{ba}$ are removed while the sum
of the weights still remains balanced, i.e.,
$\sum_{j=1}^{n}w^{q}_{ij}-w_{ab}=\sum_{i=1}^{n}w^{q}_{ij}-w_{ba}$. However,
since the weight-stochasticity does not hold anymore
$\sum_{j=1}^{n}w^{q}_{ij}-w_{ab}=\sum_{i=1}^{n}w^{q}_{ij}-w_{ba}\neq 1$, we
need to redesign the rest of the link weights to satisfy bi-stochasticity.

######  Remark 1.

The existing optimization literature mostly requires bi-stochastic weigh
design [9, 11, 12]. One drawback of such a requirement is in the case of link
removal when the network loses the bi-stochastic condition. This implies that
these existing algorithms do not converge in the case of link removal.
Therefore, the weight compensation design algorithms, e.g. [13, 14], are
needed to redesign the weights to make the network bi-stochastic again for
convergence. This adds more complexity to the solution whenever the network
topology changes. Other than the extra complexity, it is likely that the
stochastic weight design cannot be even satisfied. Note that not all network
topologies have the necessary condition for bi-stochastic weight design. More
discussions and illustrations on this are given in Section IV.

######  Lemma 1.

[22, Theorem 3] Given SC digraph $\mathcal{G}$ of size $n$, the real part of
eigenvalues of ${\overline{W}}$ is non-positive, with the algebraic
multiplicity of zero eigenvalue equal to $1$. For its irreducible $0$-$1$
structured Laplacian matrix ${\overline{W}}$, consider the right and left
eigenvectors $\mathbf{v}_{1}$ and $\mathbf{u}_{1}^{\top}$ (including the ones
associated with the $0$ eigenvalue) satisfying
$\mathbf{u}_{1}^{\top}\mathbf{v}_{1}=1$. Then, the matrix
$R=\mathbf{v}_{1}^{\top}\mathbf{u}_{1}$ is real-valued. Further, for the $0$
eigenvalue $\mathbf{v}_{1}=\frac{1}{\sqrt{n}}\mathbf{1}_{n}$.

######  Lemma 2.

For the directed graph $\mathcal{G}$ described in Assumption 2 and its
associated Laplacian $\overline{W}$, the left eigenvector
$\mathbf{u}_{1}^{\top}$ of the $0$ eigenvalue is non-negative. For such a
matrix one can assume the left eigenvector satisfies
$\sum_{i=1}^{n}u_{1,i}>0$.

The above is common in weighted-consensus literature to reach a weighted
average of the states, see [22, Corollary 2]. Recalling arguments from
consensus literature, from Assumption 2 and the disagreement-decay structure
of (4)-(5), we obtain (asymptotically):

$\displaystyle\sum_{i=1}^{n}\dot{\mathbf{y}}_{i}$
$\displaystyle=\sum_{i=1}^{n}\partial_{t}\boldsymbol{\nabla}f_{i}(\mathbf{x}_{i}),$
(14) $\displaystyle\sum_{i=1}^{n}\dot{\mathbf{x}}_{i}$
$\displaystyle=-\alpha\sum_{i=1}^{n}\mathbf{y}_{i}.$ (15)

The above directly follows from the consensus-type structure of the first
terms in Eq. (4)-(5). By setting initial values of
$\mathbf{y}(0)=\mathbf{0}_{nm}$ and simple integration, we get that
$\sum_{i=1}^{n}\dot{\mathbf{x}}_{i}$ tracks
$-\sum_{i=1}^{n}\boldsymbol{\nabla}f_{i}(\mathbf{x}_{i})$ as the gradient sum
and eventually reaches zero as discussed next. By setting
${\dot{\mathbf{x}}_{i}=0_{m}}$, the equilibrium $\underline{\mathbf{x}}^{*}$
of the dynamics (4)-(5) satisfies ${(\mathbf{1}_{n}^{\top}\otimes
I_{m})\boldsymbol{\nabla}F(\mathbf{x}^{*})=\mathbf{0}_{m}}$. This holds for
the optimal state of (1) as well [10].

######  Lemma 3.

Let $\mathbf{y}(0)=\mathbf{0}_{nm}$ and
$\mathbf{x}(0)\neq\mathbf{1}_{n}\otimes{\mathbf{x}}_{0}$, for some non-zero
$\mathbf{x}_{0}\in\mathbb{R}^{m}$. The global equilibrium of (4)-(5) includes
the states in the form
$[\mathbf{1}_{n}\otimes\overline{\mathbf{x}}^{*};\mathbf{0}_{nm}]$ (as the
invariant set) with $\overline{\mathbf{x}}^{*}\in\mathbb{R}^{m}$ and includes
the optimizer $\mathbf{x}^{*}$ of problem (1) which satisfies
$(\mathbf{1}_{n}^{\top}\otimes
I_{m})\boldsymbol{\nabla}F(\mathbf{x}^{*})=\mathbf{0}_{m}$.

###### Proof.

From the consensus-type structure of (4)-(5), the following holds for any
${\mathbf{x}=\mathbf{1}_{n}\otimes\overline{\mathbf{x}}^{*}},\overline{\mathbf{x}}^{*}\in\mathbb{R}^{m}$;
$\dot{\mathbf{x}}_{i}=\mathbf{0}_{m}$ from (4), and
$\dot{\mathbf{y}}_{i}=\partial_{t}\boldsymbol{\nabla}f_{i}(\overline{\mathbf{x}}^{*})=\boldsymbol{\nabla}^{2}f_{i}(\overline{\mathbf{x}}^{*})\dot{\mathbf{x}}_{i}=\mathbf{0}_{m}$
from (5). From (14)-(15), the following uniquely holds for any
$\mathbf{x}=\mathbf{x}^{*}$ (which is also in the form
$\mathbf{1}_{n}\otimes\overline{\mathbf{x}}^{*}$),

$\sum_{i=1}^{n}\dot{\mathbf{x}}_{i}=-\alpha(\mathbf{1}_{n}^{\top}\otimes
I_{m})\boldsymbol{\nabla}F(\mathbf{x}^{*})=\mathbf{0}_{m}.$

This shows that $[\mathbf{x}^{*};\mathbf{0}_{nm}]$ is an invariant equilibrium
under dynamics (4)-(5) and proves the lemma. ∎

Lemma 3 defines the invariant set of the networked dynamics (4)-(5) in the
form $[\mathbf{x}^{*};\mathbf{0}_{nm}]$. Next, we need to show that the
dynamics uniquely converges to this $\mathbf{x}^{*}$ as the optimal point of
(1).

###  III-B Proof of Convergence

######  Lemma 4.

[23, 24] Consider the $n$-by-$n$ matrix $P(\alpha)$ as a smooth function of
${\alpha\in\mathbb{R}_{\geq 0}}$. Let $P(0)$ has ${l<n}$ equal eigenvalues
$\lambda_{1}=\ldots=\lambda_{l}$, associated with linearly independent right
and left unit eigenvectors $\mathbf{v}_{1},\ldots,\mathbf{v}_{l}$ and
$\mathbf{u}_{1},\ldots,\mathbf{u}_{l}$. Let
${P^{\prime}=\partial_{\alpha}P(\alpha)|_{\alpha=0}}$ and
$\lambda_{i}(\alpha)$ denote the $i$-th eigenvalue as a function of $\alpha$,
which corresponds to ${\lambda_{i},i\in\\{1,\ldots,l\\}}$. Then,
$\partial_{\alpha}\lambda_{i}|_{\alpha=0}$ is the $i$-th eigenvalue of,

$\left(\begin{array}[]{ccc}\mathbf{u}_{1}^{\top}P^{\prime}\mathbf{v}_{1}&\ldots&\mathbf{u}_{1}^{\top}P^{\prime}\mathbf{v}_{l}\\\
&\ddots&\\\
\mathbf{u}_{l}^{\top}P^{\prime}\mathbf{v}_{1}&\ldots&\mathbf{u}_{l}^{\top}P^{\prime}\mathbf{v}_{l}\end{array}\right).$

######  Theorem 1.

Let Assumptions 2 and conditions in Lemma 2 hold. For sufficiently small
$\alpha$, all eigenvalues of $M$ have non-positive real-parts $\forall t>0$,
and the algebraic multiplicity of the zero eigenvalue is $m$ (as the
individual state dimension).

###### Proof.

From (13), let $M(\alpha)=M_{0}+\alpha M_{1}$ with

$\displaystyle M_{0}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cc}\overline{W}\otimes
I_{m}&\mathbf{0}_{mn\times mn}\\\ H(\overline{W}\otimes
I_{m})&\overline{A}\otimes I_{m}\end{array}\right),$ (18) $\displaystyle
M_{1}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cc}\mathbf{0}_{mn\times mn}&-{I_{mn}}\\\
{\mathbf{0}_{mn\times mn}}&-H\end{array}\right),$ (21)

Since $M_{0}$ is block triangular, its eigen-spectrum $\sigma(\cdot)$ is

$\displaystyle\sigma(M_{0})=\sigma(\overline{W}\otimes
I_{m})\cup\sigma(\overline{A}\otimes I_{m})$ (22)

From Lemma 1, both matrices $\overline{W}$ and $\overline{A}$ have one
isolated zero eigenvalue and the rest in the LHP. The $m$ sets of eigenvalues
of the matrix $M_{0}$ (of size $2mn$), associated with dimensions
$j=\\{1,\ldots,m\\}$ of states $\mathbf{x}_{i}$, then follow as,

$\operatorname{Re}\\{\lambda_{2n,j}\\}\leq\ldots\leq\operatorname{Re}\\{\lambda_{3,j}\\}<\lambda_{2,j}=\lambda_{1,j}=0,$

Next, for spectrum analysis, we consider the term $\alpha M_{1}$ as a
perturbation to $M_{0}$. Using Lemma 4, one can check how the spectrum of $M$
changes, as $M_{0}$ is perturbed by the values of $\alpha$ and $M_{1}$. In
particular, we check the variation of the zero eigenvalues $\lambda_{1,j}$ and
$\lambda_{2,j}$ of $M_{0}$ by the perturbation term $\alpha M_{1}$. Let
$\lambda_{1,j}(\alpha)$ and $\lambda_{2,j}(\alpha)$ denote the perturbed
eigenvalues associated to $M$. The right and left unit eigenvectors of
$\lambda_{1,j}$, $\lambda_{2,j}$ follow from Lemma 1 as333The normalizing
factors of the unit vectors might be ignored as in the followings we only care
about the sign of the terms in $U^{\top}M_{1}V$, not the exact values. ,

$\displaystyle V$
$\displaystyle=[V_{1}~{}V_{2}]=\frac{1}{\sqrt{n}}\left(\begin{array}[]{cc}\mathbf{1}_{n}&\mathbf{0}_{n}\\\
\mathbf{0}_{n}&\mathbf{1}_{n}\end{array}\right)\otimes I_{m}$ (25)
$\displaystyle U^{\top}$
$\displaystyle=[U_{1}~{}U_{2}]^{\top}=\left(\begin{array}[]{cc}\mathbf{u}_{1}&\mathbf{0}_{n}\\\
\mathbf{0}_{n}&\mathbf{u}_{2}\end{array}\right)^{\top}\otimes I_{m}$ (28)

These unit eigenvectors satisfy $V^{\top}V=I_{2mn}$ and $U^{\top}U=I_{2mn}$.
From the definition $\partial_{\alpha}{dM(\alpha)}|_{\alpha=0}=M_{1}$ and
Lemma 4,

$\displaystyle U^{\top}M_{1}V=\left(\begin{array}[]{cc}\mathbf{0}_{m\times
m}&\times\\\ \mathbf{0}_{m\times m}&-(\mathbf{u}_{2}\otimes
I_{m})^{\top}H(\frac{1}{\sqrt{n}}\mathbf{1}_{n}\otimes
I_{m})\end{array}\right).$ (31)

From Lemma 1, let $u_{2,i}$ denote the $i$th element of $\mathbf{u}_{2}$.

$\displaystyle-(\mathbf{u}_{2}\otimes
I_{m})^{\top}H(\frac{1}{\sqrt{n}}\mathbf{1}_{n}\otimes
I_{m})=-\sum_{i=1}^{n}\frac{u_{2,i}}{\sqrt{n}}\boldsymbol{\nabla}^{2}f_{i}(\mathbf{x}_{i})\prec
0,$ (32)

This follows from the definition of $H$, Lemma 2, and the strict convexity of
the cost function in Assumption 1. From Lemma 4,
$\partial_{\alpha}d\lambda_{1,j}|_{\alpha=0}$ and
$\partial_{\alpha}d\lambda_{2,j}|_{\alpha=0}$ follow the eigenvalues of the
triangular matrix in (31). This matrix has $m$ zero eigenvalues and, from
(32), $m$ negative eigenvalues and, thus,
$\partial_{\alpha}\lambda_{1,j}|_{\alpha=0}=0$ and
$\partial_{\alpha}\lambda_{2,j}|_{\alpha=0}<0$. This implies that $\alpha
M_{1}$ as a perturbation, pushes the $m$ zero eigenvalues
$\lambda_{2,j}(\alpha)$ of $M$ toward the LHP while the other $m$ zero
eigenvalues $\lambda_{1,j}(\alpha)$ remain unchanged. Thus, one can conclude
that for sufficiently small $\alpha$, the eigen-spectrum obeys

$\displaystyle\begin{aligned}
\operatorname{Re}\\{\lambda_{2n,j}(\alpha)\\}\leq\ldots\leq\operatorname{Re}\\{\lambda_{3,j}(\alpha)\\}\leq\lambda_{2,j}(\alpha)<\lambda_{1,j}(\alpha)=0,\end{aligned}$
(33)

which completes the proof. ∎

######  Remark 2.

Any pair of Laplacian matrices $\overline{A},\overline{W}$ whose left and
right unit eigenvectors satisfy Eq. (32) implies that the perturbed eigenvalue
moves towards LHP. The only constraint is to satisfy the consensus property on
the first terms (on the right-hand) of Eq. (4)-(5). In addition, for
homogeneous quadratic cost models, we have $H=\gamma I_{n}$ which is time-
invariant. Then, Lemma 2 can be more relaxed as Eq. (32) easily holds for any
$\mathbf{u}_{1}^{\top}$ satisfying $\sum_{i=1}^{n}u_{1,i}>0$.

The above theorem implies that, for an admissible range of $\alpha>0$, the
time-varying matrix $M$ has only one set of $m$ zero eigenvalues associated
with the eigenvectors $V_{1}$ in (25), and its null space
$\mathcal{N}(M)=\text{span}\\{[\mathbf{1}_{n};\mathbf{0}_{n}]\otimes I_{m}\\}$
is time-independent. Recall that the above proof uses the fact that the
eigenvalues are continuous functions of the matrix entries [25]. Note that
such matrix perturbation analysis as in Theorem 1 allows checking the eigen-
spectrum of hybrid systems as (10)-(13) with time-varying system matrices and
discrete jumps due to dynamic network topologies.

###  III-C Bounding $\alpha$

Next, we need to determine the upper bound on $\alpha$ such that no other
eigenvalues of $M_{0}$ under the perturbation moves to the RHP. Some analysis
is given in [2, Lemma 7] based on the notion of optimal matching distance
${d(\sigma(M),\sigma(M_{0}))=\min_{\pi}\max_{1\leq i\leq
2nm}(\lambda_{i}-\lambda_{\pi(i)}(\alpha))}$ with $\pi$ denoting all possible
permutations over $2nm$ symbols [18]. This parameter
$d(\sigma(M),\sigma(M_{0}))$ gives the smallest-radius circle centering at and
including $\lambda_{1,j},\ldots,\lambda_{2n,j}$, i.e., the farthest distance
between the eigenvalues of $M$ and $M_{0}$. The bound on $\alpha$ is defined
via the following lemma.

######  Lemma 5.

[18, Theorem 39.1] Given the system matrix ${M(\alpha)=M_{0}+\alpha M_{1}}$,
the optimal matching distance satisfies ${d(\sigma(M),\sigma(M_{0}))\leq
4(\lVert M_{0}\rVert+\lVert M\rVert)^{1-\frac{1}{nm}}\lVert\alpha
M_{1}\rVert^{\frac{1}{nm}}}$.

The sketch of the derivation is shown in Fig. 1.

Figure 1: The perturbation analysis to bound $d(\sigma(M),\sigma(M_{0}))$: the
min eigenvalue (in absolute value) of $M_{0}$ may move towards RHP if
$d(\sigma(M),\sigma(M_{0}))>\underline{\lambda}$.

In order to bound $d(\sigma(M),\sigma(M_{0}))$ by $\underline{\lambda}$ as the
minimum eigenvalue of $M_{0}$, assuming ${w_{ij},a_{ij}\geq 0}$ such that
$\sum_{j=1}^{n}w_{ij}<1$ and ${\sum_{j=1}^{n}a_{ij}<1}$ (similar to [2]), one
can find a bound on $\alpha$ for convergence as follows.

Replacing the spectral norm in Lemma 5 and after some simplification,

$\displaystyle
4(|\lambda_{2nm}|+|\lambda_{2nm}|+\underline{\lambda})^{1-\frac{1}{nm}}\alpha^{\frac{1}{nm}}\max\\{1,\gamma\\}^{\frac{1}{nm}}<\underline{\lambda},$

with $\lambda_{2nm}$ denoting the largest (in absolute value) eigenvalue of
$M_{0}$. This gives

$\displaystyle
0<{\alpha}\leq\frac{\underline{\lambda}^{nm}}{4^{nm}(2|\lambda_{2nm}|+\underline{\lambda})^{nm-1}\max\\{1,\gamma\\}}$
(34)

The next Lemma borrowed from [22], is originally known as the Gresgorin circle
Theorem [26] and gives a better perspective on the spectral localization of
the eigenvalues of $M_{0}$,

######  Lemma 6.

[22, Theorem 2] Given the Laplacian matrix $\overline{W}$ of a digraph
$\mathcal{G}$, its eigenvalues lie in the following disk (known as the
Gresgorin disk) in the complex plane:

$\displaystyle\mathcal{D}(\overline{w}_{\max},\overline{w}_{\max})=\\{z\in\mathbb{C}||z+\overline{w}_{\max}|\leq\overline{w}_{\max}\\}.$
(35)

where $\overline{w}_{\max}=\sum_{j=1,j\neq i}^{n}|\overline{w}_{ij}|$.
Following the definition of the Laplacian
$\overline{w}_{\max}=\max_{i}\\{|\overline{w}_{ii}|\\}$.

Recall that the eigenvalues of $M_{0}$ coincide with the eigenvalues of
$\overline{W}$ and $\overline{A}$, implying that in (35), the radius (and
centre) of the Gresgorin disk is
$\pm\max_{i}\\{|\overline{w}_{ii}|,|\overline{a}_{ii}|\\}$. This is used in
the next section to estimate the maximum (absolute) eigenvalue of $M(\alpha)$.

Next, we find another admissible $\alpha$ range similar to the analysis in
[11]. First, we recall from [27, Appendix] to relate the spectrum of
$M(t,\alpha)$ in (13) to $\alpha$. By proper row/column permutations in [27,
Eq. (18)], one can find $\sigma(M)$ as,

$\displaystyle\mbox{det}($ $\displaystyle\alpha
I_{mn})\mbox{det}(H(\overline{W}\otimes I_{m})+$
$\displaystyle(\overline{A}\otimes I_{m}-\alpha H-\lambda
I_{mn})(\frac{1}{\alpha})(\overline{W}\otimes I_{m}-\lambda I_{mn}))=0.$ (36)

To make the calculations simpler we let $m=1$ (but it can be extended to any
$m>1$),

$\displaystyle\mbox{det}(I_{n})\mbox{det}((\overline{A}-\lambda
I_{n})(\overline{W}-\lambda I_{n})+\alpha\lambda H)=0$ (37)

Note that $\alpha=0$ is a root of (37), and for small values $\alpha>0$ is
admissible for stability. We find the other positive root
$\alpha=\overline{\alpha}$ of (37) and, due to continuity of $\sigma(M)$ as a
function of $\alpha$ [25], we can claim that
$\alpha\in(0~{}\overline{\alpha})$ leads to stability of $M(\alpha)$. Recall
that $\alpha=0$ gives $\mbox{det}((\overline{W}-\lambda
I_{n})(\overline{A}-\lambda I_{n}))=0$ and the eigen-spectrum follows as
$\sigma(\overline{A})\cup\sigma(\overline{W})$ with two zero eigenvalues.

For any $\lambda$ in the LHP, from the diagonal structure of $\alpha H$, one
can rewrite (37) (with some abuse of notation) as

$\displaystyle\mbox{det}($ $\displaystyle(\overline{A}-\lambda
I_{n}\pm\sqrt{\alpha|\lambda|H})(\overline{W}-\lambda
I_{n}\mp\sqrt{\alpha|\lambda|H})$
$\displaystyle\pm\sqrt{\alpha|\lambda|H}(\overline{W}-\overline{A}))=0$ (38)

Letting $\overline{A}=\overline{W}$ to make the last term zero gives,

$\displaystyle\mbox{det}(\overline{W}$ $\displaystyle-\lambda
I_{n}\mp\sqrt{\alpha|\lambda|H})=$
$\displaystyle\mbox{det}(\overline{W}-\lambda(I_{n}\mp\sqrt{\frac{\alpha
H}{|\lambda|}}))=0$ (39)

with $\lambda$ as the eigenvalue of $\overline{W}$. The new perturbed
eigenvalue (in the above) is $\lambda(1\pm\sqrt{\frac{\alpha H}{|\lambda|}})$.
One can see that the minimum $\overline{\alpha}$ to make this term zero
(towards instability) for a non-zero $\lambda$ is

$\displaystyle\overline{\alpha}=\operatorname*{argmin}_{\alpha}|1-\sqrt{\frac{\alpha
H}{|\lambda|}}|\geq\frac{\min\\{|\operatorname{Re}\\{\lambda_{j}\\}|\neq
0\\}}{\max\\{H_{ii}\\}}=\frac{|\operatorname{Re}\\{\lambda_{2}\\}|}{\gamma}$
(40)

where at the last term we used $H\preceq\gamma I_{n}$. Then, from perturbation
analysis as in Theorem 1, one can claim the stability range as

$\displaystyle
0<\alpha<\overline{\alpha}:=\frac{|\operatorname{Re}\\{\lambda_{2}\\}|}{\gamma}$
(41)

for which $M(\alpha)$ has only one isolated zero eigenvalue and the rest with
negative values. Note that having $H=\gamma I_{n}$, as in homogeneous
quadratics models for CPU balancing [3] or power scheduling [4], the above
gives a tighter bound on $\alpha$. In general, for
$\overline{W}\neq\overline{A}$ matrices, the admissible range changes to,

$\displaystyle
0<\alpha<\frac{\min\\{|\operatorname{Re}\\{\lambda_{2}(\overline{W})\\}|,|\operatorname{Re}\\{\lambda_{2}(\overline{A})\\}|\\}}{\gamma}=:\overline{\alpha}$
(42)

This is because choosing the diagonal blocks of $M_{0}$ as both equal to
either $\overline{A}$ or $\overline{W}$ gives
$\overline{\alpha}:=\frac{|\operatorname{Re}\\{\lambda_{2}(\overline{A})\\}|}{\gamma}$
or
$\overline{\alpha}:=\frac{|\operatorname{Re}\\{\lambda_{2}(\overline{W})\\}|}{\gamma}$
in (42), respectively. This upper-bound on $\alpha$ is oblivious to the size
of the system and adjacency matrices (only depends on their eigen-spectrum)
and holds for any value of $m\geq 1$.

In general, the stability analysis of hybrid systems as (10)-(13) are
challenging (see some relevant discussions in [28] and also, on hybrid
consensus setups, in [22, Section IX]). The above perturbation-based analysis
is the key component of the convergence proof along with the proper choice of
the Lyapunov function as in [2, 22]. The stability is discussed next in the
main theorem of this paper.

######  Theorem 2.

Let Assumptions (1)-(2) and conditions in Lemma 2, Lemma 3, and Theorem 1
hold. For $\alpha$ satisfying (42), the dynamics (4)-(5) (and its matrix form
(10)-(13)) converges to the invariant equilibrium
$[\mathbf{x}^{*}\otimes\mathbf{1}_{n};\mathbf{0}_{nm}]$, where
$\mathbf{x}^{*}\otimes\mathbf{1}_{n}$ denotes the optimizer of the problem
(1).

###### Proof.

From Lemma 3,
$[\mathbf{x}^{*}\otimes\mathbf{1}_{n};\mathbf{0}_{nm}]\in\mathbb{R}^{2mn}$
with $\mathbf{x}^{*}\in\mathbb{R}$ is the invariant set under dynamics
(10)-(13) and belongs to $\mathcal{N}(M)$. From Theorem 1 and for an
admissible range of $\alpha\in(0~{}\overline{\alpha})$, we showed that the
eigenvalues of $M(\alpha,t)$ remain in the LHP except for the zero eigenvalue
of multiplicity one. This is irrespective of the network topology
$\mathcal{G}$ and the change in $M(\alpha,t)$ which may even allow for
possible hybrid stability analysis [29]. Following [22, Corollary 1], this
network of integrators is asymptotically globally stable. On the other hand,
its equilibrium is associated with the (eigenspace of) right eigenvectors
associated with its zero eigenvalue of multiplicity $m$, which is a set of $m$
decoupled vectors as $\\{[\mathbf{1}_{n};\mathbf{0}_{n}]\otimes I_{m}\\}$.
This simply implies that the invariant set
$[\mathbf{x}^{*}\otimes\mathbf{1}_{n};\mathbf{0}_{nm}]$ is the equilibrium
under dynamics (10)-(13). ∎

For undirected networks, since the eigenvalues of $M$ are real, one extends
the proof to hybrid setups under a switching signal $q$.

##  IV The Discrete-Time Model: Improvement over the Existing Models

###  IV-A Discrete-time version

For agents with discrete-time models, one can follow similar discussions as
for discretized consensus dynamics in [22, Section IV]. The discretized
version of the CT dynamics (10) takes the following form,

$\displaystyle\left(\begin{array}[]{c}\mathbf{x}(k+1)\\\
\mathbf{y}(k+1)\end{array}\right)=M_{d}(\eta,\alpha)\left(\begin{array}[]{c}{\mathbf{x}(k)}\\\
{\mathbf{y}(k)}\end{array}\right),$ (47)

with $k$ as the discrete time-step. The most common approximation for $M_{d}$
is by Euler Forward Method (EFM) as $M_{d}=I+\eta M$ with $\eta$ as the
(discrete) sampling step-size. This implies that the discretized version of
the node dynamics (4)-(5) is in the following form:

$\displaystyle\mathbf{x}_{i}(k+1)=\mathbf{x}_{i}(k)-$
$\displaystyle\eta\sum_{j=1}^{n}w^{q}_{ij}(\mathbf{x}_{i}(k)-\mathbf{x}_{j}(k))-\alpha\mathbf{y}_{i}(k),$
(48) $\displaystyle\mathbf{y}_{i}(k+1)=\mathbf{y}_{i}(k)-$
$\displaystyle\eta\sum_{j=1}^{n}a^{q}_{ij}(\mathbf{y}_{i}(k)-\mathbf{y}_{j}(k))$
$\displaystyle+\boldsymbol{\nabla}f_{i}(\mathbf{x}_{i}(k+1))-\boldsymbol{\nabla}f_{i}(\mathbf{x}_{i}(k)),$
(49)

As discussed in Section III-B, matrix $M$ has $m$ zero eigenvalues while all
other eigenvalues have negative real-part. This implies that $M_{d}$ given by
EFM has $m$ eigenvalues at $1$ with $m$ independent (decoupled) eigenvectors,
while for stability, its other eigenvalues need to satisfy $|\lambda_{i}|<1$
which depends on the sampling step $\eta$. The explicit upper-bound for $\eta$
such that stable CT dynamics remains stable after discretization via first-
order Euler approximation is given, for example, in [30, Table I] as,

$\displaystyle\eta<\min_{2\leq i\leq 2n,1\leq j\leq
m}\frac{2|Re\\{\lambda_{i,j}(\alpha)\\}|}{|\lambda_{i,j}(\alpha)|^{2}},$ (50)

with $\lambda_{i,j}(\alpha)$ as the eigenvalue of the time-varying matrix $M$.
For real Laplacian matrices associated with linear consensus protocols, e.g.,
$\overline{W}$ the discretized model is given as
$P_{w}=I_{n}+\eta\overline{W}$ which satisfies the following Lemma.

######  Lemma 7.

[31] Given $P_{w}=I_{n}+\eta\overline{W}$ with
$\overline{w}_{\max}=\sum_{j=1,j\neq i}^{n}|\overline{w}_{ij}|$ (from Lemma
6), $\overline{W}$ as the Laplacian matrix satisfying Assumption 2, and
$0<\eta<\frac{1}{\overline{w}_{\max}}$, the followings hold:

  * •

$P_{w}$ is row-stochastic and non-negative with one isolated eigenvalue at
$1$.

  * •

The eigenvalues of $P_{w}$ are inside the unit circle.

  * •

$P_{w}$ is primitive, i.e., it has only one eigenvalue with maximum modulus.

One can state Lemma 7 similarly for $\overline{A}$. We use this for the
stability analysis of the DT version (47) via a similar procedure as in
Section III-B by proving that $M_{d}$ has $m$ eigenvalues at $1$ with all
other eigenvalues within the unit circle. Assuming ${M_{d}=(I+\eta
M_{0})+\alpha\eta M_{1}}$, it is clear that for $\eta$ satisfying Lemma 7 all
the eigenvalues of $I+\eta M_{0}$ remain inside the unit circle except for
$2m$ eigenvalues $\lambda_{1,j}=\lambda_{2,j}=1$ (from Lemma 7). In this case,
we consider $\eta\alpha M_{1}$ as a perturbation to $I+\eta M_{0}$ and we
follow a similar analysis as in the proof of Theorem 1. Note that the
eigenvectors associated with the $1$ eigenvalues are the same as (25)-(28).
Then, from Lemma 4 and following similar equation to derive Eq. (32), we have
$\frac{d\lambda_{1,j}}{d\alpha}|_{\alpha=0}=0$ and
$\frac{d\lambda_{2,j}}{d\alpha}|_{\alpha=0}<0$. This implies that $M_{d}$ has
$m$ eigenvalue at $1$ (with the same right eigenvectors $V_{1}$) and for
sufficiently small $\alpha\eta$ all other eigenvalues remain inside the unit
circle. To bound $\alpha\eta$, we follow a similar procedure as in Section
III-C. From Lemma 5, we need $d(\sigma(M_{d}),\sigma(I+\eta
M_{0}))<1-\lambda_{\max}$ where $\lambda_{\max}=\max_{1\leq j\leq m,2\leq
i\leq n}|\lambda_{i,j}|$ as the eigenvalue (with its maximum absolute value
less than $1$) of $I+\eta M_{0}$, and the bound on $\eta\alpha$ can be
determined. From Lemma 6, for $0$-$1$ adjacency matrix $W$,
$\overline{w}_{\max}$ denotes the max node degree [22], and the algebraic
connectivity of $\overline{W}$ satisfies
$\lambda_{2}(\overline{W})\geq\frac{1}{nd_{g}}$ with $d_{g}$ as the network
diameter [32, p.571]. This gives an idea of the eigen-spectrum of $M_{0}$ and
the discretized matrix $I+\eta M_{0}$. Using perturbation analysis, one can
find an estimate on the eigen-spectrum of $M_{d}(\alpha)$ by considering the
perturbation term as $\eta\alpha M_{1}$. Then, Eq. (34) changes to the
following which gives one admissible range of $\alpha\eta$,

$\displaystyle
0<\alpha\eta\leq\frac{(1-\lambda_{\max})^{nm}}{4^{nm}(3-\lambda_{\max})^{nm-1}\max\\{1,\gamma\\}}$
(51)

Similar to the derivation resulting to Eq. (42), one can find the min
$\overline{\alpha}$ towards instability for $|\lambda|\neq 1$ as

$\displaystyle\overline{\alpha}$
$\displaystyle=\operatorname*{argmin}_{\alpha}|1-\sqrt{\frac{\alpha\eta
H}{|\lambda-1|}}|$
$\displaystyle\geq\frac{1-\lambda_{\max}}{\eta\max\\{H_{ii}\\}}=\frac{1-\lambda_{\max}}{\eta\gamma}$
(52)

and the bound in Eq. (42) changes to,

$\displaystyle
0<\alpha\eta<\frac{\min\\{1-\lambda_{\max}(\overline{A}),1-\lambda_{\max}(\overline{W})\\}}{\gamma}$
(53)

where the bound is on both $\alpha\eta$ as the perturbation parameter. For
undirected networks, i.e., symmetric $\overline{W}$, $\overline{A}$ with real
eigenvalues and eigenvectors, for simplicity assume $m=1$ where the results
can be easily extended to $m>1$. Three cases may occur for
$\lambda_{\max}(\overline{W}),\lambda_{\max}(\overline{A})$ which are largest
(in absolute value) eigenvalues inside the unit circle: (i)
$|\lambda_{\max}(\overline{W})|>|\lambda_{\max}(\overline{A})|$ with vector
$[\mathbf{v}_{w};\mathbf{0}_{n}]$ as the right eigenvector of $M_{0}$
associated to it. (ii)
$\lambda_{\max}(\overline{W})<\lambda_{\max}(\overline{A})$ with
$[\mathbf{0}_{n};\mathbf{v}_{a}]$ as its right eigenvector. (iii)
$\lambda_{\max}(\overline{W})=\lambda_{\max}(\overline{A})$ (algebraic
multiplicity of two) with
$[\mathbf{v}_{w},\mathbf{0}_{n};\mathbf{0}_{n},\mathbf{v}_{a}]$ as the right
eigenvectors. Then, from Lemma 4 and following similar equation to derive Eq.
(32), for case (i),

$\displaystyle[\mathbf{v}_{w};\mathbf{0}_{n}]^{\top}M_{1}[\mathbf{v}_{w};\mathbf{0}_{n}]=0.$
(54)

which implies that, after perturbation, we have
$\partial_{\alpha}\lambda_{\max}|_{\alpha=0}=0$. For case (ii),

$\displaystyle[\mathbf{0}_{n};\mathbf{v}_{a}]^{\top}M_{1}[\mathbf{0}_{n};\mathbf{v}_{a}]=-\sum_{i=1}^{n}\mathbf{v}^{2}_{a,i}\partial_{\mathbf{x}_{i}}f_{i}<0.$
(55)

which implies that perturbation leads to
$\partial_{\alpha}\lambda_{\max}|_{\alpha=0}<0$ and $\lambda_{\max}<1$ moves
further inside the unit circle and gets smaller in absolute value. For case
(iii) as a combination of the other two cases,

$\displaystyle\left(\begin{array}[]{cc}\mathbf{v}_{w}&\mathbf{0}_{n}\\\
\mathbf{0}_{n}&\mathbf{v}_{a}\end{array}\right)^{\top}M_{1}\left(\begin{array}[]{cc}\mathbf{v}_{w}&\mathbf{0}_{n}\\\
\mathbf{0}_{n}&\mathbf{v}_{a}\end{array}\right)=\left(\begin{array}[]{cc}0&\times\\\
0&-\mathbf{v}_{a}^{\top}H\mathbf{v}_{a}\end{array}\right).$ (62)

saying that $\partial_{\alpha}\lambda_{\max}(\overline{W})|_{\alpha=0}=0$,
$\partial_{\alpha}\lambda_{\max}(\overline{A})|_{\alpha=0}<0$. This implies
that $\lambda_{\max}(\overline{A})$ gets smaller in size, and moves further
towards inside the unit circle and $\lambda_{\max}(\overline{W})<1$ remains
constant. Since the perturbation parameter is $\alpha\eta$ in Eq. (53),
$\alpha$ needs to be decreased for larger values of the sampling period
$\eta$. Following Lemma 5, one can claim that
$|\lambda_{\max}(\alpha)|\geq|\lambda_{\max}|+d(\sigma(M(\alpha)),\sigma(M_{0}))$
with $|\lambda_{\max}|:=\max_{2\leq i\leq 2n,1\leq j\leq m}|\lambda_{i,j}|$
(for $\alpha=0$). For a given $0<\alpha<\overline{\alpha}$, then Eq. (50)
gives another bound as,

$\displaystyle\eta<\frac{1}{\max\\{\overline{w}_{\max},\overline{a}_{\max}\\}+2(\lVert
M_{0}\rVert+\lVert M(\alpha)\rVert)^{1-\frac{1}{2mn}}\lVert\alpha
M_{1}\rVert^{\frac{1}{2mn}}},$ (63)

ensures stable system mapping from CT to DT and convergence under the
discretized model (47). We summarized our proposed distributed optimization
solution in Algorithm 1.

1 Given: costs $f_{i}(\mathbf{x}_{i})$, dynamic graph $\mathcal{G}_{q}$,
weight matrices $W_{q},A_{q}$, switching signal $q$

2 Initialization: ${\mathbf{y}}_{i}(0)=\mathbf{0}_{m}$, initialize
${\mathbf{x}}_{i}(0)$ with random values ;

3 Choose $\eta$ satisfying Eq. (63) and $0<\alpha<\overline{\alpha}$
satisfying Eq. (42)

4 for  _$k <T_{termination}$_ do

5 Node $i$ receives $\mathbf{x}_{j}(k)$ and $\mathbf{y}_{j}(k)$ data from the
neighbors $j$;

6 Node $i$ finds $\mathbf{x}_{i}(k+1)$ via discrete dynamics (48);

7 Node $i$ finds $\boldsymbol{\nabla}f_{i}(\mathbf{x}_{i}(k+1))$;

8 Node $i$ finds $\mathbf{y}_{i}(k+1)$ via discrete dynamics (49);

9 Node $i$ shares $\mathbf{x}_{i}(k+1)$ and $\mathbf{y}_{i}(k+1)$ over
$\mathcal{G}_{q}$;

10 Set $k\leftarrow k+1$

11Return: final optimal state $\mathbf{x}^{*}$ and final cost
$F(\mathbf{x}^{*})$;

Algorithm 1 The proposed distributed optimization algorithm.

###  IV-B Performance Subject to Link Removal

The proposed discrete strategy shows better performance over dynamic networks
and switching topologies. This is because the WB condition for matrices $W$
and $A$ (Assumption 2) is much easier to satisfy as compared to weight-
stochastic design algorithms, e.g., in existing discrete-time literature [9,
11, 12]. Note that many existing works do not work under switching network
topologies or need weight-redesign algorithms to satisfy stochastic conditions
[13, 14] in case of link failure. Since in general weigh-balancing algorithms
[33] are computationally more efficient than stochastic design algorithms and
in some cases, it is not even possible to satisfy the weight stochasticity for
some networks (see example in Fig. 2-(Right)), the proposed strategy
outperforms the existing literature. This specifies a case where there is no
need to redesign the weights for the WB condition (for weight-symmetric
undirected networks) but needs the weights to be redesigned (if possible) for
the stochastic condition. Some other examples in distributed resource
allocation and coupling-constrained distributed optimization subject to packet
drops and link removals are given in [34]. In [34] it is claimed that such a
WB condition is more preferred in real-time applications as it requires no
need to redesign stochastic algorithms in case the network is unreliable and
dynamic, e.g., in packet loss scenarios.

######  Remark 3.

For general symmetric WB undirected graphs with unreliable links, the WB
conditions hold after link removal/failure, but the bi-stochastic condition
does not necessarily hold for the same networks after link removal. For
example, consider the illustration in Fig. 2.

Figure 2: (Left)-An undirected unreliable network of $5$ nodes which is both
stochastic and WB. The red link represents an unreliable link that might be
subject to failure. (Right)-The same network after removal/failure of the (bi-
directional) unreliable link: the network is still WB, but is not bi-
stochastic anymore. To apply the existing weight-stochastic literature, e.g.,
[9, 11, 12], one needs to redesign the weights, e.g., by weight compensation
algorithms in [13, 14], to redesign the weights if possible and make the new
topology weight-stochastic again.

##  V Simulation: Nonlinear SVM

Figure 3: This figure shows the data points scattered in 2D and classified by
the SVM hyperplane in 3D.

Figure 4: The SVM parameters ${\omega}_{i}$ and $\nu_{i}$ (at all $n=5$ nodes
in Fig. 2) under discrete dynamics (10)-(13). The top row represents the
parameters before link removal and the bottom row is for after link removal.
The gradient sum $\sum_{i=1}^{n}\boldsymbol{\nabla}f_{i}(\mathbf{x}_{i})$ is
also shown to admit the GT property. The hinge loss residual
$\overline{F}(\mathbf{x})=\sum_{i=1}^{n}f_{i}(\mathbf{x}_{i})-F^{*}$ (with
$F^{*}$ obtained via centralized SVM) is also shown.

Here, we recall an application in distributed support-vector-machines (D-SVM)
from [2]. The cost function to be optimized locally is in logarithmic hinge
loss form as

$\displaystyle
f_{i}(\mathbf{x}_{i})=\boldsymbol{\omega}_{i}^{\top}\boldsymbol{\omega}_{i}+C\sum_{j=1}^{N_{i}}\tfrac{1}{\mu}\log(1+\exp(\mu
z_{i})),$ (64)
$\displaystyle{\mathbf{x}_{i}=[\boldsymbol{\omega}_{i}^{\top};\nu_{i}]}\in\mathbb{R}^{m},~{}~{}{z_{i}=1-l_{j}(\boldsymbol{\omega}_{i}^{\top}\phi(\boldsymbol{\chi}^{i}_{j})-\nu_{i})}$
(65)

with $\phi(\cdot)$ as some proper nonlinear mapping. The simulation are done
for the 3D binary classification example with $m=4$ and $50$ random data
points. The loss function parameters are set as: $\mu=3$, $C=1.5$. The points
are not linearly separable in 2D and transformed into 3D using proper
nonlinear mapping in the form
$\phi(\boldsymbol{\chi}^{i}_{j})=[\boldsymbol{\chi}^{i}_{j}(1)^{2};\boldsymbol{\chi}^{i}_{j}(2)^{2};\sqrt{2}\boldsymbol{\chi}^{i}_{j}(1)\boldsymbol{\chi}^{i}_{j}(2)]$.
The projected 3D points are linearly separable via SVM hyperplane and proper
Kernel function. It should be clarified that every node has access to a
different part of the SVM classification dataset. In fact, only $80\%$ of
randomly chosen data points are assigned to each node. The SC graph of $n=5$
nodes is shown in Fig. 2 (before and after link removal). For the sake of this
simulation, the link weights are chosen randomly and symmetric, and one of the
links is removed to show the WB performance subject to link failure. In Fig.
4, the top row shows the algorithm evolution before link failure and the
bottom row is for after link removal. We set $\alpha=1,\eta=0.001$ in
Algorithm 1 for this simulation and choose random initialization for D-SVM
parameters. As it is clear, the proposed algorithm is still convergent under
link failure. This shows how our work advances the state-of-the-art literature
which assume stochastic link-weights. Note that stochastic weight design
suffers from change in the network topology (including link failure), since
redesign algorithms are needed to compensate for loss of stochasticity. Some
of these compensation algorithms are proposed in [13, 14, 15, 16, 17]. These
algorithms add more computational complexity to the existing bi-stochastic
weight design solutions. Therefore, such algorithms cannot handle the link
failure example in Fig. 2, while our algorithm works with no additional
computational complexity. Note that in Fig. 4 the D-SVM parameters at all
nodes reach consensus and converge to the centralized SVM parameters
asymptotically, while each node has only access to a portion of the
classification dataset. This is significant as it illustrates how our
algorithm paves the way for distributed and localized machine learning
solutions.

Figure 5: The D-SVM residual (at all $n=6$ nodes) under discrete dynamics
(10)-(13), where the weight matrices are time-varying (random-weighting at
every 100 iterations) over a 2-hop SC digraph. The hinge loss residual
$\overline{F}$ are compared for different $\alpha,\eta$ values. This example
shows that for large values of both $\alpha=5,\eta=0.005$ solution is not
stable.

We repeat the simulation for a dynamic random network with time-varying link
weights randomly changing every $100$ iterations. This is to show that our
setup can handle dynamic network topologies. The residuals
$\overline{F}(\mathbf{x})=\sum_{i=1}^{n}f_{i}(\mathbf{x}_{i})-F^{*}$ (with
$F^{*}$ obtained via centralized SVM) are shown in Fig. 5 for different values
of $\alpha,\eta$ in Algorithm 1: $\alpha=5,\eta=0.005$, $\alpha=8,\eta=0.001$,
$\alpha=1,\eta=0.01$. As it can be seen there is a trade-off between these two
variables, and by increasing one the other is decreased to satisfy the
stability. For large values of both GT tracking parameter $\alpha$ and
discrete step-size $\eta$ the stable convergence may not be achieved, e.g.,
see $\alpha=5$ and $\eta=0.005$. This is because (following the perturbation-
based analysis in Section IV) for these values of $\alpha$ and $\eta$ one or
more eigenvalues of the proposed dynamics (47) move to the RHP and causes
instability. On the other hand, for other choices of $\alpha$ and $\eta$
satisfying Eq. (53) all eigenvalues remain in the LHP and stable convergence
holds.

##  VI Conclusions and Future Directions

This work relaxes the weight constraint on the consensus optimization
algorithms, e.g., to be applied over dynamic networks. One main application is
in distributed mobile sensor networks and learning over coordinated (swarm)
robotic systems, where the connectivity mainly follows the nearest neighbour
rule [21] or (limited) broadcasting range of the mobile entities (or agents)
with the links coming and going as their formation evolves. As one direction
of future research, one can further address possible communication time-delays
or link failure (similar to the equality-constraint optimization in [35, 36]).
Also, actuation and data-transmission nonlinearities on the node dynamics,
e.g., in terms of quantization or saturation for equality-constraint
optimization [37, 3], can be extended to this current work.

## Acknowledgements

The authors would like to thank Themistoklis Charalambous and Usman A. Khan
for their help and comments.

## References

  * [1] M. Doostmohammadian, A. Taghieh, and H. Zarrabi,  “Distributed estimation approach for tracking a mobile target via formation of UAVs,”  IEEE Transactions on Automation Science and Engineering, vol. 19, no. 4, pp. 3765–3776, 2021.
  * [2] M. Doostmohammadian, A. Aghasi, T. Charalambous, and U. A. Khan,  “Distributed support vector machines over dynamic balanced directed networks,”  IEEE Control Systems Letters, vol. 6, pp. 758 – 763, 2021.
  * [3] A. I. Rikos, A. Grammenos, E. Kalyvianaki, C. N. Hadjicostis, T. Charalambous, and K. H. Johansson,  “Optimal CPU scheduling in data centers via a finite-time distributed quantized coordination mechanism,”  in 51st IEEE Conf. on Decision and Control, 2021.
  * [4] S. Kar and G. Hug,  “Distributed robust economic dispatch in power systems: A consensus \+ innovations approach,”  in IEEE Power and Energy Society General Meeting, 2012, pp. 1–8.
  * [5] L. Liu and G. Yang,  “Distributed fixed-time optimal resource management for microgrids,”  IEEE Transactions on Automation Science and Engineering, 2022.
  * [6] M. Kaheni, E. Usai, and M. Franceschelli,  “A distributed optimization and control framework for a network of constraint coupled residential BESSs,”  in IEEE 17th International Conference on Automation Science and Engineering (CASE). IEEE, 2021, pp. 2202–2207.
  * [7] M. Doostmohammadian,  “Distributed energy resource management: All-time resource-demand feasibility, delay-tolerance, nonlinearity, and beyond,”  IEEE Control Systems Letters, 2023.
  * [8] T. T. Doan and A. Olshevsky,  “Distributed resource allocation on dynamic networks in quadratic time,”  Systems & Control Letters, vol. 99, pp. 57–63, 2017.
  * [9] F. Saadatniaki, R. Xin, and U. A. Khan,  “Decentralized optimization over time-varying directed graphs with row and column-stochastic matrices,”  IEEE Transactions on Automatic Control, vol. 65, no. 11, pp. 4769–4780, 2020.
  * [10] B. Gharesifard and J. Cortés,  “Distributed continuous-time convex optimization on weight-balanced digraphs,”  IEEE Transactions on Automatic Control, vol. 59, no. 3, pp. 781–786, 2014.
  * [11] C. Xi, R. Xin, and U. A. Khan,  “ADD-OPT: accelerated distributed directed optimization,”  IEEE Trans. on Autom. Control, vol. 63, no. 5, pp. 1329–1339, 2017\.
  * [12] S. Pu, W. Shi, J. Xu, and A. Nedić,  “Push–pull gradient methods for distributed optimization in networks,”  IEEE Transactions on Automatic Control, vol. 66, no. 1, pp. 1–16, 2021.
  * [13] N. H. Vaidya, C. N. Hadjicostis, and A. D. Dominguez-Garcia,  “Robust average consensus over packet dropping links: Analysis via coefficients of ergodicity,”  in 51st IEEE Conference on Decision and Control, 2012, pp. 2761–2766.
  * [14] F. Fagnani and S. Zampieri,  “Average consensus with packet drop communication,”  SIAM Journal on Control and Optimization, vol. 48, no. 1, pp. 102–133, 2009.
  * [15] J. Xu, H. Zhang, and L. Xie,  “Consensusability of multiagent systems with delay and packet dropout under predictor-like protocols,”  IEEE Transactions on Automatic Control, vol. 64, no. 8, pp. 3506–3513, 2018.
  * [16] Z. Li and J. Chen,  “Robust consensus for multi-agent systems communicating over stochastic uncertain networks,”  SIAM Journal on Control and Optimization, vol. 57, no. 5, pp. 3553–3570, 2019.
  * [17] B. Gerencsér and J. M. Hendrickx,  “Push-sum with transmission failures,”  IEEE Transactions on Automatic Control, vol. 64, no. 3, pp. 1019–1033, 2018.
  * [18] R. Bhatia,  Perturbation bounds for matrix eigenvalues,  SIAM, 2007.
  * [19] K. Rokade and R. K. Kalaimani,  “Distributed ADMM over directed networks,” 2021,  https://arxiv.org/abs/2010.10421.
  * [20] W. Jiang and T. Charalambous,  “Distributed alternating direction method of multipliers using finite-time exact ratio consensus in digraphs,”  in 2021 European Control Conference (ECC), 2021, pp. 2205–2212.
  * [21] A. Jadbabaie, J. Lin, and A. S. Morse,  “Coordination of groups of mobile autonomous agents using nearest neighbor rules,”  IEEE Transactions on automatic control, vol. 48, no. 6, pp. 988–1001, 2003.
  * [22] R. Olfati-Saber and R. M. Murray,  “Consensus problems in networks of agents with switching topology and time-delays,”  IEEE Transactions on Automatic Control, vol. 49, no. 9, pp. 1520–1533, Sept. 2004.
  * [23] A. P. Seyranian and A. A. Mailybaev,  Multiparameter stability theory with mechanical applications, vol. 13,  World Scientific, 2003.
  * [24] K. Cai and H. Ishii,  “Average consensus on general strongly connected digraphs,”  Automatica, vol. 48, no. 11, pp. 2750–2761, 2012.
  * [25] G. W. Stewart and J. Sun,  “Matrix perturbation theory,”  1990\.
  * [26] R. A. Horn and C. R. Johnson,  Matrix Analysis,  Cambridge University Press, Cambridge, 1985.
  * [27] M. Doostmohammadian, M. Pirani, U. A. Khan, and T. Charalambous,  “Consensus-based distributed estimation in the presence of heterogeneous, time-invariant delays,”  IEEE Control Systems Letters, vol. 6, pp. 1598 – 1603, 2021.
  * [28] J. Cortés,  “Discontinuous dynamical systems,”  IEEE Control systems magazine, vol. 28, no. 3, pp. 36–73, 2008\.
  * [29] R. Goebel, R. Sanfelice, and A. Teel,  “Hybrid dynamical systems,”  IEEE control systems magazine, vol. 29, no. 2, pp. 28–93, 2009\.
  * [30] P. Axelsson and F. Gustafsson,  “Discrete-time solutions to the continuous-time differential Lyapunov equation with applications to Kalman filtering,”  IEEE Transactions on Automatic Control, vol. 60, no. 3, pp. 632–643, 2014.
  * [31] R. Olfati-Saber, J. A. Fax, and R. M. Murray,  “Consensus and cooperation in networked multi-agent systems,”  Proceedings of the IEEE, vol. 95, no. 1, pp. 215–233, 2007.
  * [32] J. L. Gross and J. Yellen,  Handbook of Graph Theory,  CRC Press, 2004.
  * [33] T. Charalambous, C.N. Hadjicostis, and M. Johansson,  “Distributed minimum-time weight balancing over digraphs,”  in 6th International Symposium on Communications, Control and Signal Processing, May 2014, pp. 190–193.
  * [34] M. Doostmohammadian, U. Khan, and A. Aghasi,  “Distributed constraint-coupled optimization over unreliable networks,”  in 10th RSI International Conference on Robotics and Mechatronics (ICRoM), Nov 2022.
  * [35] M. Doostmohammadian, A. Aghasi, M. Vrakopoulou, H. R. Rabiee, U. A. Khan, and T. Charalambous,  “Distributed delay-tolerant strategies for equality-constraint sum-preserving resource allocation,”  Systems & Control Letters, vol. 182, pp. 105657, 2023.
  * [36] F. Rahimi and R. Mahboobi Esfanjani,  “A distributed dual decomposition optimization approach for coordination of networked mobile robots with communication delay,”  in 9th RSI International Conference on Robotics and Mechatronics. 2021, pp. 18–23, IEEE.
  * [37] M. Doostmohammadian, A. Aghasi, M. Vrakopoulou, and T. Charalambous,  “1st-order dynamics on nonlinear agents for resource allocation over uniformly-connected networks,”  in IEEE Conference on Control Technology and Applications (CCTA), 2022, pp. 1184–1189.

|  Mohammadreza Doostmohammadian received his B.Sc. and M.Sc. in Mechanical
Engineering from Sharif University of Technology, Iran, respectively in 2007
and 2010, where he worked on applications of control systems and robotics. He
received his PhD in Electrical and Computer Engineering from Tufts University,
MA, USA in 2015. During his PhD in Signal Processing and Robotic Networks
(SPARTN) lab, he worked on control and signal processing over networks with
applications in social networks. From 2015 to 2017 he was a postdoc researcher
at ICT Innovation Center for Advanced Information and Communication Technology
(AICT), School of Computer Engineering, Sharif University of Technology, with
research on network epidemic, distributed algorithms, and complexity analysis
of distributed estimation methods. He was a researcher at Iran
Telecommunication Research Center (ITRC), Tehran, Iran in 2017 working on
distributed control algorithms and estimation over IoT. Since 2017 he has been
an Assistant Professor with the Mechatronics Department at Semnan University,
Iran, and he was a researcher at the School of Electrical Engineering and
Automation, Aalto University, Finland. His general research interests include
distributed estimation, control, and convex optimization over networks. He was
the chair of the robotics and control session at the ISME-2018 conference and
also the session chair at 1st Artificial Intelligent Systems Conference of
Iran, 2022.
---|---
|  Wei Jiang got his Ph. D. degree from Ecole Centrale de Lille, France,
working on cooperative control for multi-agent/robot systems and motion
planning for mobile manipulator systems. He is a Postdoctoral researcher at
the School of Electrical Engineering and Automation, Aalto University,
Finland.
---|---
|  Muwahida Liaquat received her PhD in Electrical Engineering in 2013 from
the National University of Sciences and Technology (NUST), Pakistan,
specializing in Control Systems. Her work was focused on developing a
reconstruction filter for the sampled data regulation of linear systems. She
received her MS in Electrical Engineering in 2006, and her BE in Computer
Engineering in 2004 from NUST. She is currently working as an Assistant
Professor in the Department of Electrical Engineering at NUST College of
Electrical and Mechanical Engineering, where she is part of the Control
Systems research group along with teaching specialized courses in control
engineering to both undergraduate and postgraduate levels. She has supervised
a number of MS students in the area of control engineering. She was a
Postdoctoral researcher at the School of Electrical Engineering and
Automation, Aalto University, Finland. She is currently a senior research
scientist with Finnish Geospatial Research Institute (FGI).
---|---
|  Alireza Aghasi joined the School of Electrical Engineering and Computer
Science at Oregon State University in Fall 2022. Between 2017 and 2022, he was
an assistant professor in the Department of Data Science and Analytics at the
Robinson College of Business, Georgia State University. Prior to this
position, he was a research scientist with the Department of Mathematical
Sciences, IBM T.J. Watson Research Centre, Yorktown Heights. From 2015 to 2016
he was a postdoctoral associate with the computational imaging group at the
Massachusetts Institute of Technology, and between 2012 and 2015 he served as
a postdoctoral research scientist with the compressed sensing group at Georgia
Tech. He got his PhD from Tufts University in 2012. His research fundamentally
focuses on optimization theory and statistics, with applications to various
areas of data science, artificial intelligence, modern signal processing and
physics-based inverse problems.
---|---
|  Houman Zarrabi received his PhD from Concordia University in Montreal,
Canada in 2011. Since then he has been involved in various industrial and
research projects. His main expertise includes IoT, M2M, CPS, big data,
embedded systems, and VLSI. He is the national IoT program director and
Assistant Professor at the Iran Telecommunication Research Center (ITRC).
---|---