Random sampling of bandlimited signals on graphs

Introduction

Graphs are a central modelling tool for network-structured data . Depending on the application, the nodes of a graph may represent people in social networks, brain regions in neuronal networks, or stations in transportation networks. Data on a graph, such as individual hobbies, activity of brain regions, traffic at a station, may be represented by scalars defined on each node, which form a graph signal. Extending classical signal processing methods to graph signals is the purpose of the emerging field of graph signal processing .

Natural choices of smoothness models build upon, e.g., the graph’s adjacency matrix, the combinatorial Laplacian matrix, the normalised Laplacian, or the random walk Laplacian. The sets of eigenvectors of these operators define different graph Fourier bases. Given such a Fourier basis, the equivalent of a classical $\omega$-bandlimited signal is a $k$-bandlimited graph signal whose $k$ first Fourier coefficients are non-null .

Unlike continuous time signal processing, the concept of regular sampling itself is not applicable for graph signals, apart for very regular graphs such as bipartite graphs . We are left with two possible choices for sampling: irregular or random sampling. Irregular sampling of $k$-bandlimited graph signals has been studied first by Pesenson who introduced the notion of uniqueness set associated to the subspace of $k$-bandlimited graph signals. If two $k$-bandlimited graph signals are equal on a uniqueness set, they are necessarily equal on the whole graph. Building upon this first work, and using the fact that the sampling matrix applied to the first $k$ Fourier modes should have rank $k$ in order to guarantee recovery of $k$-bandlimited signals, Anis et al. and Chen et al. showed that a sampling set of size $k$ that perfectly embeds $k$-bandlimited signals always exists. To find such an optimal set, the authors need to compute the first $k$ eigenvectors of the Laplacian, which is computationally prohibitive for large graphs. A recent work bypasses the partial diagonalisation of the Laplacian by using graph spectral proxies, but the procedure to find an optimal sampling set still requires a search over all possible subsets of nodes of a given size. This is a very large combinatorial problem. In practice, approximate results are obtained using a greedy heuristic that enables the authors to efficiently perform experiments on graphs of size up to few thousands nodes.

Several other sampling schemes exist in the literature, such as schemes based on a bipartite decomposition of the graph , on a decomposition via maximum spanning trees , on the sign of the last Fourier mode , on the sign of the Fiedler vector , or on a decomposition in communities . All these propositions are however specifically designed for graph multiresolution analysis with filterbanks, and are not suited to find optimal or close-to-optimal sets of nodes for sampling $k$-bandlimited graph signals.

In this paper, we propose a very different approach to sampling on graphs. Instead of trying to find an optimal sampling set, (i.e., a set of size $k$) for $k$-bandlimited signals, we relax this optimality constraint in order to tackle graphs of very large size. We allow ourselves to sample slightly more than $k$ nodes and, inspired by compressive sampling, we propose two random sampling schemes that ensure recovery of graph signals with high probability.

A central graph characteristic that appears from our study is the graph weighted coherence of order $k$ (see Definition 2.1). This quantity is a measure of the localisation of the first $k$ Fourier modes on the nodes of the graph. Unlike the classical Fourier modes, some graph Fourier modes have the surprising potential of being localised on very few nodes. The farther a graph is from a regular grid, the higher the chance to have a few localised Fourier modes. This particularity in graph signal processing is studied in but is still largely not understood.

First, we propose a non-adaptive sampling technique that consists in choosing a few nodes at random to form the sampling set. In this setting, we show that the number of samples ensuring the reconstruction of all $k$-bandlimited signals scales with the square of the graph weighted coherence. For regular or almost-regular graphs, i.e., graphs whose coherence is close to $\sqrt{k}$, this result shows that $O(k\log{k})$ samples selected using the uniform distribution are sufficient to sample $k$-bandlimited signals. We thus obtain an almost optimal sampling condition.

Second, for arbitrary graphs with a coherence potentially tending to $\sqrt{n}$, where $n\gg k$ is the total number of nodes, we propose a second sampling strategy that compensates the undesirable consequences of mode localisation. The technique relies on the variable density sampling strategy widely used in compressed sensing . We prove that there always exists a sampling distribution such that no more than $O(k\log{k})$ samples are sufficient to ensure exact and stable reconstruction of all $k$-bandlimited signals, whatever the graph structure. Unfortunately, computing the optimal sampling distribution requires the partial diagonalisation of the first $k$ eigenvectors of the Laplacian. To circumvent this issue, we propose a fast technique to estimate this optimal sampling distribution accurately.

Finally, we propose an efficient method to reconstruct any $k$-bandlimited signal from its samples. We prove that the method recovers $k$-bandlimited signals exactly in the absence of noise. We also prove that the method is robust to measurement noise and model errors.

Note that our sampling theorems are applicable to any symmetrical Laplacian or adjacency matrix, i.e., any weighted undirected graph. Nevertheless, the efficient recovery method we propose is specifically designed to take advantage of the semi-definite positivity of the Laplacian operator. In the following, we therefore concentrate on such symmetrical positive semi-definite Laplacians, such as the combinatorial or normalized Laplacians.

Let us acknowledge that the idea of random sampling for $k$-bandlimited graph signals is mentioned in and . In , the authors prove that the space of $k$-bandlimited graph signals can be stably embedded using a uniform sampling but for the Erdős-Rényi graph only. The idea of using a non-uniform sampling appears in . However, the authors do not prove that this sampling strategy provides a stable embedding of the space of $k$-bandlimited graph signals. We prove this result in Section 2 but also show that there always exists a sampling distribution that yields optimal results. Finally, the reconstruction methods proposed in requires a partial diagonalisation of the Laplacian matrix, unlike ours. We also have much stronger recovery guarantees than the ones presented in , which are expected recovery guarantees.

2Notations and definitions

i.e., $\mathsf{U}_{k}$ is the restriction of $\mathsf{U}$ to its first $k$ vectors. This yields the following formal definition of a $k$-bandlimited signal.

