Sampling of graph signals with successive local aggregations

I Introduction

Sampling (and subsequent interpolation) is a cornerstone problem in classical signal processing . The emergence of new fields of knowledge such as network science and big data is generating a pressing need to extend the results existing for classical time-varying signals to signals defined on graphs . This not only entails modifying the algorithms currently available for time-varying signals, but also gaining intuition on what concepts are preserved (and lost) when a signal is defined, not in the classical time grid, but in a more general graph domain.

This paper investigates the sampling and posterior recovery of signals that are defined in the nodes of a graph. The underlying assumption is that such signals admit a sparse representation in a (frequency) domain which is related to the structure of the graph where these signals reside. Most of the current efforts in this field have been focused on using the value of the signal observed at a subset of nodes to recover the signal in the entire graph . Our proposal in this paper is different. We present a new sampling method that accounts for the graph structure, can be run at a single node and only requires access to information of neighboring nodes. Moreover, we also show that the proposed method shares similarities with the classical sampling and interpolation of time-varying signals. When the graph corresponds to a directed cycle, which is the support of classical time-varying signals, our method is equivalent to classical sampling. When the graph is more general, the Vandermonde structure of the sampling matrix, which is critical to guarantee recovery in classical sampling , is preserved. Such a structure not only facilitates the interpolation process, but also helps to draw some connections between the proposed method and the sampling of time-varying signals. Sampling and interpolation are analyzed first in the absence of noise, where the conditions under which recovery is guaranteed are identified. The conditions depend both on the structure of the graph and the particular node taking the observations. They also reveal that one way to understand bandlimited graph signals is to think of signals that can be well approximated by only observing the value of the signal at a small neighborhood. We then analyze the sampling and reconstruction process when noise is present and when the specific frequencies where the signal is sparse are not known. For the noisy case, an interpolator based on the Best Linear Unbiased Estimator (BLUE) is designed and the effect on the corresponding error covariance matrix of different noise models is discussed. For the case of unknown frequency support, we also provide conditions under which the signal can be identified. This second problem falls into the category of sparse signal reconstruction where the main idea is to leverage the structure of the observation matrix to facilitate recovery. The last contribution is the design of a generalization of our sampling method that considers a subset of nodes, each of them taking multiple observations. Within that generalization, the approach of sampling a graph signal by observing the value of the signal at a subset of nodes can be viewed as a particular case. Hence, the generalization will also be useful to compare and establish relationships between existing approaches to sample signals in graphs and our proposed method.

The paper is organized as follows. Section II introduces the new aggregation sampling method, compares it to the existing selection sampling method and shows that for classical time-varying signals both methods are equivalent. Section III analyzes our sampling method in more detail and applies it to sample bandlimited graph signals. The analysis includes conditions for recovery, which are formally stated in Section III-C. Section IV investigates the effect of noise in aggregation sampling. It also discusses how to select sampling nodes and observation schemes that lead to a good recovery performance. Corresponding modifications in the interpolation in order to recover the signal when the support is not known are discussed in Section V. Section VI proposes a generalization under which the existing selection sampling and the proposed aggregation sampling can be viewed as particular cases. Several illustrative numerical results are presented in Section VII. A few concluding remarks are provided in Section VIII, which closes the paper.

II Sampling of graph signals

The graph ${\mathcal{G}}$ is endowed with a graph-shift operator ${\mathbf{S}}$ defined as an $N\times N$ matrix whose entry $(i,j)$, denoted as $S_{ij}$, can be nonzero only if $i=j$ or $(j,i)\in{\mathcal{E}}$. The sparsity pattern of the matrix ${\mathbf{S}}$ captures the local structure of ${\mathcal{G}}$ but we make no specific assumptions on the values of the nonzero entries of ${\mathbf{S}}$. Common choices for ${\mathbf{S}}$ are the adjacency matrix of the graph , the Laplacian , and its generalizations . The intuitive interpretation of ${\mathbf{S}}$ is that it represents a linear transformation that can be computed locally at the nodes of the graph. If ${\mathbf{y}}=[y_{1},\ldots,y_{N}]^{T}$ is defined as ${\mathbf{y}}={\mathbf{S}}{\mathbf{x}}$, then node $i$ can compute $y_{i}$ provided that it has access to the values of $x_{j}$ at its incoming neighbors $j\in{\mathcal{N}}_{i}$. We assume henceforth that ${\mathbf{S}}$ is diagonalizable, so that there exists a $N\times N$ matrix ${\mathbf{V}}$ and a $N\times N$ diagonal matrix $\boldsymbol{\Lambda}$ that can be used to decompose ${\mathbf{S}}$ as

In particular, (1) is true for normal matrices satisfying ${\mathbf{S}}{\mathbf{S}}^{H}={\mathbf{S}}^{H}{\mathbf{S}}$. In that case we have that ${\mathbf{V}}$ is unitary, which implies ${\mathbf{V}}^{-1}={\mathbf{V}}^{H}$, and leads to the decomposition ${\mathbf{S}}={\mathbf{V}}\boldsymbol{\Lambda}{\mathbf{V}}^{H}$.

A natural definition of sampling for a graph signal is to introduce a fat $K\times N$ selection matrix ${\mathbf{C}}$ and define the sampled signal as

If the matrix ${\mathbf{C}}$ is chosen as binary, i.e., with elements $C_{ij}\in\{0,1\}$, has a single nonzero element per row, and at most one nonzero element per column, then the signal $\bar{{\mathbf{x}}}$ is a selection of $K$ out of the $N$ elements of ${\mathbf{x}}$. In such a case, the ratio $K/N$ is the sampling rate of the signal. Uniform sampling amounts to setting ${\mathbf{C}}=[{\mathbf{e}}_{1},{\mathbf{e}}_{N/K+1},\ldots,{\mathbf{e}}_{N-N/K+1}]^{T}$ and the selection of the first $K$ elements of ${\mathbf{x}}$ is accomplished by setting ${\mathbf{C}}={\mathbf{E}}^{T}_{K}:=[{\mathbf{e}}_{1},\ldots,{\mathbf{e}}_{K}]^{T}$. We remark that, in general, it is not clear how to choose good selection matrices ${\mathbf{C}}$. This is in contrast to conventional sampling of signals in the time domain where uniform sampling is advantageous.

An equally valid, yet less intuitive, definition is to fix a node, say $i$, and consider the sampling of the signal seen by this node as the shift operator ${\mathbf{S}}$ is applied recursively. To describe this sampling methodology more clearly, define the $l$-th shifted signal ${\mathbf{y}}^{(l)}:={\mathbf{S}}^{l}{\mathbf{x}}$ and further define the $N\times N$ matrix

that groups the signal ${\mathbf{x}}$ and the result of the first $N-1$ applications of the shift operator. Associating the $i$-th row of ${\mathbf{Y}}$ with node $i$, we define the successively aggregated signal at $i$ as ${\mathbf{y}}_{i}:=({\mathbf{e}}_{i}^{T}{\mathbf{Y}})^{T}={\mathbf{Y}}^{T}{\mathbf{e}}_{i}$. Sampling is now reduced to the selection of $K$ out of the $N$ elements (rows) of ${\mathbf{y}}_{i}$, which we accomplish with a selection matrix ${\mathbf{C}}$ [cf. (2)]

We say that the signal $\bar{{\mathbf{y}}}_{i}$ samples ${\mathbf{x}}$ with successive local aggregations. This nomenclature follows from the fact that ${\mathbf{y}}^{(l)}$ can be computed recursively as ${\mathbf{y}}^{(l)}:={\mathbf{S}}{\mathbf{y}}^{(l-1)}$ and that the $i$-th element of this vector can be computed using signals associated with itself and its incoming neighbors,

We can then think of the signal ${\mathbf{y}}_{i}$ as being computed locally at node $i$ using successive variable exchanges with neighboring nodes. In fact, it is easy to show that $y_{i}^{(l)}$ can be expressed as a linear combination of the values of $x_{j}$ at nodes $j$ whose distance (number of hops) from node $i$ is less than or equal to $l$. This implies that the sampled signal $\bar{{\mathbf{y}}}_{i}$ in (4) is a selection of values that node $i$ can determine locally. An underlying idea behind the sampling in (4) is to incorporate the structure of the shift into the sampling procedure. Indeed, $\mathbf{S}$ and $\mathbf{y}^{(l)}$ play key roles in other graph-processing algorithms such as shift-invariant graph filters , where the output of the filter can be viewed a linear combination of the shifted signals $\mathbf{y}^{(l)}$.

To understand the difference between selection sampling [cf. (2)] and aggregation sampling [cf. (4)], it is instructive to consider their application to a signal defined in the time domain. We do so in the following section.

