An Example-Driven Introduction to Data Analytics on Graphs
Ljubisa Stankovic Danilo Mandic Milos Dakovic Ilya Kisil Ervin Sejdic Anthony G. Constantinides










Understanding the Basis of Graph Signal Processing via an Intuitive Example-Driven Approach
    Ljubiša Stanković, 
Danilo Mandic,
 Miloš Daković,
Ilya Kisil,
Ervin Sejdić,
Anthony G. Constantinides
    December 30, 2023
==========================================================================================================










§ SCOPE


Graphs are irregular structures which naturally account for data integrity, however, traditional approaches have been established outside Signal Processing, and largely focus on analyzing the underlying graphs rather than signals on graphs.
 Given the rapidly increasing availability of multisensor and multinode measurements, likely recorded on irregular or ad-hoc grids, it would be extremely advantageous to analyze such structured data as graph signals and thus benefit from the ability of graphs to incorporate spatial awareness of the sensing locations, physical intuition and sensor importance, and the ease of local versus global sensor association.  The aim of this Lecture Note is therefore to establish a common language between graph signals, defined on irregular signal domains, and some of the most fundamental paradigms in DSP, such as spectral analysis of multichannel signals, system transfer function, digital filter design, parameter estimation, and optimal denoising.  

This is achieved through a physically meaningful and intuitive real-world example of geographically distributed multisensor temperature estimation.
A similar spatial multisensor arrangement is already widely used in Signal Processing curricula to introduce minimum variance estimators and Kalman filters, and by adopting this framework we facilitate a seamless integration of graph theory into the curriculum of existing DSP courses. By bridging the gap between standard approaches and graph signal processing, we also show that standard methods can be thought of as special cases of their graph counterparts, evaluated on line graphs. It is hoped that our approach would not only help to demystify graph theoretic approaches in education but it would also empower practitioners  and researchers to explore a whole host of otherwise prohibitive modern applications. 







§ RELEVANCE


 In classical Signal Processing, the signal domain is determined by equidistant time instants or by a set of spatial sensing points on a uniform grid. However, increasingly the actual data sensing domain may not even be related to the physical dimensions of time and/or  space, and it typically does exhibit various forms of regularity. For example, in social or web-related networks, the sensing points and their connectivity pertain to specific objects/nodes and topology of their links. It should be noted that  even for the data acquired in well defined time and space domains,  the introduction of new relations between the signal samples, through graphs, may yield new insights into the analysis and provide enhanced data processing (e.g., based on local similarity neighborhoods). The advantage of graphs over classical data domains is that graphs account naturally for irregular data relations in the problem definition, together with the corresponding data connectivity in the analysis. 



 
  

Indeed, Graph Signal and Information Processing is particularly well suited to making sense from data acquired over irregular data domains, which can be achieved, for example, by leveraging intuitions developed on Euclidean domains,  by employing analogies with other irregular domains such as polygon meshes and manifolds, or learning the mutual connectivity from available sets of data. In many emerging applications, e.g., Big Data, this also introduces a number of new challenges:

	
  * Basic concepts must be revisited in order to accommodate structured but often incomplete information, 
	
  * New physically meaningful frameworks, specifically tailored for heterogeneous data sources, are required, and
	
  * Trade-offs between performance and numerical requirements are a prerequisite when operating in real-time.


The common language and enhanced intuition between the  graph  approaches and their standard counterparts, illuminated in this article through  the relationships between the vertex and time domains,  may be naturally generalized to address the above challenges and spur further developments in the curricula on Statistical Signal Processing, Graph Signal Processing, and Big Data.



§ PREREQUISITES

This Lecture Note assumes a basic knowledge of Linear Algebra and Digital Signal Processing.  
 


§ HISTORY OF GRAPH THEORETIC APPLICATION

Graph theory, as a branch of mathematics, has existed for almost three centuries. 
The beginning of graph theory applications in electrical engineering dates back to the mid-XIX century and the definition of Kirchoff's laws. 
Owing to their inherent “spatial awareness", graph models have since become a de facto standard for data analysis across the science and engineering areas, including chemistry,  operational research, social networks, and computer sciences. 

A systematic account of graph theory as an optimization tool can be attributed to the seminal book by Nicos Christofides of Imperial College London, published in 1975 <cit.>. Soon after gaining prominence in general optimization, it was very natural to explore the application of graph theory in signal processing and related areas <cit.>. Indeed, perhaps the first lecture course to teach graph theory to then emerging communication networks and channel coding student cohort was introduced by the author Anthony Constantinides in 1970s. This helped to establish and formalize the connections between general optimization and the topology of a communication network, and has spurred further applications in image processing <cit.>. 


