A Design Method of Distributed Algorithms via Discrete-time Blended Dynamics Theorem

We use the following convention in this paper. $\mathmybb{0}_{N}$ and $\mathmybb{1}_{N}$ denote the column vectors of size $N$ consisting of all zeros and ones, respectively. A positive vector implies a vector whose elements are all positive. The 2-norm of a vector and the induced 2-norm of a matrix are denoted by $\|\cdot\|$. For any square matrix $\mathsf{M}$, $\rho(\mathsf{M})$ denotes the spectral radius of $\mathsf{M}$. The operation defined by the symbol $\otimes$ is Kronecker product, which has the properties that, given the matrices $\mathsf{A,B,C}$ and $\mathsf{D}$ of appropriate dimensions, $(\mathsf{A}\otimes\mathsf{B})(\mathsf{C}\otimes\mathsf{D})=(\mathsf{AC})\otimes(\mathsf{BD})$ and $\|\mathsf{A}\otimes\mathsf{B}\|=\|\mathsf{A}\|\|\mathsf{B}\|$. See, e.g., (Bernstein, 2009) for more about the Kronecker product. For notational convenience, we use the convention $\mathsf{M}_{\otimes n}:=\mathsf{M}\otimes I_{n}$ for any size of matrix $\mathsf{M}$.

A communication network of $N$ agents can be represented by a directed graph ${\mathcal{G}}=({\mathcal{N}},{\mathcal{E}})$ where ${\mathcal{N}}=\left\{1,2,\dots,N\right\}$ is the node set, and ${\mathcal{E}}\subset{\mathcal{N}}\times{\mathcal{N}}$ is the edge set of ordered pairs of nodes. If agent $j$ sends information to agent $i$, then we say that the node $j$ is connected to the node $i$ by an edge $(j,i)\in{\mathcal{E}}$. A graph is said to be undirected if $(j,i)\in{\mathcal{E}}$ implies $(i,j)\in{\mathcal{E}}$. For a directed graph, in-neighbors of node $i$ are represented by ${\mathcal{N}}_{i}=\{j\in{\mathcal{N}}|(j,i)\in{\mathcal{E}}\}$ and in-degree of node $i$ is denoted by $d_{i}:=|{\mathcal{N}}_{i}|$. On the contrary, out-neighbors of node $i$ are represented by ${\mathcal{N}}_{i}^{\text{out}}:=\{j\in{\mathcal{N}}|(i,j)\in{\mathcal{E}}\}$ and out-degree of node $i$ is denoted by $d_{i}^{\text{out}}:=|{\mathcal{N}}_{i}^{\text{out}}|$. A path of length $L$ from node $i$ to node $j$ is a sequence $(i_{0},i_{1},\ldots,i_{L})$ such that $i_{0}=i$, $i_{L}=j$, $(i_{l},i_{l+1})\in{\mathcal{E}}$ for any $l\in\{0,\ldots,L-1\}$, and every $i_{l}$ is distinct. A directed graph is said to be strongly connected if there exists a path from any node to any other node. Similarly, an undirected graph is said to be connected if there exists a path between any two nodes. A cycle is a path that starts and ends at the same node and any node does not appear more than once in it. A strongly connected directed graph is said to be periodic if there exists period $\mathsf{p}>1$ that divides the length of every cycle in the graph, otherwise the graph is said to be aperiodic. Any directed graph having at least one self-connection is aperiodic. For a directed graph ${\mathcal{G}}$, its adjacency matrix $\mathcal{A}=[\alpha_{ij}]\in{\mathbb{R}}^{N\times N}$ is defined as $\alpha_{ij}>0$ if $(j,i)\in{\mathcal{E}}$ and $\alpha_{ij}=0$ otherwise. If every positive component of $\mathcal{A}$ is 1, then $\mathcal{A}$ is called a binary adjacency matrix. Meanwhile, given a non-negative square matrix $\mathsf{M}\in{\mathbb{R}}^{N\times N}$, its associated directed graph is defined as the directed graph whose adjacency matrix is $\mathsf{M}$. A non-negative square matrix $\mathsf{M}\in{\mathbb{R}}^{N\times N}$ is said to be primitive if there exists $\mathsf{n}\in{\mathbb{N}}$ such that $\mathsf{M}^{\mathsf{n}}$ is positive, i.e., every component of $\mathsf{M}^{\mathsf{n}}$ is positive. The primitive property of any non-negative square matrix can be verified by its associated graph as the following lemma.

In this subsection, we assume that the network is undirected and connected, and consider the following Metropolis-Hastings coupling weight $w_{ij}^{\text{MH}}$ in (Schwarz et al., 2014):

