Sampling and Recovery of Graph Signals
P. Di Lorenzo S. Barbarossa P. Banelli
eess.SP





CHAPTER: SAMPLING AND RECOVERY OF GRAPH SIGNALS


[1]Paolo Di Lorenzo
2]Sergio Barbarossa,
1]Paolo Banelli

[1]University of Perugia, Department of Engineering, Via G. Duranti 93, 06125, Perugia, Italy.
[2]Sapienza University of Rome, Department of Information Engineering, Electronics, and Telecommunications, via Eudossiana 18, 00184, Rome, Italy.
*[3]Corresponding: paolo.dilorenzo@unipg.it



[Abstract]
The aim of this chapter is to give an overview of the recent advances related to sampling and recovery of signals defined over graphs. First, we illustrate the conditions for perfect recovery of bandlimited graph signals from samples collected over a selected set of vertexes. Then, we describe some sampling design criteria proposed in the literature to mitigate the effect of noise and model mismatching when performing graph signal recovery. Finally, we illustrate algorithms and optimal sampling strategies for adaptive recovery and tracking of dynamic graph signals, where both sampling set and signal values are allowed to vary with time. Numerical simulations carried out over both synthetic and real data illustrate the potential advantages of graph signal processing methods for sampling, interpolation, and tracking of signals observed over irregular domains such as, e.g., technological or biological networks.


[Keywords:]
Graph signal processing sampling on graphs interpolation on graphs adaptation and learning over networks





§ INTRODUCTION



In a large number of applications involving sensor, transportation, communication, social, or biological networks, the observed data can be modeled as signals defined over graphs, or graph signals for short. As a consequence, over the last few years, there was a surge of interest in developing novel analysis methods for graph signals, thus leading to the research field known as graph signal processing (GSP), see, e.g., <cit.>. The goal of GSP is to extend classical processing tools to the analysis of signals defined over an irregular discrete domain, represented by a graph, and one interesting aspect is that such methods typically come to depend on the graph topology, see, e.g., <cit.>.


A fundamental task in GSP is to infer the values of a graph signal by interpolating the samples collected from a known set of vertices. In the GSP literature, this learning task is known as interpolation from samples, and emerges whenever cost constraints limit the number of vertices that we can directly observe. This arises in several applications such as semi-supervised learning of categorical data <cit.>, environmental monitoring <cit.>, and missing value prediction as in matrix completion problems <cit.>.
Interpolation methods on graphs rely on the implicit assumption that nodes close to each other have similar values, i.e., the graph encodes similarity
among the values observed over the vertices. For instance, in an item-item graph in a recommendation system, a user would rate two similar items with similar ratings <cit.>. In the same way, predicting the functions of proteins based on a protein network relies on some notion of closeness among the nodes <cit.>. In other words, the signals of interest must be smooth functions over the graph. In GSP, the smoothness assumption is typically formalized in terms of (approximate) bandlimitedness over a graph Fourier basis, and enables recovery of the signal after sampling over a selected subset of vertices.

A first seminal contribution to sampling theory in GSP is given by <cit.>, where a sufficient condition for unique recovery is stated for a given sampling set; the approach was then extended in <cit.>. Most of the works on graph sampling theory assume that a portion of the graph Fourier basis is explicitly known. For example, the work in <cit.> provides conditions that guarantee unique reconstruction of signals spanned over a subset of vectors composing the graph Fourier basis, proposing also a greedy method to select the sampling set in order to minimize the effect of sample noise in the worst-case. Reference <cit.> exploited a smart partitioning of the graph in local-sets, and proposed iterative methods to reconstruct bandlimited graph signal from sampled data. The work in <cit.> creates a conceptual link between uncertainty principle and sampling of graph signals, and proposes several optimality criteria (e.g., the mean-square error) to select the sampling set in the presence of noise. Another valid approach is the so-called aggregation sampling <cit.>, which involves successively shifting a signal using the adjacency matrix and aggregating the values at a given node. Greedy sampling strategies with provable performance guarantees were proposed in <cit.> in a Bayesian reconstruction setting.

If the size of the graph signal is very large as, e.g., in web-scale graphs <cit.>, complexity becomes a crucial issue, such that in many cases we cannot assume to know or efficiently compute the graph Fourier basis. Some works have then proposed sampling methods that
do not require such previous knowledge. For instance, the work in <cit.> proposes efficient methods to select the sampling set based on powers of the variation operator to approximate the bandwidth of the graph signal. There are also alternative approaches that do not consider graph spectral information and rely only on vertex-domain characteristic, e.g., maximum graph cuts <cit.> and spanning trees <cit.>. Finally, there exist randomized sampling strategies, e.g., <cit.>. The work in <cit.> provides an efficient design of sampling probability distribution over the nodes, deriving bounds on the reconstruction error in the presence of noise and/or approximatively bandlimited signals. Reference <cit.> exploits compressive sampling arguments to derive random sampling strategies with variable density, thus also proposing a fast technique to estimate the optimal sampling distribution accurately. Last, the work in <cit.> proposes a sampling strategy tailored for large-scale data based on random walks on graphs.

The sampling strategies described so far involve batch methods for sampling and recovery of graph signals. In many applications such as, e.g., transportation networks, brain networks, or communication networks, the observed graph signals are typically time-varying. This requires the development of effective methods capable to learn and track dynamic graph signals from a carefully designed, possibly time-varying, sampling set. Some previous works have considered this specific learning task, see, e.g., <cit.>. Specifically, <cit.> proposed an LMS estimation strategy enabling adaptive learning and tracking from a limited number of smartly sampled observations. The LMS method in  <cit.> was then extended to the distributed setting in <cit.>. The work in <cit.> proposed a kernel-based reconstruction framework to accommodate time-evolving signals over possibly time-evolving topologies, leveraging spatio-temporal dynamics of the observed data. Finally, reference <cit.> proposes a distributed method for tracking
bandlimited graph signals, assuming perfect observations and a fixed sampling strategy.