For a signal ${\mathbf{x}}$ defined on top of the directed cycle ${\mathcal{G}}_{dc}$, we consider selection sampling and aggregation sampling when using the shift operator ${\mathbf{S}}={\mathbf{A}}_{dc}$ and the uniform selection matrix ${\mathbf{C}}=[{\mathbf{e}}_{1},{\mathbf{e}}_{N/K+1},\ldots,{\mathbf{e}}_{N-N/K+1}]^{T}$. Illustrations of the respective sampling procedures are available in Figs. 1 and 2 for a signal with $N=6$ elements and sampling rate $K/N=1/2$.

In selection sampling we just multiply the graph signal ${\mathbf{x}}$ with the selection matrix ${\mathbf{C}}$ to obtain the sampled signal $\bar{{\mathbf{x}}}={\mathbf{C}}{\mathbf{x}}$ as indicated by (2). In aggregation sampling we consider subsequent applications of the shift matrix ${\mathbf{S}}={\mathbf{A}}_{dc}$. Each of these shifts amounts to rotating the signal clockwise so that the element at node $i$ moves to node $i+1$ for all $i<N$ and the element at node $N$ moves to node $1$. If we consider, e.g., node $i=1$, the first shift moves signal $x_{N}$ to this node so that $y_{1}^{(1)}=x_{N}$, the second shift moves signal $x_{N-1}$ to this node so that $y_{1}^{(2)}=x_{N-1}$ and so on. It follows that the aggregated signal ${\mathbf{y}}_{1}$ in (3) is given by ${\mathbf{y}}_{1}=[x_{1},x_{N},x_{N-1},\ldots,x_{2}]$. This is just a shift of the original signal ${\mathbf{x}}$, which, upon multiplication by the selection matrix ${\mathbf{C}}$ as per (4) results in a vector $\bar{{\mathbf{y}}}_{1}={\mathbf{C}}{\mathbf{y}}_{1}$ that contains the same elements that $\bar{{\mathbf{x}}}$ contains.

For the cycle graph and shift operator ${\mathbf{S}}={\mathbf{A}}_{dc}$ selection and aggregation sampling produce not only equivalent sampled signals but also reduce to conventional sampling. This is not a coincidence because both methods are designed as generalizations of conventional sampling. In general, selection sampling and aggregation sampling produce different outcomes. In selection sampling we move through nodes to collect samples at points uniquely identified by ${\mathbf{C}}$, whereas in aggregation sampling we move the signal through the graph while collecting samples at a fixed node. Observe that because aggregation sampling depends on the shift operator, it incorporates the structure of the graph into the sampling procedure. This is not true for selection sampling except for the choice of matrices ${\mathbf{C}}$ adapted to particular graphs.

III Sampling of bandlimited graph signals

Recovery of the original signal from its sampled version is possible under the assumption that the original signal admits a sparse representation. This section begins by introducing the concept of a bandlimited graph signal, which is sparse in the frequency domain, and establishing some connections with the concept of bandlimitedness in the classical time domain. Section III-B reviews briefly the recovery of a bandlimited graph signal for the case of selection sampling. Section III-C analyzes the recovery of a bandlimited graph signal for the case of aggregation sampling.

The common practice when addressing the problem of sampling signals in graphs is to suppose that the graph-shift operator $\mathbf{S}$ plays a key role in explaining the signals of interest $\mathbf{x}$. More specifically, that $\mathbf{x}$ can be expressed as a linear combination of a subset of the columns of $\mathbf{V}=[\mathbf{v}_{1},...,\mathbf{v}_{N}]$, or, equivalently, that the vector $\widehat{\mathbf{x}}=\mathbf{V}^{-1}\mathbf{x}$ is sparse. In this context, vectors $\mathbf{v}_{k}$ are interpreted as the graph frequency basis and $\widehat{x}_{k}$ as the corresponding signal frequency coefficients. To simplify exposition, it will be assumed throughout the paper that the active frequencies are the first $K$ ones, which are associated with the largest eigenvalues . Under this assumption, it holds that $\widehat{\mathbf{x}}=[\widehat{x}_{1},...,\widehat{x}_{K},0,...,0]^{T}$. However, the results presented in the paper can be applied to any set of active frequencies $\mathcal{K}$ of size $K$ provided that $\mathcal{K}$ is known. For convenience, we define $\mathbf{V}_{K}:=[\mathbf{v}_{1},...,\mathbf{v}_{K}]$ and $\widehat{\mathbf{x}}_{K}:=[\widehat{x}_{1},...,\widehat{x}_{K}]^{T}$ so that we may write $\widehat{\mathbf{x}}=[\widehat{\mathbf{x}}^{T}_{K}~{}|~{}\mathbf{0}_{1\times N-K}]^{T}$. For $\widehat{\mathbf{x}}$ to be sparse, it is reasonable to assume that $\mathbf{S}$ is involved in the generation of $\mathbf{x}$.

When $\mathcal{G}=\mathcal{G}_{dc}$, it can be easily shown that setting the shift operator either to $\mathbf{S}=\mathbf{A}_{dc}$ or to $\mathbf{S}=\mathbf{L}_{dc}:=\mathbf{I}-\mathbf{A}_{dc}$ gives rise to the Fourier basis $\mathbf{F}$. More formally, that the right eigenvectors of $\mathbf{S}$ satisfy $\mathbf{V}=\mathbf{F}$, with $F_{ij}:=\frac{1}{\sqrt{N}}e^{+\mathfrak{j}\frac{2\pi}{N}(i-1)(j-1)}$ and $\mathfrak{j}:=\sqrt{-1}$. Selecting $\mathbf{S}=\mathbf{A}_{dc}$ has the additional advantage of satisfying $\Lambda_{ii}=e^{-\mathfrak{j}\frac{2\pi}{N}(i-1)}$, i.e., the eigenvalues of the shift operator correspond to the classical discrete frequencies. Interpretations for the eigenvalues of the Laplacian matrix $\mathbf{L}_{dc}$ also exist .

III-B Selection sampling of bandlimited graph signals

Under the selection sampling approach , sampling a graph signal amounts to setting $\bar{\mathbf{x}}=\mathbf{C}\mathbf{x}$ [cf. (2)]. Since the $K\times N$ binary selection matrix $\mathbf{C}$ indexes the nodes that are observed, the issue then is how to design $\mathbf{C}$, i.e., which nodes to select, and how to recover the original signal $\mathbf{x}$ from its samples $\bar{\mathbf{x}}$.

To answer these questions, it is assumed that the signal $\mathbf{x}$ is bandlimited, so that it can be expressed as a linear combination of the $K$ principal eigenvectors in $\mathbf{V}$. The sampled signal $\bar{\mathbf{x}}$ is then $\bar{\mathbf{x}}=\mathbf{C}\mathbf{x}=\mathbf{C}\mathbf{V}_{K}\widehat{\mathbf{x}}_{K}$. Clearly, if the matrix $\mathbf{C}\mathbf{V}_{K}$ is invertible, then $\widehat{\mathbf{x}}_{K}$ can be recovered from $\bar{\mathbf{x}}$. Once the coefficients $\widehat{\mathbf{x}}_{K}$ are known, the signal in the original domain can be found as $\mathbf{x}=\mathbf{V}_{K}\widehat{\mathbf{x}}_{K}$. Combining the previous equations, we have

The expression in (6) shows how the original signal can be interpolated from its samples. For the previous equation to hold true, the matrix $\mathbf{C}\mathbf{V}_{K}$ has to be invertible. Hence, the key for guaranteeing perfect signal reconstruction is to select a subset of nodes such that the corresponding rows in $\mathbf{V}_{K}$ are linearly independent. In the classical domain of time-varying signals, the (Fourier) basis has a Vandermonde structure, both row-wise and column-wise. This readily implies that any subset of $K$ rows will give rise to a (row-wise) Vandermonde matrix and, hence, invertibility is guaranteed. However, for an arbitrary graph this is not guaranteed and algorithms to select a specific subset that guarantees recovery are required . The role of the Vandermonde structure of the sampling matrix will be analyzed in more detail in the ensuing sections.

III-C Aggregation sampling of bandlimited graph signals

As explained in (4), under the aggregation approach the sampled signal is formed by observations of the shifted signals $\mathbf{y}^{(l)}=\mathbf{S}^{l}\mathbf{x}$ taken at a given node $i$. Under this second approach, the graph-shift operator $\mathbf{S}$ plays a key role not only in explaining and recovering $\mathbf{x}$, but also in sampling $\mathbf{x}$. Another reason to consider this scheme is that the entries of $\mathbf{y}^{(l)}$ can be found by sequentially exchanging information among neighbors. This implies that: a) for setups where graph vertices correspond to nodes of an actual network, the procedure can be implemented distributedly; and b) if recovery is feasible, the observations at a single node can be used to recover the signal in the entire graph.

Mimicking the approach in the previous section, we first analyze how the bandlimitedness of ${\mathbf{x}}$ is manifested on the sampled signal. Then, we identify under which conditions recovery is feasible and describe the corresponding interpolation algorithm. For the ease of exposition, the dependence of $\mathbf{y}_{i}$ on $\widehat{\mathbf{x}}$ is given in the form of a lemma.

