Spectral Statistics of Lattice Graph Percolation Models
11 Spectral Statistics of Lattice GraphPercolation Models
Stephen Kruzick and Jos´e M. F. Moura,
Fellow, IEEE
Abstract —In graph signal processing, the graph adja-cency matrix or the graph Laplacian commonly definethe shift operator. The spectral decomposition of the shiftoperator plays an important role in that the eigenvaluesrepresent frequencies and the eigenvectors provide a spec-tral basis. This is useful, for example, in the design of filters.However, the graph or network may be uncertain dueto stochastic influences in construction and maintenance,and, under such conditions, the eigenvalues of the shiftmatrix become random variables. This paper examinesthe spectral distribution of the eigenvalues of randomnetworks formed by including each link of a D -dimensionallattice supergraph independently with identical probability,a percolation model. Using the stochastic canonical equa-tion methods developed by Girko for symmetric matriceswith independent upper triangular entries, a deterministicdistribution is found that asymptotically approximates theempirical spectral distribution of the scaled adjacencymatrix for a model with arbitrary parameters. The mainresults characterize the form of the solution to an impor-tant system of equations that leads to this deterministicdistribution function and significantly reduce the numberof equations that must be solved to find the solution fora given set of model parameters. Simulations comparingthe expected empirical spectral distributions and the com-puted deterministic distributions are provided for sampleparameters. Index Terms —graph signal processing, random net-works, random graph, random matrix, eigenvalues, spectralstatistics, stochastic canonical equations
I. I
NTRODUCTION
Much of the increasingly vast amount of data availablein the modern day exhibits nontrivial underlying struc-ture that does not fit within classical notions of signalprocessing. The theory of graph signal processing hasbeen proposed for treating data with relationships andinteractions best described by complex networks. Withinthis framework, signals manifest as functions on thenodes of a network. The shift operator used to analyze
The authors are with the Department of Electrical and Com-puter Engineering, Carnegie Mellon University, Pittsburgh, PA15213 USA (ph: (412)2686341; e-mail: [email protected];[email protected]).This work was supported by NSF grant these signals is provided by a matrix related to the net-work structure, such as the adjacency matrix, Laplacianmatrix, or normalized versions thereof [1][2]. Decompo-sition of a signal according to a basis of eigenvectorsof the shift operator serves a role similar to that of theFourier Transform in classical signal processing [3]. Inthis context, multiplication by polynomial functions ofthe chosen shift operator matrix performs shift invariantfiltering [1].The eigenvalues of the shift operator matrix playthe role of graph frequencies and are important in thedesign and analysis of graph signals and graph filters. If W is a diagonalizable graph shift operator matrix witheigenvalue λ for eigenvector v such that W v = λ v and,for example, if a filter is implemented on the networkas P ( W ) where P is a polynomial, then P ( W ) hascorresponding eigenvalue P ( λ ) for v by simultaneousdiagonalizability of powers of W [4]. The framework ofgraph signal processing regards P ( λ ) as the frequencyresponse of the filter [5]. Hence, knowledge of the eigen-values of W informs the design of the filter P when theeigenvalues of P ( W ) should satisfy desired properties.Furthermore, the eigenvalues relate to a notion of signalcomplexity known as the signal total variation, whichhas several slightly different definitions depending oncontext [2][3][5][6]. For purposes of motivation, takingthe shift operator to be the row-normalized adjacencymatrix (cid:98) A , define the l p total variation of the signal x as TV G ( x ) = (cid:13)(cid:13)(cid:13)(cid:16) I − (cid:98) A (cid:17) x (cid:13)(cid:13)(cid:13) pp = (cid:13)(cid:13)(cid:13)(cid:98) L x (cid:13)(cid:13)(cid:13) pp (1)which sums over all network nodes the p th power ofthe absolute difference between the value of the signal x at each node and the average of the value of x atneighboring nodes [5][6]. Thus, if v is a normalizedeigenvector of the row-normalized Laplacian (cid:98) L witheigenvalue λ , v has total variation | λ | p . The eigenvectorsthat have higher total variation can be viewed as morecomplex signal components in much the same waythat classical signal processing views higher frequencycomplex exponentials.As an application, consider a connected, undirected a r X i v : . [ c s . NA ] S e p etwork on N nodes with normalized Laplacian (cid:98) L withthe goal of achieving distributed average consensus viaa graph filter. For such a network, there is a simpleLaplacian eigenvalue λ ( (cid:98) L ) = 0 corresponding to theaveraging eigenvector v = and other eigenvalues < λ i ( (cid:98) L ) ≤ for ≤ i ≤ N . Any filter P suchthat P (0) = 1 and | P ( λ i ( (cid:98) L )) | < for ≤ i ≤ N will asymptotically transform an initial signal x withaverage µ to the average consensus signal ¯x = µ uponiterative application [7]. If the eigenvalues λ i ( (cid:98) L ) for ≤ i ≤ N are known, consensus can be achieved infinite time by selecting P to be the unique polynomialof degree N − with P (0) = 1 and P ( λ i ( (cid:98) L )) = 0 [4].Note that the averaging eigenvector has total variation by the above definition (1) and that all other, morecomplex eigenvectors are completely removed by thefilter. Thus, a finite time consensus filter represents themost extreme version of a non-trivial lowpass filter. Withpolynomial filters of smaller fixed degree d , knowledgeof the eigenvalues can be used to design filters of a givenlength for optimal consensus convergence rate, which isgiven by ρ (cid:16) P (cid:16)(cid:98) L (cid:17)(cid:17) /d = max ≤ i ≤ N (cid:12)(cid:12)(cid:12) P (cid:16) λ i (cid:16)(cid:98) L (cid:17)(cid:17)(cid:12)(cid:12)(cid:12) /d (2)as studied in [7]. This can also be attempted in situationsin which the graph is a random variable, leading touncertainty in the eigenvalues. For example, [7] also pro-poses two methods to accomplish this for certain randomswitching network models. One relies on the eigenvaluesof the expected iteration matrix. The other attempts tofilter over a wide range of possible eigenvalues withouttaking the model into account. Neither takes into accountany information about how the eigenvalue distributionspreads. This type of information could be relevant tograph filter design, especially when the network doesnot switch too often relative to the filter length, renderingstatic random analysis applicable.The empirical spectral distribution F W ( x ) of a Hermi-tian matrix W , which counts the fraction of eigenvaluesof W in the interval ( −∞ , x ] , describes the set of eigen-values [8]. For matrices related to random graphs, theempirical spectral distributions are, themselves, function-valued random variables. Often, it is not possible toknow the joint distribution of the eigenvalues, but thebehavior of the empirical spectral distribution can bedescribed asymptotically. For instance, the empiricalspectral distributions of Wigner matrices converge to theintegral of a semicircle asymptotically in the matrix sizeas described by Wigner’s semicircle law [9]. In othercases, a limiting distribution may or may not exist, but asequence of deterministic functions can still approximate Fig. 1: The illustration shows an example lattice graphwith three dimensions of size × × with some nodetuples labeled. Each group of circled nodes, which differby exactly one symbol, represents a complete subgraph.the empirical spectral distribution [10]. The stochasticcanonical equation techniques of Girko detailed in [8],which allow analysis of matrices with independent butnon-identically distributed elements as necessary whenstudying the adjacency matrices of many random graphmodels, provide a method to obtain such deterministicequivalents.This paper analyzes the empirical eigenvalue distribu-tions of the adjacency matrices of random graphs formedby independent random inclusion or exclusion of eachlink in a certain supergraph, a bond percolation model[11]. Specifically, the paper examines adjacency matricesof percolation models of D -dimensional lattice graphs,in which each D -tuple with M d possible symbols inthe d th dimension has an associated graph node and inwhich two nodes are connected by a link if the cor-responding D -tuples differ by only one symbol. Thesegraphs generalize other commonly encountered graphs,such as complete graphs and the cubic lattice graph[12], to an arbitrary lattice dimension. An illustrationfor the three dimensional case appears in Figure 1, inwhich a complete graph connects each group of circlednodes. The stochastic canonical equation techniques ofGirko provide the key mathematical tool employed hereto achieve our results. The resulting deterministic ap-proximations to the adjacency matrix empirical spectraldistribution can, furthermore, provide information aboutthe empirical spectral distributions of the normalizedadjacency matrices and, thus, the normalized Laplacians.This work is similar to that of [13], which analyzesthe adjacency matrices of a different random graphmodel known as the stochastic block model using similartools and suggests applications related to the study ofepidemics and algorithms based on random walks.Section II discusses relevant background informationon random graphs and the spectral statistics of randommatrices. It also introduces an important mathematical2esult from [8], presented in Therorem 1, which can beused to find deterministic functions that approximate theempirical spectral distribution of certain large randommatrices by relating it to the solution of a system ofequations (16). Section III introduces lattice graphs andaddresses application of Theorem 1 to scaled adjacencymatrices of percolation models based on arbitrary D -dimensional lattice graphs. The main contributions ofthis paper appear in Theorem 2, which derives the formof the solution to the system of equations (16) forarbitrary parameters, and also in Corollary 1, whichdescribes the particular solution for given model pa-rameters. The relationship of the empirical spectraldistribution of the adjacency matrix to those of thesymmetrically normalized adjacency matrix and the row-normalized adjacency matrix is also noted. For someexample lattice graphs, these results are used to providedeterministically calculated distribution functions thatare compared against the expected empirical spectraldistributions. Finally, Section IV provides concludinganalysis. II. B ACKGROUND
A. Random Graphs An undirected graph (network) G consists of anordered pair ( V , E ) in which V denotes a set of nodes(vertices) and E denotes a set of undirected links(edges) . These undirected links are unordered pairs { v i , v j } = { v j , v i } of nodes v i , v j ∈ V with v i (cid:54) = v j .Note that self-loops , links formed as { v, v } , are excludedfrom this definition. The graph adjacency matrix A ( G ) encapsulates the graph information with A ( G ) ij = 1 if { v j , v i } ∈ E and A ( G ) ij = 0 if { v j , v i } / ∈ E .A subgraph G sub of an undirected graph is a graph ( V sub , E sub ) with V sub ⊆ V and E sub ⊆ E , and G issaid to be a supergraph of G sub . Random undirectedgraphs , which are undirected graph-valued random vari-ables, model uncertainty in the network topology. Thereare many commonly studied random graph models,which may be specified by a probability distribution onthe collection of possible link sets P ( V × V ) between agiven number of nodes |V| = N .One of the most frequently considered random graphmodels, the Erd¨os-R´enyi model , considers a randomgraph G ER ( N, p ( N )) on N nodes formed by includingeach link between two different nodes according toindependent Bernoulli trials that succeed with proba-bility p ( N ) , which is often allowed to vary with N [14]. This distribution has important properties regardingasymptotic connectedness behavior as N increases thatcondition on the relationship between N and p ( N ) [15]. Additionally, the relationship between Erd¨os-R´enyirandom graphs and Wigner matrices leads to asymptoticspectral statistics that may be characterized [16].Similarly, starting from an arbitrary supergraph G sup ,N with N nodes, the graph-valued random variable G perc ( G sup ,N , p ( N )) that results from including eachlink of G sup ,N according to independent Bernoulli trialsthat succeed with probability p ( N ) is known as a bondpercolation model . Thus, an Erd¨os-R´enyi model is abond percolation model with a complete supergraph.Typically, study of bond percolation models concernstheir asymptotic connectivity behavior, although that isnot the purpose of this paper. Here, the terminology issimply used to refer to the type of random graph modelsconsidered in this paper. B. Spectral Statistics
The eigenvalues of matrices, such as the graph adja-cency matrix, that respect the network structure can beimportant in the design of distributed algorithms. Hence,an understanding of the eigenvalues of relevant randomsymmetric matrices is desired. For certain Hermitian ran-dom matrix models parameterized by matrix dimension N × N , results are available concerning the asymptoticbehavior of the eigenvalues as N increases through afunction called the empirical spectral distribution. Thedefinitions employed in this paper are provided below.Note that the eigenvalues of a Hermitian matrix mustbe real valued. Given a N × N Hermitian matrix-valuedrandom variable W N , order the N random eigenvaluessuch that λ ( W ) ≤ . . . ≤ λ i ( W N ) ≤ λ i +1 ( W N ) ≤ . . . ≤ λ N ( W N ) . The empirical spectral measure of W N given by µ w N ( X ) = 1 N N (cid:88) i =1 χ ( λ i ( W N ) ∈ X ) (3)assigns to each subset X ⊆ R the fraction of eigenvaluesthat appear in X using the indicator function χ . Relatedto this function, the empirical spectral distribution of W N given by [10] F W N ( x ) = µ W N (( −∞ , x ])= 1 N N (cid:88) i =1 χ ( λ i ( W N ) ≤ x ) (4)counts the number of eigenvalues in the interval ( −∞ , x ] , and the empirical spectral density of W N
3s given by [10] f W ( x ) = ddx F W N ( x )= 1 N N (cid:88) i =1 δ ( x − λ i ( W N )) . (5)indicates the locations of the eigenvalues using the Diracdelta function δ . The empirical spectral measure andempirical spectral distribution define each other and areoften referred to interchangeably. Note that each of thesedefinitions produces a function-valued random variable.From them can be defined the expected empiricalspectral measure µ exp ,W N = E [ µ W N ] , the expectedempirical spectral distribution F exp ,W N = E [ F W N ] ,and the expected empirical spectral density f exp ,W N =E [ f W N ] .Analysis of the empirical spectral distribution for Her-mitian matrices often involves the Stieltjes transform defined as below [10]. S F ( z ) = (cid:90) ∞−∞ x − z dF ( x ) , Im { z } (cid:54) = 0 (6)The Stieltjes transform can be inverted to obtain thecorresponding density by computing the following ex-pression [10]. F ( x ) = lim (cid:15) → + π (cid:90) x −∞ Im { S F ( λ + (cid:15)i ) } dλ (7) f ( x ) = lim (cid:15) → + π Im { S F ( x + (cid:15)i ) } (8)For the empirical spectral distribution F N of an N × N Hermitian matrix W N , the Stieltjes function computesto S F WN ( z ) = 1 N tr (cid:16) ( W N − zI N ) − (cid:17) , Im { z } (cid:54) = 0 (9)which is the normalized trace of the resolvent ( W N − zI N ) − [10].Given a sequence of Hermitian matrix-valued randomvariables { W N } ∞ N =1 , the sequence has limiting spectralmeasure µ lim ,W N provided µ W N converges weakly to µ lim ,W N . Similarly, it has limiting spectral distribution F lim ,W N provided F W N converges weakly to F lim ,W N .In some circumstances in which these limits may or maynot exist, the spectral statistics can still be describable bya sequence of deterministic functions. Given a sequenceof Hermitian matrix-valued random variables W N , a deterministic equivalent for W N with respect to thesequence of functionals g N is a sequence of deterministicmatrices W ◦ N such that lim N →∞ ( g N ( W N ) − g N ( W ◦ N )) = 0 (10) almost surely [10]. Additionally, the sequence of deter-ministic values g N ( W ◦ N ) is called a deterministic equiv-alent of the sequence of random values g N ( W N ) [10]. C. Stochastic Canonical Equations
The most critical theorem for this work provides amethod for finding deterministic equivalents of empiricalspectral distribution functions for symmetric randommatrices with independent entries, except for the rela-tion that the lower triangular entries are determined bysymmetry from the upper triangular entries. Providedsome regularity conditions (11), (12), and (13) hold,the Stieltjes transform (15) of a deterministic equivalentdistribution function can be found by solving a systemof equations (16) containing matrices [8]. Theorem 1provides the formal statement below. This result and acompilation of many other related results can be foundin [8].
Theorem 1 (Girko’s K1 Equation [8])
Considera family of symmetric matrix valued random vari-ables W N indexed by size N such that W N is an N × N symmetric matrix in which the entries onthe upper triangular region are independent. Thatis, (cid:110) ( W N ) ij | ≤ i ≤ j ≤ N (cid:111) are independent with ( W N ) ji = ( W N ) ij . Let W N have expectation B N =E [ W N ] and centralization H N = W N − E [ W N ] suchthat the following three conditions hold. Note that inorder to avoid cumbersome indexing, the index N willhenceforth be omitted from most expressions involving W N , B N , and H N . sup N max i N (cid:88) j =1 | B ij | < ∞ (11) sup N max i N (cid:88) j =1 E (cid:2) H ij (cid:3) < ∞ (12) lim N →∞ max i N (cid:88) j =1 E (cid:2) H ij χ ( | H ij | > τ ) (cid:3) = 0 for all τ > (13) Then for almost all x , lim N →∞ | F W N ( x ) − F N ( x ) | = 0 (14) almost surely, where F N is the distribution with Stieltjestransform S F N ( z ) = (cid:90) x − z dF N ( x )= 1 N N (cid:88) k =1 C kk ( z ) , Im { z } (cid:54) = 0 (15)4 nd the functions C kk ( z ) satisfy the canonical systemof equations C kk ( z ) = B − zI − ...... (cid:32) δ lj N (cid:88) s =1 C ss ( z )E (cid:2) H js (cid:3)(cid:33) l,j = Nl,j =1 − kk (16) for k = 1 , . . . , N . Note that the notation ( · ) l,j = Nl,j =1 indicates a matrix built from the parameterized contentsof the parentheses, such that X = ( X ij ) l,j = Nl,j =1 , and δ lj is the Kronecker delta function. There exists a uniquesolution C kk ( z ) for k = 1 , . . . , N to the canonicalsystem of equations (16) among the class L = { X ( z ) ∈ C | X ( z ) analytic , Im { z } Im { X ( z ) } > } .Furthermore, if inf i,j N E (cid:2) H ij (cid:3) ≥ c > , (17) then lim N →∞ sup x | F W N ( x ) − F N ( x ) | = 0 (18) almost surely, where F N is defined as above. In addition to the independence properties of theentries of the random matrix, Theorem 1 specifies threeconditions that must be verified. The first condition(11) bounds the absolute sum of the entries of eachrow, and the second condition (12) bounds the totalvariance of entries of each row. The third condition (13)is closely related to Lindberg’s condition of the centrallimit theorem. Under these conditions, a deterministicequivalent F N ( x ) for F W N ( x ) can be found by solvingthe system of equations (16) and computing the Stieltjestransform via (15), which can then be inverted. Under theadditional condition (17), the supremum converges [8].III. M AIN R ESULTS
This section examines, in particular, a random graphmodel formed by percolation of a D -dimensional latticegraph. Note that there are many types of graphs known aslattice graphs in other contexts, so the precise definitionused will be introduced here. A D -dimensional latticegraph with size M d along the d th dimension is a graph G lat = ( V , E ) in which the |V| = N = (cid:81) Dd =1 M d nodesare identified with the ordered D -tuples in Z [1 , M ] × . . . × Z [1 , M D ] and in which two nodes are connectedby a link if the corresponding D -tuples differ by exactlyone symbol [12]. In this way, if the node indices arefixed along D − of the lattice dimensions, the node-induced subgraph along the remaining dimension is a complete graph. Denoting the adjacency matrix of acomplete graph on M vertices by K M , in terms ofKronecker tensor products and adjacency matrices, theselattice graphs can be described by A ( G lat ) = D (cid:88) j =1 D (cid:79) d =1 X dj , X dj = (cid:26) K M d j = dI M d j (cid:54) = d (19)where the j th term in the summation contributes allcomplete graphs along the j th lattice dimension. Suchgraphs are of interest because they can be grown to alarge number of nodes by increasing the lattice dimen-sion sizes while the number of distinct eigenvalues ofthe adjacency matrix does not increase. The eigenvaluesof the lattice graph adjacency matrix are λ ( j , . . . , j D ) = D (cid:88) d =1 λ d ( j d ) ,λ d ( j d ) = (cid:26) M d − j d = 0 − j d = 1 (20)for j , . . . , j D = 0 , . Thus, the number of distincteigenvalues of a D -dimensional lattice graph dependsonly on the number of dimensions D and is at most D .The adjacency matrices of random graphs distributedaccording to percolation models G perc ( G lat , p ) of D -dimensional lattices are symmetric matrices with in-dependent entries, except for relations determined bysymmetry. In more precise terms, the matrix entries A ij ( G perc ) for i ≤ j are independent, and A ji ( G perc ) = A ij ( G perc ) . Thus, the scaled adjacency matrix eigenval-ues of percolations models formed from these latticescan be analyzed by the tools provided in Theorem 1.For a lattice with arbitrary number of dimensions D andsizes M d for d = 1 , . . . , D along with an arbitrary linkinclusion probability p , Theorem 2 derives the form ofthe unique solution C kk ( z ) for k = 1 , . . . , N to thesystem of equations (16). An outline of the methodsused can be found in the initial paragraph of the proof.This result will subsequently be used in Corollary 1 tocompute the Stieltjes transform of a deterministicallyequivalent distribution function for the empirical spectraldistribution. Theorem 2 (Solution Form for D -Lattice Percolation) Consider the D -dimensional lattice graph G lat with N nodes in which the d th dimension of the lattice had size d for d = 1 , . . . , D such that the adjacency matrix is A ( G lat ) = D (cid:88) j =1 D (cid:79) d =1 X dj , X dj = (cid:26) K M d j = dI M d j (cid:54) = d (21) N = D (cid:89) d =1 M d . (22) Form a random graph G perc ( G lat , p ) by independentlyincluding each link of G lat with probability p , and denotethe corresponding random scaled adjacency matrix W ,expectation B , and centralization HW ( G perc ) = 1 γ A ( G perc ) (23) B = E [ W ( G perc )] (24) H ( G perc ) = W ( G perc ) − E [ W ( G perc )] (25) where γ = γ ( M , . . . , M D ) = D (cid:88) d =1 ( M d − . (26) Let C kk ( z ) for k = 1 , . . . , N be the unique solutionto the system of equations (16) among the class L guaranteed to exist by Theorem 1, and write C ( z ) = B − zI N − . . .. . . (cid:32) δ lj N (cid:88) s =1 C ss ( z ) E (cid:2) H js (cid:3)(cid:33) l,j = Nl,j =1 − . (27) Note that C kk ( z ) is the k th diagonal entry of C ( z ) andthat uniqueness of C kk ( z ) implies uniqueness of C ( z ) .For some values of α i ,...,i D ( z ) for i , . . . , i D = 0 , C ( z ) = (cid:88) i ,...,i D =0 α i ,...,i D ( z ) D (cid:79) d =1 Y di d ,Y di d = (cid:26) K M d i d = 0 I M d i d = 1 . (28) Proof:
The theorem is proven primarily by use of symmetryin the random graph model and is organized in threeparts below. The first paragraph verifies that the threeconditions of Theorem 1 hold for the random matrixmodel W ( G perc ) , guaranteeing existence of a uniquesolution C kk ( z ) for k = 1 , . . . N to (16). Subsequently,the second paragraph demonstrates that C ( z ) is of thedesired form if and only if it is invariant to permu-tations that factor into Kronecker products of smallerpermutations along each lattice dimension. Finally, the third paragraph uses symmetries in the random graphmodel to show that C ( z ) , is, in fact, invariant underthese permutations, yielding the result in the theorem.In order to simplify the calculations used to verify theconditions of Theorem 1, a more explicit characterizationof the entries A ij is first introduced. Note that everyinteger ≤ x ≤ N can be uniquely represented in amixed-radix system as x = 1 + D (cid:88) d =1 β ( x, d ) d − (cid:89) j =1 M j (29)for integers ≤ β ( x, d ) ≤ M d − . Hence, node indices ( i, j ) can be described by their digit sequences β ( i ) = { β ( i, d ) } Dd =1 and β ( j ) = { β ( j, d ) } Dd =1 in this system.Let (cid:107) β ( i ) − β ( j ) (cid:107) count the differences between thesedigit sequences. Using this to describe the entries of thesupergraph adjacency matrix yields A ij ( G lat ) = (cid:26) (cid:107) β ( i ) − β ( j ) (cid:107) = 10 otherwise . (30)Because B ij = pγ A ij ( G lat ) , the following computationverifies the first condition (11) of Theorem 1. sup N max i N (cid:88) j =1 | B ij | = sup N max i (cid:88) (cid:107) β ( i ) − β ( j ) (cid:107) =1 | B ij | = sup N pγ D (cid:88) d =1 ( M d − p < ∞ (31)Noting that E [( H ij (cid:1) xy (cid:3) = p (1 − p ) γ if ( A ij ( G lat )) xy = 1 and that E [( H ij (cid:1) xy (cid:3) = 0 if ( A ij ( G lat )) xy = 0 , thefollowing computation verifies the second condition (12)of Theorem 1. sup N max i N (cid:88) j =1 E (cid:2) H ij (cid:3) = sup N max i (cid:88) (cid:107) β ( i ) − β ( j ) (cid:107) =1 E (cid:2) H ij (cid:3) = sup N p (1 − p ) γ (cid:32) D (cid:88) d =1 ( M d − (cid:33) ≤ p (1 − p ) < ∞ (32)Finally, the following computation demonstrates that the6hird condition (13) of Theorem 1 holds for any τ > . lim N →∞ max i N (cid:88) j =1 E (cid:2) H ij χ ( | H ij | > τ ) (cid:3) ≤ lim N →∞ max i (cid:88) (cid:107) β ( i ) − β ( j ) (cid:107) =1 E max (cid:16) p , (1 − p ) (cid:17) γ . . .. . . χ max ( p, − p ) γ > τ = 0 (33)This is true because lim N →∞ max ( p, − p ) γ = 0 , (34)so when γ becomes sufficiently large lim N →∞ χ (cid:18) max ( p, − p ) γ > τ (cid:19) = 0 . (35)Note that because H ij = 0 when (cid:107) β ( i ) − β ( j ) (cid:107) (cid:54) = 1 the additional condition (17) for the stronger conclusionof Theorem 1 does not hold.Given that C ( z ) is of the form in (28), let φ d be apermutation function on M d symbols with correspondingrow permutation matrix P d for d = 1 , . . . , D and P = P ⊗ . . . ⊗ P D with corresponding permutation function φ . It is clear that P C ( z ) P (cid:62) = P (cid:88) i ,...,i D =0 α i ,...,i D ( z ) D (cid:79) d =1 Y di d P (cid:62) = (cid:88) i ,...,i D =0 α i ,...,i D ( z ) D (cid:79) d =1 P d Y di d P (cid:62) d (36) = (cid:88) i ,...,i D =0 α i ,...,i D ( z ) D (cid:79) d =1 Y di d = C ( z ) . Note, again, that every integer ≤ x ≤ N can beuniquely represented as x = 1 + D (cid:88) d =1 β ( x, d ) d − (cid:89) j =1 M j (37)for integers ≤ β ( x, d ) < M d and that φ ( x ) = 1 + D (cid:88) d =1 ( φ d ( β ( x,d ) + 1) − d − (cid:89) j =1 M j . (38)If C ( z ) is of the form in (28) then C ( z ) xy = α δ β ( x, β ( y, ,...,δ β ( x,D ) β ( y,D ) ( z ) . If C ( z ) is not of theform in (28) then there are some x , y and x , y suchthat δ β ( x ,d ) β ( y ,d ) = δ β ( x ,d ) β ( y ,d ) for d = 1 , . . . , D but that C ( z ) x ,y (cid:54) = C ( z ) x ,y . Let ψ , . . . , ψ D bepermutations such that ψ d ( β ( x , d ) + 1) − β ( x , d ) and ψ d ( β ( y , d ) + 1) − β ( y , d ) with correspond-ing row permutation matrices Q d . Let Q = Q ⊗ . . . Q D be the row permutation matrix corresponding to permu-tation ψ . Then (cid:0) QC ( z ) Q (cid:62) (cid:1) x y = C ( z ) ψ ( x ) ψ ( y ) = C ( z ) x y (cid:54) = C ( z ) x y . (39)Hence, C ( z ) is of the form in (28) if and only if C ( z ) = P C ( z ) P (cid:62) for all P = P ⊗ . . . ⊗ P D where P d is arow permutation matrix on M d symbols.Let φ d be a permutation function on M d symbolswith corresponding row permutation matrix P d for d =1 , . . . , D . Also, let P = P ⊗ . . . ⊗ P D , thus forminga permutation matrix that operates on each lattice di-mension. The following equations result from permutingthe rows and columns of C ( z ) and by manipulating thematrix inverse, expectations, and transposes. P C ( z ) P (cid:62) = P B − zI N − ...... (cid:32) δ lj N (cid:88) s =1 C ss ( z )E (cid:2) H js (cid:3)(cid:33) l,j = Nl,j =1 − P (cid:62) = P BP (cid:62) − zP I N P (cid:62) − ...... (cid:32) δ lj N (cid:88) s =1 C ss ( z )E (cid:104) H φ ( j ) s (cid:105)(cid:33) l,j = Nl,j =1 − (40)Note that for the underlying lattice graph, P A ( G lat ) P (cid:62) = A ( G lat ) . Hence W ( G perc ) hasthe same distribution as P W ( G perc ) P (cid:62) as all edgeinclusions are independent and identically distributed.Therefore, B = P BP (cid:62) . (41)Furthermore, in terms of matrix entries, W φ ( x ) φ ( y ) ( G perc ) has the same distribution as W xy ( G perc ) . Hence, H φ ( j ) s has the same distribution as H jφ − ( s ) . By applying these two identities and changingthe order of summation by φ , the following equations7esult. P C ( z ) P (cid:62) = B − zI N − ...... (cid:32) δ lj N (cid:88) s =1 C ss ( z )E (cid:104) H jφ − ( s ) (cid:105)(cid:33) l,j = Nl,j =1 − = B − zI N − ...... (cid:32) δ lj N (cid:88) s =1 C φ ( s ) φ ( s ) ( z )E (cid:2) H js (cid:3)(cid:33) l,j = Nl,j =1 − (42)The entries C φ ( k ) φ ( k ) ( z ) are the k th diagonal entries of P C ( z ) P (cid:62) for k = 1 , . . . , N and, thus, a solution to thesystem of equations (16). However, because the solutionis unique, this implies that C ( z ) = P C ( z ) P (cid:62) . (43)Therefore, C ( z ) is of the form shown in (28). (cid:4) While Theorm 2 specifies the form of the uniquesolution to the system of equations (16), it does notdescribe how to find the D parameters α i ,...,i D ( z ) that correspond to a given set of lattice dimensionsand link probability parameter that are necessary to findthe Stieltjes transform of the deterministically equivalentdistribution function F N . Corollary 1 shows that theStieltjes transform of F N is equal to the diagonal elementof C ( z ) described in Theorem 2. It also derives thesystem of D nonlinear equations to which the param-eters are a solution, accomplishing this by computingthe terms of (27) and noting that all computed matricesare simultaneously diagonalizable. The matrix equationin (27) then transforms into a system of equations thatdescribes the eigenvalues, which appears in (45) ofCorollary 1. Remark 1 describes use of an iterative ap-proach to find the solution. For added clarity, please notethat i , . . . , i D = 0 , index the parameters α i ,...,i D ( z ) while j , . . . , j D = 0 , index the equations that describethe parameters. Corollary 1 (Stieltjes Transform for D -Lattice Per-colation) The Stieltjes transform of the deterministicequivalent distribution function F n specified in Theorem1 for the empirical spectral distribution of W ( G perc ) isgiven by S F N ( z ) = α ,..., ( z ) , Im ( z ) (cid:54) = 0 (44) where the D complex variables α i ,...,i D ( z ) for i , . . . , i D = 0 , solves the system (cid:88) i ,...,i D =0 α i ,...,i D ( z ) D (cid:89) d =1 λ di d ( j d ) = (cid:32) pγ (cid:32) D (cid:88) d =1 λ d ( j d ) (cid:33) − z − ...... p (1 − p ) γ (cid:32) D (cid:88) d =1 ( M d − (cid:33) α ,..., ( z ) (cid:33) − (45) of D rational equations for j , . . . , j D = 0 , where λ di d ( j d ) = M d − i d = 0 , j d = 0 − i d = 0 , j d = 11 i d = 1 . (46) Proof:
Using the Stieltjes transform expression (15) fromTheorem 1 and the derived form for C ( z ) in (28) fromTheorem 2, the Stieltjes transform of the deterministicequivalent distribution function F N specified in Theorem1 for the empirical spectral distribution of W ( G perc ) iscomputed as S F N ( z ) = 1 N N (cid:88) k =1 C kk ( z )= 1 N N (cid:88) k =1 α ,..., ( z ) = α ,..., ( z ) . (47)In order to derive the system of equations describingthe coefficients α i ,...,i D for i , . . . , i D = 0 , , the valueof the third matrix in the denominator of (27) must firstbe computed. Subsequently, the system of equations maybe found by noting that all matrices in equation (16) maybe simultaneously diagonalized by orthogonal matricesof eigenvectors. Each distinct eigenvalue of C ( z ) willresult in an equation that the coefficients must satisfy. Let β be defined as in the proof of Theorem 2. Summing thevariance over each row yields N (cid:88) s =1 C ss ( z )E (cid:2) H js (cid:3) = N (cid:88) (cid:107) β ( s ) − β ( j ) (cid:107) =1 C ss ( z )E (cid:2) H js (cid:3) = p (1 − p ) γ (cid:32) D (cid:88) d =1 ( M d − (cid:33) α ,..., ( z ) (48)8nd, therefore, that (cid:32) δ lj N (cid:88) s =1 C ss ( z ) E (cid:2) H js (cid:3)(cid:33) l,j = Nl,j =1 = p (1 − p ) γ (cid:32) D (cid:88) d =1 ( M d − (cid:33) α ,..., ( z ) I N . (49)Note that C ( z ) , B , and the matrix computed in (49)are simultaneously diagonalizable by Kronecker productproperties. Let λ di d (0) be the eigenvalue of Y di d cor-responding to eigenvector v = . Let λ di d (1) be theeigenvalue of Y di d corresponding to any other eigenvec-tor orthogonal to v = . Hence, λ di d ( j d ) = M d − i d = 0 , j d = 0 − i d = 0 , j d = 11 i d = 1 . (50)The eigenvalues of C ( z ) are indexed by j , . . . , j D =0 , and are given by (cid:88) i ,...,i D =0 α i ,...,i D ( z ) D (cid:89) d =1 λ di d ( j d ) (51)The corresponding eigenvalues of B are given by pγ (cid:32) D (cid:88) d =1 λ d ( j d ) (cid:33) (52)The corresponding eigenvalues of the matrix in (49) aregiven by p (1 − p ) γ (cid:32) D (cid:88) d =1 ( M d − (cid:33) α ,..., ( z ) . (53)Of course, the corresponding eigenvalue of zI is z .Through applying this diagonalization to the equa-tion (16), the system of D equations indexed by j , . . . , j D = 0 , that the D coefficients α i ,...,i D indexed by i , . . . , i D = 0 , satisfy is described byequation (45). (cid:4) Remark 1 (Iterative Computation of Coefficients) Com-putation of the values of α i ,...,i D for i , . . . , i D = 0 , can be approached iteratively. Let the D dimensionalvector α collect the values of α i ,...,i D for i , . . . , i D =0 , . Then the system of equations in (45) can bewritten as X α = g ( α ) where X is the invertiblematrix describing the left side of the equation and g isthe function describing the right side of the equation.Selecting an initial guess α and iterating such that α n +1 = X − g ( α n ) is a possible approach to search for the fixed point α . While convergence of this processis not proven, it appears to work well in practice.Figure 2 and Figure 3 provide example deterministicequivalents for the empirical spectral distribution ofscaled adjacency matrices for percolation models oflattice graphs with different lattice characteristics andlink inclusion parameters as calculated using Corollary 1and the methods in Remark 1. The corresponding densityfunction is also computed, and the expected empiricalspectral distributions and densities are shown for com-parison. Asymptotically, as each lattice dimension sizegrows without bound, the area of the largest region of thedensity function approaches 1, and the area of the otherregions approach 0. Thus, Theorem 1 does not guar-antee the deterministic distribution function reflects theempirical spectral distribution in the regions of smallerincrease. Nevertheless, the observed characteristics seemto provide a good approximation.Finally, Theorem 3 describes the relationship betweenthe empirical spectral distribution of the row-normalizedadjacency matrix (cid:98) A ( G perc ) = ∆ − ( G perc ) A ( G perc ) ,where ∆ ( G perc ) is the diagonal matrix of node degrees,and the empirical spectral distribution of the scaledadjacency matrix W ( G perc ) = γ A ( G perc ) . Becausethe row-normalized adjacency matrix and the symmet-rically normalized adjacancy matrix are similar, theresult also applies to both of those matrices. Withsome manipulation, it can also be used to describethe row-normalized Laplacian and the symmetricallynormalized Laplacian. The theorem shows that the L´evymetric distance d L (cid:16) F p √ γ (cid:98) A , F √ γW (cid:17) between F p √ γ (cid:98) A and F √ γW approaches asymptotically. The proof followsthat of Lemma 5 of [13] almost identically, with onlysome complications in applying Lemma 2.3 of [17] tobound the matrix entry differences. Consequently, thedeterministic equivalent F N computed in Corollary 1 canbe used to approximate the empirical spectral distributionof (cid:98) A as well. Theorem 3
Let W ( G perc ) = γ A ( G perc ) be thescaled adjacency matrix of G perc ( G lat , p ) for scale factor γ = (cid:80) Dd =1 ( M d − with empirical spectral distribution F W , and let (cid:98) A ( G perc ) = ∆ − ( G perc ) A ( G perc ) be therow-normalized adjacency matrix of G perc ( G lat , p ) withempirical spectral distribution F (cid:98) A , where ∆ ( G perc ) isthe diagonal matrix of node degress. Also let d L ( · , · ) bethe L´evy distance metric. Assume that all of the latticedimension sizes increase without bound as N → ∞ .Then, lim N →∞ d L (cid:16) F p √ γ (cid:98) A , F √ γW (cid:17) = 0 . (54)9ig. 2: For a graph percolation model G perc based on a two dimensional lattice with size
20 nodes ×
40 nodes and link inclusion probability p = . , the expected empirical spectral distribution E [ F W N ] of the scaled adjacencymatrix W N = γ A ( G perc ) was computed over Monte-Carlo trials (black dotted line, top left). The deterministicdistribution F N computed through Corollary 1 was also computed (blue line, top left). The corresponding densityfunctions appear on the top right with labeled sections of interest more closely examined on the bottom row. Proof:
Following the proof methods of Lemma 5 of [13], therow sums ∆ ii are first computed as γ ( p + (cid:15) i ) where (cid:15) = max i | (cid:15) i | = o (1) . Subsequently, the bound onthe cubed L´evy distance found in Lemma 2.3 of [18]is evaluated and shown to asymptotically approach .Consequently, the result holds. The first of these steps isless straightforward than the equivalent in [13].Let β ( x, d ) be defined as before. For a given nodeindex i and lattice dimension d , form a M d × M d matrix Q ( i, d ) from all entries A lj such that β ( l, k ) = β ( j, k ) = β ( i, k ) for k (cid:54) = d . Note that this matrix issymmetric with zeros on the diagonal and identicallydistributed Bernoulli random variables with parameter p for all entries off the diagonal. These variables areindependent, except for symmetry relations. Form a M d × M d matrix R ( i, d ) equal to Q ( i, d ) on the offdiagonal but with an identical Bernoulli random variableadded to the diagonal element. By Lemma 2.3 of [17], M d (cid:88) s =1 R l,s ( i,d ) = M d ( p + δ d,l )= ( M d − p + δ i,d,l )+( p + δ i,d,l ) (55) where max l | δ i,d,l | = o (1) . Thus, note that M d (cid:88) s =1 Q l,s ( i, d ) = ( M d −
1) ( p + κ i,d,l ) (56)where ( p + δ i,d,l − /M d ≤ κ i,d,l ≤ ( p + δ i,d,l ) /M d so max l | κ i,d,l | = o (1) . Hence, ∆ ii may be computedas ∆ ii = N (cid:88) j =1 A ij = D (cid:88) d =1 M d (cid:88) s =1 Q β ( i,d ) ,s ( i, d )= D (cid:88) d =1 ( M d − (cid:0) p + κ i,d,β ( i,d ) (cid:1) = γ ( p + (cid:15) i ) (57)where (cid:15) i = 1 γ D (cid:88) d =1 ( M d − κ i,d,β ( i,d ) . (58)Note that max i | (cid:15) i | ≤ max i max d max l | κ i,d,l | (59)so (cid:15) = max i | (cid:15) i | = o (1) . By Lemma 2.3 of [18], d L (cid:16) F p √ γ (cid:98) A , F √ γW (cid:17) ≤ n (cid:13)(cid:13)(cid:13) √ γW − p √ γ (cid:98) A (cid:13)(cid:13)(cid:13) F . (60)10ig. 3: For a graph percolation model G perc based on a two dimensional lattice with size
15 nodes ×
20 nodes ×
25 nodes and link inclusion probability p = . , the expected empirical spectral distribution E [ F W N ] of the scaledadjacency matrix W N = γ A ( G perc ) was computed over Monte-Carlo trials (black dotted line, top left).The deterministic distribution F N computed through Corollary 1 was also computed (blue line, top left). Thecorresponding density functions appear on the top right with labeled sections of interest more closely examined onthe bottom row.Note that W ij = γ A ij and that p (cid:98) A ij = p ∆ ii A ij = pγ ( p + (cid:15) i ) A ij = 1 γ A ij (1 + O ( (cid:15) i ))= 1 γ A ij (1 + O ( (cid:15) )) (61)Hence, (cid:13)(cid:13)(cid:13) √ γW − p √ γ (cid:98) A (cid:13)(cid:13)(cid:13) F = γn N (cid:88) i,j =1 γ A ij O (cid:0) (cid:15) (cid:1) = 1 nγ N (cid:88) i,j =1 A ij O (cid:0) (cid:15) (cid:1) → (62)because (cid:15) = o (1) . Therefore, it has been shown that d L (cid:16) F p √ γ (cid:98) A , F √ γW (cid:17) → . (cid:4) IV. C
ONCLUSION
The eigenvalues of matrices that encode network in-formation, such as the adjacency matrix and the Lapla-cian, provide important information that can be used, for instance, in polynomial filter design for graph signalprocessing. This paper examines the empirical spectraldistribution of the scaled adjacency matrix eigenvaluesfor a random graph model formed by independently in-cluding or excluding each link of a D -dimensional latticesupergraph, a percolation model. A stochastic canonicalequations theorem provided by Girko can be applied tofind a sequence of deterministic functions that asymptot-ically approximates the empirical spectral distribution asthe number of nodes increases. Through this theorem, theStieltjes transforms of these deterministic distributionscan be found by solving a certain system of equations.The primary contributions of this paper focus on solv-ing this system of equations for the random networkmodel studied with arbitrary lattice and link inclusionparameters. Specifically, Theorem 2 finds the form of thematrix that satisfies the system of equations, resulting ina matrix with D parameters. Subsequently, Corollary1 derives a system of D equations that describe theseparameters, significantly reducing the original system ofequations. The resulting distributions can also be relatedto the empirical spectral distribution of the symmetricallynormalized adjacency matrix, row normalized adjacencymatrix, symmetrically normalized Laplacian, and row-normalized Laplacian. Simulations for sample model11arameters validate the results. While the methods usedguarantee that the computed deterministic distributionfunction asymptotically captures the behavior of the bulkof the eigenvalues, no guarantees are made about regionswhere the fraction of eigenvalues asymptotically van-ishes. Although the observed correspondence betweenthe computed and simulated distributions appears good,future efforts could include a more precise characteri-zation. Additional topics for possible expansion couldalso include evaluation of the degree to which filterdesign benefits from the distribution information gainedas well as analysis of additional random network modelsof interest. R EFERENCES[1] A. Sandryhaila and J. M. F Moura.
Discrete Signal Processingon Graphs , IEEE Transactions on Signal Processing, vol 61, no.7, pp. 1644-1656, April 2013.[2] D. Shuman, S. Narag, P. Frossard, A. Ortega, and P. Van-dergheynst.
The Emerging Field of Signal Processing on Graphs ,IEEE Signal Processing Magazine, vol. 30, no. 3, pp. 83-98, May2013.[3] A. Sandryhaila and J. M. F. Moura.
Big Data Analysis withSignal Processing on Graphs: Representation and Processingof Massive Data Sets with Irregular Structure , IEEE SignalProcessing Magazine, vol. 31, no. 5, pp. 80-90, 2014.[4] A. Sandryhaila, S. Kar, and J. M. F. Moura.
Finite-timeDistributed Consensus through Graph Filters , Proceedings ofICAASP 2014, pp. 1080-1084, May 2014.[5] A. Sandryhaila and J. M. F. Moura.
Discrete Signal Processingon Graphs: Frequency Analysis , IEEE Transactions on SignalProcessing, vol. 62, no. 12, pp. 3042-3054, 2014.[6] S. Chen, A. Sandryhaila, J. M. F. Moura, and J. Kovaˇcevi´c.
SignalRecovery on Graphs: Variation Minimization , IEEE Transactionson Signal Processing, vol. 63, no. 17, pp. 4609-4624, Sept. 2015.[7] E. Kokiopoulou and P. Frossard. “Polynomial Filtering for FastConvergence in Distributed Consensus”.
IEEE Transactions onSignal Processing , vol. 57, pp. 342-354, Jan 2009.[8] V. Girko.
Theory of Stochastic Canonical Equations , vol. 1, pp.1-3, Springer Science+Business Media, 2001.[9] E. Wigner.
On the Distribution of the Roots of Certain SymmetricMatrices , The Annals of Mathematics, vol. 67, no. 2, pp. 325-327, Mar. 1958.[10] R. Couillet and M. Debbah.
Random Matrix Methods for Wire-less Communications , pp. 29-31, 113-115, Cambridge UniversityPress, 2011.[11] G. Grimmett.
Percolation , second edition, pp. 1-12, Springer1999.[12] R. Laskar.
Eigenvalues of the Adjacency Matrix of Cubic LatticeGraphs , Pacific Journal of Mathematics, vol. 29, no. 3, pp. 623-629, July 1969.[13] K. Avrachenkov, L. Cottatellucci, and A. Kadavankandy.
SpectralProperties of Random Matrices for Stochastic Block Model , 4thInternational Workshop on Physics-Inspired Paradigms in Wire-less Communications and Networks, pp. 537-544, May 2015.[14] B. Bollob´as.
Random Graphs , second edition, pp. 34-35, Cam-bridge University Press, 2001.[15] E. Gilbert.
Random Graphs , The Annals of Mathematical Statis-tics, vol. 30, no. 4, pp. 1141-1144, 1959.[16] X. Ding and T. Jiang.
Spectral Distributions of Adjacency andLaplacian Matrices of Random Graphs , The Annals of AppliedProbability, vol. 20, no. 6, pp. 20862117, 2010. [17] C. Bordenave, P. Caputo, and D. Chafa¨ı.
Spectrum of Large Ran-dom Reversible Markov Chains: Two Examples , Latin AmericanJournal of Probability and Mathematical Statistics, vol. 7, pp.41-64, 2010.[18] Z. Bai.
Methodologies in Spectral Analysis of Large RandomMatrices, A Review , Statistica Sinica, vol. 9, no. 3, pp. 611-677,July 1999., Statistica Sinica, vol. 9, no. 3, pp. 611-677,July 1999.