Note we use ${\rm span}(\mathsf{U}_{k})$ in our definition of $k$-bandlimited signals to handle the case where the eigendecomposition is not unique. To avoid any ambiguity in the definition of $k$-bandlimited signals, we assume that $\bm{\lambda}_{k}\neq\bm{\lambda}_{k+1}$ for simplicity.

3Outline

In Section 2, we detail our sampling strategies and provide sufficient sampling conditions that ensure a stable embedding of $k$-bandlimited graph signals. We also prove that there always exists an optimal sampling distribution that ensures an embedding of $k$-bandlimited signals for $O(k\log(k))$ measurements. In Section 3, we propose decoders able to recover $k$-bandlimited signals from their samples. In Section 4, we explain how to obtain an estimation of the optimal sampling distribution quickly, without partial diagonalisation of the Laplacian matrix. In Section 5, we conduct several experiments on different graphs to test our methods. Finally, we conclude and discuss perspectives in Section 6.

Sampling k𝑘k-bandlimited signals

In this section, we start by describing how we select a subset of the nodes to sample $k$-bandlimited signals. Then, we prove that this sampling procedure stably embeds the set of $k$-bandlimited signals. We describe how to reconstruct such signals from these measurements in Section 3.

The subset of nodes $\Omega:=\{\omega_{1},\ldots,\omega_{m}\}$ used for sampling is constructed by drawing independently (with replacements) $m$ indices from the set $\{1,\ldots,n\}$ according to the probability distribution $\bm{p}$. We thus have

Note that we discuss the case of sampling without replacement in Section 2.3.

Let us pause for a moment and highlight few important facts. First, the sampling procedure allows each node to be selected multiple times. The number of measurements $m$ includes these duplications. In practice, one can sample each selected node only once and add these duplications “artificially” afterwards. Second, the set of nodes $\Omega$ needs to be selected only once to sample all $k$-bandlimited signals on $\mathcal{G}$. One does not need to construct a set $\Omega$ each time a signal has to be sampled. Third, note that the sampling procedure is so far completely independent of the graph $\mathcal{G}$. This is a non-adaptive sampling strategy.

for all $i\in\{1,\ldots,n\}$ and $j\in\{1,\ldots,m\}$. Note that $\bm{y}=\mathsf{M}\bm{x}$. In the next section, we show that, with high probability, $\mathsf{M}$ embeds the set of $k$-bandlimited signals for a number of measurements $m$ essentially proportional to $k\log(k)$ times a parameter called the graph weighted coherence.

2The space of k𝑘k-bandlimited signals is stably embedded

Similarly to many compressed sensing results, the number of measurements required to stably sample $k$-bandlimited signals will depend on a quantity, called the graph weighted coherence, that represents how the energy of these signals spreads over the nodes. Before providing the formal definition of this quantity, let us give an intuition of what it represents and why it is important.

characterises how much the energy of $\bm{\delta}_{i}$ is concentrated on the first $k$ Fourier modes. This ratio varies between and $1$. When it is equal to $1$, this indicates that there exists $k$-bandlimited signals whose energy is solely concentrated at the $i^{\text{th}}$ node; not sampling the $i^{\text{th}}$ node jeopardises the chance of reconstructing these signals. When this ratio is equal to , then no $k$-bandlimited signal has a part of its energy on the $i^{\text{th}}$ node; one can safely remove this node from the sampling set. We thus see that the quality of our sampling method will depend on the interplay between the sampling distribution $\bm{p}$ and the quantities $\left\|\mathsf{U}_{k}^{\intercal}\bm{\delta}_{i}\right\|_{2}$ for $i\in\{1,\ldots,n\}$. Ideally, we should have $\bm{p}_{i}$ large wherever $\left\|\mathsf{U}_{k}^{\intercal}\bm{\delta}_{i}\right\|_{2}$ is large and $\bm{p}_{i}$ small wherever $\left\|\mathsf{U}_{k}^{\intercal}\bm{\delta}_{i}\right\|_{2}$ is small. The interplay between $\bm{p}_{i}$ and $\left\|\mathsf{U}_{k}^{\intercal}\bm{\delta}_{i}\right\|_{2}$ is characterised by the graph weighted coherence.

The quantity $\left\|\mathsf{U}_{k}^{\intercal}\bm{\delta}_{i}\right\|_{2}$ is called the local graph coherence at node $i$.

Let us highlight two fundamental properties of $\nu_{\bm{p}}^{k}$. First, we have

Indeed, as the columns of $\mathsf{U}_{k}$ are normalised to $1$, we have

We are now ready to introduce our main theorem which shows that $m^{-1}\mathsf{M}\mathsf{P}^{-1/2}$ satisfies a restricted isometry property on the space of $k$-bandlimited signals.

Let $\mathsf{M}$ be a random subsampling matrix constructed as in (3) with the sampling distribution $\bm{p}$. For any $\delta,\epsilon\in(0,1)$, with probability at least $1-\epsilon$,

for all $\bm{x}\in{\rm span}(\mathsf{U}_{k})$ provided that

There are several important comments to make about the above theorem.

Third, as $(\nu_{\bm{p}}^{k})^{2}\geqslant k$, we need to sample at least $k$ nodes. Note that $k$ is also the minimum number of measurements that one must take to hope to reconstruct $\bm{x}\in{\rm span}(\mathsf{U}_{k})$.

The above theorem is quite similar to known compressed sensing results in bounded orthonormal systems . The proof actually relies on the same tools as the ones used in compressed sensing. However, in our case, the setting is simpler. Unlike in compressed sensing where the signal model is a union of subspaces, the model here is a single known subspace. In the proof, we exploit this fact to refine and tighten the sampling condition. In this simpler setting and thanks to our refined result, we can propose a sampling procedure that is always optimal in terms of the number of measurements.