After a relative lull over the  next two decades, current developments in graph theory owe their prominence to the emergence of modern data sources, such as large-scale sensor and social networks, which inherently provide rich underlying physical, social, and geographic structures that require new ways to  establish statistical inference, leading to data processing on graphs, within a new fast maturing field of Graph Signal Processing, <cit.>. 






§ PROBLEM STATEMENT: AN ILLUSTRATIVE EXAMPLE

 	Graphs and Graph Signal Processing represent quite a general mathematical formalism which, albeit different from classic concepts, does admit the development of graph-domain counterparts of well established DSP paradigms. It would therefore be valuable to introduce such a general concept in an inductive and intuitive way, through a  simple, general enough and well-understood example regarding a commonly considered topic in classical DSP. 
 	
 	To this end consider a multi-sensor setup, shown in Fig. <ref>, for measuring a temperature field in a known geographical region; such a set-up is typically used in the context of minimum-variance estimators and Kalman filters. The temperature sensing locations are chosen according to the significance of a particular geographic area to local users, with N=64 sensing points in total, as shown in Fig. <ref>a). The temperature field is denoted by {x(n)} and a snapshot of its values is given in Fig. <ref>b). 
 Each measured sensor signal can then be mathematically expressed as 
    x(n)=s(n)+ε(n),

  where s(n) is the true temperature  that would have be obtained in  ideal measuring conditions and ε(n) comprises the adverse effects of the local environment on sensor readings or faulty sensor activity, and is referred to as “noise" in the sequel. For illustrative purposes, in our study ε(n) was modeled as a realization of white, zero-mean, Gaussian process, with standard deviation σ_ε=4. It was added to the signal, s(n), to yield the signal-to-noise ratio in {x(n)} of SNR_0=14.2 dB.  






Remark 1: Classical signal processing requires an arrangement of the quintessentially spatial temperature samples in Fig. <ref>a) into a line structure shown in Fig. <ref>b). Obviously, such “lexicographic" ordering is not amenable to exploiting the spatial information related to the actual sensor arrangement, dictated by the terrain. For example, this renders classical analyses of this temperature field inapplicable (or at best suboptimal), as the performance critically depends on the chosen sensor ordering scheme. This exemplifies that even a most routine temperature measurement  setup requires a more complex estimation structure than the simple  line one corresponding to the classical signal processing framework, shown in Fig. <ref>b).  

To introduce a “situation-aware" noise reduction scheme for the temperature field in Fig. <ref>, we proceed to explore  a graph-theoretic framework to this problem, starting from a local signal average operator. In classical Signal Processing this can be achieved through a moving average operator, e.g., through averaging across the neighboring data samples, or equivalently neighboring nodes, as in the  line graph in Fig. <ref>b), and  for each sensing point. Physically, such local neighborhood should indeed include close neighboring sensing points but which also exhibit similar meteorological properties defined by the distance, altitude difference, and other terrain properties. In other words, since the sensor network in Fig. <ref> measures a set of related temperatures from irregularly spaced sensors, an effective estimation strategy should include domain knowledge – not possible to achieve with standard DSP (line graph).   

Consider the local neighborhoods for the sensing points n=20, 29, 37, and 41, shown in Fig. <ref>a). The cumulative temperature for each sensing point is then given by

    y(n)=∑_m  at  and around  n   x(m),
 
so that the local average temperature for a sensing point n may be easily obtained by dividing the cumulative temperature, y(n), with the number of included sensing points.
For example, for the sensing points n=20 and n=37, presented in Fig. <ref>a), the “domain knowledge aware" local estimation takes the form

    y(20)     =x(20)+x(19)+x(22)+x(23) 
    
    y(37)     = x(37)+x(32)+x(33)+x(35)+x(61).

For convenience, the full set of relations among the sensing points can now be arranged into the matrix form, to give

    𝐲=𝐱+𝐀𝐱,
 
where the matrix 𝐀 indicates the connectivity structure of the neighboring sensing locations that should be involved in the calculation for each y(n). 
The matrix  𝐀 is therefore referred to as the connectivity or adjacency matrix of a graph. Its elements are either 1 (if the corresponding vertices are related) or 0 (if they are not related). Fig. <ref>b) shows the sensing locations with the corresponding connectivity for the temperature estimation scenario in Fig. <ref>a). From (<ref>) we can observe, for example, that the 20th row of the adjacency matrix 𝐀  will have all zero elements, except for A_20,19=1, A_20,22=1, and A_20,23=1 (for more detail see the electronic supplement). 