Define the $N\times 1$ vector $\boldsymbol{\upsilon}_{i}:=\mathbf{V}^{T}\mathbf{e}_{i}$, which collects the values of the frequency basis $\{\mathbf{v}_{k}\}_{k=1}^{K}$ at node $i$, and the $N\times N$ (column-wise) Vandermonde matrix

Then, the shifted signal $\mathbf{y}_{i}$ can be expressed as

Proof : Using the spectral decomposition of $\mathbf{S}$, signal $\mathbf{y}^{(l)}$ can be written as

Based on the definitions of $\mathbf{y}_{i}$ and $\boldsymbol{\upsilon}_{i}$, it follows that

Since the $l$-th column of matrix $\mathbf{Y}$ is $\mathbf{y}^{(l-1)}$, it can be written as $(\mathbf{V}\boldsymbol{\Lambda}^{l-1})\widehat{\mathbf{x}}$ [cf. (9)]. Hence, the $l$-th column of matrix $(\mathbf{V}^{-1}\mathbf{Y})$ can be written as $\boldsymbol{\Lambda}^{l-1}\widehat{\mathbf{x}}$ or, equivalently, as $\text{diag}(\widehat{\mathbf{x}})[\lambda_{1}^{l-1},...,\lambda_{N}^{l-1}]^{T}$. Leveraging the fact that the vector containing the $l$-th power of the eigenvalues corresponds to the row $l+1$ of matrix $\boldsymbol{\Psi}$, the shifted signal $\mathbf{y}_{i}$ can be expressed as

Notice that while in Section III-B the relationship between the sparse frequency coefficients $\widehat{\mathbf{x}}$ and the signal to be sampled was simply given by $\mathbf{x}=\mathbf{V}\widehat{\mathbf{x}}$, now it is given by $\mathbf{y}_{i}=\boldsymbol{\Psi}\text{diag}(\boldsymbol{\upsilon}_{i})\widehat{\mathbf{x}}$.

Next, we use Lemma 1 to identify under which conditions recovery is feasible. To do this, let us define the $N\times K$ matrix $\boldsymbol{\Psi}_{i}=\boldsymbol{\Psi}\text{diag}(\boldsymbol{\upsilon}_{i})\mathbf{E}_{K}$. Then, the sampled signal $\bar{\mathbf{y}}_{i}$ is

where $\mathbf{C}$ is the binary $K\times N$ selection matrix, and $\widehat{\mathbf{x}}_{K}$ the vector collecting the non-zero components of $\widehat{\mathbf{x}}$. To simplify exposition, for the time being we will assume that $\mathbf{C}\!=\!\mathbf{E}_{K}^{T}$, i.e., that the observations correspond to the original signal and the first $K\!-\!1$ shifts. This assumption can be relaxed, as discussed in Remark 1.

If matrix $\mathbf{C}\boldsymbol{\Psi}_{i}$ is invertible, then $\widehat{\mathbf{x}}_{K}$ can be recovered from $\bar{\mathbf{y}}_{i}$ [cf. (12)] and, once $\widehat{\mathbf{x}}_{K}$ is known, $\mathbf{x}$ can be found as $\mathbf{x}=\mathbf{V}_{K}\widehat{\mathbf{x}}_{K}$. Combining the previous expressions, we have [cf. (6)]

The expression in (13) shows how the original signal can be interpolated from its samples. As already stressed, for the previous equation to hold true, the matrix $\mathbf{C}\boldsymbol{\Psi}_{i}$ has to be invertible. Hence, the key for guaranteeing perfect signal reconstruction is to select samples such that the corresponding rows in $\boldsymbol{\Psi}_{i}$ are linearly independent. While for the selection sampling described in Section III-B there is no straightforward way to check the invertibility of $\mathbf{C}\mathbf{V}_{K}$ (existing algorithms typically do that by inspection ), for the aggregation sampling described in (8)-(13) the invertibility of $\mathbf{C}\boldsymbol{\Psi}_{i}$ can be guaranteed if the conditions presented in the following proposition hold.

Let $\mathbf{x}$ and $\bar{\mathbf{y}}_{i}$ be, respectively, a bandlimited graph signal with at most $K$ non-zero frequency components and the output of the sampling process defined in (12). Then, the $N$ entries of signal $\mathbf{x}$ can be recovered from the $K$ samples in $\bar{\mathbf{y}}_{i}$ if the two following conditions hold i) The first $K$ eigenvalues of the graph-shift operator $\mathbf{S}$ are distinct; i.e., $\lambda_{i}\neq\lambda_{j}$ for all $i\neq j$, $i\leq K$ and $j\leq K$. ii) The $K$ first entries of $\boldsymbol{\upsilon}_{i}$ are non-zero.

Proof : To prove the proposition it suffices to show that i) and ii) guarantee the invertibility of $\mathbf{C}\boldsymbol{\Psi}_{i}$ [cf. (13)]. Matrix $\mathbf{C}\boldsymbol{\Psi}_{i}$ can be understood as the multiplication of two matrices: matrix $(\mathbf{C}\boldsymbol{\Psi}\mathbf{E}_{K})$ and matrix $(\mathbf{E}_{K}^{T}\text{diag}(\boldsymbol{\upsilon}_{i})\mathbf{E}_{K})$. It is immediate that condition ii) guarantees that the second matrix is invertible. Moreover, condition i) guarantees invertibility of the first matrix. To see this, note that $(\boldsymbol{\Psi}\mathbf{E}_{K})$ is a $N\times K$ (column-wise) Vandermonde matrix. Hence $\mathbf{C}(\boldsymbol{\Psi}\mathbf{E}_{K})$ is a selection of the first $K$ rows of $(\boldsymbol{\Psi}\mathbf{E}_{K})$, which is also Vandermonde. Any square Vandermonde matrix has full rank provided that the basis (i.e., the eigenvalues of $\mathbf{S}$) are distinct, as required in condition i). ∎

One of the implications of the proposition is that there is no need to compute or observe the entire vector $\mathbf{y}_{i}$, since its first $K$ entries suffice to guarantee recovery.

The conditions in Proposition 1 are not difficult to check and they provide additional insights on the behavior of the sampling and interpolation procedure. Condition i) refers to the structure of the entire graph. It states that if a graph has two identical frequencies and the signal of interest is a linear combination of both of them, then the $K\times K$ matrix $(\mathbf{C}\boldsymbol{\Psi}\mathbf{E}_{K})$ cannot be inverted and the sampling procedure will fail, regardless of the chosen node. Note that this problem is not present in classical sampling of time-varying signals, because the eigenvalues of the Fourier Vandermonde matrix associated with $\mathbf{S}=\mathbf{A}_{dc}$ are always distinct. Condition ii) refers to the specific node where the samples of the shifted signal are observed. It basically states that any node in the network can be used to sample the signal provided that $(\mathbf{e}_{k}^{T}\boldsymbol{\upsilon}_{i})\neq 0$ for $k=1,\ldots,K$; i.e., that the chosen node participates in the specific frequencies on which signal $\mathbf{x}$ is expressed. It also points to the fact that if $|\mathbf{e}_{k}^{T}\boldsymbol{\upsilon}_{i}|$ are non-zero but small, selecting $i$ as the sampling node may give rise to interpolations that are potentially unstable if noise is present; see Section IV. For the particular case when $\mathbf{S}=\mathbf{A}_{dc}$, condition ii) is satisfied since all the entries of the Fourier basis are non-zero.

III-D Discussion

Suppose that we know that $\mathbf{x}$ is indeed $K$-bandlimited; i.e., that it can be expressed as a linear combination of the $K$ first frequency basis vectors $\mathbf{v}_{1},\ldots,\mathbf{v}_{K}$. Then, Proposition 1 states that a single node, say the $i$-th one, can reconstruct the entire graph signal just from its own signal $x_{i}$ and $K-1$ exchanges with its neighbors. Note that one of the consequences of this result is that linear combinations of signals at nodes that are in a neighborhood of radius $K-1$ suffice to reconstruct the entire graph signal. To be specific, suppose that $\mathbf{x}=\alpha\mathbf{v}_{1}$, which represents the extreme case of a 1-bandlimited signal. Then, it follows that $[\mathbf{y}_{i}]_{1}=x_{i}=\alpha[\boldsymbol{\upsilon}_{i}]_{1}$, from where $\alpha$ can be found – and $\mathbf{x}$ reconstructed – as long as $[\boldsymbol{\upsilon}_{i}]_{1}\neq 0$ [cf. condition ii) in Proposition 1]. This implies that a node can reconstruct a 1-bandlimited signal based solely on the value that this signal takes at the node. For the case of a 2-bandlimited signal where $\mathbf{x}=\alpha_{1}\mathbf{v}_{1}+\alpha_{2}\mathbf{v}_{2}$, Proposition 1 guarantees reconstruction based on $[\mathbf{y}_{i}]_{1}$ and $[\mathbf{y}_{i}]_{2}$, which only contain information about the signal at node $i$ and at its neighbors. Therefore, one can understand bandlimited graph signals as signals that can be identified locally by relying on observations within a given number of hops. Note that this does not necessarily imply that the variation of the signal among close-by nodes is small, it only means that the pattern of variation can be inferred just by looking at close-by nodes. This discussion will be revisited in Section V. For the recovery to be implemented locally too, the nodes need to know $\mathbf{V}_{K}$ and $\{\lambda_{k}\}_{k=1}^{K}$, i.e. the structure of the graph where the signal resides.

