Signal Recovery on Graphs: Fundamental Limits of Sampling Strategies

I Introduction

The massive amounts of data being generated from various sources, including online social networks, citation networks, biological networks, and physical infrastructures has inspired an emerging field of research for analyzing data supported on graphs . Signal processing on graphs is a theoretical framework for the analysis of high-dimensional data with complex, nonregular structures ; it extends classical discrete signal processing to signals with such structure, models it by a graph and signals by graph signals, generalizing concepts and tools from classical discrete signal processing to graph signal processing. Recent work involves sampling of graph signals , recovery of graph signals , representations for graph signals , uncertainty principles on graphs , graph dictionary construction , semi-supervised learning with graphs , graph denoising , community detection and clustering on graphs , graph-based filter banks and graph-based transforms , among others.

In this paper, we consider the classical signal processing task of sampling and recovery . As the bridge connecting sequences and functions, classical sampling theory shows that a bandlimited function can be perfectly recovered from its sampled sequence if the sampling rate is high enough. In the 60 years since Shannon, many new regular and irregular sampling and recovery frameworks have emerged to recover signals with different properties. For example, finite rate of innovation sampling considers sampling and recovery of signals that have a finite number of degrees of freedom per unit time , compressed sensing considers sampling and recovery of sparse signals or signals that can be sparsely represented by some coherent and redundant dictionary and subspace sampling considers sampling and recovery of signals from a union of subspaces .

The interest in sampling and recovery of graph signals has increased in the last few years . In , authors proposed an algorithm to recover graph signals that have small variation under uniform sampling, with an upper bound on recovery error. In , authors proposed a sampling theory for graph signals and showed perfect recovery for bandlimited graph signals under experimentally designed sampling. Extensions include a fast distributed algorithm , local weighted measurements , local aggregation and percolation from seeding nodes .

In this paper, we focus on smooth graph signals, that is, signals whose coefficients at each node are similar to the coefficients of its neighbors. We propose a new class of smooth graph signals, called approximately bandlimited and build a theoretical foundation to understand sampling and recovery of this class under uniform sampling, experimentally designed sampling and active sampling. We propose a recovery strategy to compare uniform sampling and experimentally designed sampling; this recovery strategy is an unbiased estimator for low-frequency components and achieves the optimal rate of convergence under some assumptions on the graph structures. In spirit, our work follows previous work that studied the theoretical capabilities of passive and active sampling for recovering functions from samples ; the difference is that we consider a discrete setting and deal with irregular structures. For a smooth function, active sampling, experimentally designed sampling and uniform sampling have the same performance from a statistical perspective . For approximately bandlimited graph signals, however, we will see that while active sampling achieves the same rate of convergence as experimentally designed sampling, experimentally designed sampling fundamentally outperforms uniform sampling when the graph is irregular.

To validate the recovery strategy, we test it on five specific graphs: a ring graph with $k$ nearest neighbors, an Erdős-Rényi graph, a random geometric graph, a small-world graph and a power-law graph, and show that experimental results match the proposed theory well.

Contributions. The contributions of the paper are as follows: We propose

a new class of smooth graph signals related to existing classes of smooth graph signals;

minimax lower bounds on the recovery error under three sampling strategies;

a recovery strategy that achieves optimal rates of convergence on two specific types of graphs;

a statistical analysis of graph structures; and

a comprehensive explanation of when and why experimentally designed sampling works for semi-supervised learning with graphs.

Outline of the paper. Section II briefly reviews graph signal processing; Section III reviews smooth graph signal models and formulates sampling and recovery strategies. We propose the minimax lower bounds on the recovery errors in Section IV and a recovery strategy that works for both uniform sampling and experimentally designed sampling in Section V. Section VI combines the results from the two previous sections and shows the optimal convergence rates of recovery on two types of graphs. The proposed recovery strategy is evaluated in Section VII on five known graph classes. Section VIII concludes the paper and provides pointers to future directions.

II Signal Processing on Graphs

The Jordan decomposition of $\operatorname{A}$ is

and the inverse graph Fourier transform is

where $\mathbf{v}_{k}$ is the $k$th column of $\operatorname{V}$ and $\widehat{x}_{k}$ is the $k$th component in $\widehat{\mathbf{x}}$. The vector $\widehat{\mathbf{x}}$ in (2) represents the signal’s expansion in the eigenvector basis and describes the frequency components of the graph signal $\mathbf{x}$. The inverse graph Fourier transform reconstructs the graph signal by combining graph frequency components. In general, $\operatorname{V}$ is not orthonormal; to restrict its behavior, we assume that

where $\alpha_{1},\alpha_{2}>0$, that is, $\operatorname{V}$ is a Riesz basis with stability constants $\alpha_{1},\alpha_{2}$ . When $\operatorname{A}$ represents an undirected graph, it is symmetric, we have $\operatorname{U}=\operatorname{V}^{T}$ and both $\operatorname{U}$ and $\operatorname{V}$ are orthonormal. For simplicity, we consider $\operatorname{A}$ to represent an undirected graph in this paper, that is, we assume $\alpha_{1}=\alpha_{2}=1$, but the proposed graph signal models and sampling strategies apply equally well for directed graphs. Note that there exist other versions of the graph Fourier transform based on different approaches to normalize the adjacency matrix, such as the eigenvector matrix of the graph Laplacian matrix ; our proposed methods work for all the versions of the graph Fourier transforms. Table I lists the notations used in the paper.

III Problem Formulation

We now review three classes of smooth graph signals and show connections between them. We then describe the sampling and recovery strategies of interest, connecting this work to the previous work on sampling theory on graphs.

We consider a graph signal to be smooth when coefficients at each node are similar to the coefficients of its neighbors. We have previously defined two classes of smooth graph signals in .

Denote this class of graph signals by $\operatorname{GS}_{\operatorname{A}}(\eta)$.

In the above, $\operatorname{A}\mathbf{x}$ is the shifted version of $\mathbf{x}$ and $\mathbf{x}-\operatorname{A}\mathbf{x}$ gives the first-order difference . Since we normalized the graph shift such that $|\lambda_{\max}(\operatorname{A})|=1$,

Thus, when $\eta\geq 4$, all graph signals satisfy (3).

Denote this class of graph signals by $\operatorname{BL}_{\operatorname{A}}(K)$.