for which $(\nu_{\bm{p}^{*}}^{k})^{2}=k$. The proof is simple. One just need to notice that $\sum_{i=1}^{n}\bm{p}^{*}_{i}=k^{-1}\;\sum_{i=1}^{n}\left\|\mathsf{U}_{k}^{\intercal}\bm{\delta}_{i}\right\|_{2}^{2}=k^{-1}\left\|\mathsf{U}_{k}\right\|_{\rm Frob}=k^{-1}k=1$ so that $\bm{p}^{*}$ is a valid probability distribution. Finally, it is easy to check that $(\nu_{\bm{p}^{*}}^{k})^{2}=k$. This yields the following corollary to Theorem 2.2.

Let $\mathsf{M}$ be a random subsampling matrix constructed as in (3) with the sampling distribution $\bm{p}^{*}$ defined in (6). For any $\delta,\epsilon\in(0,1)$, with probability at least $1-\epsilon$,

for all $\bm{x}\in{\rm span}(\mathsf{U}_{k})$ provided that

The sampling distribution $\bm{p}^{*}$ is optimal in the sense that the number of measurements needed to embed the set of $k$-bandlimited signals is essentially reduced to its minimum value. Note that, unlike Theorem 2.2 where the sampling is non-adaptive, the sampling distribution is now adapted to the structure of the graph and a priori requires the knowledge of a basis of ${\rm span}(\mathsf{U}_{k})$. We present a fast method that does not require the computation of a basis of ${\rm span}(\mathsf{U}_{k})$ to estimate $\bm{p}^{*}$ in Section 4.

It is important to mention that variable density sampling techniques are also popular in compressed sensing to reduce the sampling rate. We have been inspired by the works in this field to develop our sampling technique on graphs. In compressed sensing, the high efficiency of variable density sampling was first observed empirically in magnetic resonance imaging where the goal was to speed up the acquisition by reducing the amount of acquired data . Theoretical evidence of the efficiency of this technique then appeared in as well as in where additional measurement constraints and structured sparsity patterns are considered. There are similarities between the theoretical results existing in the compressed sensing literature and the ones presented in this work as we use similar proof techniques. However, we refine the proofs to take into account our specific signal model: bandlimited signals on graphs. One specificity of this setting is that there always exists a sampling distribution for which sampling $O(k\log(k))$ nodes is enough to capture all $k$-bandlimited signals. We recall that up to the log factor, one cannot hope to reduce the number of measurements much further. A second originality is that we can rapidly estimate this optimal distribution using fast filtering techniques on graphs (see Section 4).

Finally, we would like to highlight the similarity between the concept of local graph coherence and the concept of leverage scores of a matrix. Let $\mathsf{G}$ be a matrix and $\mathsf{V}$ be the left singular vectors of $\mathsf{G}$, the leverage scores related to the best rank-$r$ approximation of $\mathsf{G}$ are $l_{i}:=\left\|\mathsf{V}_{r}^{\intercal}\delta_{i}\right\|_{2}^{2}$, where $\mathsf{V}_{r}$ contains the $r$ eigenvectors of $\mathsf{G}$ with largest eigenvalues (see, e.g., for more details). The only difference between the leverage scores and the local graph coherences of $\mathsf{L}$ is thus that the latter involve the smallest eigenvalues while the leverage scores involve the largest eigenvalues. For the normalised graph Laplacian $\mathsf{L}=\mathsf{I}-\mathsf{D}^{-1/2}\mathsf{W}\mathsf{D}^{-1/2}$, one can notice that the leverage scores of $\mathsf{D}^{-1/2}\mathsf{W}\mathsf{D}^{-1/2}$ correspond exactly to the local graph coherences of $\mathsf{L}$. In machine learning, the leverage scores have been used, e.g., to improve the performance of randomised algorithms that solve overdetermined least-square problems or that compute a low-rank approximation of a given matrix . These algorithms work by building a sketch of the matrix of interest from a subset of its rows and/or columns. The leverage scores represent the importance of each row/columns in the dataset. The optimised sampling distribution for the rows/columns is then constructed by normalising the vector of leverage scores, as we do it here with the local graph coherence. Note also that fast algorithms to estimate the leverage scores have been developed in . In the future, it would be interesting to compare this method with ours, which explicitly uses the graph structure in the computations.

3Sampling without replacement

Let $\mathsf{M}$ be a random subsampling matrix constructed as in (3) with $\Omega$ built by drawing $m$ indices $\{\omega_{1},\ldots,\omega_{m}\}$ from $\{1,\ldots,n\}$ uniformly at random without replacement. For any $\delta,\epsilon\in(0,1)$, with probability at least $1-\epsilon$,

for all $\bm{x}\in{\rm span}(\mathsf{U}_{k})$ provided that

The attentive reader will notice that, unfortunately, the condition on $m$ is identical to the case where the sampling is done with replacement. This is because the theorem that we use to prove this result is obtained by “coming back” to sampling with replacement. Yet, we believe that it is still interesting to mention this result for applications where one wants to avoid any duplicated lines in the sampling matrix $\mathsf{M}$, which, for example, ensures that $\left\|\mathsf{M}\right\|_{2}=1$.

In the general case of non-uniform distributions, we are unfortunately not aware of any result allowing us to handle the case of a sampling without replacement. Yet it would be interesting to study this scenario more carefully in the future as sampling without replacement seems more natural for practical applications.

4Intuitive links between a graph’s structure and its local coherence

Recall that the local graph coherence at node $i$ reads $\left\|\mathsf{U}_{k}^{\intercal}\bm{\delta}_{i}\right\|_{2}$. We give here some examples showing how this quantity changes for different graphs.

Consider first the $d$-dimensional grid with periodic boundary conditions. In this case, the eigenvectors of its Laplacian are simply the $d$-dimensional classical Fourier modes. For simplicity, we suppose that $\lambda_{k}\neq\lambda_{k+1}$. In this case, one can show that the local coherence is independent of $i$: $\forall i\in\mathcal{V}\quad\left\|\mathsf{U}_{k}^{\intercal}\bm{\delta}_{i}\right\|_{2}=\sqrt{k/n}$. The optimal probability $\bm{p}^{*}$ is therefore uniform for the $d$-dimensional grid with periodic boundary conditions. Without the periodic boundary conditions, the optimal probability is mostly constant with an increase when going to the boundary nodes. An example in 1 dimension is shown in Fig. 4 for the path graph.

