Jeong Woo Kim

Jin Gyu Lee

Donggil Lee

Hyungbo Shim

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

Keywords: discrete-time heterogeneous multi-agent system, multi-step coupling, blended dynamics

We develop a discrete-time version of the blended dynamics theorem for the use of designing distributed computation algorithms. The blended dynamics theorem enables to predict the behavior of heterogeneous multi-agent systems. Therefore, once we get a blended dynamics for a particular computational task, design idea of node dynamics for individual heterogeneous agents can easily occur. In the continuous-time case, prediction by blended dynamics was enabled by high coupling gain among neighboring agents. In the discrete-time case, we propose an equivalent action, which we call multi-step coupling in this paper. Compared to the continuous-time case, the blended dynamics can have more variety depending on the coupling matrix. This benefit is demonstrated with three applications; distributed estimation of network size, distributed computation of the PageRank, and distributed computation of the degree sequence of a graph, which correspond to the coupling by doubly-stochastic, columnstochastic, and row-stochastic matrices, respectively.

Introduction

Over the past decades, many distributed algorithms have been actively studied for their benefits. The benefits include lessened computational burden of one node as the burden is distributed over many nodes in the network, improved reliability against faults as a fault on one node can be compensated by redundancy of many nodes, and preserved privacy as private information need not be transferred to a central node for computation. Examples that enjoy the aforementioned benefits include distributed optimization [START_REF] Nedić | Distributed subgradient methods for multi-agent optimization[END_REF][START_REF] Nedić | Distributed optimization for control[END_REF] and distributed computation of PageRank [START_REF] Ishii | Distributed randomized algorithms for the PageRank computation[END_REF][START_REF] Suzuki | Distributed randomized algorithms for pagerank based on a novel interpretation[END_REF].

⋆ This work was supported by the National Research Foundation of Korea(NRF) grant funded by the Korea government(Ministry of Science and ICT) (No. RS-2022-00165417).

Email addresses: jeongwoo.kim@hyundai.com (Jeong Woo Kim), jin-gyu.lee@inria.fr (Jin Gyu Lee), dglee@kist.re.kr (Donggil Lee), hshim@snu.ac.kr (Hyungbo Shim).

On the other hand, constructive design methods for general distributed algorithms, except the distributed optimization, are not well developed yet. One potential approach towards the constructive design is the blended dynamics approach (Lee & Shim, 2020), which is inspired by [START_REF] Kim | Robustness of synchronization of heterogeneous agents by strong coupling and a large number of agents[END_REF][START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF]. This approach is based on the blended dynamics theorem, which, briefly speaking, asserts the following. Consider a multi-agent system ẋi = f i (t,

x i ) + κ j∈Ni (x j -x i ), i ∈ N , (1)
where x i ∈ R are the states, f i : R × R → R are heterogeneous vector fields, N := {1, • • • , N } is the set of node indices, and N i is the index set of nodes that send information to node i. The individual node dynamics are represented by ẋi = f i (t, x i ), which is coupled with neighboring nodes by the coupling term j∈Ni (x j -x i ) with a common coupling gain κ. Then, under the assumption that the communication graph is undirected and connected, every agent in (1) behaves like the blended dy-namics defined as

ṡ(t) = 1 N N i=1 f i (t, s(t)) ∈ R (2)
where s is the state of the blended dynamics, if the coupling gain κ is sufficiently large and if the blended dynamics are contractive, i.e., incrementally stable [START_REF] Kim | Robustness of synchronization of heterogeneous agents by strong coupling and a large number of agents[END_REF]. More precisely, for any ϵ > 0, there exists κ min such that, if κ > κ min , lim sup

t→∞ |x i (t) -s(t)| ≤ ϵ, ∀i ∈ N . (3)
Since the blended dynamics are a simple average of individual node dynamics, it has been utilized as a design tool for many distributed algorithms; that is, one designs a desired algorithm as the blended dynamics (2) first, and then, splits it into different node dynamics (1) with couplings. This philosophy has been successfully employed in many applications such as distributed economic power dispatch problem (Yun, Shim, & Ahn, 2018), distributed state estimator [START_REF] Kim | Completely decentralized design of distributed observer for linear systems[END_REF], secure estimation by distributed median computation [START_REF] Lee | Fully distributed resilient state estimation based on distributed median solver[END_REF], distributed optimization without convexity of each node [START_REF] Lee | Blended dynamics approach to distributed optimization: Sum convexity and convergence rate[END_REF], decentralized controller design [START_REF] Kim | Decentralized design and plug-and-play distributed control for linear multi-channel systems[END_REF], and analysis of coupled oscillators [START_REF] Trummel | On optimal synchronization of diffusively coupled heterogeneous van der pol oscillators[END_REF]. See [START_REF] Lee | Design of heterogeneous multi-agent system for distributed computation[END_REF] for more comprehensive summary of these applications.

The distributed algorithms designed by the blended dynamics theorem does not require each node dynamics to be stable, as long as their average (i.e., the blended dynamics) is contractive, which yields flexibility of the design. Moreover, as long as the blended dynamics remain contractive, a new node can join the network or an existing node can leave the network during the operation, which is called as a plug-and-play feature. This is because the designed distributed algorithms are initialization-free, i.e., the proposed algorithms do not rely on any specific initial values (Lee & Shim, 2020).

While all the above results are in the continuous-time domain, it is however required to implement the designed algorithm in the discrete-time domain so that it operates on digital devices in practice. A naive idea is to use simple discretization methods such as forward difference.

For example, a discretized model of (1) becomes

x i (t + ∆ t ) = x i (t) + ∆ t f i (t, x i (t)) + κ∆ t j∈Ni (x j (t) -x i (t)), ( 4
)
where ∆ t is the sampling time. In the continuous-time case, we recall that arbitrarily small error between x i and s in (3) can be obtained by increasing the coupling gain κ, so that an emergent behavior of the multi-agent system arises with strong coupling. However, in the discretetime case of (4), increasing κ unboundedly yields instability of the network unless ∆ t is decreased with the same ratio 1 . Therefore, the discrete-time algorithm (4) is not suitable for a discrete-time version of the blended dynamics approach.

In this paper, we propose a new form of a multi-agent system (which is given by ( 5) in Section 2). We note that the meaning of using a large coupling gain κ in the continuous-time case of ( 1) is that consensus is taken more care of than the progress through the node dynamics. Based on the observation, and motivated by (Wang, Liu, Morse, & Anderson, 2019), the proposed form repeats a weighted averaging action many times before progressing through the node dynamics, which we call 'multi-step coupling.'

This approach still maintains the advantages of the continuous-time case, such as the initialization-free benefit or the plug-and-play operation, and that the individual node dynamics need not be stable as long as the blended dynamics are stable. The former benefit is useful for a distributed algorithm because, without it, the algorithm must be re-initialized whenever a change occurs in the network such as addition of new nodes or removal of existing nodes, and immediate detection of a change in a distant node is not easy for a large scale network. In particular, the initialization-free feature is beneficial for a distributed estimation of relative importance of each node in the network, i.e., the PageRank score. While there have been various studies for the distributed estimation of the PageRank scores such as [START_REF] Ishii | Distributed randomized algorithms for the PageRank computation[END_REF][START_REF] Lei | Distributed randomized PageRank algorithm based on stochastic approximation[END_REF][START_REF] Suzuki | Distributed randomized algorithms for pagerank based on a novel interpretation[END_REF], they assume the initialization process. On the contrary, the proposed approach provides a design method for initialization-free distributed algorithms. An initialization-free algorithm for PageRank will appear in Section 3.2. Now, the latter benefit that not all node dynamics have to be stable, has been successfully utilized in a few applications. Distributed state estimator in the continuous-time [START_REF] Kim | Completely decentralized design of distributed observer for linear systems[END_REF] or in the discrete-time (Wang et al., 2019) is one example because individual observers at each agent can not be stable due to lack of observability. But, with strong couplings in some sense, their collective behavior becomes a stable one. Another example can be found in Section 3.1, where individual node dynamics are not stable, but their blended dynamics are stable so that they can estimate the size of the network asymptotically.

