Structured sampling and fast reconstruction of smooth graph signals

TL;DR

This paper introduces a structured sampling method for smooth signals on graphs, enabling efficient and stable reconstruction by leveraging group-based sampling strategies and nearly piecewise constant assumptions.

Contribution

It proposes a novel group-based sampling strategy with theoretical guarantees and fast approximation methods for reconstructing smooth graph signals.

Findings

Sampling $O(k ext{log}(k))$ groups suffices for stable embedding.

Reconstruction speed is improved by reducing the problem dimension.

Numerical experiments validate the theoretical results and demonstrate practical applications.

Abstract

This work concerns sampling of smooth signals on arbitrary graphs. We first study a structured sampling strategy for such smooth graph signals that consists of a random selection of few pre-defined groups of nodes. The number of groups to sample to stably embed the set of $k$-bandlimited signals is driven by a quantity called the \emph{group} graph cumulative coherence. For some optimised sampling distributions, we show that sampling $O(k\log(k))$ groups is always sufficient to stably embed the set of $k$-bandlimited signals but that this number can be smaller -- down to $O(\log(k))$ -- depending on the structure of the groups of nodes. Fast methods to approximate these sampling distributions are detailed. Second, we consider $k$-bandlimited signals that are nearly piecewise constant over pre-defined groups of nodes. We show that it is possible to speed up the reconstruction of such signals by reducing drastically the dimension of the vectors to reconstruct. When combined with the proposed structured sampling procedure, we prove that the method provides stable and accurate reconstruction of the original signal. Finally, we present numerical experiments that illustrate our theoretical results and, as an example, show how to combine these methods for interactive object segmentation in an image using superpixels.

Equations

$${\mathsf{U}}_k:=\left({\mathbf{u}}_1,\dots ,{\mathbf{u}}_k\right)\in {\mathbb{R}}^{n\times k}.$$

$$\mathbf{x}={\mathsf{U}}_k\mathbf{\eta }.$$

$${\cup }_{\ell =1}^N{\mathcal{N}}_{\ell }$$