Note that the original definition in just requires $\widehat{\mathbf{x}}$ to be $K$-sparse, meaning that $\widehat{\mathbf{x}}$ is not necessarily lowpass. Here, instead, we fix the ordering of frequencies and requires the bandlimited graph signals to be lowpass. The following theorem details the relationship between these two classes.

For any $K\in\{0,1,\cdots,N-1\}$, $\operatorname{BL}_{\operatorname{A}}(K)$ is a subset of $\operatorname{GS}_{\operatorname{A}}(\eta)$, when $\eta\geq(1-\lambda_{K-1})^{2}$.

To show when $\operatorname{BL}_{\operatorname{A}}(K)\subseteq\operatorname{GS}_{\operatorname{A}}(\eta)$, let $\mathbf{x}\in\operatorname{BL}_{\operatorname{A}}(K)$, that is, $\mathbf{x}=\sum_{k=0}^{K-1}\widehat{x}_{k}\mathbf{v}_{k}.$ Then,

where $\mathcal{N}_{i}$ denotes the neighbors of the $i$th node (basically all $j$ for which $\operatorname{A}_{i,j}\neq 0$), (a) follows from the fact that $\mathbf{v}_{k}$ and $\lambda_{k}$ are the $k$th eigenvector and eigenvalue of $\operatorname{A}$, respectively and thus $(\operatorname{A}\mathbf{v}_{k})_{i}=\lambda_{k}(\mathbf{v}_{k})_{i}$; and (b) from the ordering of eigenvalues $\lambda_{k}\geq\lambda_{k-1}$. From the above, we see that for bandlimited graph signals, the signal coefficient at each node is close to the weighted average of all its neighbors; in other words, bandlimited graph signals are smooth locally, which implies global smoothness.

proving that the bandlimited graph signals form a subset of globally smooth graph signals. The converse is not true, however, as globally smooth graph signals can have isolated high-frequency components and are thus not bandlimited. ∎

While the recovery of globally smooth graph signals has been studied in (leading to graph signal inpainting), the criterion of global smoothness is quite loose, making it hard to provide further theoretical insight . Similarly, while the recovery of bandlimited graph signals has been studied in (leading to sampling theory on graphs), the requirement of bandlimitedness is rather restrictive, making it impractical in real-world applications.

We thus propose a third class of smooth graph signals that relaxes the requirement of bandlimitedness, but still promotes smoothness We proposed the approximately bandlimited class in ..

Denote this class of graph signals by $\operatorname{ABL}_{\operatorname{A}}(K,\beta,\mu)$.

The approximately bandlimited class allows for a tail after the first $K$ frequency components, whose shape and decay are controlled by $\mu$ and $\beta$; the smaller the $\mu$, the less energy from the high-frequency components is allowed in the tail, and the larger the $\beta$, the higher the penalty on the energy from those high-frequency components. The class of $\operatorname{ABL}_{\operatorname{A}}(K)$ is similar to the ellipsoid constraints in previous literature , where all the frequency components are considered in the constraints; in other words, $\operatorname{ABL}_{\operatorname{A}}(K)$ poses fewer restrictions on the low-frequency components. Many real graph signals exhibit the approximately bandlimited property; for example, Figures 1 and 2 show that the temperature readings across the U.S and wind speeds across Minnesota are approximately bandlimited graph signals. We have found that the approximately bandlimited class is more powerful than the bandlimited class when representing real graph signals.

The following theorem details the relationship between $\operatorname{ABL}_{\operatorname{A}}(K,\beta,\mu)$ and $\operatorname{GS}_{\operatorname{A}}(\eta)$.

With $\beta\geq 1$, $\mu,\eta\geq 0$ and $K\in\{0,1,\cdots,N-1\}$,

$\operatorname{ABL}_{\operatorname{A}}(K,\beta,\mu)$ is a subset of $\operatorname{GS}_{\operatorname{A}}(\eta)$, when

$\operatorname{GS}_{\operatorname{A}}(\eta)$ is a subset of $\operatorname{ABL}_{\operatorname{A}}(K,\beta,\mu)$, when

From Theorem 2, we see that depending on the parameters, $\operatorname{GS}_{\operatorname{A}}(\eta)$ can be a subset of $\operatorname{ABL}_{\operatorname{A}}(K,\beta,\mu)$ or vice versa.

To show when $\operatorname{ABL}_{\operatorname{A}}(K,\beta,\mu)\subseteq\operatorname{GS}_{\operatorname{A}}(\eta)$, let $\mathbf{x}\in\operatorname{ABL}_{\operatorname{A}}(K,\beta,\mu)$. Then,

where (a) follows from the triangle inequality; (b) from Theorem 1; (c) from $\operatorname{A}\mathbf{v}_{k}=\lambda_{k}\mathbf{v}_{k}$; (d) from the fact that $\left\|\operatorname{V}\alpha\right\|_{2}=\|\alpha\|_{2}$ by Parseval’s equality since $\operatorname{V}$ is orthonormal, with $\alpha_{i}=(1-\lambda_{i})\widehat{x}_{i}$, for $i=K,\,K+1,\,\ldots,\,N-1$ and otherwise; (e) from $-1\leq\lambda_{K}\leq 1$; and (f) from $\mathbf{x}\in\operatorname{ABL}_{\operatorname{A}}(K,\beta,\mu)$.

To show the second statement when $\operatorname{GS}_{\operatorname{A}}(\eta)\subseteq\operatorname{ABL}_{\operatorname{A}}(K,\beta,\mu)$, let $\mathbf{x}\in\operatorname{GS}_{\operatorname{A}}(\eta)$. Then,

where (a) follows from $\sum_{k=K}^{N-1}(1-\lambda_{k})^{2}\widehat{x}_{k}^{2}\leq\sum_{k=0}^{N-1}(1-\lambda_{k})^{2}\widehat{x}_{k}^{2}=\left\|\mathbf{x}-\operatorname{A}\mathbf{x}\right\|_{2}^{2}\leq\eta\left\|\mathbf{x}\right\|_{2}^{2}$, where the last inequality follows from Definition 1. ∎