This simple real-world example can be interpreted within the graph signal processing framework as follows: 

	
  * The sensing points where the signal is measured are designated as the graph vertices, see Fig. <ref>,
	
  * The vertex-to-vertex lines indicating the connectivity among the sensing points are called the graph edges, 
	
  * The vertices and edges form a graph, as in Fig. <ref>b), a new and very structurally rich signal domain, 
	
  * The graph, rather than a standard vector of sensing points, is then used for analyzing and processing data, as it is equipped with spatial and physical awareness,  
	
  * The measured temperatures are now interpreted as signal samples on graph, as shown in Fig. <ref>, 	
	
  * Similar to traditional signal processing, this new graph signal may have many realizations on the same graph and may include noise,
	
  * Through relation (<ref>), we have therefore introduced a simple graph system for physically and spatially aware signal averaging (a linear first-order graph system).




To emphasize our trust in a particular sensor and to model mutual sensor  relevance,  a weighting scheme may be imposed on the edges (connectivity) between the sensing points, in the form

    y(n)=x(n)+∑_m n W_nmx(m).

The weight W_nm indicates the strength of the coupling between signal values at the sensing points n and m; its value is zero if the points n and m are not related and for n=m. We have now arrived at a weighted graph, whereby each edge has an associated weight, W_nm, which adds a “mutual sensor relevance" information to the already established “spatial awareness" modeled by the edges.  In our example, a matrix form of a weighted cumulative graph signal now becomes 

    𝐲=𝐱+𝐖𝐱.

This equips graph signal models with additional flexibility. In order to produce unbiased estimates, instead of the cumulative sums in (<ref>) and (<ref>), the weighting coefficients within the estimate for each y(n) should sum up to unity. This may be achieved through a  normalized form of (<ref>), given by

    𝐲=1/2(𝐱+𝐃^-1𝐖𝐱),

where the elements of the diagonal normalization matrix, 𝐃, called the degree matrix, are D_nn=∑_m W_nm  and 𝐃^-1𝐖 is referred to as a random walk weight matrix.  
When this simple normalized first-order system is employed to filter the original noisy signal from Fig. <ref>, an improvement of  6 dB over the original signal-to-noise ratio, SNR_0=14.2 dB, is achieved. 


Another important operator for graph signal processing is the graph Laplacian, 𝐋, which is defined as

    𝐋=𝐃-𝐖.




Remark 2: A graph is fully specified  by the set of its vertices and their connectivity scheme (designated by edges). The edges may be defined by the adjacency matrix, 𝐀, with A_mn∈{0,1}, for unweighted graphs or by the “connectivity strength" weighting matrix, 𝐖, with W_mn∈ℝ^+, for weighted graphs. The degree matrix, 𝐃, and the Laplacian matrix,  𝐋, with L_mn∈ℝ,  are defined using the adjacency/weighting matrix. When the relations between all pairs of vertices are mutually symmetric, then all the matrices involved are also symmetric, and such graphs are called undirected.  If that is not the case, then the adjacency/weighting matrix is not symmetric and such graphs are called directed graphs.     




The above-introduced graph framework is quite general and admits application to many different scenarios. For example, when performing an opinion poll within a social network, the members of that social network are treated as the vertices (data acquisition points). Their friendship relations are represented by the edges which model graph connectivity while the member answers play the role of  graph signal values. 

The definition of an appropriate graph structure is a prerequisite for physically meaningful and computationally efficient
graph signal processing applications. Three important classes of problems, regarding the way how the graph topology is defined, are described in Sidebar 1.

In the following, we shall demonstrate how this simple and intuitive  concept provides a natural and straightforward platform to introduce the graph-counterparts of several fundamental signal processing algorithms.





§ SYSTEM ON A GRAPH


The signal shift operator (unit time delay) is the lynchpin in discrete-time signal processing, but it is not so obvious to define on graphs due to the rich underlying connectivity structure. Topologically, the signal shift on a graph can be viewed as the movement of a signal sample from the considered vertex along all edges connected to this vertex. The signal (backward) shift operator can then be compactly defined using the graph adjacency matrix as 𝐱_shifted=𝐀𝐱.