In this chapter, we review some of the recent advances related to sampling and recovery of signals defined over graphs. Due to space limitations, such review will be limited only to some specific contributions.

The structure of the chapter is explained in the sequel. Sec. <ref> defines the adopted notation, and recalls some background on GSP. In Sec. <ref>, we illustrate the conditions for perfect recovery of bandlimited graph signals from samples collected according to design criteria proposed to mitigate the effect of noise or model mismatching. Finally, Sec. <ref> illustrates algorithms and optimal sampling strategies for adaptive recovery and tracking of dynamic graph signals, where both sampling set and signal values are allowed to vary with time.




§ NOTATION AND BACKGROUND



In this paragraph, we first introduce the notation that we will use throughout the chapter. Then, we briefly recall some basic notions from GSP that will be instrumental for the derivations and arguments of the following sections.


Notation. We  indicate  scalars  by  normal  letters (e.g., a); vector variables with bold lowercase letters (e.g., ) and matrix variables  with  bold  uppercase  letters (e.g., ). Scalars a_i and a_ij correspond to the i-th entry of  and the ij-th entry of , respectively. We indicate by _2 and _2 the  ℓ_2 norm  and  the  spectral  norm  of  the  vector  and  matrix , respectively. If  is rectangular, we denote by σ_i() the i-th singular value of ; if  is square, λ_i() represents the i-th eigenvalue of . The trace operator of matrix  is indicated with Tr(); diag() is a diagonal matrix having  as main diagonal; rank() denotes the rank of matrix ; () represents the determinant of , whereas pdet() is the pseudo-determinant of , i.e., the product of all non-zero eigenvalues of . The superscript ^H denotes the hermitian operator, i.e., the conjugate transposition of a vector or a matrix, whereas ^† denotes the pseudo-inverse of matrix . 𝔼{·} represents the expectation operator. A set of elements is denoted by a calligraphic letter (e.g., ), and || represents the cardinality of set , i.e., the number of elements of . The symbols ∪, ∩, and ∖ denote union, intersection, and difference among sets, respectively. Given a set , we denote its complement set as ^c, i.e., =∪^c and ∩ S^c=∅. 1 denotes the vector of all ones, whereas _ is the set indicator vector, whose i-th entry is equal to one, if  i ∈, or zero otherwise. 


Background. We consider a graph = (, ) consisting of a set of N nodes = {1,2,..., N}, along with a set of weighted edges ={a_ij}_i, j ∈, such that a_ij>0, if there is a link from node j to node i, or a_ij=0, otherwise. The adjacency matrix  of a graph is the collection of all the weights a_ij, i, j = 1, …, N. The combinatorial Laplacian matrix is defined as 𝐋 =  diag(1^T)-𝐀. A signal  over a graph  is defined as a mapping from the vertex set to the set of complex numbers, i.e., : →ℂ. The graph  is endowed with a graph-shift operator  defined as an N × N matrix whose entry (i,j), denoted with S_ij, can be non-zero only if i=j or the link (j,i)∈. The sparsity pattern of matrix  captures the local structure of ; common choices for  are the adjacency matrix <cit.>, the Laplacian <cit.>, and its generalizations <cit.>. We assume that  is diagonalizable, i.e., there exists an N × N matrix =[_1,…,_N] and an N × N diagonal matrix Λ that can be used to decompose  as =Λ^-1. When  is normal, i.e., when ^H=^H, matrix  is unitary and ^-1=^H.

Recovery of a signal from its sampled version is possible under the assumption that  admits a sparse representation. The basic idea when addressing the problem of sampling graph signals is to suppose that  plays a key role in explaining the signal of interest. More specifically, we assume that  can be expressed as a linear combination of a subset of the columns of , i.e.,

    =

where ∈ℂ^N is either exactly or approximately sparse. In this context, vectors {_i}_i=1^N are interpreted as the graph Fourier basis, and {s_i}_i=1^N
are the corresponding graph signal frequency coefficients, i.e.,

    = 𝐔^H

takes the role of the Graph Fourier Transform (GFT) of signal . As an example, in many cases the graph signal exhibits clustering features, i.e., it is a smooth function over each cluster (e.g., in semi-supervised learning <cit.>), but it may vary arbitrarily from one cluster to the other. In such a case, if the columns of  are chosen to represent clusters, the only nonzero (or approximately nonzero) entries of  are the ones associated to the clusters. In the case of undirected graphs,  may be composed from the eigenvectors of the Laplacian, which have well known clustering properties <cit.>.

















The localization properties of graph signals in vertex and frequency domains will play an important role in the ensuing arguments. To introduce such properties, we first define the matrix _∈ℂ^N×||, which represents the collection of all the columns of  associated with a subset of frequency indices ⊆{1,…,N}. Then, we introduce the N× N band-limiting operator

    _ = _^H_.

The role of _ is to project a vector  onto the subspace spanned by the columns of _. Thus, we say that a vector  is perfectly localized over the frequency set  (or -bandlimited) if

    _==__,

where _ is given in (<ref>), and the second equality comes from (<ref>) where we have exploited the sparsity of , which is different from zero (and equal to _∈ℂ^||) only in the frequency support .
Similarly, given a subset of vertices ⊆, we define the N× N vertex-limiting operator

    𝐃_ =  diag{_}.

Thus, we say that a vector  is perfectly localized over the subset ⊆ (or -vertex-limited) if


