{'fraction_non_alphanumeric': 0.058021266759130836, 'fraction_numerical': 0.0146208969024503, 'mean_word_length': 3.8367651169290973, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 46, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Distributed Algorithms for Consensus and Coordination in the Presence of Packet-Dropping Communication Links - Part I: Statistical Moments Analysis Approach

This two-part paper discusses robustification methodologies for linear-iterative distributed algorithms for consensus and coordination problems in multicomponent systems, in which unreliable communication links may drop packets. We consider a setup where communication links between components can be asymmetric (i.e., component j might be able to send information to component i, but not necessarily vice-versa), so that the information exchange between components in the system is in general described by a directed graph that is assumed to be strongly connected. In the absence of communication link failures, each component i maintains two auxiliary variables and updates each of their values to be a linear combination of their corresponding previous values and the corresponding previous values of neighboring components (i.e., components that send information to node i). By appropriately initializing these two (decoupled) iterations, the system components can asymptotically calculate variables of interest in a distributed fashion; in particular, the average of the initial conditions can be calculated as a function that involves the ratio of these two auxiliary variables. The focus of this paper to robustify this double-iteration algorithm against communication link failures. We achieve this by modifying the double-iteration algorithm (by introducing some additional auxiliary variables) and prove that the modified double-iteration converges almost surely to average consensus. In the first part of the paper, we study the first and second moments of the two iterations, and use them to establish convergence, and illustrate the performance of the algorithm with several numerical examples. In the second part, in order to establish the convergence of the algorithm, we use coefficients of ergodicity commonly used in analyzing inhomogeneous Markov chains.


I. INTRODUCTION
The design of protocols and algorithms for distributed computation and control/decision tasks has attracted significant attention by the computer science, communication, and control communities (e.g., [1], [2], [3], [4], [5], [6], [7] and references therein). For example, given i) a collection of robots moving in the plane, ii) a collection of sensors in a sensor network, or iii) a collection of distributed energy resources in an electrical grid, the components may be interested in, respectively, i) agreeing on a common direction to follow (this common direction could be provided by a leader robot), ii) measurement averaging (with each sensor providing a local measurement of a global quantity), or iii) collectively providing a predetermined total amount of active power subject to the constraints of each distributed resource. In the control literature, the first and second problems are respectively known as consensus and average consensus (see, e.g., [2]), whereas the third problem can be considered as a distributed resource coordination problem [8], [9].
In this two-part paper, we consider multicomponent systems in which each component can exchange information with other components in its neighborhood in order to compute, in a distributed fashion, some quantity of interest. In our setup, communication links between components (nodes) can be asymmetric (i.e., component j might be able to send information to component i, but not necessarily vice-versa), a situation that arises in a wireless setting if the transmission power available to different nodes are also different. In this setting, the information exchange between components in the system can be described by a directed graph which is assumed to be strongly connected. Through an iterative process, nodes in the network are required to compute (using only information made available by their neighbors) the quantity of interest. In particular, we study linear-iterative algorithms in which each node j maintains a value (or a set of values) that is updated to be a weighted linear combination of node j's own previous value and the previous values of its neighboring nodes (i.e., nodes that transmit information to node j). The main focus of the paper is to develop strategies to robustify the linear-iterative algorithms described above against communication links that may drop packets.
In the context of consensus and average-consensus problems, an extensive literature in the control community focuses on the linear-iterative algorithms described above (e.g., [7], [10], [2], [11], [12], [13], [14] and references therein). These works have revealed that if the network topology satisfies certain conditions, the weights for the linear iteration can be chosen so that all the nodes asymptotically converge to the same value (even if the network connections are time-varying). Additionally, if the interconnection topology is invariant and bidirectional (i.e., if node j can send information to node i, then node i can send information to node j), simple techniques can be used to choose the weights of the linear iteration so as to ensure that, after running the linear iteration, the nodes will asymptotically reach consensus to the average of their initial values [2], [11], [12]. Other works have looked at the consensus and average-consensus problems when the interconnection topology is described by a directed graph. In particular, the authors of [15] focus on continuous-time linear iterations and state necessary and sufficient conditions for a network of integrators to asymptotically reach agreement to a common value (but not necessarily the average of their initial values). Similarly, the authors of [12] consider discretetime iterations, and provide necessary and sufficient conditions on the weights that allow the nodes to asymptotically reach consensus to the average of their initial values. Additionally, the work in [16], [17] discusses how average-consensus can be reached asymptotically with linear-iterative algorithms in which the nodes use fixed weights in their linear updates and also develops linear-iterative algorithms where the nodes adapt their weights in a distributed fashion so that asymptotically average-consensus is reached. In the context of resource coordination, there is some recent work [18], [8], [9] that also focuses on linear-iterative algorithms, similar to those used to address consensus and average-consensus problems. Other recent work has addressed the related problem of achieving consensus and averageconsensus in a multicomponent system where some nodes can exhibit malicious behavior [19], [20]. These works assume fault-free communication links, but are related to what we do in this paper in the sense that they can be used to handle unreliable nodes (as opposed to links).
In our development, we adopt a very general model for the communication modality between nodes, which allows asymmetric information structures, in the sense that if node i can transmit information to another node j, it is not necessarily true that node j can transmit information to node i. We only require that each node, apart from seeing incoming transmissions sent to it by neighboring nodes, knows the number of nodes that it can transmit information to, which in graph-theoretic terms is referred to as the out-degree of that node. In fact, in the proposed algorithm, each node will broadcast the same quantity to all receiving nodes, which simplifies the communication scheme between sending and receiving nodes (as it is not necessary for each sending node to separately communicate with each receiving node).
When the communication network is perfectly reliable (no packet drops), the collective dynamics of the linear iterations can be described by a discrete-time transition system with no inputs in which the transition matrix is column stochastic and primitive. Then, each node will run two identical copies of the linear iteration each of which, however is initialized differently depending on the problem to be solved.In this paper we mostly focus on the average consensus problem. Under proper initialization, it can be shown that each node will asymptotically calculate the desired value as a function of the outcomes of the two iterations. The details of these double-iteration approach are provided in [16], [17] for the average consensus case and in [8], [9] for the resource coordination problem. For the average-consensus problem, the double-iteration algorithm is a particular case of the algorithm in [21] (which is a generalization of the algorithm proposed in [22]), where the matrices describing each linear iteration are allowed to vary as time evolves, whereas in our setup (for the ideal case when there are no communication link failures) the transition matrix is fixed over time.
The focus of this paper is to robustify the double-iteration algorithm (informally described above and formally described in Section II) so that it can tolerate failures in communication links and converge to the average value. Our communication link reliability model assumes that at each time step, a communication link is unavailable with some probability. In other words, a packet containing information from node i to node j is dropped with some probability. Next we informally describe our robustification approach. Consider two nodes i and j, and assume that j receives information from node i but not necessarily vice-versa. Let us refer to j as the receiving node (or receiver) and i as the sending node (or sender). An important requirement is for the graph describing the communication network to be strongly connected, which implies every node must be able to act both as a sender and as a receiver. Then, for each of the two iterations each node performs, node i (the sender) will keep track of the following quantities of interest: i) its own internal state (as captured by the state variables maintained in the double iteration scheme of [21], [18]; ii) the total mass broadcasted so far (to be described in detail soon); and iii) the total received mass from each node l that sends information to node i. Similarly, for both iterations, each node j (the receiver) updates the value of its internal state to be a linear combination of its own previous internal state value (weighted by the inverse of the number of nodes that have j as a neighbor) and the difference between the two most recently received mass values from each of its neighbors (also weighted by the inverse of the number of nodes that have j as a neighbor). At time instant k, the total broadcasted mass by node j is the sum up to (and including) time step k of the sequence of values of node j's internal value, weighted by the inverse of the number of nodes that receive values from node j). Additionally, node j updates the value of the received mass from each node l that sends information to node j as follows: the received mass from node l is the total broadcasted mass sent by node l up to time k if the communication link from node l to node j is available at time step k; otherwise, the received mass remains the same as the most recently received mass from node l. An implicit assumption here is that messages broadcasted by node l are tagged with the sender's identity so that the receiving node j can determine where different packages have originated from.
Recent work that has addressed the consensus and average-consensus problems in the presence of unreliable communication links [23], [24], [25] has employed a communication link availability model similar to ours. The work in [23] assumes that the graph describing the communication network is undirected and that when a communication link fails it affects communication in both directions. Additionally, nodes have some mechanism to detect link unavailability and compensate for it by rescaling their other weights (so that the resulting transition matrix remains column stochastic). Following this strategy, the authors show asymptotic convergence to the average of initial conditions and also calculate the rate at which the variance of the total deviation from the average converges to zero. The work in [24] does not require the graph describing the communication network to be undirected and proposes two compensation methods to account for communication link failures. In the first method, the so-called biased compensation method, the receiving node compensates for the unavailability of an incoming link by adding the weight associated to the unavailable link to its own weight (so that the resulting matrix remains row stochastic). In the second method, called the balanced compensation method, the receiving node compensates for link unavailability by rescaling all the incoming link weights so that the resulting matrix remains row stochastic. The key in both methods is the fact that at each time step, the resulting weight matrix is row stochastic; the authors show that the nodes converge almost surely to the same value, but this value is not necessarily the average of the initial conditions. The work in [25], which does not require the communication graph to be undirected, proposes a correction strategy that corrects the errors in the quantity (state) iteratively calculated by each node, so that the nodes obtain the average of their initial values. This correction strategy is based on each node maintaining some auxiliary variable that accounts for the amount by which node i changes its state due to the updates from its neighbors, i.e., the nodes that can send information to node i. For their strategy to work and ensure that the nodes converge almost surely to the average consensus, the authors rely on the nodes sending acknowledgment messages and retransmitting information an appropriate number of times.
In [21], the authors proposed a gossip-based algorithm for average-consensus ver a directed graph where the transition matrices describing the nodes' collective dynamics change at every iteration step (depending on which node awakes). This scheme requires the node that is awake to perform an internal state update and send its internal state (weighted by the corresponding out-going link weight) to its neighbors. This approach results in generates a sequence of column stochastic matrices (not necessarily primitive) with the property that all the diagonal entries remain positive. The authors prove that by running two such iterations in parallel, one of them initialized with the values on which the average operation is to be performed and the other with the all-ones vector, each node will asymptotically achieve average consensus by taking the ratio of the two values in maintains. A key premise in their proof is that column stochasticity of the transition matrix is maintained over time, which requires sending nodes to know the number of nodes that are listening. This suggests that i) either the communication links are perfectly reliable, or ii) there is some acknowledgment and retransmission mechanism that ensures messages are delivered to the listening nodes at every round of information exchange. In this paper, we remove such assumptions and robustify the double-iteration algorithm against unreliable communication links using a pure broadcast-message model without any requirement for an acknowledgment/retransmission mechanism. Thus, despite the reliance of our algorithm on the ratio of two linear iterations, it is different both in the communication model we assume-a broadcast model in our case-and also in the nature of the protocol itself-our focus is on ensuring convergence in the presence of communication link failures.
An additional assumption made in [21] is that the diagonal entries of the transition matrix (at every step) remain positive. In our model, we originally consider that nodes do not drop self-packets. However, to ease the analysis, we remove this assumption and consider the case where self-packet drops are also allowed at every time step, which i) allows us to handle intermittent faults in the node processing device, and ii) removes the assumption that all diagonal entries must be positive at every step. Finally, the analysis machinery in [21] is quite different from the one used in this paper. We employ moment analysis of the two iterations to establish that they are linearly related as the number of steps goes to infinity, while [21] relies on establishing weak ergodicity of the product of the transition matrices as the number of steps goes to infinity. Finally, as it will be shown in the second part of this paper, our algorithm can be re-casted into a similar framework as the one in [21] by augmenting the dimension of the vector describing the collective dynamics to account for the packets that get dropped once there is a communication failure. Note, However, that the resulting matrices will be column stochastic but will not necessarily have strictly positive entries on their diagonals. In the second part of the paper, we provide an analysis framework to establish the convergence of our algorithm and generalizes the ideas in [21] to the case when the matrices describing the system collective dynamics do not have strictly positive diagonals. In this regard, we will show that even in the case where self-packet drops are not allowed, the resulting transition matrices might still have zero diagonal entries.
The remainder of this paper is organized as follows. Section II provides background on graph theory, introduces the communication model, and briefly describes the non-robust version of the double-iteration algorithm we use in this work. Section III describes the proposed strategy to robustify the double-iteration algorithm against communication link failures and illustrates the use/performance of the algorithms via several examples. The convergence analysis of the robustified double-iteration algorithm is provided in Section V. Concluding remarks are presented in Section VI.

