A Probabilistic Interpretation of Sampling Theory of Graph Signals

We give a probabilistic interpretation of sampling theory of graph signals. To do this, we first define a generative model for the data using a pairwise Gaussian random field (GRF) which depends on the graph. We show that, under certain conditions, reconstructing a graph signal from a subset of its samples by least squares is equivalent to performing MAP inference on an approximation of this GRF which has a low rank covariance matrix. We then show that a sampling set of given size with the largest associated cut-off frequency, which is optimal from a sampling theoretic point of view, minimizes the worst case predictive covariance of the MAP estimate on the GRF. This interpretation also gives an intuitive explanation for the superior performance of the sampling theoretic approach to active semi-supervised classification.


INTRODUCTION
Graph signal processing aims to extend the tools for analysis, approximation, denoising and interpolation of traditional signals to signals defined on graphs. The advantage of this framework is that it allows us to process the given data while taking into consideration the underlying connectivity between the data points. The graph can be inherent to the data as is the case in application areas such as social networks and sensor networks or it can be constructed using the data to capture the underlying geometry. Examples of the latter are found in image processing and machine learning (see [1,2]).
In this paper, we focus on the sampling theory of graph signals. The classical Nyquist-Shannon sampling theorem says that a signal with bandwidth f is uniquely determined by its (uniformly spaced) samples if the sampling rate is higher than 2f . Intuitively, it tells us how "smooth" the signal has to be, for perfect recovery, given the sampling density, and vice versa. Moreover, the signal can be perfectly reconstructed from the samples by a simple low pass filter. Sampling theory of graph signals similarly deals with the problem of reconstructing an unknown graph signal from its samples on a subset of nodes. Frequency domain representation of graph signals is given by the eigenvectors and eigenvalues of the Laplacian matrix associated with the graph. In order to pose a sampling theorem analogous to the Nyquist-Shannon sampling theorem, we need to find the maximum bandwidth (in the graph spectral domain) that a graph signal can have so that it is uniquely determined by its samples on the given subset of nodes. Conversely, given the bandwidth, we need to find the smallest subset of nodes, so that recovery of any signal with that bandwidth, from its samples on that subset, is unique and stable.
This work was supported in part by NSF under grant CCF-1410009.
Given that the signal is smooth enough to be uniquely represented by its samples on a subset of nodes, we need to give an efficient and stable algorithm to reconstruct the unknown samples. These questions have been answered to some extent in [3,4,5,6]. We discuss some of these results in Section 2.2.
This sampling theoretic perspective has been shown to be very useful for graph based active semi-supervised learning [7]. In this context, label prediction is considered as a graph signal reconstruction problem. The characterization of a subset of nodes given by the sampling theory, namely the associated cutoff frequency is used as a criterion function to choose the optimal set nodes to be labelled for active learning.
Sampling theoretic approaches for active and semi-supervised learning [7] are purely deterministic. However, their probabilistic interpretation is desired for the following reasons: 1. It allows us to understand them as model based methods and thus, makes it easier to include them as components of a larger probabilistic model. 2. It can also suggest a principled way to refine the model parameters (which are given by the underlying graph) as more data is observed (see [8] for an example). 3. The interpretation presented in this paper assumes a Gaussian random field model for the data. This may lead to generalizations of the sampling theory to data with non-Gaussian distributions which might be more realistic for a classification problem. 4. This interpretation also makes the relationship between the sampling theoretic approach and previously proposed semi-supervised [9] and active learning [10,11] methods more apparent as discussed in Section 5.
The main contributions of this paper are the following. We define a generative model for graph signals using a pairwise Gaussian random field (GRF) with a covariance matrix that depends on the graph. We show that, when conditions of the graph signal sampling theorem are satisfied, bandlimited reconstruction of a graph signal from a subset of its samples is equivalent to performing MAP inference on a low rank approximation of the above GRF. This learning model performs very well in classification problems, as demonstrated in the experiments, since the true data covariance matrix is expected to be close to low rank. We then show that a sampling set of given size with the largest associated cut-off frequency, which is optimal from a sampling theoretic point of view, minimizes the worst case predictive covariance of the MAP estimate on the GRF.