_=, with _ defined as in (<ref>).


We also denote by _ the set of all -bandlimited signals, and by _ the set of all -vertex-limited signals. The operators 𝐃_ and 𝐁_ are self-adjoint and idempotent, and represent orthogonal projectors onto the sets _ and _, respectively. Differently from continuous-time signals, a graph signal can be perfectly localized in both vertex and frequency domains. This property is formally stated in the following theorem <cit.>.


There is a graph signal  perfectly localized over both vertex set  and frequency set  (i.e. ∈_∩_) if and only if the operator ___ has an eigenvalue equal to one; in such a case,  is the eigenvector of ___ associated to the unit eigenvalue.

Perfect localization onto the sets  and  can be equivalently expressed in terms of the operator __ <cit.>, i.e., it holds if and only if

    ___2=____2  = 1.


























In the following section, we will illustrate the theory behind sampling and recovery of signals defined over graphs.



§ SAMPLING AND RECOVERY


Let us consider the observation of an -bandlimited graph signal over the sampling set . The observation model can be cast as:

    _=𝐏_^T =𝐏_^T𝐔__,

where _∈ℂ^|| is the observation vector over the vertex set , and 𝐏_∈ℝ^N× || is a sampling matrix whose columns are indicator functions for nodes in , and such that the orthogonal projector over _ is given by _=__^T [cf. (<ref>)]. The problem of recovering a bandlimited graph signal from its samples is then equivalent to the problem of properly selecting the sampling set , and then recover  from _ by inverting the system of equations in (<ref>). This approach is known as selection sampling and was addressed, for example, in  <cit.>, <cit.>, <cit.>, and <cit.>.


In the sequel, we will first consider the conditions for perfect recovery of bandlimited graph signals. Then, we will illustrate the effect of noise and model mismatching on the reconstruction performance. Also, since the identification of the sampling set  plays a key role in the conditions for signal recovery and in the reconstruction performance, we will illustrate optimization strategies to design the sampling set. Finally, we will illustrates results of numerical simulations carried out over synthetic and realistic data.



 §.§ Sampling and Perfect Recovery of bandlimited graph signals



We will now address the fundamental problem of assessing the conditions and the means for perfect recovery of  from _. To this aim, we introduce the operator _^c=-_, which projects onto the complement vertex set ^c = ∖. Starting from (<ref>), the necessary and sufficient conditions for perfect recovery are stated in the following Theorem <cit.>.


Any -bandlimited graph signal  can be perfectly recovered from its samples collected over the vertex set  if and only if

    _^c__2 < 1,

i.e., if there are no -bandlimited signals that are perfectly localized on ^c.



From (<ref>), a sufficient condition for signal recovery is the existence of the pseudo-inverse matrix =(_^T_)^†=(_^H__)^-1_^H_.
Since  _^H__=-_^H_^c_, we obtain that  exists if
_^H_^c__2=_^c__2<1, i.e., if (<ref>) holds true. Conversely, if _^c__2=1, there exist bandlimited signals that are perfectly localized over ^c [cf. (<ref>)]. Thus, if we sample one of such signals over , it would be impossible to recover  from those samples. This proves that condition (<ref>) is also necessary.




























Theorem <ref> and its proof also suggest the reconstruction formula:

    =_(_^T_)^†_=_(_^H__)^-1_^H__,

which guarantees reconstruction of the bandlimited graph signal  if condition (<ref>) holds true, and has computational complexity equal to O(N||^2). The above reconstruction formula is also known as consistent reconstruction <cit.> since it keeps the observed samples unchanged.

Let us consider now the implications of condition (<ref>) of Theorem <ref> on the sampling strategy. To fulfill (<ref>), we need to guarantee that there exist no signals that are perfectly localized over the vertex set ^c and the frequency set . Since, in general, we have

    _= _^T__+_^c^T __,

we need to guarantee that ___≠0̱ for any non-trivial vector _,












which requires _^T_ to be full column rank, i.e.,

    (_^T _) =  (_) = ||.

Of course, a necessary condition to satisfy (<ref>) is that

    || ≥ ||.

However, condition (<ref>) is not sufficient, because _^T_ may loose rank, depending on graph topology and samples' location. As a particular case, if the graph is not connected, the vertices can be labeled so that the Laplacian (adjacency) matrix can be written as a block diagonal matrix, with a number of blocks equal to the number of connected components. Correspondingly, each eigenvector of  can be expressed as a vector having all zero elements, except for the entries corresponding to the connected component. This implies that, if there are no samples over the vertices corresponding to the non-null entries of the eigenvectors with index included in , _^T_ looses rank.

More generally, even if the graph is connected, there may easily occur situations where matrix _^T_ is not rank-deficient, but it is ill-conditioned, depending on graph topology and samples' location. This case is particularly dangerous when the true signal is only approximately bandlimited (which is the case for most signals in practice) or when the samples are noisy. In such cases, not all sampling sets of given size are equally good, and it becomes fundamental to understand which is the best sampling set that achieves the smallest reconstruction error. Thus, in the sequel, we will first illustrate the effect of noise and model mismatching on graph signal reconstruction, and then we will describe sampling strategies satisfying several optimization criteria.



 §.§ The effect of noise and model mismatching


Let us consider first the reconstruction of bandlimited signals from noisy samples, where the observation model is given by:

    _ = _^T (  + ) = _^T __ + _^T,

where  is a zero-mean noise vector with covariance matrix _v=𝔼{^H}. To design an interpolator in the presence of noise, we consider the best linear unbiased estimator (BLUE), which is given by <cit.>:

    = _(_^H_(_^T_v_)^-1_^T_)^-1_^H_(_^T_v_)^-1_.