II. PRELIMINARIES
This section provides background of graph-theoretic notions that are used to describe the communication network and the distributed consensus/coordination setup, introduces the basic communication link availability model, and reviews a previously proposed two-iteration algorithm that can be used to solve the class of problems addressed in this paper when the communication network is perfectly reliable.

A. Network Communication Model
Let discrete time instants be indexed k = 0, 1, . . . ; then, the information exchange between nodes (components) at each time instant k can be described by a directed graph . . , n} is the vertex set (each vertex-or node-corresponds to a system component), and if node j can receive information from node i at instant k. It is assumed that E[k] ⊆ E, ∀k ≥ 0, where E is the set of edges that describe all possibly available communication links between nodes; furthermore, the graph (V, E) is assumed to be strongly connected. All nodes that can possibly transmit information to node j are called its in-neighbors, and are represented by the set N − j = {i ∈ V : (j, i) ∈ E}. Note that there are self-loops for all nodes in G (i.e., (j, j) ∈ E for all j ∈ V). The number of neighbors of j (including itself) is called the in-degree of j and is denoted by The nodes that have j as neighbor (including itself) are called its out-neighbors and are denoted by N + j = {l ∈ V : (l, j) ∈ E}; the out-degree of node j is D + j = |N + j |. The existence of a communication link from node i to node j can be described in probabilistic terms as follows. At instant k, let x ji [k], ∀i, j ∈ V be an indicator variable that takes value 1 with probability q and takes value zero with probability 1 − q, i.e., Then, for all k ≥ 0, the existence of a communication link between nodes i and j can be described be another indicator variable ℓ ji [k] defined as It follows that E[k] contains the elements of E for which ℓ ji [k] = x ji [k] = 1.