To draw distinction between the standard shift and the graph shift operator, consider the line graph in Fig. <ref>b) (bottom) and the “spatial aware" graph in Fig. <ref>a), b), and assume that the input signal is a pulse that occurs only at the sensor n=29, that is, x(n)=δ (n-29). The shifted signal in classic signal processing (line graph in the bottom of Fig. <ref>b) (bottom)) will be  x_shifted(n)=δ(n-28) and can be considered as a movement of the delta pulse along the line graph from vertex n to vertex (n-1). The same principle can be applied to the graph domain in Fig. <ref>a) whereby the delta pulse from vertex n=29 is moved to all its connected vertices, to obtain the shifted graph signal,  x_shifted(n)=δ(n-27)+δ(n-28)+δ(n-51)+δ(n-59), as shown in Fig. <ref>. 




 If the shifted signal values are also scaled by the weighting coefficients of the corresponding edges, then the shifted signal is given by 𝐖𝐱.   Since the Laplacian can also be used as a shift operator, we will adopt the symbol 𝐒 to denote a general shift operator on a graph, which yields a graph shifted signal 𝐒𝐱.


Remark 3: The standard shift operator, x(n)=x(n-1), is a “one-to-one" mapping, while the graph shift operator, 𝐱_shifted=𝐒𝐱, is a “one-to-many" mapping which accounts for the underlying physics of the sensing process (in our example), not possible to achieve with standard DSP. Moreover, it also allows us to incorporate a contextual relation between the vertices within the irregular grid trough the weighting matrix  𝐖.  Notice that the graph shift operator does not satisfy the isometry property since the energy of the shifted signal is not the same as the energy of the original signal. 


In analogy to the pivotal role of time shift in standard system theory, a system on a graph can be implemented as a linear combination of a graph signal and its graph shifted versions.  The notion of a system is used in its classical sense, as a set of physical rules (an algorithm) that transforms an input graph signal into another (output) graph signal. The output graph signal from a system on a graph can then be written as

    𝐲=h_0 𝐒^0 𝐱+h_1 𝐒^1 𝐱+…+h_M-1𝐒^M-1 𝐱=∑_m=0^M-1h_m 𝐒^m 𝐱,

where,  by definition 𝐒^0=𝐈, while h_0, h_1,  …, h_M-1 are the system coefficients to be found (see Section <ref>). Notice that for the directed and unweighted line graph in Fig. <ref>b) (bottom), the system on a graph in (<ref>) reduces to the well known standard Finite Impulse Response (FIR) filter, given by

    y(n)=h_0x(n)+h_1x(n-1)+⋯+h_M-1x(n-M+1).


Remark 4: The above established link between the classical transfer function of a physical system and its graph-theoretic counterpart may serve to promote new algorithmic approaches, which stem from signal processing, into many application scenarios that are directly considered as graphs. 

Observe that the Laplacian operator applied on a signal,  𝐋𝐱, can be considered as a combination of the scaled original signal, 𝐃𝐱, and its weighted shifted version, 𝐖𝐱, since  𝐋𝐱=𝐃𝐱-𝐖𝐱.
A system defined using the graph Laplacian is obtained from (<ref>) by replacing 𝐒=𝐋, and has the form

    𝐲=𝐋^0 𝐱+h_1 𝐋^1 𝐱+…+h_M-1𝐋^M-1 𝐱

therefore allows us to always produce an unbiased estimate of a constant c, that is, if 𝐱=𝐜 then 𝐲=𝐜, since 𝐋𝐜=0.   


From (<ref>), a simple first order system based on the graph Laplacian  can be written as

    𝐲=𝐱 + h_1 𝐋𝐱

and is amenable, with slight modifications, to being used for efficient low-pass graph filtering, see Sidebar 2.  




Remark 5:  A system on a graph is conveniently defined by  the “graph transfer function",  H(𝐒), as   

    𝐲=H(𝐒)𝐱.


For an unweighted graph, the adjacency matrix, 𝐀, is commonly used as a shift matrix, 𝐒, while the Laplacian matrix, 𝐋=𝐃-𝐖, is used to define a shift  on a weighted graph.


Properties of a system on a graph: Following the above discussion, it is now possible to link the properties of linear systems with those of systems on a graph. From equations (<ref>)-(<ref>) the system on a graph is said to be:


  * Linear, if 

    H(𝐒)(a_1𝐱_1+a_2𝐱_2)=a_1𝐲_1+a_2𝐲_2.
 

  * Shift invariant, if

    H(𝐒)(𝐒𝐱)=𝐒(H(𝐒)𝐱).