Moreover, while the continuous-time approach predicted collective synchronization behavior of the multi-agent system in [START_REF] Kim | Robustness of synchronization of heterogeneous agents by strong coupling and a large number of agents[END_REF]Lee & Shim, 2020), this discrete-time approach estimates not only emergent but also individually scaled behavior, i.e., each agent behaves similarly to the solution of the blended dynamics with an agent-wise scaling factor. For example, in Section 3.2, we will introduce an application example where each node estimates its relative importance which is possibly agentwise different so that overall nodes are not synchronized in the network.

Notation and useful facts

We use the following convention in this paper. 0 N and 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 2norm of a vector and the induced 2-norm of a matrix are denoted by ∥ • ∥. For any square matrix M, ρ(M) denotes the spectral radius of M. The operation defined by the symbol ⊗ is Kronecker product, which has the properties that, given the matrices A, B, C and

D of ap- propriate dimensions, (A ⊗ B)(C ⊗ D) = (AC) ⊗ (BD)
and ∥A ⊗ B∥ = ∥A∥∥B∥. See, e.g., [START_REF] Bernstein | Matrix Mathematics[END_REF] for more about the Kronecker product. For notational convenience, we use the convention M ⊗n := M ⊗ I n for any size of matrix M.

A communication network of N agents can be represented by a directed graph G = (N , E) where N = {1, 2, . . . , N } is the node set, and E ⊂ N × 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) ∈ E. A graph is said to be undirected if (j, i) ∈ E implies (i, j) ∈ E. For a directed graph, in-neighbors of node i are represented by

N i = {j ∈ N |(j, i) ∈ E}
:= |N out i |.
A path of length L from node i to node j is a sequence (i 0 , i 1 , . . . , i L ) such that i 0 = i, i L = j, (i l , i l+1 ) ∈ E for any l ∈ {0, . . . , 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 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 G, its adjacency matrix

A = [α ij ] ∈ R N ×N is defined as α ij > 0 if (j, i) ∈ E and α ij = 0 otherwise. If every positive component of A is 1, then A is
called a binary adjacency matrix. Meanwhile, given a non-negative square matrix M ∈ R N ×N , its associated directed graph is defined as the directed graph whose adjacency matrix is M. A non-negative square matrix M ∈ R N ×N is said to be primitive if there exists n ∈ N such that M n is positive, i.e., every component of M n is positive. The primitive property of any nonnegative square matrix can be verified by its associated graph as the following lemma.

Lemma 1 (Bullo (2020, Theorem 4.7)) For a directed graph G, its adjacency matrix A, and its binary adjacency matrix A, the following statements are equivalent: (a) G is strongly connected and aperiodic, (b) A is primitive, and (c) A is primitive.

Lemma 2 (Perron-Frobenius theorem) If a nonnegative matrix M ∈ R N ×N is primitive, then the following statements hold:

1) There exists a simple eigenvalue λ > 0 (called Perron-Frobenius eigenvalue) of M such that λ > |σ| for any other eigenvalue σ of M.

2) The right and left eigenvectors of λ are positive.

Discrete-time Blended Dynamics Theorem

For a discrete-time version of the blended dynamics theorem, we propose the following discrete-time algorithm: for each agent i ∈ N ,

x i [t k+1 ] =    f i (t k , x i [t k ]), if k = 0, (5a) j∈Ni∪{i} w ij x j [t k ], if k = 1, . . . , K -1, (5b)
where x i ∈ R n is the state, the function f i : Z×R n → R n is continuously differentiable and represents the timevarying heterogeneous node dynamics (5a), and the coefficient w ij , called coupling weight, determines the behavior of the coupling dynamics (5b). Here, t k is the symbol defined by

t k = t + k K , (6)
where K ∈ N, and we call t k by fractional discrete-time index. In particular, we call t ∈ Z as integer count and k ∈ Z as fraction count. The fraction count k varies from 0 to K -1. Keeping in mind that t K = t + K/K = (t + 1) + 0/K = (t + 1) 0 , we see that the fractional discrete-time t k advances as 0

0 , 0 1 , • • • , 0 K-1 , 1 0 , 1 1 , • • • .
The time t 0 will often be written as t for convenience. The fractional discrete-time has nothing to do with real time, and can be implemented in practice just as a sequential order in an algorithm.

We will choose K sufficiently large, which determines how many times the coupling dynamics (5b) are executed before the next node dynamics (5a) are executed.

It will be shown that, in this way, the effect of strong coupling κ in the continuous-time blended dynamics theorem can be similarly reflected in discrete-time. To emphasize the difference, we call this type of coupling in (5) by multi-step coupling.

The coupling weights w ij in (5b) have the property:

w ij > 0, j ∈ N i ∪ {i}, = 0, otherwise. (7)
Now, we assume that the matrix W := [w ij ] ∈ R N ×N , which we call a weight matrix, satisfies the following.

Assumption 1 The spectral radius of W is 1.

Note that the communication protocols in many discrete-time multi-agent systems in the literature have the form of linear combination like in (5b) and their weight matrices satisfy Assumption 1. Examples include [START_REF] Ishii | Distributed randomized algorithms for the PageRank computation[END_REF][START_REF] Lei | Distributed randomized PageRank algorithm based on stochastic approximation[END_REF][START_REF] Saber | Consensus and cooperation in networked multi-agent systems[END_REF][START_REF] Ren | Consensus seeking in multiagent systems under dynamically changing interaction topologies[END_REF], in which the weight matrices are given by stochastic matrices whose spectral radius is 1.

Meanwhile, the communication network under consideration is represented by the directed graph G, which does not have self-connection at any node by definition, and we assume the following.

Assumption 2 The network G is strongly connected.