The estimator in (<ref>) minimizes the least square error and, if noise is Gaussian in (<ref>), it coincides with the minimum variance unbiased estimator, which attains the Cramér-Rao lower bound. It is immediate to see that (<ref>) is an unbiased estimator, i.e., 𝔼{}=. Furthermore, the mean square error (MSE) is given by <cit.>:

    MSE=𝔼-^2= Tr{(_^H_(_^T_v_)^-1_^T_)^-1}.

As a particular case, if noise is spatially uncorrelated, i.e., _v= diag{r_1^2,…,r_N^2}, and letting _,i^H be the i-th row of matrix _, we obtain:

    MSE= Tr{(_^H__v^-1_)^-1}= Tr{(∑_i∈_,i_,i^H/r_i^2)^-1}.

This illustrates how, in the presence of uncorrelated noise, the design of the sampling set should minimize the trace of the inverse of matrix _^H__v_.




















So far we assumed that the true signal  is perfectly bandlimited, i.e., ∈_. However, in most applications, the signals are only approximately bandlimited. In such a case, the recovery formula in (<ref>) applied to such signals leads to a reconstruction error, which is analyzed next. In general, an approximately bandlimited graph signal can be expressed as

    =_+Δ,

where _=_ is the bandlimited component, whereas Δ=_^c represents the model mismatch. Sampling the signal over the vertex set  and using (<ref>) as a recovery formula, then an upper bound (i.e., a worst-case) on the reconstruction error is given by <cit.>:

    -≤Δ/cos(θ_max),

where θ_max represents the maximum angle between the subspaces _ and _, which is defined as:

    cos(θ_max)=      =1inf  𝐃__2
         subject to  𝐁_=.

In particular, from (<ref>), it is easy to see that cos(θ_max)>0 if condition (<ref>) holds true. Intuitively, the bound in (<ref>) says that, for the worst-case error to be minimum, the sampling and reconstruction subspaces should be as aligned as possible. Therefore, for approximatively bandlimited signals, an optimal sampling set should be selected in order to maximize the smallest maximum angle between the subspaces _ and _. Interestingly, from (<ref>) and (<ref>), it appears clear that

    cos(θ_max)=σ_min(__)

Thus, in the presence of model mismatching, the design of the sampling set should maximize the minimum singular value of matrix __ or, equivalently, the minimum eigenvalue of matrix _^H__.

In the next section, we will illustrate the strategies used to optimize the selection of the sampling set.



 §.§ Sampling strategies



As previously mentioned, when sampling graph signals, besides choosing the right number of samples, whenever possible it is also fundamental to have a strategy indicating where to sample, as the samples' location plays a key role in the performance of reconstruction algorithms. In principle, in the ideal case (<ref>), any sampling set  that satisfies condition (<ref>) enables unique reconstruction through the interpolation formula in (<ref>). However, in the presence of noise or model mismatching, from (<ref>) and (<ref>)-(<ref>), it is clear that the quality of reconstruction is strongly affected by a careful design of the sampling set .  Different costs can then be defined to measure the reconstruction error and are based on optimal design of experiments <cit.>. For instance, if we seek for the optimal sampling set ^ opt of size M, as the set that minimizes the mean squared error in (<ref>), we have:

    ^ A-opt=||=Mmin   Tr{(_^H__v^-1_)^-1}.

This is analogous to the so-called A-optimal design <cit.>, and is equivalent to the one proposed in <cit.>. Similarly, if we aim to design the optimal sampling set of size M to minimize the worst-case reconstruction error in the presence of model mismatching [cf. (<ref>)-(<ref>)], we have:

    ^ E-opt=||=Mmax  σ_min(__),

which is equivalent to the so-called E-optimal design <cit.>. The above criterion is equivalent to the one proposed in <cit.> and, in general, it is useful to find a stable sampling set that satisfies condition (<ref>). To select the optimal sampling set, we should solve one of the problems in (<ref>) or (<ref>), which entail the selection of an M-element subset of 𝒱 that optimizes the adopted design criterion. This is a finite combinatorial optimization problem (which is known to be NP-hard <cit.>), whose solution in general requires an exhaustive search over all the possible combinations. Since the number of possible subsets grows factorially as |𝒱| increases, a brute force approach quickly becomes infeasible also for graph signals of moderate dimensions. To cope with this issue, in the sequel, we will introduce lower complexity methods based on: i) greedy approaches, and ii) convex relaxations.



  §.§.§ Greedy Sampling

In this section, we  will consider a numerically efficient, albeit sub-optimal, greedy algorithm to tackle the problem of selecting the sampling set. The greedy approach is described in Algorithm 1.

The simple idea underlying such method is to iteratively add to the sampling set those vertices of the graph that lead to the largest increment of an adopted performance metric, i.e., a specific set function f():2^𝒱→ℝ. We will set f()=- Tr{(_^H__v^-1_)^-1} if we use an A-optimality design as in (<ref>), or f()=σ_min(__) if we consider an E-optimality design as in (<ref>). In fact, since Algorithm 1 starts from the empty set, when ||<||, matrix _^H__v^-1_ is inevitably rank deficient, and its inverse does not exist. In this case, considering an A-optimality criterion, we can use f()=- Tr{(_^H__v^-1_)^†}, which becomes equivalent to (<ref>) when condition (<ref>) is satisfied.

In general, the performance of the greedy strategy will be sub-optimal with respect to an exhaustive search procedure. Nevertheless, if the set function f() satisfies some structural properties, the greedy Algorithm 1 can  be proved to be close to optimality. In  particular, submodularity plays a similar role in combinatorial optimization to convexity in continuous optimization and shares other features of concave functions <cit.>. 