Remark 6: A system on a graph, defined by 

    H(𝐒)=h_0 𝐒^0+h_1 𝐒^1+…+h_S-1𝐒^M-1

is linear and shift invariant, since the matrix multiplication of the square weighting matrices is associative ( 𝐒(𝐒𝐒)=(𝐒𝐒)𝐒), that is  𝐒𝐒^m=𝐒^m𝐒.  






§ GRAPH FOURIER TRANSFORM


While classic spectral analysis is performed in the Fourier domain,  spectral representations of graph signals employ either the adjacency/weighting matrix or the graph Laplacian eigenvalue decomposition. For the latter case we have

    𝐋=𝐔Λ𝐔^-1,

where 𝐔 is an orthonormal matrix of the eigenvectors, 𝐮_k, of the graph Laplacian matrix, 𝐋, (in its columns), and Λ is a diagonal matrix of the corresponding eigenvalues, λ_k. These eigenvectors may then be used for the spectral-based clustering of graph vertices, see Sidebar 3.




The graph Fourier transform, 𝐗, of a graph signal, 𝐱, is then defined as

    𝐗=𝐔^-1𝐱.


Physically, since 𝐔^-1=𝐔^T,  the element X(k) of  a graph Fourier transform, 𝐗, represents a projection of the graph signal, 𝐱, onto the k-th eigenvector, 𝐮_k∈𝐔, that is

    X(k)=∑_n=1^Nx(n)u_k(n).



 
The  inverse graph Fourier transform is then straightforwardly obtained as

    𝐱=𝐔 𝐗

or

    x(n)=∑_k=1^NX(k)u_k(n).


Remark 7: In analogy to the  classic Fourier transform where the signal is projected onto a set of harmonic orthogonal bases, 𝐗=𝐔^-1𝐱, where 𝐔 is the matrix of harmonic bases 𝐮_k=[1, e^j2π k/N, …, e^jπ (N-1)k/N]^T/√(N),  the graph  Fourier transform can be understood as a signal decomposition onto the set of eigenvectors of the graph Laplacian (or the adjacency matrix) that serve as orthonormal basis functions. In the case of a circular graph, the graph Fourier transform reduces to the standard discrete Fourier transform (DFT). For this reason, the transform in (<ref>) is referred to as the Graph  Fourier transform (GFT).

Classic spectral analysis can thus be considered as a special case of graph signal spectral analysis, with the adjacency matrix defined on an unweighted   circular directed graph (a line graph with the  connected last and first vertex), when 𝐮_k=[1, e^j2π k/N, …, e^jπ (N-1)k/N]^T/√(N). This becomes obvious by recognizing that the eigenvalues of a directed unweighted circular graph,  λ_k=e^-j2π k /N, are easily obtained as a solution of the eigenvalue/eigenvector (EVD) relation 𝐀𝐮_k=λ_k𝐮_k. For a vertex n, this relation is of the form u_k(n-1)=λ_k u_k(n). The previous vector elements u_k(n) and eigenvalues λ_k are the solutions of this difference equation.  It can be shown that the eigenvectors of the graph Laplacian of a line graph are real-valued harmonic functions, whose combinations can produce the standard complex-valued DFT basis functions, in an indirect way. The standard signal representation in Fig <ref>b) therefore corresponds to a signal whose domain is a line graph.


As is common in signal processing, the true temperature was simulated through a linear combination of several graph Laplacian eigenvectors (serving as basis functions) in the form 𝐱=-160𝐮_1+16𝐮_2-8𝐮_3-40𝐮_4+ 16𝐮_5-24𝐮_6+ε(n), where the random Gaussian noise, ε(n), had standard deviation  σ_ε=4.



§ SPECTRAL DOMAIN OF A SYSTEM ON GRAPHS

Consider a system on a graph, as in (<ref>), defined by its Laplacian matrix, given by

    𝐲=∑_m=0^M-1h_m 𝐋^m 𝐱.

Upon employing the eigen-domain (graph spectral) representation of the Laplacian matrix, 𝐋=𝐔Λ𝐔^-1, we have

    𝐲=∑_m=0^M-1h_m 𝐔Λ^m 𝐔^-1 𝐱
    =
    𝐔  H(Λ)𝐔^-1 𝐱,

where

    H(Λ)=
    ∑_m=0^M-1h_m Λ^m