$$\displaystyle\mathsf{N}^{(\ell)}_{ij}:=\left\{\begin{array}[]{ll}1&\text{for }j=n_{i}^{(\ell)},\\ 0&\text{otherwise}.\end{array}\right.$$

$$\sum _{\ell =1}^N{{\mathsf{N}}^{(\ell )}}^{⊺}{\mathsf{N}}^{(\ell )}=\mathsf{I}.$$

$$\displaystyle\mathsf{N}^{(\ell)}_{ij}:=\left\{\begin{array}[]{ll}\beta_{n_{i}^{(\ell)}}^{-1/2}&\text{for }j=n_{i}^{(\ell)},\\ 0&\text{otherwise},\end{array}\right.$$

$$\mathbb{P}({\omega }_j=\ell )={\mathbf{p}}_{\ell },\forall j\in \{1,\dots ,s\}\text{ and }\forall \ell \in \{1,\dots ,N\}.$$

$$m:=\sum _{j=1}^s\left|{\mathcal{N}}_{{\omega }_j}\right|.$$

$$\displaystyle\mathsf{M}:=\left(\begin{array}[]{c}\mathsf{N}^{(\omega_{1})}\\ \vdots\\ \mathsf{N}^{(\omega_{s})}\end{array}\right),$$

$$\mathbf{y}:=\mathsf{M}\mathbf{x}.$$

$${\mathsf{P}}^{(\ell )}:={\mathbf{p}}_{\ell }^{-1/2}\mathsf{I}\in {\mathbb{R}}^{|{\mathcal{N}}_{\ell }|\times |{\mathcal{N}}_{\ell }|},$$

$$\mathsf{P}:=\mathrm{diag}\left({\mathsf{P}}^{({\omega }_1)},\dots ,{\mathsf{P}}^{({\omega }_s)}\right).$$

$${\Vert {\mathsf{N}}^{(\ell )}{\mathsf{U}}_k\Vert }_2=\underset{\mathbf{\eta }\in {\mathbb{R}}^k:{\Vert \mathbf{\eta }\Vert }_2=1}{\sup }{\Vert {\mathsf{N}}^{(\ell )}{\mathsf{U}}_k\mathbf{\eta }\Vert }_2$$

$${\nu }_{\mathbf{p}}:=\underset{1⩽\ell ⩽N}{\max }\left\{{\mathbf{p}}_{\ell }^{-1/2}{\Vert {\mathsf{N}}^{(\ell )}{\mathsf{U}}_k\Vert }_2\right\}.$$

$$1$$

$$=\underset{1⩽\ell ⩽N}{\max }\left\{\frac{{\Vert {\mathsf{N}}^{(\ell )}{\mathsf{U}}_k\mathbf{\eta }\Vert }_2^2}{{\mathbf{p}}_{\ell }}\right\}⩽\underset{1⩽\ell ⩽N}{\max }\left\{\frac{{\Vert {\mathsf{N}}^{(\ell )}{\mathsf{U}}_k\Vert }_2^2}{{\mathbf{p}}_{\ell }}\right\}={\nu }_{\mathbf{p}}^2.$$

$$(1-\delta ){\Vert \mathbf{x}\Vert }_2^2⩽\frac{1}{s}{\Vert \mathsf{PM}\mathbf{x}\Vert }_2^2⩽(1+\delta ){\Vert \mathbf{x}\Vert }_2^2$$

$$s⩾\frac{3}{{\delta }^2}{\nu }_{\mathbf{p}}^2\log \left(\frac{2k}{\xi }\right).$$

$${\mathbf{p}}_{\ell }^{\ast }:=\frac{{\Vert {\mathsf{N}}^{(\ell )}{\mathsf{U}}_k\Vert }_2^2}{{\sum }_{{\ell }^{\prime }=1}^N{\Vert {\mathsf{N}}^{({\ell }^{\prime })}{\mathsf{U}}_k\Vert }_2^2},\text{ for all }\ell \in \{1,\dots ,N\},$$

$${\nu }_{{\mathbf{p}}^{\ast }}^2=\sum _{\ell =1}^N{\Vert {\mathsf{N}}^{(\ell )}{\mathsf{U}}_k\Vert }_2^2.$$

$${\nu }_{{\mathbf{p}}^{\prime }}^2⩾{{\mathbf{p}}_{{\ell }^{\prime }}^{\prime }}^{-1}{\Vert {\mathsf{N}}^{({\ell }^{\prime })}{\mathsf{U}}_k\Vert }_2^2>{{\mathbf{p}}_{{\ell }^{\prime }}^{\ast }}^{-1}{\Vert {\mathsf{N}}^{({\ell }^{\prime })}{\mathsf{U}}_k\Vert }_2^2={\nu }_{{\mathbf{p}}^{\ast }}^2,$$

$${\nu }_{{\mathbf{p}}^{\ast }}^2=\sum _{\ell =1}^N{\Vert {\mathsf{N}}^{(\ell )}{\mathsf{U}}_k\Vert }_2^2⩽\sum _{\ell =1}^N{\Vert {\mathsf{N}}^{(\ell )}{\mathsf{U}}_k\Vert }_F^2=k.$$

$${\nu }_{{\mathbf{p}}^{\ast }}^2=\sum _{\ell =1}^N{\Vert {\mathsf{N}}^{(\ell )}{\mathsf{U}}_k\Vert }_2^2⩽N,$$

$${\bar{\nu }}_{\mathbf{p}}:=\underset{1⩽\ell ⩽N}{\max }\left\{{\mathbf{p}}_{\ell }^{-1/2}{\Vert {\mathsf{N}}^{(\ell )}{\mathsf{U}}_k\Vert }_F\right\}.$$

$$k$$

$$s⩾\frac{3}{{\delta }^2}{\bar{\nu }}_{\mathbf{p}}^2\log \left(\frac{2k}{\xi }\right).$$

$${\mathbf{q}}_{\ell }^{\ast }:=\frac{{\Vert {\mathsf{N}}^{(\ell )}{\mathsf{U}}_k\Vert }_F^2}{k},\text{ for all }l\in \{1,\dots ,N\},$$

$${\bar{\nu }}_{{\mathbf{q}}^{\ast }}^2=k.$$

$$\underset{\mathbf{z}\in {\mathbb{R}}^n}{\min }{\Vert \mathsf{P}(\mathsf{M}\mathbf{z}-\mathbf{y})\Vert }_2^2+\gamma {\mathbf{z}}^{⊺}g(\mathsf{L})\mathbf{z},$$

$$g(\mathsf{L}):=\sum _{i=0}^d{\alpha }_i{\mathsf{L}}^i=\sum _{i=0}^d{\alpha }_i\mathsf{U}{\mathsf{\Lambda }}^i{\mathsf{U}}^{⊺},$$

$${\mathbf{a}}^{(\ell )}:=\frac{1^{⊺}{\mathsf{N}}^{(\ell )}}{{\left|{\mathcal{N}}_{\ell }\right|}^{1/2}},$$

Full text

G. Puy and P. Pérez are with Technicolor, 975 Avenue des Champs Blancs, 35576 Cesson-Sévigné, France..

1Introduction

This work was initially inspired by studying developments in edit propagation for interactive image or video manipulations, where the goal is to propagate operations made by a user in some parts of the image to the entire image, e.g., propagate foreground/background scribbles for object segmentation [1, 2, 3, 4, 5] or propagate manual color/tone modifications [3, 4, 6, 7, 5]. First, a graph $\mathcal{G}$ modelling the similarities between the pixels of the image is built. Second, the user-edits specify values at some nodes (pixels) of $\mathcal{G}$. Finally, the complete signal is estimated by assuming that it is smooth on $\mathcal{G}$. The quality of the propagation depends on the structure of $\mathcal{G}$ and on the location of annotated nodes. If a part of the image is weakly connected to the rest, the user-edits do not propagate well to this region unless this region is edited directly. Therefore, highlighting beforehand which regions or groups of nodes are important to edit to ensure a good propagation would be a useful feature to facilitate user interactions. Furthermore, designing a fast reconstruction/propagation method is also important so that the user can visualise immediately the effect of his inputs. To this end, “superpixels” – small groups of connected pixels where the image varies only little – are sometimes used among other things to speed up computations, e.g., in [3]. We address these problems from a graph signal processing point-of-view [8]. More precisely, we view them as a sampling problem where we should select a structured set of nodes (regions to edit), “measure” the signal on this set (user-edits), and reconstruct the signal on the entire graph. We believe that the sampling strategy and the fast reconstruction method proposed here provide useful tools to optimize user’s inputs and accelerate computations.

As mentionned above, the user-edits are propagated to the entire image by assuming that the global signal is smooth on $\mathcal{G}$. In the context of graph signal processing, smooth or $k$-bandlimited signal is a widely used and studied model. Several sampling methods have been designed to sample such signals. Pesenson introduced the notion of uniqueness set (of nodes) for $k$-bandlimited graph signals in [9, 10]. Two different $k$-bandlimited signals are also necessarily different when restricted to a uniqueness set. Therefore, one can sample all $k$-bandlimited signals on a uniqueness set. Then, Anis et al., [11, 12] and Chen et al., [13, 14] proved that one can always find a uniquess set of $k$ nodes to sample all $k$-bandlimited signals. Finding this set is however computationally expensive. In [15], graph spectral proxies are used to find such a set more efficiently. Yet the combinatorial problem that needs to be solved to find such a set makes the method still difficult to use for large graphs. Other authors used the idea of random sampling to be able to handle large graphs [14, 16, 17]. Recently, Puy et al. [17] proved that there always exists a random sampling strategy for which sampling $O(k\log(k))$ nodes is sufficient to stably embed the set of bandlimited signals. They also designed a fast and scalable algorithm to estimate the optimal sampling distribution.

In this paper, we study first a random sampling strategy for $k$-bandlimited signals where we sample few groups of nodes instead of sampling individual nodes. We introduce the concept of local group graph coherence that quantifies the importance of sampling each group. Second, in order to build a fast reconstruction technique for $k$-bandlimited signals, we use the intuition that a smooth graph signal is a signal that varies slowly from one node to its connected nodes. If we group few connected nodes together, we usually expect a bandlimited signal to be essentially constant on this set of nodes; as long as we do not group together too many weakly connected nodes. We propose here to use this property to accelerate the reconstruction of such $k$-bandlimited signals, i.e., $k$-bandlimited signals nearly piecewise constant over (pre-defined) groups of nodes. When combined with the proposed sampling technique, we prove that this fast method provides stable and accurate reconstructions of the signals of interest. Finally, we illustrate how to use these results for interactive object segmentation in an image, with required node groups being superpixels.

1.1Contributions

The random sampling strategy that we propose generalises the method proposed in [17]. We use here a structured sampling strategy. Let us already acknowledge here that such strategies are also studied in the field of compressed sensing [18, 19, 20, 21] and that some of our solutions are directly inspired by these works.

First, in this structured sampling setting, we show that the number of groups to sample is directly linked to a quantity called the group graph cumulative coherence. This quantity generalises the concept of graph cumulative coherence introduced in [17] and characterises how much the energy of $k$-bandlimited can stay concentrated in each group of nodes.

Second, we can then choose to sample the groups non-adaptively or to optimise the sampling distribution to minimise the number of groups to sample. With this optimised sampling distribution, our result shows that, in the worst case, sampling $O(k\log(k))$ groups of nodes is sufficient to ensure the reconstruction of all $k$-bandlimited signals. As each group can contain many nodes, we might have to sample a large number of nodes. This is the potential price to pay when sampling the nodes by groups. Fortunately, a smaller number of groups – down to $O(\log(k))$ – might already be sufficient if the groups are well designed.

Third, we describe a method to estimate the optimal sampling distribution without the need of computing the graph Fourier matrix and which is thus able to handle large graphs.

Fourth, estimating the optimal sampling distribution may still be too slow when a large number of groups are involved. We thus also present a sufficient recovery condition that involves a relaxed version of group graph cumulative coherence. The sampling distribution that minimises this relaxed coherence is fast to estimate for large graphs and large number of groups. With this sampling strategy, we prove that sampling $O(k\log(k))$ groups is always sufficient to ensure the reconstruction of all $k$-bandlimited signals. This strategy is mainly interesting at small $k$.

Finally, we propose a fast reconstruction method for $k$-bandlimited signals that are also nearly piecewise constant over pre-defined groups of nodes. We show that we can reduce drastically the dimension of the reconstruction problem for such signals. When the above sampling procedure is used to sample them, we prove that the proposed method provides accurate and stable recovery of this type of signals.

1.2Applications

The proposed sampling methods can have interest in several applications. For example, if one needs to develop sensors in a large scale network, it might be easier to deploy and install these sensors at nodes around few spatial locations instead of scattering them all over the network. Finding the best regions to monitor is thus important. In a social network, one might be interested in monitoring a signal defined over communities and thus should find which groups of users are the most important to sample. In semi-supervised learning, it may be easier to label jointly some similar nodes than individual nodes. Let us now detail such an example for semi-supervised classification.

The task we consider is interactive object segmentation in an image where, similarly to what is done in [22] for instance, we build a graph $\mathcal{G}$ that models similarities between the pixels of the image, ask the user to label some regions depending on whether they are part of the object or not, diffuse the result on the complete image by supposing that indicator vectors of the object is smooth on $\mathcal{G}$. To propose the regions to label, we view this segmentation problem as a sampling problem of a smooth signal on $\mathcal{G}$ where the regions to label are chosen to ensure the recovery of a $k$-bandlimited signal. To choose the regions to label, we start by partitioning the image into superpixels, e.g., with SLIC technique [23]. We see in Fig. 1 that the superpixels divide the image into homogeneous regions. The superpixels follow the edges and most of them thus belong to either the tiger or the background but rarely both. One interest of dividing the image into superpixels is that it facilitates user interactions. It is easier to determine if a superpixel belongs to the object of interest than if a pixel belongs to the tiger, especially at the boundaries, as recently exploited for segmentation on touch-screens [24]. Using the sampling method that we developed we can propose to the user a small number of superpixels to label. Furthermore, this proposition is adapted to the structure of the graph. Another advantage of using superpixels is that the indicator function of the object, beyond being smooth on $\mathcal{G}$, is also approximately piecewise constant on the superpixels. We can thus use our reconstruction method to estimate rapidly the segmentation result from the labelled superpixels.

1.3Notations and definitions

We consider that $\mathcal{G}=\{\mathcal{V},\mathcal{E},\mathsf{W}\}$ is an undirected weighted graph, where $\mathcal{V}$ is the set of $n$ nodes, $\mathcal{E}$ is the set of edges, and $\mathsf{W}\in\mathbb{R}^{n\times n}$ is the weighted adjacency matrix with nonnegative entries. We denote the graph Laplacian by $\mathsf{L}\in\mathbb{R}^{n\times n}$. We assume that $\mathsf{L}$ is real, symmetric, and positive semi-definite, e.g., the combinatorial graph Laplacian $\mathsf{L}:=\mathsf{D}-\mathsf{W}$, or the normalised one $\mathsf{L}:=\mathsf{I}-\mathsf{D}^{-1/2}\mathsf{W}\mathsf{D}^{-1/2}$. The matrix $\mathsf{D}\in\mathbb{R}^{n\times n}$ is the diagonal degree matrix and $\mathsf{I}\in\mathbb{R}^{n\times n}$ is the identity matrix [25]. The diagonal degree matrix $\mathsf{D}$ has entries $D_{ii}:=\sum_{i\neq j}\mathsf{W}_{ij}$.

We denote by $\mathsf{U}\in\mathbb{R}^{n\times n}$ the orthonormal eigenvectors of $\mathsf{L}$ and by $0=\lambda_{1}\leqslant\ldots\leqslant\lambda_{n}$ the ordered real eigenvalues of $\mathsf{L}$. We have $\mathsf{L}=\mathsf{U}\mathsf{\Lambda}\mathsf{U}^{\intercal}$, where $\mathsf{\Lambda}:={\rm diag}(\lambda_{1},\ldots,\lambda_{n})\in\mathbb{R}^{n\times n}$. The matrix $\mathsf{U}$ is the graph Fourier basis [8]. For any signal $\bm{x}\in\mathbb{R}^{n}$ defined on the nodes of the graph $\mathcal{G}$, $\hat{\bm{x}}=\mathsf{U}^{\intercal}\bm{x}$ contains the Fourier coefficients of $\bm{x}$ ordered in increasing frequencies. This work deals with $k$-bandlimited signals $\bm{x}\in\mathbb{R}^{n}$ on $\mathcal{G}$, i.e., signals whose Fourier coefficients $\hat{\bm{x}}_{k+1},\ldots,\hat{\bm{x}}_{n}$ are null. Let $\mathsf{U}_{k}$ be the restriction of $\mathsf{U}$ to its first $k$ vectors:

[TABLE]

Definition 1.1** ($k$-bandlimited signal on $\mathcal{G}$).**

A signal $\bm{x}\in\mathbb{R}^{n}$ defined on the nodes of the graph $\mathcal{G}$ is $k$-bandlimited with $k\in\mathbb{N}\setminus\{0\}$ if $\bm{x}\in{\rm span}(\mathsf{U}_{k})$, i.e., there exists $\bm{\eta}\in\mathbb{R}^{k}$ such that

[TABLE]

This definition was also used in [14, 11, 17]. We assume that $\lambda_{k}\neq\lambda_{k+1}$ to avoid any ambiguity in the definition of $k$-bandlimited signals.

Finally, for any matrix $\mathsf{X}\in\mathbb{R}^{n_{1}\times n_{2}}$, $\left\|\mathsf{X}\right\|_{2}$ denotes its spectral norm and $\left\|\mathsf{X}\right\|_{F}$ its Frobenius norm; when $n_{1}=n_{2}$, $\lambda_{\rm max}(\mathsf{X})$ denotes its largest eigenvalue and $\lambda_{\rm min}(\mathsf{X})$ its smallest eigenvalue. For any vector $\bm{x}\in\mathbb{R}^{n_{1}}$, $\left\|\bm{x}\right\|_{2}$ denotes the Euclidean norm of $\bm{x}$. Depending on the context, $\bm{x}_{j}$ may represent the $j^{\text{th}}$ entry of the vector $\bm{x}$ or the $j^{\text{th}}$ column-vector of the matrix $\mathsf{X}$. The entry on the $i^{\text{th}}$ row and $j^{\text{th}}$ column of $\mathsf{X}$ is denoted by $\mathsf{X}_{ij}$. The identity matrix is denoted by $\mathsf{I}$ (its dimensions are determined by the context).

We present in Fig. 2 a representation of the important variables and processes involved in this paper in order to facilitate the understanding of the different results.

2Sampling using groups of nodes

In this section, we explain our sampling strategies, starting with the definition of the groups of nodes.

2.1Grouping the nodes

We consider that the $n$ nodes of $\mathcal{G}$ are divided into $N$ different groups $\mathcal{N}_{1},\ldots,\mathcal{N}_{N}\subseteq\{1,\ldots,n\}$. The size of the $\ell^{\text{th}}$ group is denoted $\left|\mathcal{N}_{\ell}\right|$. We suppose that these groups form a partition of $\{1,\ldots,n\}$, so that each node belongs to exactly one group. We have

[TABLE]

For the object segmentation application discussed in the introduction, these groups represent the superpixels. However, we do not impose the groups to be made of neighbouring nodes in the graph; they can be made of nodes “far” from each other.

For each group $\mathcal{N}_{\ell}=\{n_{1}^{(\ell)},\ldots,n_{\left|\mathcal{N}_{\ell}\right|}^{(\ell)}\}$, we associate a matrix $\mathsf{N}^{(\ell)}~{}\in~{}\mathbb{R}^{\left|\mathcal{N}_{\ell}\right|\times n}$ that restricts a graph signal to the nodes appearing in $\mathcal{N}_{\ell}$, i.e.,

[TABLE]

Note that

[TABLE]

The case of overlapping groups can be handled by changing the definition of $\mathsf{N}^{(\ell)}$ to

[TABLE]

where $1\leqslant\beta_{i}\leqslant N$, $i=1,\ldots,n$, is the number of times node $i$ appears in the different groups $\mathcal{N}_{1},\ldots,\mathcal{N}_{N}$. Equation (7) also holds in this case. All results presented in Section 2 are valid for overlapping groups with this definition of $\mathsf{N}^{(\ell)}$.

2.2Sampling the groups

The sampling procedure consists in selecting $s$ groups out of the $N$ available ones. In the application of Section 1.2, it corresponds to the selection of the superpixels to label. We select these groups at random using a sampling distribution on $\{1,\ldots,N\}$ represented by a vector $\bm{p}\in\mathbb{R}^{N}$. The probability of selecting the $\ell^{\text{th}}$ group is $\bm{p}_{\ell}$. We assume that $\bm{p}_{\ell}>0$ for all $\ell=1,\ldots,N$, so that all groups may be selected with a non-zero probability. We obviously have $\sum_{\ell=1}^{N}\bm{p}_{\ell}=1$.

The indices $\Omega:=\{\omega_{1},\ldots,\omega_{s}\}$ of the selected groups are obtained by drawing independently – thus with replacements – $s$ indices from the set $\{1,\ldots,N\}$ according to $\bm{p}$, i.e.,

[TABLE]

The selected groups are $\mathcal{N}_{\omega_{1}},\ldots,\mathcal{N}_{\omega_{s}}$ and the total number of selected nodes is

[TABLE]

Once the groups are selected, we build the sampling matrix $\mathsf{M}\in\mathbb{R}^{m\times n}$ that satisfies

[TABLE]

and which restricts any signal to the nodes belonging to the selected groups. For a signal $\bm{x}\in\mathbb{R}^{n}$ defined on the nodes of $\mathcal{G}$, its sampled version $\bm{y}\in\mathbb{R}^{m}$ satisfies

[TABLE]

Our goal now is to determine what number $s$ is enough to ensure that all $k$-bandlimited signals can be reconstructed from its sampled version obtained with $\mathsf{M}$. To conduct this study, we need to define few more matrices. First, we associate the matrix

[TABLE]

to each group $\mathcal{N}_{\ell}$. Then, once the groups are drawn, we construct the block diagonal matrix $\mathsf{P}\in\mathbb{R}^{m\times m}$

[TABLE]

This matrix takes into account the probability of sampling each group and will be used to rescale $\bm{y}$ for norm preservation. This matrix ensures that $s^{-1}\,\mathbb{E}_{\Omega}\left\|\mathsf{P}\bm{y}\right\|_{2}^{2}=s^{-1}\,\mathbb{E}_{\Omega}\left\|\mathsf{P}\mathsf{M}\bm{x}\right\|_{2}^{2}=\left\|\bm{x}\right\|_{2}^{2}$.111This property is a consequence of (A), proved in Appendix A. Both matrices $\mathsf{P}$ and $\mathsf{M}$ depend on $\Omega$ and are random.

2.3Group graph coherence

The sampling procedure in [17] is similar to the ones proposed here at the difference that the nodes are sampled individually and not by groups. It was proved there that the number of nodes to sample is driven by a quantity called the graph coherence. This quantity measures how the energy of $k$-bandlimited signals spreads over the nodes. Similarly, we prove here that the number of groups to sample is driven by a quantity that measures the energy of $k$-bandlimited signals spreads over the groups. We now introduce this quantity.

The matrix $\mathsf{N}^{(\ell)}\mathsf{U}_{k}$ is the matrix that restricts a $k$-bandlimited signal to the nodes belonging to $\mathcal{N}_{\ell}$. Therefore,

[TABLE]

measures the energy on the nodes $\mathcal{N}_{\ell}$ of the normalised $k$-bandlimited signal that is most concentrated on $\mathcal{N}_{\ell}$. This energy varies between [math] and $1$. When this energy is close to $1$, there exists a $k$-bandlimited signal whose energy is essential concentrated on $\mathcal{N}_{\ell}$. This signal lives only on the nodes in $\mathcal{N}_{\ell}$ and does not spread elsewhere. On the contrary, when this energy is close to [math], there is no $k$-bandlimited signal living only on $\mathcal{N}_{\ell}$.

The sampling distribution $\bm{p}$ is adapted to the graph and the structure of the groups if: whenever $\|\mathsf{N}^{(\ell)}\mathsf{U}_{k}\|_{2}$ is high, $\bm{p}_{\ell}$ is high; whenever $\|\mathsf{N}^{(\ell)}\mathsf{U}_{k}\|_{2}$ is small, $\bm{p}_{\ell}$ is small. In other words, the ratio between $\|\mathsf{N}^{(\ell)}\mathsf{U}_{k}\|_{2}$ and $\bm{p}_{\ell}$ should be as constant as possible. This ensures that the groups where some $k$-bandlimited signals are concentrated are sampled with higher probability. Similarly to what was done in [17] with individual nodes, we define the group graph weighted coherence as the largest ratio between $\|\mathsf{N}^{(\ell)}\mathsf{U}_{k}\|_{2}$ and $\bm{p}_{\ell}^{-1/2}$.

Definition 2.1** (Group graph cumulative coherence).**

The group graph cumulative coherence of order $k$ is

[TABLE]

The quantity $\|\mathsf{N}^{(\ell)}\mathsf{U}_{k}\|_{2}$ is called the local group graph coherence.

In the extreme case where the groups $\mathcal{N}_{1},\ldots,\mathcal{N}_{N}$ reduce all to singletons, we recover the definition of the graph weighted coherence introduced in [17].

It is easy to prove that $\nu_{\bm{p}}$ is lower bounded by $1$. Indeed, for any $\bm{\eta}\in\mathbb{R}^{k}$ with $\left\|\bm{\eta}\right\|_{2}=1$, we have

[TABLE]

We have $\nu_{\bm{p}}=1$ in, for example, the degenerated case where $\mathcal{N}_{1}=\{1,\ldots,n\}$.

2.4Stable embedding

We now have all the tools to present our main theorem that shows that sampling $O(\nu_{\bm{p}}^{2}\log(k))$ groups is sufficient to stably embed the whole set of $k$-bandlimited signals. Hence, it is possible to reconstruct any $\bm{x}\in{\rm span}(\mathsf{U}_{k})$ from its measurements $\bm{y}=\mathsf{M}\bm{x}$.

Theorem 2.2** (Restricted isometry property - RIP).**

Let $\mathsf{M}$ be a random subsampling matrix constructed as in (16) using the groups $\mathcal{N}_{1},\ldots,\mathcal{N}_{N}$ and the sampling distribution $\bm{p}$. For any $\delta,\xi\in(0,1)$, with probability at least $1-\xi$,

[TABLE]

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

[TABLE]

Proof.

See Appendix A. ∎

In the above theorem, we recall that $s$ is the number of selected groups, each of them containing several nodes. We thus control the number of groups to sample and not directly the number of nodes. As the lower bound on $\nu_{\bm{p}}$ is $1$, sampling $O(\log(k))$ groups might already be sufficient if the groups and the sampling distribution are well-designed.

The number of groups to sample is driven by $\nu_{\bm{p}}$, which itself depends on the structure of the groups $\mathcal{N}_{1},\ldots,\mathcal{N}_{N}$ and on the sampling distribution $\bm{p}$. To reduce the number of samples to measure, we might optimise the structure of the groups and the sampling distribution. For example, if we were able to construct $n/L$ groups ($L\geqslant 1$) such that $\|\mathsf{N}^{(\ell)}\mathsf{U}_{k}\|_{2}\approx(L/n)^{1/2}$ – i.e., no $k$-bandlimited signals have more than $100\cdot(L/n)^{1/2}$ percent of its energy concentrated in each group – then setting $\bm{p}_{\ell}=L/n$, $\ell=1,\ldots,n/L$, would yield $\nu\approx 1$. In this case, sampling one group is enough to embed the set of $k$-bandlimited signals. However, it is not obvious how we can construct such groups in practice and we might not even have the flexibility to modify the structure of the groups. In such a case, the only possibility to reduce the number of measurements is to optimise the sampling distribution $\bm{p}$ to minimise $\nu_{\bm{p}}$.

The sampling distribution minimizing the coherence $\nu_{\bm{p}}$ is the distribution $\bm{p}^{*}\in\mathbb{R}^{N}$ that satisfies

[TABLE]

and for which

[TABLE]

Indeed, let $\bm{p}^{\prime}\neq\bm{p}^{*}$ be another sampling distribution. As the entries of $\bm{p}^{\prime}$ and $\bm{p}^{*}$ are nonnegative and sum to $1$, we necessarily have $\bm{p}^{\prime}_{\ell^{\prime}}<\bm{p}_{\ell^{\prime}}^{*}$ for some $\ell^{\prime}>0$. Then,

[TABLE]

where the last equality is obtained by replacing $\bm{p}_{\ell^{\prime}}^{*}$ with its value. Therefore, $\nu_{\bm{p}^{\prime}}>\nu_{\bm{p}^{*}}$ for any $\bm{p}^{\prime}\neq\bm{p}^{*}$. As similar proof can be found in, e.g., [18] where the authors derive the optimal sampling distribution for a compressive system.

We notice that

[TABLE]

Hence, by using this distribution, (24) shows that sampling $O(k\log(k))$ groups is always sufficient to sample all $k$-bandlimited signals. The exact number is proportional to $\nu_{\bm{p}^{*}}^{2}\log(k)$. This is not in contradiction with the fact that at least $k$ measurements are required in total as one group contains at least one node. We also have

[TABLE]

as $\left\|\mathsf{N}^{(\ell)}\mathsf{U}_{k}\right\|_{2}^{2}\leqslant 1$. Therefore, in any case, the bound never suggests to sample much more than $N$ groups, as one would expect it.

We recall that the results in [17] proves that it is always sufficient to sample $O(k\log(k))$ nodes to embed the set of $k$-bandlimited signals. When sampling the nodes by groups, it is the number of groups to sample that should be $O(k\log(k))$. This can be a large number of individual nodes but this is the potential price to pay when sampling the nodes by groups.

Variable density sampling [26, 27, 19] and structured sampling [18, 20, 21] are also important topics in compressed sensing. The method proposed here is closely inspired by these studies, especially by [20, 21]. Our results thus share several similarities with these works. We however benefit from a simpler signal model and take advantage of the graph structure to refine the results, propose simpler decoders to reconstruct the signal, and design efficient algorithms to estimate $\bm{p}^{*}$.

2.5A more practical result at small $k$’s

Optimising the sampling distribution reduces to estimating the spectral norm of the matrices $\mathsf{N}^{(\ell)}\mathsf{U}_{k}$. We will present in Section 4 a method avoiding the computation of $\mathsf{U}_{k}$. However, the method might still be too slow when a large number of groups is involved as it requires an estimation of a spectral norm for each group separately. It is thus interesting to characterise the performance of the proposed method using other quantities easier to compute.

In this section, we present results involving the following quantity

[TABLE]

The only difference between $\nu_{\bm{p}}$ and $\bar{\nu}_{\bm{p}}$ is that we substituted the Frobenius norm for the spectral norm. As $\left\|\mathsf{X}\right\|_{2}\leqslant\left\|\mathsf{X}\right\|_{F}$ for any matrix $\mathsf{X}$, we have $\nu_{\bm{p}}\leqslant\bar{\nu}_{\bm{p}}$; hence the results involving $\bar{\nu}_{\bm{p}}$ will be more pessimistic than those involving $\nu_{\bm{p}}$. We have $\bar{\nu}_{\bm{p}}\geqslant k$. Indeed,

[TABLE]

As with $\nu_{\bm{p}}$, the lower bound is also attained in, for example, the degenerated case where $\mathcal{N}_{1}=\{1,\ldots,n\}$. As $\nu_{\bm{p}}\leqslant\bar{\nu}_{\bm{p}}$, we have the following corollary to Theorem 2.2.

Corollary 2.3**.**

Let $\mathsf{M}$ be a random subsampling matrix constructed as in (16) using the groups $\mathcal{N}_{1},\ldots,\mathcal{N}_{N}$ and the sampling distribution $\bm{p}$. For any $\delta,\xi\in(0,1)$, with probability at least $1-\xi$, (23) holds for all $\bm{x}\in{\rm span}(\mathsf{U}_{k})$ provided that

[TABLE]

Proof.

As $\bar{\nu}_{\bm{p}}\geqslant\nu_{\bm{p}}$, (32) implies (24). Theorem 2.2 then proves that (23) holds with probability at least $1-\xi$ for all $\bm{x}\in{\rm span}(\mathsf{U}_{k})$. ∎

The sufficient condition (32) can be much more pessimistic than (24). Indeed, as we have $\bar{\nu}_{\bm{p}}^{2}\geqslant k$, Condition (32) suggests to always sample more than $O(k\log(k))$ groups, while we know that sampling $O(\log(k))$ can be enough. The interest of this result is thus in the regime where $k$ is small.

As we have done it with $\nu_{\bm{p}}$, we can also optimise $\bm{p}$ to minimise $\bar{\nu}_{\bm{p}}$. The sampling distribution that minimises $\bar{\nu}_{\bm{p}}$ is the distribution $\bm{q}^{*}$ satisfying

[TABLE]

and for which

[TABLE]

With this distribution, (32) proves that sampling $O(k\log(k))$ groups is always enough to sample all $k$-bandlimited signals. This result is particularly interesting at small $k$’s because estimating $\bm{q}^{*}$ is much easier and faster than estimating $\bm{p}^{*}$ (see Section 4).

In some particular cases, this result can be interesting also at large $k$’s. Indeed, (32) is a pessimistic bound. In reality, we may have $\nu_{\bm{q}^{*}}^{2}\leqslant k$ so that, according to (24), fewer samples than $O(k\log(k))$ are actually sufficient when using $\bm{q}^{*}$. Furthermore, $\bm{q}^{*}$ might actually be close to $\bm{p}^{*}$, in which case we would reach almost optimal results with $\bm{q}^{*}$. The occurrence of these events is however quite difficult to predict and depends on the structure of the groups and of the graph.

Note that we need the knowledge of the truncated Fourier matrix $\mathsf{U}_{k}$ to compute the distributions $\bm{p}^{*}$ and $\bm{q}^{*}$. Computing $\mathsf{U}_{k}$ can be intractable for large graphs. We will present a method to overcome this issue in Section 4. We continue first by explaining how to estimate $\bm{x}$ from its measurements.

3Fast reconstruction

Once the signal has been sampled, we shall also be able to reconstruct it. In [17], the authors propose to estimate the original signal by solving

[TABLE]

where $\gamma>0$ and $g$ is a nonnegative nondecreasing polynomial. We have

[TABLE]

where $\alpha_{0},\ldots,\alpha_{d}\in\mathbb{R}$ are the coefficients of the polynomial and $d\in\mathbb{N}$ its degree. These polynomial’s parameters can be tuned to improve the quality of the reconstruction. The role of $g$ can be viewed as a filter on $\mathcal{G}$ and should ideally be a high-pass filter. The matrix $\mathsf{P}$ is introduced to account for the RIP. The advantage of this method is that it can be solved efficiently even for large graph for example by conjugate gradient. Indeed, each step of this algorithm can be implemented by using only matrix-vector multiplications with $\mathsf{P}\mathsf{M}$ and $\mathsf{L}$, which are both sparse matrices. The matrix $g(\mathsf{L})$ does not have to be computed explicitly.

For brevity, we do not recall the theorem proving that (35) provides accurate and stable recovery of all $k$-bandlimited signals. However, this theorem applies as soon as the restricted isometry property (23) holds, and thus applies here when (24) holds. The reconstruction quality is controlled by $g(\lambda_{k})$ and the ratio ${g(\lambda_{k})}/{g(\lambda_{k+1})}$. One should seek to design a filter $g$ such that these quantities are as close as possible to zero to improve the reconstruction quality.

We propose now a method to obtain a faster estimation of the original signal when it is nearly piecewise constant.

3.1Piecewise constant graph signals

Before continuing, we want to stress that we consider non-overlapping groups $\mathcal{N}_{1},\ldots,\mathcal{N}_{N}$ in this rest of Section 3.

If a graph signal is nearly piecewise constant over the groups $\mathcal{N}_{1},\ldots,\mathcal{N}_{N}$ then reconstructing the mean values of this signal for each group is enough to obtain a good approximation of the original signal. Instead of estimating $n$ unknowns, we reduce the estimation to $N$ unknowns. When $N\ll n$, this is a large reduction of dimension yielding a significant speed up and memory reduction.

The fact that a signal $\bm{x}$ is piecewise constant over the groups $\mathcal{N}_{1},\ldots,\mathcal{N}_{N}$ is characterized as follows. We construct the averaging row-vectors $\bm{a}^{(\ell)}\in\mathbb{R}^{1\times n}$ that satisfy

[TABLE]

and the matrix

[TABLE]

As the groups do not overlap, we have

[TABLE]

hence $\left\|\mathsf{A}\right\|_{2}=1$. Applying $\mathsf{A}$ to $\bm{x}$ provides $N$ values, each one of them corresponding to the sum of the values of $\bm{x}$ within the group $\mathcal{N}_{\ell}$, scaled by $\left|\mathcal{N}_{\ell}\right|^{-1/2}$. Then, applying $\mathsf{A}^{\intercal}$ to $\mathsf{A}\bm{x}$ gives an approximation of $\bm{x}$ where the values in the vector $\mathsf{A}^{\intercal}\mathsf{A}\bm{x}$ are constant within each group; this is a piecewise constant vector over the groups. The value of $\mathsf{A}^{\intercal}\mathsf{A}\bm{x}$ appearing within the group $\mathcal{N}_{\ell}$ is exactly the average of $\bm{x}$ within $\mathcal{N}_{\ell}$. Saying that $\bm{x}$ is nearly constant within each group corresponds to assume that

[TABLE]

where $\epsilon\geqslant 0$ is a small value. The signal model of interest in this section is thus

[TABLE]

3.2Reducing the dimension

To build a fast algorithm exploiting the above property, we use a reconstruction method similar to (35) but involving vectors of smaller dimension. We define the averaged vector $\tilde{\bm{x}}:=\mathsf{A}\bm{x}\in\mathbb{R}^{N}$ of dimension $N$. As $\bm{x}\in\mathcal{A}_{\epsilon}$, estimating $\tilde{\bm{x}}$ is enough to get a good approximation of $\bm{x}$ – we just need to multiply it with $\mathsf{A}^{\intercal}$. Furthermore, as $\bm{x}$ is nearly piecewise constant over the groups $\mathcal{N}_{1},\ldots,\mathcal{N}_{N}$, by construction of the matrix $\mathsf{M}$, the measurement vector $\bm{y}=\mathsf{M}\bm{x}$ is also almost piecewise constant over the sampled groups $\mathcal{N}_{\omega_{1}},\ldots,\mathcal{N}_{\omega_{s}}$. We thus average $\bm{y}$ over these groups by multiplying it with the matrix $\widetilde{\mathsf{A}}\in\mathbb{R}^{s\times m}$ that satisfies

[TABLE]

We obtain

[TABLE]

We now have to link $\tilde{\bm{y}}$ to $\tilde{\bm{x}}$. We create the matrix $\widetilde{\mathsf{M}}\in\mathbb{R}^{s\times N}$ that restricts the $N$ mean value of $\tilde{\bm{x}}$ to the $s$ mean value of the selected groups, i.e.,

[TABLE]

We have therefore

[TABLE]

The goal is now to estimate $\tilde{\bm{x}}$ from $\tilde{\bm{y}}$. To ensure that the reconstruction method is stable to measurement noise, we do not consider the perfect scenario above but instead the scenario where

[TABLE]

and $\tilde{\bm{n}}\in\mathbb{R}^{s}$ models noise. We now need a regularisation term to estimate $\tilde{\bm{x}}$. We obtain this term by reducing the dimension of the regularisation involving the Laplacian $\mathsf{L}$ in (35). We compute

[TABLE]

Note that $\widetilde{\mathsf{L}}$ is a symmetric positive definite matrix. Like $g(\mathsf{L})$, it can thus be used as a regularisation. We thus propose to estimate $\tilde{\bm{x}}$ by solving

[TABLE]

where $\gamma>0$ and $\widetilde{\mathsf{P}}\in\mathbb{R}^{s\times s}$ is the diagonal matrix with entries satisfying

[TABLE]

Let $\tilde{\bm{x}}^{*}\in\mathbb{R}^{N}$ be a solution of (55). We finally obtain an estimation of $\bm{x}$ by computing $\mathsf{A}^{\intercal}\tilde{\bm{x}}^{*}$.

In the particular case where $g(\cdot)$ is the identity, one can notice that

[TABLE]

is non-zero only if there is at least one edge in $\mathcal{E}$ joining the groups $\mathcal{N}_{\ell}$ and $\mathcal{N}_{\ell^{\prime}}$. The Laplacian $\widetilde{\mathsf{L}}_{\ell\ell^{\prime}}$ preserves the connections present in the original graph represented by $\mathsf{L}$.

The dimension of the unknown vector in (55) is $N$, which can be much smaller than $n$. This leads to a large gain in memory and computation time when either $\widetilde{\mathsf{L}}$ or matrix-vector multiplications with $\widetilde{\mathsf{L}}$ can be computed rapidly. If $g$ has a small degree, then $\widetilde{\mathsf{L}}$ can be computed explicitly in a short amount of time. In such a case, we can solve (55) faster than (35) as it involves matrices of (much) smaller dimensions. In general, it is however not always straightforward to find an a efficient implementation for matrix-vector multiplications with $\widetilde{\mathsf{L}}$ without temporary going back to the signal domain of dimension $n$, i.e., multiplying the vector $\tilde{\bm{z}}$ with $\mathsf{A}^{\intercal}$, filtering the high dimensional signal $\mathsf{A}^{\intercal}\tilde{\bm{z}}$, and downsampling the result. Even though solving (55) might still be faster than solving (35) in this situation, we loose part of the efficiency by working temporarily in the high dimensional domain. We thus have less flexibility in the choice of $g$ with this reconstruction technique.

Let us also mention that (55) can be used to initialise the algorithm used to solve (35) with a good approximate solution. As in multigrid approaches to solve linear systems of equations, see, e.g., [28, 29].

The following theorem bounds the error between the signal recovered by (55) and the original vector.

Theorem 3.1**.**

Let $\Omega=\{\omega_{1},\ldots,\omega_{s}\}$ be a set of $s$ indices selected independently from $\{1,\ldots,N\}$ using a sampling distribution $\bm{p}\in\mathbb{R}^{N}$, $\mathsf{M},\mathsf{P},\widetilde{\mathsf{M}},\widetilde{\mathsf{P}}$ be the associated matrices constructed respectively in (16), (19), (51) and (56), and $M_{\rm max}>0$ be a constant such that $\left\|\mathsf{P}\mathsf{M}\right\|_{2}\leqslant M_{\rm max}$. Let $\xi,\delta\in(0,1)$ and suppose that $s$ satisfies (24). With probability at least $1-\xi$, the following holds for all $\bm{x}\in\mathcal{A}_{\epsilon}$, all $\tilde{\bm{n}}\in\mathbb{R}^{s}$, all $\gamma>0$, and all nonnegative nondecreasing polynomial functions $g$ such that $g(\lambda_{k+1})>0$.

Let $\tilde{\bm{x}}^{*}$ be the solution of (55) with $\tilde{\bm{y}}=\widetilde{\mathsf{M}}\mathsf{A}\bm{x}+\tilde{\bm{n}}$. Define $\bm{\alpha}^{*}:=\mathsf{U}_{k}\mathsf{U}_{k}^{\intercal}\,{\mathsf{A}}^{\intercal}\tilde{\bm{x}}^{*}$ and $\bm{\beta}^{*}:=(\mathsf{I}-\mathsf{U}_{k}\mathsf{U}_{k}^{\intercal})\,{\mathsf{A}}^{\intercal}\tilde{\bm{x}}^{*}$. Then,

[TABLE]

and

[TABLE]

Proof.

See Appendix B. ∎

The vector $\bm{\alpha}^{*}$ is the orthogonal projection of ${\mathsf{A}}^{\intercal}\tilde{\bm{x}}^{*}$ onto ${\rm span}(\mathsf{U}_{k})$. The vector $\bm{\beta}^{*}$ is the projection of ${\mathsf{A}}^{\intercal}\tilde{\bm{x}}^{*}$ onto the orthogonal complement of ${\rm span}(\mathsf{U}_{k})$. There are several remarks to make about the above theorem:

• •

Theorem 3.1 shows that the result obtained via (55) is similar to the one we would have obtained by solving (35) – see [17] for the error bounds – with additional errors controlled by $\epsilon$. We recall that $\epsilon$ characterises how far $\bm{x}$ is from a piecewise constant signal. As expected, the smaller $\epsilon$, the better the reconstruction.

• •

The reconstruction quality improves when $g(\lambda_{k})$ and the ratio ${g(\lambda_{k})}/{g(\lambda_{k+1})}$ go to [math], and when ${g(\lambda_{n})}/{g(\lambda_{k+1})}$ tends to $1$. We recall that we have $g(\lambda_{n})\geqslant g(\lambda_{k+1})>g(\lambda_{k})$ by assumption.

• •

The effect of the noise $\tilde{\bm{n}}$ decreases when $g(\lambda_{k+1})$ increases, and, obviously, $\gamma$ should be adapted to the signal-to-noise ratio.

Let us mention that the idea of “coarsening” a graph and the signals that live on it using a partition into different groups of nodes can also be found in [30], where a multiresolution analysis method of graph signals is proposed. The coarsening method is however different than the one used here.

4Optimal sampling estimation

In this section, we come back to the sampling process of Section 2 and leave the reconstruction problem. We explain how to estimate the sampling distributions $\bm{p}^{*}$ and $\bm{q}^{*}$ of Section 2 without computing the truncated Fourier matrix $\mathsf{U}_{k}$, as this computation is intractable for large graphs. The methods below only involve matrix-vector multiplications with the sparse Laplacian matrix $\mathsf{L}$ and are thus computationally tractable even for large $n$.

4.1Estimation of $\bm{p}^{*}$

The distribution $\bm{p}^{*}$, defined in (25), that minimises the coherence $\nu_{\bm{p}}$ is entirely defined by the values of $\|\mathsf{N}^{(\ell)}\mathsf{U}_{k}\|_{2}$ for $\ell=1,\ldots,N$, which are thus the quantities we need to evaluate.

We recall that, in graph signal processing, a filter is represented by a function $h:\mathbb{R}\rightarrow\mathbb{R}$, and that the signal $\bm{x}$ filtered by $h$ is

[TABLE]

where $\hat{\bm{h}}=(h(\lambda_{1}),\ldots,h(\lambda_{n}))^{\intercal}\in\mathbb{R}^{n}$. To filter the signal $\bm{x}$ without actually computing the graph Fourier transform of $\bm{x}$, we can approximate the function $h$ by a polynomial

[TABLE]

of degree $d$, and compute $\bm{x}_{r}$ instead of $\bm{x}_{h}$. The filtered signal $\bm{x}_{r}$ is computed rapidly using the formula

[TABLE]

that involves only matrix-vector multiplications with the sparse Laplacian matrix $\mathsf{L}$. We let the reader refer to [31] for more information on this fast filtering technique. For any polynomial function $r$ of the form above and any matrix $\mathsf{A}\in\mathbb{R}^{n\times n}$, we define

[TABLE]

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

Let $i_{\lambda_{k}}:\mathbb{R}\rightarrow\mathbb{R}$ be the ideal low-pass filter at cutoff frequency $\lambda_{k}$, i.e., the filter that satisfies

[TABLE]

We have $\mathsf{U}_{k}\mathsf{U}_{k}^{\intercal}=i_{\lambda_{k}}\left(\mathsf{L}\right)$. Then, we notice that

[TABLE]

We recall that $\mathsf{N}^{(\ell)}$ is the matrix that restricts the signal to the nodes belonging to $\mathcal{N}_{\ell}$. The matrix appearing on the right hand side of the last equality corresponds to the linear operator that 1) extends a vector on the complete graph by inserting [math] in all groups $\ell^{\prime}\neq\ell$, 2) low-pass filters the extended signal, 3) restricts the result to the group $\mathcal{N}_{\ell}$. This process can be approximated by replacing the ideal low-pass filter $i_{\lambda_{k}}$ with a polynomial approximation $\tilde{i}_{\lambda_{k}}$ of $i_{\lambda_{k}}$ and