Let us now consider a graph made of $k$ disconnected components of size $n_{1},\ldots,n_{k}$. We have $\sum_{j=1}^{k}n_{j}=n$. Considering the combinatorial graph Laplacian, one can show that a basis of ${\rm span}{(\mathsf{U}_{k})}$ is the concatenation of the indicator vectors of each component. Moreover, ${\rm span}{(\mathsf{U}_{k})}$ is the eigenspace associated to the eigenvalue . The local coherence associated to node $i$ in component $j$ is $1/\sqrt{n_{j}}$, and the probability to choose this node reads $p_{i}^{*}=1/(kn_{j})$. If all components have the same size, i.e., $n_{1}~{}=~{}\ldots~{}=~{}n_{k}$, the optimal sampling is the uniform sampling. If the components have different sizes, the smaller is a component, the larger is the probability of sampling one of its nodes. With the optimal sampling distribution, each component is sampled with probability $1/k$, no matter the size of the component. The probability that each component is sampled at least once - a necessary condition for perfect recovery - is thus higher than when using uniform sampling. One may relax this strictly disconnected component example into a loosely defined community-structured graph, where $k$ sets of nodes (forming a partition of the $n$ nodes) are more connected with themselves than with the rest of the graph. In this case, one also expects that the probability to sample a node is inversely proportional to the size of the community it belongs to.

Signal recovery

In a situation where one knows a basis of ${\rm span}{(\mathsf{U}_{k})}$, the standard method to estimate $\bm{x}$ from $\bm{y}$ is to compute the best approximation to $\bm{y}$ from ${\rm span}(\mathsf{U}_{k})$, i.e., to solve

Note that we introduced a weighting by the matrix $\mathsf{P}^{-1/2}_{\Omega}$ in (7) to account for the fact that $m^{-1}\,\mathsf{P}^{-1/2}_{\Omega}\mathsf{M}=m^{-1}\,\mathsf{M}\mathsf{P}^{-1/2}$ satisfies the RIP, not $\mathsf{M}$ alone. The following theorem proves that the solution of (7) is a faithful estimation of $\bm{x}$.

Let $\bm{x}^{*}$ be the solution of Problem (7) with $\bm{y}=\mathsf{M}\bm{x}+\bm{n}$. Then,

We notice that in the absence of noise $\bm{x}^{*}=\bm{x}$, as desired. In the presence of noise, the upper bound on the error between $\bm{x}^{*}$ and $\bm{x}$ increases linearly with $\|\mathsf{P}^{-1/2}_{\Omega}\bm{n}\|_{2}$. For a uniform sampling, we have $\|\mathsf{P}^{-1/2}_{\Omega}\bm{n}\|_{2}=\sqrt{n}\;\|\bm{n}\|_{2}$. For a non-uniform sampling, we may have $\|\mathsf{P}^{-1/2}_{\Omega}\bm{n}\|_{2}\gg\sqrt{n}\;\left\|\bm{n}\right\|_{2}$ for some particular draws of $\Omega$ and noise vectors $\bm{n}$. Indeed, some weights $\bm{p}_{\omega_{i}}$ might be arbitrarily close to . Unfortunately, one cannot in general improve the upper bound in (8) as proved by the second part of the theorem with (9). Non-uniform sampling can thus be very sensitive to noise unlike uniform sampling. However, this is a worst case scenario. First, it is unlikely to draw an index $\omega_{i}$ where $\bm{p}_{\omega_{i}}$ is small by construction of the sampling procedure. Second,

so that $\|\mathsf{P}^{-1/2}_{\Omega}\bm{n}\|_{2}$ is not too large on average over the draw of $\Omega$. Furthermore, in our numerical experiments, we noticed that we have $\min_{i}\bm{p}_{i}=1/(\alpha^{2}\,n)$, where $\alpha>1$ is a small constantIn the numerical experiments presented below, we have $\alpha$ smaller or equal to $3$ in all cases tested with the optimal sampling distribution for the graphs presented in Fig. 1., for the optimal sampling distributions $\bm{p}=\bm{p}^{*}$ obtained in practice. This yields $\|\mathsf{P}^{-1/2}_{\Omega}\bm{n}\|_{2}\leqslant\alpha\,\sqrt{n}\,\left\|\bm{n}\right\|_{2}$, which shows that non-uniform sampling is just slightly more sensitive to noise than uniform sampling in practical settings, with the advantage of reducing the number of measurements. Non-uniform sampling is thus still a beneficial solution.

We have seen a first method to estimate $\bm{x}$ from its measurements. This method has however a major drawback: it requires the estimation of a basis of $\mathsf{U}_{k}$, which can be computationally very expensive for large graphs. To overcome this issue, we propose an alternative decoder which is computationally much more efficient. This algorithm uses techniques developed to filter graph signals rapidly. We thus briefly recall the principle of these filtering techniques.

2Fast filtering on graphs

According to the above definition, filtering a priori requires the knowledge of the matrix $\mathsf{U}$. To avoid the computation of $\mathsf{U}$, one can approximate the function $h$ by a polynomial

of degree $d$ and compute $\bm{x}_{\kappa}$, which will approximate $\bm{x}_{h}$. This computation can be done rapidly as it only requires matrix-vector multiplications with $\mathsf{L}$, which is sparse in most applications. Indeed,

We let the reader refer to for more information on this fast filtering technique.

Remark that $\kappa(\mathsf{L})=\mathsf{U}\,\kappa(\mathsf{\Lambda})\,\mathsf{U}^{\intercal}$.

3Efficient decoder