B. Double-Iteration Algorithm Formulation in Perfectly Reliable Communication Networks
When the communication network of a multi-component system is perfectly reliable, i.e., Pr{ℓ ji [k] = 1} = 1, ∀(j, i) ∈ E, ∀k ≥ 0, it was shown in [17], [9] that the components of the multi-component system can asymptotically solve average consensus and resource coordination problems in a distributed fashion by running two separate appropriately initialized linear iterations of the form where D + j (D + i ) is the out-degree of node j (i). A requirement in all cases is that the underlying communication graph (G, E) is strongly connected.

1) Average Consensus Problem:
In this problem, the nodes aim to obtain the average of the values v j , j = 1, . . . , n, they each posses. In [17], it was shown that if the initial conditions in (3) (referred to as iteration 1) are set to y j [0] = v j , and the initial conditions in (4) (referred to as iteration 2) are set to z j [0] = 1, then the nodes can asymptotically calculate v : by running the two iterations in (3) and (4).

2) Resource Coordination Problem:
In this problem, each node j can contribute a certain amount π j ≥ 0 of a given resource, which is upper and lower bounded by known capacity limits π max j and π min j respectively. The challenge is to coordinate the components so that they collectively provide a pre-determined total amount ρ d = n j=1 π j of the resource 1 as specified by an external "leader." In [9], it was shown that i) if the initial conditions in (3) are set to y j [0] = ρ d /m − π min j if j is an out-neighbor of the leader (where m ≥ 1 is the number of nodes contacted initially by the external leader) and y j [0] = −π min j otherwise, and ii) if the initial conditions in (4) are set to z j [0] = π max j − π min j , then the 1 In the development in [18], [9], it is assumed that n j=1 π min j ≤ ρ d ≤ n j=1 π max j ; this is not a restrictive assumption because in the proposed algorithms, all nodes will be able to know if ρ d < n j=1 π min j or if ρ d > n j=1 π max j (which means that no feasible solution exists). nodes can asymptotically calculate their own resource contribution π j as which satisfies In this paper, we start with a double iteration of the form in (3)-(4) that is used for either average consensus or coordination, and develop systematic methodologies to handle packet drops in the communication links.