We may decompose the interpolator $\mathbf{V}_{K}(\mathbf{C}\boldsymbol{\Psi}_{i})^{-1}$ in (13) into three factors $\mathbf{V}_{K}(\mathbf{E}_{K}^{T}\text{diag}(\boldsymbol{\upsilon}_{i})\mathbf{E}_{K})^{-1}(\mathbf{C}\boldsymbol{\Psi}\mathbf{E}_{K})^{-1}$ to reveal that it can be computed in closed-form, since a non-zero diagonal matrix can be trivially inverted and closed-form expressions for the inverse of a Vandermonde matrix exist . Moreover, notice that one of these three factors is related to the structure of the graph and the support where the signal is bandlimited $\mathbf{V}_{K}$; one is related to the structure of the graph, the support of the signal and the subset of observations $(\mathbf{C}\boldsymbol{\Psi}\mathbf{E}_{K})^{-1}$; and the third one depends on the specific node where the samples are taken $(\mathbf{E}_{K}^{T}\text{diag}(\boldsymbol{\upsilon}_{i})\mathbf{E}_{K})^{-1}$.

IV Sampling and interpolation in the presence of noise

When sampling in the absence of noise, two main questions are how to recover the signal from its samples and the conditions under which recovery is feasible. When the samples are noisy, perfect reconstruction is, in general, unfeasible and new issues arise. In Section IV-A we estimate the noisy signal through interpolation via the Best Linear Unbiased Estimator (BLUE) for a general noise model. In Section IV-B, we specify noise models that are likely to arise in graph domains. Then, in Section IV-C, we discuss the effect on the interpolation error of selecting the sampling node and the rows of the selection matrix.

Key to design the interpolator in the presence of noise is to notice that the relation between the observed samples $\bar{\mathbf{z}}_{i}$ and the original signal $\mathbf{x}$ is given by

The BLUE estimator of $\widehat{\mathbf{x}}_{K}$, which minimizes the least squares error, is given by

provided that the inverse in (16) exists. Additionally, for the particular case of Gaussian noise in (14), the estimator in (16) coincides with the Minimum Variance Unbiased (MVU) estimator which attains the Cramér-Rao lower bound. In this case, it also holds true that the inverse of the error covariance matrix associated with (16) corresponds to the Fisher Information Matrix (FIM) . Clearly the larger the number of rows in (14), the better the estimation is. When the selection matrix $\mathbf{C}$ selects exactly $K$ rows (and not more), (16) reduces to

After obtaining $\hat{\widehat{\mathbf{x}}}_{K}^{(i)}$ – either via (16) or (17) –, the time signal recovered at the $i$-th node $\hat{\mathbf{x}}^{(i)}$ can be found as

Note that the error covariance matrix $\mathbf{R}_{e}^{(i)}$ depends on the noise model, the frequencies of the graph (eigenvalues of the shift operator), the node taking the observations, and the sample-selection scheme adopted (cf. Remark 1).

The error covariance matrix enables us to assess the performance of the estimation. Smaller errors give rise to better estimators. However, there exist multiple alternatives to quantify the error, as analyzed by the theory of optimal design of experiments . The most common approach is to find an estimator which minimizes the trace of the covariance matrix

which corresponds to the minimization of the Mean Square Error (MSE). Other common error metrics based on the error covariance matrix are the largest eigenvalue

and the inverse of the trace of its inverse

Notice that the error metrics $e_{3}$ and $e_{4}$ are computed based on the error covariance matrix for the frequency estimator $\widehat{\mathbf{R}}_{e}^{(i)}$ instead of the time estimator since $\mathbf{R}_{e}^{(i)}$ is a singular matrix [cf. (20)].

IV-B Noise models

The results presented so far consider a general error covariance matrix $\mathbf{R}_{w}^{(i)}$, so that they can be used regardless of the color of the noise. In this section, we present three particular examples of interest.

White noise in the observed signal $\mathbf{z}_{i}$. This implies that $\mathbf{w}_{i}$ is white and therefore $\mathbf{R}_{w}^{(i)}=\sigma^{2}\mathbf{I}$, with $\sigma^{2}$ denoting the noise power. In this case, the $K\!\times\!K$ matrix $\bar{\mathbf{R}}_{w}^{(i)}$ is given by

White noise in the original signal $\mathbf{x}$. With $\mathbf{w}$ denoting the white additive noise present in $\mathbf{x}$, we can use the linear observation model to write $\mathbf{w}_{i}=\boldsymbol{\Psi}\text{diag}(\boldsymbol{\upsilon}_{i})\mathbf{V}^{-1}\mathbf{w}$. Then, the $N\times N$ error correlation matrix is simply given by $\mathbf{R}_{w}^{(i)}=\sigma^{2}\boldsymbol{\Psi}\text{diag}(\boldsymbol{\upsilon}_{i})\mathbf{V}^{-1}(\mathbf{V}^{-1})^{H}\text{diag}(\boldsymbol{\upsilon}_{i})^{H}\boldsymbol{\Psi}^{H}$. When the shift is a normal matrix, $\mathbf{V}$ is unitary and the previous expression reduces to $\mathbf{R}_{w}^{(i)}=\sigma^{2}\boldsymbol{\Psi}|\text{diag}(\boldsymbol{\upsilon}_{i})|^{2}\boldsymbol{\Psi}^{H}$. As before, the $K\times K$ error correlation matrix is obtained just by selecting the rows and columns of the former,

The previous expressions show not only that the noise is correlated, but also that the correlation depends on the graph structure (eigenvalues and eigenvectors of $\mathbf{S}$), the node collecting the observations, and the specific selection of observations.

White noise in the active frequency coefficients $\widehat{\mathbf{x}}_{K}$. With $\widehat{\mathbf{w}}_{K}$ denoting the white additive noise present in $\widehat{\mathbf{x}}_{K}$, we can use the linear observation model to write $\mathbf{w}_{i}=\boldsymbol{\Psi}\text{diag}(\boldsymbol{\upsilon}_{i})\mathbf{E}_{K}\widehat{\mathbf{w}}_{K}\!=\!\boldsymbol{\Psi}_{i}\widehat{\mathbf{w}}_{K}$. It follows that the $N\!\times\!N$ and $K\!\times\!K$ error covariance matrices are $\mathbf{R}_{w}^{(i)}=\sigma^{2}\boldsymbol{\Psi}_{i}\boldsymbol{\Psi}_{i}^{H}$ and

This model can be appropriate for scenarios where the signal of interest is the output of a given “graph process” – e.g., a diffusion process – and the noise is present in the input of that process. This noise model can also arise when the signal to be sampled has been previously processed with a low-pass graph filter .

There are many other noise models that can be of interest in graph setups. For example, the error covariance can be a linear combination of the previous covariance matrices (noise is present in both the original signal and the observation process). Alternatively, the noise at a specific node can be also rendered dependent on the number of neighbors. This last situation would be reasonable, for example, in distributed setups where the information of neighboring nodes is exchanged via noisy channels.

IV-C Selection of the sampling set

The two elements that define the set of samples to be interpolated are: the sampling node, i.e., the node $i$ which aggregates the information; and the sample-selection scheme, i.e., the elements of $\mathbf{y}_{i}$ selected by $\mathbf{C}$.

The recovery results in Section III-C show that any node $i$ can be used to sample and recover the entire graph signal, provided that the entries of $\boldsymbol{\upsilon}_{i}$ corresponding to the active frequencies in $\widehat{\mathbf{x}}$ are non-zero. However, when noise is present, the error covariance matrix $\mathbf{R}_{e}^{(i)}$ – which is the key element to evaluate the quality of the interpolation – is different for each $i$. In this context, it is reasonable to select as a sampling node one leading to a small error. Note that selecting the best one will only require the computation of $N$ closed-form expressions which involve matrix inversions. In scenarios where computational complexity is a limiting factor, the structure of the noise correlation, as well as the structure of the interpolation matrix, can be exploited to reduce considerably the computational burden. E.g., for the case where white noise is present in the active frequency coefficients, when substituting (27) into (19) and (20), it follows that