where $\mu\in(0,1)$. Note that $w_{ij}^{\text{MH}}$ depends on both $d_{i}$ and $d_{j}$, so that each agent should additionally exchange its degree information with neighbors to update its coupling weights online, and this will enable the plug-and-play operation to the multi-step coupling framework.

It can be easily seen that $w_{ij}^{\text{MH}}$ satisfies (7) and Assumption 1 holds for $W^{\text{MH}}:=[w_{ij}^{\text{MH}}]$ because it is doubly-stochastic. Hence, we can choose $p=\mathmybb{1}_{N}$ and $q=(1/N)\mathmybb{1}_{N}$ from Lemma 2.1. Then, the blended dynamics is obtained as a simple average of $f_{i}$s as follows:

Since all $p_{i}$s are chosen evenly as 1, by Theorem 1 or Corollary 1, the behavior of every trajectory ${x}_{i}[t]$ is approximately synchronized to the solution $s$ of the blended dynamics (21). It should be emphasized that this collective synchronized behavior comes from $p=\mathmybb{1}_{N}$ and this choice of $p$ is always possible for any row-stochastic (not necessarily to be doubly-stochastic) weight matrices. In Section 3.3, we will consider the row-stochastic weight matrix whose column-sums are not 1.

As an application example for the Metropolis-Hastings coupling, we can design a distributed algorithm for network size estimation as follows. In fact, many distributed algorithms such as (Ishii et al., 2012; Nedic & Ozdaglar, 2009) are often assumed to know the network size $N$. The insight of the proposed algorithm is to make its blended dynamics converge to $N$. For example, if the blended dynamics is designed as the following scalar dynamics

such that it has the stable equilibrium point at $N$, then each state $x_{i}[t]$ will also approach to $N$ under the Metropolis-Hastings coupling. Thus, by increasing $K$ until the synchronization error $\epsilon$ in (10) or (19) gets smaller than 0.5, each agent can find the exact network size through the round-off to the nearest integer.

One idea to design the heterogeneous dynamics $f_{i}$s whose average becomes (22) comes from

this is, one agent has $f_{i}=1$ and all others have $f_{i}(s)=s+1$. For this, we intentionally add one specific node which does not leave the network during the operation of the algorithm. Without loss of generality, let an index of this node be 1 and it runs the following dynamics:

On the other hand, all the other nodes of $i=2,\ldots,N$ run the following dynamics:

Note that, even though the individual dynamics of $(N-1)$ nodes for $i=2,\ldots,N$ are (marginally) unstable, the overall networked system becomes stable and the trajectories of individual agent approach close to $N$ (less than the distance of 0.5 with sufficiently large $K$). In addition, this distributed algorithm can be applied even when some agents might join or leave the network during the process of the algorithm, because it does not rely on the initial condition of agents. This idea is motivated by (Lee et al., 2018) which proposed continuous-time distributed network size estimation algorithm.

Now, we turn our attention to a different type of the weight matrix whose column-sums are all one. In this subsection, we consider the multi-step coupling framework whose coupling dynamics is the following iterative power method of PageRank (Brin & Page, 1998):

where $x_{i}\in{\mathbb{R}}$ is the state, the parameter $m\in(0,1)$ is typically chosen as 0.15, ${\mathcal{N}}_{i}$ is the in-neighbors of node $i$, and $d_{j}^{\text{out}}$ is the out-degree of node $j$. In fact, PageRank score provides an information on relative importance of each node in the network, so it has been widely utilized in diverse areas such as informatics (Chen et al., 2007), bibliometrics (Liu et al., 2005), and biology (Zaki et al., 2012).

Then, the coupling weight $w_{ij}^{\text{PR}}$ is given by

where $a_{ij}$ is the $ij$-th element of the binary adjacency matrix $A$. It can be easily seen that $w_{ij}^{\text{PR}}$ satisfies (7) and its weight matrix $W^{\text{PR}}:=[w_{ij}^{\text{PR}}]$ is obtained as

where $D_{\text{out}}$ is a diagonal matrix whose diagonal components are $d_{1}^{\text{out}},\cdots,d_{N}^{\text{out}}$ in sequence. Since $W^{\text{PR}}$ is column-stochastic, it has the spectral radius of 1 with the left eigenvector $q=\mathmybb{1}_{N}$. By Lemma 2.1, there exists a positive right eigenvector $p\in{\mathbb{R}}^{N}$ for the eigenvalue 1 such that

Here, each element $p_{i}$ of $p$ is called as PageRank score of node $i$ and it represents the relative importance of each node in the network (Brin & Page, 1998). From this, the blended dynamics under PageRank coupling (23) is written as, with $s\in{\mathbb{R}}$,

By the results in Section 2, the $i$-th agent’s trajectory over the integer count $t$ is approximated by $p_{i}s[t]$, which PageRank-scaled solution $s$ of the blended dynamics (24). Therefore, if one is interested in solving the PageRank score of each node, then the blended dynamics can be designed to have a stable equilibrium of 1.