[TABLE]

To estimate $\bm{p}^{*}$, we estimate the spectral norm of the matrix appearing on the right hand side for which matrix-vector multiplication is fast. The quality of the approximation depends obviously on the choice of the polynomial $\tilde{i}_{\lambda_{k}}$.

Estimating $\|\mathsf{N}^{(\ell)}\;\tilde{i}_{\lambda_{k}}\left(\mathsf{L}\right)\;{\mathsf{N}^{(\ell)}}^{\intercal}\|_{2}$ amounts to computing the largest eigenvalue of $\mathsf{N}^{(\ell)}\;\tilde{i}_{\lambda_{k}}\left(\mathsf{L}\right)\;{\mathsf{N}^{(\ell)}}^{\intercal}$ which can be done, e.g., by using the power method. This method requires matrix-vector multiplication only with $\mathsf{N}^{(\ell)}$ and $\tilde{i}_{\lambda_{k}}\left(\mathsf{L}\right)$ and is thus fast. Finally, the approximation $\bar{\bm{p}}\in\mathbb{R}^{N}$ of $\bm{p}^{*}$ satisfies

[TABLE]

Note that an estimation of $\lambda_{k}$ is required beforehand to define the filter $\tilde{i}_{\lambda_{k}}$. We estimate this value using the dichotomy method presented in [17].