Instead of solving (7), we propose to estimate $\bm{x}$ by solving the following problem

We argue that solving (11) is computationally efficient because $\mathsf{L}$ is sparse in most applications. Therefore, any method solving (11) that requires only matrix-vector multiplications with $g(\mathsf{L})$ can be implemented efficiently, as it requires multiplications with $\mathsf{L}$ only (recall the definition of $g(\mathsf{L})$ in (10)). Examples of such methods are the conjugate gradient method or any gradient descend methods. Let us recall that one can find a solution to (11) by solving

The next theorem bounds the error between the original signal $\bm{x}$ and the solution of (11).

Let $\bm{x}^{*}$ be the solution of (7) with $\bm{y}=\mathsf{M}\bm{x}+\bm{n}$. Then,

where $\bm{\alpha}^{*}:=\mathsf{U}_{k}\mathsf{U}_{k}^{\intercal}\,\bm{x}^{*}$ and $\bm{\beta}^{*}:=(\mathsf{I}-\mathsf{U}_{k}\mathsf{U}_{k}^{\intercal})\,\bm{x}^{*}$.

In the above theorem, $\bm{\alpha}^{*}$ is the orthogonal projection of $\bm{x}^{*}$ onto ${\rm span}(\mathsf{U}_{k})$ and $\bm{\beta}^{*}$ onto the orthogonal complement of ${\rm span}(\mathsf{U}_{k})$. To obtain a bound on $\left\|\bm{x}^{*}-\bm{x}\right\|_{2}$, one can simply use the triangle inequality and the bounds (3.2) and (14).

In the presence of noise, for a fixed function $g$, the upper bound on the reconstruction error is minimised for a value of $\gamma$ proportional to $\|\mathsf{P}^{-1/2}_{\Omega}\bm{n}\|_{2}/\left\|\bm{x}\right\|_{2}$. To optimise the result further, one should seek to have $g(\bm{\lambda}_{k})$ as small as possible and $g(\bm{\lambda}_{k+1})$ as large as possible.

Decoder (11) has close links to several “decoders” used in the semi-supervised learning litterature that attempt to estimate the label of unlabeled data from a small number of labeled data, by supposing that the label functions are smooth either 1) in the data space or 2) in a suitable transformed space –using similarity kernels that define graphs modeling the underlying manifold for instance– , or 3) in both . In our work, smoothness is defined solely using the graph (case 2) which we suppose given; there is no equivalent of a data space (case 1) on which to define another smoothness constraint. Nevertheless, other types of smoothness could be considered instead of the Laplacian smoothness $\bm{z}^{\intercal}g(\mathsf{L})\bm{z}$. For instance, one could decide to use an $l_{1}$ penalisation of the graph difference operator, as in , to allow the signal to depart from the smoothness prior at some nodes. The closest semi-supervised framework to our naive decoder (7) is found in , where the authors constrain the solution to be in ${\rm span}{(\mathsf{U}_{k})}$ without specifying precisely the value of $k$; and the closest technique to our efficient decoder (11) is found in even though their cost function has an additional term of the form $\sum_{i\not\in\Omega}\bm{z}_{i}^{2}$ compared to ours. Another decoding method may be found in , where the authors have a similar cost function to ours. However, they work in the data space (case 1 above) and try to optimise this cost function directly in this space, i.e., without explicitly constructing and using a graph.

Even though we use similar decoders than in the semi-supervised learning literature, let us stress that an important difference is that we choose beforehand which nodes to sample/label. In this sense, our work also has connections with the literature in active learning , more precisely with the works that concentrate on the offline (all nodes to label are chosen from the start), single-batch (the nodes to sample are drawn simultaneously) selection problem , such as in . Yet another connection may be found in the area of Gaussian Random Fields . The originality of our method compared to these works comes from our specific smoothness model (${\rm span}{(\mathsf{U}_{k})}$) for which we devise an original sampling scenario which ensures stable reconstruction when coupled with the decoder (11) (see Theorem 3.2).

Estimation of the optimal sampling distribution

In this section, we explain how to estimate the optimal sampling distribution $\bm{p}^{*}$ efficiently. This distribution is entirely defined by the values $\left\|\mathsf{U}_{k}^{\intercal}\bm{\delta}_{i}\right\|_{2}^{2}$, $i=1,\ldots,n$ (see (6)). In order to be able to deal with large graphs and potentially large $k$, we want to avoid the computation of a basis of ${\rm span}(\mathsf{U}_{k})$ to estimate this distribution. Instead, we take another route that consists in filtering a small number of random signals. Note that the idea of filtering few random signals to estimate the number of eigenvalues of a Hermitian matrix in a given interval is already proposed and studied in . We show here that this technique can be used to estimate $\bm{p}^{*}$.

The estimation of the optimal sampling distribution is based on the following property. The $i^{\text{th}}$ entry of $\bm{r}_{b_{\bm{\lambda}_{k}}}$ is

and the mean of $(\bm{r}_{b_{\bm{\lambda}_{k}}})_{i}^{2}$ satisfies

This shows that $(\bm{r}_{b_{\bm{\lambda}_{k}}})_{i}^{2}$ is an unbiased estimation of $\|\mathsf{U}_{k}^{\intercal}\bm{\delta}_{i}\|_{2}^{2}$, the quantity we want to evaluate. Therefore, a possibility to estimate the optimal sampling distribution consists in filtering $L$ random signals $\bm{r}^{1},\ldots,\bm{r}^{L}$ with the same distribution as $\bm{r}$ and average $(\bm{r}_{b_{\bm{\lambda}_{k}}}^{1})_{i}^{2},\ldots,(\bm{r}_{b_{\bm{\lambda}_{k}}}^{L})_{i}^{2}$ for each $i\in\{1,\ldots,n\}$. The next theorem shows that if $\bm{\lambda}_{k}$ is known, then $L\geqslant O(\log(n))$ random vectors are sufficient to have an accurate estimation of $\|\mathsf{U}_{k}^{\intercal}\bm{\delta}_{i}\|_{2}^{2}$.