Definition 1: A set function
f:2^𝒱→ℝ is submodular if and only if the derived set functions
f_a:2^𝒱∖{a}→ℝ

    f_a()=f(∪{a})-f()

are monotone decreasing, i.e., if for all subsets a,𝒜,ℬ⊆𝒱 it holds that

    𝒜⊆ℬ ⇒  f_a(𝒜)≥ f_a(ℬ).


Intuitively, submodularity is a diminishing returns property where  adding  an  element to a  smaller  set  gives  a  larger  gain than  adding  one  to  a  larger set. The maximization of monotone increasing submodular functions is still NP-hard, but the greedy heuristic can be used to obtain a solution that is provably close to optimality, with a solution having objective value within 1-1/e of the optimal combinatorial solution <cit.>.













Unfortunately, both set functions in (<ref>) and (<ref>) are not submodular functions [Interestingly, in a Bayesian recovery setting <cit.>, the negative of the MSE function was proved to be approximatively submodular.] <cit.>. Thus, even if the design criteria in (<ref>) and (<ref>) are useful to minimize the effect of noise [cf. (<ref>)] and model mismatching [cf. (<ref>)-(<ref>)], respectively, we do not have theoretical performance guarantees when applying Algorithm 1 to solve such problems. Nevertheless, in the literature of experimental design, a further design criterion is often considered as a surrogate for (<ref>) [or (<ref>)], which writes as:

    ^ D-opt    = ||=Mmax  log(_^H__v^-1_)
        = ||=Mmax  log(∑_i∈_,i_,i^H/r_i^2).

This is analogous to the so-called D-optimal design <cit.>, and is equivalent to one of the methods proposed in <cit.> for graph signals sampling. This design strategy aims at maximizing the volume of the parallelepiped built with the selected rows {_,i^H}_i∈ of matrix _ (weighted by the inverse of the noise variances {r_i^2}_i∈), and the rationale is to design a well suited basis for the graph signal that we want to estimate. Interestingly, the set function f()=log(_^H__v^-1_) is a monotone increasing submodular function <cit.>. Thus, in this case, the greedy approach in Algorithm 1 can be used to solve (<ref>) with provable performance guarantees. In the implementation of Algorithm 1, when ||<|| and matrix _^H__v^-1_ is rank-deficient, we can use f()=log pdet(_^H__v^-1_), which is equivalent to (<ref>) when sampling condition (<ref>) on perfect recovery is satisfied.



  §.§.§ Convex Relaxation


Another possible algorithmic solution to problems like (<ref>), (<ref>), (<ref>), is to resort to convex relaxation techniques, see, e.g., <cit.>. To this aim, let us introduce the indicator vector ={d_i}_i=1^N, such that the i-th entry is binary and given by d_i=1 if node i belongs to the sampling set , and d_i=0 otherwise. Using the indicator vector , we can build a general sampling design problem that can be cast as:

    min_   f() 
       s.t.    1^T=M, 
               ∈{0,1}^N,

where f()= Tr{(_^H diag()_v^-1_)^-1} for the A-optimal design [cf. (<ref>)], f()=-σ_min( diag()_) for the E-optimal design [cf. (<ref>)], and f()=-log(_^H diag()_v^-1_) for the D-optimal design [cf. (<ref>)]. Problem (<ref>) has still combinatorial complexity, due to the integer nature of the optimization variable . Nevertheless, we can simple relax the indicator variable  to be a real vector belonging to the hypercube [0,1]^N, thus leading to the following formulation:

    ∈[0,1]^Nmin   f() 
       s.t.      1^T=M.

It is now easy to check that problem (<ref>) is convex for all objective functions f() defined by the design criteria in
(<ref>), (<ref>), (<ref>), and its global solution can be found using efficient numerical methods <cit.>. Of course, since (<ref>) is a relaxed version of (<ref>), its real solution ^* might need a further selection/thresholding step in order to generate a valid integer vector, as required by (<ref>). For instance, a possible solution is to select the M sampling nodes as the ones associated with the M largest entries of ^*. Finally, one can also formulate the sampling design problem in the opposite way with respect to (<ref>). In particular, we might be interested in searching for the optimal indicator vector  that minimizes the number of collected samples, i.e., the ℓ_0 norm of the vector , under a performance requirement on the function f(), e.g, the MSE in (<ref>). This category of design problems takes the name of sparse sensing <cit.> and, using similar relaxation arguments as before, such criteria lead to convex optimization problems.

In the next section, we will illustrate some numerical results aimed at assessing the performance of the described sampling and recovery strategies.





 §.§ Numerical Results


In the sequel, we consider the application of the described sampling and recovery methods to two real graphs: a power network, and a road network.



Sampling over power grids. The first example involves the IEEE  118  Bus  Test  Case,  i.e.,  a  portion  of  the American Electric Power System (in the Mid-western US) as of December 1962. The graph is composed of 118 nodes (i.e., buses), its  topology  (i.e.,  transmission  lines  connecting buses) is depicted in Fig. <ref>(a) <cit.>, and the color of each node encodes the entries of the eigenvector of the Laplacian matrix associated to the second smallest eigenvalue (these entries highlight the presence of three distinct clusters in the network). As illustrated in <cit.>, the dynamics of the power generators give rise to smooth graph signals, so that the bandlimited assumption is justified, although in approximate sense. In our example, we randomly generate a lowpass signal with ||=12 and we take a number of samples equal to ||=12. The green squares correspond to the samples selected using the greedy Algorithm 1 and the A-optimality design in (<ref>). It is interesting to see how the method distributes samples over the clusters, and puts the samples, within each cluster, quite far apart from each other. Finally, we compare the reconstruction performance obtained by the considered greedy sampling strategies [cf. (<ref>), (<ref>), and (<ref>)] and by random sampling. To this aim we consider graph signal recovery, in the presence of an uncorrelated, zero mean Gaussian random noise with unit variance, and considering ||=20. Thus, Fig. <ref>(b) reports the MSE in (<ref>) versus the number of samples collected over the graph. As expected, the MSE decreases as the number of  samples increases. We can also notice how random sampling performs quite poorly, whereas the A-optimal design (and the D-optimal design) outperforms all other strategies.