4.2Estimation of $\bm{q}^{*}$

Computing $\bar{\bm{p}}$ requires the estimation of $N$ eigenvalues. Even though these estimations can be done in parallel, this process might still be too slow for certain applications. As explained before, when $k$ is small, we can use the sampling distribution $\bm{q}^{*}$ in (33) that minimises $\bar{\nu}_{\bm{p}}$. This distribution is faster to compute than $\bar{\bm{p}}$.

We start by noticing that we have

[TABLE]

for each group $\ell=1,\ldots,N$. The vector $\bm{\delta}_{i}\in\mathbb{R}^{n}$ is the unit vector that is null on all nodes except at node $i$. Hence, we only need an estimation of $\|\mathsf{U}_{k}^{\intercal}\bm{\delta}_{i}\|_{2}^{2}$, $i=1,\ldots,n$, to estimate $\bm{q}^{*}$. An algorithm was already proposed in [17] to estimate these values. We let the reader refer to Algorithm $1$ in [17] for the details of the method. We just recall that this estimation is obtained by filtering $O(\log(n))$ random signals with a polynomial approximation of $i_{\lambda_{k}}$. Finally, our estimation $\bar{\bm{q}}\in\mathbb{R}^{N}$ of $\bm{q}^{*}$ has entries

[TABLE]