Then, under Assumptions 1 and 2, the following is wellknown (but we put its proof for readers' convenience).

Lemma 3 Let λ i , i ∈ N , be the eigenvalues of W such that |λ 1 | ≥ |λ 2 | ≥ . . . ≥ |λ N |.
Under Assumptions 1 and 2, λ 1 = 1, λ 1 > |λ j | for all j = 2, . . . , N , and there exist positive vectors p, q ∈ R N such that

W p = p, q ⊤ W = q ⊤ , q ⊤ p = 1. (8)
PROOF. From ( 7), the associated graph of W not only contains all edges of G but also has a self-connection for every node because all diagonal entries of W are positive. Thus, the associated graph is aperiodic as well as strongly connected. This implies that W is primitive by Lemma 1, and, by Assumption 1 and Lemma 2, W has the simple Perron-Frobenius eigenvalue 1, i.e., λ 1 = 1 and λ 1 > |λ j | for all j = 2, . . . , N , with positive right and left eigenvectors p and q, respectively. Scaling p and q yields that q ⊤ p = 1. ■

With p and q from Lemma 3 at hand, we now introduce discrete-time blended dynamics, which are defined as a weighted average of node dynamics:

s[t + 1] = N i=1 q i f i (t, p i s[t]) =: f s (t, s[t]) ∈ R n (9)
where p i and q i is the i-th element of p and q, respectively, and t is the integer count of the fractional time (i.e., t = t 0 ). In particular, we assume the blended dynamics are stable in the sense of [START_REF] Lohmiller | On contraction analysis for non-linear systems[END_REF][START_REF] Tran | Convergence properties for discrete-time nonlinear systems[END_REF] as follows.

Assumption 3 The blended dynamics (9) are contractive; i.e., there exist a (symmetric) positive definite matrix H ∈ R n×n and a positive constant γ < 1 such that

∂f s ∂s (t, s) ⊤ H 2 ∂f s ∂s (t, s) ≤ γH 2 , ∀s ∈ R n , t ∈ Z.
Remark 1 Assumption 3 does not ask each node dynamics

x i [t + 1] = f i (t, x i [t]
) to be stable. Rather it allows unstable node dynamics whose instability can be compensated by other node dynamics so that the blended dynamics become stable in the sense of Assumption 3. For example, when there are four agents with

f 1 (t, x) = f 2 (t, x) = 0.1x and f 3 (t, x) = f 4 (t, x) = 1.5x
, the agents 1 and 2 have stable node dynamics while the agents 3 and 4 have unstable ones. If the weight matrix has the vectors p = 1 4 and q = 1 4 /4, then Assumption 3 holds because f s (s) = 0.8s.

We will see that the blended dynamics (9) allow to predict the behavior of (5) when K is large. To make the prediction effective from any initial conditions globally in the state-space, we impose the following assumption.

Assumption 4

The function f i (t, x) is uniformly bounded in t and globally Lipschitz with respect to x uniformly in t: i.e., ∃ a non-decreasing continuous function M : R → R and a constant L ≥ 0 such that, ∀x, y ∈ R n , t ∈ Z, and i ∈ N ,

∥f i (t, x)∥ ≤ M (∥x∥) , ∥f i (t, x) -f i (t, y)∥ ≤ L∥x -y∥.
Theorem 1 Under Assumptions 1-4, for any ϵ > 0, there exists K min such that, for all K > K min , the solution x i of (5) and the solution s of (9) with arbitrary initial conditions satisfy

lim sup t→∞ x i [t] -p i s[t] ≤ ϵ, ∀i ∈ N . ( 10
)
In addition, for each k ∈ {1, 2, . . . , K -1} and i ∈ N ,

lim sup t→∞ x i [t k ] -p i s[t + 1] ≤ ϵ 2 1 + 1 |λ N | K-k . (11)
Theorem 1 states that, with sufficiently large number of steps for the coupling (5b), the behavior of node dynamics (5a), which is represented by the state x i at the integer count t, can be approximately predicted by the solution s of the blended dynamics with the scaling factor p i , and the approximation error can be made arbitrarily small by increasing K. Moreover, the behavior of x i over the fraction counts is also bounded with respect to the scaled trajectory of s.

Remark 2 From proof of Theorem 1 in Appendix B, K min can be explicitly defined as the smallest integer such that

K min ≥ log |λ2| 1 - √ γ 2 M1 max 1 L , ϵ 2M (∥p∥M s ) √ N ,
where M1 and M s are some constants which will be induced from (B.5) and (B.6). This expression yields a reasonable interpretation in the sense that K min increases as the second largest eigenvalue λ 2 of the weight matrix W approaches to 1, the performance index ϵ decreases, or the Lipschitz constant L and the network size N get larger, while it decreases as the degree of stability of the blended dynamics, (1 -√ γ), gets larger. Note that the network topology implicitly affects K min through λ 2 . In practice, a suitable K is chosen by repeated simulations, or by a maximal value if an upper bound of the right-hand side for all possible cases can be computed (like in [START_REF] Lee | Distributed algorithm for the network size estimation: Blended dynamics approach[END_REF]).

Remark 3 Selection of the eigenvectors p and q as (8) is not unique, but the result of Theorem 1 remains the same. To see this, we note that different selection of p ′ and q ′ from p and q, respectively, should satisfy p ′ = cp and q ′ = (1/c)q for some c > 0 because of the Perron-Frobenius theorem (Lemma 2). In addition, we note that the new blended dynamics become

s ′ [t + 1] = (1/c) N i=1 q i f i (t, cp i s ′ [t]) =: f s ′ (t, s ′ [t]). Comparing it with (9), it is seen that s ′ [t] = (1/c)s[t],
and thus, we have ∥x

i [t] -p ′ i s ′ [t]∥ = ∥x i [t] -cp i (1/c)s[t]∥ = ∥x i [t] -p i s[t]∥.
Finally, it is also seen that the new blended dynamics satisfy Assumption 3 because

(∂f s ′ /∂s ′ )(t, s ′ ) = (∂f s /∂s)(t, s) with s = cs ′ . Remark 4 In (11), the state x i [t k ] is compared not with s[t] but with s[t + 1].
One may find this is natural considering the behavior of the overall system. At each integer time t = t 0 , each x i obeys the heterogeneous node dynamics (5a), which potentially updates

x i [t 1 ] in different directions from the updated p i s[t + 1] (even if x i [t 0 ] is close to p i s[t]). Instead, repeated execution of (5b) drives x i [t k ] to p i s[t + 1], which is well reflected in (11).
We now present intuitive explanations for Theorem 1, whose rigorous proof continues in the Appendix. For simplicity, define x :=

[x ⊤ 1 , . . . , x ⊤ N ] ⊤ ∈ R nN . Then, (5) is simply written as x[t k ] = W k-1 ⊗n      f 1 (t 0 , x 1 [t 0 ]) . . . f N (t 0 , x N [t 0 ])      =: W k-1 ⊗n F (t 0 , x[t 0 ]),
(12) for k = 1, . . . , K -1 and t ∈ Z, where W ⊗n = W ⊗ I n , and, since t K = (t + 1) 0 = t + 1, we have

x[t + 1] = W K-1 ⊗n F (t, x[t]). ( 13
)
Similar to (13), the blended dynamics ( 9) are written as

s[t + 1] = N i=1 q i f i (t, p i s[t]) = q ⊤ ⊗n F (t, p ⊗n s[t]). ( 14
)
By Lemma 3, there exist R, Z ∈ R N ×(N -1) such that

W = p R 1 0 0 Λ q ⊤ Z ⊤ , ( 15
)
Z ⊤ R = I N -1 , and Z ⊤ p = R ⊤ q = 0 N -1 , (16)
where Λ ∈ R (N -1)×(N -1) is a matrix whose eigenvalues are λ 2 , . . . , λ N . With them, we consider the coordinate transformation

ξ = ξ 1 ξ = q ⊤ ⊗n Z ⊤ ⊗n x ( 17
)
whose inverse is

x = p ⊗n ξ 1 + R ⊗n ξ.
The overall dynamics (13) at each integer time t become

ξ 1 [t + 1] = q ⊤ ⊗n W K-1 ⊗n F (t, p ⊗n ξ 1 [t] + R ⊗n ξ[t]) = q ⊤ ⊗n F (t, p ⊗n ξ 1 [t] + R ⊗n ξ[t]), ξ[t + 1] = Z ⊤ ⊗n W K-1 ⊗n F (t, p ⊗n ξ 1 [t] + R ⊗n ξ[t]) = (Λ K-1 Z ⊤ ) ⊗n F (t, p ⊗n ξ 1 [t] + R ⊗n ξ[t]). (18)
Define the error variable e := ξ 1 -s. Then the above dynamics are rewritten by

e[t + 1] = q ⊤ ⊗n F (t, p ⊗n (e[t] + s[t]) + R ⊗n ξ[t]) -q ⊤ ⊗n F (t, p ⊗n s[t]), ξ[t + 1] = (Λ K-1 Z ⊤ ) ⊗n F (t, p ⊗n (e[t] + s[t]) + R ⊗n ξ[t]).
Since the spectral radius ρ(Λ) = |λ 2 | < 1, it may be inferred that ∥ ξ∥ gets small if K is sufficiently large. On the other hand, if ξ happens to be identically zero, then e[t] converges to zero as t tends to infinity, which follows from the following lemma.

Lemma 4 Under Assumption 3,

∥H{f s (t, s 2 ) -f s (t, s 1 )}∥ ≤ √ γ∥H(s 2 -s 1 )∥ for all t ∈ Z and s 1 , s 2 ∈ R n .
In fact, if ξ ≡ 0 n(N -1) , then,

He[t + 1] = ∥H{q ⊤ ⊗n F (t, p ⊗n (e[t] + s[t])) -q ⊤ ⊗n F (t, p ⊗n s[t])}∥ = ∥H{f s (t, e[t] + s[t]) -f s (t, s[t])}∥ ≤ √ γ∥He[t]∥.
This implies that ξ 1 and s tend to get closer as time goes on. When this happens, x = p ⊗n ξ 1 + R ⊗n ξ tends to p ⊗n s, which motivates Theorem 1.

However, ξ does not become identically zero in general, and the above arguments need rigorous analysis, which is done in the Appendix with a Lyapunov function.

In Theorem 1, the approximation of x i [t] by p i s[t] is stated in an asymptotic format, i.e., using lim sup. If one questions when the approximation becomes effective, the following corollary assures that it can be very early if K is sufficiently large. In particular, if the initial conditions of x i are in a known compact set, then K min can be computed (see the proof in the Appendix).

Corollary 1 Under Assumptions 1-4, for any ϵ > 0 and compact set C ⊂ R n , there exists K min > 0 such that, for all K > K min , the solution x i of (5) with x i [0] ∈ C, and the solution s of (9) with s

[1] = N i=1 q i f i (0, x i [0]) satisfy x i [t] -p i s[t] ≤ ϵ, ∀t ≥ 1, i ∈ N . ( 19
)
In addition, for all k ∈ {1, 2, . . . , K -1}, i ∈ N , and t ≥ 1,

x i [t k ] -p i s[t + 1] ≤ ϵ 2 1 + 1 |λ N | K-k . ( 20
)
Remark 5 In Corollary 1, the solution s[t] of the blended dynamics is initiated not at t = 0 but at t = 1. One may find this is natural because the state

s ′ [1] = f s (0, s ′ [0]) with s ′ [0] = N i=1 q i x i [0] (i.e., s ′ [1] = N i=1 q i f i (0, p i N i=1 q i x i [0])) can be very differ- ent from the considered state s[1] = N i=1 q i f i (0, x i [0]).
In fact, the states x j [0 1 ], j ∈ N , may be very different from each other, but they converge with sufficiently large K towards N i=1 q i x i [0 1 ] with the scaling factor p j , which is s[1] (not s ′ [1]).

Network Synthesis with Examples

As mentioned before, the proposed approach is useful as a design method for distributed algorithms by first designing suitable blended dynamics such that they behave as desired and then synthesizing each heterogeneity of the multi-agent system such that it has the predesigned blended dynamics. Moreover, comparing with the continuous-time approach, the discrete-time version handles more general protocols as long as the spectral radius of its weight matrix is 1. In fact, many studies on discrete-time multi-agent system including PageRank [START_REF] Ishii | Distributed randomized algorithms for the PageRank computation[END_REF][START_REF] Lei | Distributed randomized PageRank algorithm based on stochastic approximation[END_REF] or consensus [START_REF] Saber | Consensus and cooperation in networked multi-agent systems[END_REF][START_REF] Ren | Consensus seeking in multiagent systems under dynamically changing interaction topologies[END_REF]) have used the communication protocol which can be represented as a linear combination of agents' state information like (5b) with the weight matrix of unit spectral radius. Thus, in this section, some of the protocols are chosen to be considered as the coupling dynamics (5b) and the behaviors of corresponding multi-step coupling systems are illustrated using the results in previous section. Based on this, the design process for each application example is also provided.