for all $i\in\{1,\ldots,n\}$, provided that

The above theorem indicates that if $\lambda\in[\bm{\lambda}_{k},\bm{\lambda}_{k+1})$ and $\hat{e}$ is null, then $\sum_{l=1}^{L}\;(\bm{r}_{c_{\lambda_{k}}}^{l})_{i}^{2}$ estimates $\left\|\mathsf{U}_{k}^{\intercal}\bm{\delta}_{i}\right\|_{2}^{2}$ with an error at most $\delta$ on each entry $i\in\{1,\ldots,n\}$. Recalling that the optimal sampling distribution has entries

approximates the optimal sampling distribution. If we know $\bm{\lambda}_{k}$ and $\bm{\lambda}_{k+1}$, we can thus approximate $\bm{p}^{*}$. In order to complete the method, we now need a solution to estimate $\bm{\lambda}_{j}$ with $j=k$ or $j=k+1$.

𝑘1\lambda_{k+1} Let $\lambda\in(0,\bm{\lambda}_{n})$. Theorem 4.1 shows that, with probability $1-\epsilon$,

when using the filter $b_{\lambda}$. Noticing that

as the columns of $\mathsf{U}$ are normalised, yields

In other words, the total energy of the filtered signals is tightly concentrated around $j^{*}$, which is the largest integer such that $\bm{\lambda}_{j^{*}}\leqslant\lambda$. Therefore, the total energy of the filtered signals provides an estimation of the number of eigenvalues of $\mathsf{L}$ that are below $\lambda$.

Using this phenomenon, one can obtain, by dichotomy, an interval $(\underline{\lambda},\bar{\lambda})$ such that $k-1$ eigenvalues are below $\underline{\lambda}$ and $k$ eigenvalues are below $\bar{\lambda}$ and thus obtain an estimation of $\bm{\lambda}_{k}$. The same procedure can be used to estimate $\bm{\lambda}_{k+1}$. Note that we cannot filter the signals using an ideal low-pass filter in practice, so that an additional error will slightly perturb the estimation.

3The complete algorithm

Experiments

In this section, we run several experiments to illustrate the above theoretical findings. First we show how the sampling distribution affects the number of measurements required to ensure that the RIP holds. Then, we show how the reconstruction quality is affected with the choice of $g$ and $\gamma$ in (11).

All our experiments are done using five different types of graph, all available in the GSP toolbox and presented in Fig. 1. We use a) different community-type graphs of size $n=1000$, b) the graph representing the Minnesota road network of size $n=2642$, c) the graph of the Stanford bunny of size $n=2503$, d) the unweighted path graph of size $n=1000$ and e) a binary tree of depth $9$ and size $n=1023$. We recall that each node in the unweighted path graph is connected to its left and right neighbours with weight $1$, except for the two nodes at the boundaries which have only one neighbour. We use the combinatorial Laplacian in all experiments. All samplings are done in the conditions of Theorem 2.2, i.e., with replacement. Finally, the reconstructions are obtained by solving (12) using the mldivide function of Matlab. For the graphs and functions $g$ considered, we noticed that it was faster to use this function than solving (12) by conjugate gradient.

We conduct a first set of experiments using five types of community graph, denoted by $\mathcal{C}_{1},\ldots,\mathcal{C}_{5}$. They all have $10$ communities. To study the effect of the size of the communities on the sampling distribution, we choose to build these graphs with $9$ communities of (approximately) equal size and reduce the size of last community:

the graphs of type $\mathcal{C}_{1}$ have $10$ communities of size $100$;

the graphs of type $\mathcal{C}_{2}$ have $1$ community of size $50$, $8$ communities of size $105$, and $1$ community of size $110$;

the graphs of type $\mathcal{C}_{3}$ have $1$ community of size $25$, $8$ communities of size $108$, and $1$ community of size $111$;

the graphs of type $\mathcal{C}_{4}$ have $1$ community of size $17$, $8$ communities of size $109$, and $1$ community of size $111$;

the graphs of type $\mathcal{C}_{5}$ have $1$ community of size $13$, $8$ communities of size $109$, and $1$ community of size $115$.

for different numbers of measurements $m$. Note that to compute $\underline{\delta}_{10}$, one just needs to notice that

For the uniform distribution $\bm{\pi}$, the first figure from the left in Fig. 2 indicates the value of $(\nu_{\bm{\pi}}^{10})^{2}(j)$, $j=1,\ldots,5$ for the five different types of graph. We have $(\nu_{\bm{\pi}}^{10})^{2}(1)\leqslant\ldots\leqslant(\nu_{\bm{\pi}}^{10})^{2}(5)$ and $m_{1,\bm{\pi}}^{*}\leqslant m_{2,\bm{\pi}}^{*}\leqslant\ldots\leqslant m_{5,\bm{\pi}}^{*}$, in accordance with Theorem 2.2.

For the optimal sampling distribution $\bm{p}^{*}$, we have $(\nu_{\bm{p}^{*}}^{10})^{2}=10$. Therefore $m_{j,\bm{p}^{*}}^{*}$ must be identical for all graph-types, as observed in the second panel of Fig. 2.

1.2Using the Minnesota and bunny graphs

To confirm the results observed above, we repeat the same experiments but using four other graphs: the Minnesota, the bunny and path graphs, and the binary tree. For the first three graphs, the experiments are performed for $k$-bandlimited signals with band-limits $10$ and $100$, i.e., we compute $\underline{\delta}_{k}$ - defined as in (16) - with $\mathsf{U}_{k}$. For the binary tree, we set the band-limits at $16$ and $64$. These choices are due to the fact that some eigenvalues have a multiplicity larger than $1$ for this last graph. These choices ensure that $\bm{\lambda}_{k}<\bm{\lambda}_{k+1}$, as required in our assumptions.