Consequently, for this particular noise model, the estimator performance is independent of the node choice. This is true for every error metric [cf. (21)-(24)]. The result is intuitive: given that the noise and the signal are present in the same frequencies, it is irrelevant if a node amplifies or attenuates a particular frequency. Differently, if the white noise is present in the observed signal, we can substitute (25) into (19) to obtain

Thus, if we are interested in minimizing, e.g., the error metric $e_{4}$ [cf. (24)], our objective may be reformulated as finding the optimal node $i^{*}$ such that

For a selection matrix of the form $\mathbf{C}=\mathbf{C}_{K}(n_{0},N_{0})$ (cf. Remark 1), the $k$-th diagonal element of the matrix in (30) can be written as $|[\boldsymbol{\upsilon}_{i}]_{k}|^{2}\sum_{m=0}^{K-1}|\lambda_{k}|^{2\,(n_{0}+mN_{0})}$. The trace is simply the sum of those elements, so that, using the closed form for a geometric sum, (30) can be rewritten as

Thus, the optimal sampling node $i^{*}$ will be one with large values of $|[\boldsymbol{\upsilon}_{i}]_{k}|$ for the active frequencies $k\leq K$. The relative importance of frequency $k$ is given by the fraction in (31), which depends on the modulus of the associated eigenvalue and the structure of the selection matrix $\mathbf{C}$ (values of $n_{0}$ and $N_{0}$).

IV-C2 Design of the sample-selection scheme

The error covariance matrix, and hence the different error metrics presented in (21)-(24), depend on the selection matrix $\mathbf{C}$. By changing $\mathbf{C}$ one can tradeoff the quality of a given sample and the detrimental effect of the corresponding noise. The specific set of samples that minimizes the error will in general depend on the error metric chosen.

Recall that any matrix in the set of admissible selection matrices $\mathcal{C}_{K}$ defined in Remark 1 is guaranteed to lead to a feasible recovery according to the conditions stated in Proposition 1. If $\mathbf{C}$ is not constrained to belong to $\mathcal{C}_{K}$, the number of candidate matrices is $N$ choose $K$. However, $\mathcal{C}_{K}$ has a much smaller cardinality: $N_{0}$ can take at most $N/K$ values, and $n_{0}$ at most $(N-N_{0}(K-1))$. Moreover, as it was the case for the sampling node selection, in some cases the noise structure can be exploited to readily determine the optimal observation strategy. E.g., for the case where white noise is present in the active frequencies, it is immediate to see that the performance is independent of the sample-selection scheme [cf. (28)]. For the case where white noise is present in the observed signal, let us assume that the selection matrix is given by $\mathbf{C}=\mathbf{C}_{K}(n_{0},N_{0})$ where $N_{0}$ is fixed and we want to design $n_{0}$.

For the first equality we have used that the product of diagonal matrices is commutative and for the second one that right and left multiplying by the canonical matrix amounts to selecting the columns and rows of the multiplied matrix. Using (32), we have that

which results in the following optimal strategy for the solution of $e_{3}$: if $\prod_{k=1}^{K}|\lambda_{k}|^{2}\leq 1$ then $n_{0}^{*}=1$, otherwise $n_{0}^{*}$ should be as large as possible; see Remark 2. Equivalently, the optimal strategy states that if an application of the shift operator has an overall effect of amplification in the active frequencies, then we should aim to apply it as many times as possible, whereas if the opposite is true, we should avoid its application.

Needless to say, one can also look at selection matrices that are not always guaranteed to lead to a Vandermonde structure, i.e., matrices not in $\mathcal{C}_{K}$. In that case, the problem can be formulated as a binary optimization over $\mathbf{C}$, which is typically challenging. If the size of the space search ($N$ choose $K$) is not too large, the problem can be solved by exhaustive search – first by checking that the matrix guarantees recovery and then evaluating the corresponding error covariance. For more general cases, a reasonable approach is to formulate the problem, relax it, and exploit the problem structure to find a good approximate solution efficiently. The problem formulation and the subsequent relaxation will depend on the specific optimality criteria. Although of interest, developing approximate algorithms to design the selection matrix $\mathbf{C}$ is out of the scope of this paper and is left as future work.

It is worth stressing that the sample-selection scheme that minimizes the error does not have to be the same for all nodes. Both the selection of the sampling node and the sampling shifts can be combined to obtain the best local reconstruction across all nodes in the graph.

Designing $\mathbf{C}$ entails the selection of a subset of $K$ entries out of the $N$ entries in $\mathbf{y}_{i}$. However, $\mathbf{y}_{i}$ has only $N$ entries because $\mathbf{Y}$ has only $N$ columns [cf. (3)]. Strictly speaking, there is no need to impose this restriction and more columns could be added to $\mathbf{Y}$. As a matter of fact, if for a given noisy graph signal the application of the shift operator $\mathbf{S}$ attenuates the noise while amplifying the signal, the sampling procedure will benefit from further applications of $\mathbf{S}$, even beyond the size of the graph $N$. In practice, the maximum number of applications will be limited by the computational and signaling cost associated with the application of the shift.

V Identifying the support of the graph signal

In the previous sections, it has been assumed that the frequency support of the bandlimited signal corresponded to the $K$ principal eigenvectors, which are the ones associated with the largest eigenvalues. However, the results presented also hold true as long as the basis support, i.e., the frequencies that are present in $\mathbf{x}$, are known. To be specific, let $\mathcal{K}:=\{k_{1},\ldots,k_{K}\}$ denote the set of indexes where the signal $\mathbf{x}$ is sparse and, based on it, define the $N\times K$ matrices $\mathbf{V}_{\mathcal{K}}:=[\mathbf{v}_{k_{1}},\ldots,\mathbf{v}_{k_{K}}]$ and $\mathbf{E}_{\mathcal{K}}:=[\mathbf{e}_{k_{1}},\ldots,\mathbf{e}_{k_{K}}]$. Then, all the results presented so far hold true if $\mathbf{V}_{K}$ is replaced with $\mathbf{V}_{\mathcal{K}}$, and $\mathbf{E}_{K}$, when used to select the active frequencies, is replaced with $\mathbf{E}_{\mathcal{K}}$.

A related but more challenging problem is to design the sampling and interpolation procedures when the frequency support $\mathcal{K}$ is not known. Generically, this problem falls into the class of sparse signal reconstruction . However, the particularities of our setup can be exploited to achieve stronger results. In particular, for the sampling procedure proposed in this paper, the so-called sensing matrix – the one relating the signal of interest to the observed samples – has a useful Vandermonde structure that can be exploited.

Consider the noiseless aggregation sampling of Section III-C, where we know that $\widehat{\mathbf{x}}$ is $K$-sparse but we do not know the support of the $K$ non-zero entries [cf. (12)]

For the case where the support is known, it was shown that a selection matrix $\mathbf{C}$ that picks the first $K$ rows of $\boldsymbol{\Psi}$ is enough for perfect reconstruction (cf. Proposition 1).

If we reformulate the recovery problem as

for the unknown support case, there is no guarantee that the solution $\widehat{\mathbf{x}}^{*}$ coincides with the $K$-sparse representation of the observed signal. When the frequency support is known, and provided that the selection matrix satisfies the conditions in Remark 1, selecting $K$ rows of $\boldsymbol{\Psi}\text{diag}(\boldsymbol{\upsilon}_{i})\mathbf{E}_{\mathcal{K}}$ leads to a (one-to-one) invertible transformation. When the support is unknown, guaranteeing identifiability requires selecting a higher number of rows (samples) . The following proposition states this result formally. To simplify notation, the proposition assumes that $K\leq N/2$, but the result holds true also when that is not the case.

Let $\mathbf{x}$ and $\mathbf{C}$ be, respectively, a bandlimited graph signal with at most $K$ non-zero frequency components and a selection matrix with $2K$ rows of the form $\mathbf{C}=\mathbf{C}_{2K}(n_{0},N_{0})$ (cf. Remark 1). Then, if all the entries in $\boldsymbol{\upsilon}_{i}$ are non-zero and all the eigenvalues of $\mathbf{S}$ are non-zero and satisfy that $\lambda_{k}^{N_{0}}\neq\lambda_{k^{\prime}}^{N_{0}}$ for all $k\neq k^{\prime}$, it holds that i) the solution to (35) is unique; and ii) the original graph signal can be recovered as $\mathbf{x}=\mathbf{V}\widehat{\mathbf{x}}^{*}$.

Proof : The proof proceeds into two steps. The first step is to show that any selection of $2K$ columns of the $2K\times N$ matrix $\mathbf{M}:=\mathbf{C}\boldsymbol{\Psi}\text{diag}(\boldsymbol{\upsilon}_{i})$ has rank $2K$ and, hence, it leads to an invertible $2K\times 2K$ matrix. To prove this, let $\mathcal{F}=\{f_{1},\ldots,f_{2K}\}$ be a set with cardinality $2K$ containing the indexes of the selected columns and define the $N\times 2K$ canonical matrix $\mathbf{E}_{\mathcal{F}}=[\mathbf{e}_{f_{1}},\ldots,\mathbf{e}_{f_{2K}}]$. Using this notation, the matrix containing the columns of $\mathbf{M}$ indexed by $\mathcal{F}$ is $\mathbf{M}\mathbf{E}_{\mathcal{F}}$, which can be alternatively written as