Distributed Network Size Estimation with Metropolis-Hastings Coupling

In this subsection, we assume that the network is undirected and connected, and consider the following Metropolis-Hastings coupling weight w MH ij in [START_REF] Schwarz | On the convergence of average consensus with generalized Metropolis-Hasting weights[END_REF]:

w MH ij :=            1 -µ max {d i , d j } , (j, i) ∈ E and i ̸ = j, 0, (j, i) / ∈ E and i ̸ = j, 1 - l̸ =i w MH il , i = j,
where µ ∈ (0, 1). Note that w MH ij 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 MH ij satisfies (7) and Assumption 1 holds for W MH := [w MH ij ] because it is doubly-stochastic. Hence, we can choose p = 1 N and q = (1/N )1 N from Lemma 3. Then, the blended dynamics are obtained as a simple average of f i s as follows:

s[t + 1] = 1 N N i=1 f i (t, s[t]). (21)
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 = 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 [START_REF] Ishii | Distributed randomized algorithms for the PageRank computation[END_REF][START_REF] Nedić | Distributed subgradient methods for multi-agent optimization[END_REF] 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 are designed as the following scalar dynamics

s[t + 1] = 1 - 1 N s[t] + 1, (22)
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 ϵ 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

1 N 1 + N i=2 (s[t] + 1) = 1 - 1 N s[t] + 1,
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:

x 1 [t k+1 ] = 1, if k = 0, j∈N1∪{1} w MH 1j x j [t k ], otherwise. (23)
On the other hand, all the other nodes of i = 2, . . . , N run the following dynamics:

x i [t k+1 ] = x i [t k ] + 1, if k = 0, j∈Ni∪{i} w MH ij x j [t k ], otherwise. (24)
Note that, even though the individual dynamics of (N -1) nodes for i = 2, . . . , 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 [START_REF] Lee | Distributed algorithm for the network size estimation: Blended dynamics approach[END_REF] which proposed continuous-time distributed network size estimation algorithm.

Initialization-free Distributed PageRank Estimation with PageRank Coupling

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 are the following iterative power method of PageRank [START_REF] Brin | The anatomy of a large-scale hypertextual web search engine[END_REF]:

x i [t k+1 ] = mx i [t k ] + (1 -m) j∈Ni x j [t k ] d out j , (25)
where x i ∈ R is the state, the parameter m ∈ (0, 1) is typically chosen as 0.15, N i is the in-neighbors of node i, and d out j 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 [START_REF] Chen | Finding scientific gems with Google PageRank algorithm[END_REF], bibliometrics [START_REF] Liu | Co-authorship networks in the digital library research community[END_REF], and biology (Zaki, Berengueres, & Efimov, 2012).

Then, the coupling weight w PR ij is given by

w PR ij =    m, i = j, (1 -m) a ij d out j , i ̸ = j,
where a ij is the ij-th element of the binary adjacency matrix A. It can be easily seen that w PR ij satisfies (7) and its weight matrix W PR := [w PR ij ] is obtained as

W PR = mI + (1 -m)AD -1 out ,
where D out is a diagonal matrix whose diagonal components are d out 1 , • • • , d out N in sequence. Since W PR is column-stochastic, it has the spectral radius of 1 with the left eigenvector q = 1 N . By Lemma 3, there exists a positive right eigenvector p ∈ R N for the eigenvalue 1 such that