We present in Fig. 3 the probability that $\underline{\delta}_{k}$ is less than $0.995$, estimated over $500$ draws of $\mathsf{M}$, as a function of $m$.

1.3Examples of optimal and estimated sampling distributions

For illustration, we present some examples of sampling distributions in Fig. 4 for five of the graphs used above. The top panels in Fig. 4 show the optimal sampling distribution computed with $\mathsf{U}_{k}$. The bottom panels show the estimated sampling distribution obtained with Algorithm 1.

It is interesting to notice that for the path graph, the optimal sampling distribution is essentially constant except for the nodes at the boundaries that are less connected than the other nodes and need to be sampled with higher probability. For the binary tree, the probability of sampling a node is only determined by its depth in the tree - all nodes at a given depth are equally important - and the deeper the node is, the higher the probability of sampling this node should be - it is easier to predict the value at one node from the values bore by its children than its parents.

2Reconstruction of k𝑘k-bandlimited signals

We present the mean reconstruction errors obtained in the absence of noise on the measurements in Fig. 5. In this set of experiments, we reconstruct the signals using $g(\mathsf{L})=\mathsf{L}$, then $g(\mathsf{L})=\mathsf{L}^{2}$, and finally $g(\mathsf{L})=\mathsf{L}^{4}$. Before describing these results, we recall that the ratio $g(\bm{\lambda}_{10})/g(\bm{\lambda}_{11})$ decreases as the power of $\mathsf{L}$ increases. We observe that all reconstruction errors, $\left\|\bm{x}^{*}-\bm{x}\right\|_{2}$, $\left\|\bm{\alpha}^{*}-\bm{x}\right\|_{2}$ and $\left\|\bm{\beta}^{*}\right\|_{2}$ decrease when the ratio $g(\bm{\lambda}_{k})/g(\bm{\lambda}_{k+1})$ in the range of small $\gamma$, as predicted by the upper bounds on these errors in Theorem 3.2.

We present the mean reconstruction errors obtained in the presence of noise on the measurements in Fig. 6. In this set of experiments, we reconstruct the signals using $g(\mathsf{L})=\mathsf{L}^{4}$. As expected the best regularisation parameter $\gamma$ increases with the noise level.

3Illustration: sampling of a real image

We finish this experimental section with an example of image sampling using the developed theory.

where $n=190816$. Each column of $\mathsf{X}$ represents a color patch of the original image at a given position.

We present in Fig. 7 the sampled image, where all non-sampled pixels appear in black. We remark that the regions where many patches are similar (sky, lake, snow) are very sparsely sampled. This can be explained as follows. The patches in such a region being all similar, one can fill this region by copying a single representative patch. In practice this is done via the Laplacian matrix, which encodes the similarities between the patches, by solving (11).

Conclusion and perspectives

We proposed two efficient sampling procedures for $k$-bandlimited signals defined on the nodes of a graph $\mathcal{G}$. The performance of these sampling techniques is governed by the graph weighted coherence, which characterises the interaction between the sampling distribution and the localisation of the first $k$ Fourier modes over the nodes of $\mathcal{G}$. For regular graph with non-localised Fourier modes and a uniform sampling distribution, we proved that $O(k\log{k})$ samples are sufficient to embed the set of $k$-bandlimited signals. For arbitrary graphs, uniform sampling might perform very poorly. In such cases, we proved that it is always possible to adapt the sampling distribution to the structure of the graph and reach optimal sampling conditions. We designed an algorithm to estimate the optimal sampling distribution rapidly. Finally, we proposed an efficient decoder that provides accurate and stable reconstruction of $k$-bandlimited signals from their samples.

We believe that the sampling method developed in this work can be used to speed up computations in multiple applications using graph models. Let us take the example of the fast robust PCA method proposed in . In this work, the authors consider the case where one has access to two graphs $\mathcal{G}_{1}$ and $\mathcal{G}_{2}$ that respectively model the similarities between the rows and the columns of a matrix $\mathsf{X}$. In this context, they propose an optimisation technique that provides a low-rank approximation of $\mathsf{X}$. We denote this low-rank approximation by $\mathsf{X}^{*}$. The intuition is that the left singular vectors and the right singular vectors of $\mathsf{X}^{*}$ live respectively in the span of the first eigenvectors of $\mathsf{L}_{1}$ and $\mathsf{L}_{2}$, the Laplacians associated to $\mathcal{G}_{1}$ and $\mathcal{G}_{2}$. Therefore, the singular vectors of $\mathsf{X}^{*}$ can be drastically subsampled using our sampling method. The low-rank matrix $\mathsf{X}^{*}$ can be reconstructed from a subset of its rows and columns. Instead of estimating $\mathsf{X}^{*}$ from the entire matrix $\mathsf{X}$, one could thus first reduce the dimension of the problem by selecting a small subset of the rows and columns of $\mathsf{X}$.

In semi-supervised learning, a small subset of nodes are labeled and the goal is to infer the label of all nodes. Advances in sampling of graph signals give insight on which nodes should be preferentially observed to infer the labels on the complete graphs. Similarly, in spectral graph clustering, cluster assignments are well approximated by $k$-bandlimited signals and can therefore be heavily subsampled. This leaves the possibility to initially cluster a small subset of the nodes and infer the clustering solution on the complete graphs afterwards, as we have proposed in .

Sensor networks provide other applications of our sampling methods. Indeed, if signals measured by a network of sensors are smooth, one can deduce beforehand from the structure of the network which sensors to sample in priority in an active sampling strategy, using the optimal or estimated sampling distribution.

Appendix A - Proof of the theorems in Section 2

We start with the proof of Theorem 2.2. For this proof, we need the following result obtained by Tropp in .

Consider a finite sequence $\{\mathsf{X}_{i}\}$ of independent, random, self-adjoint, positive semi-definite matrices of dimension $d\times d$. Assume that each random matrix satisfies

We also need the following facts. For all $\delta\in$, we have

