Gilles Puy

Patrick Pérez

Structured sampling and fast reconstruction of smooth graph signals

come

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, e.g., propagate foreground/background scribbles for object segmentation [START_REF] Boykov | Interactive graph cuts for optimal boundary and region segmentation of objects in n-d images[END_REF][START_REF] Grady | Random walks for image segmentation[END_REF][START_REF] Chen | Manifold preserving edit propagation[END_REF][START_REF] Chen | Image matting with local and nonlocal smooth priors[END_REF][START_REF] Xu | A sparse control model for image and video editing[END_REF] or propagate manual color/tone modifications [START_REF] Chen | Manifold preserving edit propagation[END_REF][START_REF] Chen | Image matting with local and nonlocal smooth priors[END_REF][START_REF] Xu | A sparse control model for image and video editing[END_REF][START_REF] Levin | Colorization using optimization[END_REF][START_REF] Lischinski | Interactive local adjustment of tonal values[END_REF]. First, a graph G modelling the similarities between the pixels of the image is built. Second, the user-edits specify values at some nodes (pixels) of G. Finally, the complete signal is estimated by assuming that it is smooth on G. The quality of the propagation depends on the structure of 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 [START_REF] Chen | Manifold preserving edit propagation[END_REF]. We address these problems from a graph signal processing point-of-view [START_REF] Shuman | The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains[END_REF]. 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 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 [START_REF] Pesenson | Sampling in Paley-Wiener spaces on combinatorial graphs[END_REF][START_REF] Pesenson | Sampling, Filtering and Sparse Approximations on Combinatorial Graphs[END_REF]. 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., [START_REF] Anis | Towards a sampling theorem for signals on arbitrary graphs[END_REF][START_REF] Anis | Asymptotic justification of bandlimited interpolation of graph signals for semi-supervised learning[END_REF] and Chen et al., [START_REF] Chen | Sampling theory for graph signals[END_REF][START_REF] Chen | Discrete signal processing on graphs: Sampling theory[END_REF] 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 [START_REF] Anis | Efficient sampling set selection for bandlimited graph signals using graph spectral proxies[END_REF], 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 [START_REF] Chen | Discrete signal processing on graphs: Sampling theory[END_REF][START_REF] Chen | Signal recovery on graphs: Random versus experimentally designed sampling[END_REF][START_REF] Puy | Random sampling of bandlimited signals on graphs[END_REF]. Recently, Puy et al. [START_REF] Puy | Random sampling of bandlimited signals on graphs[END_REF] 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.

Contributions

The random sampling strategy that we propose generalises the method proposed in [START_REF] Puy | Random sampling of bandlimited signals on graphs[END_REF]. We use here a structured sampling strategy. Let us already acknowledge here that such strategies are also studied in the field of compressed sensing [START_REF] Chauffert | Variable density compressed sensing in mri. theoretical vs heuristic sampling strategies[END_REF][START_REF] Adcock | Breaking the coherence barrier: a new theory for compressed sensing[END_REF][START_REF] Boyer | Compressed sensing with structured sparsity and structured acquisition[END_REF][START_REF] Bigot | An analysis of block sampling strategies in compressed sensing[END_REF] 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 [START_REF] Puy | Random sampling of bandlimited signals on graphs[END_REF] 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.

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 [START_REF] Chen | Manifold preserving edit propagation[END_REF] for instance, we build a graph 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 G. To propose the regions to label, we view this segmentation problem as a sampling problem of a smooth signal on 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 [START_REF] Achanta | Slic superpixels compared to state-of-the-art superpixel methods[END_REF]. 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 [START_REF] Korinke | Superselect: An interactive superpixel-based segmentation method for touch displays[END_REF]. 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 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.

Notations and definitions

We consider that G = {V, E, W} is an undirected weighted graph, where V is the set of n nodes, E is the set of edges, and W ∈ R n×n is the weighted adjacency matrix with nonnegative entries. We denote the graph Laplacian by L ∈ R n×n . We assume that L is real, symmetric, and positive semi-definite, e.g., the combinatorial graph Laplacian L := D -W, or the normalised one L := I -D -1/2 WD -1/2 . The matrix D ∈ R n×n is the diagonal degree matrix and I ∈ R n×n is the identity matrix [START_REF] Chung | Spectral graph theory[END_REF]. The diagonal degree matrix D has entries

D ii := i =j W ij .
We denote by U ∈ R n×n the orthonormal eigenvectors of L and by 0 = λ 1 . . . λ n the ordered real eigenvalues of L. We have L = UΛU , where Λ := diag(λ 1 , . . . , λ n ) ∈ R n×n . The matrix U is the graph Fourier basis [START_REF] Shuman | The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains[END_REF]. For any signal x ∈ R n defined on the nodes of the graph G, x = U x contains the Fourier coefficients of x ordered in increasing frequencies. This work deals with k-bandlimited signals x ∈ R n on G, i.e., signals whose Fourier coefficients xk+1 , . . . , xn are null. Let U k be the restriction of U to its first k vectors:

U k := (u 1 , . . . , u k ) ∈ R n×k . (1) Definition 1.1 (k-bandlimited signal on G). A signal x ∈ R n defined on the nodes of the graph G is k-bandlimited with k ∈ N \ {0} if x ∈ span(U k ), i.e., there exists η ∈ R k such that x = U k η. (2)
This definition was also used in [START_REF] Anis | Towards a sampling theorem for signals on arbitrary graphs[END_REF][START_REF] Chen | Discrete signal processing on graphs: Sampling theory[END_REF][START_REF] Puy | Random sampling of bandlimited signals on graphs[END_REF]. We assume that λ k = λ k+1 to avoid any ambiguity in the definition of k-bandlimited signals.