III. ROBUSTIFICATION OF DOUBLE-ITERATION ALGORITHM
In this section, the algorithm described in Section II-B is modified so as to make it robust against communication link failures. As in (3)-(4), each node will run two iterations (which we refer to as iterations 1 and 2) to calculate quantities of interest and eventually solve the average consensus or coordination problems.
Consider the setup described in the previous section: we are given a (possible directed) strongly connected graph (G, E) representing a multicomponent system and its communication links between its components. For the sake of generality, let us refer to j as the receiving node (or receiver) and i as the sending node (or sender). For each of the two iterations, node i (the sender) will calculate several quantities of interest, which we refer to as: i) internal state; ii) total broadcasted mass; and iii) total received mass from each in-neighbor l of node i, i.e., for each node l ∈ N − i . For both iterations 1 and 2, each node j updates the value of its internal state to be a linear combination of its own previous internal state value (weighted by the inverse of the number of nodes that have j as a neighbor, i.e., 1/D + j ) and the sum (over all its in-neighbors) of the difference between the two most recently received mass values. At instant time k, the total broadcasted mass is the sum up to (and including) step k of the weighted value of node j's internal state (used to update the internal state of node j). Additionally, node j (the receiver) updates the value of the received mass from node l to be either the total broadcasted mass sent by node i if the communication link from i to j is available at instant k, or the most recently received mass value from node i, otherwise. An implicit assumption here is that messages broadcasted by node i are tagged with the sender's identity so that the receiving node j can determine where messages originated from.
For iteration 1, let y j [k] denote node j's internal state at time instant k, µ lj [k] denote the mass broadcasted from node j to each of its out-neighbors l (this is a single value ad it is the quantity is the same for each out-neighbor l of node j, i.e., for each l ∈ N − j ), and ν ji [k] denote the total mass received at node j from node i ∈ N − j . Similarly, let z j [k] denote node j's internal state takes at time instant k, σ lj [k] denote node j's broadcasted mass for each out-neighbor l, l ∈ N + j , and τ ji [k] denote the total mass received from node i ∈ N − j . Then, the progress of iteration 1 is described by where Recall that D + j (D + i ) is the number of nodes that node j (i) can transmit information to. Similarly, the progress of iteration 2 is described by where As mentioned earlier, for solving the average consensus problem, the initial conditions in (8) are set to y j [0] = v j , whereas the initial conditions in (9) are set to z j [0] = 1. Similarly, for solving the resource coordination problem, the initial conditions in (8)

Main Result:
We shall argue that with the proposed robustification strategy, despite the presence of unreliable communication links (at each time step, each link (j, i) ∈ E, fails independently from other links and independently between time steps, with some probability 1 − q ji ), nodes can asymptotically estimate the exact solution v to the average consensus by calculating, whenever z j [k] > 0 the ratio Similarly, exact solution to the resource coordination problem can be obtained as whenever z j [k] > 0. In both cases, we run the iterations in (8) and (9) and using the corresponding initial conditions as outlined above. In particular, we will show that, for every j, Additionally, we will show that z j [k] > 0 occurs infinitely often.

A. Examples
We now illustrate how the proposed algorithm works for the case of average consensus in the presence of packet-dropping communication links. We start with the rather small network shown in Fig. 1 and assume that the packets on each link (including the self-links which are not drawn in the figure 2 ) can be dropped with probability 1 − q, independently between different links and independently between different iterations. We also assume that the initial values of the five nodes are given by v = [−4, 5, 6, −3, 1] T , with their average equal to 1. Thus, in the iterations (8) and (9) We run the iterations in (8) and (9) and plot the ratio z j [k] as a function of the iteration step k for each node j (j = 1, 2, 3, 4, 5). Figure 2 shows the typical behavior that we observe for q = 0.99 (i.e., for a probability of a packet drop equal to 0.01). As can be seen in the figure, the ratio at each node quickly converges to the correct average, though the individual values for y j [k] and z j [k] do not converge.
In Figs. 3 and 4 we show typical behaviors of the same multicomponent system for q = 0.5 and q = 0.1. The behavior remains similar to the one observed before: even though y j [k] and z j [k] do not 2 We make this assumption later in the paper for the purposes of simplifying notation.    converge (in fact, they seem to behave more radically with decreasing q), the ratio z j [k] does converge to the average of the initial values. Note that the plot of the ratio in Fig. 4 is quite different than the rest: in this case, q is small enough so that z j [k] (and simultaneously y j [k]) can take the value zero (e.g., when all packets destined to node j are dropped at iteration k); thus, the ratio in such cases is not defined and is not plotted, resulting in a discontinuous set of points in the plot. Nevertheless, we can see that when packets are received (which happens frequently enough for each node), the ratio has the correct value. This is a point addressed later in the paper.
An example of what happens in larger graphs is shown in Fig. 5. Here we consider a graph with 50 nodes, randomly generated by choosing a directed edge from node i to node j, 1 ≤ i, j ≤ 50, i = j, independently with probability 1/2, and ensuring that the resulting graph is strongly connected. As can be seen the behavior remains similar to what we observed for the smaller graph: the ratio

IV. FIRST AND SECOND MOMENT ANALYSIS
In this section, we obtain recurrence relations that describe the first and second moment of the iterations after (8) and (9); this analysis is used in Section V to establish the claims in (10) and (11). In order to ease the moment calculations, the expressions in (8)-(9) will be rewritten more compactly in vector form. Also, in order to facilitate notation, we will allow each node j to drop the packet carrying its own previous value when updating its value. This way, node j handles its own value in the same way as its neighbors' values and notation is simplified significantly.