is the transfer function of the graph system.

From (<ref>),
𝐔^-1𝐲=
 H(Λ)𝐔^-1 𝐱,
or in terms of the graph Fourier transform of the input and output signal
 
    𝐘=
     H(Λ) 𝐗.



The classic spectral transfer function for (<ref>) is then obtained by using the adjacency matrix of an unweighed directed circular graph whose eigenvalues are  λ_k=e^-j2π k /N. 




§ SPECTRAL DOMAIN FILTER DESIGN


Consider a desired  graph transfer function, G(Λ). Like in classic signal processing, a system with this transfer function can be implemented either in the spectral domain or in the vertex domain. 

The spectral domain implementation is straightforward and can be performed in the following three steps: 


  * Calculate the  GFT of the input graph signal 𝐗=𝐔^-1𝐱, 

  * Multiply the  GFT of the input graph signal with transfer function  G(Λ) to obtain  𝐘=G(Λ)𝐗, and 

  * Calculate the output graph signal as the
inverse graph Fourier transform of 𝐘 to yield  𝐲=𝐔𝐘. 


Notice that this procedure may be computationally very demanding for large graphs where it may be easier to implement the desired filter (or its close approximation) in the vertex domain, in analogy to the time domain in the classical approach. This means that we have to find the coefficients, h_0,h_1,…,h_M-1 in (<ref>), such that its spectral representation, H(Λ),  is equal (or at least as close as possible) to the desired G(Λ). 




In other words, the transfer function of the vertex domain system in (<ref>), given by
H(λ_k)=h_0 +h_1 λ_k^1+… h_M-1λ_k^M-1, should be equal to the desired transfer function, G(λ_k), for each spectral index, k. This condition leads to a system of linear equations

    h_0 +h_1 λ_1^1+… h_M-1λ_1^M-1=G(λ_1) 
    
    	h_0 +h_1 λ_2^1+… h_M-1λ_2^M-1=G(λ_2) 
    ⋮
    
    	h_0 +h_1 λ_N^1+… h_M-1λ_N^M-1=G(λ_N).

The matrix form of this system is given by 

    𝐕_λ 𝐡=𝐠,

where 𝐕_λ is a Vandermonde matrix formed of the eigenvalues, λ_k, while
𝐡=[h_0, h_1,…,h_M-1]^T
is the vector of system coefficients that we wish to estimate, and   

    𝐠=[G(λ_1), G(λ_2),…,G(λ_N)]^T=diag(G(Λ)).


The system order M is typically significantly lower than the number of equations, N, in (<ref>).
For such an overdetermined case, the least-squares approximation  of 𝐡 is obtained by minimizing the squared error, 
e^2=‖𝐕_λ𝐡-𝐠‖ _2 ^2.
Like in standard least-squares, the solution is obtained by a direct minimization, ∂ e ^2 /∂𝐡^T=0, to yield 

    ĥ=(𝐕_λ^T 𝐕_λ)^-1𝐕_λ^T 𝐠=pinv(𝐕_λ)𝐠.
 
The so obtained solution, ĥ, therefore represents the mean square error minimizer for 𝐕_λ𝐡=𝐠. Notice that this solution may not satisfy 𝐕_λ𝐡=𝐠, in which case the coefficients ĝ (its spectrum Ĝ(Λ)) may be used, that is

    𝐕_λĥ=ĝ.

Such a  solution, in general, differs from the desired system coefficients 𝐠 (its spectrum G(Λ)).


Example: Consider the graph signal from Fig. <ref>. The task is to design a graph filter whose frequency response is g(λ_k)=exp(-λ_k) and to then filter the graph signal using this spectral domain graph filter. For M=4,  the corresponding system coefficients can be found to be
h_0=0.9606,
h_1=-0.7453,
h_2=0.1936, and
h_3=-0.0162. Upon signal filtering using the so defined graph transfer function,  the output signal-to-noise ratio was SNR=21.74  dB, that is a 7.54 dB improvement over the original  signal-to-noise ratio  SNR_0=14.2  dB.    

More detail on the solution of the system in (<ref>) and (<ref>) is provided in Sidebar 4. 





§ OPTIMAL DENOISING


Consider a measurement, as in the temperature measurement scenario in Fig. <ref>, which is composed of a slow-varying desired signal, 𝐬, and a superimposed fast changing disturbance, ε, to give

    𝐱=𝐬+ε.