As the $i^{\text{th}}$ row-vector of $\mathsf{M}\mathsf{P}^{-1/2}\mathsf{U}_{k}$ is $\bm{\delta}_{\omega_{i}}^{\intercal}\mathsf{U}_{k}/\sqrt{\bm{p}_{\omega_{i}}}$, we have

The expected value of each $\mathsf{X}_{i}$ is

Furthermore, for all $i=1,\ldots,n$, we have

Lemma A.1 yields, for any $\delta\in(0,1)$,

Therefore, for any $\delta\in(0,1)$, we have, with probability at least $1-\epsilon$,

for all $\bm{x}\in{\rm span}(\mathsf{U}_{k})$, terminates the proof. ∎

The proof of Theorem 2.4 is based on the following results, also obtained by Tropp.

Let $\mathcal{X}$ be a finite set of positive-semidefinite matrices of dimension $d\times d$, and suppose that

Sample $\{\mathsf{X}_{1},\ldots,\mathsf{X}_{l}\}$ uniformly at random from $\mathcal{X}$ without replacement. Compute

Using Lemma A.1, one can notice that the above probability bounds would be identical if the matrices $\{\mathsf{X}_{1},\ldots,\mathsf{X}_{l}\}$ were sampled uniformly at random from $\mathcal{X}$ with replacement. It is thus not necessary to detail the complete proof which is entirely similar to the one of Theorem 2.2, at the exception of the sampling procedure.

Appendix B - Proof of the theorems in Section 3

We recall that $\bm{x}^{*}$ is a solution to (7). By optimality of $\bm{x}^{*}$, we have

for any $\bm{z}\in{\rm span}(\mathsf{U}_{k})$. In particular for $\bm{z}=\bm{x}$, we obtain

Then, the triangle inequality and (4) yields

In the second step, we used the fact that $\mathsf{M}\mathsf{P}^{-1/2}=\mathsf{P}_{\Omega}^{-1/2}\mathsf{M}$. Combining (18) and (19) directly yields (8), the first bound in Theorem 3.1.

To prove the second bound, let us choose $\bm{n}_{0}=\mathsf{M}\bm{z}_{0}$ with $\bm{z}_{0}\in{\rm span}(\mathsf{U}_{k})$. Therefore, $\bm{y}=\mathsf{M}(\bm{x}+\bm{z}_{0})$ and $\bm{x}^{*}=\bm{x}+\bm{z}_{0}$ is an obvious solution to (7) in this case. To finish the proof, we use (4) which yields

As $\bm{x}^{*}$ is a solution to (11), we have

Choosing $\bm{z}=\bm{x}$ in (20) and using the facts that $\bar{\mathsf{U}}_{k}^{\intercal}\bm{\alpha}^{*}=\bm{0}$, $\mathsf{U}_{k}^{\intercal}\bm{\beta}^{*}=\bm{0}$, $\bar{\mathsf{U}}_{k}^{\intercal}\bm{x}=\bm{0}$, and that $g(\mathsf{L})=\mathsf{U}\,g(\mathsf{L})\,\mathsf{U}^{\intercal}$, we obtain

where we used the fact that $\left\|\bar{\mathsf{U}}_{k}^{\intercal}\bm{\beta}^{*}\right\|_{2}=\left\|\bm{\beta}^{*}\right\|_{2}$ and $\left\|\mathsf{U}_{k}^{\intercal}\bm{x}\right\|_{2}=\left\|\bm{x}\right\|_{2}$. As the left hand side of the last inequality is a sum of two positive quantities, we also have

Inequality (22) proves (14), the second inquality in Theorem 3.2. It remains to prove (3.2). To prove this inequality, we continue by using (4), which yields

Finally, combining (21), (22) and (23) gives

Appendix C - Proof of the theorem in Section 4

We use the classical technique to prove the Johnson-Lindenstrauss lemma (see, e.g., ).

Each filtered signal $\bm{r}_{\hat{c}_{\lambda}}^{l}$, $l\in\{1,\ldots,L\}$, satisfies

where $\mathsf{C}_{\lambda}:={\rm diag}(\hat{c}_{\lambda}(\bm{\lambda}_{1}),\ldots,\hat{c}_{\lambda}(\bm{\lambda}_{n}))$. Let $i$ be fixed for the moment. We have

The expected value of this sum is $\left\|\mathsf{C}_{\lambda}\mathsf{U}^{\intercal}\bm{\delta}_{i}\right\|_{2}^{2}$. Indeed,

This is a sum of $L$ independent centered random variables. Furthermore, as each $\bm{r}^{l}$ is a zero-mean Gaussian random vector with covariance matrix $L^{-1}\,\mathsf{I}$, the variables $(\bm{r}_{\hat{c}_{\lambda}}^{l})_{i}$ are subgaussian with subgaussian bounded by $C\,L^{-1/2}\,\left\|\mathsf{C}_{\lambda}\mathsf{U}^{\intercal}\bm{\delta}_{i}\right\|_{2}$, where $C\geqslant 1$ is an absolute constant. We let the reader refer to, e.g., for more information on the definition and properties of subgaussian random variables. Using Lemma $5.14$ and Remark $5.18$ in , one can prove that each summand of $X$ is a centered subexponential random variable with subexponentinal norm bounded by $4C^{2}\,L^{-1}\,\left\|\mathsf{C}_{\lambda}\mathsf{U}^{\intercal}\bm{\delta}_{i}\right\|_{2}^{2}$. Corollary $5.17$ in shows that there exists an absolute contant $c>0$ such that

for all $t\in(0,4C^{2}\,L^{-1}\,\left\|\mathsf{C}_{\lambda}\mathsf{U}^{\intercal}\bm{\delta}_{i}\right\|_{2}^{2})$, or, equivalently, that

This proves that, with probability at least $1-\epsilon$,

for all $i\in\{1,\ldots,n\}$, provided that

To finish the the proof, one just needs to remark that

by definition of $\hat{c}_{\lambda}$ (see (15)) and use the triangle inequality. ∎

References