where each $\|\mathsf{U}_{k}^{\intercal}\bm{\delta}_{i}\|_{2}^{2}$ are estimated by Algorithm $1$ in [17].

This estimation is faster than for $\bar{\bm{\bm{p}}}$ because the power method is an iterative method that involves one filtering at each iteration. Furthermore, the power method is run independently for each group $\mathcal{N}_{\ell}$. In total, (much) more than $\mathcal{N}_{\ell}$ filterings are thus required. On the contrary, for $\bar{\bm{q}}$, we just need to filter $O(\log(n))$ signals to obtain the estimation. In most situations, we already have $O(\log(n))\leqslant\mathcal{N}_{\ell}$ and computing $\bar{\bm{q}}$ is thus faster than computing $\bar{\bm{p}}$.

5Experiments

In this last section, we first test our sampling strategies on two different graphs to illustrate the effect of the different sampling distributions on the minimum number of samples required to ensure that the RIP holds. Then, we apply our sampling strategy for user-guided object segmentation. In this application, we also test the different recovery techniques proposed in Section 3.

5.1Sampling distributions

We perform experiments on two different graphs: the Minnesota graph of size $n=2642$ and the bunny graph of size $n=2503$. Both graphs are presented in Fig. 3 and are available in the GSP toolbox [32]. For each graph, we group the nodes using the spatial coordinates associated to each node. For the Minnesota graph, we divide the space into $100$ cells and group the nodes that fall in the same cell. After removing empty cells, we obtain the $N=73$ groups represented in Fig. 3. For the bunny graph, we obtain $N=213$ groups with a similar procedure (see Fig. 3).