Finally, for any matrix X ∈ R n1×n2 , 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 x ∈ R n1 , x 2 denotes the Euclidean norm of x. Depending on the context, x j may represent the j th entry of the vector x or the j th column-vector of the matrix X. The entry on the i th row and j th 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. 2 a representation of the important variables and processes involved in this paper in order to facilitate the understanding of the different results.

x 2 R n x 2 R N ỹ 2 R s à A y 2 R m M e M L 2 R n⇥n e L 2 R N ⇥N Averaging Sampling p 2 R N Figure 2.
Summary of the main variables and processes involved in the presented sampling strategies. All variables are defined when needed hereafter.

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 n nodes of G are divided into N different groups N 1 , . . . , N N ⊆ {1, . . . , n}. The size of the th group is denoted |N |. We suppose that these groups form a partition of {1, . . . , n}, so that each node belongs to exactly one group. We have

∪ N =1 N = {1, . . . , n}, and N ∩ N = ∅. (3)
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 N = {n

ij := 1 for j = n ( ) i , 0 otherwise. (4) Note that N =1 N ( ) N ( ) = I. ( 5
)
The case of overlapping groups can be handled by changing the definition of N ( ) to

N ( ) ij := β -1/2 n ( ) i for j = n ( ) i , 0 otherwise, (6)
where 1 β i N , i = 1, . . . , n, is the number of times node i appears in the different groups N 1 , . . . , N N . Equation ( 5) also holds in this case. All results presented in Section 2 are valid for overlapping groups with this definition of N ( ) .

Sampling 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, . . . , N } represented by a vector p ∈ R N . The probability of selecting the th group is p . We assume that p > 0 for all = 1, . . . , N , so that all groups may be selected with a non-zero probability. We obviously have

N =1 p = 1.
The indices Ω := {ω 1 , . . . , ω s } of the selected groups are obtained by drawing independently -thus with replacementss indices from the set {1, . . . , N } according to p, i.e., P(ω j = ) = p , ∀j ∈ {1, . . . , s} and ∀ ∈ {1, . . . , N }. [START_REF] Lischinski | Interactive local adjustment of tonal values[END_REF] The selected groups are N ω1 , . . . , N ωs and the total number of selected nodes is

m := s j=1 N ωj . (8)
Once the groups are selected, we build the sampling matrix M ∈ R m×n that satisfies

M :=    N (ω1)
. . .

N (ωs)    , (9)
and which restricts any signal to the nodes belonging to the selected groups. For a signal x ∈ R n defined on the nodes of G, its sampled version y ∈ R m satisfies y := Mx. [START_REF] Pesenson | Sampling, Filtering and Sparse Approximations on Combinatorial Graphs[END_REF] 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 M. To conduct this study, we need to define few more matrices. First, we associate the matrix

P ( ) := p -1/2 I ∈ R |N |×|N | , (11)
to each group N . Then, once the groups are drawn, we construct the block diagonal matrix P ∈ R m×m P := diag P (ω1) , . . . , P (ωs) . [START_REF] Anis | Asymptotic justification of bandlimited interpolation of graph signals for semi-supervised learning[END_REF] This matrix takes into account the probability of sampling each group and will be used to rescale y for norm preservation. This matrix ensures that s 1 Both matrices P and M depend on Ω and are random.

-1 E Ω Py 2 2 = s -1 E Ω PMx 2 2 = x 2 2 .

Group graph coherence

The sampling procedure in [START_REF] Puy | Random sampling of bandlimited signals on graphs[END_REF] 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 kbandlimited 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 N ( ) U k is the matrix that restricts a k-bandlimited signal to the nodes belonging to N . Therefore,

N ( ) U k 2 = sup η∈R k : η 2 =1 N ( ) U k η 2 (13)
measures the energy on the nodes N of the normalised k-bandlimited signal that is most concentrated on N . This energy varies between 0 and 1. When this energy is close to 1, there exists a k-bandlimited signal whose energy is essential concentrated on N . This signal lives only on the nodes in N and does not spread elsewhere. On the contrary, when this energy is close to 0, there is no k-bandlimited signal living only on N . The sampling distribution p is adapted to the graph and the structure of the groups if: whenever

N ( ) U k 2 is high, p is high; whenever N ( ) U k 2 is small, p is small.
In other words, the ratio between N ( ) U k 2 and p 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 [START_REF] Puy | Random sampling of bandlimited signals on graphs[END_REF] with individual nodes, we define the group graph weighted coherence as the largest ratio between N ( ) U k 2 and p -1/2 . Definition 2.1 (Group graph cumulative coherence). The group graph cumulative coherence of order k is

ν p := max 1 N p -1/2 N ( ) U k 2 . ( 14
)
The quantity N ( ) U k 2 is called the local group graph coherence.

In the extreme case where the groups N 1 , . . . , N N reduce all to singletons, we recover the definition of the graph weighted coherence introduced in [START_REF] Puy | Random sampling of bandlimited signals on graphs[END_REF].

It is easy to prove that ν p is lower bounded by 1. Indeed, for any η ∈ R k with η 2 = 1, we have

1 = U k η 2 2 = N =1 N ( ) U k η 2 2 = N =1 p • N ( ) U k η 2 2 p p 1 • max 1 N N ( ) U k η 2 2 p = max 1 N N ( ) U k η 2 2 p max 1 N N ( ) U k 2 2 p = ν 2 p . ( 15
)
We have ν p = 1 in, for example, the degenerated case where N 1 = {1, . . . , n}.