From Theorem 2, we see that $\operatorname{ABL}_{\operatorname{A}}(K,\beta,\mu)$ is not only more general than $\operatorname{BL}_{\operatorname{A}}(K)$, but describes $\operatorname{GS}_{\operatorname{A}}(\eta)$ in a more controlled fashion. In this paper, we focus on $\operatorname{ABL}_{\operatorname{A}}(K,\beta,\mu)$, and study the recovery performance of this class under three sampling strategies.

III-B Sampling & Recovery Strategies

We consider three different sampling strategies:

uniform sampling means that sample indices are chosen from $\{0,1,\ldots,N-1\}$ independently and randomly;

experimentally designed sampling means that sample indices can be chosen beforehand; and

active sampling means that sample indices can be chosen as a function of the sample points and the samples collected up to that instance, that is, $\mathcal{M}_{i}$ depends only on $\{\mathcal{M}_{j},y_{j}\}_{j<i}$.

It is clear that uniform sampling is a subset of experimentally designed sampling, which is again a subset of active sampling.

The goal in this paper is to study the fundamental limitations of these three sampling strategies when recovering approximately bandlimited graph signals. This study is related to many real-world applications. For example, in semi-supervised learning, datasets are modeled as a graph with data samples as nodes and similarities between those data samples as edges. Features and labels associated with data samples form approximately bandlimited graph signals. We aim to select data samples as labeled data and recover the labels for the unlabeled data. The sampling strategy helps select the most informative data samples and minimizes the recovery error. Some other applications include route planning based on wind speed , sensor position selection and compressive spectral clustering .

IV Fundamental Limits of Sampling Strategies

In this section, we study the fundamental limitations of the three sampling strategies for recovering $\operatorname{ABL}_{\operatorname{A}}(K,\beta,\mu)$ by showing minimax lower bounds. We do this by following the minimax decision rule and finding tight lower bounds for the minimax risk over all possible recovery strategies . In other words, we try to minimize the recovery error in the worst case and use a tight lower bound to describe this minimax error.

The lower bounds we present will be of the form

The upper bounds on the maximum risk are usually obtained through explicit recovery strategies, as will be shown in Section V. The upper bounds we present will be of the form

where $C>0$ is a constant. If (9) and (10) both hold, then $\phi_{m}$ is said to be the optimal rate of convergence and there is no recovery strategy that asymptotically outperforms the proposed one. When talking about optimal rates of convergence, we bound the sequence by a polynomial and make statements of the form: Given $\gamma_{1}<\gamma<\gamma_{2}$, a rate of convergence $\phi_{m}^{2}$ is equivalent to $m^{-\gamma}$ if and only if $m^{-\gamma_{2}}<\phi_{m}^{2}<m^{-\gamma_{1}}$ for $n$ large enough.

Note that the general bounds are for arbitrary graphs and thus involve parameters that depend on the graph structure; given the graph structure, we can specify the parameters and show explicit rates of convergence, as we will do in Section VI.

Let $\operatorname{V}_{(2,K)}$ be the sub-matrix of $\operatorname{V}$, consisting of the $K$th to the $(2K-1)$th columns of $\operatorname{V}$.

Given the class $\operatorname{ABL}_{\operatorname{A}}(K,\beta,\mu)$: (1) Under uniform sampling,

where $c_{1}>0$ , $0<c<1$, and $\Theta_{\rm u}$ denotes the set of all recovery strategies based on uniform sampling; (2) Under experimentally designed sampling,

where $c_{1}>0$ , $0<c<1$, and $\Theta_{\rm e}$ denotes the set of all recovery strategies based on experimentally designed sampling; (3) Under active sampling,

where $c_{1}>0$ , $0<c<1$, and $\Theta_{\rm a}$ denotes the set of all recovery strategies based on active sampling.

See Appendix A in the supporting document for the proof of these results. From Theorem 3, we see that experimentally designed sampling has the same minimax lower bound as active sampling, which means that collecting the feedback before taking samples does not fundamentally improve the recovery performance. We also see that the three minimax lower bounds depend on the properties of $\operatorname{V}_{(2,\kappa_{0})}$, which depend on the graph structure. When each row of $\operatorname{V}_{(2,\kappa_{0})}$ has roughly similar energy, $\left\|\operatorname{V}_{(2,\kappa_{0})}\right\|_{F}^{2}$ and $N\left\|\operatorname{V}_{(2,\kappa_{0})}\right\|_{\infty,2}^{2}$ are similar; when the energy is concentrated in a few rows, $N\left\|\operatorname{V}_{(2,\kappa_{0})}\right\|_{\infty,2}^{2}$ is much larger than $\left\|\operatorname{V}_{(2,\kappa_{0})}\right\|_{F}^{2}$; in other words, the minimax lower bound for experimentally designed sampling is tighter than that for uniform sampling. This happens in many real-world graphs that have complex, irregular structure. The minimax lower bounds thus show the potential advantage of experimentally designed sampling and active sampling over uniform sampling. We will elaborate on this in Sections VI and VII.

V Recovery Strategy

In the previous section, we presented the minimax lower bounds for each of the three sampling strategies and showed that active sampling cannot fundamentally perform better than experimentally designed sampling. We now propose a recovery strategy for both uniform sampling and experimentally designed sampling and evaluate its statistical properties.

To analyze uniform sampling and experimentally designed sampling in a similar manner, we consider sampling score-based sampling that unifies both sampling strategies. Sampling score-based sampling means that the sample indices are chosen from an importance sampling distribution that is proportional to some sampling score. Let $\{\pi_{i}\}_{i=1}^{N}$ be a set of sampling scores, where $\pi_{i}$ denotes the probability to choose the $i$th sample in each random trial. When the sampling score for each node is the same, we get uniform sampling; when the sampling score for each node is designed based on the graph structure, we get experimentally designed sampling.

For $\operatorname{ABL}_{\operatorname{A}}(K)$, most of the energy is concentrated in the first $K$ frequency components and the graph signal can be approximately recovered by using those. We consider the following recovery strategy to estimate those frequency components by projecting samples onto the bandlimited space spanned by the first $K$ frequency components.

where $\mathbf{y}$ is a vector representation of the samples (7).

It is also intuitive to consider the following recovery strategy ,

where we fit the sampled elements of the recovered graph signal to the samples by solving the least squares problem and $(\cdot)^{\dagger}$ is the pseudo-inverse operator . The corresponding expected recovered graph signal is