Traffic flow prediction over road networks. The second example considers sampling of a portion of the road network in the neighborhood of Mazzini square, which is in the city of Rome, Italy. We have placed landmarks (nodes of the graph) over the streets in a regular fashion, and connected adjacent landmarks on the same lane and at the junctions, thus obtaining the graph topology depicted in Fig. <ref>(a). The signal lying on the vertices of the graph represents the flow (number of vehicles per unit of time) of cars passing through the landmarks during a period of 30 seconds, and was obtained using a realistic simulator of urban mobility, namely, SUMO <cit.>. The similarity of values of the signal over adjacent nodes makes the signal to be smooth, but only approximatively bandlimited. In this sense, there is a mismathcing between the observed signal and the bandlimited model used for processing. The goal is to infer the traffic situation over all the road network from a small number of collected samples. Thus, we consider a bandwidth equal to ||=30, and we take a number of samples equal to ||=40. The green squares in Fig. <ref>(a) correspond to the samples selected using the convex relaxation in (<ref>) and the E-optimality design in (<ref>). It is interesting to see how the method distributes almost uniformly the samples over the streets and the junctions. Finally, we compare the reconstruction performance obtained by the sampling strategies based on convex relaxation [cf. (<ref>), (<ref>), and (<ref>)] in the presence of model mismatching. To this aim, Fig. <ref>(b) reports the normalized MSE (NMSE), i.e., NMSE=-^2/^2, versus the graph signal bandwidth, and selecting ||=||. From Fig. <ref>(b), as expected, the NMSE decreases if we use a larger bandwidth. Furthermore, in this case, the E-optimal design outperforms all other strategies, as expected from (<ref>)-(<ref>).


 

 §.§ ℓ_1-norm reconstruction of graph signals


Let us consider now a different observation model, where a bandlimited graph signal ∈_ is observed everywhere, but a subset of nodes
 is strongly corrupted by noise, i.e.,

    = +_,

where the noise is arbitrary but bounded, i.e., _1<∞. This model is relevant, for example, in sensor networks, where a subset of sensors can be damaged
or  highly  interfered.  The  problem  in  this  case  is  whether  it is possible to recover the graph signal  exactly, i.e., irrespective of noise. Even though this is not a sampling problem, the solution is  still  related  to  sampling  theory.  Clearly,  if  the  signal 
is bandlimited  and  if  the  indexes  of  the  noisy  observations  are known, the answer is simple: ∈_
can be perfectly recovered from the noisy-free observations, i.e., by completely discarding the noisy observations, if the sampling theorem condition (<ref>)
holds true. But of course, the challenging situation occurs when the location of the noisy observations is not known. In such a
case, we may resort to an ℓ_1-norm minimization, by formulating the problem as follows <cit.>:

    = min_∈-_1.

We provide next some theoretical bounds on the cardinality of  and  enabling perfect recovery of the bandlimited graph signal using (<ref>). To this purpose, we recall the following lemma from <cit.>.

Let us define

μ := max_j ∈
 i ∈u_j (i),

where u_j (i) is the i-th entry of the j-th vector of the graph Fourier basis. If for some unknown , we have

    < 1/2 μ^2,

then the ℓ_1-norm reconstruction method (<ref>) recovers ∈ perfectly, i.e. =, for any arbitrary noise  present on at most  vertices.




An example of ℓ_1 reconstruction based on (<ref>) is useful to
grasp some interesting features. We consider the IEEE 118 bus graph in Fig. <ref>(a). The signal
is assumed to be bandlimited, with a spectral content limited to the first
|| eigenvectors of the Laplacian matrix. In Fig. <ref>, we report the behavior of the NMSE associated to the ℓ_1-norm
estimate in (<ref>) , versus the number of noisy samples, considering different values of bandwidth
||. As we can notice from
Fig. <ref>, for any value of ||, there exists a threshold value such that, if the number of noisy samples is lower than the threshold, the reconstruction of the signal is error free. As expected, a smaller signal bandwidth allows perfect reconstruction with a larger number of noisy samples.



§ ADAPTIVE SAMPLING AND RECOVERY



In this section, we consider processing methods capable to learn and track dynamic graph signals from a carefully designed, possibly time-varying, sampling set. To this aim, let us assume that, at each time n, noisy samples of the signal are taken over a (randomly) time-varying subset of vertices, according to the following model:

    [n] =    _[n](+[n])= _[n]__+_[n][n]

where _[n]= diag{d_i[n]}_i=1^N∈ℝ^N× N [cf. (<ref>)], with d_i[n] denoting a random sampling binary coefficient, which is equal to 1 if i∈[n], and 0 otherwise (i.e., [n] represents the instantaneous, random sampling set at time n); and [n]∈ℂ^N is zero-mean, spatially and temporally independent observation noise, with covariance matrix _v={r_1^2,…,r_N^2}. The estimation task consists in recovering the vector  (or, equivalently, its GFT _) from the noisy, streaming, and partial observations [n] in (<ref>). Following an LMS approach <cit.>, from (<ref>), the optimal estimate for _ can be found as the vector that solves the following optimization problem:

    min_s 𝔼 _[n]([n]-_)^2