Stable embedding

We now have all the tools to present our main theorem that shows that sampling O(ν 2 p log(k)) groups is sufficient to stably embed the whole set of k-bandlimited signals. Hence, it is possible to reconstruct any x ∈ span(U k ) from its measurements y = Mx. Theorem 2.2 (Restricted isometry property -RIP). Let M be a random subsampling matrix constructed as in (9) using the groups N 1 , . . . , N N and the sampling distribution p. For any δ, ξ ∈ (0, 1), with probability at least 1ξ,

(1 -δ) x 2 2 1 s PM x 2 2 (1 + δ) x 2 2 (16)
for all x ∈ span(U k ) provided that

s 3 δ 2 ν 2 p log 2k ξ . ( 17
)
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 ν 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 ν p , which itself depends on the structure of the groups N 1 , . . . , N N and on the sampling distribution 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 1) such that N ( ) U k 2 ≈ (L/n) 1/2i.e., no k-bandlimited signals have more than 100 • (L/n) 1/2 percent of its energy concentrated in each group -then setting p = L/n, = 1, . . . , n/L, would yield ν ≈ 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 p to minimise ν p .

The sampling distribution minimizing the coherence ν p is the distribution p * ∈ R N that satisfies

p * := N ( ) U k 2 2 N =1 N ( ) U k 2 2
, for all ∈ {1, . . . , N }, [START_REF] Chauffert | Variable density compressed sensing in mri. theoretical vs heuristic sampling strategies[END_REF] and for which

ν 2 p * = N =1 N ( ) U k 2 2 . ( 19
)
Indeed, let p = p * be another sampling distribution. As the entries of p and p * are nonnegative and sum to 1, we necessarily have p < p * for some > 0. Then,

ν 2 p p -1 N ( ) U k 2 2 > p * -1 N ( ) U k 2 2 = ν 2 p * , (20)
where the last equality is obtained by replacing p * with its value. Therefore, ν p > ν p * for any p = p * . As similar proof can be found in, e.g., [START_REF] Chauffert | Variable density compressed sensing in mri. theoretical vs heuristic sampling strategies[END_REF] where the authors derive the optimal sampling distribution for a compressive system.

We notice that

ν 2 p * = N =1 N ( ) U k 2 2 N =1 N ( ) U k 2 F = k. (21)
Hence, by using this distribution, [START_REF] Puy | Random sampling of bandlimited signals on graphs[END_REF] shows that sampling O(k log(k)) groups is always sufficient to sample all k-bandlimited signals. The exact number is proportional to ν 2 p * 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

ν 2 p * = N =1 N ( ) U k 2 2 N, (22) as N ( ) U k 2 2
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 [START_REF] Puy | Random sampling of bandlimited signals on graphs[END_REF] 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 [START_REF] Adcock | Breaking the coherence barrier: a new theory for compressed sensing[END_REF][START_REF] Puy | On variable density compressive sampling[END_REF][START_REF] Krahmer | Stable and robust sampling strategies for compressive imaging[END_REF] and structured sampling [START_REF] Chauffert | Variable density compressed sensing in mri. theoretical vs heuristic sampling strategies[END_REF][START_REF] Boyer | Compressed sensing with structured sparsity and structured acquisition[END_REF][START_REF] Bigot | An analysis of block sampling strategies in compressed sensing[END_REF] are also important topics in compressed sensing. The method proposed here is closely inspired by these studies, especially by [START_REF] Boyer | Compressed sensing with structured sparsity and structured acquisition[END_REF][START_REF] Bigot | An analysis of block sampling strategies in compressed sensing[END_REF]. 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 p * .

A more practical result at small k's

Optimising the sampling distribution reduces to estimating the spectral norm of the matrices N ( ) U k . We will present in Section 4 a method avoiding the computation of 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 νp := max

1 N p -1/2 N ( ) U k F . (23)
The only difference between ν p and νp is that we substituted the Frobenius norm for the spectral norm. As X 2 X F for any matrix X, we have ν p νp ; hence the results involving νp will be more pessimistic than those involving ν p . We have νp k. Indeed,

k = U k 2 F = N =1 N ( ) U k 2 F = N =1 p • N ( ) U k 2 F p p 1 • max 1 s N ( ) U k 2 F p = ν2 p . ( 24
)
As with ν p , the lower bound is also attained in, for example, the degenerated case where N 1 = {1, . . . , n}. As ν p νp , we have the following corollary to Theorem 2.2.

Corollary 2.3. Let M be a random subsampling matrix constructed as in (9) using the groups N 1 , . . . , N N and the sampling distribution p. For any δ, ξ ∈ (0, 1), with probability at least 1ξ, (16) holds for all x ∈ span(U k ) provided that

s 3 δ 2 ν2 p log 2k ξ . (25)
Proof. As νp ν p , (25) implies [START_REF] Puy | Random sampling of bandlimited signals on graphs[END_REF]. Theorem 2.2 then proves that [START_REF] Chen | Signal recovery on graphs: Random versus experimentally designed sampling[END_REF] holds with probability at least 1ξ for all x ∈ span(U k ).

The sufficient condition (25) can be much more pessimistic than [START_REF] Puy | Random sampling of bandlimited signals on graphs[END_REF]. Indeed, as we have ν2 p k, Condition (25) 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 ν p , we can also optimise p to minimise νp . The sampling distribution that minimises νp is the distribution q * satisfying q