W PR p = p, 1 ⊤ N p = 1.
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 [START_REF] Brin | The anatomy of a large-scale hypertextual web search engine[END_REF]. From this, the blended dynamics under PageRank coupling (25) are written as, with s ∈ R,

s[t + 1] = 1 ⊤ N F (t, ps[t]) = N i=1 f i (t, p i s[t]) . (26)
By the results in Section 2, the i-th agent's trajectory over the integer count t is approximated by p i s[t], i.e., PageRank-scaled solution s of the blended dynamics (26). 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 [START_REF] Ishii | Distributed randomized algorithms for the PageRank computation[END_REF][START_REF] Lei | Distributed randomized PageRank algorithm based on stochastic approximation[END_REF][START_REF] Suzuki | Distributed randomized algorithms for pagerank based on a novel interpretation[END_REF]. 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 multistep coupling of ( 25). Thus, we first design the blended dynamics which have a stable equilibrium point at 1 as

s[t + 1] = νs[t] + (1 -ν), (27)
where ν ∈ (0, 1) is a design parameter. Since 27) to each node by proposing the following algorithm

νs[t] + (1 - ν) = N i=1 {νp i s[t] + (1 -ν)/N }, we can divide (
x i [t k+1 ] =        νx i [t k ] + 1 -ν N , if k = 0, mx i [t k ] + (1 -m) j∈Ni x j [t k ] d out j , otherwise. (28)
Indeed, the proposed algorithm has the blended dynamics of ( 27). As stated in Theorem 1 or Corollary 1, the proposed distributed algorithm does not rely on a particular initialization. The algorithm (28) uses a global information of N , but it can be distributively estimated by the result in Section 3.1.

Distributed Degree Sequence Estimation with Average Coupling

Average consensus protocol has been widely utilized in many discrete-time consensus problems including [START_REF] Saber | Consensus and cooperation in networked multi-agent systems[END_REF][START_REF] Ren | Consensus seeking in multiagent systems under dynamically changing interaction topologies[END_REF]. Thus, in this subsection, we consider a multi-step coupling framework whose coupling dynamics (5b) are the following average con-sensus protocol

x i [t k+1 ] = θx i [t k ] + 1 -θ |N i | j∈Ni x j [t k ],
where θ ∈ (0, 1) is a parameter which represents weight between its own state and the average of the neighbors.

Then, the coupling weight w avg ij is obtained as

w avg ij = θ, i = j, (1 -θ) a ij d i , otherwise,
and the weight matrix W avg := [w avg ij ] is given by

W avg = θI + (1 -θ)D -1 A,
where D is the diagonal matrix whose diagonal components are d 1 , . . . , d N sequentially. Note that w avg ij satisfies ( 7) and Assumption 1 holds because W avg is rowstochastic matrix. Thus, we can choose p = 1 N and find a positive vector q by Lemma 3 such that q ⊤ W avg = q ⊤ , q ⊤ 1 N = 1. Now, the blended dynamics are given by

s[t + 1] = N i=1 q i f i (t, s[t]). (29)
Meanwhile, if the network under consideration is undirected, q is easily obtained as q = (1/d sum )d for p = 1 N where

d sum := N j=1 d j and d = [d 1 , . . . , d N ] ⊤ ∈ R N because d ⊤ D -1 A = 1 ⊤ N A = [d out 1 , . . . , d out N ] = [d 1 , . . . , d N ] = d ⊤ .
From this, the blended dynamics (29) are rewritten as

s[t + 1] = 1 d sum N i=1 d i f i (t, s[t]). ( 30
)
Similarly with Section 3.1, the overall trajectories of x i [t] are approximately synchronized to the solution of blended dynamics (29) (or ( 30)) 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 (29) whose weights q i s are uneven in general.

For undirected graphs, a non-increasing sequence of all degrees is called as degree sequence [START_REF] Diestel | Graph Theory[END_REF]. 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, [START_REF] Viger | Efficient and simple generation of random simple connected graphs with pre[END_REF] realized the given degree sequence by a simple graph (realization problem) and [START_REF] Harary | Graphical Enumeration[END_REF] 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 X i ∈ Z be the id of the agent i and suppose 1 < X i for all i ∈ N without loss of generality.

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

x i [t k+1 ] = (1 -1/d i ) x i [t k ] + N Xi , if k = 0, j∈Ni∪{i} w avg ij x j [t k ], otherwise,
where x i ∈ R is the state. From (30), the blended dynamics are obtained as

s[t + 1] = 1 d sum N i=1 d i 1 - 1 d i s[t] + N i=1 d i N Xi = 1 - N d sum s[t] + 1 d sum N i=1 d i N Xi .
For any connected network, d sum ≥ N , and this guarantees the contraction stability of the above blended dynamics. In addition, its equilibrium point s * is obtained as

s * = 1 N N i=1 d i N Xi = N i=1 d i N Xi-1 .
Meanwhile, it is easily seen that d i < N for all i ∈ N because the network has no self-connection. Thus, the number s * can be regarded as a representation of the numeral system with the base N : where B = max i∈N X i and δ b represents the (b + 1)-th right-most digit such that

[s * ] N = [δ (B-1) δ (B-2) • • • δ 1 δ 0 ] N , (31)
δ b = d i if b = X i -1 and δ b = 0 otherwise.
As a result, since every state x i [t] is approximately synchronized to s * of (31), each agent can estimate the degree sequence by removing zero digits in the N -base numeral representation of its state and ordering the rest digits. Here, we suppose that K is properly chosen such that the synchronization error ϵ < 1 because this guarantees that the error does not affect even in the rightmost digit δ 0 .

Simulations

In this section, we present simulation results2 to validate the algorithms proposed in Section 3 for estimating network size and PageRank scores. In particular, we consider an undirected network that changes over time and experiences the following scenarios:

1) Initially, the network has five agents.

2) At t = 100, five new agents join the network.

3) At t = 200, two agents leave the network while the others remain.

The network for each time interval is demonstrated in Fig. 1. The agents in the network can be classified into three types: red circle agents, which always remain in the network; blue rectangle agents, which join the network later; and green triangle agents, which stay for a while but eventually leave the network.

Network Size Estimation

To estimate the network size through the algorithm in Section 3.1, we first select an agent from the red circle agents and denote it as i = 1. Then, we allocate the node dynamics as f i = 1 for the selected agent and f i (x i ) = x i + 1 for all the other agents. Note that, under the Metropolis-Hastings coupling, the blended dynamics of (5) become (22). Here, K is set to be 25 for multi-step coupling and µ to be 0.1 for the Metropolis-Hastings coupling weight. Each state converges to the total number of connected agents with a deviation no greater than 0.5, regardless of agents joining or leaving from the network. This enables every agent to obtain the network size accurately. In contrast, due to initial condition dependency, the IAM and OSM algorithms' estimation values converge to the network size only when new agents join the network and they are no longer accurate when some agents leave the network as shown in Fig. 2 (b). Meanwhile, Fig. 3 demonstrates the estimation results of the proposed algorithm with different values of K. As K gets larger, the convergence error to the number of total agents decreases.

PageRank Estimation

As presented in (28), the implementation of proposed algorithm for PageRank estimation requires a global information of the network size N . We utilize the network size estimation algorithm in Section 3.1 to replace the use of N in (28), so that the algorithm becomes fully distributed.