The expression reveals that $\mathbf{M}\mathbf{E}_{\mathcal{F}}$ is invertible because it can be written as the product of two $2K\times 2K$ invertible matrices. The latter is true because: a) conditions $\mathbf{C}=\mathbf{C}_{2K}(n_{0},N_{0})$, $\lambda_{k}^{N_{0}}\neq\lambda_{k^{\prime}}^{N_{0}}$ for all $k\neq k^{\prime}$, and $\lambda_{k}\neq 0$ for all $k$ guarantee that $(\mathbf{C}\boldsymbol{\Psi}\mathbf{E}_{\mathcal{F}})$ is invertible because it is a product of a diagonal and a full rank Vandermonde matrix (cf. Remark 1) and b) condition $[\boldsymbol{\upsilon}_{i}]_{k}\neq 0$ for all $k$ guarantees that $(\mathbf{E}_{\mathcal{F}}^{T}\text{diag}(\boldsymbol{\upsilon}_{i})\mathbf{E}_{\mathcal{F}})$ is an invertible diagonal matrix. This is true for any $\mathcal{F}$. The second step is to show that $2K$ observations are enough to guarantee identifiability. To see why this is the case, assume that two different feasible solutions $\widehat{\mathbf{x}}_{A}$ and $\widehat{\mathbf{x}}_{B}$ exist. This would imply that $\mathbf{M}(\widehat{\mathbf{x}}_{A}-\widehat{\mathbf{x}}_{B})=0$. Nevertheless, the vector $(\widehat{\mathbf{x}}_{A}-\widehat{\mathbf{x}}_{B})$ has, at most, $2K$ non-zero components and any choice of $2K$ columns of $\mathbf{M}$ generates a full rank square matrix which forces $\widehat{\mathbf{x}}_{A}=\widehat{\mathbf{x}}_{B}$, contradicting the assumption of multiple solutions. Finally, it is worth noting that although the proposition requires all the eigenvalues to be non-zero and distinct, only the ones associated with $\mathcal{K}$ need to satisfy those requirements. Note that the previous proof amounts to say that the matrix $\mathbf{C}\boldsymbol{\Psi}\text{diag}(\boldsymbol{\upsilon}_{i})$ has full spark and, hence, the claims in the proposition coincide with those in for the -norm recovery. ∎

It is worth stressing that the proof for the joint recovery and identification support leverages once more the fact that $\boldsymbol{\Psi}$ is a Vandermonde matrix, which is a distinct feature of the aggregation sampling scheme proposed in this paper. To gain more intuition about the result, we revisit the discussion provided after Proposition 1 and suppose that we know that $\mathbf{x}$ is a bandlimited signal with only one non-zero frequency component. This means that $K=1$ and that the graph signal can be written as $\mathbf{x}=\alpha\mathbf{v}_{k}$. If the value of $k$ is known, which amounts to say that the support where the signal is sparse is known, then node $i$ can interpolate the entire signal $\mathbf{x}$ using $x_{i}$ (cf. Section III-D). If the support is not known, Proposition 2 establishes $2K=2$ samples are needed. Clearly, one sample is not enough because $x_{i}=\hat{\alpha}[\mathbf{v}_{k}]_{i}$ admits $N$ different solutions, one per $k$. To see why two samples suffice, note that the first shift corresponds to $[\mathbf{y}_{i}]_{1}=x_{i}$ and the second to $[\mathbf{y}_{i}]_{2}=S_{ii}x_{i}+\sum_{j\in\mathcal{N}_{i}}S_{ij}x_{j}=\lambda_{\hat{k}}x_{i}$. Then, node $i$ can identify first the active frequency by finding the frequency index $\hat{k}$ whose associated eigenvalue is $[\mathbf{y}_{i}]_{2}/[\mathbf{y}_{i}]_{1}$. For the identification to succeed, the eigenvalues need to be distinct, as required by Proposition 2. Once the active frequency is known, the corresponding frequency coefficient can be estimated as before by setting $\hat{\alpha}=x_{i}/[\mathbf{v}_{\hat{k}}]_{i}$ and then the entire graph signal is $\hat{\mathbf{x}}=\hat{\alpha}\mathbf{v}_{\hat{k}}$. This discussion provides additional support to the idea that bandlimited graph signals can be understood as signals that can be inferred from local information.

From a computational perspective, the presence of the -norm in (35) renders the optimization non-convex, thus challenging to solve. A straightforward way to convexify it is to replace the -norm with a $1$-norm. Note that if such a process finds a feasible solution, call it $\widehat{\mathbf{x}}_{1}^{*}$, such that $||\widehat{\mathbf{x}}_{1}^{*}||_{0}=K$, then it holds that $\widehat{\mathbf{x}}^{*}=\widehat{\mathbf{x}}_{1}^{*}$. Conditions under which this process is guaranteed to identify the support can be found by analyzing the coherence and the restricted isometry property (RIP) of the matrix $\mathbf{C}\boldsymbol{\Psi}\text{diag}(\boldsymbol{\upsilon}_{i})$ . Unfortunately, determining the conditioning of all submatrices of a deterministic matrix (and, hence, the RIP) is challenging . The coherence of the matrix $\mathbf{C}\boldsymbol{\Psi}\text{diag}(\boldsymbol{\upsilon}_{i})$, denoted as $\mu_{i}(\mathbf{C})$, is easier to find and it depends on the most similar pair of eigenvalues of $\mathbf{S}$. However, the sparsity bound given by the matrix coherence, which requires K\leq\frac{1}{2}\big{(}1+\frac{1}{\mu_{i}(\mathbf{C})}\big{)} , is oftentimes too restrictive. A better alternative in that case is to use the concept of $t$-averaged mutual coherence and apply the results in for deterministic sensing matrices.

V-B Noisy joint recovery and support identification

If noise is present and the frequency support of the signal is unknown, the ($K$-sparse) least squares estimate of $\widehat{\mathbf{\mathbf{x}}}$ can be found as the solution to the following optimization problem

where the matrix multiplication $(\bar{\mathbf{R}}_{w}^{(i)})^{-1/2}$ in the objective accounts for the fact of the noise being colored. As in the noiseless case, a straightforward approach to convexify the problem is to replace the -norm with the $1$-norm and solve the problem \widehat{\mathbf{x}}_{1}^{*}:=\arg\min_{\widehat{\mathbf{x}}}\|(\bar{\mathbf{R}}_{w}^{(i)})^{-1/2}\big{(}\mathbf{C}\mathbf{y}_{i}-\mathbf{C}\boldsymbol{\Psi}\text{diag}(\boldsymbol{\upsilon}_{i})\widehat{\mathbf{x}}\big{)}\|_{2}^{2}+\gamma||\widehat{\mathbf{x}}||_{1} for different values of the parameter $\gamma$.

The challenges for support identification and the penalty paid in terms of error performance are related to those in the previous sections . If the conditioning of matrix $\mathbf{C}\boldsymbol{\Psi}\text{diag}(\boldsymbol{\upsilon}_{i})$ is poor, which depends heavily on how similar the eigenvalues in $\boldsymbol{\Lambda}$ are, the performance will be bad. Bounds can be found using the coherence of $\mathbf{C}\boldsymbol{\Psi}\text{diag}(\boldsymbol{\upsilon}_{i})$, which is tractable, or by analyzing its RIP. The results in for deterministic matrices can also be used here. An alternative to have performance guarantees in this case is to consider that the matrix $\mathbf{C}\boldsymbol{\Psi}\text{diag}(\boldsymbol{\upsilon}_{i})$ is random. This can be the case if, for example, $\mathbf{C}$ is designed as random or if there is noise in the application of the shift operator $\mathbf{S}$.

VI Space-shift sampling of graph signals

This section presents an alternative – more general – sampling setup that combines the selection sampling presented in Section III-B with the aggregation sampling proposed in Section III-C.

Based on this, $\underset{\bar{}}{\mathbf{z}}$ can be written as