* := N ( ) U k 2 F
k , for all l ∈ {1, . . . , N }, [START_REF] Puy | On variable density compressive sampling[END_REF] and for which

ν2 q * = k. ( 27
)
With this distribution, [START_REF] Chung | Spectral graph theory[END_REF] 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 q * is much easier and faster than estimating p * (see Section 4).

In some particular cases, this result can be interesting also at large k's. Indeed, ( 25) is a pessimistic bound. In reality, we may have ν 2 q * k so that, according to [START_REF] Puy | Random sampling of bandlimited signals on graphs[END_REF], fewer samples than O(k log(k)) are actually sufficient when using q * . Furthermore, q * might actually be close to p * , 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 U k to compute the distributions p * and q * . Computing 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 x from its measurements.

Fast reconstruction

Once the signal has been sampled, we shall also be able to reconstruct it. In [START_REF] Puy | Random sampling of bandlimited signals on graphs[END_REF], the authors propose to estimate the original signal by solving min z∈R n P(Mzy)

2 2 + γ z g(L)z, (28)
where γ > 0 and g is a nonnegative nondecreasing polynomial. We have

g(L) := d i=0 α i L i = d i=0 α i UΛ i U , ( 29
)
where α 0 , . . . , α d ∈ R are the coefficients of the polynomial and d ∈ 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 G 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 PM and L, which are both sparse matrices. The matrix g(L) does not have to be computed explicitly.