A. Vectorized Description of Double-Iteration Algorithm
Using the definition for the indicator variable x ji [k] given in (1) and the resulting indicator variable ℓ ji [k] given in (2), which describes the successful transmission of information from node i to node j over an existing, unreliable communication link, iterations (8) and (9) can be rewritten as and where Let A • B denote the Hadamard (entry-wise) product of a pair of matrices A and B of identical size. Then, for all k ≥ 0, iteration (12)- (14) can be rewritten in matrix form as where P = [p ji ] ∈ R n×n , with p ji = 1 D + i , ∀j ∈ N + i and p ji = 0 otherwise; M −1 = N −1 = 0; y k = y[k]; U ∈ R n×n , with [U ji ] = 1, ∀i, j; diag(y k ) is the diagonal matrix that results by having the entries of y k on the main diagonal; and e = [1, 1, . . . , 1] ′ (note that U = ee T ). Similarly, for k ≥ 0, (15)-(17) can be rewritten in matrix form as where S −1 = T −1 = 0, z k = z[k], and diag(z k ) is the diagonal matrix that results by having the entries of z k on the main diagonal. By defining A k := M k − N k−1 and B k := S k − T k−1 , iteration (18)- (20) can be rewritten more compactly as and iteration (21)-(23) as For analysis purposes, each matrix in (24)- (25) and (26)-(27) will be rewritten in vector form by stacking up the corresponding columns. 3 Then, (24)- (25) and (26)-(27) can be rewritten in vector form as follows. Let F = [I n I n . . . I n ] ∈ R n×n 2 , where I n is the n × n identity matrix, andP = [E 1 P T E 2 P T . . . E n P T ] T ∈ R n 2 ×n , where E i ∈ R n×n has E i (i, i) = 1 and all other entries equal zero.
[The entries of E i P T ∈ R n×n (P E T i = P E i ) are all zero except for the i th row (column) entries, which are those of the i th row (column) of matrix P T (P ).] Then, (24)-(25) can be rewritten as where a k ∈ R n 2 , x k ∈ R n 2 , and x k−1 ∈ R n 2 result from stacking the columns of matrices A k , X k , and X k−1 , respectively. Similarly, (26)- (27) can be rewritten as where b k ∈ R n 2 results from stacking the columns of matrix B k . Remark 1: It is important to note that matrices A k and B k , and their corresponding vectors a k and b k , have some entries that remain at zero for all k ≥ 0. Specifically, the (j, i) entry of matrices A k and B k (and their corresponding entries in a k and b k ) remain zero if there is no communication link from node i to node j, i.e., (j, i) / ∈ E. The reason we keep these entries (despite the fact they are zero and do not play a role in the analysis) is because they facilitate matrix notation and calculations in subsequent developments.
Since it will appear later at several points of the analysis, it is worth noting that when premultiplying P by F , we recover the matrix P , i.e.,

B. First Moment Analysis
In this section, we describe the first moment dynamics of (28)-(31) via discrete-time transition systems with no inputs, where (as shown below) the corresponding transition matrices (which are obtained from P and q) are column stochastic and primitive. In both iterations, the sum of the entries of the first moment vectors for y k and z k is shown to remain constant over time and be respectively equal to the sum of q i y 0 (i) and q i z 0 (i). Furthermore, both first moments E[y k ] and E[z k ] are shown to reach a steady-state value as k goes to infinity. The above discussion is formalized in the following lemma.
Lemma 1: Let a k , b k , y k , and z k be described by the recurrence relations in (28)- (29), and (30)-(31) respectively. Let the first moments of a k , y k , b k , and z k (i.e., E[a k ], E[y k ], E[b k ], and E[z k ]) be denoted by a k , y k , b k , and z k respectively. Then the evolution of a k , y k , b k , and z k , ∀k ≥ 1, is governed by where I m is the m × m identity matrix, with a 0 =P y 0 , y 1 = qP y 0 , b 0 =P z 0 , and z 1 = qP z 0 . Proof: Since the development for obtaining a k and y k is parallel to that for obtaining b k and z k , our analysis focuses on the first case. For k = 0 in (28)-(29), by taking expectations of both sides and noting that packet drops at time step k = 0 are independent of the initial values for a 0 , it follows that Substituting (37) into (38), we obtain y 1 = qFP y 0 = qP y 0 . For k ≥ 1 in (28)- (29), noting that packet drops at time step k are independent of previous packet drops and the initial values of a 0 , it follows, by taking expectations on both sides, that Substituting (40) into (39), we obtain Similarly, substituting (39) into (40), we have where I n is the n × n identity matrix.

C. Second Moment Analysis
In the order to calculate the second moment dynamics for (28)-(31), we utilize in the following lemma.
Lemma 2: Let x, c and d be random vectors of dimension n. Furthermore, assume that the entries of x are Bernoulli i.i.d. random variables such that Pr{x i = 1} = q and Pr{x i = 0} = 1 − q, ∀i = 1, 2, . . . n, and are independent from c and d. Then where diag(cd T ) is a diagonal matrix with the same diagonal as matrix cd T . Proof: The (i, j), i = j, entry of S can be obtained as follows: Since x i and x j are pairwise independent, and independent from c and d, it follows that . . , n, we obtain the corresponding entry of S as From these observations, the result in (46) follows. Similarly, the (i, j), i = j, entry of T can be obtained as follows: Since x i and (1 − x j ) are independent, it follows that For i = j, and observing that E[x i (1 − x i )] = 0, ∀i = 1, . . . , n, the corresponding entry of T can obtained as follows; The result in (47) follows from (52) and (53).
The following lemma establishes that the evolution of , and can be expressed as linear iterations with identical dynamics but different initial conditions. Similarly, the evolution of E[y k y T k ], E[z k z T k ], and E[y k z T k ] can also be expressed as linear iterations with identical dynamics but different initial conditions. Lemma 3: Consider the second moments of a k , y k , b k , and z k , and let , and E[y k z T k ]) be denoted by Γ k , Φ k , Ψ k , Λ k , Ξ k , and Υ k respectively. Then, the evolutions of Γ k , Φ k , Ψ k , Λ k , Ξ k , Υ k , ∀k ≥ 1, are described by the following iterations (where all I denote n 2 × n 2 identity matrices): with initial conditions Proof: The derivation of (54), (55), (60), and (61), is the same as the derivation of (56), (57), (62), and (63), thus the developments in the proof will only address the former. For k = 0, it follows from Lemma 1 and (28) that a 0 =P y 0 . Then, and Applying the results in Lemmas 1 and 2 to (67), it follows that where we used the fact that FP = P refer to Eq. (32) . For k ≥ 1, and taking into account that y k = F (a k−1 • x k−1 ), it follows that Then, from Lemma 2, (69) can be rewritten as By re-arranging terms in (70) and observing that Γ k−1 = E a k−1 a T k−1 and diag(Γ k−1 ) = E diag(a k−1 a T k−1 ) , the result in (54) follows.
Additionally, from Lemma 2, it follows that To obtain the iterations for Ξ k and Υ k , the developments are very similar to the ones above. For k = 0, and where again we used the fact that FP = P refer to Eq. (32) . For k ≥ 1, from Lemma 2 and (28), and taking into account that y k+1 = F (a k • x k ) and z k+1 = By re-arranging terms in (74) and observing that , the result in (58) follows. Finally, which completes the proof.
Although omitted in the statement of Lemma 3, it is easy to see that the dynamics of ∆ k = E[b k a T k ] and Θ k = E[z k y T k ] can also be obtained by noting that ∆ k = Ψ T k and Θ k = Υ T k .