where $\underset{\bar{}}{\mathbf{w}}$ is a vector of length $N^{2}$ obtained by concatenating the noise vectors $\mathbf{w}_{i}$ for all nodes $i$. This implies that (39) is a system of $N^{2}$ linear equations with $K<N$ variables. Thus, our objective is to pick $K$ of these equations in order to estimate $\widehat{\mathbf{x}}_{K}$ – and, hence, $\mathbf{x}$ through (18) – while minimizing the error introduced by the noise $\underset{\bar{}}{\mathbf{w}}$. Notice that if we restrict ourselves to pick $K$ equations out of the $N$ equations in positions $(i-1)\,N+1$ to $i\,N$ for some node $i\in\{1,2,\ldots,N\}$, then the problem reduces to local aggregation sampling at node $i$ as developed in Section III-C. Similarly, if we restrict ourselves to pick the $K$ equations out of the $N$ equations in positions $[1,1+N,1+2\,N,\ldots,1+(N-1)\,N]$, then the problem reduces to selection sampling as presented in Section III-B. In this sense, the formulation in (39) is more general. To implement the selection of the $K$ equations out of the $N^{2}$ options in (39), we use a binary selection matrix ${\mathbf{C}}$ as done in previous sections but, in this case, the size of $\mathbf{C}$ is $K\times N^{2}$. The reduced square system of linear equations can then be written as [cf. (14)]

The correlation matrices of the frequency $\widehat{\mathbf{R}}_{\underset{\bar{}}{e}}$ and time $\mathbf{R}_{\underset{\bar{}}{e}}$ errors of the estimator computed as the solution of (40) are [cf. (19) and (20)]

where $\mathbf{R}_{\underset{\bar{}}{w}}$ is the covariance matrix of the stacked vector of noise $\underset{\bar{}}{\mathbf{w}}$. In this aggregated case, the same noise models introduced in Section IV-B can be present. For the white noise in observations, we have that $\mathbf{R}_{\underset{\bar{}}{w}}=\sigma^{2}\mathbf{I}$, for the white noise in the original signal, we have that $\mathbf{R}_{\underset{\bar{}}{w}}=\sigma^{2}\left(\mathbf{I}\otimes\boldsymbol{\Psi}\right)\boldsymbol{\Upsilon}\boldsymbol{\Upsilon}^{H}\left(\mathbf{I}\otimes\boldsymbol{\Psi}\right)^{H}$, and for the white noise in the active frequency coefficients we have that $\mathbf{R}_{\underset{\bar{}}{w}}=\sigma^{2}\left(\mathbf{I}\otimes(\boldsymbol{\Psi}\mathbf{E}_{K})\right)\bar{\boldsymbol{\Upsilon}}\bar{\boldsymbol{\Upsilon}}^{H}\left(\mathbf{I}\otimes(\boldsymbol{\Psi}\mathbf{E}_{K})\right)^{H}$.

In the previous discussion, no structure was assumed in the selection matrix $\mathbf{C}$. A case of particular interest is when the sampling schemes are implemented in a distributed manner using message passing. Suppose that the sampling is performed at node $i$. To compute $y_{i}^{(l)}$, the node $i$ needs to have access to $y_{j}^{(l^{\prime})}$ for all $j\in\mathcal{N}_{i}$ and $l^{\prime}<l$. To simplify notation, and without loss of generality, we will assume that the sampling node is $i=1$ and that the neighbors of $i\!=\!1$ are $i=2,\ldots,N_{1}+1$. Suppose also that node $i\!=\!1$ computes $L_{1}$ shifts, from $y_{1}^{(0)}$ up to $y_{1}^{(L_{1})}$. This implies that node $i=1$ has access to $L_{1}+1$ of its own samples and to $L_{1}$ samples of each of its $N_{1}$ neighbors. The selection matrix $\mathbf{C}$ can be written as

Matrix $\mathbf{C}$ has $1+L_{1}(1+N_{1})$ rows, one per observation. The first $1+L_{1}$ rows correspond to the samples at node $i\!=\!1$ and the remaining $L_{1}N_{1}$ to the samples at its neighbors. Note also that matrix $\mathbf{C}\left(\mathbf{I}\otimes(\boldsymbol{\Psi}\mathbf{E}_{K})\right)\bar{\boldsymbol{\Upsilon}}$ is not full (row) rank. The reason is that all the samples obtained at node $i\!=\!1$, except for the first one, are linear combinations of the samples at its neighbors. This implies that the number of frequencies that can be recovered using (43) is, at most, $1+L_{1}N_{1}$.

Structured observation models different from the one in (43) can be also of interest. For example, one can consider setups where nodes from different parts of the graph take a few samples each and forward those samples to a central fusion center. In such a case, since the nodes gathering data need not be neighbors, the problem associated with some of the samples being a linear combination of the others will not necessarily be present.

VII Numerical experiments

We start by illustrating the perfect recovery of synthetic noiseless graph signals, both when the frequency support is known and when it is not (Section VII-A). We then present results for real-world graph signals corresponding to the exchange among the different sectors of the economy of the United States. These are used to test recovery under the presence of noise (Section VII-B) as well as to illustrate the space-shift sampling method (Section VII-C).

Consider the 20-node undirected graph $\mathcal{G}$ depicted in Fig. 3(a), which corresponds to a realization of a symmetric Erdõs-–Rényi graph with edge probability $0.20$ . With $\mathbf{A}=\mathbf{V}\boldsymbol{\Lambda}_{A}\mathbf{V}^{H}$ denoting the adjacency matrix of $\mathcal{G}$, three different graph-shift operators are considered: $\mathbf{S}_{1}=\mathbf{A}$, $\mathbf{S}_{2}=\mathbf{I}-\mathbf{A}$, and $\mathbf{S}_{3}=0.5\mathbf{A}^{2}$. Notice that, even though the support of $\mathbf{S}_{3}$ differs from that of $\mathbf{S}_{1}$ and $\mathbf{S}_{2}$, the graph-shift operator $\mathbf{S}_{3}$ still preserves the notion of locality as defined by a two-hop neighborhood. Note also that the three shift operators share the same set of eigenvectors $\mathbf{V}$, but they have a different set of eigenvalues.

Let $\mathbf{x}$ be a graph signal supported on $\mathcal{G}$. This signal is represented in Fig. 3(a). To facilitate interpretation, the value of the signal at a given node is written explicitly next to the node and also coded by the color of the node. Although seemingly random in the node domain, the structure of the signal $\mathbf{x}$ is highly determined by the graph. To illustrate this, Fig. 3(b) presents the frequency components $\widehat{\mathbf{x}}$ of signal $\mathbf{x}$, where the graph frequency basis are given by the columns of $\mathbf{V}$. The figure reveals that $\mathbf{x}$ has a bandwidth of $K=3$. Since $\mathbf{V}$ is the basis for $\mathbf{S}_{1}$, $\mathbf{S}_{2}$ and $\mathbf{S}_{3}$, the frequency representation $\widehat{\mathbf{x}}$ and the bandwidth $K$ are the same for any of the three operators. As a result, the procedure described in Section III-C allows recovering the whole signal using three aggregated samples, no mater which operator is chosen for the aggregation.

Suppose that we select node $i=4$ as sampling node, which is circled in red in Fig. 3(a). If the shift is $\mathbf{S}_{1}$, the $3$ first observations taken by that node are $\mathbf{y}_{4}=[-0.55,1.27,-2.94]^{T}$. The first observation corresponds to the value of the signal at node 4, the second one to the aggregated signal at its neighbors and the third observation corresponds to a linear combination of the signal values within its two-hop neighborhood. Since $K=3$, Proposition 1 guarantees recovery if: i) the $3$ first eigenvalues of the shift operator are distinct and ii) the $3$ first values of $\boldsymbol{\upsilon}_{4}$ are non-zero. It turns out that for $\mathbf{S}_{1}$ and node 4 these two conditions hold true and, hence, the interpolation in (13) yields the original signal in Fig. 3(a). In fact, for the network at hand, these two conditions are satisfied for all nodes and shift operators considered. This implies that perfect reconstruction in a noiseless setting is achieved independently of which node aggregates the information and which shift operator – among the three presented – is picked. To better asses the conditions in Proposition 1, we build 10,000 different random graphs where the edge probability is randomly chosen from the interval $[0.15,0.25]$. Realizations that do not give rise to a connected graph are discarded. We vary the number of nodes from 10 to 30 and the active frequencies of the graph signals from 1 to 5. For each random graph and signal defined on it, we test for perfect signal recovery on every node. The simulations show that in $99.89\%$ of the cases the signal is successfully recovered.