Preliminaries and Notation
We consider a connected, undirected and weighted graph G = (V, E). The nodes V in the graph are indexed by {1, 2, . . . , N }. S c denotes the complement of S in V, i.e., S c = V \ S. The edge set E is given by {(i, j, wij)}, where i, j ∈ V and wij ∈ R + .
(i, j, wij) denotes an edge with weight wij connecting nodes i and j. The connectivity information given by E is encoded by the adjacency matrix W of size N × N with W(i, j) = wij. The degree matrix D is a diagonal matrix diag{d1, . . . dN }, where di = j wij is the degree of node i. The Laplacian matrix is defined as L = D − W. The symmetric normalized form of the Laplacian is given by L = D −1/2 LD −1/2 . A graph signal f : V → R is a mapping which takes a real value on each node of the graph. It can be represented as f = (f1, . . . , fN ) ∈ R N . For x ∈ R N , xS denotes a sub-vector of x consisting of its components indexed by S. Similarly, for A ∈ R N ×N , AS 1 S 2 is the sub-matrix of A with rows indexed by S1 and columns indexed by S2. For simplicity, we denote ASS by AS . We use λmax[.] and λmin[.] to denote the largest and the smallest eigenvalue of a matrix, respectively. tr(.) denotes the trace of a matrix. A + is used to denote the pseudo-inverse of A. 1 and 0 denote vectors or matrices of ones and zeros, respectively.
It can be shown that L and L are positive semi-definite. Hence, L has real eigenvalues 0 = λ1 < λ2 ≤ . . . ≤ λN and a corresponding orthogonal set of eigenvectors {u 1 , u 2 , . . . , u N }. It can be diagonalized as L = UΛU , where U = (u 1 , . . . , u N ) and Λ = diag{λ1, . . . , λN }. Variation in the eigenvectors of L over the graph (as captured by u Lu = i,j wij(ui − uj) 2 ) increases as the corresponding eigenvalues increase. Thus, these eigenvectors allow us a to define a graph dependent notion of frequency for the graph signals. The so-called Graph Fourier Transform (GFT) 1 is defined asfi = f , u i (or in an equivalent matrix formf = U f ), wherefi is the GFT coefficient corresponding to frequency λi. An ω-bandlimited signal has its GFT supported on [0, ω], i.e.,fi = 0 for λi > ω. Conversely, such a signal is said to have a bandwidth equal to ω. If {λ1, . . . , λr} are the eigenvalues less than ω, then any ω-bandlimited signal can be written as a linear combination of corresponding eigenvectors where a is the coefficient vector. The space of ω-bandlimited signals is called a Paley-Wiener space P Wω(G).

Sampling Theorem and Bandlimited Reconstruction
Sampling theory deals with the problem reconstructing an ωbandlimited signal f from its samples fS on the nodes in S ⊆ V.
There are three important questions that need to be answered in this context: 1. Given S, what is the maximum bandwidth ω that f can have so that it is uniquely determined by fS ? 2. Which is the best sampling set Sopt of a given size m? 3. Given that f is uniquely determined by fS , how to find the unknown samples fSc ? We briefly review some of the results related to each of the above problems. Let L2(S c ) be the space of signals which are identically zero on S but can have non-zero samples on S c , i.e., gS = 0 ∀g ∈ L2(S c ). It is easy to see that for all signals in P Wω(G) to be uniquely determined by their samples on S, we need P Wω(G) ∩ L2(S c ) = {0}. This observation leads to the following theorem.
Theorem 1 (Sampling Theorem [6]). Any signal in P Wω(G) can be uniquely reconstructed from its samples on a subset of nodes S if and only if ω < inf where ω(.) denotes the bandwidth of a signal. If the above condition is satisfied, then S is said to be a uniqueness set for P Wω(G).
To ensure unique recovery of a signal from its samples on S, its bandwidth has to be less than inf g∈L 2 (S c ) ω(g). This is called the cut-off frequency associated with the subset S and is denoted by ω(S). An estimate of the cut-off frequency is given by [6] Ω k (S) = λmin (L k )Sc It can be shown that Ω k (S) ≤ ω(S) and we get closer to ω(S) as k increases.
A larger cut-off frequency estimate Ω k (S) implies that a bigger space of signals can be perfectly recovered from their samples on S. Therefore, Ω k (S) can be used as a criterion function to be maximized for choosing the optimal sampling set Sopt of given size m, i.e., The above problem is combinatorial and NP-hard. A greedy algorithm for finding an approximate solution is proposed in [6].
Consider a signal f ∈ P Wω(G) with ω < ω(S). Using the representation of a bandlimited signal in (1), we get that fS = USRa. Since f is uniquely sampled on S, USR must have full column rank so that the least squares solution a of the above system of equations is unique. The unknown samples can then be reconstructed by: A faster, iterative method for bandlimited reconstruction is proposed in [5], which does not need the computation of eigenvectors. These sampling theory based algorithms for subset selection and signal reconstruction have been applied to graph based active semisupervised learning and are shown to perform better than many state of the art approaches [7].