In Fig. 4, the estimation results obtained by the proposed algorithm are compared with those obtained by other algorithms, including the Randomized Algorithm (RA) by [START_REF] Ishii | Distributed randomized algorithms for the PageRank computation[END_REF], and Synchronous Algorithm based on New Interpretation (SANI) by [START_REF] Suzuki | Distributed randomized algorithms for pagerank based on a novel interpretation[END_REF]. Fig. 4 (a) shows the estimation result obtained by the proposed algorithm, where each state trajectory of the agent at every integer count is displayed.

We set the parameter ν as 0.9 for the blended dynamics ( 27) and m as 0.2 for the coupling weight w PR ij . Despite changes in network configuration, each state of the agent i converges to its respective PageRank score p i . On the other hand, the RA and SANI algorithms' estimation values converge to the correct value at first, but they fail to handle changes in the network, as shown in Fig. 4 (b).

Conclusion

In this paper, we introduced a discrete-time blended dynamics theorem, which inherits all the benefits of the continuous-time blended dynamics theorems in (Lee & Shim, 2020). This was achieved by the proposed multistep coupling in the multi-agent algorithm (5). The proposed approach does not require stability of individual agents as long as the blended dynamics are stable, the plug-and-play operation is easily achieved. To illustrate utility of the proposed method as a design tool, three application examples are included. scribed degree sequence. In Proceedings of the International Computing and Combinatorics Conference, pp. 400-449. Berlin: Springer. Wang, L., Liu, J., Morse, A. S., & Anderson, B. D. (2019)

:= {f s (t, s 2 ) -f s (t, s 1 )} ⊤ H 2 × {f s (t, cs 2 + (1 -c)s 1 ) -f s (t, s 1 )}. (c)
Since the function g t,s1,s2 : R → R is continuously differentiable, by the mean-value theorem, there exists c ∈ [0, 1] such that g t,s1,s2 (1)-g t,s1,s2 (0) = g ′ t,s1,s2 (c)(1-0), which is equivalent to

∆(t, s 1 , s 2 ) = {f s (t, s 2 ) -f s (t, s 1 )} ⊤ H 2 × ∂f s ∂s (t, s)(s 2 -s 1 ) (A.2)
where s = cs 2 + (1 -c)s 1 . Using this, we have

∆(t, s 1 , s 2 ) = ∥H{f s (t, s 2 ) -f s (t, s 1 )}∥ 2 ≤ ∥H{f s (t, s 2 ) -f s (t, s 1 )}∥ H ∂f s ∂s (t, s)(s 2 -s 1 ) , which in turn implies ∥H{f s (t, s 2 ) -f s (t, s 1 )}∥ ≤ ∥H(∂f s /∂s)(t, s)(s 2 -s 1 )∥.
From this, the claim (A.1) is justified. Finally, it follows from Assumption 3 that

∥H{f s (t, s 2 ) -f s (t, s 1 )}∥ 2 = ∆(t, s 1 , s 2 ) ≤ (s 2 -s 1 ) ⊤ ∂f s ∂s (t, s) ⊤ H 2 ∂f s ∂s (t, s)(s 2 -s 1 ) ≤ γ(s 2 -s 1 ) ⊤ H 2 (s 2 -s 1 ) = γ∥H(s 2 -s 1 )∥ 2 ,
which proves Lemma 4.

B Proof of Theorem 1

The proof begins with the claim that the solution s[t] of the blended dynamics is bounded by Assumptions 3-4.

Lemma B.1 Under Assumption 3, the solutions of the blended dynamics (9), initiated at time t = t 0 , are bounded as

∥Hs[t]∥ ≤ √ γ t-t 0 ∥Hs[t 0 ]∥ + sup τ ≥t 0 ∥Hf s (τ, 0 n )∥ 1 - √ γ . (B.1)
Additionally, with Assumption 4, we have

lim sup t→∞ ∥s[t]∥ ≤ √ N ∥q∥∥H -1 ∥∥H∥M (0) 1 - √ γ =: M s (B.2)
PROOF. We prove (B.1) by showing that

∥Hs[t]∥ ≤ w[t], ∀t ≥ t 0 (B.3)
where w ∈ R is the solution of

w[t + 1] = √ γw[t] + ∥Hf s (t, 0 n )∥, w[t 0 ] = ∥Hs[t 0 ]∥. (B.4) Indeed, since ∥Hs[t 0 ]∥ ≤ w[t 0 ], let us suppose ∥Hs[τ ]∥ ≤ w[τ ] for some integer τ ≥ t 0 . Then, ∥Hs[τ + 1]∥ = ∥H{f s (τ, s[τ ]) -f s (τ, 0 n ) + f s (τ, 0 n )}∥ ≤ ∥H{f s (τ, s[τ ]) -f s (τ, 0 n )}∥ + ∥Hf s (τ, 0 n )∥ ≤ √ γ∥Hs[τ ]∥ + ∥Hf s (τ, 0 n )∥ ≤ √ γw[τ ] + ∥Hf s (τ, 0 n )∥ = w[τ + 1],
where the second inequality comes from Lemma 4. This justifies (B.3). Meanwhile, the solution w of (B.4) is

w[t] = √ γ t-t 0 w[t 0 ] + t-1 τ =t 0 √ γ t-τ -1 ∥Hf s (τ, 0 n )∥,
which yields (B.1), and, Assumption 4 and (B.1) implies lim sup

t→∞ ∥s[t]∥ ≤ ∥H -1 ∥ lim sup t→∞ ∥Hs[t]∥ ≤ ∥H -1 ∥∥H∥ sup τ ≥t 0 ∥q ⊤ ⊗n F (t, p ⊗n 0 n )∥ 1 - √ γ ≤ ∥H -1 ∥∥H∥ ∥q∥ √ N M (0) 1 - √ γ
which completes the proof. ■ Now, we analyze the behavior of ( 18) with ( 14), which describes the evolution of the overall system at every integer time t. For this, with ξ 1 and ξ from (17), we introduce a Lyapunov function

V (ξ 1 , ξ, s) = ∥H(ξ 1 -s)∥ + η∥ ξ∥
where η > L∥q∥∥R∥∥H∥/ √ γ. Then, by recalling (18),

V [t + 1] = ∥H(ξ 1 [t + 1] -s[t + 1])∥ + η∥ ξ[t + 1]∥ = H{q ⊤ ⊗n F (t, p ⊗n ξ 1 [t]) -q ⊤ ⊗n F (t, p ⊗n s[t]) + q ⊤ ⊗n F (t, p ⊗n ξ 1 [t] + R ⊗n ξ[t]) -q ⊤ ⊗n F (t, p ⊗n ξ 1 [t])} + η (Λ K-1 Z ⊤ ) ⊗n {F (t, p ⊗n s[t]) + F (t, p ⊗n ξ 1 [t]) -F (t, p ⊗n s[t]) + F (t, p ⊗n ξ 1 [t] + R ⊗n ξ[t]) -F (t, p ⊗n ξ 1 [t])} ≤ H{q ⊤ ⊗n F (t, p ⊗n ξ 1 [t]) -q ⊤ ⊗n F (t, p ⊗n s[t])} + Hq ⊤ ⊗n {F (t, p ⊗n ξ 1 [t] + R ⊗n ξ[t]) -F (t, p ⊗n ξ 1 [t])} + |λ 2 | K-1 η ∥Z ⊤ ⊗n {F (t, p ⊗n ξ 1 [t]) -F (t, p ⊗n s[t])}∥ + ∥Z ⊤ ⊗n {F (t, p ⊗n ξ 1 [t] + R ⊗n ξ[t]) -F (t, p ⊗n ξ 1 [t])}∥ + ∥Z ⊤ ⊗n F (t, p ⊗n s[t])∥ .
The above inequality can be simplified by the following three inequalities. First, by Lemma 4,

