Structured sampling and fast reconstruction of smooth graph signals
Gilles Puy Patrick Pérez


Analysis of the coupled Navier-Stokes/Biot problem
    Aycil Cesmelioglu
    December 30, 2023
==================================================





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 -bandlimited signals is driven by a quantity called the group graph cumulative coherence. For some optimised sampling distributions, we show that sampling  O(log()) groups is always sufficient to stably embed the set of -bandlimited signals but that this number can be smaller – down to O(log()) – depending on the structure of the groups of nodes. Fast methods to approximate these sampling distributions are detailed. Second, we consider -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.

Graph signal processing, compressive sampling, bandlimited graph signals.










§ INTRODUCTION


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, , propagate foreground/background scribbles for object segmentation <cit.> or propagate manual color/tone modifications  <cit.>.
First, a graph  modelling the similarities between the pixels of the image is built. Second, the user-edits specify values at some nodes (pixels) of . Finally, the complete signal is estimated by assuming that it is smooth on . The quality of the propagation depends on the structure of  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, , in <cit.>. We address these problems from a graph signal processing point-of-view <cit.>. 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 . In the context of graph signal processing, smooth or -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 <cit.>. Two different -bandlimited signals are also necessarily different when restricted to a uniqueness set. Therefore, one can sample all -bandlimited signals on a uniqueness set. Then, Anis et al., <cit.> and Chen et al., <cit.> proved that one can always find a uniquess set of  nodes to sample all -bandlimited signals. Finding this set is however computationally expensive. In <cit.>, 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 <cit.>. Recently, Puy et al. <cit.> proved that there always exists a random sampling strategy for which sampling O(log()) 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 -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 -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 -bandlimited signals, , -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.   





 §.§ Contributions


The random sampling strategy that we propose generalises the method proposed in <cit.>. We use here a structured sampling strategy. Let us already acknowledge here that such strategies are also studied in the field of compressed sensing <cit.> 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 <cit.> and characterises how much the energy of -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(log()) groups of nodes is sufficient to ensure the reconstruction of all -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()) – 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(log()) groups is always sufficient to ensure the reconstruction of all -bandlimited signals. This strategy is mainly interesting at small .

Finally, we propose a fast reconstruction method for -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.





 §.§ Applications



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 <cit.> for instance, we build a graph  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 . To propose the regions to label, we view this segmentation problem as a sampling problem of a smooth signal on  where the regions to label are chosen to ensure the recovery of a -bandlimited signal. To choose the regions to label, we start by partitioning the image into superpixels, , with SLIC technique <cit.>. We see in Fig. <ref> 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 <cit.>. 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 , 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.










 §.§ Notations and definitions
 

We consider that = {𝒱, ℰ, W} is an undirected weighted graph, where 𝒱 is the set of  nodes, ℰ is the set of edges, and W∈^× is the weighted adjacency matrix with nonnegative entries. We denote the graph Laplacian by ∈^×. We assume that  is real, symmetric, and positive semi-definite, , the combinatorial graph Laplacian := D - W, or the normalised one := I - D^-1/2WD^-1/2. The matrix D∈^× is the diagonal degree matrix and I∈^× is the identity matrix <cit.>. The diagonal degree matrix D has entries D_ii := ∑_i ≠ jW_ij.

We denote by ∈^× the orthonormal eigenvectors of  and by 0 = _1 ≤…≤_n the ordered real eigenvalues of . We have = ^, where := (_1, …, _n) ∈^×. The matrix  is the graph Fourier basis <cit.>. For any signal ∈^ defined on the nodes of the graph , = ^ contains the Fourier coefficients of  ordered in increasing frequencies. This work deals with k-bandlimited signals ∈^ on , , signals whose Fourier coefficients _k+1, …, _ are null. Let _ be the restriction of  to its first  vectors:


    _ := ( _1, …, _) ∈^×.



 A signal ∈^ defined on the nodes of the graph  is k-bandlimited with ∈∖{0} if ∈(U_), , there exists η⃗∈^k such that


    = _η⃗.




This definition was also used in <cit.>. We assume that _≠_+1 to avoid any ambiguity in the definition of -bandlimited signals.

Finally, for any matrix X∈^_1 ×_2, X_2 denotes its spectral norm and X_F its Frobenius norm; when n_1=n_2, λ_ max(X) denotes its largest eigenvalue and λ_ min(X) its smallest eigenvalue. For any vector ∈^_1, _2 denotes the Euclidean norm of . Depending on the context, _j may represent the j^þ entry of the vector  or the j^þ column-vector of the matrix X. The entry on the i^þ row and j^þ column of X is denoted by X_ij. The identity matrix is denoted by I (its dimensions are determined by the context).

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










§ SAMPLING USING GROUPS OF NODES



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





 §.§ Grouping the nodes