where $\operatorname{P}=\left(\operatorname{D}^{2}\Psi^{T}\Psi\operatorname{V}_{(K)}\right)^{\dagger}\operatorname{D}$ $\Psi^{T}\Psi=\left(\operatorname{U}_{(K)}\Psi^{T}\Psi\operatorname{D}^{2}\Psi^{T}\Psi\operatorname{V}_{(K)}\right)^{-1}$ $\operatorname{U}_{(K)}\Psi^{T}\Psi\operatorname{D}^{2}\Psi^{T}\Psi$, $\operatorname{V}_{(-K)}$ chooses the last $K$ columns of $\operatorname{V}$, and $\widehat{\mathbf{x}}_{(-K)}$ chooses the last $K$ columns of $\mathbf{x}$. When $\mathbf{x}$ is bandlimited, $\mathbf{x}^{*}_{\rm LS}$ perfectly recovers $\mathbf{x}$; when $\mathbf{x}$ is approximately bandlimited, however, $\mathbf{x}^{*}_{\rm LS}$ is a mixture of low-frequency and high-frequency components and it is hard to show its statistical properties. Moreover, it is less computationally efficient to compute $\mathbf{x}^{*}_{\rm LS}$ than $\mathbf{x}^{*}_{\rm SP}$ because of the presence of the inverse term.

Since the definitions of the bandlimited class and the approximately bandlimited class, and the proposed sampling strategies are all based on the graph Fourier transform instead of on graph shift, all the proposed methods work for other versions of the graph Fourier transform as well.

V-B Statistical Analysis

We study the statistical properties of Algorithm 1 by providing the bias, covariance, mean square error (MSE), and an optimal sampling distribution.

The following lemma shows that the sampled projection estimator is an unbiased estimator of the first $K$ frequency components for any sampling scores.

The sampled projection estimator with bandwidth $K$ and arbitrary sampling scores is an unbiased estimator of the first $K$ frequency components, that is,

where $\mathbf{x}^{*}$ is the solution of Algorithm 1.

Lemma 2 gives the exact covariance of the sampled projection estimator.

The covariance of sampled projection estimator $\mathbf{x}^{*}$ has the following property:

where Tr($\cdot$) is the trace operator and $\operatorname{W}_{C}$ is a diagonal matrix with $\left(\operatorname{W}_{C}\right)_{i,i}=(x_{i}^{2}+\sigma^{2})/(m\pi_{i})$.

Theorem 4 shows the exact MSE of sampled projection estimator and an upper bound.

For $\mathbf{x}\in\operatorname{ABL}_{\operatorname{A}}(K,\beta,\mu)$, let $\mathbf{x}^{*}$ be the sampled projection estimator with bandwidth $\kappa\geq K$. Then,

We merge the proofs of Lemmas 1, 2 and Theorem 4 in Appendix B in the supporting document. The main idea follows from the bias-variance tradeoff. The bias term $\mu/(1+\kappa^{2\beta})\left\|\mathbf{x}\right\|_{2}^{2}$ comes from the high-frequency components and ${\rm Tr}\left(\operatorname{U}_{(\kappa)}\operatorname{W}_{C}\operatorname{V}_{(\kappa)}\right)$ comes from the covariance. The last inequality in Theorem 4 comes from relaxing the bias term and omitting a constant, which has nothing to do with the sampling scores. Thus, optimizing the upper bound over sampling scores is equivalent to optimizing the exact MSE.

We next study the MSEs of sampled projection estimator based on both uniform sampling and experimentally designed sampling.

The upper bound of MSE of uniform sampling is

The upper bound just specifies that the sampling score be $\pi_{i}=1/N$ for all $i$. For experimentally designed sampling, we are allowed to study the graph structure and design sample indices. In the following corollary, we show a set of optimal sampling scores of the sampled projection estimator by minimizing the upper bound of MSE.

The optimal sampling score of the sampled projection estimator with bandwidth $\kappa\geq K$ is

To obtain the optimal sampling score for the estimator, we minimize the MSE given in Theorem 4 and solve the following optimization problem.

The objective function is the variance term of the MSE derived in Theorem 4. Since the bias term has nothing to do with the sampling score, minimizing the variance term is equivalent to minimizing the MSE. The constraints require $\{\pi_{i}\}_{i=1}^{N}$ to be a valid probability distribution. The Lagrangian function is then

where $\eta_{0},\eta_{i}$ are Lagrangian multipliers. We set the derivative of the Lagrangian function to zero,

Substituting the optimal sampling score $\pi_{i}$ to the upper bound of the MSE (13),

where (a) follows from substituting the optimal sampling score (14) into the upper bound of the MSE (13). ∎

We see that the optimal sampling score includes a trade-off between signal and noise. In practice, we cannot access $\mathbf{x}$ and thus need to approximate the ratio between each $x_{i}$ and $\sigma^{2}$. For active sampling, we can collect feedback to approximate each signal coefficient $x_{i}$; for experimentally designed sampling, we approximate beforehand; one way is to use the graph structure to sketch the shape of $\mathbf{x}$. We first use the Cauchy-Schwarz inequality to bound $x_{i}$,

where $\mathbf{v}_{i}$ is the $i$th row of $\operatorname{V}$. Thus, in the upper bound, we have a tradeoff between signal and noise: when the signal-to-noise ratio (SNR) $\left\|\mathbf{x}\right\|_{2}/\sigma^{2}$ is small, the approximate optimal sampling score is

which is the square root of the leverage score of $\operatorname{V}_{(K)}$.

When the SNR $\left\|\mathbf{x}\right\|_{2}/\sigma^{2}$ is large, the approximate optimal sampling score is

For approximately bandlimited signals, when $\beta$ is large and $\mu$ is small, the main energy concentrates in the first $K$ frequency components, $\left\|\mathbf{v}_{i}\right\|_{2}$ is concentrated in $\left\|\mathbf{v}_{i,(K)}\right\|_{2}$, where $\mathbf{v}_{i,(K)}$ is the first $K$ elements in $\mathbf{v}_{i}$. Specifically,

In this case, when the SNR $\left\|\mathbf{x}\right\|_{2}/\sigma^{2}$ is large, the approximate optimal sampling score is

which is the leverage score of $\operatorname{V}_{(K)}$.

V-C Relation to the Sampling Theory on Graphs