∥H{q ⊤ ⊗n F (t, p ⊗n ξ 1 [t]) -q ⊤ ⊗n F (t, p ⊗n s[t])}∥ = ∥H{f s (t, ξ 1 [t]) -f s (t, s[t])}∥ ≤ √ γ∥H(ξ 1 [t] -s[t])∥ Second, by Assumption 4, ∥F (t, p ⊗n ξ 1 [t] + R ⊗n ξ[t]) -F (t, p ⊗n ξ 1 [t])∥ ≤ N i=1 ∥f i (t, p i ξ 1 [t] + (R i ⊗ I n ) ξ[t]) -f i (t, p i ξ 1 [t])∥ 2 1/2 ≤ N i=1 L 2 ∥(R i ⊗ I n ) ξ[t]∥ 2 1/2 = L∥R ⊗n ξ[t]∥ ≤ L∥R∥∥ ξ[t]∥,
where R i is the i-th row of R, and similarly the third one is

∥F (t, p ⊗n ξ 1 [t]) -F (t, p ⊗n s[t])∥ ≤ L∥p∥∥ξ 1 [t] -s[t]∥ = L∥p∥∥H -1 H(ξ 1 [t] -s[t])∥ ≤ L∥p∥∥H -1 ∥∥H(ξ 1 [t] -s[t])∥.
Then, we have that

V [t + 1] ≤ √ γ∥H(ξ 1 [t] -s[t])∥ + ∥H∥∥q∥L∥R∥∥ ξ[t]∥ + |λ 2 | K-1 η∥Z∥L∥p∥∥H -1 ∥∥H(ξ 1 [t] -s[t])∥ + |λ 2 | K-1 η∥Z∥L∥R∥∥ ξ[t]∥ + |λ 2 | K-1 η∥Z∥∥F (t, p ⊗n s[t])∥ ≤ √ γV [t] + |λ 2 | K-1 ηL∥Z∥M 1 V [t] + |λ 2 | K-1 η∥Z∥∥F (t, p ⊗n s[t])∥,
where M 1 := max{∥p∥∥H -1 ∥, ∥R∥/η}.

For the given ϵ, let K min be a positive integer such that

|λ 2 | K min ηLM 1 ∥Z∥ ≤ 1 - √ γ 2 (B.5) |λ 2 | K min 2ηM 1 M (∥p∥M s ) √ N ∥Z∥ 1 - √ γ ≤ ϵ 2 . (B.6)
Then, for all K > K min ,

V [t + 1] -V [t] ≤ - (1 - √ γ) 2 V [t] + |λ 2 | K-1 η∥Z∥∥F (t,
p i (ξ 1 [t] -s[t]) + (R i ⊗ I n ) ξ[t] ≤ max ∥p∥∥H -1 ∥, ∥R∥ η lim sup t→∞ V [t] ≤ |λ 2 | K-1 2η∥Z∥ 1 - √ γ √ N M (∥p∥M s )M 1 ≤ ϵ (B.9)
where we used x = p ⊗n ξ 1 + R ⊗n ξ. This completes the proof for (10) of Theorem 1. Now, in order to inspect the behavior of the system over the fractional time, let us apply the transformation of (17) to (12), which yields, for k = 1,

• • • , K -1, ξ 1 [t k ] = q ⊤ ⊗n W k-1 ⊗n F (t 0 , x[t 0 ]) = q ⊤
⊗n F (t 0 , x[t 0 ]) = q ⊤ ⊗n x[(t + 1) 0 ] = ξ 1 [(t + 1) 0 ] (B.10) in which, the third equality can also be seen from ( 13). Similarly, for k = 1,

• • • , K -1, ξ[t k ] = Z ⊤ ⊗n W k-1 ⊗n F (t 0 , x[t 0 ]) = (Λ k-1 Z ⊤ ) ⊗n F (t 0 , x[t 0 ]) = Λ k-K ⊗n (Λ K-1 Z ⊤ ) ⊗n F (t 0 , x[t 0 ]) = Λ k-K ⊗n Z ⊤ ⊗n W K-1 ⊗n F (t 0 , x[t 0 ]) = (Λ -1 ) K-k
N | K-k ≤ ∥p∥∥H -1 ∥ M 1 ϵ 2 + ∥R∥/η M 1 ϵ 2|λ N | K-k ≤ ϵ 2 1 + 1 |λ N | K-k ,
where the last two inequalities come from (B.12) and the definition of M 1 , and this completes the proof.

C Proof of Corollary 1

In this proof, we show that, after K -1 times execution of (5b), the solution of the overall system from the initial condition enters a positively invariant set, in which ( 19) holds. For this, let us first construct a few sets as C s 1 = {q ⊤ ⊗n F (0, x) :

x i ∈ C, i ∈ N } ⊂ R n , C s t+1 = {q ⊤ ⊗n F (t, p ⊗n s) : s ∈ C s t } ⊂ R n , ∀t ≥ 1, C s ∞ = ∞ t=1 C s t ⊂ R n ,
in which, C s t is the set of all possible s[t] for each t ≥ 1, which is bounded. The set C s ∞ is also bounded by Lemma B.1. Now, for the overall state x, consider two more sets: On the other hand, the set C is reached within K -1 executions of (5b) from the initial time 0 0 . Indeed, s[1] = q ⊤ ⊗n F (0, x[0]) ∈ C s 1 ⊂ C s ∞ , ∥ ξ[1]∥ = ∥(Λ K-1 Z ⊤ ) ⊗n F (0, x[0])∥ ≤ |λ 2 | K-1 ∥Z∥∥F (0, x[0])∥ ≤ δ, and ∥H(ξ 1 [1] -s[1])∥ = ∥H{q ⊤ ⊗n W K-1 ⊗n F (0, x[0])q ⊤ ⊗n F (0, x[0])}∥ = 0. Therefore, (x[t], s[t]) remains in C for all t ≥ 1.

C ′ = {x : ∥H(ξ 1 -s)∥ ≤ ϵ 0 , s ∈ C s ∞ , ∥ ξ∥ ≤ δ} ∪ C N ⊂ R nN , C = {(x, s) ∈ C ′ × C s ∞ : ∥H(ξ 1 -s)∥ ≤ ϵ 0 , ∥ ξ∥ ≤ δ} ⊂ R nN × R n ,
Finally, ( 19) is proved because, in the set C,

∥x i [t] -p i s[t]∥ = ∥p i (ξ 1 [t] -s[t]) + (R i ⊗ I n ) ξ[t]∥ ≤ max ∥p∥∥H -1 ∥, ∥R∥ (∥H(ξ 1 [t] -s[t])∥ + ∥ ξ[t]∥) ≤ ϵ.
To show (20), we note that (B.10) and (B.11) still hold. Therefore, for all k = 1, • • • , K -1, i ∈ N , and t ≥ 1,

∥x i [t k ] -p i s[t + 1]∥ = ∥p i (ξ 1 [t k ] -s[t + 1]) + (R i ⊗ I n ) ξ[t k ]∥ ≤ max ∥p∥∥H -1 ∥, ∥R∥ × (∥H(ξ 1 [t k ] -s[t + 1])∥ + ∥ ξ[t k ]∥) ≤ max ∥p∥∥H -1 ∥, ∥R∥ × (∥H(ξ 1 [t + 1] -s[t + 1])∥ + ∥ ξ[t + 1]∥/|λ N | K-k ) ≤ max ∥p∥∥H -1 ∥, ∥R∥ ϵ 0 + δ |λ N | K-k ≤ ϵ 2 1 + 1 |λ N | K-k ,
which completes the proof of (20).

Fig. 1 .

 1 Fig. 1. The network is initially given by the graph (a) and fixed for t ∈ [0, 100). Following the arrival of five new agents, the network is then represented by the graph (b) for t ∈ [100, 200). After two green triangle agents depart from the network, the network changes to the graph (c).

Fig. 2 .

 2 Fig. 2. The values N [t] ± 0.5 over time are indicated by the gray dashed lines. Colored solid curves represent the network size estimation results obtained by the proposed algorithm [(a): all the state trajectories of agents at each integer count are depicted] and by IAM and OSM [(b): the state trajectory of selected agent is depicted].

Fig. 2

 2 Fig. 2 illustrates the comparison of the network size estimation results obtained by the proposed algorithm and other algorithms including the Initial condition Averaging Method (IAM) by (Shames, Charalambous, Hadjicostis, & Johansson, 2012) and Order Statistics Method (OSM) by[START_REF] Varagnolo | Distributed cardinality estimation in anonymous networks[END_REF]. Specifically, Fig.2(a) depicts the estimation result obtained by the proposed algorithm, where all the state trajectories, at each integer count, i.e., x i [t 0 ], are shown. Here, K is set to be 25 for multi-step coupling and µ to be 0.1 for the Metropolis-Hastings coupling weight. Each state converges to the total number of connected agents with a deviation no greater than 0.5, regardless of agents joining or leaving from the network. This enables every agent to obtain the network size accurately. In contrast,

Fig. 4 .

 4 Fig. 4. Each PageRank score pi[t] over time is indicated by the gray dashed lines. Colored solid curves represent the PageRank score estimation results obtained by the proposed algorithm [(a): all the state trajectories of agents at each integer count are depicted] and by RA and SANI [(b): the state trajectory of selected agent is depicted].

  where C N is N -ary Cartesian power of C, i.e.,C N = C × C × • • • × C K min such that |λ 2 | K min ∥Z∥ sup x∈C ′ ∥F (t, x)∥ ≤ δ, ∀t ≥ 0.We claim that the set C is positively invariant for any K > K min . To see this, suppose that, for any τ ≥ 1,s[τ ] ∈ C s τ ⊂ C s ∞ , ∥ ξ[τ ]∥ ≤ δ, and ∥H(ξ 1 [τ ] -s[τ ])∥ ≤ ϵ 0 , so that (x[τ ], s[τ ]) ∈ C. Then, it follows that s[τ + 1] ∈ C s τ +1 ⊂ C s ∞ , ∥ ξ[τ + 1]∥ ≤ |λ 2 | K-1 ∥Z∥∥F (τ, x[τ ])∥∥ ≤ δ,and∥H(ξ 1 [τ + 1] -s[τ + 1])∥ ≤ ∥H{ξ 1 [τ + 1] -q ⊤ ⊗n F (τ, p ⊗n ξ 1 [τ ])}∥ + ∥H{q ⊤ ⊗n F (τ, p ⊗n ξ 1 [τ ]) -s[τ + 1])}∥ ≤ L∥q∥∥R∥∥H∥δ + ∥H{f s (τ, ξ 1 [τ ]) -f s (τ, s[τ ])}∥ ≤ L∥q∥∥R∥∥H∥δ + √ γ∥H(ξ 1 [τ ] -s[τ ])∥ ≤ ϵ 0in which, we used∥ξ 1 [τ + 1] -q ⊤ ⊗n F (τ, p ⊗n ξ 1 [τ ])∥ = ∥q ⊤ ⊗n {F (τ, p ⊗n ξ 1 [τ ] + R ⊗n ξ[τ ]) -F (τ, p ⊗n ξ 1 [τ ])}∥ ≤ L∥q∥∥R∥∥ ξ[τ ]∥ ≤ L∥q∥∥R∥δ.