For each graph, we compute the combinatorial Laplacian and $\mathsf{U}_{k}$ for different values of $k$. Then, we compute the lower RIP constant, i.e., the constant $\underline{\delta}_{k}>0$ that satisfies

[TABLE]

This constant is the smallest value that $\delta$ can take such that the left-hand side of the RIP (23) holds. Remark that

[TABLE]

We estimate $\underline{\delta}_{k}$ for $500$ independent draws of the set $\Omega$, which defines the matrices $\mathsf{P}\mathsf{M}$, and different numbers of selected groups $s$. All samplings are done in the conditions of Theorem 2.2 using the sampling distributions $\bm{u},\bm{p}^{*},\bar{\bm{p}}$ and $\bar{\bm{q}}$. The vector $\bm{u}$ denotes the uniform distribution over $\{1,\ldots,N\}$. When conducting this experiment with the estimated distributions $\bar{\bm{p}}$ and $\bar{\bm{q}}$, we re-estimate these distributions at each of the $500$ trials. These distributions are estimated using Jackson-Chebychev polynomials of order $50$ [33]. For the Minnesota graph, we consider the bandlimits $k=5,10,20$. For the bunny graph, we consider the bandlimits $k=10,25,50$.

We present the probability that $\underline{\delta}_{k}$ is less than $0.995$, estimated over the $500$ draws of $\Omega$, as a function of $s$ in Fig. 3. For the Minnesota graph, the performance is better when using the optimal distribution $\bm{p}^{*}$ than when using the uniform distribution $\bm{u}$ for all $k$, which is in line with the theory. The estimated $\bar{\bm{p}}$ and $\bar{\bm{q}}$ yield performance equivalent to $\bm{p}^{*}$. This confirms that we can achieve similar sampling performance without having to compute the Fourier matrix $\mathsf{U}_{k}$, which, we recall, is intractable for large graphs. This also shows that $\bar{\bm{q}}$ can lead to nearly optimal results. For the bunny graph, all sampling distributions yield essentially the same results at all bandlimits. We notice a slight improvement at $k=50$ when using $\bar{\bm{p}}$, $\bar{\bm{q}}$ or $\bm{p}^{*}$ instead of $\bm{u}$.

