Chebyshev Polynomial Approximation for Distributed Signal Processing
Unions of graph Fourier multipliers are an important class of linear operators for processing signals defined on graphs. We present a novel method to efficiently distribute the application of these operators to the high-dimensional signals collected by sensor networks. The proposed method features approximations of the graph Fourier multipliers by shifted Chebyshev polynomials, whose recurrence relations make them readily amenable to distributed computation. We demonstrate how the proposed method can be used in a distributed denoising task, and show that the communication requirements of the method scale gracefully with the size of the network.
Chebyshev polynomial approximation, denoising, distributed optimization, regularization, signal processing on graphs, spectral graph theory, wireless sensor networks
Wireless sensor networks are now prevalent in applications such as environmental monitoring, target tracking, surveillance, medical diagnostics, and manufacturing process flow. The sensor nodes are often deployed en masse to collectively achieve tasks such as estimation, detection, classification, and localization. While such networks have the ability to collect large amounts of data in a short time, they also face a number of resource constraints. First, they are energy constrained, as they are often expected to operate for long periods of time without human intervention, despite being powered by batteries or energy harvesting. Second, they may have limited communication range and capacity due to the need to save energy. Third, they may have limited on-board processing capabilities. Therefore, it is critical to develop distributed algorithms for in-network data processing that help balance the trade-offs between performance, communication bandwidth, and computational complexity.
Due to the limited communication range of wireless sensor nodes, each sensor node in a large network is likely to communicate with only a small number of other nodes in the network. To model the communication patterns, we can write down a graph with each vertex corresponding to a sensor node and each edge corresponding to a pair of nodes that communicate. Moreover, because the communication graph is a function of the distances between nodes, it often captures spatial correlations between sensors’ observations as well. That is, if two sensors are close enough to communicate, their observations are more likely to be correlated. We can further specify these spatial correlations by adding weights to the edges of the graph, with higher weights associated to edges connecting sensors with closely correlated observations. For example, it is common to construct the graph with a thresholded Gaussian kernel weighting function based on the physical distance between nodes, where the weight of edge connecting nodes and that are a distance apart is
for some parameters and .
In this paper, we consider signals collected by a sensor network whose nodes can only send messages to their local neighbors (i.e., they cannot communicate directly with a central entity). While much of the literature on distributed signal processing (see, e.g., - and references therein) focuses on coming to an agreement on simple features of the observed signal (e.g., consensus averaging, parameter estimation), we are more interested in processing the full function in a distributed manner, with each node having its own objective. Some example tasks under this umbrella include:
Distributed denoising – In a sensor network of sensors, a noisy -dimensional signal is observed, with each component of the signal corresponding to the observation at one sensor location. Using the prior knowledge that the denoised signal should be smooth or piecewise smooth with respect to the underlying weighted graph structure, the sensors’ task is to denoise each of their components of the signal by iteratively passing messages to their local neighbors and performing computations.
Distributed semi-supervised learning / binary classification – A binary label (-1 or 1) is associated with each sensor node; however, only a small number of nodes in the network have knowledge of their labels. The cooperative task is for each node to learn its label by iteratively passing messages to its local neighbors and performing computations.
These and similar tasks have been considered in centralized settings in the relatively young field of signal processing on graphs. For example, - consider general regularization frameworks on weighted graphs; - present graph-based semi-supervised learning methods; and - consider regularization and filtering on weighted graphs for image and mesh processing. In distributed settings,  considers denoising via wavelet processing and  presents a denoising algorithm that projects the measured signal onto a low-dimensional subspace spanned by smooth functions. References - consider different distributed regression problems.
Our main contributions in this paper are i) to show that a key component of many distributed signal processing tasks is the application of linear operators that are unions of graph Fourier multipliers; and ii) to present a novel method to efficiently distribute the application of the graph Fourier multiplier operators to the high-dimensional signals collected by sensor networks.
To elaborate a bit, graph Fourier multiplier operators are the graph analog of filter banks, one of the most commonly used tools in digital signal processing. Multiplying a signal on the graph by one of these matrices is analogous to reshaping the signal’s frequencies by multiplying it by a filter in the Fourier domain in classical signal processing. The crux of our novel distributed computational method is to approximate each graph Fourier multiplier by a truncated Chebyshev polynomial expansion. In a centralized setting,  shows that the truncated Chebyshev polynomial expansion efficiently approximates the application of a spectral graph wavelet transform, which is a specific example of a union of graph Fourier multipliers. In this paper, we extend the Chebyshev polynomial approximation method to the general class of unions of graph Fourier multiplier operators, and show how the recurrence properties of the Chebyshev polynomials also enable distributed application of these operators. The communication requirements for distributed computation using this method scale gracefully with the number of sensors in the network (and, accordingly, the size of the signals).
The remainder of the paper is as follows. In the next section, we provide some background from spectral graph theory. In Section III, we introduce graph Fourier multiplier operators and show how they can be efficiently approximated with shifted Chebyshev polynomials in a centralized setting. We then discuss the distributed computation of quantities involving these operators in Section IV, and provide some application examples in Section V. Section VI concludes the paper.
Ii Spectral Graph Theory
Before proceeding, we introduce some basic notations and definitions from spectral graph theory . We model the sensor network with an undirected, weighted graph , which consists of a set of vertices , a set of edges , and a weight function that assigns a non-negative weight to each edge. We assume the number of sensors in the network, , is finite, and the graph is connected. The adjacency (or weight) matrix for a weighted graph is the matrix with entries , where
The degree of each vertex is the sum of the weights of all the edges incident to it. We define the degree matrix to have diagonal elements equal to the degrees, and zeros elsewhere. The non-normalized graph Laplacian is defined as . For any on the vertices of the graph, satisfies
where indicates vertices and are connected.
As the graph Laplacian is a real symmetric matrix, it has a complete set of orthonormal eigenvectors. We denote these by for , with associated real, non-negative eigenvalues satisfying . Zero appears as an eigenvalue with multiplicity equal to the number of connected components of the graph . Without loss of generality, we assume the eigenvalues of the Laplacian of the connected graph to be ordered as
Just as the classical Fourier transform is the expansion of a function in terms of the eigenfunctions of the Laplace operator
the graph Fourier transform of any function on the vertices of is the expansion of in terms of the eigenfunctions of the graph Laplacian. It is defined by
where we adopt the convention that the inner product be conjugate-linear in the first argument. The inverse graph Fourier transform is given by
Iii Chebyshev Polynomial Approximation of Graph Fourier Multipliers
In this section, we introduce graph Fourier multiplier operators, unions of graph Fourier multiplier operators, and a computationally efficient approximation to unions of graph Fourier multiplier operators based on shifted Chebyshev polynomials. All methods discussed here are for a centralized setting, and we extend them to a distributed setting in Section IV.
Iii-a Graph Fourier Multiplier Operators
For a function defined on the real line, a Fourier multiplier operator or filter reshapes the function’s frequencies through multiplication in the Fourier domain:
Equivalently, denoting the Fourier and inverse Fourier transforms by and , we have
We can extend this straightforwardly to functions defined on the vertices of a graph (and in fact to any group with a Fourier transform) by replacing the Fourier transform and its inverse in (4) with the graph Fourier transform and its inverse, defined in (2) and (3). Namely, a graph Fourier multiplier operator is a linear operator that can be written as
We refer to as the multiplier. A high-level intuition behind (III-A) is as follows. The eigenvectors corresponding to the lowest eigenvalues of the graph Laplacian are the “smoothest” in the sense that is small for neighboring vertices and . At the extreme is , which is a constant vector ( for all and ). The inverse graph Fourier transform (3) provides a representation of a signal as a superposition of the orthonormal set of eigenvectors of the graph Laplacian. The effect of the graph Fourier multiplier operator is to modify the contribution of each eigenvector. For example, applying a multiplier that is 1 for all below some threshold, and 0 for all above the threshold is equivalent to projecting the signal onto the eigenvectors of the graph Laplacian associated with the lowest eigenvalues. This is analogous to low-pass filtering in the continuous domain. Section V contains further intuition about and examples of graph Fourier multiplier operators. For more properties of the graph Laplacian eigenvectors, see .
Iii-B Unions of Graph Fourier Multiplier Operators
In order for our distributed computation method of the next section to be applicable to a wider range of applications, we can generalize slightly from graph Fourier multipliers to unions of graph Fourier multiplier operators. A union of graph Fourier multiplier operators is a linear operator () whose application to a function can be written as (see also Figure 1)
where for every , is a graph Fourier multiplier operator with multiplier , and
Iii-C The Chebyshev Polynomial Approximation
Exactly computing requires explicit computation of the entire set of eigenvectors and eigenvalues of , which becomes computationally challenging as the size of the network, , increases, even in a centralized setting. As discussed in detail in [20, Section 6], a computationally efficient approximation of can be computed by approximating each multiplier by a truncated series of shifted Chebyshev polynomials. Doing so circumvents the need to compute the full set of eigenvectors and eigenvalues of . We summarize this approach below.
For , the Chebyshev polynomials are generated by
These Chebyshev polynomials form an orthogonal basis for
. So every function on that is square integrable with respect to the measure can be represented as , where is a sequence of Chebyshev coefficients that depends on . For a detailed overview of Chebyshev polynomials, including the above definitions and properties, see -.
By shifting the domain of the Chebyshev polynomials to via the transformation , we can represent each multiplier as
For , the shifted Chebyshev polynomials satisfy
Thus, for any , we have
where and the element of is given by
Now, to approximate the operator , we can approximate each multiplier by the first terms in its Chebyshev polynomial expansion (7). Then, for every and , we have
To recap, we propose to compute by first computing the Chebyshev coefficients according to (III-C), and then computing the sum in (III-C). The computational benefit of the Chebyshev polynomial approximation arises in (III-C) from the fact the vector can be computed recursively from and according to (9). The computational cost of doing so is dominated by the matrix-vector multiplication of the graph Laplacian , which is proportional to the number of edges, . Therefore, if the underlying communication graph is sparse (i.e., scales linearly with the network size ), it is far more computationally efficient to compute than . Finally, we note that in practice, setting the Chebyshev approximation order to around 20 results in approximating very closely in all of the applications we have examined.
Iv Distributed Computation
In the previous section, we showed that the Chebyshev polynomial approximation to a union of graph Fourier multipliers provides computational efficiency gains, even in a centralized computation setting. In this section, we discuss the second benefit of the Chebyshev polynomial approximation: it is easily distributable.
Iv-a Distributed Computation of
We consider the following scenario. There is a network of nodes, and each node begins with the following knowledge:
, the component of the signal
The identity of its neighbors, and the weights of the graph edges connecting itself to each of its neighbors
An upper bound on , the largest eigenvalue of the graph Laplacian. This bound need not be tight, so we can precompute a bound such as , where is the degree of node 
The task is for each network node to compute
by iteratively exchanging messages with its local neighbors in the network and performing some computations.
As a result of (III-C), for node to compute the desired sequence in (12), it suffices to learn . Note that and for all nodes that are not neighbors of node . Thus, to compute , sensor node just needs to receive from all neighbors . So once all nodes send their component of the signal to their neighbors, they are able to compute their respective components of . In the next step, each node sends the newly computed quantity to all of its neighbors, enabling the distributed computation of according to (9). The iterative process of computation and local communication continues for rounds until each node has computed the required sequence . In all, messages of length 1 are required for every node to compute its sequence of coefficients in (12) in a distributed fashion. This distributed computation process is summarized in Algorithm 1.
Inputs at node : , , ,
Outputs at node :
An important point to emphasize again is that although the operator and its approximation are defined through the eigenvectors of the graph Laplacian, the Chebyshev polynomial approximation helps the sensor nodes apply the operator to the signal without explicitly computing (individually or collectively) the eigenvalues or eigenvectors of the Laplacian, other than the upper bound on its spectrum. Rather, they initially communicate their component of the signal to their neighbors, and then communicate simple weighted combinations of the messages received in the previous stage in subsequent iterations. In this way, information about each component of the signal diffuses through the network without direct communication between non-neighboring nodes.
Iv-B Distributed Computation of
The application of the adjoint of the Chebyshev polynomial approximate operator can also be computed in a distributed manner. Let
where . Then it is straightforward to show that
We assume each node starts with knowledge of for all . For each , the distributed computation of the corresponding term on the right-hand side of (13) is done in an analogous manner to the distributed computation of discussed above. Since this has to be done for each , messages, each a vector of length , are required for every node to compute .
Iv-C Distributed Computation of
V Application Examples
In this section, we provide more detailed explanations of how the Chebyshev polynomial approximation of graph Fourier multipliers can be used in the context of specific distributed signal processing tasks.
V-a Distributed Smoothing
Perhaps the simplest example application is distributed smoothing with the heat kernel as the graph Fourier multiplier. One way to smooth a signal is to compute , where, for a fixed , . clearly satisfies our definition of a graph Fourier multiplier operator (with ). In the context of a centralized image smoothing application,  discusses in detail the heat kernel, , and its relationship to classical Gaussian filtering. Similar to the example at the end of Section III-A, the main idea is that the multiplier acts as a low-pass filter that attenuates the higher frequency (less smooth) components of .
Now, to perform distributed smoothing, we just need to compute in a distributed manner according to Algorithm 1, where is the shifted Chebyshev polynomial approximation to the graph Fourier multiplier operator .
V-B Distributed Regularization
Regularization is a common signal processing technique to solve ill-posed inverse problems using a priori information about a target signal to recover it accurately. Here we use regularization to solve the distributed denoising task discussed in Section I, starting with a noisy signal defined on a graph of sensors. The prior belief we want to enforce is that the target signal is smooth with respect to the underlying graph topology. The class of regularization terms we consider is for , and the resulting regularization problem has the form
To see intuitively why incorporating such a regularization term into the objective function encourages smooth signals (with as an example), note that if and only if is constant across all vertices, and, more generally
so is small when the signal has similar values at neighboring vertices with large weights (i.e., it is smooth).
We now show how our novel method is useful in solving this distributed regularization problem.
The objective function in (14) is convex in . Differentiating it with respect to , any solution to
From the real, symmetric nature of and the definition of the Laplacian eigenvectors (), we have:
Finally, taking the inverse graph Fourier transform of (18), we have
So, one way to do distributed denoising is to compute , the Chebyshev polynomial approximation of , in a distributed manner via Algorithm 1. We show this now with a numerical example. We place 500 sensors randomly in the square. We then construct a weighted graph according to the thresholded Gaussian kernel weighting (1) with and , so that two sensor nodes are connected if their physical separation is less than 0.075. We create a smooth 500-dimensional signal with the component given by , where and are node ’s and coordinates in . One instance of such a network and signal are shown in Figure 2, and the eigenvectors of the graph Laplacian are shown in Figure 3.
Next, we corrupt each component of the signal with uncorrelated additive Gaussian noise with mean zero and standard deviation 0.5. Then we apply the graph Fourier multiplier operator , the Chebyshev polynomial approximation to from Proposition 1, with . The multiplier and its Chebyshev polynomial approximations are shown in Figure 4, and the denoised signal is shown in Figure 5. We repeated this entire experiment 1000 times, with a new random graph and random noise each time, and the average mean square error for the denoised signals was 0.013, as compared to 0.250 average mean square error for the noisy signals.
We conclude this section by returning to the distributed binary classification task discussed in the introduction. In , Belkin et al. show that the regularizer also works well in graph-based semi-supervised learning. One approach to distributed binary classification is to let be the labels (-1 or 1) of those nodes who know their labels, and 0 otherwise. Then the nodes compute in a distributed manner via Algorithm 1, and each node sets it label to 1 if and -1 otherwise. We believe our approach to distributedly applying graph Fourier multipliers can also be used for more general learning problems, but we leave this for future work.
V-C Distributed Wavelet Denoising
In this section, we consider an alternate method of distributed denoising that may be better suited to signals that are piecewise smooth on the graph, but not necessarily globally smooth. The setup is the same as in Section V-B, with a noisy signal , and each sensor observing . Instead of starting with a prior that the signal is globally smooth, we start with a prior belief that the signal is sparse in the spectral graph wavelet domain . The spectral graph wavelet transform, , defined in , is precisely of the form of in (6). Namely, it is composed of one multiplier, , that acts as a low-pass filter to stably represent the signal’s low frequency content, and wavelet operators, defined by , where is a set of scales and is the wavelet multiplier that acts as a band-pass filter.
The most common way to incorporate a sparse prior in a centralized setting is to regularize via a weighted version of the least absolute shrinkage and selection operator (lasso) , also called basis pursuit denoising :
where . The optimization problem in (20) can be solved for example by iterative soft thresholding . The initial estimate of the wavelet coefficients is arbitrary, and at each iteration of the soft thresholding algorithm, the update of the estimated wavelet coefficients is given by
where is the step size and is the shrinkage or soft thresholding operator
We now turn to the issue of how to implement the above algorithm in a distributed fashion by sending messages between neighbors in the network. One option would be to use the distributed lasso algorithm of , which is a special case of the Alternating Direction Method of Multipliers [31, p. 253]. In every iteration of that algorithm, each node transmits its current estimate of all the wavelet coefficients to its local neighbors. With a transform the size of the spectral graph wavelet transform, this requires total messages at every iteration, with each message being a vector of length . A method where the amount of communicated information does not grow with (beyond the number of edges, ) would be highly preferable.
The Chebyshev polynomial approximation of the spectral graph wavelet transform allows us to accomplish this goal. Our approach is to approximate by , and use the distributed implementation of the approximate wavelet transform and its adjoint to perform iterative soft thresholding. In the first soft thresholding iteration, each node must learn at all scales , via Algorithm 1. These coefficients are then stored for future iterations. In the iteration, each node must learn the coefficients of centered at , by sequentially applying the operators and in a distributed manner via the methods of Sections IV-B and IV-A, respectively. Finally, when a stopping criterion for the soft thresholding is satisfied, the adjoint operator is applied again in a distributed manner to the resulting coefficients , and node ’s denoised estimate of its signal is .
We now examine the communication requirements of this approach. Recall from Section IV that messages of length 1 are required to compute in a distributed fashion. Distributed computation of , the other term needed in the iterative thresholding update (V-C), requires messages of length and messages of length . The final application of the adjoint operator to recover the denoised signal estimates requires another messages, each a vector of length . Therefore, the Chebyshev polynomial approximation to the spectral graph wavelet transform enables us to iteratively solve the weighted lasso in a distributed manner where the communication workload only scales with the size of the network through , and is otherwise independent of the network dimension .
Vi Concluding Remarks and Future Work
We presented a novel method to distribute a class of linear operators called unions of graph Fourier multiplier operators. The main idea is to approximate the graph Fourier multipliers by Chebyshev polynomials, whose recurrence relations make them readily amenable to distributed computation in a sensor network. Key takeaways from the discussion and application examples include:
A number of distributed signal processing tasks can be represented as distributed applications of unions of graph Fourier multiplier operators (and their adjoints) to signals on weighted graphs. Examples include distributed smoothing, denoising, and semi-supervised learning.
The graph Fourier multiplier operators are the graph analog of filter banks, as they reshape functions’ frequencies through multiplication in the Fourier domain.
The amount of communication required to perform the distributed computations only scales with the size of the network through the number of edges of the communication graph, which is usually sparse. Therefore, the method is well suited to large-scale sensor networks.
Our ongoing work includes extending the scope and depth of our application examples. In addition to considering more applications and larger size networks, we plan a more thorough empirical comparison of the computation and communication requirements of the approach described in this paper to alternative distributed optimization methods. The second major line of ongoing work is to analyze robustness issues that arise in real networks. For instance, we would like to incorporate quantization and communication noise into the sensor network model, in order to see how these propagate when using the Chebyshev polynomial approximation approach to distributed signal processing tasks. It is also important to analyze the effects of a sensor node dropping out of the network or communicating nodes losing synchronicity to ensure that the proposed method is stable to these disturbances.
-  M. Rabbat and R. Nowak, “Distributed optimization in sensor networks,” in Proc. Int. Symp. Inf. Process. Sensor Netw., Berkeley, CA, Apr. 2004, pp. 20–27.
-  J. B. Predd, S. R. Kulkarni, and H. V. Poor, “Distributed learning in wireless sensor networks,” IEEE Signal Process. Mag., vol. 23, pp. 56–69, Jul. 2006.
-  R. Olfati-Saber, J. Fax, and R. Murray, “Consensus and cooperation in networked multi-agent systems,” Proc. IEEE, vol. 95, no. 1, pp. 215–233, Jan. 2007.
-  A. G. Dimakis, S. Kar, J. M. F. Moura, M. G. Rabbat, and A. Scaglione, “Gossip algorithms for distributed signal processing,” Proc. IEEE, vol. 98, no. 11, pp. 1847–1864, Nov. 2010.
-  A. J. Smola and R. Kondor, “Kernels and regularization on graphs,” in Proc. Ann. Conf. Comp. Learn. Theory, ser. Lect. Notes Comp. Sci., B. Schölkopf and M. Warmuth, Eds. Springer, 2003, pp. 144–158.
-  D. Zhou and B. Schölkopf, “A regularization framework for learning from graph data.” in Proc. ICML Workshop Stat. Relat. Learn. and Its Connections to Other Fields, Jul. 2004, pp. 132–137.
-  ——, “Regularization on discrete spaces,” in Pattern Recogn., ser. Lect. Notes Comp. Sci., W. G. Kropatsch, R. Sablatnig, and A. Hanbury, Eds. Springer, 2005, vol. 3663, pp. 361–368.
-  X. Zhu and Z. Ghahramani, “Semi-supervised learning using Gaussian fields and harmonic functions,” in Proc. Int. Conf. Mach. Learn., Washington, D.C., Aug. 2003, pp. 912–919.
-  M. Belkin, I. Matveeva, and P. Niyogi, “Regularization and semi-supervised learning on large graphs,” in Learn. Theory, ser. Lect. Notes Comp. Sci. Springer-Verlag, 2004, pp. 624–638.
-  D. Zhou, O. Bousquet, T. N. Lal, J. Weston, and B. Schölkopf, “Learning with local and global consistency,” in Adv. Neural Inf. Process. Syst., S. Thrun, L. Saul, and B. Schölkopf, Eds. MIT Press, 2004, pp. 321–328.
-  S. Bougleux, A. Elmoataz, and M. Melkemi, “Discrete regularization on weighted graphs for image and mesh filtering,” in Scale Space Var. Methods Comp. Vision, ser. Lect. Notes Comp. Sci., F. Sgallari, A. Murli, and N. Paragios, Eds. Springer, 2007, vol. 4485, pp. 128–139.
-  A. Elmoataz, O. Lezoray, and S. Bougleux, “Nonlocal discrete regularization on weighted graphs: a framework for image and manifold processing,” IEEE Trans. Image Process., vol. 17, pp. 1047–1060, Jul. 2008.
-  F. Zhang and E. R. Hancock, “Graph spectral image smoothing using the heat kernel,” Pattern Recogn., vol. 41, pp. 3328–3342, Nov. 2008.
-  G. Peyré, S. Bougleux, and L. Cohen, “Non-local regularization of inverse problems,” in Proc. ECCV’08, ser. Lect. Notes Comp. Sci., D. A. Forsyth, P. H. S. Torr, and A. Zisserman, Eds. Springer, 2008, pp. 57–68.
-  R. Wagner, V. Delouille, and R. Baraniuk, “Distributed wavelet de-noising for sensor networks,” in Proc. IEEE Int. Conf. Dec. and Contr., San Diego, CA, Dec. 2006, pp. 373–379.
-  S. Barbarossa, G. Scutari, and T. Battisti, “Distributed signal subspace projection algorithms with maximum convergence rate for sensor networks with topological constraints,” in Proc. IEEE Int. Conf. Acc., Speech, and Signal Process., Taipei, Apr. 2009, pp. 2893–2896.
-  C. Guestrin, P. Bodik, R. Thibaux, M. Paskin, and S. Madden, “Distributed regression: an efficient framework for modeling sensor network data,” in Proc. Int. Symp. Inf. Process. Sensor Netw., Berkeley, CA, Apr. 2004, pp. 1–10.
-  J. B. Predd, S. R. Kulkarni, and H. V. Poor, “A collaborative training algorithm for distributed learning,” IEEE Trans. Inf. Theory, vol. 55, no. 4, pp. 1856–1871, Apr. 2009.
-  G. Mateos, J.-A. Bazerque, and G. B. Giannakis, “Distributed sparse linear regression,” IEEE Trans. Signal Process., vol. 58, no. 10, pp. 5262–5276, Oct. 2010.
-  D. K. Hammond, P. Vandergheynst, and R. Gribonval, “Wavelets on graphs via spectral graph theory,” Appl. Comput. Harmon. Anal., vol. 30, no. 2, pp. 129–150, Mar. 2011.
-  F. K. Chung, Spectral Graph Theory. Vol. 92 of the CBMS Regional Conference Series in Mathematics, AMS Bokstore, 1997.
-  T. Bıyıkoğlu, J. Leydold, and P. F. Stadler, Laplacian Eigenvectors of Graphs. Lecture Notes in Mathematics, vol. 1915, Springer, 2007.
-  J. C. Mason and D. C. Handscomb, Chebyshev Polynomials. Chapman and Hall, 2003.
-  G. M. Phillips, Interpolation and Approximation by Polynomials. CMS Books in Mathematics, Springer-Verlag, 2003.
-  T. J. Rivlin, Chebyshev Polynomials. Wiley-Interscience, 1990.
-  W. N. Anderson and T. D. Morley, “Eigenvalues of the Laplacian of a graph,” Linear Multilinear Algebra, vol. 18, no. 2, pp. 141–145, 1985.
-  R. Tibshirani, “Regression shrinkage and selection via the Lasso,” J. Royal. Statist. Soc. B, vol. 58, no. 1, pp. 267–288, 1996.
-  S. Chen, D. Donoho, and M. Saunders, “Atomic decomposition by basis pursuit,” SIAM J. Sci. Comp., vol. 20, no. 1, pp. 33–61, Aug. 1998.
-  I. Daubechies, M. Defrise, and C. De Mol, “An iterative thresholding algorithm for linear inverse problems with a sparsity constraint,” Commun. Pure Appl. Math., vol. 57, no. 11, pp. 1413–1457, Nov. 2004.
-  P. L. Combettes and V. R. Wajs, “Signal recovery by proximal forward-backward splitting,” Multiscale Model. Sim., vol. 4, no. 4, pp. 1168–1200, Nov. 2005.
-  D. P. Bertsekas and J. N. Tsitsiklis, Parallel and Distributed Computation: Numerical Methods. Prentice-Hall, 1989.