degree of node i is denoted by d out i

	and in-degree of node i is de-
	noted by d i := |N i |. On the contrary, out-neighbors of node i are represented by N out i := {j ∈ N |(i, j) ∈ E}
	and out-

  . A distributed observer for a discrete-time linear system. In Proceedings of the IEEE 58th Conference on Decision and Control (CDC), pp. 367-372.Before proving Lemma 4, we first claim that, for each t ∈ Z and s 1 , s 2 ∈ R n , there exists s ∈ R n in the line connecting s 1 and s 2 such that {f s (t, s 2 ) -f s (t, s 1 )} ⊤ H 2 {f s (t, s 2 ) -f s (t, s 1 )}

	Yun, H., Shim, H., & Ahn, H. S. (2018). Initialization-
	free privacy-guaranteed distributed algorithm for eco-
	nomic dispatch problem. Automatica, 102, 86-93.
	Zaki, N. N., Berengueres, J. J., & Efimov, D. D. (2012).
	Detection of protein complexes using a protein rank-
	ing algorithm. Proteins: Structure, Function, and
	Bioinformatics, 80 (10), 2459-2468.
	A Proof of Lemma 4
	≤ (s 2 -s 1 ) ⊤ ∂f s ∂s	(t, s) ⊤ H 2 ∂f s ∂s	(t, s)(s 2 -s 1 ).	(A.1)
	It can be proved by the mean-value theorem with a trick.
	Let the left-hand side of the equality (A.1) as ∆(t, s 1 , s 2 )
	for convenience. With a variable c ∈ R, define
	g t,s1,s2

  p ⊗n s[t])∥. (B.7)

	By Assumption 4 and by Lemma B.1,
		√
	lim sup	∥F (t, p ⊗n s[t])∥ ≤	N M (∥p∥ lim sup	∥s[t]∥)
	t→∞	√			t→∞
		≤	N M (∥p∥M s )
	Using this and (B.7), the ultimate bound of V is obtained
	as
	lim sup t→∞	V [t] ≤ |λ 2 | K-1 2η∥Z∥ 1 -√ γ	√	N M (∥p∥M s ). (B.8)
	Therefore, for each agent i ∈ N and K > K min ,
	lim sup

t→∞ ∥x i [t] -p i s[t]∥ = lim sup t→∞

  ⊗n ξ[(t + 1) 0 ].

	Hence,
	∥ ξ[t k ]∥ ≤	∥ ξ[(t + 1) 0 ]∥ |λ N | K-k	(B.11)
	for each k = 1, • • • , K -1.
	On the other hand, by (B.8) and (B.6), we have
	lim sup
	≤ lim sup t→∞	V [t] ≤	2M 1 ϵ	.	(B.12)
	Therefore, for each k = 1, • • • , K -1,
	lim sup
			+	∥R∥ η	η∥ ξ[t k ]∥
	≤ lim sup

t→∞ max{∥H(ξ 1 [t] -s[t])∥, η∥ ξ[t]∥} t→∞ ∥x i [t k ] -p i s[t + 1]∥ = lim sup t→∞ ∥p i (ξ 1 [t k ] -s[t + 1]) + (R i ⊗ I n ) ξ[t k ]∥ ≤ lim sup t→∞ ∥p∥∥H -1 ∥∥H(ξ 1 [t k ] -s[t + 1])∥ t→∞ ∥p∥∥H -1 ∥∥H(ξ 1 [t + 1] -s[t + 1])∥ + ∥R∥ η η∥ ξ[t + 1]∥ |λ

Codes are at https://github.com/donggil-lee/multi-agent.