For illustration, we present in Fig. 4 examples of computed sampling distributions $\bm{p}^{*}$, $\bar{\bm{p}}$ and $\bar{\bm{q}}$. All sampling distributions exhibit similar structures, which explains why they all yield about the same performance in our experiments.

5.2Object segmentation

5.2.1Protocol

We now test our method for interactive object segmentation. We consider the image of size $n=321\times 481$ presented in Fig. 1 for which our goal is to segment the tiger. The ground truth segmentation map $\bm{x}\in\{0,1\}^{n}$ is presented in Fig. 5. The value $1$ (white) indicates the presence of the tiger and the value [math] (black) stands for the foreground. The original image and the ground truth image are part of the dataset available222http://www.ntu.edu.sg/home/asjfcai/Benchmark_Website/benchmark_index.html. in [34]. Our objective is to recover the original map $\bm{x}$ from few user-inputs. To facilitate the interactions with the user, we divide the original image into the $N=600$ superpixels showed in Fig. 1 and computed with SLIC [23], choose a small number of superpixels at random and ask the user to label these superpixels: $1$ if the superpixel belongs to the tiger; [math] otherwise.

The graph $\mathcal{G}$ used to propagate the user-labels to the complete image is constructed as follows. We build a feature vector for each pixel by extracting a color RGB patch of size $3\times 3$ around the pixel, transform this patch in vector form, and augment this vector with the absolute 2D coordinates of the pixels in this extracted patch. This yields $n$ feature vectors $\bm{g}_{i}\in\mathbb{R}^{45}$, $i=1,\ldots,n$. We then connect each feature vector to its $9$ nearest neighbours (in the Euclidean sense), which gives a set of $9n$ edges $\mathcal{E}$. The adjacency matrix $\mathsf{W}\in\mathbb{R}^{n\times n}$ satisfies $\mathsf{W}_{ij}:=\exp[-{\|\bm{g}_{i}-\bm{g}_{j}\|_{2}^{2}}/{\sigma^{2}}],$ where $\sigma>0$ is the $25^{\text{th}}$ percentile of the set $\{\|\bm{g}_{i}-\bm{g}_{j}\|_{2}:(i,j)\in\mathcal{E}\}$. We finally symmetrise the matrix $\mathsf{W}$ and compute the combinatorial Laplacian $\mathsf{L}\in\mathbb{R}^{n\times n}$.

We study three strategies to choose the superpixels. The first strategy consists in choosing the superpixels uniformly at random, i.e., using the sampling distribution $\bm{u}$. The second and third strategies consist in choosing the superpixels with respectively the optimised distributions $\bar{\bm{q}}$ and $\bar{\bm{p}}$, which we evaluate at $k_{0}=50$ using Jackson-Chebychev polynomials of order $75$ [33]. For illustration, we present in Fig. 5 the estimated values $\|\mathsf{N}^{(\ell)}\mathsf{U}_{k_{0}}\|_{F}^{2}$ and $\|\mathsf{N}^{(\ell)}\mathsf{U}_{k_{0}}\|_{2}^{2}$, which define the optimised distributions $\bar{\bm{q}}$ and $\bar{\bm{p}}$, respectively. Both distributions indicate that one should label more superpixels around the tiger. The distribution $\bar{\bm{q}}$ “focuses” however more on specific regions, like the head of tiger. The distribution $\bar{\bm{p}}$ spreads the measurements over the entire tiger more uniformly.

We emulate user-interactions as follows. For each chosen superpixel, we compute the mean of the ground truth map $\bm{x}$ within this superpixel. If the mean value is larger than $0.5$, we label the superpixel as part of the tiger. Otherwise, we label the superpixel as part of the background. This strategy obviously introduces noise if some superpixels cover part of the background and of the tiger. Once the labelling is done, we have access to the measurement vector $\tilde{\bm{y}}\in\mathbb{R}^{m}$ from which we want to reconstruct $\bm{x}$. We repeat this procedure for $s\in\{50,70,\ldots,250\}$. For each $s$, we also repeat the experiments $50$ times with independent draws of the superpixels. We draw the superpixels with replacements in all cases.

To reconstruct the original map $\bm{x}$, we first use the fast reconstruction method (55) and then refine the solution at the pixel level with (35), using the solution of the first minimisation problem as initialisation to solve the second minimisation problem. We choose $g(\mathsf{L})=\mathsf{L}$ and solve both problems in the limit where $\gamma\rightarrow 0$. In this limit, the problems (55) and (35) become

[TABLE]

and

[TABLE]

respectively. Both problems are solved using FISTA [35]. The same stopping criteria are used for all experiments.

5.2.2Results

We present in the top panel of Fig. 6 the reconstruction snr obtained with the different methods. The reconstruction snr is defined as $-20\log_{10}(-\left\|\bm{x}-\bm{x}^{*}\right\|_{2}/\left\|\bm{x}\right\|_{2})$, where $\bm{x}^{*}$ is the reconstructed signal. We notice that the snr attained with the fast decoder (55) is very similar to the snr attained with (35). We also remark that the optimised distributions $\bar{\bm{p}}$ and $\bar{\bm{q}}$ yield better reconstructions than the uniform distribution $\bm{u}$. The mean reconstruction snr is slightly better with $\bar{\bm{p}}$ than with $\bar{\bm{q}}$ at $s\geqslant 150$.

We present the computation time of each method in the bottom panel of Fig. 6. We notice that solving (55) is much faster than solving (35), while they yield almost the same quality. This highlight the interest of the fast reconstruction technique. It is also interesting to note that it is faster to solve (35) when the measurements are drawn with $\bar{\bm{p}}$ or with $\bar{\bm{q}}$ than with $\bm{u}$. The reason is probably a better initialisation of (35) or a better “quality” of the measurements with the optimised distributions than with the uniform distribution.

Finally, we present in Fig. 5 some examples of reconstructions from $s=150$ sampled superpixels for each method. We notice that the optimised sampling distributions improve the reconstruction of $\bm{x}$ around the head and tail of the tiger, i.e., where the optimised distributions have higher values. With a uniform distribution, the structure of the graph makes it difficult to reconstruct $\bm{x}$ around the head and tail from the values of other superpixels. The optimised sampling distribution compensate this issue by favouring this area when selecting the measurements.

6Discussion and Conclusion

We presented a structured sampling strategy for $k$-bandlimited signals where the nodes are selected by groups. We proved that the local group graph cumulative coherence quantifies the importance of sampling each group to ensure a stable embedding of all $k$-bandlimited signals. Finally, we presented a fast reconstruction technique for $k$-bandlimited signals which are also nearly piecewise-constant over pre-defined groups of nodes.

Among the possible applications of these methods, we believe that they can also be useful to accelerate the compressive spectral clustering method proposed in [36]. After having computed some features vectors for each nodes, this compressive method works by downsampling the set of features of vectors, performing $k$-means on this reduced set to find $k$ clusters, and interpolating the clustering results on all nodes by solving (35). To accelerate the method, one could 1) pre-group similar nodes to form $N$ groups such that $k\leqslant N\ll n$, e.g., by running few iterations of the $k$-means algorithm; 2) subsample this set of $N$ groups; 3) cluster this subset to find $k$ clusters; and 4) solve (55) to cluster all nodes. If the overhead of computing the $N$ groups is small, this method has the potential to be faster than the original compressive spectral clustering method.

Finally, we would like to discuss two limitations in the proposed methods. First, the optimal sampling distribution depends on the parameter $k$. In some applications, the final result may change a lot depending on the value of $k$ which was chosen to compute this distribution. Finding a range of values of $k$ which give acceptable and stable results is thus an important step in every application. Second, the estimation of the optimal sampling distribution depends on the quality of the polynomial approximation of the ideal low-pass filter $\lambda_{\lambda_{k}}$. It is sometimes necessary to use a polynomial of large degree to get a correct estimation, which limits the computational efficiency of the proposed methods. In such cases, it would be especially useful to find more efficient alternatives to estimate the distributions $\bm{p}^{*}$ and $\bm{q}^{*}$.

Appendix A - Proof of the Theorem 2.2

As done in [17], the proof is obtained by applying the following lemma obtained by Tropp in [37].

Lemma A.1** (Theorem $1.1$, [37]).**

Consider a finite sequence $\{\mathsf{X}_{j}\}$ of independent, random, self-adjoint, positive semi-definite matrices of dimension $d\times d$. Assume that each random matrix satisfies $\lambda_{\rm max}(\mathsf{X}_{j})\leqslant R$ almost surely. Define

[TABLE]

Then

[TABLE]

We will also use the facts that, for all $\delta\in[0,1]$,

[TABLE]

Proof of Theorem 2.2.

We start by noticing that

[TABLE]

We define

[TABLE]

The matrix $\mathsf{X}$ is thus a sum of $s$ independent, random, self-adjoint, positive semi-definite matrices. We are in the setting of Lemma A.1. We continue by computing $\mathbb{E}\,\mathsf{X}_{j}$ and $\lambda_{\rm max}(\mathsf{X}_{j})$.

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

[TABLE]

Therefore, $\lambda_{\rm min}\left(\sum_{j}\mathbb{E}\,\mathsf{X}_{j}\right)=1$ and $\lambda_{\rm max}\left(\sum_{j}\mathbb{E}\,\mathsf{X}_{j}\right)=1.$ Furthermore, for all $j=1,\ldots,s$, we have