The aim is to design a graph filter for disturbance suppression (denoising), the output of which is denoted by 𝐲, <cit.>.

The optimal denoising task can then be defined through a minimization of the cost function

    J=1/2‖𝐲-𝐱‖_2^2 + α𝐲^T 𝐋𝐲.

The minimization of the first term, 1/2‖𝐲-𝐱‖_2^2, enforces the output signal, 𝐲, to be as close as possible, in terms of the minimum residual disturbance power, to the available observations, 𝐱. As mentioned before, the second term, 𝐲^T 𝐋𝐲, represents a measure of smoothness of the graph filter output,  𝐲. For more detail on promoting smoothness of a graph signal, see Sidebar 2. The parameter α models a balance between the closeness of the output, 𝐲, to the observed data, 𝐱,  and the smoothness of output estimate 𝐲. While the problem in (<ref>) could be expressed through a constrained Lagrangian  optimization, whereby we choose to focus more on the graph theoretic issues and hence we adopt a simpler option whereby the mixing parameter α is chosen empirically.

The solution to this minimization problem follows from

    ∂ J/∂𝐲^T = 𝐲-𝐱 +  2α𝐋𝐲 =0

and results in a smoothing optimal denoiser in the form

    𝐲= (𝐈 + 2α𝐋)^-1𝐱.

The Laplacian spectral domain form of this relation is

    𝐘= (𝐈 + 2αΛ)^-1𝐗,

with the corresponding graph filter transfer function

    H(λ_k)=1/1+2αλ_k.

For a small α, H(λ_k) ≈ 1 and 𝐲≈𝐱, while for a large α, H(λ_k) ≈δ(k) and 𝐲≈ const., which enforces 𝐲 to be maximally smooth (a constant, without any variation). Using α=4, the obtained output signal-to-noise ratio for the graph signal from  Fig. <ref>  was SNR=26 dB, a 11.8 dB improvement over the original SNR_0=14.2 dB. 


Remark 8: There are many cases when the graph topology is unknown, so that the graph structure, i.e.,
the Laplacian (graph edges and their weights) is also unknown. 
To this end, we may employ a class of methods for graph topology learning, based on the minimization of the cost function in (<ref>) with respect to both the Laplacian, 𝐋, and the output signal, 𝐲, with additional (commonly sparsity) constraints imposed on the Laplacian values.





§ CURRENT GRAPH SIGNAL PROCESSING CHALLENGES



Current research is mainly focused on graphs themselves, like for example,  on reducing the complexity of calculation in very large graphs, including downsampling, multirate analysis, compressive sensing, graph segmentation,  non-linear GSP, robust GSP, deep learning architectures for graph signals, multidimensional graph signals, and vertex-varying and vertex-frequency analysis. 





§ WHAT WE HAVE LEARNED


Natural signals (speech, biomedical, video) reside over irregular domains and are, unlike the signals in communications, not adequately processed using, e.g., standard harmonic analyses. 
While Data Analytics are heavily dependent on advances in DSP, neither the EE graduates worldwide nor practical data analysts are yet best prepared to employ graph algorithms in their future jobs. Our aim has been to fill this void by providing an example-driven platform to introduce graphs and their properties through the well understood notions of transfer functions, Fourier transform, and digital filtering. 

While both a graph with N vertices and a classical discrete time signal with N samples can be viewed as N-dimensional vectors, structured graphs are much richer irregular domains which convey information about both the signal generation and propagation mechanisms. This allows us to employ intuition and our know-how from Euclidean domains to revisit basic dimensionality reduction operations, such as coarse graining of graphs (cf. standard downsampling). In addition, in the vertex domain a number of different distances (shortest-path, resistance, diffusion) have useful properties which can be employed to maintain data integrity throughout the processing, storage, communication and analysis stages, as the connectivities and edge weights are either dictated by the physics of the problem at hand or are inferred from the data. This particularly facilitates maintaining control and intuition  over distributed operations throughout the processing chain. 


It is our hope that this lecture note has helped to
demystify graph signal processing for students and educators, together with empowering practitioners with enhanced intuition in graph-theoretic  design and optimization. This material may also serve as a vehicle to seamlessly merge curricula in Electrical Engineering and Computing. 
The generic and physically meaningful nature of this example-driven Lecture Note is also likely to promote intellectual curiosity and serve as a platform to explore the numerous opportunities in manifold applications in our ever-growing interconnected world, facilitated by the Internet of Things.