V. CONVERGENCE ANALYSIS OF ROBUSTIFIED DOUBLE-ITERATION ALGORITHM
The previous Section established that the iterations governing the evolution of Γ k , Ψ k and Ξ k are identical except for the initial conditions. We will show next that the steady-state solutions of these iterations are also identical up to a multiplicative constant. To see this, we will rewrite (54), (56), and (58) in vector form using Kronecker products. For given matrices C, A, and B of appropriate dimensions, the matrix equation C = AXB (where X is an unknown matrix) can be rewritten as a set of linear equations of the form (B T ⊗ A)x = c, where x and c are the vectors that result from stacking the columns of matrices X and C respectively, and ⊗ denotes the Kronecker product 4 of matrices [26]. Let γ k be the vector that results from stacking the columns of Γ k andγ k the vector that results from stacking the columns of diag(Γ k ). Then, it can be easily seen that (54) can be rewritten as Let G be a diagonal matrix with entries G (l − 1)n 2 + l, (l − 1)n 2 + l = 1, ∀l = 1, 2, . . . , n 2 , and zero otherwise. Then, the second term on the right hand side of (76) can be written as which leads us to Let ψ k and ξ k and δ k be the vectors that result from stacking the columns of Ψ k , Ξ k and ∆ k respectively. Then, it is easy to see that the same recurrence relation as in (78) governs the evolution of ψ k and ξ k . Theorem 1: Let P ∈ R n×n be a column stochastic and primitive weight matrix associated with a directed graph G = {V, E}, with V = {1, 2, . . . , n} and E ⊆ V × V. Let F = [I n I n . . . I n ] ∈ R n×n 2 , where I n is the n × n identity matrix, andP = [E 1 P T E 2 P T . . . E n P T ] T ∈ R n 2 ×n , where each E i ∈ R n×n , i ∈ {1, 2, . . . , n}, satisfies E i (i, i) = 1 and has all other entries equal to zero. Then, for any q, 0 < q ≤ 1, the matrix Π defined as is column stochastic, and it has a single eigenvalue of maximum magnitude at value one. Proof: We show first column stochasticity of matrix Π. Let C = qP F + (1 − q)I and D = I −P F , so that Π = C ⊗C +q(1−q)(D ⊗D)G. We will establish that C ⊗C is column stochastic and also show 4 The Kronecker product of matrices A = [aij ] ∈ R m×n and B = [bij ] ∈ R p×q is defined (see, e.g., [26]) as the block matrix that the column sums of D ⊗D are all zero. By construction, the entries of the i th column ofP ∈ R n 2 ×n are all zero, with the possible exception of the ones indexed by (i − 1)n + j, i , i, j = 1, 2, . . . , n, each of which corresponds to the (j, i) entry of matrix P . Then, it follows that n 2 l=1P li = n j=1 P ji = 1, ∀i = 1, 2, . . . , n 2 . The matrixP F ∈ R n 2 ×n 2 is also column stochastic by construction, as it results from horizontally concatenating n times the matrixP , i.e.,P F = [PP . . .P ]; therefore, the matrix C is also column stochastic. The kronecker product of C with itself, results in an n 4 × n 4 block matrix of the form Then, it follows that the sum of the entries of the l th column of C j is n 4 m=1 C j (m, l) = ( n 2 i=1 c ij )( n 2 r=1 c rl ). Since n 2 i=1 c ij and n 2 r=1 c rl are the sum of the entries of the j th and l th columns of C = qP F + (1 − q)I n 4 (which is column stochastic), it follows that n 4 m=1 C j (m, l) = 1; therefore, C ⊗ C is also column stochastic. SinceP F is column stochastic, the column-sums of D = I −P F are zero. The kronecker product of D with itself is of the form Using similar arguments as above, it follows that n 4 m=1 D j (m, l) = ( n 2 i=1 d ij )( n 2 r=1 d rl ) = 0, which implies that the column-sums of D ⊗ D are zero. The only thing left to establish that Π is column stochastic is to show that all entries of Π are nonnegative (from where it immediately follows that Π = C ⊗ C + q(1 − q)(D ⊗ D)G is column stochastic). We argue nonnegativity of Π as follows: due to the sparsity structure of G in (79), the only nonzero entries of (D ⊗ D)G will be in columns (k−1)n 2 +k, k = 1, 2, . . . , n 2 ; thus except for entries in these columns, the entries of Π will be identical to the corresponding entries in C ⊗ C. From the structure ofP F , entries of C ⊗ C and q(1 − q)(D ⊗ D) can, respectively, take one of the following three forms: or or where p ij ≥ 0 and p lm ≥ 0 are the (i, j) and (l, m) entries of matrix P . For (80) and (81), the corresponding entry of Π is of the form and satisfies 0 ≤ qp ij p lm + (1 − q) ≤ 1 For (82) and (83), the corresponding entry of Π is of the form and satisfies 0 ≤ qp ij p lm ≤ 1. For (84) and (85), the corresponding entry of Π is of the form and satisfies 0 ≤ qp ij p lm ≤ 1.
To prove the second assertion, we will show first that matrixP F can be written via a permutation of its indices in the form where U is an irreducible column stochastic matrix and lim k→∞ W k = 0. SinceP F is column stochastic, we can assume that it corresponds to the weight matrix of some graphG = {Ṽ,Ẽ}. We will show that this graph has a single recurrent class plus a few transient states, from which the decomposition ofP F in (89) follows. Let From the structure ofP F , it follows that for any node (i, j) ∈Ṽ, one-step transitions out of (i, j) are to nodes of the form (m, i), with i ∈ N − m , where N − m is the set in-neighbors of node m in the graph G (with weight matrix P ). From the structure ofP F , it also follows that there are possibly several rows ofP F with all entries equal to zero, which means that a node (i, j) that is associated with such row cannot be reached from any other node; however, as already argued, from nodes of the form (i, j), it is possible to reach nodes of the form (m, i), where i ∈ N − m . Clearly, the nodes corresponding to rows with all entries being zero are transient. Note that the possibility of individual nodes that cannot be reached from any other node being disconnected is ruled out as it is easy to see the only nonzero diagonal entries ofP F correspond to diagonal entries of P , which are strictly smaller than one.
Next we will show that from a node (i, j) whose corresponding row inP F has some nonzero entries one can reach any other node (m, l) whose corresponding row inP F has some nonzero entries. This means that all non-transient nodes form a single recurrent class (as already argued all nonzero diagonal entries are strictly smaller than one which means there cannot be absorbing nodes). This follows from the fact that the graph G is strongly connected, which means that for any j, l ∈ V, there exists a path between j and l. Let i 1 , i 2 , . . . , i t denote the nodes traversed along the path between j and l. We will show next that for any two non-transient nodes (i, j), (r, l) ∈Ṽ there exists a path. As already argued, from (i, j) one can reach in a single hop any node of the form (m, i), where m is a neighbor of node i in the graph G. Since i 1 is the first node traversed in the path between j and l, it follows that (i 1 , i) ∈Ṽ can be reached in one step from (i, j). By repeatedly using this argument, it follows that the sequence of nodes (i 1 , i), (i 2 , i 1 ), . . . , (i t , i t−1 ), (r, i t ) forms a path between (i, j) and (r, l), which means that any non-transient node can be reached by any other non-transient node; thus, the set of non-transient nodes forms a single recurrent class. Clearly, the vertex setṼ can be decomposed into a single recurrent class and possibly several transient nodes. By re-ordering the nodes, it follows thatP F can be rewritten as in (89) (see, e.g., [27, p. 126]). Furthermore, since Q in (89) is irreducible, it follows that qQ + (1 − q)I (where I is the identity matrix) is primitive. It follows that C = qP F + (1 − q)I has a unique largest eigenvalue of value one, i.e., λ 1 = 1, and 1 > |λ 2 | ≥ · · · ≥ |λ n 2 |. Let σ(C) = {λ 1 , λ 2 , . . . , λ n 2 }. Then, σ(C ⊗ C) = {λ i λ j , i = 1, . . . , n, j = 1, . . . , n}, including algebraic multiplicities in both cases [26, p. 245]. Since λ 1 = 1 is unique (multiplicity one) and |λ i | < 1, i = 2, . . . , n 2 , it follows that the eigenvalue of C ⊗ C = [qP F + (1 − q)I] ⊗ [qP F + (1 − q)I] of largest magnitude also takes value 1 and is unique. Since C ⊗ C is column stochastic and λ 1 = 1 is unique, we know that either C ⊗ C is also primitive or it can be decomposed following a permutation of indices to the form [27, p. 126]: where L is a primitive matrix and lim k→∞ N k = 0.
We will show next that Π = C ⊗ C + q(1 − q)(D ⊗ D)G has exactly the same nonzero entries as C ⊗ C and therefore can be decomposed following the same permutation of indices to the form in (91).
As argued before, due to the sparsity structure of G in (79), the only nonzero entries of (D ⊗ D)G will be in columns (k − 1)n 2 + k, k = 1, 2, . . . , n 2 , thus except for entries in the aforementioned columns, the nonzero entries of Π will be the same as those in C ⊗ C. For all other columns in Π (that include nonzero entries in (D ⊗ D)G), it was shown in (86)-(88) that the nonzero entries of Π are strictly positive, from where it follows that Π has the same sparsity structure as C ⊗ C, which means that Π can also be decomposed in the form of (91) (for some matrices L ′ , M ′ , N ′ ), and the resulting upper-right block is also a primitive matrix. Therefore, Π has a unique largest eigenvalue at one.
The following two lemmas establish that the first and second moments of a k and b k , and y k and z k converge to the same solution up to a scalar multiplication. These two lemmas will be used to show that as k → ∞, the random vector v k = z k − αy k , for α = n j=1 z 0 (j) n j=1 y 0 (j) , will converge almost surely to v = 0. This suggests that, as k → ∞, and whenever z k is nonzero, each node i can obtain an estimate of α = n j=1 z 0 (j) n j=1 y 0 (j) by calculating the ratio y k (i)/z k (i). We will also show that, in fact, z k will be larger than some threshold infinitely often.