Sampling theory on graphs considers a bandlimited graph signal under the experimentally designed sampling . It shows that for a noiseless bandlimited graph signal, experimentally designed sampling guarantees perfect recovery while uniform sampling cannot, which also implies that active sampling cannot perform better than experimentally designed sampling. The recovery strategy is to solve the following optimization problem,

where $\Psi$ is the sampling operator (8) and $\mathbf{y}$ is a vector representation of the samples (7). When the original graph signal is bandlimited, the estimator (15) is unbiased and its MSE comes solely from the variance term caused by noise, that is,

In , the authors propose an optimal sampling operator based on minimizing $\left\|(\Psi\operatorname{V}_{(k)})^{\dagger}\right\|_{2}^{2}$. The sample set is deterministic and the procedure is computationally efficient when the sample size is small. Compared with sampling score-based sampling, however, the optimal sampling operator is computationally inefficient when the sample size is large. When the original graph signal is not bandlimited, similarly to $\mathbf{x}^{*}_{\rm LS}$, the solution of (15) is biased, that is,

We see that the high-frequency components $\operatorname{V}_{(-K)}\widehat{\mathbf{x}}_{(-K)}$ are projected on the low-frequency space spanned by $\operatorname{V}_{(K)}$, which causes aliasing.

VI Optimal Rates of Convergence

In this section, we introduce two types of graphs and show that the proposed recovery strategies achieve the optimal rates of convergence on both. Type-1 graphs model regular graphs, where the corresponding graph Fourier bases are not sparse and elements have similar magnitudes; examples are circulant graphs and nearest-neighbor graphs . Since the energy is evenly spread over the graph Fourier basis, each element contains similar amount of information; for such graphs, experimentally designed sampling performs similarly to uniform sampling. Type-2 graphs model irregular graphs, where the corresponding graph Fourier bases are sparse and elements do not have similar magnitudes; examples are small-world graphs and scale-free graphs . Since the energy is concentrated in a few elements in the graph Fourier basis, these elements are more informative and should be selected. For such graphs, experimentally designed sampling outperforms uniform sampling.

For a Type-1 graph, elements in $\operatorname{V}$ have roughly similar magnitudes, that is, the energy is evenly spread over $\operatorname{V}$.

Based on Theorem 4, we can specify the parameters for a Type-1 graph and show the following results.

Let $\mathbf{x}^{*}$ be the sampled projection estimator with the bandwidth $\kappa\geq K$ and uniform sampling; we have

where $C$ is a positive constant, $m$ is the number of samples (8), $\beta$ is the spectral decay factor in the approximately bandlimited class (6), and the rate is achieved when $\kappa$ is of the order of $m^{1/(2\beta+1)}$ and upper bounded by $N$;

Let $\mathbf{x}^{*}$ be the sampled projection estimator with the bandwidth $\kappa\geq K$ and the optimal sampling score in Corollary 2; we have

where $C$ is a positive constant, and the rate is achieved when $\kappa$ is of the order of $m^{1/(2\beta+1)}$ and upper bounded by $N$.

When $m\gg N$, we set $\kappa=N$, the bias term is then zero, and both upper bounds are actually $C\,m^{-1}$. We see that uniform sampling and optimal sampling score based sampling have the same convergence rate, that is, experimentally designed sampling does not perform better than uniform sampling for the Type-1 graphs.

Based on Theorem 3 and Corollary 3, we conclude the following.

where constant $C>c>0$, and the rate is achieved when $\kappa$ is of the order of $m^{1/(2\beta+1)}$ and upper bounded by $N$;

where constant $C>c>0$, and the rate is achieved when $\kappa$ is of the order of $m^{1/(2\beta+1)}$ and upper bounded by $N$.

We merge the proofs of Corollaries 3 and 4 in Appendix C in the supporting document.

We see that under both random and experimentally designed sampling, the lower and upper bounds have the same rate of convergence, which achieves the optimum. In addition, random and experimentally designed sampling have the same optimal rate of convergence and we can thus conclude that experimentally designed sampling does not perform asymptotically better than uniform sampling for the Type-1 graphs. The sampled projection estimator attains the optimal rate of convergence.

VI-B Type-2 Graphs

where $\mathbf{v}_{i,T}$ is the $T$th largest elements in the $i$th column of $\operatorname{V}$ and $T\ll N$ is some constant.

A Type-2 graph requires that each column vector of $\operatorname{V}$ be approximately sparse. When we take a few columns from $\operatorname{V}$ to form a submatrix, the energy in the submatrix concentrates in a few rows of the submatrix. This is equivalent to the sampling scores being approximately sparse. Simulations show that star graphs, scale-free graphs and small-world graphs approximately fall into this type of graphs.

Based on Theorem 4, we can specify the parameters for a Type-2 graph and show the following results.

Let $\mathbf{x}^{*}$ be the sampled projection estimator with the bandwidth $\kappa\geq K$ and uniform sampling; we have

where $C$ is a positive constant, and the rate is achieved when $\kappa$ is of the order of $m^{1/(2\beta+\gamma)}$, $\gamma=\log(N)/\log(\kappa)\geq 1$;

Let $\mathbf{x}^{*}$ be the sampled projection estimator with the bandwidth $\kappa\geq K$ and optimal sampling score based sampling; we have

where $C$ is a positive constant, the rate is achieved when $\kappa$ is of the order of $m^{1/(2\beta+1)}$ and upper bounded by $N$.

Based on Theorem 3 and Corollary 5, we conclude the following.

where constant $C>c>0$, and the rate is achieved when $\kappa$ is of the order of $m^{1/(2\beta+\gamma)}$ and $\gamma=\log(N)/\log(\kappa)\geq 1$;

Under experimentally designed sampling, there exists a $\gamma>1$,

where $C$ is a positive constant, the rate is achieved when $\kappa$ is of the order of $m^{1/(2\beta+1)}$ and upper bounded by $N$.

We merge the proofs of Corollaries 5 and 6 in Appendix D in the supporting document.

We see that under both uniform and experimentally designed sampling, the lower and upper bounds have the same rate of convergence, which achieves the optimum. However, experimentally designed sampling has a larger optimal rate of convergence, and we can thus conclude that experimentally designed sampling exhibits asymptotically better performance than uniform sampling for a Type-2 graph. The sampled projection estimator attains the optimal rate of convergence.

VII Experimental Results