where in (<ref>) we have exploited the fact that _[n] is an idempotent matrix for any fixed n [cf. (<ref>)]. An LMS-type solution optimizes (<ref>) by means of a stochastic steepest-descent procedure, relying only on instantaneous information. Thus, letting [n] be the current estimates of vector , the LMS algorithm for graph signals evolves as illustrated in Algorithm 2 <cit.>.



In the sequel, we illustrate how the design of the sampling strategy affects the reconstruction capability of Algorithm 2. To this aim, let us denote the expected sampling set by ={i=1,…,N  |  p_i>0}, i.e., the set of nodes of the graph that are sampled with a probability p_i=𝔼{d_i[n]} strictly greater than zero. Also, let 𝒮_c be the complement set of . Then, the following results illustrates the conditions for adaptive recovery of graph signals <cit.>.


Any -bandlimited graph signal can be reconstructed via the adaptive Algorithm 2 if and only if

    _ 𝒮_c__2 < 1,

i.e., if there are no -bandlimited signals that are perfectly localized on 𝒮_c.



Differently from batch sampling and recovery of graph signals, see, e.g.,
<cit.>, condition (<ref>) depends on the expected sampling set. In particular, it implies that there are no -bandlimited signals that are perfectly localized over the set _c. As a consequence, the adaptive Algorithm 2 with probabilistic sampling does not need to collect all the data necessary to reconstruct one-shot the graph signal at each iteration, but can learn acquiring the needed information over time. The only important thing required by condition (<ref>) is that a sufficiently large number of nodes is sampled in expectation.

We now illustrate the mean-square performance of Algorithm 2. The main results are summarized in the following Theorem <cit.>.


Assume spatial and temporal independence of the random variables extracted by the sampling process {d_i[n]}_i,n. Then, for any initial condition, Algorithm 2 is mean-square error stable if the sampling probability vector  and the step-size μ satisfy (<ref>) and


    0< μ < 2λ_min(_^H  diag()_)/λ^2_max(_^H  diag()_).


Furthermore, under a small step-size assumtpion, the MSE writes as:

    MSE   =lim_n→∞ 𝔼[n]-[n]^2 
       = μ/2  Tr[(_^H  diag()_)^-1_^H  diag()_v_]

and the convergence rate α is well approximated by

    α=
     1-2μλ_min(_^H  diag()_).



The results of Theorem 5 are instrumental to devise optimal probabilistic sampling strategies for Algorithm 2, which are described in the sequel.





 §.§ Probabilistic Sampling Strategies


We consider a sampling design that seeks for the probability vector  that minimizes the total sampling rate over the graph, i.e., 1^T, while guaranteeing a target performance in terms of MSE in (<ref>) and of convergence rate in (<ref>) <cit.>.

Then, the optimization problem can be cast as:

    min_p  1^T
       s.t.    λ_min(_^H  diag()_)≥1-α̅/2μ, 
    
    
                Tr[(_^H  diag()_)^-1_^H  diag()_v_]≤2γ/μ,
               0≤≤^max.

The first constraint imposes that the convergence rate of the algorithm is larger than a desired value, i.e., α in (<ref>) is smaller than a target value, say, e.g., α̅∈(0,1). Furthermore, as illustrated in <cit.>, the first constraint on the convergence rate also guarantees adaptive signal reconstruction, i.e., condition (<ref>) holds true. The second constraint guarantees a target mean-square performance, i.e., the MSE in (<ref>) must be less than or equal to a prescribed value, say, e.g., γ>0. Finally, the last constraint limits the probability vector to lie in the box p_i∈[0,p^max_i], for all i, with 0≤ p^max_i≤ 1 denoting an upper bound on the sampling probability at each node that might depend on external factors such as, e.g., limited energy, processing, and/or communication resources, node or communication failures, etc.

Unfortunately, problem (<ref>) is non-convex, due to the presence of the non-convex constraint on the MSE. To handle the non-convexity of (<ref>), we exploit an upper bound of the MSE function in (<ref>), given by:

    MSE() ≤  MSE()≜μ/2  Tr(_^H  diag()_v_)/λ_min(_^H  diag()_),   for all ∈ℝ^N.

Of course, replacing the MSE function with the bound (<ref>), the second constraint in (<ref>) is always satisfied. Furthermore, the function in (<ref>) has the nice property to be pseudo-convex, since it is the ratio between a convex and a concave function, which are both differentiable and positive for all  satisfying the other constraints <cit.>. Thus, exploiting the upper bound (<ref>), we can formulate a surrogate optimization problem for the selection of the probability vector , which can be cast as:

    min_p  1^T
       subject to
        λ_min(_^H  diag()_)≥1-α̅/2μ, 
         Tr(_^H  diag()_v_)/λ_min(_^H  diag()_)≤2γ/μ,
        0≤≤^max.

Since the sublevel sets of pseudo-convex functions are convex sets <cit.>, it is straightforward to see that the approximated problem (<ref>) is convex, and its global solution can be found using efficient numerical tools <cit.>.



 §.§ Distributed Adaptive Recovery