The graph signal $\mathbf{x}$ can be sampled and recovered even when the frequency support is unknown, i.e., when we know that $\widehat{\mathbf{x}}=\mathbf{V}^{-1}\mathbf{x}$ contains $K=3$ non-zero components, but we do not know the indices of these $K$ active frequencies. In this case, however, $2K=6$ samples are needed to guarantee identifiability (cf. Proposition 2). By solving problem (35), the signal can be recovered at every node and using any of the three shift operators, as in the previous case. However, when solving a relaxed version of problem (35), accurate signal recovery depends on the specific network, signal and node selected for reconstruction. Moreover, the recovery rate depends on the choice of the graph-shift operator $\mathbf{S}$. For example, for the signal in Fig. 3(a), solving a 1-norm relaxation of the problem (35) yields the original graph signal $\mathbf{x}$ if $\mathbf{S}=\mathbf{A}$ and $i=4$, but fails if $\mathbf{S}=\mathbf{I}-\mathbf{A}$ and $i=5$ where node $i=5$ is the right neighbor of node $i=4$. To assess recovery better, Fig. 5 plots the success rate – fraction of realizations for which the actual signal was recovered – for graph-shifts $\mathbf{S}_{1}$, $\mathbf{S}_{2}$ and $\mathbf{S}_{3}$, and different number of observations. Each point in the plots represents an average across all nodes in the network, 5 signal realizations and 10 random graph realizations. The three plots correspond to symmetric Erdõs-–Rényi graphs generated using different edge probabilities: 0.15, 0.20, and 0.25. The recovery rate for $\mathbf{S}_{3}=0.5\mathbf{A}^{2}$ is consistently higher than for the other shift operators considered. This is not surprising: when squaring the adjacency matrix to generate $\mathbf{S}_{3}$, the dissimilarity between any pair of eigenvalues is increased, which reduces the matrix coherence $\mu_{i}(\mathbf{C})$ associated with $\mathbf{S}_{3}=0.5\mathbf{A}^{2}$ and facilitates sparse recovery (cf. last paragraph in Section V-A). Nonetheless, if success rate is the main concern, there exist relaxations of the 0-norm that give better results than the 1-norm used .

VII-B Recovery in the presence of noise

In Sections VII-B1 to VII-B3 we consider the bandlimited signal $\mathbf{x}_{4}$ as noiseless and add different types of Gaussian noise to analyze the interpolation performance at different nodes. Differently, in Section VII-B4 we interpret the whole graph signal $\mathbf{x}$ as a noisy version of $\mathbf{x}_{4}$ and analyze the reconstruction error when interpolating $\mathbf{x}$ from just 4 samples.

We perform aggregation sampling of multiple noisy versions of $\mathbf{x}_{4}$ via successive applications of the graph-shift $\mathbf{S}$ at different economic sectors (nodes). The noisy versions of $\mathbf{x}_{4}$ are generated by adding noise to the observed signal as described in (25). The power of the white noise $\sigma^{2}$ is the same for all nodes and is computed so that, averaging across nodes, the linear signal to noise ratio (SNR) for the first, second, third and fourth observations in each node is 2, 10, 50, and 250, respectively. This increase in SNR is attributable to the fact that successive applications of the shift $\mathbf{S}$ increase the signal magnitude while we keep the noise power $\sigma^{2}$ constant. In Fig. 6(b) we plot the empirical average reconstruction error at different nodes across 1,000 noisy realizations of $\mathbf{x}_{4}$ and compare it with the theoretical average error, i.e., the trace of $\mathbf{R}_{e}^{(i)}$ in (20) [cf. (21)]. We first observe that the computed theoretical error indeed coincides with the average empirical error across realizations. Moreover, notice that the reconstruction performance is highly node dependent. The error is minimized for the reconstruction based on the artificial sectors AV and FU. This is reasonable since these two nodes – unlike other sectors – are closely related to every other sector of the economy (cf. Fig. 4). Furthermore, the sectors achieving the worst reconstruction errors are ‘Publishing Industries’ and ‘Ground Passenger Transportation’ corresponding to nodes 34 and 31. This can be explained by analyzing the vectors $\bar{\boldsymbol{\upsilon}}_{34}=\boldsymbol{\upsilon}_{34}\mathbf{E}_{4}$ and $\bar{\boldsymbol{\upsilon}}_{31}$ (cf. Lemma 1). Even though both vectors have all four components different from zero, which guarantees perfect reconstruction in the noiseless case (cf. Proposition 1), they possess an element whose absolute value is in the order of $10^{-4}$, increasing the sensitivity of the reconstruction in the presence of noise. For all other nodes the smallest element of $\bar{\boldsymbol{\upsilon}}_{i}$ is at least one order of magnitude larger. Fig. 6(d) presents the reconstruction obtained by aggregation sampling in node 46 corresponding to ‘Professional Services’ – best among real economic sectors, i.e., excluding AV and FU – which achieves an error of 0.26.

VII-B2 White noise in the original signal

Similarly to the analysis performed in the previous section, we investigate the reconstruction performance of aggregation sampling at different nodes. However, in this case, the noise is added to the original signal, following the model described in (26). The power of the white noise $\sigma^{2}$ is set to induce a linear SNR of $10^{2}$. As was the case in the previous section, the average empirical error (across 1,000 realizations) matches closely our theoretical estimates; see Fig. 6(c). Moreover, the specific nodes that lead to a good (bad) interpolation performance are very similar to those in the previous noise model. Indeed, sectors 34 and 31 have the highest reconstruction error whereas AV and FU attain the best reconstructions. Fig. 6(d) shows the best reconstruction – excluding AV and FU – which amounts to an error of 0.001 and corresponds to the sector ‘Professional Services’ at node 46.

VII-B3 White noise in the active frequencies

We consider a third category of noisy versions of $\mathbf{x}_{4}$ obtained by adding white noise only to the four active frequencies, as described in (27). The power of the white noise $\sigma^{2}$ is set to induce a linear SNR of $10^{2}$. The empirical reconstruction error associated with each node – averaged over 1,000 noisy realizations of $\mathbf{x}_{4}$ – is the same among nodes. This validates the analysis in (28), which stated that, for this noise model, the quality of the reconstruction is node independent. In Fig. 6(d) we present an example of such reconstruction, achieving an error of 0.01.

VII-B4 Real-world noisy signal

We interpret the graph signal $\mathbf{x}$ as a noisy realization of a signal of bandwidth 4. Hence, our goal is to obtain the best reconstruction of $\mathbf{x}$ based on 4 observations. As described in (20) and shown before, interpolation performance is highly node dependent. Indeed, the reconstruction error when keeping the first 4 observations at each node spans 5 orders of magnitude depending on the sampling node, although for most nodes it is contained between $10^{-3}$ and $10^{-1}$; see Fig. 6(e). The best reconstruction among the real sectors is achieved by ‘Insurance Carriers’ (node 40). The best and the median reconstructions are acceptable, attaining errors of 0.0035 and 0.019, respectively. Fig. 6(f) depicts the best reconstruction.

VII-C Space-shift sampling

In Section VII-B4 we analyzed the accuracy of interpolating the U.S. economic activity after aggregation sampling in different economic sectors. The minimum and median reconstruction errors are presented in the first row of Table I, where the reconstruction error is quantified as the ratio between the energy of the error and that of the original signal. An alternative approach is to implement selection sampling, i.e. to sample the signal $\mathbf{x}$ in 4 different sectors – excluding the artificial sectors AV and FU – and interpolate the whole signal from these 4 observations, as explained in Section III-B. Recall that reconstruction is not guaranteed for every subset of 4 nodes since we must have invertibility of $(\mathbf{C}\mathbf{V}_{K})$ [cf. (6)]. By analyzing the minimum and median reconstruction errors – see the two first rows in Table I – it is clear that the node aggregation sampling outperforms the node selection sampling. This is intuitive since most of the energy of the signal is contained in the two first frequencies [cf. Fig. 6(a)(top)], which are associated with the largest eigenvalues. Hence, after successive implementations of the graph-shift, the error in estimating these frequencies is reduced, resulting in a smaller error in the interpolation of the whole signal.

As developed in Section VI, more general sampling strategies can be implemented. For example, we can sample the value of the signal at 4 nodes after the application of one, two or three graph-shifts. The results – listed in rows 3, 4 and 5 of Table I– reveal that reduction in the median error after each graph-shift application is conspicuous, especially when going from no applications – median error of 4.2 – to one application – median error of 0.03. A different alternative is a sampling strategy that selects the original signal and the signal after one shift in two different sectors. The results, listed in the last row of Table I, show that this configuration leads to a very good reconstruction performance: 0.0035 minimum error and 0.039 median error. Note that with this sampling configuration, the two sectors are only required to compute the aggregated activity of their one-hop neighbors.

VIII Conclusions

A novel scheme for sampling bandlimited graph signals – that admit a sparse representation in the frequency domain – was proposed. The scheme was based on the aggregation of local information at a single node after successive applications of the graph-shift operator. This contrasted to most existing works, which focus on sampling the value of the signal observed at a subset of nodes. Our scheme was shown to be equivalent to classical sampling for directed cycle graphs whereas, for more general graphs, the Vandermonde structure of the sampling matrix was exploited to determine the conditions for perfect reconstruction in the absence of noise. Reconstruction under correlated noise was analyzed, and design criteria to pick the sampling node and shifts leading to optimal noisy reconstruction were discussed. Scenarios where the specific set of frequencies present in the bandlimited signal is not known were also investigated and connections with sparse signal reconstruction were drawn. Finally, a more general sampling scheme was presented which contained, as particular cases, the selection sampling as well as our local aggregation approach. The various sampling and interpolation scenarios were illustrated through numerical experiments in both synthetic and real-world graph signals.

References