Graph Signal Processing: Filter Design and Spectral Statistics

Graph Signal Processing: Filter Design and Spectral Statistics Stephen Kruzick and Jos ´ e M. F . Moura Carnegie Mellon University , Department of Electrical Engineering 5000 Forbes A venue, Pittsbur gh, Pennsylv ania 15213 Abstract —Graph signal pr ocessing analyzes signals supported on the nodes of a graph by defining the shift operator in terms of a matrix, such as the graph adjacency matrix or Laplacian matrix, related to the structur e of the graph. With respect to the graph shift operator , polynomial functions of the shift matrix perform filtering. An application considered in this paper , con vergence acceleration filters for distributed av erage consensus may be viewed as lowpass graph filters periodically applied to the states. Design of graph filters depends on the shift matrix eigendecomposition. Consequently , random graphs pr esent a challenge as this inf ormation is often difficult to obtain. Nevertheless, the asymptotic behavior of the shift matrix empirical spectral distrib ution pro vides a substitute for suitable random matrix models. This paper employs deterministic approximations for empirical spectral statistics from other works to propose optimization criteria for consensus acceleration filters, evaluating the r esults through simulation. Index T erms —graph signal processing, distributed average consensus, filter design, Chebyshev appr oximation, random graphs, random matrices, spectral statistics I . I N T RO D U C T I O N Signal processing applications increasingly benefit from analysis methods that account for underlying structure in data, often modeled by a graph [1, 2]. Graph signal processing ana- lyzes data in terms of this structure by representing signals as functions on the graph nodes and by defining the shift operator as a matrix that respects the graph structure, such as the adjacency matrix or the Laplacian matrix [1, 2]. Under these definitions, left multiplication of polynomials in the graph shift matrix with the graph signal vector performs filtering [1]. When the graph shift matrix is diagonalizable, decomposition of the signal vector in the basis of shift matrix eigenv ectors provides an analogy to the Fourier transform, where the eigen vectors serve as pure frequency signal components [3]. Because the shift matrix eigen values relate to the total v ariation of the eigen vectors, the y provide a notion of frequenc y that can be used to order the eigen vectors [3]. F or a shift matrix W G , the response of a graph filter p ( W G ) to an eigenv ector v of W G with corresponding eigen value λ is p ( W G ) v = p ( λ ) v [4]. Thus, the eigen values of the shift matrix play a critical role in graph filter design. Random graphs and matrices result Stephen Kruzick (skruzick@andrew .cmu.edu) and Dr. Jos ´ e M. F . Moura (moura@ece.cmu.edu) are members of the Department of Electrical and Computer Engineering at Carnegie Mellon University in Pittsburgh, P A, USA. This work was supported by NSF grant #CCF1513936. in random eigenv alues, which complicates the filter design process. This paper examines design of graph filters used to improv e the con vergence rate of the distrib uted av erage consensus algo- rithm for large scale random netw orks. In distrib uted a verage consensus, the network nodes must reach agreement on the mean of data distributed among the nodes through an iterativ e algorithm only using local communications [5], a problem relev ant to applications such as sensor data fusion [6], flocking of multiagent systems [7], and processor load balancing [8]. In each iteration of the algorithm, the network nodes update a state variable, initially set to the node data, by computing a linear combination of neighbor node states. Thus, distributed av erage consensus on a graph G is described by the dynamic system x n +1 = W G x n (1) where x n collects the node states at time n , x 0 collects the initial node data, and the weight matrix W G describes the state update, which respects the graph structure G [9]. If the weight matrix W G satisfies W G 1 = 1 , ` > W G = ` > , ρ ( W G − J ` ) < 1 (2) where ρ is the spectral radius and J ` = 1 ` > / ` > 1 is the weighted av erage consensus matrix, then x n will asymptot- ically approach a weighted average of the initial data J ` x 0 = ( ` > x 0 / ` > 1 ) 1 at rate closely related to ln ρ ( W G − J ` ) [9]. Note that if W G is a doubly stochastic matrix, ` = 1 so the unweighted average is produced. Con vergence rate acceleration for distributed av erage con- sensus may be accomplished by periodically applying a filter to previous node states. F or a filter of de gree d , the modified algorithm performs the state update in (1) and additionally sets x n := d X k =0 a k x n − d + k , n ≡ 0 ( mod d ) (3) on ev ery d th iteration, where the filter coefficients form the polynomial p ( W G ) = P d k =0 a k W k G [10]. If the filter satisfies p ( W G ) 1 = 1 , ` > p ( W G ) = ` > , ρ ( p ( W G ) − J ` ) < 1 (4) the state conv erges to the weighted average consensus at rate closely related to 1 /d ln ρ ( p ( W G ) − J ` ) [10]. Thus, consensus acceleration filters should be designed to reduce ρ ( p ( W G ) − J ` ) . For known network topologies and weight matrices with K distinct eigenv alues, finite time consensus filters use polynomials with zeros at the K − 1 distinct eigen values λ 6 = 1 to reach consensus in a finite number of iterations, representing an extreme e xample potentially requiring high filter degree [11]. For filters of lower fixed degree 1 ≤ d ≤ K − 1 and kno wn weight matrices, [10] formulates a semidefinite program that yields the optimal solution. For random switching network topologies, [10] also attempts to design consensus acceleration filters by applying the semidefinite program to the mean weight matrix. Howe ver , this can lead to suboptimal results or e ven di ver gence when the eigen values of the mean matrix do not suf ficiently approximate the random weight matrix eigen values. Therefore, consensus acceleration filter design should employ , when possible, a more complete understanding of the weight matrix spectral statistics. Because the consensus weight matrix respects the graph structure and distributed av erage consensus asymptotically produces a constant vector , consensus acceleration filters can be vie wed as lo wpass graph filters. The methods presented in this paper combine deterministic approximations of the empirical eigen value distribution of large scale random ma- trices with linear programming for Chebyshe v approximation to optimize the conv ergence rate for large scale constant (with respect to time iterations) random networks. Related literature features contrasting approaches including filters based on Chebyshev polynomials of increasing degree [12], ARMA filter response specifications selected independently from the graph [13], spectral clustering based on the smallest and sec- ond largest eigenv alue modulus [14], and different asymptotic methods [15]. In this paper , section II describes background information from random matrix theory used to describe the weight matrix spectral statistics. Section III proposes the filter optimization problem and discusses numerical practicalities. Section IV supports the proposed design method with simula- tion results and also discusses choice of weight matrix. Finally , Section V provides concluding remarks. For an extended version of this paper , refer to [16]. I I . B AC K G RO U N D : R A N D O M M AT R I X T H E O RY Filter design for signal processing on graphs depends on knowledge of the graph shift matrix eigen values. As previously noted, the value of the filter polynomial at each eigen value determines the response to the corresponding eigen vector . Therefore, spectral information should inform design criteria for the filter response. Howe ver , for contexts in which the graph is subject to stochastic influences, the corresponding shift matrix and associated eigenv alues also become random variables. W ith few e xceptions [17], the joint eigenv alue distri- bution typically proves analytically elusi ve. Ne vertheless, for some matrices of large size [18]–[22], useful information can be obtained through the asymptotic behavior of the empirical eigen value distribution defined below . For an N × N matrix W N with real eigenv alues λ 1 , . . . , λ N , the empirical spectral distribution function, defined by F W N ( λ ) = 1 N N X i =1 χ ( λ i ≤ λ ) (5) Fig. 1: Example empirical spectral distribution (blue shaded) and deterministic approximation (black curve) for row normal- ized Laplacian of 2 D lattice stochastic block model network from Figure 5 described in Section IV where χ is the indicator function, counts the number of eigen- values on the interval ( −∞ , λ ] [21]. Likewise, the empirical spectral density function, defined by f W N ( λ ) = 1 N N X i =1 δ ( λ − λ i ) (6) where δ is the Dirac delta function, indicates eigen value locations [21]. While these are function v alued random vari- ables due to the random eigen values, the empirical spectral distribution and density may sometimes be approximated by deterministic functions for random matrices of large size. For a family W N of random matrices of dimensions N × N parameterized by N , the sequence F W N of empirical spectral distributions may approach a limiting spectral distribution F lim . W ell kno wn examples include the W igner semicircular law [18], the Marchenko-P astur Law [21], and the Girko circu- lar law [19]. Similarly , a sequence of deterministic distribution functions F N that asymptotically approximate the empirical spectral distribution sequence is known as a deterministic equiv alent for the sequence [21]. The stochastic canonical equation methods of Girko provide one approach to obtaining these deterministic equiv alents [19]. For random matrix models that satisfy certain re gularity con- ditions, the Stieltjes transform of a deterministic equiv alent distribution can be computed by solving a system of equations that depends on the random matrix model parameters [19]. These methods allow analysis of matrices with elements that are independent but not necessarily identically distributed, as would arise from many random graph models, such as percolations of non-complete supergraphs [19]. Use of Girko’ s methods to approximate empirical spectral distributions of graph adjacency matrices, normalized adjacency matrices, and normalized Laplacian matrices has been e xamined in [23]– [25]. This paper employs the deterministic equiv alents com- puted using these methods along with the filter design criteria proposed in Section III to produce the simulation results appearing in Section IV. I I I . F I LT E R D E S I G N M E T H O D In order to design graph filters that accelerate the distributed av erage consensus process, note that the worst case consensus con vergence error at time n = md is k p ( W ) m − J ` k 2 where d is the degree of filter polynomial p ( W ) = P d k =0 a k W k and ` > W = ` > . When the weight matrix can be diagonalized by matrix V , it follows that ρ ( p ( W ) − J ` ) m ≤ k p ( W ) m − J ` k 2 ≤ k V k 2   V − 1   2 ρ ( p ( W ) − J ` ) m . (7) Thus, bounds for the worst case con ver gence rate can be optimized by designing the filter to minimize ρ ( p ( W ) − J ` ) . When comparing performance of dif ferent length filters, the per iteration con ver gence rate 1 /d ln ρ ( p ( W ) − J ` ) should be used. Given a random N × N weight matrix model with deterministic approximation f N for the empirical spectral density , the proposed optimization (8) solves this problem in the space of polynomials with degree at most d . The condition that p ( W ) 1 = 1 and the condition that W 1 = 1 impose the equality constraint p (1) = 1 . min p ∈ P d max λ ∈ Λ κ,τ | p ( λ ) | s.t. p (1) = 1 Λ κ,τ = { λ < 1 − κ | f N ( λ ) > τ } (8) The set Λ κ,τ for small constants κ and τ defines the spectral region of interest, where κ guarantees a transition region around the equality constraint and τ specifies where the density function f N has negligible value. Intuitively , (8) minimizes the worst case graph filter frequency response at eigen values of W captured as the support of the deterministic approximation f N to the empirical eigenv alue distribution. While a loss of robustness is possible if an eigen value falls outside the approximate spectral distribution support, the true support is well approximated asymptotically large networks with spectral conv ergence behavior and reasonable filter de- grees, and additional constraints could be added to e xplicitly prev ent this. This minimax polynomial optimization can be understood in the context of Chebyshe v approximation and produces equirip- ple behavior . While a Remez algorithm could be employed, it can also be formulated as a linear program for simplicity . In practice, the substitution in (9)-(10) may be utilized to improve the conditioning of the resulting linear program by eliminating the equality constraint, where the { φ n } denote a basis of scaled and shifted Chebyshev polynomials of the first kind chosen as in (13). p ( λ ) = 1 + (1 − λ ) q ( λ ) (9) q ( λ ) = d − 1 X n =0 a n φ n ( λ ) (10) This results in the modified optimization problem in (11), which can be used to obtain p ( λ ) by finding q ( λ ) . min q ∈ P d − 1 max λ ∈ Λ κ,τ | (1 − λ ) ( q ( λ ) + 1 / (1 − λ )) | Λ κ,τ = { λ < 1 − κ | f N ( λ ) > τ } (11) The linear program in (12) solves the problem by minimizing the bound on the objective function absolute value at rep- resentativ e sample points Λ S ⊂ Λ κ,τ (e.g., sev eral hundred uniformly spaced points). min { a n } ,> 0  s.t. (1 − λ i ) d − 1 X n =0 a n φ n ( λ i ) − 1 1 − λ i ! <  − (1 − λ i ) d − 1 X n =0 a n φ n ( λ i ) − 1 1 − λ i ! <  for all λ i ∈ Λ S (12) For good numerical performance, the basis of Chebyshev polynomials { φ n } should be scaled to the spectral region of interest Λ κ,τ as follows, where T n is the degree n standard Chebyshev polynomial of the first kind [26]. φ n ( λ ) = T n  λ − β α  α = λ max − λ min 2 , β = λ max + λ min 2 λ max = max (Λ κ,τ ) , λ min = min (Λ κ,τ ) (13) I V . S I M U L A T I O N S In order to demonstrate the performance of the filters designed using the proposed methods, this section provides supporting simulation results for constant random networks. The first pair of simulations compares the relativ e conv ergence rates of filters for weight matrices based on the unnormalized Laplacian matrix and weight matrices based on the row- normalized Laplacian matrix. The second pair of simulations compares the proposed method to other filter design methods av ailable in the literature. The simulations presented in this section cover two random network models. An Erd ˝ os-R ´ enyi network on N nodes de- scribes a random graph model in which each pair of nodes connects according to independent Bernoulli trials with link probability θ [27]. A D -dimensional lattice stochastic block model consists of N 1 × · · · × N D populations, each with M nodes. The populations correspond to D -tuples and collec- tiv ely form a lattice in the sense of [28]. Nodes connect to other nodes according to independent Bernoulli trials if their population tuples dif fer by at most one symbol, with link probability θ 0 within populations and θ k between nodes in populations along lattice dimension k . The adjacency matrix empirical spectral distributions of both of these models are amenable to analysis by the methods of Girk o, as done in [23]– [25]. The second simulation uses the spectral distrib ution approximation results from [23]–[25] without repeating the deriv ation. The proposed optimization problem does not specify a particular scheme for choosing the weight matrix W G for graph G , requiring only that W G satisfy the weighted consensus conditions (2). A common choice for the weight matrix, W G = I − αL ( G ) depends on the unnormalized graph Laplacian matrix L ( G ) = D ( G ) − A ( G ) where D ( G ) is the diagonal degree matrix and A ( G ) is the graph adjacency Fig. 2: Consensus conv ergence rates for an Erd ˝ os-R ´ enyi network with 2000 nodes and connection probability θ = 0 . 03 using filters based on weights W = I − αL and W = I − α b L R Fig. 3: Consensus con vergence rates for a 2 D lattice stochastic block model network with 3 × 7 populations of 100 nodes and connection probabilities ( θ 0 , θ 1 , θ 2 ) = (0 . 15 , 0 . 09 , 0 . 06) using filters based on weights W = I − α L and W = I − α b L R Fig. 4: Consensus conv ergence rates for an Erd ˝ os-R ´ enyi network with 1000 nodes and connection probability θ = 0 . 05 using filters of degree d = 1 , . . . , 10 based on b L R for listed methods Fig. 5: Consensus con vergence rates for a 2 D lattice stochastic block model network with 3 × 4 populations of 100 nodes and connection probabilities ( θ 0 , θ 1 , θ 2 ) = (0 . 10 , 0 . 10 , 0 . 10) using filters of degree d = 1 , . . . , 10 based on b L R for listed methods matrix. The scale parameter α must be chosen to satisfy the spectral radius constraint. Note that W G = I − αL ( G ) is a doubly stochastic matrix, so the unweighted average is produced. Howe ver , the spectral statistics of the Laplacian are typically not approachable using Girko’ s methods. In the case of Erd ˝ os-R ´ enyi networks, the limit of the Laplacian empirical spectral distribution has density given by the free con volution [22] of a Gaussian distribution and a semicircular distribution [29], but results for other models are typically inaccessible. Alternativ ely , the weight matrix scheme W G = I − α b L R ( G ) where b L R ( G ) = I − D ( G ) − 1 A ( G ) is the row-normalized Laplacian can be selected. This is a row-stochastic matrix with left eigen vector ` = d , the vector of node degrees, correspond- ing to λ = 1 . Although this leads to a weighted a verage, the unweighted a verage can be produced through premultiplication by the corrective transform  d > 1  /  1 > 1  D − 1 , which can be applied at each node using the av erage node degree and the local node degree. Because the row-normalized Laplacian is more amenable to analysis using the methods of Girk o [23]– [25], W G = I − α b L R ( G ) presents an appealing choice for use with (8). Figures 2 and 3 show the relativ e conv ergence rates for each of these weight matrices using the simulated expected empirical densities for an Erd ˝ os-R ´ enyi network and a 2 D lattice stochastic block model with 1 /α chosen as the approximate center of the distribution support. Note that W G = I − α b L R ( G ) outperforms W G = I − αL ( G ) further recommending use of the weight matrix based on the row- normalized Laplacian. Intuitively , this observation is expected due to the conv olution-like composition of the Laplacian empirical spectral density resulting in less compact support. The second group of simulations compares the results obtained from the proposed design method (8) to those from the mean matrix semidefinite program method described in [10] and from the Newton interpolating polynomial method described in [10] for various filter degrees. A deterministic approximation for the empirical spectral density was com- puted using Girko’ s stochastic canonical equation theorem as in [23]–[25] to be used with the proposed optimization method (8). The results for the proposed optimization method compare fa vorably to these methods as seen in Figures 4 and 5, which plot the per iteration con vergence rates for filters of de- gree d = 1 , . . . , 10 . Note that because the polynomial based on the semidefinite program only has a unique solution for filter degree less than the number of distinct mean weight matrix eigen values, it only appears for d = 1 in Figure 4 and d ≤ 4 in Figure 5. Furthermore, the results were compared with filters designed using the proposed optimization method (8) with an oracle for the true empirical spectral distribution. The proposed method using the deterministic equi valent distrib ution achiev es nearly equal results, demonstrating good performance. V . C O N C L U S I O N Filter design for signal processing on random graphs re- quires the ability to obtain information about the random shift matrix eigen values. Thus, methods from random matrix theory that capture asymptotic deterministic structure in the empirical spectral distributions of suitable matrices provide useful tools for filter design in the context of large scale ran- dom graphs. This paper proposed an optimization problem to deriv e con vergence acceleration filters for distributed a verage consensus, which can be understood as lowpass graph filters. In practical terms, these filters enable improv ed accuracy over a fixed number of iterations or a giv en lev el of accuracy in fewer iterations. The proposed method combines Chebyshev approximation techniques with deterministic equiv alents for shift matrices derived in other papers using Girko’ s stochastic canonical equations. Simulation results demonstrate that the filters derived perform well on constant random networks, comparing fav orably to other tested methods. Consideration of weight matrices based on the unnormalized and ro w- normalized Laplacians suggest faster conv ergence can be achiev ed using random row-normalized Laplacians. Further work will focus on analysis of time-varying random networks and additional graph signal processing applications. R E F E R E N C E S [1] A. Sandryhaila and J. M. F . Moura, “Discrete signal processing on graphs, ” IEEE T ransactions on Signal Pr ocessing , vol. 61, no. 7, pp. 1644–1656, April 2013. [2] D. Shuman, S. Narag, P . Frossard, A. Ortega, and P . V ander gheynst, “The emerging field of signal processing on graphs, ” IEEE Signal Processing Magazine , vol. 30, no. 3, pp. 83–98, May 2013. [3] A. Sandryhaila and J. M. F . Moura, “Discrete signal processing on graphs: Frequency analysis, ” IEEE T ransactions on Signal Processing , vol. 62, no. 12, pp. 3042–3054, June 2014. [4] A. Sandryhaila and J. M. F . Moura, “Big data analysis with signal processing on graphs: Representation and processing of massive data sets with irregular structure, ” IEEE Signal Processing Magazine , vol. 31, no. 5, pp. 80–90, Sept. 2014. [5] R. Olfati-Saber, A. Fax, and R. M. Murray, “Consensus and cooperation in networked multi-agent systems, ” Pr oceedings of the IEEE , vol. 98, no. 7, pp. 1354–1355, June 2010. [6] L. Xiao and S. Boyd, “Fast linear iterations for distributed averaging, ” Systems and Control Letter s , vol. 53, pp. 65–78, Feb. 2004. [7] R. Olfati-Saber, “Flocking for multi-agent dynamic systems: Algorithms and theory , ” IEEE T ransactions on Automatic Control , vol. 51, no. 3, pp. 401–420, Mar . 2006. [8] G. Cybenko, “Dynamic load balancing for distributed memory multi- processors, ” Journal of P arallel and Distrib uted Computing , v ol. 7, no. 2, pp. 279–301, Oct. 1989. [9] S. Kar and J. M. F . Moura, “Consensus+innovations distributed inference over networks: Cooperation and sensing in networked systems, ” IEEE Signal Pr ocessing Magazine , vol. 30, no. 3, pp. 99–109, May 2013. [10] E. Kokiopoulou and P . Frossard, “Polynomial filtering for fast con- ver gence in distributed consensus, ” IEEE T ransactions on Signal Pr ocessing , vol. 57, pp. 342–354, Jan. 2009. [11] A. Sandryhaila, S. Kar, and J. M. F . Moura, “Finite-time distributed consensus through graph filters, ” Pr oceedings of ICAASP 2014 , pp. 1080–1084, May 2014. [12] E. Montijano, J. I. Montijano, and C. Sagues, “Chebyshev polynomials in distributed consensus applications, ” IEEE T ransactions on Signal Pr ocessing , vol. 61, no. 3, pp. 693–706, Feb. 2013. [13] A. Loukas, A. Simonetto, and G. Leus, “Distributed autoregressiv e moving a verage graph filters, ” IEEE Signal Processing Letter s , vol. 22, no. 11, pp. 1931–1935, Nov . 2015. [14] S. Apers and A. Sarlette, “ Accelerating consensus by spectral clustering and polynomial filters, ” IEEE T ransactions on Contr ol of Network Systems , vol. 4, no. 3, pp. 544–554, Sept. 2016. [15] F . Gama and A. Ribeiro, “W eak law of large numbers for stationary graph processes, ” 2017 IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP 2017) , pp. 4124–4128, March 2017. [16] S. Kruzick and J. M. F . Moura, “Optimal filter design for signal processing on random graphs: Accelerated consensus, ” Submitted , June 2016. [17] C. A. Trace y and H. Widom, “The distrib ution of the largest eigen- value in the Gaussian ensembles, ” Caloger o-Moser-Sutherland Models, CRM Series in Mathematical Ph ysics, vol. 4, pp. 461–472, Springer Science+Business Media, 2000. [18] E. W igner, “On the distribution of the roots of certain symmetric matric- es, ” The Annals of Mathematics , vol. 67, no. 2, pp. 325–327, Mar. 1958. [19] V . Girko, Theory of Stoc hastic Canonical Equations , vol. 1, pp. 1-3, 365-366, Springer Science+Business Media, 2001. [20] Z. Bai, “Methodologies in spectral analysis of large random matrices, a revie w , ” Statistica Sinica , vol. 9, no. 3, pp. 611–677, July 1999. [21] R. Couillet and M. Debbah, Random Matrix Methods for Wir eless Communications , Cambridge Uni versity Press, 2011. [22] A. M. Tulino and S. V erd ´ u, F oundations and T r ends in Communications and Information Theory: Random Matrix Theory and W ir eless Commu- nications , vol. 1, 2004. [23] K. A vrachenkov , L. Cottatellucci, and A. Kadavankandy , “Spectral prop- erties of random matrices for stochastic block model, ” 4th International W orkshop on Physics-Inspir ed P aradigms in W ireless Communications and Networks , pp. 537–544, May 2015. [24] S. Kruzick and J. M. F . Moura, “Spectral statistics of lattice graph percolations, ” Arxiv: https://arxiv .org/abs/1611.02655 , Sept. 2016. [25] S. Kruzick and J. M. F . Moura, “Spectral statistics of lattice graph structured, non-uniform percolations, ” 42nd IEEE International Confer - ence on Acoustics, Speech, and Signal Pr ocessing (ICASSP 2017) , pp. 5930–5934, Sept. 2016. [26] J. P . Boyd, Chebyshev and F ourier Spectral Methods , Dov er , second edition, 2001. [27] B. Bollob ´ as, Random Graphs , Cambridge University Press, 2001. [28] R. Laskar, “Eigen values of the adjacency matrix of cubic lattice graphs, ” P acific Journal of Mathematics , v ol. 29, no. 3, pp. 623–629, July 1969. [29] X. Ding and T . Jiang, “Spectral distributions of adjacency and Laplacian matrices of random graphs, ” The Annals of Applied Probability , vol. 20, no. 6, pp. 2086–2117, 2010.