Polynomial Filtering for Fast Convergence in Distributed Consensus

P olynomia l Filtering for F ast Co n v ergence i n Dis tributed C onsensus Effrosyni Kokiop oulou and P ascal F ro ssard ∗ Ecole P olytec hnique F´ ed ´ erale de Lausanne (EPFL) Signal Pro cessing Lab ora tory (L TS4) Lausanne - 10 15, Switzerland { effrosyni.k okiopoulou,pascal.frossard } @epfl.ch Abstract In the past few years, the problem of distributed consensus has received a lot of atten tion, particularly in the framework of ad ho c senso r net works. Mo s t metho ds prop ose d in the liter- ature address the consensus a veraging problem b y distributed linear iterative algorithms, with asymptotic conv erge nce o f the consensus s olution. The conv erg ence rate of such distributed algorithms typically dep ends on the netw or k top o lo gy and the weights g iven to the edg es b e- t ween neighboring sens o rs, as describ ed by the netw ork matrix. In this pap er, we prop o se to accelerate the conv er gence rate for given ne tw o rk ma trices by the use of p o lynomial filtering algorithms. The main idea o f the prop o sed metho dolog y is to apply a p oly no mial filter on the net work ma trix that will sha p e its sp ectrum in order to increase the c o nv ergenc e rate. Such an alg o rithm is equiv alent to p erio dic up dates in each of the sensors by a ggre g ating a few o f its previo us estimates . W e formulate the co mputation of the co efficients of the optimal p oly- nomial as a semi- de finite prog ram that can be efficient ly and globally solved fo r both sta tic and dynamic net work top olog ies. W e fina lly provide simulation r esults that demonstra te the effectiveness of the pr op osed solutions in acc e lerating the c onv erge nce of distributed co nsensus av erag ing problems. 1 In tro d uction W e consider the p roblem of distribu ted consensus [1] that h as b ecome recen tly v ery interesti n g esp ecially in the con text of ad ho c sens or netw orks. In particular, the pr ob lem of distribu ted a v erage consensus has attracte d a lot of researc h efforts due to its n umerous applications in dive r se areas. A few examples in clude distribu ted estimatio n [2], distribu ted compression [3], co ordination of net wo r ks of autonomous agen ts [4] and computation of a v erages and least-squares in a distribu ted fashion (see e.g., [5, 6, 7, 8] and references therein). In general the main goal of distr ibuted consensus is to r eac h a global solution using only lo cal computation and comm un ication while sta ying r ob u st to c hanges in the netw ork top ology . Giv en the initial v alues at the sensors, the p r oblem of d istributed a ve r aging is to compu te their a ve r age at e ach sensor using d istributed linear iterations. Each distrib uted iteration inv olv es lo cal comm u nication among the sensors. In particular, eac h sensor up dates its own lo cal estimate of th e a verage b y a w eight ed linear com bination of the corresp onding estimates of its neigh b ors. The we ights that are ∗ This work has b een partly supp orted by the Swiss NSF, und er NCCR IM2. 1 represent ed in a netw ork weigh t matrix W t ypically d riv e the imp ortance of the measurements of the different neigh b ors. One of the imp ortan t c h aracteristics of the distrib uted consensus algorithms is the rate of con v er gence to the asymptotic solution. In man y cases, the a v erage consensus solution can b e reac hed by s u ccessiv e multiplicatio n s of W with the v ector of initial sensor v alues. F urthermore, it has b een sho wn in [5] that in the case of fixed n etw ork top ology , the con ve r gence r ate dep ends on the s econd largest eigen v alue of W , λ 2 ( W ) . I n particular, the conv ergence is faster when the v alue of λ 2 ( W ) is small. Simila r con ve r gence results ha ve b een prop osed recen tly in the case of random net wo r k top ology [9, 10], where the conv ergence rate is gov erned by the exp ected v alue of the λ 2 ( W ) , E [ λ 2 ( W ) ]. The main researc h direction so far fo cuses on the computation of the optimal weigh ts W that yield the fastest conv ergence rate to the consensu s solution [5, 6, 7]. In this work, we div erge from metho ds that are b ased on successiv e m u ltiplications of W , and we r ather allo w the sensors to u se their previous estimates, in order to accelerate th e conv ergence r ate. This is similar in s p irit to the w orks prop osed [12, 13] that reac h the consensus solution in a fin ite num b er of steps. They use r esp ectiv ely extrap olation metho ds and linear dynamical sy s tem form ulation for fixed netw ork top ologies. In order to add ress more generic netw ork top ologies, we prop ose here to u se a matrix p olynomial p app lied on the w eight matrix W in order to shap e its sp ectrum. Give n the fact that the conv ergence rate is drive n by λ 2 ( W ) , it is ther efore p ossible to impact on the con ve r gence rate b y careful design of the p olynomial p . In the implemen tation viewp oin t, working with p ( W ) is equiv alen t to eac h sen s or aggregating its v alue p erio d ically usin g its o wn pr evious estimates. W e further form ulate the problem of the computation of the optimal co efficien ts for b oth static and dynamic netw ork top ologies. W e show th at this problem can b e efficien tly and globally solved in b oth cases by the definition of a semi-defin ite program (SDP). The rest of this pap er is organized as follo ws. In Section 2 we review the main conv ergence results of a v erage consensus in b oth fix ed and dyn amic rand om netw ork top ologies. Next, in Section 3 we introdu ce the p olynomial filtering metho dology and discuss its im p lemen tation for d istributed consensus problems. W e compute the optimal p olynomial filter in Section 4 for b oth static and dynamic netw ork top ologies. In S ection 5 we p ro vide sim ulation results that v erify the v alidit y and the effectiv eness of our m etho d. Related work is finally p resen ted in S ection 6. 2 Con v ergence in Distributed Consensus Av eraging Let us first defi ne formally the pr oblem of distribu ted consensu s a ve r aging. Assum e that in itially eac h sensor i rep orts a scalar v alue x 0 ( i ) ∈ R . W e d enote by x 0 = [ x 0 (1) , . . . , x 0 ( n )] ⊤ ∈ R n the v ector of initial v alues on the n et wo r k. Denote by µ = 1 n n X i =1 x 0 ( i ) (1) the a verag e of the initial v alues of the sensors. Ho w ever, one rarely h as a complete view of the net wo r k. The p roblem of distribu ted a verag in g therefore b ecomes t yp ically to compute µ at e ach sensor b y d istributed linear iterations. In what follo ws we review the main con vergence results for distributed consensus algorithms on b oth fix ed and random net work top ologies. 2 2.1 Static netw ork top ology W e mo del the static n et wo r k top ology as an undir ected graph G = ( V , E ) with no d es V = { 1 , . . . , n } corresp ondin g to sen sors. An edge ( i, j ) ∈ E is d ra wn if and only if sensor i can communicat e with sensor j . W e den ote the set of n eigh b ors for no d e i as N i = { j | ( i, j ) ∈ E } . Unless otherwise stated, w e assume that eac h graph is sim p le i.e., no lo ops or m u ltiple edges are allo wed. In this work, we consider distributed linear iterations of the follo wing form x t +1 ( i ) = W ii x t ( i ) + X j ∈N i W ij x t ( j ) , (2) for i = 1 , . . . , n , where x t ( j ) represent s the v alue computed by sensor j at iteration t . Since the sensors comm un icate in eac h iteration t , w e assu me th at they are synchronized. The p arameters W ij denote the edge we ights of G . Since eac h sensor comm un icates only with its d irect neigh b ors, W ij = 0 w hen ( i, j ) / ∈ E . The ab ov e iteration can b e compactly written in the follo wing form x t +1 = W x t , (3) or more generally x t =  t − 1 Y i =0 W  x 0 = W t x 0 . (4) W e call the matrix W that gathers the edge weig hts W ij , as the weigh t matrix. Note that W is a sparse matrix whose sparsit y p attern is dr iv en b y the netw ork top ology . W e assume that W is symmetric, and we denote its eigen v alue decomp osition as W = Q Λ Q ⊤ . Th e (real) eigen v alues can furth er b e arr anged as follo ws: λ 1 ( W ) ≥ λ 2 ( W ) ≥ . . . ≥ λ n − 1 ( W ) ≥ λ n ( W ) . (5) The distribu ted linear iteration giv en in eq. (3) con verge s to the a verag e if and only if lim t →∞ W t = 11 ⊤ n , (6) where 1 is the vec tor of ones [5]. In deed, notice that in this case x ∗ = lim t →∞ x t = lim t →∞ W t x 0 = 11 ⊤ n x 0 = µ 1 . It has b een shown that for fix ed netw ork top ology the con vergence rate of eq. (3) dep en d s on the magnitude of th e second largest eigen v alue λ 2 ( W ) [5]. The asymp totic con vergence factor is defined as r asym ( W ) = sup x 0 6 = µ 1 lim t →∞  k x t − µ 1 k 2 k x 0 − µ 1 k 2  1 /t , (7) and the p er-step con v ergence factor is written as r step ( W ) = sup x 0 6 = µ 1 k x t +1 − µ 1 k 2 k x t − µ 1 k 2 . (8) F ur thermore, it has b een sho wn that the con ve r gence rate relates to the sp ectrum of W , as giv en by the follo wing theorem [5]. 3 Theorem 1 The c onver genc e given by e q. (6) is gu ar ante e d if and only if 1 ⊤ W = 1 ⊤ , (9) W 1 = 1 , (10) ρ  W − 11 ⊤ n  < 1 , (11) wher e ρ ( · ) denotes the sp e ctr al r adius of a matrix. F urthermor e, r asym ( W ) = ρ  W − 11 ⊤ n  , (12) r step ( W ) = k W − 11 ⊤ n k 2 . (13) According to the ab o ve theorem, 1 √ n 1 is a left and r igh t eigen v ector of W asso ciated with the eigen v alue one, and the magnitude of all other eigen v alues is strictly less th an one. Note finally , that since W is symm etric, the asymptotic con v ergence factor coincides with the p er-step conv ergence factor, which implies that th e relations (12) and (13) are equiv alen t. W e give no w an alternate pr o of of the ab o ve theorem that illustrates the imp ortance of the second largest eigen v alue in the conv ergence rate. W e expand the initial state vect or x 0 to the orthogonal eigen basis Q of W ; th at is, x 0 = ν 1 1 √ n + n X j =2 ν j q j , where ν 1 = h 1 √ n , x 0 i and ν j = h q j , x 0 i . W e f u rther assume that ν 1 6 = 0. Then , eq. (4) implies that x t = W t x 0 = W t  ν 1 √ n 1 + n X j =2 ν j q j  = ν 1 √ n 1 + n X j =2 ν j W t q j = ν 1 √ n 1 + n X j =2 ν j λ t j q j . Observe no w that if | λ j | < 1 , ∀ j ≥ 2, then in the limit, the second term in the ab ov e equation deca ys and x ∗ = lim t →∞ x t = ν 1 √ n 1 = µ 1 . W e see th at the smaller the v alue of λ 2 ( W ) , the faster the con verge n ce rate. An alogous con ver- gence results h old in the case of dyn amic n et wo r k top ologies discuss ed n ext. 2.2 Dynamic netw ork top ology Let u s consider no w netw orks w ith rand om link failures, w here the state of a link changes ov er the iterations. In particular, we u se the r andom net wo r k mo d el pr op osed in [9, 10]. W e assume 4 that the netw ork at any arbitrary iteration t is G ( t ) = ( V , E ( t )), where E ( t ) denotes the edge set at iteration t , or equiv alen tly at time instan t t . Since the n et wo r k is dynamic, the edge set c hanges o v er the iterations, as links f ail at r andom. W e assume that E ( t ) ⊆ E ∗ , where E ∗ ⊆ V × V is the set of realizable ed ges when there is no link failure. W e also assu me th at eac h link ( i, j ) fails w ith a probabilit y (1 − p ij ), indep en d en tly of the other links. Two r andom edge sets E ( t 1 ) and E ( t 2 ) at d ifferen t iterations t 1 and t 2 are indep endent. The probabilit y of forming a particular E ( t ) is thus giv en by Q ( i,j ) ∈E ( t ) p ij . W e defin e the matrix P as P ij = ( p ij if ( i, j ) ∈ E ∗ and i 6 = j, 0 otherwise . (14) The matrix P is symmetric and its d iagonal elemen ts are zero, since it corresp onds to a simple graph. It r epresent s the probabilities of edge formation in the n et wo r k, and th e edge set E ( t ) is therefore a random subs et of E ∗ driv en by th e P matrix. Finally , the weigh t matrix W b ecomes dep end en t on the edge set since only the weigh ts of existing edges can tak e non zero v alues. In the dynamic case, the distribu ted linear iteration of eq. (2) b ecomes x t +1 ( i ) = W ii ( t ) x t ( i ) + X j ∈N i W ij ( t ) x t ( j ) (15) or in compact form, x t +1 = W t x t , (16) where W t denotes the we ight matrix corresp ond ing to th e graph realization G ( t ) of iteration t . Th e iterativ e relation giv en by eq. (16) can b e wr itten as x t =  t − 1 Y i =0 W i  x 0 . Clearly , x t no w represen ts a sto chastic pro cess sin ce the edges are dr a wn rand omly . The conv ergence rate to the consensu s solution therefore d ep end s on the b eha vior of the pro d uct Q t − 1 i =0 W i . W e say that the algorithm conv erges if ∀ x 0 ∈ R n , lim t →∞ E k x t − µ 1 k = 0 . (17) W e review no w some con verge n ce results from [10], whic h fi rst shows that Lemma 1 F or any x 0 ∈ R n , k x t +1 − µ 1 k ≤  t Y i =0 ρ  W i − 11 ⊤ n   k x 0 − µ 1 k . It leads to the f ollo wing con ve r gence th eorem [10 ] for dynamic netw orks. Theorem 2 If E [ ρ  W − 11 ⊤ n   ] < 1 , the ve ctor se quenc e { x t } ∞ t =0 c onver ges in the sense of e q. (17). 5 W e define the conv ergence factor in dynamic net work top ologies as r d ( W ) = E h ρ  W − 11 ⊤ n  i . This factor dep end s in general on the sp ectral p rop erties of the ind u ced n et wo r k matrix and drives the con v ergence r ate of eq. (16). The authors in [14] sho w that | λ 2 ( E [ W ]) | < 1 is also a necessary and suffi cien t condition for asymptotic (almost sure) conv ergence of th e consensu s algorithm in the case of random net works, wh ere b oth netw ork top ology and weigh ts are random (in particular i.i.d and indep endent ov er time). Finally , it is interesting to note that th e consensus p roblem in a r andom net work r elates to gossip algorithms. Distributed a v eraging und er th e synchronous gossip constr aint implies that m ultiple no d e p airs ma y communicate sim u ltaneously only if these n o de pairs are disjoint. In other w ords , the set of links implied by the activ e no de pairs f orms a matc hin g of the graph . T herefore, the d istr ibuted a verag ing problem describ ed ab o ve is closely related to the distr ibuted syn c hr on ou s algorithm und er the gossip constrain t that h as b een p rop osed in [15, Sec. 3.3.2]. It has b een sh o wn in this case that the av eraging time (or conv ergence rate) of a gossip algorithm d ep ends on th e second largest eigen v alue of a doub ly sto c hastic net wo r k matrix. 3 Accelerated consensus with p olynomial filtering 3.1 Exploiting memory As we ha ve seen ab o ve, the con ve r gence rate of the distributed consens u s algorithms d ep end s in general on the sp ectral prop erties of an induced net wo r k matrix. Th is is the case for b oth fixed and random net work top ologies. Most of the r esearc h wo r k has b een d ev oted to find ing w eight matrix W for accelerating the conv ergence to the consensus solution when sen s ors only u se their cur rent estimates. W e c h o ose a different approac h where we exploit the memory of sens ors , or the v alues of previous estimates in ord er to augment to con ve r gence r ate. Therefore, we ha ve prop osed in our previous work [12] the Scalar Epsilon Algorithm (SEA) for accelerating the con v ergence rate to th e consensus solution. SEA b elongs to th e family of extrap olation metho d s for accele r ating vecto r sequences, such as eq. (3). Th ese metho ds exploit the fact that th e fixed p oin t of the sequence b elongs to the su bspace spann ed by any ℓ + 1 consecutiv e terms of it, wh ere ℓ is the degree of the minimal p olynomial of the sequ ence generator matrix (for more d etails, s ee [12] and references therein). SEA is a lo w complexit y algorithm, w hic h is ideal for sensor netw orks and it is kno wn to reac h the consensu s solution in 2 ℓ steps. How ev er, ℓ is unknown in practice, so one ma y u se all th e a v ailable terms of the ve ctor sequence. Hence, the memory requirement s of S EA are O ( T ), where T is the n u m b er of terms. Moreo v er, SEA assumes that the sequence generator matrix (e.g., W in the case of eq. (3)) is fixed, so that it d o es n ot adapt easily to dynamic net work top ologies. In this pap er, w e p rop ose a more flexible algorithm based on the p olynomial filtering tec hniqu e. P olynomial filtering p ermits to “sh ap e” the sp ectrum of a certain symmetric w eight matrix, in order to accelerate th e conv ergence to the consensus s olution. Similarly to S EA, it allo ws the sensors to us e the v alue of their pr evious estimates. Ho we ver, the p olynomial filtering metho dology in tro duced b elo w p resen ts thr ee main adv anta ges: (i) it is r obust to dyn amic top ologies (ii) it has explicit con trol on the conv ergence rate and (iii) its memory requ iremen ts can b e adju sted to the memory constrain ts imp osed by the sensor. 6 3.2 P olynomial filtering Starting from a giv en (p ossibly optimal) w eight matrix W , we p r op ose the app lication of a p olyno- mial filter on the sp ectrum of W in order to impact the m agnitude of λ 2 ( W ) that mainly drives th e con v er gence rate. Denote by p k ( λ ) the p olynomial filter of degree k that is applied on the sp ectrum of W , p k ( λ ) = α 0 + α 1 λ + α 2 λ 2 + . . . + α k λ k . (18) Accordingly , the matrix p olynomial is giv en as p k ( W ) = α 0 I + α 1 W + α 2 W 2 + . . . + α k W k . (19) Observe no w that p k ( W ) = p k ( Q Λ Q ⊤ ) = α 0 I + α 1 ( Q Λ Q ⊤ ) + . . . + α k ( Q Λ k Q ⊤ ) = Qp k (Λ) Q ⊤ , (20) whic h implies that the eigenv alues of p k ( W ) are simply the p olynomial filtered eigen v alues of W i.e., p k ( λ i ( W )) , i = 1 , . . . , n . In the implemen tation lev el, w orkin g on p k ( W ) implies a p erio dic up d ate of the current sensor’s v alue with a linear com bin ation of its p revious v alues. T o see w h y th is is true, we obs erv e th at: x t + k + 1 = p k ( W ) x t = α 0 x t + α 1 W x t + . . . + α k W k x t = α 0 x t + α 1 x t +1 + . . . + α k x t + k . (21) A careful d esign of p k ma y impact the con vergence r ate dr amatically . Then, eac h sensor t yp ically applies p olynomial fi ltering for distributed consensus by follo wing the m ain steps tabulated in Algorithm 1. Note that, for b oth fixed and rand om net work top ology cases, the α k ’s are compu ted off-line assuming th at W and resp ectiv ely E [ W ] are kno wn a p riori. In what follo ws, w e pr op ose d ifferen t approac hes for computing the co efficients α i of the filter p k . 3.3 Newton’s interpolating p olynomial One simple and rather intuitiv e approac h for the design of the p olynomial p k ( λ ) is to use Hermite in terp olation. Recall that the ob jectiv e is to damp en the smallest eigen v alues λ 2 , . . . , λ n of W , while ke epin g the eigen v alue one intact . Therefore, w e assume that the sp ectrum of W lies in an in terv al [ a, 1] and w e imp ose s mo othness constrain ts of p k at the left endp oin t a . In particular, the p olynomial p k ( λ ) : [ a, 1] → R that we seek, will b e d etermined by the follo w ing constrain ts: p k ( a ) = 0 , p ( i ) k ( a ) = 0 , i = 1 , . . . , k − 1 (22) p k (1) = 1 , (23) where p ( i ) k ( a ) denotes the i -th deriv ativ e of p k ( λ ) ev aluated at a . Th e damp ening is ac hieve d b y imp osing smo othness constraint s of the p olynomial on the left end p oint of the in terv al. The 7 Algorithm 1 Po lyn omial Filtered Distribu ted Consens us 1: Input : p olynomial co efficien ts α 0 , . . . , α k +1 , tolerance ǫ . 2: Output : a v erage estimate ¯ µ . 3: Initializat ion : 4: ¯ µ 0 = x 0 ( i ). 5: S et the iteration in dex t = 1. 6: rep eat 7: if mo d ( t, k + 1) ==0 then 8: x t ( i ) = α 0 x t − k − 1 ( i ) + α 1 x t − k ( i ) + α 2 x t − k + 1 ( i ) + . . . + α k x t − 1 ( i ) { p olynomial fi ltered up date } 9: x t ( i ) = W ii x t ( i ) + P j ∈N i W ij x t ( j ). 10: else 11: x t ( i ) = W ii x t − 1 ( i ) + P j ∈N i W ij x t − 1 ( j ). 12: end if 13: Increase the iteration ind ex t = t + 1. 14: ¯ µ t = x t ( i ). 15: un t il ¯ µ t − ¯ µ t − 1 < ǫ 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 2 4 6 8 10 Figure 1: Newton’s p olynomial of v arious d egrees k . computed p olynomial will ha ve degree equal to k . Finally , the co efficien ts of p k that satisfies the ab o ve constrain ts can b e computed by the Newton’s divided differences table. Figure 1 shows the form of p k ( λ ) f or a = 0 and different v alues of the d egree k . As k increases, the damp ening of the sm all eigen v alues b ecomes more effectiv e. I t is w orth mentio n ing that the design of Newton’s p olynomial do es not d ep end on the net w ork size or the size of the w eight matrix W . What is only needed is a left endp oint a , which encloses the sp ectrum of W , as well as the desired degree k , wh ic h moreo ver ma y b e imp osed by memory constrain ts. This f eature of Newton’s p olynomial is v ery interesting and it is p articularly app ealing in the case of dynamic net wo r k top ologies. Ho wev er, the ab o ve p olynomial design is m ostly driv en b y intuitiv e arguments, whic h tend to obtain sm all eigen v alues for faster con verge n ce. In the follo wing section, w e pr o vide an alternativ e tec hniqu e for computing the p olynomial filter that optimizes the con v er gence rate. 8 4 Optimal p olynomial filters 4.1 P olynomial filtering for static netw ork top ologies W e are no w in terested in fin ding th e p olynomial that leads to the fastest con ve r gence of linear iteration describ ed in eq. (3), for a giv en we ight matrix W and a certain d egree k . F or notational ease, we call W p = p k ( W ) = P k i =0 α i W i . According to T heorem 1, the optimal p olynomial is the one that minimizes th e second largest eigen v alue of W p , i.e., ρ  W p − 11 ⊤ n  . Therefore, we n eed to solv e an optimization problem wh ere the optimization v ariables are the k + 1 p olynomial co efficient s α 0 , . . . , α k and th e ob jective f unction is the sp ectral radius of W p − 11 ⊤ n . Optimization pr oblem: OPT1 min α ∈ R k +1 ρ  P k i =0 α i W i − 11 ⊤ n  sub j ect to  P k i =0 α i W i  1 = 1 . In terestingly , the optimization pr oblem OPT1 is con vex. Firs t, its ob jectiv e function is con vex, as stated in Lemma 2 . Lemma 2 F or a giv e n symmetric weight matrix W and de gr e e k , ρ  P k i =0 α i W i − 11 ⊤ n  is a c onvex function of the p olynomial c o e fficients α i ’s. Pr o of : Let β , γ ∈ R k +1 and 0 ≤ θ ≤ 1. Since W is symmetric, P k i =0 α i W i is also sy m metric. Hence, the sp ectral radius is equal to the matrix 2-norm. Thus, w e ha ve ρ  k X i =0 ( θ β i + (1 − θ ) γ i ) W i − 11 ⊤ n  =    k X i =0 ( θ β i + (1 − θ ) γ i ) W i − 11 ⊤ n    2 =    θ  k X i =0 β i W i − 11 ⊤ n  + (1 − θ )  k X i =0 γ i W i − 11 ⊤ n     2 ≤ θ    k X i =0 β i W i − 11 ⊤ n    2 + (1 − θ )    k X i =0 γ i W i − 11 ⊤ n    2 ≤ θ ρ  k X i =0 β i W i − 11 ⊤ n  + (1 − θ ) ρ  k X i =0 γ i W i − 11 ⊤ n  , whic h p ro ves the Lemma. In addition, the constrain t of OPT1 is linear whic h implies that the set of f easible α i ’s is conv ex. As OPT1 minimizes a con v ex function o ve r a con vex set, the optimization problem is indeed conv ex. In order to solv e OPT1, we use an auxiliary v ariable s to b ound the ob jectiv e fu n ction, and then we express the sp ectral radiu s constrain t as a linear matrix inequality (LMI). Thus, we need to solv e the follo wing optimization problem. 9 Optimization pr oblem: OPT2 min s ∈ R ,α ∈ R k +1 s sub j ect to − sI  P k i =0 α i W i − 11 ⊤ n  sI ,  P k i =0 α i W i  1 = 1 . Recall that since W is sym metric, W p = P k i =0 α i W i will b e symmetric as well. Hence, th e constraint W p 1 = 1 is su fficien t to ensure that 1 w ill b e also a left eigen v ector of W p . The s p ectral r adius constrain t, − sI  W p − 11 ⊤ n  sI ensures that that all the eigen v alues of W p other than the fi rst one, are less or equal to s . D u e to the L MI, the ab o ve optimization p r oblem b ecomes equiv alen t to a semi-definite pr ogram (SDP) [16]. S DPs are conv ex problems and can b e globally and efficien tly solve d . The solution to OPT 2 is therefore computed efficien tly in pr actice, wher e the SDP only h as a mo d erate num b er of k + 2 unknowns (including s ). 4.2 P olynomial filtering for dynamic netw ork top ologies W e extend now the idea of p olynomial filtering to dynamic net work top ologies. Th eorem 2 suggests that the con v ergence rate in the ran d om net work top ology case is gov erned by E h ρ  W − 11 ⊤ n  i . Since W d ep ends on a d ynamic edge set, p k ( W ) no w b ecomes sto c hastic. F ollo wing the same in tuition as ab ov e, we could form an optimization problem, similar to OPT1, whose ob j ectiv e function w ould b e E [ ρ  P k i =0 α i W i − 11 ⊤ n  ]. Although this ob jectiv e fun ction can b e shown to b e con v ex, its ev aluation is hard and typica lly requires seve r al Mon te Carlo simulati ons steps. Recall that the con ve r gence rate of eq. (16) is related to the second largest eigen v alue of E [ W ], whic h is muc h easier to ev aluate. Let W denote the a v er age we ight matrix E [ W ]. W e then observe that E h ρ  W − 11 ⊤ n  i ≥ ρ  W − 11 ⊤ n  , whic h is du e to Lemma 2 and Jensen’s inequ alit y . Th e ab o ve inequalit y implies that the sp ectral radius ρ  W − 11 ⊤ n  shall b e small in order E h ρ  W − 11 ⊤ n  i to b e small to o. Additionally , the auth ors pro vid e exp erimenta l evidence in [10], whic h indicates th at ρ  W − 11 ⊤ n  seems to b e closely related to the conv ergence rate of eq. (16). Based on the ab o ve facts, we prop ose to build our p olynomial fi lter based on W . Hence, w e form u late the follo wing optimization problem for computing the p olynomial co efficien ts α i ’s in the random net work top ology case. Optimization pr oblem: OPT3 min α ∈ R k +1 ρ  P k i =0 α i W i − 11 ⊤ n  sub j ect to  P k i =0 α i W i  1 = 1 . 10 OPT3 could b e viewed as the analog of OPT 1 f or the case of dynamic netw ork top ology . The main difference is that we work on W , whose eigen v alues can b e easily obtained. Using again the auxiliary v ariable s , we reac h the follo w in g f orm ulation for obtaining the α i ’s. Optimization pr oblem: OPT4 min s ∈ R ,α ∈ R k +1 s sub j ect to − sI  P k i =0 α i W i − 11 ⊤ n  sI ,  P k i =0 α i W i  1 = 1 . Once W h as b een compu ted, this optimization pr oblem is solv ed efficient ly by a S DP similarly to the case of static net works. 5 Sim ulation Results 5.1 Setup In this section, w e pr o vide simulati on results wh ic h sh o w the effectiv eness of the p olynomial filtering metho dology . First w e in tro duce a few w eight matrices that hav e b een extensiv ely used in the distributed a verag ing literature. Supp ose that d ( i ) denotes the degree of the i -th sensor. It h as b een shown in [5, 6] that iterating eq. (3) with the follo wing matrices leads to conv ergence to µ 1 . • Maximum-de gr e e wei ghts. Th e Maxim um-degree weig ht matrix is W ij =      1 n if ( i, j ) ∈ E , 1 − d ( i ) n i = j, 0 otherwise . (24) • Metr op olis weig hts. Th e Metrop olis w eigh t matrix is W ij =      1 1+max { d ( i ) ,d ( j ) } if ( i, j ) ∈ E , 1 − P ( i,k ) ∈E W ik i = j, 0 otherwise . (25) • L apla c i an w eights. S upp ose that A is the adjacency matrix of G and D is a diagonal matrix whic h holds the v ertex d egrees. The Lap lacian matrix is defin ed as L = D − A and the Laplacian weig ht matrix is d efined as W = I − γ L, (26) where the scalar γ m ust satisfy γ < 1 /d max [5]. The sensor net works are built using the random geographic graph mo del [19]. In p articular, w e place n no d es u niformly distributed on th e 2-dimensional un it area. T wo no des are adjacen t if their Eu clidean distance is smaller than r = q log( N ) N in ord er to guarantee connectedness of the graph with high p robabilit y [19]. Finally , th e S DP pr ograms for optimizing the p olynomial fi lters are solv ed in Matlab using the SeDuMi [18] solv er 1 . 1 Publically av ailable at: http://sedumi.m cmaster.ca/ 11 0.65 0.7 0.75 0.8 0.85 0.9 0.95 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 (a) SDP p olynomial filter p k ( λ ) , λ min ≤ λ ≤ 1 0 10 20 30 40 50 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 eig(W) eig(P(W)) (b) Effect on the sp ectrum λ ( W ) Figure 2: Effect of p olynomial filtering on the sp ectrum of the Maximum-degree matrix of eq. (24). 5.2 Static netw ork top ologies W e illustrate fi rst the effect of p olynomial filtering on the sp ectrum of W . W e build a net work of n = 50 sensors and we apply p olynomial fi ltering on the Maxim u m-degree weig ht matrix W , give n in (24). W e use k = 4 and we solv e the optimization p roblem OPT2 using the Maxim um-degree matrix W as input. Figure 2(a) shows the obtained p olynomial filter p k ( λ ), when λ ∈ [ λ min ( W ) , 1]. Next, w e apply th e p olynomial on W and Figure 2(b) sho w s the sp ectrum of W b efore (star-solid line) and after (circle-soli d lin e) p olynomial filtering, v ersu s the vec tor index. O b serv e that p olynomial filtering dr amatically in creases the sp ectral gap | 1 − λ 2 | , whic h further leads to accelerating the distributed consensus, as w e sho w in th e simulat ions that follo w. Then we compare the p erformance of the different distribu ted consensus algorithms, with all the aforemen tioned weigh t matrices; that is, Maxim u m-degree, Metrop olis and Laplacian w eight matrices f or distributed a v eraging. W e compare b oth Newton’s p olynomial and the SDP p olynomial (obtained from the solution of O PT2) with the standard iterativ e metho d, w hic h is based on successiv e iterations of eq. (3). F or the sak e of completeness, we also pro vid e the resu lts of the scalar epsilon algorithm (SE A) that uses all p revious estimates [12]. First, we explore th e b eha vior of p olynomial filtering m etho ds un der v ariable degree k from 2 to 6 with step 2. W e use the Lap lacian w eight matrix for this exp eriment. Figures 3(a) and 3(b) illustrate the ev olution of the absolute error k x t − µ 1 k 2 v ersus the iteration in dex t , for p olynomial filtering with SDP and Newton’s p olynomials resp ectiv ely . W e also pro vide th e curve of the standard iterativ e metho d as a b aseline. Ob s erv e fir st that b oth p olynomial filtering metho ds outp erform the standard metho d by exhib iting faster con verge n ce rates, across all v alues of k . Notice also, that the degree k go v ern s the con vergence r ate, since larger k implies more effectiv e filtering and ther efore faster conv ergence. Finally , the stagnati on of the con vergence p ro cess of the SDP p olynomial filtering and large v alues of k is due to th e limited accuracy of the SDP solv er. Next, we sho w th e results obtained with the other t wo weig ht matrices on the same sensor net- w ork. Figures 4(a) and 4(b) show the con vergence b eha vior of all metho ds for the Maxim um-degree and Metrop olis matrices resp ectiv ely . In b oth p olynomial filtering m etho ds we use a representat ive 12 0 50 100 150 200 250 10 −15 10 −10 10 −5 10 0 10 5 Number of iterations ||x − µ || 2 Iter k=2 k=4 k=6 SEA (a) SDP p olynomial filtering 0 50 100 150 200 250 10 −10 10 −5 10 0 10 5 Number of iterations ||x − µ || 2 Iter k=2 k=4 k=6 SEA (b) Newton p olyn omial filtering Figure 3: Beha vior of p olynomial fi ltering for v ariable degree k on fixed top ology u sing the Laplacian w eight matrix. 0 50 100 150 200 250 10 −10 10 −5 10 0 10 5 Number of iterations ||x − µ || 2 Iter SDP−PF Newton−PF SEA (a) Max-degree wei ght matrix 0 50 100 150 200 250 10 −15 10 −10 10 −5 10 0 10 5 Number of iterations ||x − µ || 2 Iter SDP−PF Newton−PF SEA (b) Metropolis wei ght matrix Figure 4: Comparison b et we en Newton’s p olynomial and SDP p olynomial ( k = 4) on fi xed top ology . v alue of k , namely 4. Notice again that p olynomial filtering accelerates th e con verge n ce of the standard iterativ e metho d (solid line). As exp ected, the optimal p olynomial compu ted with SDP outp erforms Newton’s p olynomial, which is based on int u itiv e argumen ts only . Finally , we can see from Figur es 3 and 4 that in s ome cases the con verge n ce rate is comparable for SEA and SDP p olynomial filtering. Note ho wev er that the former uses all previous iterates, in con trast to the latter that uses only the k + 1 most recent ones. Hence, the memory requir ements are smaller for p olynomial fi ltering, since they are directly drive n by k . This moreo ver allo ws more direct con trol on the con verge n ce r ate, as we ha ve seen in Fig. 3. Inte r estingly , we see that the 13 0 50 100 150 200 250 10 −15 10 −10 10 −5 10 0 10 5 Number of iterations ||x − µ || 2 Iter SDP−PF Newton−PF (a) k = 1, q = 1 0 50 100 150 200 250 10 −15 10 −10 10 −5 10 0 10 5 Number of iterations ||x − µ || 2 Iter SDP−PF Newton−PF (b) k = 2, q = 0 . 1 0 50 100 150 200 250 10 −15 10 −10 10 −5 10 0 10 5 Number of iterations ||x − µ || 2 Iter SDP−PF Newton−PF (c) k = 2, q = 0 . 3 0 50 100 150 200 250 10 −20 10 −10 10 0 10 10 10 20 10 30 Number of iterations ||x − µ || 2 Iter SDP−PF Newton−PF (d) k = 2, q = 0 . 8 Figure 5: Random net w ork top ology sim ulations. q denotes the probabilit y that the n et wo r k c hanges at eac h iteration. con v er gence pro cess is smo other with p olynomial filtering, w h ic h f urther p ermits easy extension to dynamic net work top ologies. 5.3 Dynamic netw ork top ologies W e study n o w the p erformance of p olynomial filtering for d ynamic netw orks top ologies. W e build a sequence of r andom net wo r k s of n = 50 sensors , an d we assum e that in eac h iteration the net wo r k top ology changes indep endently from the p revious iterations, with probabilit y q and with probabilit y 1 − q it remains the same as in the previous iteration. W e compare all metho d s f or differen t v alues of the probabilit y q . W e use the Laplacian we ight matrix (26). In the SDP p olynomial fi ltering metho d , w e solve the S DP program OP T4 (see Sec. 4.2). Fig. 5 shows the a v erage p erf ormance of p olynomial filtering for some representa tive v alues of the degree k and the 14 probabilit y q . The av erage p erformance is computed using the median o v er the 100 exp eriments. W e ha ve not rep orted the p erform ance of the SEA algorithm, since it is n ot robu st to c h anges of the net work top ology . Notice th at wh en k = 1 (i.e., eac h sensor uses only its cur ren t v alue and the r ight pr evious one) p olynomial filtering accelerat es the con verge n ce o v er the s tandard metho d. A t the same time, it sta ys robust to net work top ology c hanges. Also, obser ve that in this case, the SDP p olynomial outp erforms Newton’s p olynomial. Ho we ver, wh en k = 2, the r oles b et ween the t wo p olynomial filtering metho ds c h an ge as the p r obabilit y q increases. F or instance, when q = 0 . 8, the SDP metho d ev en diverges. This is exp ected if we think that the co efficients of Newton’s p olynomial are computed using Hermite in terp olation in a giv en inte r v al and they do not dep end on the sp ecific realizatio n of the u n derlying weig ht matrix. Thus, they are m ore generic th an those of the SDP p olynomial th at takes W into accoun t, and therefore less sensitiv e to the actual top ology realization. Algorithms based on optimal p olynomial filtering b ecome inefficient in a h ighly dynamic n et wo rk , whose top ology changes very frequent ly . 6 Related w ork Sev eral w orks hav e studied the conv ergence rate of d istributed consen s us algorithms. In p articular, the authors in [5] and [9, 10] ha ve sho wn that the con ve r gence rate dep end s on the second largest eigen v alue of the n et wo r k weig ht matrix, for fixed and random n et wo r ks, r esp ectiv ely . Th ey b oth use semi-definite p r ograms to compute the optimal w eight matrix, and the optimal top ology . Other works ha ve addressed the consensus p roblem, and we mention here only the most relev an t ones. A. O lshevsky and J. N. Tsitsiklis in [8] pr op ose t wo consensus algorithms f or fixed net work top ologies, which b uild on the “agreemen t algorithm”. The prop osed alg orithm s make use of spanning trees and the authors b oun d their worst-ca se con vergence rate. F or d ynamic net wo r k top ologies, they prop ose an algorithm wh ic h build s on a previously known distrib uted load b alancing algorithm. In this case, the auth ors sh o w that the algo r ithm has a p olynomial b ound on the con v er gence time ( ǫ -con v ergence). The authors in [25] stud y the con verge n ce prop erties of agreemen t o ve r random n et wo r ks fol- lo wing the Erd˝ os and R ´ enyi random graph mo del. According to th is m o del, eac h edge of the graph exists w ith probability p , in dep en d en tly of other edges and the v alue of p is the s ame for all edges. By agreemen t, we consid er the case where all no des of the graph agree on a particular v alue. The authors employ resu lts from sto c hastic stabilit y in order to establish con verge n ce of agreemen t o ver random net wo r ks. Also, it is sho w n that the rate of con verge n ce is go v ern ed by the exp ectation of an exp onen tial factor, whic h in vol ves the second sm allest eigen v alue of th e Laplacian of the graph . Gossip algorithms hav e also b een applied successfully to solving d istributed a veragi n g prob- lems. In [15] p ro vide con ve r gence r esu lts on randomized gossip algorithm in b oth synchronous an d async hr on ou s settings. Based on the obtained results, they optimize the net wo r k top ology (edge formation pr obabilities) in order to maximize th e conv ergence rate of randomized gossip. This optimization pr oblem is also formulated as a semi-definite program (SDP). In a r ecen t study , th e authors in [20] ha ve b een able to improv e the standard gossip proto cols in cases w here the sen s ors kno w th eir geometric p ositions. The main idea is to exp loit geographic routing in ord er to aggregate v alues among random no des th at are far a wa y in the net work. Under the same assumption of knowing the geometric p ositions of th e sensors, the authors in [21] pr op ose a fast consen s us algorithm for geographic r andom graphs. In particular, they utilize 15 lo cation in formation of the s en sors in order to constru ct a nonreve r sible lifted Mark o v c hain that mixes faster th an corresp ondin g rev ersible chains. The main idea of lifting is to distinguish the graph n o des f rom the states of the Mark o v chai n and to “split” the states into virtual states that are conn ected in such a w ay that p ermits faster mixing. The lifted graph is then “pro j ected” bac k to the original graph, wh ere the dynamics of the lifted Mark o v c hain are s imulated s ub j ect to the original graph top ology . Ho wev er, the prop osed algorithm is not applicable in the case where the no des’ geographic lo cation is not av ailable. In [22] th e auth ors pr op ose a cluster-based distributed a v eraging algorithm, app licable to b oth fixed linear iteration and random gossiping. Th e induced ov erla y graph that is constructed b y clustering the no des is b etter connected relativ ely to the original graph; hence, the random walk on the o verla y graph mixes faster than the corresp onding walk on the original graph. Along th e same lines, K . J ung et al. in [23], ha ve u sed n onrev ersib le lifted Mark o v c h ains to accelerate consensus. They us e the lifting scheme of [24] and they prop ose a deterministic gossip algorithm b ased on a set of disjoint maximal matc h ings, in order to simulate the dynamics of the lifted Mark o v chain. Finally , ev en if w e hav e mostly considered sync hr onous algorithms in this pap er, it is wo rth men tioning that th e authors in [26] pr op ose t wo asynchronous algorithms for distribu ted a veragi n g. The fir st algorithm is based on b lo c king (that is, w hen tw o no des u p d ate their v alues they blo c k unt il the u p d ate has b een completed) and the other algorithm d rops the blo c kin g assumption. The authors show the con vergence of b oth algorithms und er very general async hr onous timing as- sumptions. Moreo v er, the authors in [27] prop ose c onsensus pr op agation , w hic h is an async hr onous distributed p r oto col that is a sp ecial case of b elief pr opagation. In the case of singly-connected graphs (i.e., connected with n o lo op s ), sync h ronous consensus propagation con v erges in a num b er of iterations th at is equal to the diameter of the graph. The authors provi d e con vergence analysis for regular graphs . 7 Conclusions In th is pap er, we prop osed a p olynomial filtering metho dology in order to accelerate distribu ted av- erage consensus in b oth fixed and rand om n etw ork top ologies. The main idea of p olynomial filtering is to sh ap e th e sp ectrum of the p olynomial weigh t matrix in order to minimize its second largest eigen v alue and subsequent ly increase the conv ergence rate. W e ha ve constructed semi-defin ite programs to compute th e optimal p olynomial co efficien ts in b oth static and d ynamic net works. Sim u lation results with sev eral common w eight matrices ha ve sho wn that the con ve r gence rate is m uch higher than for state-of-the-a rt algorithms in most scenarios, except in the s p ecific case of highly dynamic net works and small m emory sensors. 8 Ac kno wledgemen ts The first author w ould like to thank prof. Y ou s ef Saad for the v aluable and insigh tfu l discussions on p olynomial fi ltering. References [1] D. Bertsek as and J. N. Tsitsiklis, P ar al lel and Distribu te d Computation: Numeric al Metho ds . Pren tice Hall, 1989. 16 [2] I. D. Sc hizas, A. Rib eiro, and G. B. Giannakis, “Consensus in Ad Ho c WSNs with Noisy Links - Pa r t I: Distributed Estimation of Deterministic Signals,” IEEE T r ansactions on Singal Pr o c essing , v ol. 56, no. 1, pp. 350–364, Jan uary 2008. [3] M. R ab b at, J. Haupt, A. Singh, and R. No wak, “Decen tralized compression and pred istribution via r an d omized gossiping,” 5th ACM Int. Conf. on Information Pr o c essing in Sensor Networks (IPSN) , pp. 51 – 59, April 2006. [4] V. D. Blondel, J. M. Hendric kx, A. Olshevsky , and J. N. Tsitsiklis, “Conv ergence in multiagen t co ordination, consensus and flo c kin g,” IEEE Conf. on D e cision and Contr ol, and the Eu r op e an Contr ol Confer enc e , pp. 2996–300 0, Decem b er 2005. [5] L. Xiao and S. Bo yd, “F ast linear iterations for distribu ted a v eraging,” Systems and Contr ol L etters , n o. 53, pp. 65–78, F ebru ary 2004. [6] L. Xiao, S. Bo yd , and S. Lall, “A sc h eme for robust distributed sensor fus ion based on a verag e consensus,” Int. Conf. on Information Pr o c essing in Se nsor Networks , p p. 63–70, Apr il 2005, los Angeles. [7] ——, “Distributed av erage consensu s with time-v arying metrop olis w eights,” Automat i c a , June 2006, submitted. [8] A. O lshevsky and J. Tsitsiklis, “Con vergence rates in distribu ted consens us and a v eraging,” IEEE Confer enc e on De cision and Contr ol , Decem b er 2006, san Diego, CA. [9] S. Kar and J. M. F. Moura, “Distributed a verag e consensus in sensor n et wo r ks with random link failures,” IEEE Int. Conf. on A c oustics, Sp e e ch and Signal P r o c e ssing (ICASSP) , Apr il 2007. [10] ——, “Sen s or net works with random links: T op ology design for distribu ted consensus,” IEE E T r ansactio ns on Signal Pr o c essing , April 2007, su bmitted. [11] E. K okiop oulou and Y. S aad, “P olynomial fi ltering in laten t seman tic indexin g for informa- tion r etriev al,” 27th ACM-SIGIR Confer e nc e on R ese ar ch and Development in Information R e triev al , 2004. [12] E. Kokiop ou lou and P . F r ossard, “Accelarating distributed consensu s using extrap olation,” IEEE Signal Pr o c essing L etters , v ol. 14, no. 10, p p. 665–668 , Octob er 2007. [13] S. Sun daram and C. N. Hadjicostis, “Distributed consensus and linear functional calculatio n in n et wo r ks: An obser v abilit y p ersp ectiv e,” 6th ACM Int. Conf. on Information Pr o c essing i n Sensor Networks (IPSN) , Apr il 25-27 2007. [14] A. T ah baz-Salehi and A. J ad b abaie, “On consensu s o v er r andom net wo r ks,” In Pr o c e e dings of 44th annual Al lerton Confer enc e on Communic ation, Contr ol and Computing , pp. 1315–132 1, Septem b er 2006. [15] B. P . S . Bo yd, A. Ghosh and D. Shah, “Randomized gossip algorithms,” IEE E T r ansactio ns on Information The ory , vo l. 52, pp . 2508–2530 , Ju n e 2006. [16] S. Bo yd and L. V andenberghe, Convex Optimization . Cam br idge Unive r sit y Press, 2004. 17 [17] Y. Saad, Iter ative metho ds for sp arse line ar systems , 2nd ed. SIAM, 2003. [18] J. F. Sturm , “Imp lemen tation of in terior p oin t metho ds for mixed semidefin ite and second order cone optimization problems,” Ec onPap ers 73 , August 2002, Tilbu r g Un iv ersity , Center for Economic Researc h. [19] P . Gup ta and P . R. Kumar, “The capacit y of w ir eless netw orks,” IEEE T r ans. on Information The ory , v ol. 46, no. 2, pp. 388–4 04, Marc h 2000. [20] A. Dimakis, A. D. S arwate, and M. J. W ain w righ t, “Geographic gossip: Efficien t a veragi n g for sensor net works,” IEEE T r ans. on Signal Pr o c essing , to app ear. [21] W. L i and H. Dai, “Lo cation-ai d ed fast distribu ted consensus,” IEEE T r ansactions on Infor- mation The ory , June 2007, su bmitted. [22] ——, “Cluster-based fast distrib uted consensus,” IEEE Int. Conf. on A c oustics, Sp e e ch and Signal Pr o c essing (ICASSP) , April 2007. [23] K. Jung and D. Shah, “F ast gossip via nonrev ers ib le random wal k,” IE E E Information The ory Workshop (ITW) , Marc h 2006. [24] F. Ch en , L. Lov´ asz, and I. P ak, “Lifting mark ov chains to sp eed u p mixing,” STOC ’ 99: Pr o c e e dings of the thirty-first annual ACM symp osium on The ory of c omputing , pp . 275–281 , 1999. [25] Y. Hatano and M. Mesbahi, “Agreement ov er random net wo r ks,” IEE E Conf. on De cision and Contr ol , vo l. 2, pp. 2010–2015 , Decem b er 2004. [26] M. Meh ya r , D. Sp anos, J. P ongsa japan , S . H. Lo w, and R. M. Murray , “Asyn c hronous dis- tributed av eraging on comm unication net works,” IEEE/ACM T r ansactio ns on N etworking , v ol. 15, no. 3, pp . 512–520, June 2007. [27] C. C . Moallemi and B. V. Roy , “Consensu s propagation,” IEEE T r ansactions on Information The ory , v ol. 52, no. 11, p p. 4753–4 766, No v ember 2006. 18