Design of Sampling Set for Bandlimited Graph Signal Estimation

DESIGN OF SAMPLING SET FOR B ANDLIMITED GRAPH SIGNA L ESTIMA TION Xuan Xie, Hui F eng , Junlian Jia, Bo Hu Research Center of Smart Networks and Systems, Fudan Uni versity , Shanghai 200433, China Emails: { xxie15, hfeng, jljia10, bohu } @fudan.edu.cn ABSTRA C T It is of p ar ticular interest to reconstruct or estimate bandlimited graph signals, which are smoothly varying sig- nals d efined over gr ap hs, fro m partial noisy measurements. Howe ver , choosing an optimal subset of nodes to sample is NP-hard. W e for m ularize the problem as the experimental design o f a linear regression model if we allow mu ltiple measurements on a single node. By relaxing it to a conv ex optimization pr oblem, we get the p roportion o f sample f o r each node giv en the budget of total sample size. Then, we use a probab ilistic quantization to get the number of each node to be sampled. Moreover , we analyze how the sample size influences whether our object fun ction is well- defined by perturbation analysis. Finally , we demonstrate the perfo rmance of the pro posed approach throug h various numerical experiments. Index T erms — Graph S ignal, Sampling Theory , Conv ex Optimization, Experimental Design, Pertur b ation Analysis 1. INTR ODUCTION Graph signals (GS) are rath er useful for describin g sig- nals and information in irregular dom ain s, such as social, sensor and economic networks [1]. Graph signal pro cessing (GSP) th eory extends and gener alizes the classic discrete signal processing th eo ry to graphs by introducing graph Fourier transform [2]–[15], wavel et [16]–[20], etc. T wo main approaches h av e been considered to co nstruct a grap h Fourier analysis framework: algebraic graph theory linked to the adjenacy m atr ix [2] , [ 4 ]–[8] and spectral graph theory based o n th e g r aph Laplacian matrix [3], [8]–[1 5]. GSP theory hav e been developed to graph filter [ 13], [17], [20], community detection and clustering on graphs [11], [19], sampling and interpolation [3], [ 6], [7], [9], [10] and corre- sponding distrib uted algorithm [8], [12]. Recently , stationary processes to the gr ap h domain has b een inv estigated in [1 4], [15]. Sampling theor y for GS deals with the p roblem of re- covering a b an dlimited signal fro m its samples on a subset of nodes o f the graph. Bandlimited GS widely exist in most real-world networks due to the fact that the signals on adjacent nodes won’t differ dramatically in reality . T o for- mulate a sampling theor y fo r bandlimited GS, the fo llow ing questions arises: How to choose the best sampling set which can achieve the smallest reconstruction er ror? Generally , the signals can be bandlimited or appro ximate bandlimited and the samples can be noise-free o r noisy as shown in T able. 1. In this paper, we focus o n choosing the sampling set f or bandlimited GS from no isy measurements. T able 1 . Literature based on different assumptions Noise-free Noisy Bandlimited [6], [7] [3], [7]–[10], [12] Approximate bandlimit ed [6] [2], [3] In this work, we con nect the design of sampling set for noisy bandlimited GS with experimental design problem, which is closely related to the op timal experimental d esign problem in [21], [22]. Th is connection with experimental design has also been noticed in [3], [8] for GS sampling but not fully in vestigated. The original form of experimental design problem with given measurement size is a com bi- natorial problem and NP-hard. There are two main kinds of solutions. One is heuristic algorithms, which is time- consuming and can hardly produ ce an optimal solution, such as Fedorovs exchange algorithm [23] and W yn n s algorithm [24]. Th e other is the conv ex relaxation proposed in [25, Sec. 7.5] which tries to obtain an optimal solution. I t has been also claimed to giv e excellent results in m any app lications [26]. In this p aper , we fo llow the latter . In this work, we con sider the case when any node is allow ed to be sampled multiple times. Fro m the view of experimental d esign , the op timal solution u sually achiev es when we allow mu ltiple-time sampling. This is very feasible for many real-world networks such as sensor n etwork and social network. W e n ot o nly d ecide the nod es to be sampled but also get the prop ortion of sample for each nod e of a gi ven sample size by relaxing the combinational problem to con vex optimization. In order to get the number of each nod e to be sampled, we use probab ilistic quantization [27] to u n biased quantify the solution o f the relaxed prob lem to integers to get a subo ptimal solution of the combination al pro blem. Different from the heur istic algo r ithm [3], [8] whose perfor- mance can not be evaluated, we take the quantization erro r as the per tu rbation and analyze how sample size influence the performance o f our algorithm by pertu r bation analysis. Moreover , we find a lower bound of the sample size to ensure the object function of our algorithm well-defined, which can provide a refer en ce to practical p roblems. Finally , the perf o rmance of ou r method is analyzed and shown to hav e smaller reconstruction error an d is more robust against additiv e noise. 2. SYSTEM MODEL Consider an N -vertex und irected conn ected graph G = ( V , E , W ) , wher e V is the verte x set, E is the edge set, and W is the weigh ted adjacency matrix. If there is an ed ge e = ( i, j ) between vertices i and j , then the entry W i,j represents th e weight of the ed ge; otherwise W i,j = 0 . A signal f defined on the vertices of a graph can b e represented as a vector f ∈ R N , an d its element f i represents the signal value at the i th vertex in V . The graph L ap lacian is defined as L = D − W , where the degree matrix D = diag ( 1 W ) . Since the Laplacian matrix is real symmetric, it has a complete eigenbasis and the spectral decomposition L = V Λ V T , (1) where the eigenv ectors { v k } 0 ≤ k ≤ N − 1 of L form the columns of V , and Λ ∈ C N × N is a diagonal matr ix of eigen values 0 = λ 0 ≤ λ 1 ≤ · · · ≤ λ N − 1 of L . Th e g raph Fourier transfo r m corresponds to the basis expansion o f a signal. The eigenv ectors of the grap h Laplacian are regarded as the Fourier basis and the eigen values are regarded as frequencies [5]. The expansion coeffi cients of a graph signal f in terms of eigen vectors are defined as ˆ f , so that a graph signal can be expressed as f = V ˆ f . (2) A grap h signal is called band limited when there exis ts a K ∈ { 0 , 1 , · · · , N − 1 } such that its g raph Fourier transform ˆ f satisfies ˆ f k = 0 f o r all k ≥ K [ 7]. Let f be a signal with bandwidth K , then it satisfies f = V K ˆ f K , (3) where V K denotes the first K column s o f V an d ˆ f K denotes the first K coeffici ents of ˆ f . Suppose that we sample M measurements from the graph signal f ∈ C N to pr o duce a sampled signal f S ∈ C M , usually M ≤ N , where S = ( S 1 , · · · , S M ) denotes the sequence of sampled indices, and S i ∈ { 1 , 2 , · · · , N } . The sampling operator Ψ : C N 7→ C M is defined as Ψ i,j =  1 , j = S i ; 0 , otherwise. (4) Let w ∈ C M be the i.i .d. noise with zero mean and unit variance in troduced during sampling. T h en the samp les are g i ven by f S = Ψ f , and the observation mo del is y S = Ψ f + w . Considering ( 3), fo r bandlimited GS, the observation mo del can b e expressed as y S = Ψ V K ˆ f K + w . Let V M K = Ψ V K , then y S = V M K ˆ f K + w . (5) The best linear unbiased estimation [ 28] o f ˆ f K from ob- served samples y S is, ˆ f ′ K = V † M K y S , (6) where V † M K = ( V T M K V M K ) − 1 V T M K is the pseudo -in verse of V M K . By (2), th e estimation error of f i s e = f ′ − f = V K V † M K w . (7) The covariance m atr ix of estimation err or is E = E [ ee T ] = V K  V T M K V M K  − 1 V T K . (8) Our main problem is to cho ice an optimal Ψ that minimize the erro r covariance E in cer tain scalarization f orms. Se veral scalarizations of (8) hav e been proposed for the con venience of minimization as fo llow : D - optimal f ( E ) = log det ( E ) , (9) E - optimal f ( E ) = k E k 2 , (10) A - optimal f ( E ) = T r ( E ) . (11) It is obvious that f ( E ) = f   V T M K V M K  − 1  for all the three scalarizations above, since V K is orthogo nal. 3. ALGORITHM In this section, we consider the sampling set design prob- lem b y estimating f from measuremen ts. Let u T 1 , · · · , u T N be the rows of V K . Recall from V M K = Ψ V K that the the rows of V M K , which characterize the measurements, can be cho sen among N possible test vectors u T 1 , · · · , u T N . Our goal o f op timal experimental design is to make f ( E ) as small as possible. Let m i denotes th e number of experiments f or which u i is chosen, and assume that the sample size is M , so we have m 1 + · · · + m N = M . (12) The scalarization o f er r or cov ariance matrix can b e ex- pressed as f ( E ) = f   N X i =1 m i u i u T i ! − 1   . (13) It shows that the error cov ariance depen ds on ly o n th e numbers of each u i is ch o sen. The basic e xperimental d esign problem is as fo llow , minimize m i f ( E ) subject to m i ≥ 0 , m 1 + · · · + m N = M m i ∈ Z . (14) 3.1. THE RELAXED PROBLEM The basic experimen tal design problem (14) is an in- tractable co m binatorial pro blem. W e relax the con straint that the m i are integers following Boyd’ s method [25, Sec. 7.5]. Let p i = m i / M , which indicates the propo rtion of experiment i , an d relax the constraint that each p i is an integer multiple of 1 M , we obtain the relaxed experimental design problem minimize p i f ( E ) = f   N X i =1 p i u i u T i ! − 1   subject to p  0 , 1 T p = 1 . (15) This is a conve x op timization which can be solved by any optimization tool like interior -point methods [25]. In the rest of this paper , we consider only the relaxed experimental problem. The optimal v alue of the relaxed problem (15) provides a lower b ound on the optimal value of the com b i- natorial one since the combinatorial problem has an addition constraint. After solving (1 5), each en try o f p need to be quantified to an integer multiple of 1 M . Dif ferent form Boyd’ s method [25, Sec. 7.5], we u se probabilistic qu antization [27] instead of the simple roundin g to ensure an unbiased q uantization in mean. The pr obabilistic quan tization Q : p i → Q ( p i ) is defined as follow: p i ∈ [0 , 1] is equally divide into M − 1 sub-intervals. The quantization points is defined as { 1 / M , 2 / M , . . . , 1 } . Then, for p i ∈ [ k / M , ( k + 1) / M ] , k ∈ { 0 , 1 , . . . , M − 1 } , Q ( p i ) is a random variable defined by Q ( p i ) = ( k M with prob ability ( k +1 M − p i ) M k +1 M with prob ability ( p i − k M ) M . (16) Clearly we have | p i − Q ( p i ) | ≤ 1 / (2 M ) . So when M is large enough, we h a ve p ≈ Q ( p ) , wh ich implies the error cov ariance matrix associated with p and Q ( p ) are closed. As a result, we can use Q ( p ) to gener ate a sub o ptimal sampling set: let m i = M Q ( p i ) be the sample quota of the i th node. 3.2. PER TURBA TION AN AL YSIS OF QU ANTIZA- TION By relaxing the experimental design problem, we do not need to solve the combinatorial problem directly . Mean- while, the suboptimal solution brings a new problem of how to ensure the objective fu nction o f (15) in vertible when we using Q ( p i ) to replace p i , since the object function will be ill-condition if P N i =1 Q ( p i ) u i u T i is no t inv ertible. Let A = P N i =1 p i u i u T i and ˆ A = P N i =1 Q ( p i ) u i u T i , o u r goal is to ensure ˆ A , perturb ed from A , in vertible. Let the perturbation erro r ∆ p i = Q ( p i ) − p i , and the p erturbation is given by δ A = ˆ A − A = N X i =1 ∆ p i u i u T i . (17) According to the [29, Th 2 .1], known that A is nonsin- gular , the relativ e l 2 -norm distance fro m A to the nearest singular matrix is min  k ( δ A ) k 2 k A k 2 : A + δ A singular  = 1 k A k 2 k A − 1 k 2 . (18) Therefore, to en sure ˆ A inv ertible, the f ollowing con dition needs to be satisfied, k δ A k 2 < 1 k A − 1 k 2 = σ min ( A ) . (19) Recall f rom (17) that k δ A k 2 is related to ∆ p i . Thus, ∆ p i affects the probability o f (19) bein g held. Theorem 1. Suppose that ∆ p 1 , . . . , ∆ p N are independen t, and each Q ( p i ) is obtained b y the q uantization g iven in (16). The probability of (19) is P ( k δ A k 2 < σ min ( A )) > N Y i =1  1 − V ar [∆ p i ] ( σ min ( A )) 2  . (20) Pr oo f: Let ∆ P = d iag (∆ p 1 , ∆ p 2 , . . . , ∆ p N ) , then k δ A k 2 =   V T k ∆ P V k   2 ≤   V T k   2 k ∆ P k 2 k V k k 2 = max | ∆ p i | . (21) Since ∆ p 1 , . . . , ∆ p N are independent, we can get the fol- lowi ng ineq uation, P ( k δ A k 2 < σ min ( A )) ≥ N Y i =1 P ( | ∆ p i | < σ min ( A )) . (22) According to Chebycheff i nequality , for any p i ∈ [ k / M , ( k + 1) / M ] , the f ollowing ineq uation holds, P ( | ∆ p i | < σ min ( A )) > 1 − V ar [∆ p i ] ( σ min ( A )) 2 . (23) Thus pr oving our claim. Assume that M is large enough, then any p i ∈ [ k / M , ( k + 1) / M ] is ap proximated to be a un iform distribu- tion: p i ∼ U [ k / M , ( k + 1 ) / M ] . According to (16) , for every p i ∈ [ k / M , ( k + 1) / M ) , E [ Q ( p i )] = p i , which means Q ( p i ) is an u nbiased rep resentation of p i , so V ar [∆ p i ] = Z 1 2 M − 1 2 M ∆ p i 2 p (∆ p i ) d ∆ p i ≈ 5 192 M 3 . (24) W e do n ot need the p robability in (20) to be 1 most of the time, it can be reduced to η (0 < η < 1) according to actual needs. This leads to the following corollary . Corollary 1. ˆ A = P N i =1 Q ( p i ) u i u T i is in vertible with probability η if the sample size M satisfi es: M ≥ &  5 192(1 − N √ η )( σ min ( A )) 2  1 3 ' , (25) where ⌈·⌉ denotes the ceiling o peration. 10 12 14 16 18 20 Bandwidth 1 2 3 4 5 6 7 8 Reconstruction error 10 -4 Proposed M1 M2 M3 (a) Signal F1 0 2 4 6 8 10 SNR(dB) 0 1 2 3 4 5 6 Reconstruction error 10 -3 Proposed M1 M2 M3 (b) Signal F2 Fig. 1 . Recon struction results o f G1 fo r different signals. 4. SIMULA TION W e now numerically e val uate the per f ormance of the proposed work. The experiments compares the average re- construction error k e k 2 of different sample set selection algorithms versus different b andwidth of the true signal an d different SNR of the samples. W e compare our appr oach with the follow ing methods: M1 [7] and M2 [ 3], which are both greedy algorithm , as well as M3 [26], which uses relaxed conv ex optimization to calculate the probability to sample each n ode and select M no des with largest pr o ba- bility . W e give some simulation results on the fo llow ing sim- ulated undirected g r aphs: Small world grap h (G1) [30] of 1000 nodes rewiring with probab ility 0.1. Random geometric graph ( G2 ) [3 1] with 500 nodes placed randomly in th e unit square and edges are placed b etween any nodes within 0.6. The edge weights are assigned via a Gaussian kernel. For each of the above graphs, we con sider the pro blem in the following scenarios: ( F1) The bandwidth of the true signal varying fro m 10 to 20 and non-zero GFT coefficients are generated from N (1 , 0 . 5 2 ) . The samples are noisy with additiv e i.i.d. Gaussian no ise such that the SNR equ als 10 dB. (F2) The true signal is exactly with the bandwidth of 15 and non- zero GFT coefficients are generated from N (1 , 0 . 5 2 ) . The samples are n o isy with additiv e i.i.d. Gaus- sian noise and the SNR varies from 0dB to 10dB. 10 12 14 16 18 20 Bandwidth 2 4 6 8 10 12 14 Reconstruction error 10 -4 Proposed M1 M2 M3 (a) Signal F1 0 2 4 6 8 10 SNR(dB) 0 1 2 3 4 5 6 7 Reconstruction error 10 -3 Proposed M1 M2 M3 (b) Signal F2 Fig. 2 . Recon struction results o f G2 fo r different signals. W e generate 200 signals from each of the two signal models on each of the g r aphs and set the sample size to 4 times of the bandwidth. Th e parameter k of M2 that controls how closely the estimate the ban dwidth of a signal is set to 8 in ou r simulation. For no isy signal mod el F1, our method has a better per formance in reco nstruction error than the o thers with the same sample size in most case as d emonstrated in Fig. (1a) and Fig. (2a). This is becau se M1 and M2 are heuristic methods, which can not ensure a global optimal solution. M3 formularizes the p roblem as an experimen tal design p r oblem like what we do but not allo ws multiple-time sampling on any no de. So their solution is a suboptimal solution from the view o f measurement, while we try to find th e optimal one. The q uantization step is a major factor that influences the p erformance of ou r metho d . The smaller sample size M is, the larger the quantization step is. As a result, the corresponding reconstruction p erformance will be poorer . Our m ethod also has a better performance for signal model F2, especially when the SNR of samples is low as demonstrated in Fig. (1b) and Fig. (2 b). This also proves that our method is more robust against additiv e n oise. 5. CONCLUSION In this paper, an alg orithm is proposed to obtain the best sampling set for band limited GS estimation f or noisy samples. By the method of experimental design, we find out the sample q uota for each node. Meanwh ile, a loo se lower bound of sample size is given to ensure the object fun ction in our method is well-defined. Acknowledgment This work is supp orted by the NSF of China (Gran t No. 6150112 4). 6. REFERENCES [1] David I Shuman, Sun il K Narang, Pascal Frossard, An- tonio Ortega, and Pierre V andergheynst, “The emerging field of signal processing on graphs: Extending high- dimensional data analysis to networks and other irreg- ular domains, ” IEEE Signa l Pr ocess. Mag. , v ol. 30, no. 3, pp . 83–98, 2013. [2] Siheng Chen, Rohan V arma, Aarti Singh, and Jelena K ov a ˇ cevi ´ c, “Signal recovery on g r aphs: Fun damental limits of sampling strategies, ” IEEE T rans. Sign al Inf. Pr ocess. Netw . , vol. 2 , no. 4, pp . 539–55 4, 2016. [3] Aamir Anis, Akshay Gadde, and Antonio Ortega, “Ef - ficient sampling set selection for bandlimited grap h signals using graph spectral prox ies, ” I EEE T rans. Signal Pr ocess. , vol. 6 4, no. 14, pp. 3775–37 8 9, 201 6. [4] Aliaksei Sandryh aila and Jos ´ e MF Mour a, “Discrete signal processing on g raphs, ” IEEE T rans. Signal Pr ocess. , vol. 61, no. 7 , pp. 1644–1 6 56, 2013. [5] Aliaksei Sandryhaila and Jos ´ e MF Moura, “Dis- crete signal processing on graphs: Frequency analysis, ” IEEE T rans. Signal Pr ocess. , vol. 62, no. 12, pp. 3042– 3054, 20 14. [6] S. K. Narang, A. Gadde, and A. Ortega, “Signal pro- cessing techniqu es for interpolation in grap h structured data, ” in Pr oc. IEEE Int. Conf. Acoust., Speech, S ignal Pr ocess. (ICA SSP) , 2013, p p . 5445–5 449. [7] Siheng Chen , Rohan V arma, Aliaksei Sandryhaila, and Jelena K ov a ˇ cevi ´ c, “Discrete signal pro cessing o n graphs: Sampling theory , ” IEEE T rans. Signal Pr ocess. , vol. 63, no. 24, pp . 6510–6 523, 2015. [8] Antonio G Marques, Santiago Segarra, Geert L eu s, and Alejandro Ribeiro, “Sampling of graph signals with successi ve local aggregations, ” IEEE T rans. Signa l Pr ocess. , vol. 64, no. 7 , pp. 1832–1 8 43, 2016. [9] Gilles Puy , Nicolas T remblay , R ´ emi Gribonv al, and Pierre V andergheynst, “Ran dom sampling of bandlim- ited signals on grap hs, ” Ap pl. Compu t. Ha rmon. An al. , 2016. [10] Mikhail Tsitsvero, Sergio Barb arossa, and Paolo Di Lo renzo, “Sign als on graphs: Uncertainty principle and sampling, ” IEEE T rans. Sig nal Pr ocess. , vol. 64, no. 1 8, pp. 48 45–4860 , 20 1 6. [11] Xiao wen Don g, Pascal Frossard, Pierre V an d ergheynst, and Nikolai Nefedov , “C lustering on multi-layer gr aphs via subspace analysis on grassmann m anifolds, ” IEEE T rans. S ignal Pr ocess. , vol. 6 2, no. 4, pp. 905–918, 2014. [12] Xiaohan W ang, Pengfei Liu, and Y uantao Gu, “L o cal- set-based graph signal reconstruction, ” IEEE T rans. Signal Pr ocess. , vol. 63 , no. 9, p p. 2432– 2444, 2015. [13] Xuesong Shi, Hui Feng, Muyuan Zhai, T ao Y ang, and Bo Hu, “Infinite impulse response graph filters in wireless sensor networks, ” I EEE Signal Pr ocess. Lett. , vol. 22, no. 8, pp. 111 3–1117, 2015. [14] Benjamin Girault, “Stationary g raph signals using an isometric graph translation, ” in Pr oc. Eur opean S ignal Pr ocessing Co nf. (EUSIPCO) , 2015 , pp. 151 6–1520. [15] Nathana ¨ el Perraudin and Pierre V andergheynst, “Sta- tionary signal pro cessing on graphs, ” IEEE T rans. Signal Pr ocess. , vol. 6 5, no . 13, pp. 3462– 3477, 2017. [16] Da vid K Hamm ond, Pierre V and ergheynst, and R ´ emi Gribonv al, “W avel ets on graphs via spectral grap h theory , ” App l. C omput. Harmon. Anal. , v ol. 30, no. 2, pp . 129–15 0, 2011. [17] Sunil K Narang and Antonio Ortega, “Com pact sup - port biorthogonal wa velet filterbanks for ar bitrary und i- rected graphs, ” IEEE T rans. S ignal Pr ocess. , vol. 61, no. 19, pp. 4 673–468 5, 20 13. [18] Da vid I Shuman, Christoph Wiesme yr , Nicki Ho- lighaus, and Pierre V andergheynst, “Spectrum-adapted tight graph wa velet and vertex-frequency frames, ” IEEE T rans. Sign al Pr ocess. , vol. 63, no. 1 6, pp. 4223– 4235, 20 15. [19] Nicolas T rem blay and Pierre B o rgnat, “Graph wa velets for multiscale community mining, ” IEEE T rans. Signal Pr ocess. , vol. 62, no. 20, pp. 5 2 27–5239 , 20 1 4. [20] Sunil K Narang and Anton io Ortega, “Perfect recon- struction two-channel wa velet filter banks for grap h structured d ata, ” IEEE T rans. Signal Pr ocess. , vol. 60, no. 6, pp. 27 86–2799 , 201 2. [21] V alerii V ad imovich Fedorov , Theory of o ptimal exper- iments , Elsevier , 1 972. [22] Friedrich Pukelsheim, Optimal design o f experiments , SIAM, 2006 . [23] Alan J Miller and Nam-K y Nguyen , “ Algorithm as 2 9 5: A fed orov exchange algor ithm for d -optimal design, ” Appl. Statist. , vol. 43, no . 4, p p. 669–6 77, 1994. [24] Henry P W y n n, “R esults in the theory and construction of d- optimum experimental designs, ” J . Roy . Statist. Soc. Series B (Methodological) , vol. 34, no. 2, pp. 133– 147, 197 2. [25] Stephen Boyd an d Lieven V andenberghe, Conve x op - timization , Cambr idge univ ersity press, 2004. [26] Siddharth Joshi and Stephen Boyd, “Sensor selection via con vex optimization, ” IEEE T rans. Signal Pr ocess. , vol. 57, no. 2, pp. 451 –462, 2009. [27] T un cer C A ysal, Mark C oates, and Michael Rab- bat, “Distrib uted average consensus u sing probab ilistic quantization, ” in Pr oc. I EEE/SP W orkshop on Statist. Signal Pr ocess. W orkshop (S SP) , 2007, pp . 640–64 4. [28] Ste ven M Kay , Fundamentals o f statistical signal pr ocessing, volume I: Es timation Theory , Prentice Hall, 1993. [29] James W Demmel, Applied numerical linear algebra , SIAM, 1997 . [30] Duncan J W atts an d Ste ven H Strogatz, “Collective dynamics of small-worldnetworks, ” n atur e , vol. 393, no. 6684 , pp. 440 –442, 1998. [31] Nathana ¨ el Perrau din, Johan Paratte, David Shuman, Li- onel Martin, V assilis Kalofolias, Pierre V andergheynst, and Da vid K Hammond, “Gspbox: A toolbo x for signal processing on gr ap hs, ” arXiv prepri nt arXiv:1408.57 81 , 2 014.