GRF MODEL FOR GRAPH SIGNALS
In order to give a probabilistic interpretation of the graph signal processing framework, we define a generative model for the signal using a pairwise Gaussian Random Field (GRF) based on the graph G. A random signal f = (f1, . . . , fN ) is assumed to be drawn from the following distribution: where I denotes an identity matrix of size N × N . Let K be the covariance matrix of the the GRF. Then, from the above equation, the inverse covariance matrix (also known as the precision matrix) can be written as: Note that K has the same eigenvectors as L, while the corresponding eigenvalues are σi = 1 λ i +δ . Thus, K can be diagonalized as where Σ = diag{σ1, . . . , σN }. The advantage of introducing the parameter δ is that it leads to a non-singular precision matrix and thus, allows us to have a proper covariance matrix. σ1 = 1/δ can be thought of as the variance of the DC component of f since u 1 = 1.

SAMPLING THEORY AND INFERENCE OVER GRF
Consider a signal f generated using the GRF defined in (6) with covariance matrix K = (L + δI) −1 . As in the sampling problem, we observe the samples of f on a subset S of nodes. Our goal is to estimate the unknown samples. It is well known that the conditional distribution of fSc given fS equals N (µ S c |S , K S c |S ), where µ S c |S = KScS (KS ) + fS and (9) are the MAP estimate and the predictive covariance matrix of fSc given fS , respectively [9,12].

Bandlimited Reconstruction as MAP Inference
Let λr be the largest eigenvalue of L which is less than ω. We definê K to be a low rank approximation of K which only contains the spectral components corresponding to {λ1, . . . , λr}, i.e., Consider the problem of reconstructing a random signal generated using a GRF with covarianceK, from its samples on S. The following theorem shows that, if conditions of the sampling theorem are satisfied, then the error of bandlimited reconstruction is zero.
Theorem 2. Let f be a random graph signal generated using the GRF with covarianceK given by (11). LetfSc be the bandlimited reconstruction of fSc obtained from its samples on S, where S is a uniqueness set for P Wω(G). Then, fSc −fSc = 0.
Before proving the above theorem, we show, in the following lemma, that bandlimited reconstruction is equivalent to MAP inference on the GRF with covarianceK. Lemma 1. Let S ⊆ V be a uniqueness set for P Wω(G). Then the MAP estimate of fSc given fS in a GRF with covariance matrixK is equal to the bandlimited reconstruction given by (5).
Proof. Under a permutation which groups together nodes in S c and S, we can writeK as the following block matrix Therefore, we can write the MAP estimate obtained with covariancê K as,μ Because ω < ω(S), we have that USR has full column rank and equivalently, U SR has full row rank. Therefore, we can write (USRΣRU SR ) + = (U SR ) + Σ + R U + SR and U + SR = (U SR USR) −1 U SR . Simplifying (13) using these equalities leads tof which is equal to the least squares solution given in (5).
Proof of Theorem 2. From Lemma 1,fSc =μ S c |S . Therefore, Using the block form ofK in (12), and the fact that USR has full column rank, it is easy to show thatK S c |S = 0, which implies E( fSc −fSc 2 ) = 0. But since, fSc −fSc ≥ 0, we get fSc −fSc = 0.

Cut-off Frequency and Estimation Error
If the true covariance matrix is only approximately low rank, then MAP inference withK gives a non-zero reconstruction error. The best sampling set in this case is the one which minimizes the predictive covariance. According to the sampling theory of graph signals, the optimal sampling set of given size is the one which has the largest associated cut-off frequency. We show that finding a sampling set S which maximizes a crude estimate of the cut-off frequency Ω1(S) is equivalent to minimizing the maximum eigenvalue of the predictive covariance of fSc given fS . Proof. Consider a block matrix representation of K similar to (12).
Using the block matrix inversion formula, we can write K −1 as are the Schur complements of KS and KSc respectively. LSc = (K −1 )Sc − δISc = S −1 K S − δISc . Note that S K S = K S c |S . Thus, the estimated cut-off frequency corresponding to the subset S of nodes can be written in terms of the conditional covariance matrix The result readily follows from this.
A sampling set with the largest estimated cut-off frequency Ω1(S) also minimizes the worst case prediction error of the MAP estimate on a GRF with K = (L + δI) −1 . However, as shown in Lemma 1, bandlimited signal reconstruction is equivalent to MAP estimation with a low rank approximation of K. Intuitively, a better estimate of the predictive covariance, in this model of signal reconstruction, can be obtained with by ((K k ) S c |S ) 1/k with larger values of k as it gives more weight to the principal components with larger variance. This justifies the use of Ω k (S) with k > 1 as a criterion for active learning.