For brevity, we do not recall the theorem proving that (28) provides accurate and stable recovery of all k-bandlimited signals. However, this theorem applies as soon as the restricted isometry property (16) holds, and thus applies here when [START_REF] Puy | Random sampling of bandlimited signals on graphs[END_REF] holds. The reconstruction quality is controlled by g(λ k ) and the ratio g(λ k )/g(λ 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.

Piecewise constant graph signals

Before continuing, we want to stress that we consider non-overlapping groups N 1 , . . . , N N in this rest of Section 3.

If a graph signal is nearly piecewise constant over the groups N 1 , . . . , 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 n, this is a large reduction of dimension yielding a significant speed up and memory reduction.

The fact that a signal x is piecewise constant over the groups N 1 , . . . , N N is characterized as follows. We construct the averaging row-vectors a ( ) ∈ R 1×n that satisfy

a ( ) := 1 N ( ) |N | 1/2 , (30)
and the matrix

A :=   
a (1) . . .

a (N )    ∈ R N ×n . ( 31
)
As the groups do not overlap, we have

AA = I, (32)
hence A 2 = 1. Applying A to x provides N values, each one of them corresponding to the sum of the values of x within the group N , scaled by |N | -1/2 . Then, applying A to Ax gives an approximation of x where the values in the vector A Ax are constant within each group; this is a piecewise constant vector over the groups. The value of A Ax appearing within the group N is exactly the average of x within N . Saying that x is nearly constant within each group corresponds to assume that

A Ax -x 2 x 2 , ( 33
)
where 0 is a small value. The signal model of interest in this section is thus

A := {x ∈ span(U k ) | (A A -I) x 2 x 2 } . (34)

Reducing the dimension

To build a fast algorithm exploiting the above property, we use a reconstruction method similar to (28) but involving vectors of smaller dimension. We define the averaged vector x := Ax ∈ R N of dimension N . As x ∈ A , estimating x is enough to get a good approximation of x -we just need to multiply it with A . Furthermore, as x is nearly piecewise constant over the groups N 1 , . . . , N N , by construction of the matrix M, the measurement vector y = Mx is also almost piecewise constant over the sampled groups N ω1 , . . . , N ωs . We thus average y over these groups by multiplying it with the matrix A ∈ R s×m that satisfies

A ji := N ωj -1/2 for j-1 j =1 N ω j i j =1 N ω j , 0 otherwise. ( 35
)
We obtain ỹ := Ay = AMx ∈ R s . [START_REF] Tremblay | Compressive spectral clustering[END_REF] We now have to link ỹ to x. We create the matrix M ∈ R s×N that restricts the N mean value of x to the s mean value of the selected groups, i.e.,

M ji := 1 if i = ω j , 0 otherwise. ( 37
)
We have therefore

ỹ = M x. ( 38
)
The goal is now to estimate x 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 ỹ = M x + ñ, (39) and ñ ∈ R s models noise. We now need a regularisation term to estimate x. We obtain this term by reducing the dimension of the regularisation involving the Laplacian L in [START_REF] Borzi | Multigrid methods for pde optimization[END_REF]. We compute

L := A g(L) A = (AU) g(Λ) (AU) ∈ R N ×N . (40)
Note that L is a symmetric positive definite matrix. Like g(L), it can thus be used as a regularisation. We thus propose to estimate x by solving min

z∈R N P( M z -ỹ) 2 2 + γ z L z, ( 41
)
where γ > 0 and P ∈ R s×s is the diagonal matrix with entries satisfying

P jj := p -1/2 ωj . ( 42
)
Let x * ∈ R N be a solution of (41). We finally obtain an estimation of x by computing A x * .

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

L = (i,j)∈N ×N L ij |N | 1/2 |N | 1/2 (43)
is non-zero only if there is at least one edge in E joining the groups N and N . The Laplacian L preserves the connections present in the original graph represented by L.

The dimension of the unknown vector in (41) is N , which can be much smaller than n. This leads to a large gain in memory and computation time when either L or matrix-vector multiplications with L can be computed rapidly. If g has a small degree, then L can be computed explicitly in a short amount of time. In such a case, we can solve (41) faster than [START_REF] Borzi | Multigrid methods for pde optimization[END_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 L without temporary going back to the signal domain of dimension n, i.e., multiplying the vector z with A , filtering the high dimensional signal A z, and downsampling the result. Even though solving (41) might still be faster than solving [START_REF] Borzi | Multigrid methods for pde optimization[END_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 (41) can be used to initialise the algorithm used to solve [START_REF] Borzi | Multigrid methods for pde optimization[END_REF] with a good approximate solution. As in multigrid approaches to solve linear systems of equations, see, e.g., [START_REF] Borzi | Multigrid methods for pde optimization[END_REF][START_REF] Dreyer | Multigrid optimization in applications[END_REF].

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

Theorem 3.1. Let Ω = {ω 1 , . . . , ω s } be a set of s indices selected independently from {1, . . . , N } using a sampling distribution p ∈ R N , M, P, M, P be the associated matrices constructed respectively in (9), ( 12), ( 37) and (42), and M max > 0 be a constant such that PM 2 M max . Let ξ, δ ∈ (0, 1) and suppose that s satisfies [START_REF] Puy | Random sampling of bandlimited signals on graphs[END_REF]. With probability at least 1ξ, the following holds for all x ∈ A , all ñ ∈ R s , all γ > 0, and all nonnegative nondecreasing polynomial functions g such that g(λ k+1 ) > 0.

Let x * be the solution of (41) with ỹ = MAx + ñ. Define α * := U k U k A x * and β * := (I -

U k U k ) A x * . Then, α * -x 2 1 s(1 -δ) • 2 + M max γg(λ k+1 ) P ñ 2 + M max g(λ k ) g(λ k+1 ) + γg(λ k ) x 2 + 2M max + M max g(λ n ) g(λ k+1 ) + γg(λ n ) x 2 , ( 44
) and β * 2 1 γg(λ k+1 ) P ñ 2 + g(λ k ) g(λ k+1 ) x 2 + g(λ n ) g(λ k+1 ) x 2 . (45) Proof. See Appendix B.
The vector α * is the orthogonal projection of A x * onto span(U k ). The vector β * is the projection of A x * onto the orthogonal complement of span(U k ). There are several remarks to make about the above theorem:

• Theorem 3.1 shows that the result obtained via (41) is similar to the one we would have obtained by solving (28) -see [START_REF] Puy | Random sampling of bandlimited signals on graphs[END_REF] for the error bounds -with additional errors controlled by . We recall that characterises how far x is from a piecewise constant signal. As expected, the smaller , the better the reconstruction. • The reconstruction quality improves when g(λ k ) and the ratio g(λ k )/g(λ k+1 ) go to 0, and when g(λ n )/g(λ k+1 ) tends to 1. We recall that we have g(λ n ) g(λ k+1 ) > g(λ k ) by assumption. • The effect of the noise ñ decreases when g(λ k+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 [START_REF] Tremblay | Subgraph-based filterbanks for graph signals[END_REF], 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 2 and leave the reconstruction problem. We explain how to estimate the sampling distributions p * and q * of Section 2 without computing the truncated Fourier matrix U k , as this computation is intractable for large graphs. The methods below only involve matrix-vector multiplications with the sparse Laplacian matrix L and are thus computationally tractable even for large n.

Estimation of p *

The distribution p * , defined in [START_REF] Chauffert | Variable density compressed sensing in mri. theoretical vs heuristic sampling strategies[END_REF], that minimises the coherence ν p is entirely defined by the values of N ( ) U k 2 for = 1, . . . , 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 : R → R, and that the signal x filtered by h is

x h := U diag( ĥ) U x ∈ R n , ( 46
)
where ĥ = (h(λ 1 ), . . . , h(λ n )) ∈ R n . To filter the signal x without actually computing the graph Fourier transform of x, we can approximate the function h by a polynomial

r(t) := d i=0 α i t i ≈ h(t) (47)
of degree d, and compute x r instead of x h . The filtered signal x r is computed rapidly using the formula

x r = d i=0 α i U diag(λ i 1 , . . . , λ i n ) U x = d i=0 α i L i x, (48)
that involves only matrix-vector multiplications with the sparse Laplacian matrix L. We let the reader refer to [START_REF] Hammond | Wavelets on graphs via spectral graph theory[END_REF] for more information on this fast filtering technique. For any polynomial function r of the form above and any matrix A ∈ R n×n , we define

r(A) := d i=0 α i A i . (49) Note that r(L) = U r(Λ) U .
Let i λ k : R → R be the ideal low-pass filter at cutoff frequency λ k , i.e., the filter that satisfies

i λ k (t) = 1 if t λ k , 0 otherwise. (50) We have U k U k = i λ k (L). Then, we notice that N ( ) U k 2 2 = N ( ) U k U k N ( ) 2 = N ( ) i λ k (L) N ( ) 2 . ( 51
)
We recall that N ( ) is the matrix that restricts the signal to the nodes belonging to N . 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 N . This process can be approximated by replacing the ideal low-pass filter i λ k with a polynomial approximation ĩλ k of i λ k and

N ( ) i λ k (L) N ( ) 2 ≈ N ( ) ĩλ k (L) N ( ) 2 . ( 52
)
To estimate 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 ĩλ k .

Estimating N ( ) ĩλ k (L) N ( ) 2 amounts to computing the largest eigenvalue of N ( ) ĩλ k (L) N ( ) which can be done, e.g., by using the power method. This method requires matrix-vector multiplication only with N ( ) and ĩλ k (L) and is thus fast. Finally, the approximation p ∈ R N of p * satisfies

p := λ max (N ( ) ĩλ k (L) N ( ) ) N =1 λ max (N ( ) ĩλ k (L) N ( ) ) . ( 53
)
Note that an estimation of λ k is required beforehand to define the filter ĩλ k . We estimate this value using the dichotomy method presented in [START_REF] Puy | Random sampling of bandlimited signals on graphs[END_REF].

Estimation of q *

Computing 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 q * in (26) that minimises νp . This distribution is faster to compute than p.

We start by noticing that we have

N ( ) U k 2 F = i∈N U k δ i 2 2 (54)
for each group = 1, . . . , N . The vector δ i ∈ R n is the unit vector that is null on all nodes except at node i. Hence, we only need an estimation of U k δ i 2 2 , i = 1, . . . , n, to estimate q * . An algorithm was already proposed in [START_REF] Puy | Random sampling of bandlimited signals on graphs[END_REF] to estimate these values. We let the reader refer to Algorithm 1 in [START_REF] Puy | Random sampling of bandlimited signals on graphs[END_REF] 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 λ k . Finally, our estimation q ∈ R N of q * has entries

q := i∈N U k δ i 2 F N =1 i∈N U k δ i 2 F , (55)
where each U k δ i 2 2 are estimated by Algorithm 1 in [START_REF] Puy | Random sampling of bandlimited signals on graphs[END_REF]. This estimation is faster than for 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 N . In total, (much) more than N filterings are thus required. On the contrary, for q, we just need to filter O(log(n)) signals to obtain the estimation. In most situations, we already have O(log(n)) N and computing q is thus faster than computing p.

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

Sampling 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 [START_REF] Perraudin | Gspbox: A toolbox for signal processing on graphs[END_REF]. 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). Other columns: probability that δ k is less that 0.995 as a function of s. The curves in black, red, blue and green are obtained using the sampling distributions u, p * , q and p, respectively. The top row shows the result for the Minnesota graph. The bottom row shows the result for the bunny graph. The bandlimit k is indicated on top of each curve.

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

δ k = 1 - 1 s inf x∈span(U k ) x 2 =1 PM x 2 2 .
This constant is the smallest value that δ can take such that the left-hand side of the RIP ( 16) holds. Remark that

δ k = 1 - 1 s λ min U k M P 2 MU k . (56)
We estimate δ k for 500 independent draws of the set Ω, which defines the matrices PM, and different numbers of selected groups s. All samplings are done in the conditions of Theorem 2.2 using the sampling distributions u, p * , p and q. The vector u denotes the uniform distribution over {1, . . . , N }. When conducting this experiment with the estimated distributions p and q, we re-estimate these distributions at each of the 500 trials. These distributions are estimated using Jackson-Chebychev polynomials of order 50 [START_REF] Napoli | Efficient estimation of eigenvalue counts in an interval[END_REF]. 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 δ k is less than 0.995, estimated over the 500 draws of Ω, as a function of s in Fig. 3. For the Minnesota graph, the performance is better when using the optimal distribution p * than when using the uniform distribution u for all k, which is in line with the theory. The estimated p and q yield performance equivalent to p * . This confirms that we can achieve similar sampling performance without having to compute the Fourier matrix U k , 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 p, q or p * instead of u.

For illustration, we present in Fig. 4 examples of computed sampling distributions p * , 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 n = 321×481 presented in Fig. 1 for which our goal is to segment the tiger. The ground truth segmentation map x ∈ {0, 1} n is presented in Fig. 5. 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 available2 in [START_REF] Li | A benchmark for semantic image segmentation[END_REF]. Our objective is to recover the original map 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 [START_REF] Achanta | Slic superpixels compared to state-of-the-art superpixel methods[END_REF], 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 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 × 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 g i ∈ R 45 , i = 1, . . . , n. We then connect each feature vector to its 9 nearest neighbours (in the Euclidean sense), which gives a set of 9n edges E. The adjacency matrix W ∈ R n×n satisfies W ij := exp[-g ig j 2 2 /σ 2 ], where σ > 0 is the 25 th percentile of the set { g ig j 2 : (i, j) ∈ E}. We finally symmetrise the matrix W and compute the combinatorial Laplacian L ∈ R n×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 u. The second and third strategies consist in choosing the superpixels with respectively the optimised distributions q and p, which we evaluate at k 0 = 50 using Jackson-Chebychev polynomials of order 75 [START_REF] Napoli | Efficient estimation of eigenvalue counts in an interval[END_REF]. For illustration, we present in Fig. 5 the estimated values N ( ) U k0 2 F and N ( ) U k0 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 p spreads the measurements over the entire tiger more uniformly. [START_REF] Borzi | Multigrid methods for pde optimization[END_REF]. The continuous curves are obtained by solving (41). Right: mean computation time in seconds as a function of s. The dash-dotted curves are obtained by solving [START_REF] Borzi | Multigrid methods for pde optimization[END_REF]. The continuous curves are obtained by solving (41) using the result of (28) as initial point. In both graphs, the errorbars are at one standard deviation. The black curves are obtained with the uniform sampling distribution u. The blue curves are obtained with q. The red curves are obtained with p.

Results

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 (-xx * 2 / x 2 ), where x * is the reconstructed signal. We notice that the snr attained with the fast decoder (41) is very similar to the snr attained with [START_REF] Borzi | Multigrid methods for pde optimization[END_REF]. We also remark that the optimised distributions p and q yield better reconstructions than the uniform distribution u. The mean reconstruction snr is slightly better with p than with q at s 150.

We present the computation time of each method in the bottom panel of Fig. 6. We notice that solving (41) is much faster than solving (28), 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 [START_REF] Borzi | Multigrid methods for pde optimization[END_REF] when the measurements are drawn with p or with q than with u. The reason is probably a better initialisation of (28) 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 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 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.

Discussion 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 [START_REF] Tremblay | Compressive spectral clustering[END_REF]. 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 [START_REF] Borzi | Multigrid methods for pde optimization[END_REF]. To accelerate the method, one could 1) pre-group similar nodes to form N groups such that k N 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 (41) 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 λ λ 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 p * and q * . Appendix A -Proof of the Theorem 2.2

As done in [START_REF] Puy | Random sampling of bandlimited signals on graphs[END_REF], the proof is obtained by applying the following lemma obtained by Tropp in [START_REF] Tropp | User-friendly tail bounds for sums of random matrices[END_REF].

Lemma A.1 (Theorem 1.1, [START_REF] Tropp | User-friendly tail bounds for sums of random matrices[END_REF]). Consider a finite sequence {X j } of independent, random, selfadjoint, positive semi-definite matrices of dimension d × d. Assume that each random matrix satisfies λ max (X j ) R almost surely. Define The matrix 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 E X j and λ max (X j ).

The expected value of each X j is

E X j = E 1 s U k N (ωj ) P (ωj ) P (ωj ) N (ωj ) U k = 1 s U k N =1 p N ( ) P ( ) P ( ) N ( ) U k = 1 s U k N =1 N ( ) N ( ) U k = 1 s U k U k = 1 s I. (65)
Therefore, λ min j E X j = 1 and λ max j E X j = 1. Furthermore, for all j = 1, . . . , s, we have λ max (X j ) = X j 2 max Therefore, for any δ ∈ (0, 1), we have, with probability at least 1ξ, 1δ λ min (X) and λ max (X) 1 + δ (69) provided that (17) holds. Noticing that (69) implies that (1δ) α

Figure 1 .

 1 Figure 1. One image (left) and its partition into superpixels with SLIC technique [23] (right).

  , we associate a matrix N ( ) ∈ R |N |×n that restricts a graph signal to the nodes appearing in N , i.e., N ( )

Figure 3 .

 3 Figure 3. First column from the left: Minnesota graph (top); bunny graph (bottom). The groups N1, . . . , NN are indicated by different colors.Other columns: probability that δ k is less that 0.995 as a function of s. The curves in black, red, blue and green are obtained using the sampling distributions u, p * , q and p, respectively. The top row shows the result for the Minnesota graph. The bottom row shows the result for the bunny graph. The bandlimit k is indicated on top of each curve.

Figure 4 .

 4 Figure 4. Example of sampling distributions. Top panels: p * (left), p (middle), and q (right) for the Minnesota graph at k = 10. Bottom panels: p * (left), p (middle), and q (right) for the bunny graph at k = 25.

Figure 6 .

 6 Figure 6. Left: the curves represent the mean reconstruction snr as a function of the number of sampled superpixels s. The dash-dotted curves are obtained by solving[START_REF] Borzi | Multigrid methods for pde optimization[END_REF]. The continuous curves are obtained by solving (41). Right: mean computation time in seconds as a function of s. The dash-dotted curves are obtained by solving[START_REF] Borzi | Multigrid methods for pde optimization[END_REF]. The continuous curves are obtained by solving (41) using the result of (28) as initial point. In both graphs, the errorbars are at one standard deviation. The black curves are obtained with the uniform sampling distribution u. The blue curves are obtained with q. The red curves are obtained with p.

  use the facts that, for all δ ∈ [0, 1], e -δ (1δ) Proof of Theorem 2.2. We start by noticing that1 s U k M PPMU k = 1 s s j=1U k N (ωj ) P (ωj ) P (ωj ) N (ω l ) U k . (ωj ) P (ωj ) P (ωj ) N (ωj ) U k and X :=

  1 yields, for any δ ∈ (0, 1), P {λ min (X)(1 -

This property is a consequence of (65), proved in Appendix A.

http://www.ntu.edu.sg/home/asjfcai/Benchmark_Website/benchmark_index.html.

F evaluated with the Frobenius norm (right). Middle row: segmentation map estimated from s = 50 sampled superpixels obtained by solving (41) and using a uniform sampling (left), p (middle), q (right). Bottom row: segmentation map estimated from s = 50 sampled superpixels obtained by solving (28) and using a uniform sampling (left), p (middle), q (right).

We emulate user-interactions as follows. For each chosen superpixel, we compute the mean of the ground truth map 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 ỹ ∈ R m from which we want to reconstruct x. We repeat this procedure for s ∈ {50, 70, . . . , 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 x, we first use the fast reconstruction method (41) and then refine the solution at the pixel level with (28), using the solution of the first minimisation problem as initialisation to solve the second minimisation problem. We choose g(L) = L and solve both problems in the limit where γ → 0. In this limit, the problems (41) and (28) become min

respectively. Both problems are solved using FISTA [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF]. The same stopping criteria are used for all experiments.

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

Indeed, for any z ∈ R n , the j th entry of P MAz is

Then, the j th entry of APMz is the scaled sum of the values in the j th sampled group appearing in PMz, which is p -1/2 ωj N (ωj ) z. From the definition of A, the sum is scaled by N ωj -1/2 . Therefore, ( P MAz) j = ( APMz) j for all j ∈ {1, . . . , s}, which terminates the proof. The second property is

The third property is

for all z ∈ R N . To prove this property, we remark that PMA z = A P M z. Indeed, the entries in PMA z corresponding to the first selected group N ω1 are all equal to (p ω1 |N ω1 |) -1/2 zω1 . There are |N ω1 | such entries. One can also notice that the first |N ω1 | entries of A P M z are also equal to (p ωj |N ω1 |) -1/2 zω1 . Repeating this reasoning for all the sampled groups proves the equality. On the one side, we thus have PMA z 2 = A P M z 2 = P M z 2 , where we used (73). On the other side, we have APMA z 2 = P MAA z 2 = P M z 2 , where we used (71) and [START_REF] Perraudin | Gspbox: A toolbox for signal processing on graphs[END_REF]. This terminates the proof.

Proof of Theorem 3.1. As x * is a minimiser of (41), we have

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

[Bounding the terms in (75) involving L]. We define the matrices

By definition of α * and β * , A x * = α * + β * with α * ∈ span(U k ) and β * ∈ span( Ūk ). We recall that L = (AU) g(Λ) (AU) . Therefore, we obtain

In the first step, we used the facts that Ū k α * = 0 and U k β * = 0. The second step follows form the fact that Ū k β * 2 = β * 2 . We also have

The second inequality follows from the facts that U k 2 = 1 and x = Ax. To obtain the third inequality, we used A 2 = 1 and the triangle inequality. For the last step, notice that Ū k x 2 = 0 (as x ∈ span(U k )), Ūk 2 = 1 and use [START_REF] Napoli | Efficient estimation of eigenvalue counts in an interval[END_REF].

[Bounding the terms in (75) involving M]. By definition of ỹ, it is immediate that

For the other term involving M, the triangle inequality yields

Then, we have

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

We continue by lower bounding ÃPMA x * -ÃPMA Ax 2 and upper bounding ÃPMA Ax-ÃPMx 2 separately.

[Lower bound on ÃPMA x * -ÃPMA Ax 2 ]. Equality (74) yields

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

The restricted isometry property [START_REF] Chen | Signal recovery on graphs: Random versus experimentally designed sampling[END_REF] then yields

We used the equality A x * = α * + β * . We thus proved that

We used the facts that à 2 = 1 (see (73)) and that x ∈ A in the last inequality.

Using the inequalities (88) and ( 89) in (84), we arrive at We only rearranged the term in the last step. This proves (44) and terminates the proof.


Colin Laville

Yannick Tillier

Yannick Tillier Mechanical

may come   MECHANICAL CHARACTERISATION OF AORTIC VALVE TISSUES

Introduction

Because of their high prevalence in degenerative valve diseases, aortic valves (AVs) have been widely studied during the last decades. In particular finite element analysis is becoming increasingly popular to support the design of enhanced prostheses with mechanical properties as close as possible to those of natural tissues. The mechanical behaviors of the AV leaflets under a wide range of physiological biaxial loading conditions have been already investigated in the literature [1]. AVs have an oriented network of collagen fibers embedded in an elastin matrix, which is responsible for their hyperelastic and anisotropic behavior. Non-linear transverse isotropic constitutive equations are often used for their modeling, as far as a macroscopically preferred fiber direction can be identified. Dozens of mechanical tests with various loading conditions must generally be carried out to identify those parameters. In this study, a method is proposed in order to simultaneously estimate relevant material and structural characteristics of AV while reducing the number of experiments. An inverse analysis procedure based on the finite element computation of biaxial tensile tests coupled with full-field measurements was used to set up this reduced protocol. In a second time, this protocol was applied to porcine aortic tissues to identify real parameters.

Methods

Sensitivity analysis: AV leaflet tissue behavior was modeled using a modified hyperelastic transversely isotropic HGO [2] material law. This model was implemented in the finite element code FORGE® using the strain energy function (1):

[ ( ) ( ) ] ( ) (1)
where are material parameters, and respectively the first invariant and fourth pseudoinvariant of the right Cauchy-Green tensor, the volume ratio and a Lagrange multiplier. Biaxial tensile tests with several loading conditions were modeled using an arbitrary set of parameters. Then a sensitivity analysis was performed on the set of curves obtained for each loading condition in order to determine an optimal number of experimental conditions (i.e. by removing unnecessary experiments) that enable to identify accurately the initial set of parameters used for creating the reference curves. Experimental measurement: Fresh native AV valves collected from porcine hearts were used. Rectangular samples excised from the center of the leaflets (in the "belly" area) were submitted to optimal biaxial loading conditions that were numerically identified (Fig. 1). The biaxial tensile tests were coupled with full-field measurements in order to get local strain measurements to be compared to numerical sensors.

Results

The sensitivity analysis carried out on the artificial set of curves showed that both local strain measurements and force/displacement curves obtained with at least two different loading conditions are necessary if one wants to accurately determine structural and material parameters. This conclusion was assessed and confirmed thanks to the biaxial tensile tests coupled with digital image correlation (DIC) data that were performed on AV leaflets.

Discussion

In this study we considered heart valve samples as monolayer homogenous materials with only one family of fibers. The same fiber orientation is assigned to each element. Second-harmonic imaging microscopy, using a confocal microscope, is planned to assess this assumption and to estimate the error that was made.

Figure 1 :

 1 Figure 1: the "belly" area on a AV leaflet (left) and the biaxial tensile device (right)

Acknowledgements

The authors thank Carnot M.I.N.E.S for support through grant 20531.