Lemma 4:
The first moments of a k and b k (also y k and z k asymptotically converge to the same solution up to scalar multiplication: where α = n j=1 z 0 (j) n j=1 y 0 (j) .
Proof: In Lemma 1, it was shown that y k+1 = qP +(1−q)I y k and z k+1 = qP +(1−q)I z k with y 1 = qy 0 , and z 1 = qz 0 . Since P is column stochastic and primitive, it follows that [qP +(1−q)I is also column stochastic and primitive. Thus, lim k→∞ z k = α lim k→∞ y k , where from the column stochasticity property it follows that n j=1 z k (j) = q( n j=1 z 0 (j)) and n j=1 y k (j) = q n j=1 y 0 (j) , ∀k ≥ 1; this implies that α = n j=1 z 0 (j) n j=1 y 0 (j) , which establishes (92). By noting that a 0 =P y 0 , and b 0 =P y 0 , and using the fact thatP is column stochastic, it follows that n j=1 a k (j) = n j=1 a 0 (j) = n j=1 y 0 (j) and n j=1 b k (j) = n j=1 b 0 (j) = n j=1 z 0 (j). Since qP F + (1 − q)I (i.e., the matrix that governs the dynamics of a k and b k ) is column stochastic and, as shown in the proof of Theorem 1, has a single largest eigenvalue at value 1, a similar development to the one above can be used to show (93).
Proof: The result follows from the first Borel-Cantelli lemma [28,Theorem 7.3.10]. For all k ≥ 0 and all ǫ > 0, define the event E k (ǫ) = { v k ∞ > ǫ}. We will first establish an upper bound on ∞ k=0 Pr{E k (ǫ)} by noting that Pr{E k (ǫ) ])) 1/2 = (trace(E[z k z T k ]) + α 2 trace(E[y k y T k ]) − 2αtrace(E[y k z T k ]) 1/2 . We will next show that E [ v k 2 ] → 0 as k → ∞ geometrically fast. Using Lemma 3, it can be established that E[v k v T k ] = E[z k z T k ]+α 2 E[y k y T k ]− α(E[y k z T k ] + E[z k y T k ])) = F q 2 X k−1 + q(1 − q)diag(X k−1 ) F T where X k−1 = E[w k w T k ] as defined in Lemma 5, thus the evolution of E[v k v T k ] is governed by the evolution of X k−1 or by χ k−1 (the vector that results from stacking the columns of X k−1 ). In Theorem 1, we showed that Π has a unique eigenvector (with all entries strictly positive) associated to the largest eigenvalue λ 1 = 1. Then, the solution of (94) is unique and equal to this eigenvector (up to scalar multiplication). Since Π is a column stochastic matrix, and Lemma 5 established that n 4 l=1 χ 0 (l) = 0, it follows that n 4 l=1 χ k (l) = 0, k ≥ 0, and therefore lim k→∞ χ k (l) = 0, ∀l. Additionally, it is well-known that the convergence of (94) is geometric with a rate of convergence given by the eigenvalue λ 2 of Π with the second largest modulus, which satisfies |λ 2 | < λ 1 = 1 (see, e.g., [27]). Thus, we have established that χ k (l) → 0, ∀l, geometrically fast, from where it follows that all the entries of E[v k v T k ] go to zero also geometrically fast. Therefore, the trace(E[v k v T k ]) also goes to zero geometrically fast, so that E [ v k 2 ] also goes to 0 geometrically fast. It immediately follows that ∞ k=0 E [ v k 2 ] < ∞ and therefore ∞ k=0 Pr{ v k ∞ ≥ ǫ} < ∞. Then, from the first Borel-Cantelli lemma Pr{ v k ∞ ≥ ǫ infinitely often} = 0 (or Pr{ v k ∞ ≥ ǫ i.o.} = 0).
[This follows from the fact that M 0 = n j=1 z j [0] = n and the observation that which is equal to M k .] The definition of M k+1 in (99) involves the summation of n + m nonnegative quantities, namely, z j [k + 1] for j = 1, 2, ..., n and m ji [k + 1] := (σ ji [k] − τ ji [k − 1])(1 − x ji [k]) for (j, i) ∈ E. We can think of these quantities as follows: z j [k + 1] is the mass at node j, whereas m ji [k + 1] is the mass waiting to get transferred to node j from node i. Since all of these quantities are nonnegative, at least one of them is larger or equal to n n+m . Regardless of whether this quantity is associated with a node (say node j * ) or a link (say link (j * , i * )), this mass has at least one way of reaching any node i of interest in graph G via a path of length at most n − 1 (because the graph G is strongly connected): in particular, there is at least one path of length at most n − 1 from node j * to node i and all the links in this path have weight at least 1 D + max . If all these links are activated, which occurs with probability q n−1 (q n in the case of link (j * , i * ) because the mass needs to first transfer to j * ), then a fraction ( 1 D + max ) n−1 of the mass will transfer to node i in at most n steps. Then, since for every j, z j [k] ≥ 0, k ≥ 0, independently of the values of z j [ln], l = 1, 2, . . . , k, Pr{z j [(k + 1)n] ≥ C | ζ j k , ζ j k−1 , . . . , ζ j 1 } ≥ q n obtains, whatever ζ 1 , ζ 2 , . . . , ζ k−1 . Finally, for every j , Pr{z j [(k + 1)n] ≥ C | ζ j k , ζ j k−1 , . . . , ζ j 1 } = 1 − Pr{z j [(k + 1)n] < C | ζ j k , ζ j k−1 , . . . , ζ j 1 } ≤ 1−P r{z j [(k +1)n] = 0 | ζ j k , ζ j k−1 , . . . , ζ j 1 } ≤ 1−q D − j , where D − j is the in-degree of node j.
Given a sequence of events E 1 , E 2 , . . . , E n , . . . defined on some probability space, the next theorem (which we do not prove) states the 1912 Borel criterion for establishing whether the event that infinitely many of the E k occur, denoted by {E k i.o}, will occur with probability one or zero (see, e.g., [29], [30]). This result, together with the result in Theorem 3 will be used to establish that, for every j, the event E j k = {z j [kn] ≥ C}, k ≥ 1 occurs infinitely often. Theorem 4: Let {E k }, k = 1, 2, . . . , be a sequence of events defined on some probability space. Let