Justification for the Sampling Theoretic Approach to Active Semi-supervised Classification
MAP estimation is optimal for reconstructing signals generated using a GRF with a full rank covariance matrix, because it minimizes the mean squared error of estimation. Moreover, since the estimation error equals tr(K S c |S ), an optimal sampling set of size m is given by arg min |S|=m tr(K S c |S ). Indeed, this is the so-called V -optimality criterion for active learning proposed in [10].
However, in a classification problem, data points in the same class are highly correlated whereas data points in different classes have very small correlation. Since the number of classes is typically very small compared to the number of data points, we expect the (unknown) "true" covariance matrix to be very well-approximated by a low rank matrix [13]. Thus, bandlimited interpolation is a better model for signal reconstruction in this context, since it is equivalent to MAP estimation with a low rank covariance matrix. Maximizing the cut-off frequency is a natural set selection criterion for this learning model.

RELATED WORK
Different criteria have been proposed for batch mode active learning on Gaussian random fields. The approach presented in [14] selects the points to label such the mutual information between the labelled and unlabelled data points is maximized. Our sampling theoretic approach (4) is more similar to the methods proposed in [10,11]. These methods use MAP estimation on GRF [9] as their model for label prediction. As stated before, [10] chooses the sampling set S by minimizing tr(K S c |S ). The method in [11], on the other hand, tries to minimize ij (K S c |S )ij (also known as Σ-optimality criterion). This is equivalent to minimizing the risk of the surveying problem [15] (which is the problem of determining the proportion of nodes belonging to one class). All the above methods are closely related to the optimal design of experiments [16]. Experiment design deals with the problem of estimating a vector from a set of linear measurements. The goal is to choose the optimal set of m measurements so that the estimation error is minimized. Different error measures lead to different optimality criteria. For example, minimizing the trace of estimation covariance leads to A-optimal design whereas minimizing its determinant gives the D-optimal design. The sampling theoretic approach is closer to the so-called E-optimal design which minimizes the worst case prediction error given by the maximum eigenvalue of the predictive covariance matrix.

EXPERIMENTS
To demonstrate the effectiveness of the framework of sampling theory, we first apply it to the problem of graph based active semisupervised classification. In our experiment, we use a subset of the USPS handwritten digit dataset containing 100 16 × 16 images each of digits 0 to 9. We construct a weighted K-NN graph of 1000 nodes with K = 10 and the similarities given by wij = exp − x i −x j 2 σ 2 . The problem is to choose the nodes to be labelled and then predict the unknown labels from the queried labels. We consider different combinations of active learning criteria and learning models. As expected from the discussion in Section 4.3, selecting the sampling set by maximizing the cutoff frequency and then performing bandlimited reconstruction outperforms Σ and Voptimality criteria used in conjunction with MAP estimation (see Figure 1(a)). Even if the learning model is fixed to bandlimited interpolation, the sampling theoretic approach gives better results as seen in Figure 1(b)). This is because maximizing the cutoff frequency is a more suitable set selection criterion under this model.
On the other hand, if we consider the problem of regression of a random real valued graph signal generated using a covariance matrix that is not low rank, a V -optimal set is expected to give a better SNR of reconstruction. This is demonstrated in Figure 2 where we reconstruct a random real valued signal generated with the covariance matrix obtained using the graph from the previous example.

CONCLUSION AND FUTURE WORK
In this paper, we gave a probabilistic interpretation for the sampling theory of graph signals. We showed that if the data is generated using a Gaussian random field whose precision matrix equals the graph Laplacian, then bandlimited reconstruction is equivalent to the MAP inference on an approximation of this GRF which has a low rank covariance matrix. Moreover, an optimal sampling set obtained via sampling theory minimizes the worst case predictive covariance of MAP estimation on the GRF.
A probabilistic interpretation allows us to view graph signal sampling theory as a model based method. It would be interesting to consider it as part of a larger probabilistic model which refines the covariance matrix as more data is observed. This interpretation also suggests a generalization of the sampling theory to non-Gaussian models which might be more realistic for some applications.