In fact, when the network has a large number of agents, these PageRank scores are not easy to be computed in a centralized manner. Thus, the distributed PageRank algorithms have been proposed in (Ishii & Tempo, 2010; Ishii et al., 2012; Lei & Chen, 2014). Unfortunately, most of them commonly assume an initialization. However, when nodes are added to or removed from the network during the process of the algorithm, the whole algorithm must be re-initialized whenever a change occurs in the network, and this is not easy to be achieved in a distributed manner.

On the contrary, we can design an initialization-free distributed PageRank estimation algorithm by employing the proposed multi-step coupling framework. It can be easily inferred that, if the solution of the blended dynamics simply converges to 1, then every sampled state $x_{i}[t]$ will approach to its PageRank score $p_{i}$ under the multi-step coupling of (23). Thus, we first design the blended dynamics which has a stable equilibrium point at 1 as

where $\nu\in(0,1)$ is a design parameter. Since $\nu s[t]+(1-\nu)=\sum_{i=1}^{N}\left\{\nu p_{i}s[t]+(1-\nu)/N\right\}$, we can divide (25) to each node by proposing the following algorithm

Indeed, the proposed algorithm has the blended dynamics of (25). As stated in Theorem 1 or Corollary 1, the proposed distributed algorithm does not rely on a particular initialization. The algorithm (26) uses a global information of $N$, but it can be distributively estimated by the result in Section 3.1.

Average consensus protocol has been widely utilized in many discrete-time consensus problems including (Saber et al., 2007; Ren & Beard, 2005). Thus, in this subsection, we consider a multi-step coupling framework whose coupling dynamics (5b) is the following average consensus protocol

where $\theta\in(0,1)$ is a parameter which represents weight between its own state and the average of the neighbors.

Then, the coupling weight $w_{ij}^{\text{avg}}$ is obtained as

and the weight matrix $W^{\text{avg}}:=[w_{ij}^{\text{avg}}]$ is given by

where $D$ is the diagonal matrix whose diagonal components are $d_{1},\ldots,d_{N}$ sequentially. Note that $w_{ij}^{\text{avg}}$ satisfies (7) and Assumption 1 holds because $W^{\text{avg}}$ is row-stochastic matrix. Thus, we can choose $p=\mathmybb{1}_{N}$ and find a positive vector $q$ by Lemma 2.1 such that

Now, the blended dynamics is given by

Meanwhile, if the network under consideration is undirected, $q$ is easily obtained as $q=(1/d_{\text{sum}})d$ for $p=\mathmybb{1}_{N}$ where $d_{\text{sum}}:=\sum_{j=1}^{N}d_{j}$ and $d=[d_{1},\ldots,d_{N}]^{\top}\in{\mathbb{R}}^{N}$ because $d^{\top}D^{-1}A=\mathmybb{1}_{N}^{\top}A=[d_{1}^{\text{out}},\ldots,d_{N}^{\text{out}}]=[d_{1},\ldots,d_{N}]=d^{\top}$. From this, the blended dynamics (27) is rewritten as

Similarly with Section 3.1, the overall trajectories of $x_{i}[t]$ are approximately synchronized to the solution of blended dynamics (27) (or (28)) for sufficiently large $K$. The difference between the doubly-stochastic and row-stochastic weight matrix is that the former has the blended dynamics as the simple average of $f_{i}$s like (21), while the latter has the blended dynamics as the weighted average like (27) whose weights $q_{i}$s are uneven in general.

For undirected graphs, a non-increasing sequence of all degrees is called as degree sequence (Diestel, 2017). Since the degree sequence does not uniquely identify a graph, there has been much attention to obtain information of the graph structure from the given degree sequence. For example, Viger & Latapy (2005) realized the given degree sequence by a simple graph (realization problem) and Harary & Palmer (2014) estimated the number of graphs with the given degree sequence (graph enumeration). If each agent can predict possible structures of the network with the degree sequence, it could obtain global information such as the algebraic connectivity. To achieve this, a distributed algorithm to estimate the degree sequence is required, and we proposed one implementation by employing the proposed framework as follows.

In the proposed algorithm, we additionally assume that each agent knows the network size $N$ and has its unique id. Again, $N$ can be distributively estimated using the application example in Section 3.1. Moreover, the assumption on the unique id for each agent is quite natural in the sense that, in practice, every communication device has its own identifier such as mac address. Let ${\mathcal{X}}_{i}\in{\mathbb{Z}}$ be the id of the agent $i$ and suppose $1<{\mathcal{X}}_{i}$ for all $i\in{\mathcal{N}}$ without loss of generality.

Under the above assumptions, an arbitrary $i$-th agent runs the following dynamics