On Asymptotic Consensus Value in Directed Random Networks

Distributed consensus algorithms have attracted a significant amount of attention in the past few years. Besides their wide range of applications in distributed and parallel computation [1], distributed control [2], [3] and robotics [4], they have also been used as models of opinion dynamics and belief formation in social networks [5], [6]. The central focus in this vast body of literature is to study whether a group of agents in a network, with local communication capabilities can reach a global agreement, using simple, deterministic information exchange protocols. 1 In recent years, there has also been some interest in understanding the behavior of consensus algorithms in random settings [7]- [12]. The randomness can be either due to the choice of a randomized network communication protocol or simply caused by the potential unpredictability of the environment in which the distributed consensus algorithm is implemented [13]. It is recently shown that consensus algorithms over i.i.d. random networks lead to a global agreement on a possibly random value, as long as the network is connected in expectation [10]. While different aspects of consensus algorithms over random switching networks, such as conditions for convergence [7]- [10] and the speed of convergence [13], have been widely studied, a characterization of the distribution of the asymptotic consensus value has attracted little attention. Two notable exceptions are Boyd et al. [14], who study the asymptotic behavior of the random consensus value in the special case of symmetric networks, and Tahbaz-Salehi and Jadbabaie [15], who compute the mean and variance of the consensus value for general i.i.d. graph processes. Nevertheless, a complete characterization of the distribution of the asymptotic value for general asymmetric random consensus algorithms remains an open problem. In this paper, we study asymptotic properties of consensus algorithms over a general class of switching, directed random graphs. More specifically, building on the results of [15], we derive closed-form expressions for the mean and an upper-bound for the variance of the asymptotic consensus value, when the underlying network evolves according to an i.i.d. directed random graph process. In our model, at each time period, a directed communication link is established between two agents with some exogenously specified probability. Due to the potential asymmetry in pairwise communications between different agents, the asymptotic value of consensus is not guaranteed to be the average of the initial conditions. Instead, agents will asymptotically agree on some random value in the convex hull of the initial conditions. Furthermore, our closed-form characterization of the variance provides a quantitative measure of how dispersed the random agreement point is around the average of the initial conditions in terms of the fundamentals of the model, namely, the structure of the network, the exogenous probabilities of communication, and the initial conditions. The rest of the paper is organized as follows. In the next section, we describe our model of random consensus algorithms. In Sections III and IV, we derive an explicit expression for the mean and an upper bound for the variance of the limiting consensus value over switching directed random graphs, respectively. Section V contains simulations of our results and Section VI concludes the paper. Consider the discrete-time linear dynamical system where k ∈ {1, 2, . . . } is the discrete time index, x(k) ∈ R n is the state vector at time k, and {W k } ∞ k=1 is a sequence of stochastic matrices. We interpret (1) as a distributed scheme where a collection of agents, labeled 1 through n, update their state values as a convex combination of the state values of their neighbors at the previous time step. Given this interpretation, x i (k) corresponds to the state value of agent i at time k, and W k captures the neighborhood relation between different agents at time k: the i j element of W k is positive only if agent i has access to the state of agent j. For the remainder of the paper, we assume that the weight matrices W k are randomly generated by an independent and identically distributed matrix process. We say that the dynamical system (1) reaches consensus asymptotically on some path {W k } ∞ k=1 , if along that path, there exists x * ∈ R such that x i (k) → x * for all i as k → ∞. We refer to x * as the asymptotic consensus value. It is wellknown that for i.i.d. random networks, the dynamical system (1) reaches consensus on almost all paths if and only if the graph corresponding to the communications between agents is connected in expectation. More precisely, it was shown in [10] that W k . . .W 2 W 1 -→ 1d T almost surelywhere d is some random vectorif and only if the second largest eigenvalue modulus of EW k is subunit. Clearly, under such conditions, the dynamical system in (1) reaches consensus with probability one where the consensus value is a random variable equal to x * = d T x(0), where x(0) is the vector of initial conditions. A complete characterization of the random consensus value x * is an open problem. However, it is possible to compute its mean and variance in terms of the first two moments of the i.i.d. weight matrix process. In [15], the authors prove that the conditional mean of the random consensus value is given by the random consensus value are given by and its conditional variance is equal to where v 1 (•) denotes the normalized left eigenvector corresponding to the unit eigenvalue, and ⊗ denotes the Kronecker product. In [16], the authors derive an explicit expression for the mean and variance of x * in the particular case of a switching Erdös-Rényi random graph process. In the following, we shall use (3) to extend these results and derive an explicit expression for the mean and an upper-bound for the variance of the asymptotic consensus value over a wider class of switching, directed random graph processes. Consider a connected undirected graph G c = (V, E c ) with a fixed set of vertices V = [n] , and unweighted edges (no self-loops allowed). Each undirected edge in E c represents a potential communication channel between nodes i and j, where this channel can be used to send information in both directions. In this paper, we focus on directed communications, i.e., the event of node i sending information towards node j is independent from the even of node j sending information towards i. In this context, it is convenient to interpret an undirected edge {i, j} ∈ E c as the union of two independent directed edges, {(i, j) , ( j, i)}, where the ordered pair (i, j) represents a directed link from node i to node j. In this paper, we study randomized time-switching consensus processes. In particular, in each discrete time slot k ≥ 1, we construct a random directed graph G k = (V, E k ), with E k ⊆ E c , such that the existence of a directed edge (u, v) ∈ E k is determined randomly and independently of all other directed edges (including the reciprocal edge (v, u)) with a probability p uv ∈ (0, 1) for (u, v) ∈ E c , and p uv = 0 for (u, v) ∈ E p . In other words, in each time slot, we randomly select a subset E k of directed links chosen from a set of candidate (directed) links in E c . We are specially interested in the case in which the probability p uv of existence of a directed link (u, v) depends exclusively on the node that receives information via that link, i.e., Pr where p v ∈ (0, 1). In this setting, we can model the ability of a node to 'listen' to their neighboring nodes. For example, in the context of opinion dynamics in social networks [5], [6], [12], the probability p v can represent the tendency of the individual at node v to take into account the opinion of her neighbors (which could depend, for example, on how many acquaintances the individual has). Let us denote by A c the symmetric adjacency matrix of the graph G c , where entries a i j = 1 if {i, j} ∈ E c , and a i j = 0 otherwise. We define the degree of node i as d i = ∑ n j=1 a ji , and the associated degree matrix as D c = diag(d i ). We also denote the random (nonsymmetric) adjacency matrix associated with We denote the in-degrees of uv , and the indegree matrix as Dk = diag( We describe the consensus dynamics in (1) via a sequence of stochastic matrices {W k } ∞ k=1 associated to the sequence of random directed graphs k=1 as follows: Notice that adding the identity matrix to the adjacency in ( 5) is equivalent to adding a self-loop to every vertex in V . We include these self-loops to avoid singularities associated with the existence of isolated nodes in G k (for which di = 0, and Dk would not invertible). Since G c is connected and p v > 0 for all v ∈ V , the expected communications graph is connected and the stochastic dynamical system in (1) reaches consensus on almost all paths, although the asymptotic consensus value x * is a random variable (not the initial average). In the following subsections, we first derive closed-form expressions for the mean of x * , and an upper-bound for the variance var (x * ) in terms of the following elements: 1) The set of initial conditions, {x u (0)} u∈V , 2) the set of nodes properties, {(p u , d u )} u∈V , and 3) the network topology, via the eigenvalues of the expected matrix EW k . As we shall show in Section IV, our expression for the variance has a nice interpretation, since it separates the influence of each one of the above elements into three multiplicative terms. In the next subsection we provide the details regarding our analysis of the expectation of x * . We use (2) to study the mean of the consensus value. We first derive an expression for EW k , and then study its dominant left eigenvector v 1 (EW k ). For notational convenience, we define the random variable z i 1/ di + 1 where di ∼ Ber (d i , p i ), and denote its first and second moments as The diagonal entries of EW k are then given by Ew ii = E 1/ di + 1 , which present the following explicit expression (see Appendix I for details): where q i = 1p i . Furthermore, the off-diagonal entries of EW k are equal to (see Appendix I for details): Taking ( 6) and ( 7) into account, we can write EW k as follows: where Σ diag M (1) i . As expected, it is easy to check that EW k is a stochastic matrix, i.e., (EW k ) (8) we can write E (x * ) explicitly in terms of d i and p i , as follows: Theorem 1: Consider the random adjacency matrix Ãk in (4) and the associated (random) stochastic matrix W k in (5). The expectation of the asymptotic consensus value of ( 1) is given by where . Proof: Our proof is based on computing v 1 (EW k ) and applying (2). Let us define v v 1 (EW k ) and w (I -Σ)v. From (8), we have that the eigenvalue equation corresponding to the dominant left eigenvector of This last equation can be written as w T D -1 c A T c = w T . The solution to this equation is the stationary distribution of the Markov chain with transition matrix where we include the normalizing parameter Hence, from (2) we have (1) i in ( 6), we reach (9) via simple algebraic simplifications. In general, the asymptotic mean E (x * ) does not coincide with the initial average x0 = 1 n ∑ i x i (0). There is a simple technique, based on Theorem 1, that allows us to make the expected consensus value to be equal to the initial average. This technique consists of using y i (0) = (ρw i ) -1 x i (0) as initial conditions in (1). Hence, one can easily check that the asymptotic consensus value E (y * ) equals the initial average x0 . In this section, we derive an expression that explicitly relates the variance var (x * ) with the three elements that influences it, namely, the set of initial conditions {x u (0)} u∈V , the nodes properties {(p u , d u )} u∈V , and the network structure (via the eigenvalues of the expected matrix EW k ). For simplicity in notation, we denote E (W k ⊗ W k ) and (EW k ⊗ EW k ) by R and Q, respectively. Our analysis starts in expression (3), which can be rewritten as Hence, we can upper-bound the variance of the asymptotic consensus value as follows: From the rules of Kronecker multiplication, we can write the first factor in terms of the initial conditions as In the following, we derive an upper bound for the second factor v 1 (R)v 1 (Q) ∞ in terms of the nodes properties and the network structure. Our approach to bound v 1 (R)v 1 (Q) ∞ is based on the observation that both R and Q are n 2 × n 2 stochastic matrices, and the dominant left eigenvectors v 1 (R) and v 1 (Q) are stationary distributions of the Markov chains with transition matrices R and Q. We denote these distributions by v 1 (R) π and v 1 (Q) π, respectively. In this setting, we can apply the following lemma from [17] which studies the sensitivity of the stationary distribution of Markov chains: Lemma 2: Consider two Markov chains with transition matrices Q and R, and stationary distributions π and π, respectively. We define G = I -Q, and denote its pseudoinverse by where κ s = max i, j g † i j is called the condition number of the chain described by the transition matrix Q. In the next subsections, apply the above lemma to bound the factor v 1 (R)v 1 (Q) ∞ . In the first subsection we compute an explicit expression for the norm of the perturbation R -Q ∞ as a function of the properties of the nodes. In the second subsection, we study the coefficient κ s in terms of the eigenvalues of EW k . Our approach is based on studying the entries of the and compare them with the entries of Q = EW k ⊗ EW k . The entries of Q and R are of the form E(w i j )E(w rs ) and E(w i j w rs ), respectively, with i, j, r, and s ranging from 1 to n. These entries can be classified into seven different cases depending on the relations between the indices. In the Appendices II and III, we present explicit expressions for each one of these cases. A key observation in our approach is to notice the pattern that the above entries induces in the matrices R and Q. For sake of clarity, we illustrate this pattern for n = 3 in Fig. 1, where the numbers in parenthesis correspond to each one of the seven cases identified in the Appendices. Since the entries of R follow the same pattern as the entries of Q (although these entries are different), the perturbation matrix ∆ = R -Q also follows the same pattern as R and Q. Hence, comparing the entries of Q and R we can easily deduce the following seven cases for the entries of ∆: , where M (2) i E 1/( di + 1) 2 can be written as a hypergeometric function that depends on p i and d i (see Appendix III). From the above entries, we can compute via simple, but tedious, algebraic manipulations the following expression for the infinity norm of the perturbation: where Note that S i is a function that depends exclusively on the sequence of nodes properties {(p u , d u )} u∈V ; hence, R -Q ∞ depends exclusively on the set of nodes degrees and probabilities, but it is independent on how these nodes interconnect. In the following subsection, we show how the pattern of interconnection among nodes influences the upper-bound of the variance in (10) via the condition number κ s in (12). In this subsection we derive an explicit relationship between the condition number κ s and the network structure via the spectral properties of EW k = Σ + (I -Σ)D -1 c A T c . This result will then be used to bound the variance of x * in (10). We base our analysis on a the following bound, derived by Meyer in [18], relating κ s with the eigenvalues of Q, denoted as Before we present our main result, we need some notation. Denote by {λ i } n i=1 and µ j n 2 j=1 the set of eigenvalues of EW k and EW k ⊗ EW k , respectively. The ordering of the eigenvalue sequences is determined by their distance to 1, i.e., |1 -λ i | ≤ 1 -λ j for i ≤ j. Hence, our result can be stated as follows: Theorem 3: The variance of the asymptotic consensus value of (1) can be upper-bounded by where {λ i } n i=1 are the eigenvalues of EW k , and the product Proof: We start our proof from ( 10) where we have used Lemma 2 in inequality (a), the expression for R -Q ∞ in ( 13) in equality (b), and the upper bound in (14) in step (c). We obtain the statement of the theorem by applying the following standard property of the Kronecker product: Remark 4: The bound in ( 15) separates the variance into three multiplicative terms representing each one of the following elements (for convenience, we have underlined each one of these terms in ( 15)): (A) This first term exclusively depends on the initial condition as indicated by ( 11). (B) The second term depends solely on the nodes properties d i and p i . (C) The last term represents the influence from the overall graph structure via the spectral properties of EW k . Remark 5: It is specially interesting to study the implications of Term (C) in the asymptotic variance. For example, given the sequences of degrees and probabilities, {d i } n i=1 and {p i } n i=1 , the influence of the network structure on the variance is given via Term (C). Since the eigenvalues of EW k are key in this term, it is interesting to briefly describe the homogeneous Markov chain with transition matrix P EW k = Σ+ (I -Σ) D -1 c A T c . This Markov chain presents a selfloop in each one of its n states with transition probability i , as well as a link from i to j with transition probability . From an analysis proposed by Meyer in [18], we conclude that Term (C) is primarily governed by how close the subdominant eigenvalues of EW k are from 1. In particular, the further the subdominant eigenvalues of EW k are from 1, the lower the upper bound in (15). In the next section, we illustrate the relationship between spectral properties of EW k and the variance of the asymptotic consensus value with several numerical simulations. In this section we present several numerical simulations illustrating our results. In the first subsection, we numerically verify the result in Theorem 1. In the second subsection we present some examples to illustrate the implications of Theorem 3 in the variance of the asymptotic consensus value. and3, with a set of followers. Connections between leaders are represented as pairs of directed links with existence probabilities p uv and p vu , with u,v ∈ {1,2,3}. We have also included in this figure the connection of each leader and one of its followers (drawn as a smaller circle on top of each leader). In this case, the directed link incoming the follower (leader) exists with probability p l f (p f l ). The rest of followers are represented under each leader, where we represent each pair of directed links as a single undirected edge for clarity. In our first simulation we take a graph G c composed by 3 stars connected in a chain (see Fig. 2). This graph is intended to represent, in a social context, three leaders with a set of followers. In Fig. 2, we represent the leaders using large circles marked as 1,2, and 3, and the followers using smaller circles. We assume that each follower only listen to one leader (the center of a star) and nobody else. In this particular example, the first, second and third leaders have 4, 8, and 16 followers, respectively. In each time step, a directed random graph G k is built by choosing a set of directed communication links from G c . We have fixed the probability of existence of a directed link coming into followers to be equal to p l f = 1/2, for all such links (see Fig. 2). Also, the probability of existence of a directed link coming into a leader is inversely proportional to the degree of the leader. More specifically, we have chosen p uv = 1/d v , for v = 1, 2, and 3, where p uv represents the existence of a directed communication link from the u-th node to the v-th leader. In this example, we compute the asymptotic consensus value for 100 realizations of a random consensus algorithm with initial conditions x i (0) = i/N. We represent the histogram of these realizations in Fig. 3, where the empirical average of the 100 realizations equals 0.5077. The theoretical expectation for the asymptotic consensus value applying Theorem 1 equals E (x * ) = ∑ i ρw i (i/N) =0.4595, which is in great accordance with the empirical value. For future reference, we have also computed the empirical standard deviation of the 100 realizations to be 0.0298, and the three eigenvalues of EW k closest to 1 are {0.9823, 0.9449, 0.75}. In this subsection, we numerically illustrate some of the implications of Theorem 3 in the variance of the asymptotic consensus value. Although the upper bound stated in this Theorem is not tight, there are some qualitative implications that are consistent throughout our numerical experiments. For example, as we mentioned in Remark 5, given a set of initial conditions and nodes properties, the upper bound in ( 15) is primarily governed by how close the subdominant eigenvalues of EW k are from 1. In particular, the further the subdominant eigenvalues of EW k are from 1, the lower the upper bound. We illustrate in the following numerical experiments that the lower the upper bound, the lower one should expect the variance to be. In the first experiment, we slightly modify the network described in the previous subsection and study the influence of this modification on the eigenvalues of EW k and the variance of x * . Our first modification is a change in the probabilities of existence of a directed link from the uth node to the v-th leader without changing the network topology. In particular, we choose the new probabilities to be p uv = 3/d v , for v = 1, 2, and 3. In the modified network, the eigenvalues of EW k closest to 1 are {0.9789, 0.9372, 0.75}. Hence, since the effect of our modification on the eigenvalues is to move them away from 1, we should expect the variance of x * to be reduced according to Remark 5. In fact, running 100 random consensus algorithms with the same initial conditions, x i (0) = i/N, using our new probabilities, we obtain a standard deviation 0.0286 (which is less than the 0.0298 obtained before). In the second experiment, we illustrate how the larger the gap between the eigenvalues of EW k and 1, the smaller the variance of x * . In this case, apart from keeping the new set of probabilities p uv = 3/d v , we convert the 3-chain of leaders into a 3-ring of leaders, as depicted in Fig. 4. In this case, the three eigenvalues of EW k closest to 1 are {0.9577, 0.9212, 0.75}, which are even further away from 1 than in the second example. Hence, as expected, the standard deviation of 100 random consensus algorithms is even smaller than in the second example, in particular 0.0274. In conclusion, our simulations are consistent with Theorem 3 and with the qualitative behavior described in Remark 5. We have studied the asymptotic properties of the consensus value in distributed consensus algorithms over switching, directed random graphs. While different aspects of consensus algorithms over random switching networks, such as conditions for convergence and the speed of convergence, have been widely studied, a characterization of the distribution of the asymptotic consensus for general asymmetric random consensus algorithms remains an open problem. In this paper, we have derived closed-form expressions for the expectation of the asymptotic consensus value as a function of the set of initial conditions, {x u (0)} u∈V , and the set of nodes properties, {(p u , d u )} u∈V , as stated in Theorem 1. We have also studied the variance of the asymptotic consensus value in terms of several elements that influence it, namely, (i) the initial conditions, (ii) nodes properties, and (iii) the network topology. In Theorem 3, we have derived an upper bound for the variance of the asymptotic consensus value that explicitly describes the influence of each one of these elements. We also provide an interpretation of the influce of the network topology on the variance in terms of the eigenvalues of the expected matrix EW k . From our analysis we conclude that, in most cases, the variance of x * is primarily governed by how close the subdominant eigenvalues of EW k are from 1. We have checked the validity of our predictions with several numerical simulations. We start by computing the entries of EW k . The diagonal entries of EW k are given by: Similarly, the non-diagonal entries of EW k result in: We now compute the possible entries in Q = EW k ⊗ EW k . The entries of the Kronecker matrix EW k ⊗ EW k , with EW k = Σ + (I -Σ) D -1 A T , present entries that can be classified into seven different cases depending on the relations between the indices. These are the cases, where we assume that all four indices i, j, r, and s are distinct: . We now turn to the computation of the elements of E[W k ⊗ W k ], which are of the form E(w i j w rs ). Again, we classify the entries into seven different cases depending on the relations between the subindices: Similarly, we also have: R 2 E (w ii w j j ) = M (1) i M (1) j , R 3 E (w ii w i j ) = a ji d i M (1) i -M (2) i , R 4 E (w i j w i j ) = a ji d i M (1) i -M (2) i , R 5 E (w ii w rs ) = a sr d r M (1) i (1) r , R 6 E (w i j w is ) = a ji a si d i (d i -1) 1 + 2M (2) i -3M (1) i , for d i > 1, 0, for d i = 1, R 7 E(w i j w ji ) = E(w i j w js ) = E(w i j w ri ) = E(w i j w r j ) = E(w i j w rs ) = a ji a sr d i d r 1 -M (1) i (1) r .