[TABLE]

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

[TABLE]

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

[TABLE]

provided that (24) holds. Noticing that (85) implies that

[TABLE]

for all $\bm{\alpha}\in\mathbb{R}^{k}$, which is equivalent to (23) for all $\bm{x}\in{\rm span}(\mathsf{U}_{k})$, terminates the proof. ∎

Appendix B - Proof of Theorem 3.1

In order to prove Theorem 3.1, we need to establish few properties between the different matrices used in this work.

The first useful property is

[TABLE]

Indeed, for any $\bm{z}\in\mathbb{R}^{n}$, the $j^{\text{th}}$ entry of $\widetilde{\mathsf{P}}\widetilde{\mathsf{M}}\mathsf{A}\bm{z}$ is

[TABLE]

Then, the $j^{\text{th}}$ entry of $\widetilde{\mathsf{A}}\mathsf{P}\mathsf{M}\bm{z}$ is the scaled sum of the values in the $j^{\text{th}}$ sampled group appearing in $\mathsf{P}\mathsf{M}\bm{z}$, which is $\bm{p}_{\omega_{j}}^{-1/2}\mathsf{N}^{(\omega_{j})}\bm{z}$. From the definition of $\widetilde{\mathsf{A}}$, the sum is scaled by $\left|\mathcal{N}_{\omega_{j}}\right|^{-1/2}$. Therefore, $(\widetilde{\mathsf{P}}\widetilde{\mathsf{M}}\mathsf{A}\bm{z})_{j}=(\widetilde{\mathsf{A}}\mathsf{P}\mathsf{M}\bm{z})_{j}$ for all $j\in\{1,\ldots,s\}$, which terminates the proof.

The second property is

[TABLE]

which implies $\|\widetilde{\mathsf{A}}\|_{2}=1$.

The third property is

[TABLE]

for all $\tilde{\bm{z}}\in\mathbb{R}^{N}$. To prove this property, we remark that $\mathsf{P}\mathsf{M}\mathsf{A}^{\intercal}\tilde{\bm{z}}=\widetilde{\mathsf{A}}^{\intercal}\widetilde{\mathsf{P}}\widetilde{\mathsf{M}}\tilde{\bm{z}}.$ Indeed, the entries in $\mathsf{P}\mathsf{M}\mathsf{A}^{\intercal}\tilde{\bm{z}}$ corresponding to the first selected group $\mathcal{N}_{\omega_{1}}$ are all equal to $(\bm{p}_{\omega_{1}}\,\left|\mathcal{N}_{\omega_{1}}\right|)^{-1/2}\tilde{\bm{z}}_{\omega_{1}}$. There are $\left|\mathcal{N}_{\omega_{1}}\right|$ such entries. One can also notice that the first $\left|\mathcal{N}_{\omega_{1}}\right|$ entries of $\widetilde{\mathsf{A}}^{\intercal}\widetilde{\mathsf{P}}\widetilde{\mathsf{M}}\tilde{\bm{z}}$ are also equal to $(\bm{p}_{\omega_{j}}\,\left|\mathcal{N}_{\omega_{1}}\right|)^{-1/2}\tilde{\bm{z}}_{\omega_{1}}$. Repeating this reasoning for all the sampled groups proves the equality. On the one side, we thus have $\left\|\mathsf{P}\mathsf{M}\mathsf{A}^{\intercal}\tilde{\bm{z}}\right\|_{2}=\|\widetilde{\mathsf{A}}^{\intercal}\widetilde{\mathsf{P}}\widetilde{\mathsf{M}}\tilde{\bm{z}}\|_{2}=\|\widetilde{\mathsf{P}}\widetilde{\mathsf{M}}\tilde{\bm{z}}\|_{2},$ where we used (89). On the other side, we have $\|\widetilde{\mathsf{A}}\mathsf{P}\mathsf{M}\mathsf{A}^{\intercal}\tilde{\bm{z}}\|_{2}=\|\widetilde{\mathsf{P}}\widetilde{\mathsf{M}}\mathsf{A}\mathsf{A}^{\intercal}\tilde{\bm{z}}\|_{2}=\|\widetilde{\mathsf{P}}\widetilde{\mathsf{M}}\tilde{\bm{z}}\|_{2},$ where we used (87) and (42). This terminates the proof.

Proof of Theorem 3.1.

As $\bm{x}^{*}$ is a minimiser of (55), we have

[TABLE]

To prove the theorem, we need to lower and upper bound the left and right hand sides of (91), respectively. We start with the bound involving $\widetilde{\mathsf{L}}$ and then with the ones involving $\widetilde{\mathsf{M}}$.

[Bounding the terms in (91) involving $\widetilde{\mathsf{L}}$]. We define the matrices

[TABLE]

By definition of $\bm{\alpha}^{*}$ and $\bm{\beta}^{*}$, $\mathsf{A}^{\intercal}\tilde{\bm{x}}^{*}=\bm{\alpha}^{*}+\bm{\beta}^{*}$ with $\bm{\alpha}^{*}\in{\rm span}(\mathsf{U}_{k})$ and $\bm{\beta}^{*}\in{\rm span}(\bar{\mathsf{U}}_{k})$. We recall that $\widetilde{\mathsf{L}}=(\mathsf{A}\mathsf{U})\;g(\mathsf{\Lambda})\;(\mathsf{A}\mathsf{U})^{\intercal}$. Therefore, we obtain

[TABLE]

In the first step, we used the facts that $\bar{\mathsf{U}}_{k}^{\intercal}\bm{\alpha}^{*}=\bm{0}$ and $\mathsf{U}_{k}^{\intercal}\bm{\beta}^{*}=\bm{0}$. The second step follows form the fact that $\left\|\bar{\mathsf{U}}_{k}^{\intercal}\bm{\beta}^{*}\right\|_{2}=\left\|\bm{\beta}^{*}\right\|_{2}$. We also have

[TABLE]

The second inequality follows from the facts that $\left\|\mathsf{U}_{k}\right\|_{2}=1$ and $\tilde{\bm{x}}=\mathsf{A}\bm{x}$. To obtain the third inequality, we used $\left\|\mathsf{A}\right\|_{2}=1$ and the triangle inequality. For the last step, notice that $\left\|\bar{\mathsf{U}}_{k}^{\intercal}\bm{x}\right\|_{2}=0$ (as $\bm{x}\in{\rm span}(\mathsf{U}_{k})$), $\left\|\bar{\mathsf{U}}_{k}\right\|_{2}=1$ and use (43).

[Bounding the terms in (91) involving $\widetilde{\mathsf{M}}$]. By definition of $\tilde{\bm{y}}$, it is immediate that

[TABLE]

For the other term involving $\widetilde{\mathsf{M}}$, the triangle inequality yields

[TABLE]

Then, we have

[TABLE]

The first equality follows from $\mathsf{A}\mathsf{A}^{\intercal}=\mathsf{I}$, the second from (87), and the triangle inequality was used in the last step. To summarise, we are at

[TABLE]

We continue by lower bounding $\|\tilde{\mathsf{A}}\mathsf{P}\mathsf{M}\mathsf{A}^{\intercal}\tilde{\bm{x}}^{*}-\tilde{\mathsf{A}}\mathsf{P}\mathsf{M}\mathsf{A}^{\intercal}\mathsf{A}\bm{x}\|_{2}$ and upper bounding $\|\tilde{\mathsf{A}}\mathsf{P}\mathsf{M}\mathsf{A}^{\intercal}\mathsf{A}\bm{x}-\tilde{\mathsf{A}}\mathsf{P}\mathsf{M}\bm{x}\|_{2}$ separately.

[Lower bound on $\|\tilde{\mathsf{A}}\mathsf{P}\mathsf{M}\mathsf{A}^{\intercal}\tilde{\bm{x}}^{*}-\tilde{\mathsf{A}}\mathsf{P}\mathsf{M}\mathsf{A}^{\intercal}\mathsf{A}\bm{x}\|_{2}$]. Equality (90) yields

[TABLE]

Using the triangle inequality and the fact that $\bm{x}\in\mathcal{A}_{\epsilon}$, we obtain

[TABLE]

The restricted isometry property (23) then yields

[TABLE]

We used the equality $\mathsf{A}^{\intercal}\tilde{\bm{x}}^{*}=\bm{\alpha}^{*}+\bm{\beta}^{*}$. We thus proved that

[TABLE]

[Upper bounding $\|\tilde{\mathsf{A}}\mathsf{P}\mathsf{M}\mathsf{A}^{\intercal}\mathsf{A}\bm{x}-\tilde{\mathsf{A}}\mathsf{P}\mathsf{M}\bm{x}\|_{2}$]. We have

[TABLE]

We used the facts that $\|\tilde{\mathsf{A}}\|_{2}=1$ (see (89)) and that $\bm{x}\in\mathcal{A}_{\epsilon}$ in the last inequality.

Using the inequalities (104) and (105) in (100), we arrive at

[TABLE]

[Finishing the proof]. Remark that (91) implies

[TABLE]

Using (95), (B), (97) in (108), we obtain

[TABLE]

which implies (59) in Theorem 3.1. It remains to prove (3.1) to finish the proof.

Using (B), (97) and (106) in (107), we obtain

[TABLE]

Using (59) to bound $\left\|\bm{\beta}^{*}\right\|_{2}$ on the right hand side, we have

[TABLE]

We only rearranged the term in the last step. This proves (3.1) and terminates the proof. ∎