The implementation of Algorithm 2 would require to collect all the data {y_i[n]}_i:d_i[n]=1, for all n, in a single processing unit that performs the computation. In many practical systems, data are collected in a distributed network, and sharing local information with a central processor might be either unfeasible or not efficient, owing to the large volume of data, time-varying network topology, and/or privacy issues. Motivated by these observations, in this section we extend the LMS strategy in Algorithm 2 to a distributed setting, where the nodes perform the reconstruction task via online in-network processing, only exchanging data between neighbors defined over a sparse (but connected) communication network, which is described by the graph G_c=(𝒱,ℰ_c). Proceeding as in <cit.> to derive distributed solution methods for problem (<ref>), let us introduce local copies {_i}_i=1^N of the global variable , and recast problem (<ref>) in the following equivalent form:

    min_{s_i}_i=1^N  ∑_i=1^N  𝔼|d_i[n](y_i[n]-_,i^H_i)|^2
    
           subject to  _i=_j   for all i=1,…,N,    j∈𝒩_i,

where _,i^H is the i-th row of matrix _ (supposed to be known at node i, or computable in distributed fashion, see, e.g., <cit.>), and 𝒩_i={j|a_ij>0} is the local neighborhood of node i. To solve problem (<ref>), we consider an Adapt-Then-Combine (ATC) diffusion strategy <cit.>, and the resulting algorithm is reported in Algorithm 3. The first step in (<ref>) is an adaptation step, where the intermediate estimate ψ_i[n] is updated adopting the current observation taken by node i, i.e. y_i[n].  The second step is a diffusion step where the intermediate estimates ψ_j[n], from the (extended) spatial neighborhood 𝒩_i=𝒩_i⋃{i}, are combined through the weighting coefficients {w_ij}. Several possible combination rules have been proposed in the literature, such as the Laplacian or the Metropolis-Hastings weights, see, e.g. <cit.>, <cit.>, <cit.>. Finally, given the estimate _i[n] of the GFT at node i and time n, the last step produces the estimate x_i[n+1] of the graph signal value at node i [cf. (<ref>)].

Here, we assume that graphs G (i.e., the one used for GSP) and G_c (i.e., the one describing the communication pattern among nodes) might have in general distinct topologies. We remark that both graphs play an important role in the proposed distributed processing strategy (<ref>). First, the processing graph determines the structure of the regression data _,i used in the adaptation step of (<ref>). In fact, {_,i^H}_i are the rows of the matrix _, whose columns are the eigenvectors of the Laplacian matrix associated with the set of support frequencies . Then, the topology of the communication graph determines how information is spread all over the network through the diffusion step in (<ref>). This illustrates how, when reconstructing graph signals in a distributed manner, we have to take into account both the processing and communication aspects of the problem.



 §.§ Numerical Results


In this section, we first illustrate the performance of the probabilistic sampling method in (<ref>) over the IEEE 118 bus graph. Then, we consider an application to dynamic inference of brain activity.




Optimal probabilistic sampling. As a first example, let us consider an application to the IEEE 118 Bus Test Case in Fig. <ref>(a). The spectral content of the graph signal is assumed to be limited to the first ten eigenvectors of the Laplacian matrix of the graph. The observation noise in (<ref>) is zero-mean, Gaussian, with a diagonal covariance matrix _v, where each element is illustrated in Fig. <ref> (bottom). The other parameters are: μ=0.1, and γ=10^-3. Then, in Fig. <ref> (top and middle), we plot the optimal probability vector obtained solving (<ref>), for two different values of α̅. In all cases, the constraints on the MSE and convergence rate are attained strictly. From Fig. <ref> (top and middle), we notice how the method increases the sampling rate if we require a faster convergence (i.e., a smaller value of α̅); it also finds a very sparse probability vector and usually avoids to assign large sampling probabilities to nodes having large noise variances. Interestingly, with the proposed formulation, sparse sampling patterns are obtained thanks to the optimization of the sampling probabilities, which are already real numbers, without resorting to any relaxation of complex integer optimization problems [cf. (<ref>)].
















Inference of brain activity. The last example  presents  test  results  on  Electrocorticography (ECoG) data,  captured through  experiments  conducted  in  an epilepsy study <cit.>. Data were collected over a  period  of  five days, where the electrodes recorded 76 ECoG time series, consisting of voltage levels measured in different regions of the brain. Two  temporal intervals of interest were picked for analysis, namely, the preictal and ictal intervals. In the sequel, we focus on the ictal interval. Further  details  about  data  acquisition and  pre-processing are provided in <cit.>. The GFT matrix _ is learnt from the first 200 samples of ictal data, using the method proposed in <cit.>, and imposing a bandwidth equal to ||=30. In Fig. <ref>, we illustrate the true behavior of the ECoG present at an unobserved electrode chosen at random, over the first 400 samples of ictal data, along with estimate carried out using Algorithm 2 (with μ=1.5). The sampling set is fixed over time (i.e., p_i=1 for all i), and chosen according to the E-optimal design in (<ref>), selecting 32 samples. As we can notice from Fig. <ref>, the method is capable to efficiently infer and track the unknown dynamics of ECoG data at unobserved regions of the brain.






§ CONCLUSIONS


In this chapter, we have reviewed some of the methods recently proposed to sample and interpolate signals defined over graphs. First, we have recalled the conditions for perfect recovery under a bandlimited assumption. Second, we have illustrated sampling strategies, based on greedy methods or convex relaxations, aimed at reducing the effect of noise or aliasing on the recovered signal. Then, we considered ℓ_1-norm reconstruction, which allows perfect recovery of graph signals in the presence of a strong impulsive noise over a limited number of nodes. Finally, adaptive methods based on (possibly distributed) LMS strategies were illustrated to enable tracking and recovery of time-varying signals over graphs. Several interesting problems need further investigation, e.g., sampling and recovery in the presence of directed/switching topologies, sampling adaptation in time-varying scenarios, distributed implementations, and the extension of GSP methods to incorporate multi-way relationships among data, e.g., under the form of hypergraphs or simplicial complexes.



§ REFERENCES

Vancouver-Numbered-Style_3_