We consider that the  nodes of  are divided into  different groups _1, …, _⊆{1, …, }. The size of the ℓ^þ group is denoted _ℓ. We suppose that these groups form a partition of {1, …, }, so that each node belongs to exactly one group. We have


    ∪_ℓ=1^N  _ℓ     =  {1, …, },
     and _ℓ ∩ _ℓ' = ∅.


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 _ℓ = {_1^(ℓ), …, __ℓ^(ℓ)}, we associate a matrix ^(ℓ) ∈ ^_ℓ× that restricts a graph signal to the nodes appearing in _ℓ, ,


    ^(ℓ)_ij := 
    {[                1 for  j = _i^(ℓ),;                0       otherwise. ].


Note that


    ∑_ℓ=1^^(ℓ)^^(ℓ) = I.



The case of overlapping groups can be handled by changing the definition of ^(ℓ) to


    ^(ℓ)_ij := 
    {[    β__i^(ℓ)^-1/2 for  j = _i^(ℓ),;                0       otherwise, ].


where 1 ≤β_i≤, i = 1, …,, is the number of times node i appears in the different groups _1, …, _. Equation (<ref>) also holds in this case. All results presented in Section <ref> are valid for overlapping groups with this definition of ^(ℓ).





 §.§ Sampling the groups


The sampling procedure consists in selecting  groups out of the  available ones. In the application of Section <ref>, it corresponds to the selection of the superpixels to label. We select these groups at random using a sampling distribution on { 1, …, } represented by a vector ∈^. The probability of selecting the ℓ^þ group is _ℓ. We assume that _ℓ > 0 for all ℓ = 1, …,, so that all groups may be selected with a non-zero probability. We obviously have ∑_ℓ=1^_ℓ = 1. 


The indices Ω := {ω_1, …, ω_} of the selected groups are obtained by drawing independently – thus with replacements –  indices from the set {1, …, } according to , ,


    (ω_j = ℓ) = _ℓ,  ∀ j ∈{1, …, } and ∀ℓ∈{1, …, }.


The selected groups are _ω_1, …, _ω_ and the total number of selected nodes is


    := ∑_j=1^_ω_j.


Once the groups are selected, we build the sampling matrix ∈^× that satisfies


    := 
    (
    [ ^(ω_1);      ⋮;  ^(ω_) ]),


and which restricts any signal to the nodes belonging to the selected groups. For a signal ∈^n defined on the nodes of , its sampled version ∈^ satisfies


    := .



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


    P^(ℓ) := _ℓ^-1/2 I∈^_ℓ×_ℓ,


to each group _ℓ. Then, once the groups are drawn, we construct the block diagonal matrix P∈^×


    P := ( P^(ω_1), …, P^(ω_)).


This matrix takes into account the probability of sampling each group and will be used to rescale  for norm preservation. This matrix ensures that ^-1 _ΩP_2^2 = ^-1 _ΩP_2^2 = _2^2.[This property is a consequence of (<ref>), proved in Appendix <ref>.] Both matrices P and M depend on Ω and are random.





 §.§ Group graph coherence


The sampling procedure in <cit.> 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 -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 -bandlimited signals spreads over the groups. We now introduce this quantity.

The matrix ^(ℓ)_ is the matrix that restricts a -bandlimited signal to the nodes belonging to _ℓ. Therefore, 


    ^(ℓ)__2 = sup_η⃗∈^k : η⃗_2=1^(ℓ)_η⃗_2


measures the energy on the nodes _ℓ of the normalised -bandlimited signal that is most concentrated on _ℓ. This energy varies between 0 and 1. When this energy is close to 1, there exists a -bandlimited signal whose energy is essential concentrated on _ℓ. This signal lives only on the nodes in _ℓ and does not spread elsewhere. On the contrary, when this energy is close to 0, there is no -bandlimited signal living only on _ℓ.

The sampling distribution  is adapted to the graph and the structure of the groups if: whenever ^(ℓ)__2 is high, _ℓ is high; whenever ^(ℓ)__2 is small, _ℓ is small. In other words, the ratio between ^(ℓ)__2 and _ℓ should be as constant as possible. This ensures that the groups where some -bandlimited signals are concentrated are sampled with higher probability. Similarly to what was done in <cit.> with individual nodes, we define the group graph weighted coherence as the largest ratio between ^(ℓ)__2 and _ℓ^-1/2.


The group graph cumulative coherence of order  is


    _ := max_1 ≤ℓ≤{_ℓ^-1/2^(ℓ)__2 }.


The quantity ^(ℓ)__2 is called the local group graph coherence.


In the extreme case where the groups _1, …, _ reduce all to singletons, we recover the definition of the graph weighted coherence introduced in <cit.>. 

It is easy to prove that _ is lower bounded by 1. Indeed, for any η⃗∈^ with η⃗_2 = 1, we have


    1 
        = _k η⃗_2^2 
    = ∑_ℓ=1^N^(ℓ)_k η⃗_2^2
    = ∑_ℓ=1^N_ℓ·^(ℓ)_k η⃗_2^2/_ℓ≤_1 ·max_1 ≤ℓ≤{^(ℓ)_k η⃗_2^2/_ℓ}
        = max_1 ≤ℓ≤{^(ℓ)_k η⃗_2^2/_ℓ}≤max_1 ≤ℓ≤{^(ℓ)_k_2^2/_ℓ}
    = _^2.


We have _ = 1 in, for example, the degenerated case where _1 = {1, …, }. 





 §.§ Stable embedding


We now have all the tools to present our main theorem that shows that sampling O(_^2 log()) groups is sufficient to stably embed the whole set of -bandlimited signals. Hence, it is possible to reconstruct any ∈(_k) from its measurements =.




Let  be a random subsampling matrix constructed as in (<ref>) using the groups _1, …, _ and the sampling distribution . For any δ, ξ∈ (0, 1), with probability at least 1-ξ,


    (1 - δ) _2^2 ≤P _2^2 ≤ (1 + δ) _2^2


for all ∈(_) provided that


    ≥3/δ^2 _^2  log( 2/ξ).






See Appendix <ref>.



In the above theorem, we recall that  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 _ is 1, sampling O(log()) groups might already be sufficient if the groups and the sampling distribution are well-designed.

The number of groups to sample is driven by _, which itself depends on the structure of the groups _1, …, _ and on the sampling distribution . 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 /L groups (L≥1) such that ^(ℓ)_k_2 ≈ (L/)^1/2 – , no -bandlimited signals have more than 100 · (L/)^1/2 percent of its energy concentrated in each group – then setting _ℓ = L/, ℓ = 1, …, /L, would yield ≈ 1. In this case, sampling one group is enough to embed the set of -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  to minimise _.

The sampling distribution minimizing the coherence _ is the distribution ^* ∈^ that satisfies


    ^*_ℓ := ^(ℓ)_k_2^2/∑_ℓ'=1^^(ℓ')_k_2^2,  for all ℓ∈{1, …, },


and for which


    _^*^2 = ∑_ℓ=1^^(ℓ)_k_2^2.


Indeed, let p⃗' ≠^* be another sampling distribution. As the entries of p⃗' and ^* are nonnegative and sum to 1, we necessarily have p⃗'_ℓ' < _ℓ'^* for some ℓ'>0. Then, 


    _p⃗'^2
    ≥p⃗'_ℓ'^-1^(ℓ')_k_2^2
    > 
    _ℓ'^*^-1^(ℓ')_k_2^2
    =
    _^*^2,


where the last equality is obtained by replacing _ℓ'^* with its value. Therefore, _p⃗' > _^* for any p⃗' ≠^*. As similar proof can be found in, , <cit.> where the authors derive the optimal sampling distribution for a compressive system.

We notice that


    _^*^2 = ∑_ℓ=1^^(ℓ)_k_2^2 ≤∑_ℓ=1^^(ℓ)_k_F^2 = .


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


    _^*^2 = ∑_ℓ=1^^(ℓ)_k_2^2 ≤,


as ^(ℓ)_k_2^2 ≤ 1. Therefore, in any case, the bound never suggests to sample much more than  groups, as one would expect it.

We recall that the results in <cit.> proves that it is always sufficient to sample O(log()) nodes to embed the set of -bandlimited signals. When sampling the nodes by groups, it is the number of groups to sample that should be O(log()). 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 <cit.> and structured sampling <cit.> are also important topics in compressed sensing. The method proposed here is closely inspired by these studies, especially by <cit.>. 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 ^*.





 §.§ A more practical result at small k's


Optimising the sampling distribution reduces to estimating the spectral norm of the matrices ^(ℓ)_. We will present in Section <ref> a method avoiding the computation of _. 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


    _ := max_1 ≤ℓ≤{_ℓ^-1/2^(ℓ)__F }.


The only difference between _ and _ is that we substituted the Frobenius norm for the spectral norm. As X_2 ≤X_F for any matrix X, we have _≤_; hence the results involving _ will be more pessimistic than those involving _. We have _≥. Indeed,


    k 
        = _k_F^2 
    = ∑_ℓ=1^N^(ℓ)_k_F^2
    = ∑_ℓ=1^N_ℓ·^(ℓ)_k_F^2/_ℓ≤_1 ·max_1 ≤ℓ≤{^(ℓ)_k_F^2/_ℓ}
    = _^2.


As with _, the lower bound is also attained in, for example, the degenerated case where _1 = {1, …, }. As _≤_, we have the following corollary to Theorem <ref>.


Let  be a random subsampling matrix constructed as in (<ref>) using the groups _1, …, _ and the sampling distribution . For any δ, ξ∈ (0, 1), with probability at least 1-ξ, (<ref>) holds for all ∈(_) provided that


    ≥3/δ^2 _^2  log( 2/ξ).







As _≥_, (<ref>) implies (<ref>). Theorem <ref> then proves that (<ref>) holds with probability at least 1-ξ for all ∈(_). 



The sufficient condition (<ref>) can be much more pessimistic than (<ref>). Indeed, as we have _^2 ≥, Condition (<ref>) suggests to always sample more than O(log()) groups, while we know that sampling O(log()) can be enough. The interest of this result is thus in the regime where  is small.

As we have done it with _, we can also optimise  to minimise _. The sampling distribution that minimises _ is the distribution q⃗^* satisfying


    q⃗^*_ℓ := ^(ℓ)_k_F^2/,  for all  l ∈{1, …, },


and for which


    _q⃗^*^2 = .


With this distribution, (<ref>) proves that sampling O(log()) groups is always enough to sample all -bandlimited signals. This result is particularly interesting at small 's because estimating q⃗^* is much easier and faster than estimating ^* (see Section <ref>).

In some particular cases, this result can be interesting  also at large 's. Indeed, (<ref>) is a pessimistic bound. In reality, we may have _q⃗^*^2 ≤ so that, according to (<ref>), fewer samples than O(log()) are actually sufficient when using q⃗^*. Furthermore, q⃗^* might actually be close to ^*, in which case we would reach almost optimal results with 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 _k to compute the distributions p⃗^* and q⃗^*. Computing _k can be intractable for large graphs. We will present a method to overcome this issue in Section <ref>. We continue first by explaining how to estimate  from its measurements.







§ FAST RECONSTRUCTION



Once the signal has been sampled, we shall also be able to reconstruct it. In <cit.>, the authors propose to estimate the original signal by solving


    min_z⃗∈^P(z⃗ - )_2^2 +  z⃗^ g() z⃗,


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


    g() := ∑_i=0^d α_i  ^i = ∑_i=0^d α_i  ^i ^,


where α_0, …, α_d∈ are the coefficients of the polynomial and d ∈ 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  and should ideally be a high-pass filter. The matrix 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 P and , which are both sparse matrices. The matrix g() does not have to be computed explicitly.

For brevity, we do not recall the theorem proving that (<ref>) provides accurate and stable recovery of all -bandlimited signals. However, this theorem applies as soon as the restricted isometry property (<ref>) holds, and thus applies here when (<ref>) holds. The reconstruction quality is controlled by g(_) and the ratio g(_)/g(_+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.





 §.§ Piecewise constant graph signals


Before continuing, we want to stress that we consider non-overlapping groups _1, …, _ in this rest of Section <ref>.

If a graph signal is nearly piecewise constant over the groups _1, …, _ 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  unknowns, we reduce the estimation to  unknowns. When ≪, this is a large reduction of dimension yielding a significant speed up and memory reduction. 

The fact that a signal  is piecewise constant over the groups _1, …, _ is characterized as follows. We construct the averaging row-vectors ^(ℓ)∈^1 × that satisfy


    ^(ℓ) := 1⃗^^(ℓ)/_ℓ^1/2,


and the matrix


    := 
    (
    [ ^(1);    ⋮;  ^() ]) ∈^×.


As the groups do not overlap, we have 

    ^ = I,

hence _2 = 1. Applying  to  provides  values, each one of them corresponding to the sum of the values of  within the group _ℓ, scaled by _ℓ^-1/2. Then, applying ^ to  gives an approximation of  where the values in the vector ^ are constant within each group; this is a piecewise constant vector over the groups. The value of ^ appearing within the group _ℓ is exactly the average of  within _ℓ. Saying that  is nearly constant within each group corresponds to assume that 


    ^ - _2 ≤ϵ_2,


where ϵ≥ 0 is a small value. The signal model of interest in this section is thus


    := {∈(_)   |  (^ - I)  _2 ≤ϵ_2 }.







 §.§ Reducing the dimension


To build a fast algorithm exploiting the above property, we use a reconstruction method similar to (<ref>) but involving vectors of smaller dimension. We define the averaged vector := ∈^ of dimension . As ∈, estimating  is enough to get a good approximation of  – we just need to multiply it with ^. Furthermore, as  is nearly piecewise constant over the groups _1, …, _, by construction of the matrix , the measurement vector = is also almost piecewise constant over the sampled groups _ω_1, …, _ω_. We thus average  over these groups by multiplying it with the matrix ∈^× that satisfies


    _ji := 
    {[                                 _ω_j^-1/2 for ∑_j'=1^j-1_ω_j'≤  i   ≤∑_j'=1^ℓ_ω_j',;                                         0                                otherwise. ].


We obtain


    :=  
    = ∈^.


We now have to link  to . We create the matrix ∈^× that restricts the  mean value of  to the  mean value of the selected groups, ,


    _j i := 
    {[            1 if  i = ω_j,;            0   otherwise. ].


We have therefore


    =  .



The goal is now to estimate  from . To ensure that the reconstruction method is stable to measurement noise, we do not consider the perfect scenario above but instead the scenario where


    =  + ,


and ∈^ models noise. We now need a regularisation term to estimate . We obtain this term by reducing the dimension of the regularisation involving the Laplacian  in (<ref>). We compute


    :=   g()  ^ = ()   g()   ()^∈^×.


Note that  is a symmetric positive definite matrix. Like g(), it can thus be used as a regularisation. We thus propose to estimate  by solving


    min_z̃⃗̃∈^P (z̃⃗̃ - )_2^2 +  z̃⃗̃^  z̃⃗̃,


where >0 and P∈^× is the diagonal matrix with entries satisfying


    P_jj := _ω_j^-1/2.


Let ^* ∈^ be a solution of (<ref>). We finally obtain an estimation of  by computing ^^*.


In the particular case where g(·) is the identity, one can notice that 


    _ℓℓ' = ∑_(i,j) ∈_ℓ×_ℓ'_ij/_ℓ^1/2_ℓ'^1/2


is non-zero only if there is at least one edge in E joining the groups _ℓ and _ℓ'. The Laplacian _ℓℓ' preserves the connections present in the original graph represented by . 

The dimension of the unknown vector in (<ref>) is , which can be much smaller than . This leads to a large gain in memory and computation time when either  or matrix-vector multiplications with  can be computed rapidly. If g has a small degree, then  can be computed explicitly in a short amount of time. In such a case, we can solve (<ref>) faster than (<ref>) 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  without temporary going back to the signal domain of dimension , , multiplying the vector z̃⃗̃ with ^, filtering the high dimensional signal ^z̃⃗̃, and downsampling the result. Even though solving (<ref>) might still be faster than solving (<ref>) 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 (<ref>) can be used to initialise the algorithm used to solve (<ref>) with a good approximate solution. As in multigrid approaches to solve linear systems of equations, see, , <cit.>.

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




Let Ω = {ω_1, …, ω_} be a set of  indices selected independently from { 1, …, } using a sampling distribution ∈^, , P, , P be the associated matrices constructed respectively in (<ref>), (<ref>), (<ref>) and (<ref>), and M_ max>0 be a constant such that P_2 ≤ M_ max. Let ξ, δ∈ (0, 1) and suppose that  satisfies (<ref>). With probability at least 1-ξ, the following holds for all ∈, all ∈^, all > 0, and all nonnegative nondecreasing polynomial functions g such that g(_+1) > 0.

Let ^* be the solution of (<ref>) with =  +. Define α⃗^* := __^ ^^* and β⃗^* := (I - __^)  ^^*. Then,


    α⃗^* - _2
     ≤ √( (1-δ)) ·[    ( 2 +  M_ max/√( g(_+1)))P_2
    + ( M_ max√(g(_)/g(_+1)) + √( g(_))) _2 
       . +  ϵ( 2 M_ max + M_ max√(g(_)/g(_+1)) + √( g(_))) _2 ],


and


    β⃗^*_2 
    ≤ 1/√( g(_+1))P_2   +  √(g(_)/g(_+1))_2   +  ϵ √(g(_)/g(_+1))_2.






See Appendix <ref>.



The vector α⃗^* is the orthogonal projection of ^^* onto (_). The vector β⃗^* is the projection of ^^* onto the orthogonal complement of (_). There are several remarks to make about the above theorem:




  * Theorem <ref> shows that the result obtained via (<ref>) is similar to the one we would have obtained by solving (<ref>) – see <cit.> for the error bounds – with additional errors controlled by ϵ. We recall that ϵ characterises how far  is from a piecewise constant signal. As expected, the smaller ϵ, the better the reconstruction.


  * The reconstruction quality improves when g(_) and the ratio g(_)/g(_+1) go to 0, and when g(_)/g(_+1) tends to 1. We recall that we have g(_) ≥ g(_+1) > g(_) by assumption.


  * The effect of the noise  decreases when g(_+1) increases, and, obviously, γ 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 <cit.>, where a multiresolution analysis method of graph signals is proposed. The coarsening method is however different than the one used here.







§ OPTIMAL SAMPLING ESTIMATION



In this section, we come back to the sampling process of Section <ref> and leave the reconstruction problem. We explain how to estimate the sampling distributions ^* and q⃗^* of Section <ref> without computing the truncated Fourier matrix _, as this computation is intractable for large graphs. The methods below only involve matrix-vector multiplications with the sparse Laplacian matrix  and are thus computationally tractable even for large . 





 §.§ Estimation of ^*


The distribution ^*, defined in (<ref>), that minimises the coherence _  is entirely defined by the values of ^(ℓ)_k_2 for ℓ = 1, …,, which are thus the quantities we need to evaluate.

We recall that, in graph signal processing, a filter is represented by a function h: →, and that the signal  filtered by h is


    _h :=  (ĥ⃗̂)  ^∈^,


where ĥ⃗̂ = (h(_1), …, h(_))^∈^. To filter the signal  without actually computing the graph Fourier transform of , we can approximate the function h by a polynomial 


    r(t) := ∑_i=0^d α_i   t^i ≈ h(t)


of degree d, and compute _r instead of _h. The filtered signal _r is computed rapidly using the formula


    _r = ∑_i=0^d α_i   (_1^i, …, _^i)  ^ = ∑_i=0^d α_i  ^i ,


that involves only matrix-vector multiplications with the sparse Laplacian matrix . We let the reader refer to <cit.> for more information on this fast filtering technique. For any polynomial function r of the form above and any matrix A∈^×, we define


    r(A) := ∑_i=0^d α_i  A^i .


Note that r() =   r()  ^.


Let i__k : → be the ideal low-pass filter at cutoff frequency _k, , the filter that satisfies


    i__k(t)
    =
    {[          1 if  t ≤_k,;          0 otherwise. ].


We have _k _k^ = i__k(). Then, we notice that


    ^(ℓ)_k_2^2
    =
    ^(ℓ)_k _k^^(ℓ)^_2
    =
    ^(ℓ)  i__k()  ^(ℓ)^_2.


We recall that ^(ℓ) is the matrix that restricts the signal to the nodes belonging to _ℓ. 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 0 in all groups ℓ' ≠ℓ, 2) low-pass filters the extended signal, 3) restricts the result to the group _ℓ. This process can be approximated by replacing the ideal low-pass filter i__k with a polynomial approximation ĩ__k of i__k and 


    ^(ℓ)  i__k()  ^(ℓ)^_2
    ≈^(ℓ) ĩ__k()  ^(ℓ)^_2.


To estimate ^*, 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 ĩ__k.

Estimating ^(ℓ) ĩ__k()  ^(ℓ)^_2 amounts to computing the largest eigenvalue of ^(ℓ) ĩ__k()  ^(ℓ)^ which can be done, , by using the power method. This method requires matrix-vector multiplication only with ^(ℓ) and ĩ__k() and is thus fast.
Finally, the approximation p̅⃗̅∈^ of ^* satisfies


    _ℓ := λ_ max(^(ℓ) ĩ__k()  ^(ℓ)^)/∑_ℓ'=1^λ_ max(^(ℓ') ĩ__k()  ^(ℓ')^).



Note that an estimation of _ is required beforehand to define the filter ĩ__k. We estimate this value using the dichotomy method presented in <cit.>.





 §.§ Estimation of q⃗^*


Computing  requires the estimation of  eigenvalues. Even though these estimations can be done in parallel, this process might still be too slow for certain applications. As explained before, when  is small, we can use the sampling distribution q⃗^* in (<ref>) that minimises _p⃗. This distribution is faster to compute than .

We start by noticing that we have


    ^(ℓ)_k_F^2 = ∑_i ∈_ℓ_^δ⃗_i_2^2


for each group ℓ = 1, …,. The vector δ⃗_i ∈^ is the unit vector that is null on all nodes except at node i. Hence, we only need an estimation of _^δ⃗_i_2^2, i = 1, …,, to estimate q⃗^*. An algorithm was already proposed in <cit.> to estimate these values. We let the reader refer to Algorithm 1 in <cit.> for the details of the method. We just recall that this estimation is obtained by filtering O(log()) random signals with a polynomial approximation of i__. Finally, our estimation q̅⃗̅∈^ of q⃗^* has entries


    q̅⃗̅_ℓ := ∑_i ∈_ℓ_^δ⃗_i_F^2/∑_ℓ=1^∑_i ∈_ℓ_^δ⃗_i_F^2,


where each _^δ⃗_i_2^2 are estimated by Algorithm 1 in <cit.>.

This estimation is faster than for  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 _ℓ. In total, (much) more than _ℓ filterings are thus required. On the contrary, for q̅⃗̅, we just need to filter O(log()) signals to obtain the estimation. In most situations, we already have O(log()) ≤_ℓ and computing q̅⃗̅ is thus faster than computing .








§ EXPERIMENTS



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 <ref>.





 §.§ Sampling distributions






We perform experiments on two different graphs: the Minnesota graph of size = 2642 and the bunny graph of size = 2503. Both graphs are presented in Fig. <ref> and are available in the GSP toolbox <cit.>. 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 = 73 groups represented in Fig. <ref>. For the bunny graph, we obtain = 213 groups with a similar procedure (see Fig. <ref>).

For each graph, we compute the combinatorial Laplacian and _k for different values of k. Then, we compute the lower RIP constant, , the constant δ_k>0 that satisfies


    δ_k 
    = 1 -   inf_∈(_) 
    x_2 = 1 P _2^2.


This constant is the smallest value that δ can take such that the left-hand side of the RIP (<ref>) holds. Remark that


    δ_k = 1 -  λ_ min( _k^^P^2 _k).


We estimate δ_k for 500 independent draws of the set Ω, which defines the matrices P, and different numbers of selected groups . All samplings are done in the conditions of Theorem <ref> using the sampling distributions u⃗, p⃗^*, and q̅⃗̅. The vector u⃗ denotes the uniform distribution over {1, …, }. When conducting this experiment with the estimated distributions  and q̅⃗̅, we re-estimate these distributions at each of the 500 trials. These distributions are estimated using Jackson-Chebychev polynomials of order 50 <cit.>. For the Minnesota graph, we consider the bandlimits = 5, 10, 20. For the bunny graph, we consider the bandlimits = 10, 25, 50.

We present the probability that δ_k is less than 0.995, estimated over the 500 draws of Ω, as a function of  in Fig. <ref>. For the Minnesota graph, the performance is better when using the optimal distribution ^* than when using the uniform distribution u⃗ for all k, which is in line with the theory. The estimated  and q̅⃗̅ yield performance equivalent to ^*. This confirms that we can achieve similar sampling performance without having to compute the Fourier matrix _, which, we recall, is intractable for large graphs. This also shows that 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 , q̅⃗̅ or ^* instead of u⃗.

For illustration, we present in Fig. <ref> examples of computed sampling distributions p⃗^*,  and q̅⃗̅. All sampling distributions exhibit similar structures, which explains why they all yield about the same performance in our experiments.








 §.§ Object segmentation





  §.§.§ Protocol





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

The graph  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 × 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  feature vectors g⃗_i ∈^45, i=1, …,. We then connect each feature vector to its 9 nearest neighbours (in the Euclidean sense), which gives a set of 9 edges E. The adjacency matrix W∈^× satisfies
W_ij := exp [- g⃗_i  - g⃗_j_2^2/σ^2],
where σ > 0 is the 25^þ percentile of the set {g⃗_i  - g⃗_j_2 : (i, j) ∈E}. We finally symmetrise the matrix W and compute the combinatorial Laplacian ∈^×.

We study three strategies to choose the superpixels. The first strategy consists in choosing the superpixels uniformly at random, , using the sampling distribution u⃗. The second and third strategies consist in choosing the superpixels with respectively the optimised distributions q̅⃗̅ and , which we evaluate at _0 = 50 using Jackson-Chebychev polynomials of order 75 <cit.>. For illustration, we present in Fig. <ref> the estimated values ^(ℓ)__0_F^2 and ^(ℓ)__0_2^2, which define the optimised distributions q̅⃗̅ and p̅⃗̅, respectively. Both distributions indicate that one should label more superpixels around the tiger. The distribution q̅⃗̅ “focuses” however more on specific regions, like the head of tiger. The distribution  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  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 ∈^ from which we want to reconstruct . We repeat this procedure for ∈{50, 70, …, 250}. For each , 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 , we first use the fast reconstruction method (<ref>) and then refine the solution at the pixel level with (<ref>), using the solution of the first minimisation problem as initialisation to solve the second minimisation problem. We choose g() = and solve both problems in the limit where → 0. In this limit, the problems (<ref>) and (<ref>) become


    min_z̃⃗̃∈^ z̃⃗̃^  z̃⃗̃   subject to   z̃⃗̃ =


and


    min_z⃗∈^ z⃗^ g() z⃗   subject to   z⃗ = ,


respectively. Both problems are solved using FISTA <cit.>. The same stopping criteria are used for all experiments.

 



  §.§.§ Results


We present in the top panel of Fig. <ref> the reconstruction snr obtained with the different methods. The reconstruction snr is defined as -20log_10(- - ^*_2/_2), where ^* is the reconstructed signal.  We notice that the snr attained with the fast decoder (<ref>) is very similar to the snr attained with (<ref>). We also remark that the optimised distributions  and q̅⃗̅ yield better reconstructions than the uniform distribution u⃗. The mean reconstruction snr is slightly better with  than with q̅⃗̅ at s ≥ 150.

We present the computation time of each method in the bottom panel of Fig. <ref>. We notice that solving (<ref>) is much faster than solving (<ref>), 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 (<ref>) when the measurements are drawn with  or with q̅⃗̅ than with u⃗. The reason is probably a better initialisation of (<ref>) or a better “quality” of the measurements with the optimised distributions than with the uniform distribution.

Finally, we present in Fig. <ref> some examples of reconstructions from s=150 sampled superpixels for each method. We notice that the optimised sampling distributions improve the reconstruction of  around the head and tail of the tiger, , where the optimised distributions have higher values. With a uniform distribution, the structure of the graph makes it difficult to reconstruct  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.











§ DISCUSSION AND CONCLUSION


We presented a structured sampling strategy for -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 -bandlimited signals. Finally, we presented a fast reconstruction technique for -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 <cit.>. 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  clusters, and interpolating the clustering results on all nodes by solving (<ref>). To accelerate the method, one could 1) pre-group similar nodes to form  groups such that ≤≪, , by running few iterations of the k-means algorithm; 2) subsample this set of  groups; 3) cluster this subset to find  clusters; and 4) solve (<ref>) to cluster all nodes. If the overhead of computing the  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 . In some applications, the final result may change a lot depending on the value of  which was chosen to compute this distribution. Finding a range of values of  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 _λ_. 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 ^* and q⃗^*.







§ 
[hang]Appendix 0mm[]







§  - PROOF OF THE THEOREM <REF>



§ PROOF OF THE THEOREM <REF>



As done in <cit.>, the proof is obtained by applying the following lemma obtained by Tropp in <cit.>.




Consider a finite sequence {X_j } of independent, random, self-adjoint, positive semi-definite matrices of dimension d × d. Assume that each random matrix satisfies
λ_ max(X_j) ≤ R
almost surely. Define


    μ_ min := λ_ min( ∑_j  X_j )
     and μ_ max := λ_ max( ∑_j  X_j ).


Then


    {λ_ min( ∑_j X_j )  ≤ (1 - δ) μ_ min}   ≤ d  [ ^-δ/(1-δ)^1-δ]^μ_ min/R for δ∈ [0, 1],
    and  {λ_ max( ∑_j X_j )  ≥ (1 + δ) μ_ max}   ≤ d  [ ^δ/(1+δ)^1+δ]^μ_ max/R for δ≥ 0.





We will also use the facts that, for all δ∈ [0, 1],


    [ ^-δ/(1-δ)^1-δ]^μ_ min/R ≤ exp( - δ^2 μ_ min/3   R)
     and [ ^δ/(1+δ)^1+δ]^μ_ max/R ≤ exp( - δ^2 μ_ max/3   R).





We start by noticing that


    _^^PP_ = ∑_j = 1^(_^^(ω_j)^P^(ω_j)) (P^(ω_j)^(ω_l)_).


We define


    X_j := 1/(_^^(ω_j)^P^(ω_j)) (P^(ω_j)^(ω_j)_)
       and   X := ∑_j=1^X_j =  _^^P^2 _.


The matrix X is thus a sum of  independent, random, self-adjoint, positive semi-definite matrices. We are in the setting of Lemma <ref>. We continue by computing X_j and λ_ max(X_j).

The expected value of each X_j is


    X_j    = 
    [ 1/(_^^(ω_j)^P^(ω_j)) (P^(ω_j)^(ω_j)_) ]
    = 
     _^(∑_ℓ=1^_ℓ(^(ℓ)^P^(ℓ)) (P^(ℓ)^(ℓ))  ) _
        =  
     _^(∑_ℓ=1^^(ℓ)^^(ℓ)) _
    =
     _^_
    =
     I.


Therefore, 
λ_ min( ∑_j X_j) = 1
and
λ_ max( ∑_j X_j) = 1.
Furthermore, for all j = 1, …,, we have


    λ_ max (X_j)
    =
    X_j_2 
    ≤max_1 ≤ℓ≤P^(ℓ)^(ℓ)_/_2^2
    =
    1/ max_1 ≤ℓ≤{^(ℓ)__2^2/_ℓ}
    = _^2/.


Lemma <ref> yields, for any δ∈ (0, 1),


    {λ_ min( X)  ≤ (1 - δ) }    ≤ ·[ ^-δ/(1-δ)^1-δ]^/_^2 ≤  exp( - δ^2 /3  _^2) 
    and  {λ_ max( X)  ≥ (1 + δ) }    ≤ ·[ ^δ/(1+δ)^1+δ]^/_^2 ≤  exp( - δ^2 /3  _^2).


Therefore, for any δ∈ (0, 1), we have, with probability at least 1 - ξ,


    1 - δ≤λ_ min( X) 
      and  λ_ max( X) ≤ 1+δ


provided that (<ref>) holds. Noticing that (<ref>) implies that


    (1-δ) α⃗_2^2
    ≤P_α⃗_2^2 
    ≤ 
    (1+δ) α⃗_2^2,


for all α⃗∈^k, which is equivalent to (<ref>) for all ∈(_k), terminates the proof.









§  - PROOF OF THEOREM <REF>



§ PROOF OF THEOREM <REF>



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

The first useful property is


    P = P.


Indeed, for any z⃗∈^, the j^þ entry of Pz⃗ is


    ( Pz⃗)_j
    =
    1⃗^^(ω_j)z⃗/( _ω_j _ω_j)^1/2.


Then, the j^þ entry of Pz⃗ is the scaled sum of the values in the j^þ sampled group appearing in Pz⃗, which is _ω_j^-1/2^(ω_j)z⃗. From the definition of , the sum is scaled by _ω_j^-1/2. Therefore, ( Pz⃗ )_j = ( Pz⃗ )_j for all j ∈{1, …, }, which terminates the proof.

The second property is 


    ^ = I,


which implies _2 = 1.

The third property is


    P^z̃⃗̃_2 
    = 
    P^z̃⃗̃_2,


for all z̃⃗̃∈^. To prove this property, we remark that
P^z̃⃗̃
= 
^Pz̃⃗̃.
Indeed, the entries in P^z̃⃗̃ corresponding to the first selected group _ω_1 are all equal to (_ω_1 _ω_1)^-1/2z̃⃗̃_ω_1. There are _ω_1 such entries. One can also notice that the first _ω_1 entries of ^Pz̃⃗̃ are also equal to (_ω_j _ω_1)^-1/2z̃⃗̃_ω_1. Repeating this reasoning for all the sampled groups proves the equality. On the one side, we thus have
P^z̃⃗̃_2
= 
^Pz̃⃗̃_2
=
Pz̃⃗̃_2,
where we used (<ref>). On the other side, we have
P^z̃⃗̃_2
= 
P^z̃⃗̃_2
=
Pz̃⃗̃_2,
where we used (<ref>) and (<ref>). This terminates the proof.



As ^* is a minimiser of (<ref>), we have


    P(^* - )_2^2 +   (^*)^  ^*
     ≤ P( - )_2^2 +  ^  .


To prove the theorem, we need to lower and upper bound the left and right hand sides of (<ref>), respectively. We start with the bound involving  and then with the ones involving .

[Bounding the terms in (<ref>) involving ]. We define the matrices


    _    := ( _+1, …, _) ∈^× ( - ),
      G_    := ( g(_1), …, g(_) ) ∈^×,
    G̅_   := ( g(_+1), …, g(_) ) ∈^( - ) × ( - ).


By definition of α⃗^* and β⃗^*, ^^* =  α⃗^* + β⃗^* with α⃗^* ∈(_) and β⃗^* ∈(_). We recall that = ()   g()   ()^. Therefore, we obtain 


    (^*)^  ^*
    =
    (_^α⃗^*)^ G_  (_^α⃗^*)
    + (_^β⃗^*)^ G̅_  (_^β⃗^*)
    ≥
    g(_+1) β⃗^*_2^2.


In the first step, we used the facts that _^α⃗^* = 0⃗ and _^β⃗^* = 0⃗. The second step follows form the fact that _^β⃗^* _2 = β⃗^* _2. We also have


    ^      =
    (^^)^  g()   (^^)
     ≤  
    g(_) _^_2^2
    + 
    g(_) _^^_2^2 
       ≤
    g(_) ^_2^2
    + 
    g(_) _^^_2^2
     ≤ 
    g(_) _2^2
    + 
    g(_) [ _^_2^2
    + 
    _^ (^ - )_2^2 ] 
       ≤
    g(_) _2^2
    + 
    ϵ^2 g(_) _2^2.


The second inequality follows from the facts that _k_2 = 1 and =. To obtain the third inequality, we used _2 = 1 and the triangle inequality. For the last step, notice that _^_2 = 0 (as ∈(_)), _k_2 = 1 and use (<ref>). 

[Bounding the terms in (<ref>) involving ]. By definition of , it is immediate that


    P( - )_2^2
    =
    P_2^2.


For the other term involving ,
the triangle inequality yields


    P^* - P_2 
       ≥P^* - P_2 - P_2.


Then, we have


    P^* - P_2
        =
    P^^* - P_2
    = 
    P^^* - P_2 
       ≥P^^* - P^_2
    - P^ - P_2.


The first equality follows from ^ = I, the second from (<ref>), and the triangle inequality was used in the last step. To summarise, we are at


    P^* - P_2
    ≥P^^* - P^_2
    - P^ - P_2 - P_2.


We continue by lower bounding P^^* - P^_2 and upper bounding P^ - P_2 separately.

[Lower bound on P^^* - P^_2]. Equality (<ref>) yields


    P^^* - P^_2
    = 
    P^^* - P^_2.


Using the triangle inequality and the fact that ∈, we obtain

    P^^* - P^_2
       ≥P^^* - P_2 - P - P^_2 
       ≥P^^* - P_2 - P_2 ^ - _2 
       ≥P^^* - P_2 - ϵ  M_ max_2.


The restricted isometry property (<ref>) then yields


    P^^* - P_2
        =
    P  (α⃗^* - ) + Pβ⃗^*_2
    ≥P (α⃗^* - )_2 - Pβ⃗^*_2 
       ≥√( (1-δ)) α⃗^* - _2 - M_ maxβ⃗^*_2.


We used the equality ^^* = α⃗^* + β⃗^*. We thus proved that


    P^^* - P^_2
    ≥√( (1-δ)) α⃗^* - _2 - M_ maxβ⃗^*_2 - ϵ M_ max_2.



[Upper bounding P^ - P_2]. We have


    P^ - P_2
    ≤_2
    P_2
    ^ - _2
    ≤ϵ  M_ max_2.


We used the facts that _2 = 1 (see (<ref>)) and that ∈ in the last inequality.

Using the inequalities (<ref>) and (<ref>) in (<ref>), we arrive at


    P^* - P_2
    ≥√( (1-δ)) α⃗^* - _2 - M_ maxβ⃗^*_2
    - 2ϵ M_ max_2 - P_2.



[Finishing the proof]. Remark that (<ref>) implies


    P(^* - )_2^2 
    ≤P( - )_2^2 +  ^  ,
    and    (^*)^  ^*
    ≤P( - )_2^2 +  ^  .


Using (<ref>), (<ref>), (<ref>) in (<ref>), we obtain


    g(_+1) β⃗^*_2^2
    ≤ P_2^2 
    +    g(_) _2^2
    + ϵ^2    g(_) _2^2,


which implies (<ref>) in Theorem <ref>. It remains to prove (<ref>) to finish the proof.

Using (<ref>), (<ref>) and (<ref>) in (<ref>), we obtain


    √( (1-δ))α⃗^* - _2
    ≤
    2 P_2 + M_ maxβ⃗^*_2
    + (√( g(_)) + ϵ√( g(_)) + 2 ϵ M_ max) _2.


Using (<ref>) to bound β⃗^*_2 on the right hand side, we have


    √( (1-δ))α⃗^* - _2
        ≤  
    2 P_2 + M_ max/√( g(_+1))P_2
    + M_ max√(g(_)/g(_+1))_2 
        + ϵ   M_ max√(g(_)/g(_+1))_2
    + (√( g(_)) + ϵ√( g(_)) + 2 ϵ M_ max) _2  
        =
    ( 2 +  M_ max/√( g(_+1)))P_2
    + ( M_ max√(g(_)/g(_+1)) + √( g(_))) _2 
        + ϵ( 2 M_ max + M_ max√(g(_)/g(_+1)) + √( g(_))) _2.


We only rearranged the term in the last step. This proves (<ref>) and terminates the proof.







IEEEtran
imaiai