In this section, we validate the proposed recovery strategy on five specific graphs: a ring graph, an Erdős-Rényi graph, a random geometric graph, a small-world graph and a power-law graph. Based on the graph structure, we roughly label each as a Type-1 or Type-2, and then, for each, we compare the empirical performance of the proposed recovery strategy based on uniform and experimentally designed sampling. For experimentally designed sampling, we use both the leverage score of $\operatorname{V}_{(K)}$ (approximately optimal in the noiseless case) and the square root of the leverage score of $\operatorname{V}_{(K)}$ (approximately optimal in the noisy case). In the experiments, the graph Fourier basis is the eigenvector matrix of the adjacency matrix. We find similar results when the graph Fourier basis is the eigenvector matrix of the graph Laplacian matrix.

We now introduce the five graphs, each with $10,000$ nodes, used in our experiments.

Ring graph with $k$-nearest neighbors. A ring graph is a graph where each node connects to its $k$-nearest neighbors. We use a ring graph where each node connects to exactly four nearest neighbors. The eigenvectors are similar to the discrete cosine transform and the energy evenly spreads to each element in $\operatorname{V}$ ; this is thus a Type-1 graph. Figure 3 illustrates some properties of the ring graph: Figure 3(a) shows the graph plot (for easier visualization, only $20$ nodes are shown and with enough zoom, one can clearly see that each node connects to exactly four nearest neighbors; Figure 3(b) shows the histogram of the degrees that concentrate on 4, as expected; and Figure 3(c) shows the histogram of the leverage scores of $\operatorname{V}_{(20)}$, which are the optimal sampling scores when the SNR is large. Note that we set the bandwidth $K=20$ to show the low-frequency band of the graph Fourier transform matrix. We see that the leverage scores concentrate around $10^{-4}$; this means that each node has the same probability to be chosen and uniform sampling is approximately the optimal sampling strategy.

Erdős-Rényi graph. An Erdős-Rényi graph is a random graph where each pair of nodes is connected with some probability . We use a graph where each pair of nodes is connected with probability of $0.01$, that is, each node has 100 neighbors in expectation. Figure 4 illustrates some properties of the Erdős-Rényi graph. Figure 4(a) shows the graph plot (for easier visualization, only 100 nodes are shown); Figure 4(b) shows the histogram of the degrees that concentrate around 100, as expected; and Figure 4(c) shows the histogram of the leverage scores of $\operatorname{V}_{(20)}$, which are the optimal sampling scores when the SNR is large. We see that the leverage scores concentrate around $10^{-4}$; this means that each node has the same probability to be chosen and uniform sampling is approximately the optimal sampling strategy, similarly to the ring graph. Based on the above, an Erdős-Rényi graph is approximately a Type-1 graph.

Random geometric graph. A random geometric graph is a spatial graph where each of the nodes is assigned random coordinates in the box $^{d}$ and an edge appears when the distance between two nodes is in a given range . We used a graph lying in the box $^{2}$, and two nodes are connected when the Euclidean distance is less than $0.03$. Figure 5 illustrates some properties of the random geometric graph: Figure 5(a) shows the graph plot; Figure 5(b) shows the histogram of the degrees that concentrate around 30, as expected (this matches previous assertion that given proper parameters, the degree distribution of a random geometric graph is the same as that of an Erdős-Rényi graph ); Figure 5(c) shows the histogram of the leverage scores of $\operatorname{V}_{(20)}$, which are the optimal sampling scores when the SNR is large; and Figure 5(d) shows the histogram of the leverage score on a log-scale. We see that the histogram of the leverage scores is skewed; this means that a few nodes are more important than other nodes during sampling and have much higher probabilities to be chosen. In , the authors show that the cluster properties are different for a random geometric graph and an Erdős-Rényi graph. The proposed sampling scores capture these cluster properties through the decomposition of the graph adjacency matrix. Based on the above, a random geometric graph is approximately a Type-2 graph.

Small-world graph. A small-world graph is a graph where most nodes are not neighbors of one another, but can be reached from any other node by a small number of hops (steps) . It is well known that many real-world graphs, including social networks, the connectivity of the Internet, and gene networks show small-world graph characteristics. We use a graph generated from the Watts–Strogatz model, where a ring graph is first built, followed by the rewiring of the edges with probability $0.01\%$. Figure 6 illustrates some properties of the small-world graph: Figure 6(a) shows the graph plot; Figure 6(b) shows the histogram of the degrees that concentrate around 2 (a few nodes have 6 neighbors because of rewiring); Figure 6(c) shows the histogram of the leverage scores of $\operatorname{V}_{(20)}$, which are the optimal sampling scores when the SNR is large; and Figure 6(d) shows the histogram of the leverage scores on a log-scale. We see that the histogram of the leverage scores is skewed; this means that a few nodes are more important than others during sampling and have much higher probabilities to be chosen, similarly to the random geometric graph. Based on the above, a small-world graph is approximately a Type-2 graph.

Power-law graph. A power-law graph is a graph where the more connected a node is, the more likely it is to receive new links, known as a preferential attachment graph . It is well known that the degree distribution of a preferential attachment graph follows the power law. We use a graph generated from the Barabási-Albert model, where new nodes are added to the network one at a time. Each new node is connected to one existing node with a probability that is proportional to the number of links that the existing nodes already have. Figure 7 illustrates some properties of the small-world graph: Figure 7(a) shows the graph plot; Figure 7(b) shows the histogram of the degrees that is skewed, which clearly follows the power law; Figure 7(c) shows the histogram of the leverage scores of $\operatorname{V}_{(20)}$, which are the optimal sampling scores when the SNR is large; and Figure 7(d) shows the histogram of the leverage scores on a log-scale. We see that the histogram of the leverage scores is skewed; this means that a few nodes are more important than others during sampling and have much higher probabilities to be chosen, similarly to the random geometric graph and the small-world graph. Based on the above, a power-law graph is approximately a Type-2 graph.

Types. Based on our observation of the graph Fourier transform matrices of each graph, the ring graph and the Erdős-Rényi graph approximately satisfy the requirement to be the Type-1 graph, while the random geometric graph, the small-world graph and the power-law graph approximately satisfy the requirement to be the Type-2 graph. We thus expect that the experimentally designed sampling has similar performance to uniform sampling for the ring graph and the Erdős-Rényi graph, while it outperforms uniform sampling for the random geometric graph, the small-world graph and the power-law graph.

VII-B Simulated Graph Signals

For each graph $\operatorname{A}$, we generate 1,000 graph signals through the following two steps: We first generate the graph frequency components as

We then normalize $\widehat{\mathbf{x}}$ to have unit norm, and obtain $\mathbf{x}=\operatorname{V}\widehat{\mathbf{x}}$. It is clear that $\mathbf{x}\in\operatorname{ABL}_{\operatorname{A}}(K,\beta,\mu)$, where $K=10$ and $\beta$ varies as $0.5$ and $1$. During sampling, we simulate noise $\epsilon\sim\mathcal{N}(0,\sigma^{2})$, vary the sample size $m$ from $1,000$ to $20,000$, and vary $\sigma^{2}$ from low noise level $10^{-4}$ ($\left\|\mathbf{x}\right\|_{2}/\left\|\epsilon\right\|_{2}=100$) to high noise level $0.02$ ($\left\|\mathbf{x}\right\|_{2}/\left\|\epsilon\right\|_{2}=0.5$). During recovery, we set the bandwidth $K=\max(10,m^{1/2\beta+1})$ as suggested in Corollaries 4 and 6.

VII-C Results

We compare four sampling strategies, including uniform sampling (in blue), leverage score based sampling (in orange), square root of the leverage score based sampling (in purple), and degree based sampling (in red). Note that the last three sampling strategies all belong to experimentally designed sampling because they are designed based on the structure of the graph. As shown in Section V-B, leverage score based sampling is approximately optimal when the SNR is large, square root of the leverage score based sampling is approximately the optimal when the SNR is small. We also use the degrees as the sampling scores because previous works show that the largest eigenvectors of adjacency matrices often have most of their mass localized on high-degree nodes , which implies the high correlation between degree and leverage score. We evaluate the recovery performance by using the MSE,

where $\mathbf{x}^{*}$ is the recovered graph signal and $\mathbf{x}$ is the original graph signal. The simulation results for the ring graph, the Erdős-Rényi graph, the random geometric graph, the small-world graph and the power-law graph are shown in Figures 8, 9, 10, 11 and 12, respectively. We summarize the important points below.

All of the sampling strategies perform similarly on the ring graph and the Erdős-Rényi graph, matching Corollary 4.

Experimentally designed sampling outperforms uniform sampling on the random geometric graph, the small-world graph and the power-law graph, matching Corollary 6. Especially for the small-world graph and the power-law graph, uniform sampling is much worse than experimentally designed sampling.

Leverage score based sampling outperforms all other sampling strategies when the noise level is low.

Square root of the leverage score based sampling outperforms all other sampling strategies when the noise level is high.

Degree based sampling outperforms uniform sampling because of its correlation to the leverage score based sampling, but for the small-world graph, degree based sampling is still much worse than leverage score based sampling.

When $\beta$ is larger, the recovery performance is better because less energy is concentrated in the high-frequency band for approximately bandlimited graph signals.

When $\sigma^{2}$ is smaller, the recovery performance is better because of less noise.

The degree distribution is not a reliable indicator of when experimentally designed sampling outperforms uniform sampling. The degree distributions of the Erdős-Rényi graph and the random geometric graph are similar, but experimentally designed sampling only outperforms uniform sampling on the random geometric graph. This implies that the first-order information provided by the degree is not sufficient in designing samples.

VII-D Discussion

Graph Fourier basis is critical for understanding the graph structure. For example, given the graph Fourier basis, Theorem 3 shows that active sampling does not fundamentally outperform experimentally designed sampling while Corollary 6 shows that experimentally designed sampling fundamentally outperforms uniform sampling on type-2 graphs. In other words, graph Fourier basis is critical for choosing samples.

For irregular (Type-2) graphs, using experimentally designed sampling to choose anchor points can fundamentally aid semi-supervised learning (not true for regular, Type-1 graphs), a technique for training classifiers with both labeled and unlabeled data. Semi-supervised learning assumes that unlabeled data can provide distribution information to build a stronger classifier . Many algorithms for semi-supervised learning are based on graphs that are constructed from a given dataset , often by modeling each node as a data sample and connecting two nodes by an edge if the distance between their features is in a given range, which is similar to the construction of random geometric graphs. Based on the assumption that adjacent nodes have similar labels, semi-supervised learning diffuses label probabilities from labeled data to unlabeled data along the graph structure and classifies unlabeled data according to those label probabilities. While in some works , training data is selected uniformly and randomly, in others, algorithms are designed adapted to the structure , which is essentially equivalent to experimentally designed sampling we propose. In other words, experimentally designed sampling is used implicitly without being able to articulate when and why it works; this paper, on the other hand, provides a comprehensive explanation of why experimentally designed sampling helps semi-supervised learning through showing the lower and upper bounds of three sampling strategies.

VIII Conclusions

We build a theoretical foundation for the recovery of a newly proposed class of smooth graph signals, approximately bandlimited signals, which generalizes the class of bandlimited graph signals, under uniform sampling, experimentally designed sampling and active sampling. We show that experimentally designed sampling and active sampling have the same fundamental limitations, and can outperform uniform sampling on irregular graphs. We propose a recovery strategy and analyze its statistical properties. We show that the proposed recovery strategy attains the optimal rates of convergence on two specific types of graphs. To validate the recovery strategy, we test it on five specific types of graphs: a ring graph with $k$ nearest neighbors, an Erdős-Rényi graph, a random geometric graph, a small-world graph and a power-law graph, and show that experimental results match the proposed theory well. This work also gives a comprehensive explanation for why experimentally designed sampling works for semi-supervised learning with graphs and shows the critical role of the graph Fourier basis in analyzing graph structures.

IX Acknowledgment

We gratefully acknowledge support from the NSF through awards 1130616,1421919, the University Transportation Center grant (DTRT12-GUTC11) from the US Department of Transportation, and the CMU Carnegie Institute of Technology Infrastructure Award. We also thank the editor and the reviewers for comments that led to improvements in the manuscript. Initial parts of this work were presented at SampTA 2015 .

References

Appendix A Appendices

We aim to construct a typical set of vectors in $\mathcal{F}$, and use the Fano’s method. Let $\mathbf{v}_{k}$ be the $k$th column of $\operatorname{V}$, $\mathbf{v}^{(i)}$ be the $i$th row of $\operatorname{V}$, $\mathcal{X}$ be a pruned hypercube and

where $\kappa_{0}$ is no smaller than the bandwidth $K$,

and $0<c<1$. It is easy to check that $\mathcal{F}^{\prime}\subseteq\mathcal{F}$. Let $d(\mathbf{x},\mathbf{y})=\left\|\mathbf{x}-\mathbf{y}\right\|_{2}$; we thus have

where $(a)$ follows from $k\leq 2K$, and $(b)$ from the Varshamov-Gilbert lemma. To use the Fanno’s method, we need to bound the Kullback-Leibler divergence,

where $\delta=\sum_{k,k^{\prime}=\kappa_{0},k\neq k^{\prime}}^{2\kappa_{0}-1}(-1)^{k+k^{\prime}}\frac{c\mu\left\|\mathbf{x}\right\|_{2}^{2}w_{k}w_{k^{\prime}}\operatorname{V}_{ik}\operatorname{V}_{ik^{\prime}}}{\kappa_{0}\sqrt{1+k^{2\beta}}\sqrt{1+k^{\prime 2\beta}}}$ is small because of the cross signs. For uniform sampling, the sampling set $\mathcal{M}$ is chosen randomly, thus, we have

For experimentally designed sampling, we can choose the sampling set $\mathcal{M}$ to maximize $\sum_{i\in\mathcal{M}}\sum_{k=\kappa_{0}}^{2\kappa_{0}-1}\operatorname{V}_{ik}^{2}$; we thus have

For active sampling, we cannot get more benefit from signal coefficients, so the KL divergence is the same of the experimentally designed sampling.r By Fanno’s lemma, we finally have the lower bounds for three sampling strategies. $\blacksquare$

A-B Proof of Theorem 4

We aim to bound the MSE by splitting to a bias term and a variance term.

where the first term is bias and the second term is variance. For each element in the bias term, we have

where $(a)$ follows from the independence of each sample. For all the elements, we have

We next bound the variance term by splitting into two parts, with and without noise. For each element in the variance term, we have

where $\Delta^{(1)}$ is the variance from noise and $\Delta^{(2)}$ is the variance from sampling. To bound $\Delta^{(1)}_{i}$, we have

We combine the bounds for both $\Delta^{(1)}_{i}$ and $\Delta^{(2)}_{i}$, and obtain the bounds for the variance term,

which leads to Lemma 2. Finally, we obtain the MSE by combine the bias term and the variance term,

A-C Proof of Corollaries 3 and 4

For a Type-1 graph, we assume that each element in an approximately bandlimited signal has a similar magnitude. Since all the elements in $\operatorname{V}$ have the same magnitude, each element of a graph signal, $x_{i}=\mathbf{v}_{i}^{T}\widehat{\mathbf{x}}$, should have a similar magnitude. In other words, $N\max_{i}x_{i}^{2}$ and $\left\|\mathbf{x}\right\|_{2}^{2}$ are of the same order.

For uniform sampling, based on Corollary 1, we have

where $(a)$ follows from $\operatorname{U}$ being orthornormal, $\kappa$ being of the order of $m^{\frac{1}{2\beta+1}}$, and $C>0$ some constant. Since at least the sampled projection estimator satisfies this rate of convergence, we have

We next show the lower bound. Based on Theorem 3, we have

where $\kappa_{0}$ is of the order of $m^{\frac{1}{2\beta+1}}$.

For optimal sampling scores based sampling, based on Corollary 2, we have

where $(a)$ follows from that, based on Definition 5, $\left\|\operatorname{U}_{(\kappa)}\right\|_{2,1}^{2}=\left(N\sqrt{\kappa(\frac{c}{\sqrt{N}})^{2}}\right)^{2}=O(N\kappa)$, which is of the same order of $N\left\|\operatorname{U}_{(\kappa)}\right\|_{F}^{2}$. Since at least the sampled projection estimator satisfies this rate of convergence, we have

Based on Definition 5, $\left\|\operatorname{V}_{(2,\kappa_{0})}\right\|_{\infty,2}^{2}=\kappa_{0}(\frac{c}{\sqrt{N}})^{2}=c^{2}\kappa_{0}/N$, which is of the same order of $\left\|\operatorname{V}_{(2,\kappa_{0})}\right\|_{F}^{2}/N$, we have

where $\kappa_{0}$ is of the order of $m^{\frac{1}{2\beta+1}}$. $\blacksquare$

A-D Proof of Corollaries 5 and 6

For a Type-2 graph, we assume that a few elements in an approximately bandlimited signal have much higher magnitudes than the others. The intuition is that, since $\operatorname{V}$ is sparse, the energy concentrates in $O(\kappa)$ rows of $\operatorname{V}_{(\kappa)}$, thus, $O(\kappa)$ components of an approximately bandlimited signal, $x_{i}\approx\mathbf{v}_{i,(\kappa)}^{T}\widehat{\mathbf{x}}_{(\kappa)}$, have much higher magnitudes than the others. In other words, $\kappa\max_{i}x_{i}^{2}$ and $\left\|\mathbf{x}\right\|_{2}^{2}$ are of the same order.

For uniform sampling, based on Corollary 1, we have

where $\gamma$ varies with $\kappa$ to satisfy $\kappa^{\gamma}\leq N$ ($\gamma>1$) and $\kappa$ is of the order of $m^{\frac{1}{2\beta+\gamma}}$, and $C>0$ is some constant. Since at least Algorithm 1 satisfies this rate of convergence, we thus have

We next show the lower bound. Based on Theorem 3, we have

where $\kappa_{0}$ is of the order of $m^{\frac{1}{2\beta+\gamma}}$.

For optimal sampling scores based sampling, based on Corollary 2, we have

where $(a)$ follows from the energy concentrating in $O(\kappa)$ columns of $\operatorname{U}_{(\kappa)}$ as shown in Definition 6 and the upper bound reaching the minimum when $\kappa$ is of the order of $m^{\frac{1}{2\beta+1}}$. Since at least Algorithm 1 satisfies this rate of convergence, we have

Based on Definition 6, $\left\|\operatorname{V}_{(2,\kappa_{0})}\right\|_{\infty,2}^{2}=c$, we have