§ ACKNOWLEDGMENTS
 We are privileged to have had the help and advice of one of the pioneers in Graph Theory Professor Nicos Christofides. We are grateful for his time, his incisive comments and valuable advice. We would also like to express our  sincere gratitude to the students in our respective postgraduate courses, for their feedback on the material taught based on this Lecture Note.



§ AUTHORS


Ljubiša Stanković, FIEEE, (ljubisa@ac.me) is professor at the 
University of Montenegro. His research interests include time-frequency analysis, compressive sensing, and graph signal processing.  He is a vice-president of the National Academy of Sciences and Arts of Montenegro (CANU) and a member of the European Academy of Sciences and Arts. Prof. Stanković is a recipient of the 2017 EURASIP Best Journal Paper Award.


Danilo P. Mandic, FIEEE, (d.mandic@imperial.ac.uk) 
is a professor of signal processing at Imperial College London, United Kingdom. He is a member of the IEEE Signal Processing Society Education Technical Committee, and has received President’s Award for Excellence in Postgraduate Supervision at Imperial College. He is a recipient of the 2018 Best Paper Award in IEEE Signal Processing Magazine.

Miloš Daković (milos@ac.me) is professor at
the University of Montenegro. His research interests include
graph signal processing, and time-frequency analysis. 

Ilya Kisil (i.kisil15@imperial.ac.uk) is a Ph.D. candidate at Imperial College London. His research interests include tensor decompositions, big data, efficient software for large scale problems, and graph signal processing. 

Ervin Sejdić, SMIEEE, (esejdic@ieee.org) is an assistant
professor at the University of Pittsburgh, USA. His research interests
include biomedical signal processing, rehabilitation
engineering, and neuroscience. He received the USA
Presidential Early Career Award for Scientists and Engineers in
2016.

Anthony G. Constantinides, LFIEEE, (a.constantinides @imperial.­ac.uk) 
is emeritus professor of signal processing in the Department of Electrical and Electronic Engineering at Imperial College London, United Kingdom.


9                                                                              

NC N. Christofides, “Graph theory: An algorithmic approach", Academic Press, 1975.

TC F. Afrati and A. G. Constantinides, “The use of graph theory in binary block code construction", in Proceedings of the International Conference on Digital Signal Processing, pp. 228-233, Florence, Italy, 310 August - 2 September, 1978

AC O. J. Morris, M. de J. Lee, and A. G.  Constantinides, “Graph theory for image analysis: An approach based on the shortest spanning tree", IEE Proceedings F-Communications, Radar and Signal Processing, vol. 133, no. 2, pp. 146-152, 1986.

 
 JK S. Chen, R. Varma, A. Sandryhaila, and J. Kovačević, “Discrete
 signal processing on graphs: Sampling theory,"  IEEE Trans. on Signal
 	Processing, vol. 63, no. 24, pp. 6510-6523, Dec.15, 2015.
 
 Maura1A. Sandryhaila and J. M. F. Moura, “Discrete signal processing
 on graphs," IEEE Transactions on Signal Processing, vol. 61, no. 7, pp.
 1644–1656,
 Apr. 2013.
 
 Maura2A. Sandryhaila and J. M. F. Moura, “Discrete signal processing
 on graphs: Frequency analysis," IEEE Transactions on Signal Processing,
 vol.
 62, no. 12, pp. 3042–3054, Jun. 2014.
 
 ChuG. Cheung, E. Magli, Y. Tanaka, and M. K. Ng,  “Graph spectral image processing,” Proceedings of the IEEE, vol. 106(5), pp. 907-930, May 2018.
 
 shuman2013may
 A. Ortega, P. Frossard, J. Kovačević, J. M. F. Moura, P. Vandergheynst, “Graph signal processing: Overview, challenges, and applications”, Proceedings of the IEEE,  vol. 106(5), pp. 808–28, May 2018.
 
  DM S. Saito, H. Suzuki, and D. P. Mandic, “Hypergraph p-Laplacian: A differential geometry view”,  in Proceedings of the the Thirty Second AAAI Conference on Artificial Intelligence (AAAI-18), pp. 3984–3991, 2018.
 
 VF L. Stanković and E. Sejdić,
 “Vertex-frequency analysis of graph signals,”, Springer Nature, 2019.
 
 Opt S. Segarra, A. G. Marques,  and A. Ribeiro,  “Optimal graph-filter design and applications to distributed linear network operators”, IEEE Transactions on Signal Processing, 65(15), pp. 4117–4131, 2017.
 

 



