Fabienne Castell
email: fabienne.castell@univ-amu.fr

Luca Avena
email: l.avena@math.leidenuniv.nl

Alexandre Gaudilliere
email: alexandre.gaudilliere@math.cnrs.fr

Clothilde Melot
email: clothilde.melot@univ-amu.fr

AND Alexandre Gaudillière

Clothilde Mélot

Intertwining wavelets or Multiresolution analysis on graphs through random forests

Keywords: 2010 Mathematics Subject Classification. 94A12, 05C81, 05C85, 15A15, 60J20, 60J28 Graph signal processing, multiresolution analysis, wavelet basis, intertwining, Markov processes, random spanning forests, determinantal processes

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Graphs provide a flexible representation of geometric structures of irregular domains, and they are now a commonly used tool in numerous applications including neurosciences, social sciences, biology, transport, communications, etc. One edge between two vertices models an interaction between them, and one can also associate a weight to each edge to quantify levels of interaction, leading to a weighted graph. A graph signal associates data to each vertex of the graph. As a measure of the brain activity in distinct functional regions, functional magnetic resonance images are an instance of a such a graph signal. Therefore, signal processing on weighted graphs has raised significant interest in recent years (see [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 a review on the subject). This paper is concerned with the problem of defining a multiresolution analysis for graph signals. To introduce the problem, let us recall the now classical method to perform such a multiresolution analysis on a regular domain.

1.1. Multiresolution analysis on regular grids. Let us consider a discrete periodic function f : Z n = Z/nZ → R, viewed as a vector in R n . The multiresolution analysis of f is based on wavelet analysis and the so-called multiresolution scheme. Performing the wavelet analysis amounts to compute an "approximation" f 1 ∈ R n/2 and a "detail" component g 1 ∈ R n/2 through classical operations in signal processing such as "filtering" and "downsampling". Roughly speaking one wishes the approximation to give the main trends present in f whereas the detail would contain more refined information. This is done by splitting the frequency content of f into two components: f 1 focuses on the low frequency part of f whereas the high frequencies in f are contained in g 1 . In signal processing and image processing, "filtering" (i.e computing the convolution f k for some well chosen kernel k) allows to perform such frequency splittings. In our case K l (f ) = f k l yields K l (f ) as a low frequency component of f and K h (f ) = f k h yields K h (f ) as a high frequency version of f . The vectors f 1 and g 1 are "downsampled" versions of K l (f ) and K h (f ) by a factor of 2, which means that one keeps one coordinate of K l (f ) and K h (f ) out of two, to build f 1 and g 1 respectively. Thus the total length of the concatenation of the two vectors [f 1 , g 1 ] is exactly n, the length of f . To sum up we have

f 1 (x) = K l (f )(x) = ϕ x , f g 1 (x) = K h (f )(x) = ψ x , f
where • x belongs to the set of downsamples Z n isomorphic to Z/ n 2 Z. • {ϕ x , x ∈ Z n } is the set of functions such that the equality between linear forms ϕ x , • = K l (•)(x) holds for all x ∈ Z n . • In the same way {ψ x , x ∈ Z n } is such that ψ x , • = K h (•)(x) holds for all x ∈ Z n .

The choice of k l and k h is clearly crucial and is done in such a way that perfect reconstruction of f is possible from f 1 and g 1 , so that there is no lost information in the representation [f 1 , g 1 ]. One can iterate the scheme to produce a sequence f N 0 ∈ R n/2 N 0 , g N 0 ∈ R n/2 N 0 , ..., g 1 ∈ R n/2 such that the total length of the concatenated vectors [f N 0 , g N 0 , ..., g 1 ] is exactly n. This is the reason why this scheme yields a multiresolution representation of f . Remark that the perfect reconstruction condition amounts to have B = {ϕ x , ψ x , x ∈ Z n } a basis for the signals f on Z n . A famous construction by Ingrid Daubechies [START_REF] Daubechies | Ten lectures on wavelets[END_REF] derives several families of orthonormal compactly supported such basis B. These families combine localization in space around the point x and localization properties in frequency due to the filtering step they have been built from. Using this space-frequency localization one can derive key properties of the wavelet analysis of a signal which rely on the deep links between the local regularity properties of f and the behavior and decay properties of detail coefficients. Let us finish by saying that wavelet type algorithms led to developments of many powerful algorithms of compression, signal restauration, classification, and other kind of time-frequency transforms etc. We will not go into more detail and refer the interested reader to one of the numerous books on wavelet methods and their applications such as [START_REF] Daubechies | Ten lectures on wavelets[END_REF] or [START_REF] Mallat | A wavelet tour of signal processing[END_REF]. 1.2. Signal processing on graphs. Our aim is to define a multiresolution analysis on a generic weighted finite graph G = (X , w, µ), where:

• X is the set of vertices of G with cardinality |X | = n; • w : X × X → R + is a weight function, which associates to each pair (x, y) a positive weight w(x, y); • µ is a measure on X such that w is reversible w.r.t. µ, meaning that ∀x, y ∈ X , µ(x)w(x, y) = µ(y)w(y, x) .

(

The weight function w provides a graph structure on X , (x, y) being an edge of the graph iff w(x, y) > 0. On such a graph, the operations of translation and downsampling are no longer canonically defined and several attempts to tackle these issues and generalize the wavelet constructions have been proposed: see [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 a review on this subject and [START_REF] Hammond | Wavelets on graphs via spectral graph theory[END_REF] for one of the most popular methods. We will describe some of them in section 1.5, but let us stress that when going to the case of a generic weighted graph, there are three main problems one has to face:

(Q1): What kind of downsampling should we use? What is the meaning of "keep one point out of two"? (Q2): On which weighted graph should the approximation f 1 be defined to iterate the procedure? (Q3): Which kind of filters should one use? What is a good "local" mean?

This paper proposes a construction of a multiresolution analysis on a weighted graph which is based on a random downsampling method to answer question (Q1), and on Markov processes interwining to answer questions (Q2) and (Q3). Markov processes enter naturally into the game through the graph Laplacian L defined on 2 (X , µ) by

Lf (x) := y∈X w(x, y)(f (y) -f (x)), (2)
where f : X → R is an arbitrary function, and w(x, y) is now viewed as the transition rate from x to y. For x ∈ X , let w(x) := y∈X \{x} w(x, y) .

Note that L acts on functions as the matrix, still denoted by L:

L(x, y) = w(x, y) for x = y ; L(x, x) = -w(x) .

L is the generator of a continuous time Markov process X = (X(t), t ≥ 0), which jumps from x to a neighboring site y at random exponential times (see section 2.3). The law of X starting from x is denoted by P x , and E x is the expectation w.r.t. P x . Apart from the symmetry assumption (1), we will assume throughout the paper that

L is irreducible. (3)
Under this assumption, µ is the unique (up to multiplicative constants) invariant measure of the process (X(t), t ≥ 0) (i.e µL = 0), and µ(x) > 0 for any x ∈ X . Since X is finite, we can assume without loss of generality that µ is a probability measure on X . 1.3. Our approach for downsampling, weighting and filtering procedures on graphs. The first step in our proposal of a multiresolution scheme consists in answering (Q1) and constructing a random subset X of X , whose main feature is to be "well spread" on X . In this respect, we use an adaptation of Wilson's algorithm ( [START_REF] Wilson | Generating random spanning trees more quickly than the cover time[END_REF]) studied in [START_REF] Avena | Two applications of random spanning forests[END_REF], whose output is a random rooted spanning forest Φ on G. Its law π q depends on the weight function w and on a positive real q. This algorithm and the properties of the corresponding random forest will be described in section 3.3. The set of roots of Φ, denoted by ρ(Φ), has the nice property of being a determinantal process on X with kernel given by the Green kernel:

K q (x, y) := q(qId -L) -1 (x, y) = P x [X(T q ) = y] , (4)
where T q is an exponential random variable with parameter q, independent of the process X. Being a determinantal process, the roots of Φ tend to repulse each other, and this repulsiveness property results in the noteworthy fact that the mean time needed by the process (X(t), t ≥ 0) to reach ρ(Φ), does not depend on the starting point in X \ ρ(Φ). Therefore, this set of roots is a natural choice for X , especially as the parameter q can be tuned to control the mean number of points in ρ(Φ). See Section 3 for more details.

Once the downsampling X is fixed, it remains to answer questions (Q2) and (Q3), i.e. we have to construct a weight function w on X × X , and to define a filtering method which gives the approximation f 1 and the detail g 1 of a given function f 0 in 2 (X , µ). This will be achieved by looking for a solution ( L, Λ) to the intertwining relation

LΛ = ΛL , (5)
where • L is a Markov generator on X ;

• Λ : X × X → R + is a positive rectangular matrix. An intertwining relation gives a natural link between two Markov processes living on different state spaces. It appeared in the context of diffusion processes in the paper by Rogers and Pitman [START_REF] Rogers | Markov functions[END_REF], as a tool to state identities in laws, and was later successfully applied to many other examples (see for instance [START_REF] Carmona | Beta-gamma random variables and intertwining relations between certain Markov processes[END_REF], [START_REF] Matsumoto | An analogue of Pitman's 2M -X theorem for exponential Wiener functionals. I. A time-inversion approach[END_REF], etc). In the context of Markov chains, intertwining was used by Diaconis and Fill [START_REF] Diaconis | Strong stationary times via a new form of duality[END_REF] to construct strong stationary times, and to control the rate of convergence to equilibrium.

At the time being, applications of intertwining include random matrices [START_REF] Donati-Martin | Some properties of the Wishart processes and a matrix extension of the Hartman-Watson laws[END_REF], particle systems [START_REF] Warren | Dyson's Brownian motions, intertwining and interlacing[END_REF], etc. For our purpose, the intertwining relation can be very useful since

• L provides a natural choice for the graph structure to be put on X ; • each row ν x := Λ(x, •) (x ∈ X ) defines a positive measure on X , which can serve as a "local mean" around x. f 0 will be approximated by the function f 1 defined on X by

∀x ∈ X , f 1 (x) := ν x(f 0 ) = x∈X Λ(x, x)f 0 (x) ;
• it will serve as a basis for getting a Jackson-type inequality (see Proposition 18).

To be of any use in signal processing, the "filters" (ν x, x ∈ X ) have to be well localized in space and frequency. In the graph context, frequency localization means that the filters belong to an eigenspace of the graph laplacian L. Hence, we are interested in solutions to [START_REF] Cohen | Numerical analysis of wavelet methods[END_REF] such that the measures (ν x, x ∈ X ) are linearly independent measures tending to be non-overlapping (space localization), and contained in eigenspaces of L. In addition, in order to iterate the procedure, we also need L to be reversible on X . Note that saying that (Λ, L) is an exact solution to [START_REF] Cohen | Numerical analysis of wavelet methods[END_REF] implies that the linear space spanned by the measures (ν x, x ∈ X ) is stable by L, and is therefore a direct sum of eigenspaces of L, so that these measures provide filters which are frequency localized. Hence the error in the intertwining relation is a measure of frequency localization: the smaller the intertwining error, the better the frequency localization. Finding solutions to [START_REF] Cohen | Numerical analysis of wavelet methods[END_REF] is the purpose of [START_REF] Avena | Approximate and exact solutions of intertwining equations through random forests[END_REF], where an exact linearly independent solution to (5) is provided. However this solution tends to be overlapping. [START_REF] Avena | Approximate and exact solutions of intertwining equations through random forests[END_REF] provides also approximate solutions to [START_REF] Cohen | Numerical analysis of wavelet methods[END_REF] with small overlapping, and it is one of these approximate solutions we use in our multiresolution analysis scheme. In order to describe it, we first consider the trace process of (X(t), t ≥ 0) on X , i.e. the process (X(t), t ≥ 0) sampled at the passage times on X . This process is a Markov process on X , whose generator L defines a weight function w on X , symmetric with respect to the measure µ restricted on X . We will recall some general facts about Markov processes, and make the previous statements more precise in section 2.3.

Turning to the rectangular matrix Λ, it depends on a parameter q and involves the Green kernel K q (see definition [START_REF] Coifman | Diffusion wavelets[END_REF]). The rectangular matrix Λ is just the restriction of K q to X × X . Note that • when q goes to 0, for any x ∈ X , K q (x, y) goes to µ(y) so that ( 5) is clearly satisfied.

µ being the left-eigenvector of L corresponding to the eigenvalue 0, the K q (x, y) are well frequency localized. However, the vectors (K q (x, •), x ∈ X ) become linearly dependent and very badly space localized. • when q goes to ∞, K q (x, •) goes to δ x. Hence, the space localization is perfect. However, the frequency localization is lost, and the error in (5) tends to grow.

Hence, a compromise has to be made concerning the choice of q , and we will discuss this point later.

To sum up our proposal, two parameters q and q being fixed,

• take X := ρ(Φ), Φ being sampled with π q . Set X = X \ X .

• L is the generator of the trace process of X on X . It can be shown (see Lemma 8 section 4.1) that L is irreducible and reversible w.r.t. the probability measure µ X , which is the measure µ conditioned to the set X : µ X (A) = µ(A ∩ X )/µ( X ). L is actually the Schur complement in L, of L restricted to X , which was already used for instance in [START_REF] Hammond | Wavelets on graphs via spectral graph theory[END_REF] . For any x, ȳ in X , such that x = ȳ, the weight function w on X × X is then defined by w(x, ȳ) = L(x, ȳ).

• Λ(x, •) := K q (x, •).
f (=f 0 ) being a function in 2 (X , µ), we define the approximation of f as the function f (=f 1 ) defined on X by

∀x ∈ X , f (x) = Λf (x) = K q f (x) , (6)
and its detail function as the function f (=g 1 ) defined on X by

∀x ∈ X , f (x) = (K q -Id)f (x) . (7)
We proved in [START_REF] Avena | Approximate and exact solutions of intertwining equations through random forests[END_REF] that some localization property of the (ν x = Λ(x, •), x ∈ X ) defined by ( 6) followed from the fact that X is a determinantal process with kernel K q . Since X = X \ X is also a determinantal process, with kernel Id -K q , this suggested the detail definition [START_REF] Diaconis | Strong stationary times via a new form of duality[END_REF], in which the sign convention is chosen to have a self-adjoint analysis operator

U : f → ( f , f ) in 2 (X , µ).
The process can then be iterated with ( X , L, µ X ) in place of (X , L, µ). This leads to a multiresolution scheme.

1.4. Description of the results. We give now a description of our main results. Note that whatever the choice of a subset X of X , one can still define L, Λ, f and f as we just did it previously. Our first set of results assumes that X is any subset of X , and provides a reconstruction formula, and bounds on various operators including approximation, detail and analysis operators. They are expressed in terms of the hitting time of X :

H X := inf t ≥ 0, X(t) ∈ X .
To state them, we introduce some notation. First of all, we define the maximal rate α > 0 by

α = max x∈X w(x) . (8)
Hence, the matrix P = Id + L/α is a stochastic matrix on X . X being any subset of X , we define the two positive real numbers β and γ by

1 γ := max x∈ X E x [H X ] , 1 β := max x∈ X z∈X P (x, z)E z [H X ] . (9)
ᾱ denotes the maximal rate of L:

ᾱ := max x∈ X -L(x, x) .
We stress the obvious fact that γ, β and ᾱ are functions of the subset X . Finally, when (M (u, v), u ∈ U, v ∈ V) is a rectangular matrix, when U is a subset of U, and V a subset of V, M U V denotes the submatrix of M obtained by restricting the entries of M to U and V .

To begin with our results, we get a reconstruction formula:

Reconstruction formula (Proposition 10.) For any f ∈ 2 (X , µ), let f ∈ 2 ( X , µ) and f ∈ 2 ( X , µ) be defined by (6) and [START_REF] Diaconis | Strong stationary times via a new form of duality[END_REF]. Then,

f = R f + Ȓ f ,
where R are Ȓ are rectangular matrices indexed respectively by X × X and X × X , whose block decompositions are

R = Id X -1 q L (-L X X ) -1 L X X , and Ȓ = L X X (-L X X ) -1 -Id X -q (-L X X ) -1 . ( 10
)
We provide also upper bounds on norms of various operators: the approximation operator R, the detail operator Ȓ, and the interwining error operator LΛ -ΛL. These bounds are given here when these operators are seen as operators from one ∞ space to another one. They are stated in greater generality in Propositions 14, 11 and 16.

Approximation and detail operator norms (Propositions 14 and 11).

Let X be any proper subset of X (i.e. ∅ X X ), X = X \ X , and let R and Ȓ be the operators defined in [START_REF] Gavish | Multiscale wavelets on trees, graphs and high dimensional data: theory and applications to semi supervised learning[END_REF]. For any f ∈ ∞ ( X , µ X ) and any

f ∈ ∞ ( X , µ X ), R f ∞,X ≤ 1 + 2 ᾱ q f ∞, X , and Ȓ f ∞,X ≤ max α β ; 1 + q γ f ∞, X . (11)
Intertwining error norm (Proposition 16).

Let X be any proper subset of X , X = X \ X . For any f ∈ ∞ (X , µ),

LΛ -ΛL f ∞, X ≤ 2q α β f ∞,X . (12)
Note that these bounds reflect the competition between LΛ -ΛL f ∞, X (frequency localization) and R f ∞,X (space localization). Actually, the term ᾱ/q appearing in ( 11) is a decreasing function of q and an increasing function of X , while the term q /β in ( 12) is increasing in q and decreasing in X .

Other results are stated in the paper, including a bound on the norm of the detail f in terms of the norm of Lf (Proposition 15), and a Jackson's type inequality on the approximation error after K steps of the mutiresolution scheme (Proposition 18).

The other set of results focuses on the case where X is the set of roots of the random forest Φ. They provide estimates on its cardinality X , and on the quantities ᾱ, β and γ involved in the previous statements. X being random, all these quantities are random ones, and the estimates we get are averaged ones. To state them, we will assume that the random variable Φ is defined on some probability space (Ω f , A f , P q ). The corresponding expectation will be denoted by E q .

Cardinality estimates (Proposition 6). Let X = ρ(Φ) and X = X \ X . Let r ∈]0; 1[, and define R r = {x ∈ X ; w(x) ≥ rα} the set of "rapid" points. Then,

|X | q q + α ≤ E q X ; E q X ∩ R r ≤ |R r | q + (2 -r)α q + 2α .
These estimates can be used to tune the parameter q in order to target a given proportion of kept points. Next, we obtain bounds on the quantities ᾱ, β and γ when X = ρ(Φ). Unfortunately, these bounds are just lower bounds, and are moreover averaged ones. Nevertheless, since they are expressed in terms of the cardinality of ρ(Φ), they can easily be estimated by Monte Carlo methods, and can serve as a guide for the choice of q and q . These are then the key estimates where our particular random choice for X plays a central role.

Estimates on ᾱ, β, γ (Propositions 20, 21 and 22). Let X = ρ(Φ) and X = X \ X .

E q [ᾱ] ≥ qE q   X X + 1   , E q 1 β ≥ E q   X α X   , E q 1 γ ≥ 1 q E q   |X | X + 1   . ( 13
)
Discussion on the choice of q and q Based on these estimates, we argue in section 6.4 that the parameter q should be chosen in [θ 1 α; θ 2 α] in order to ensure that the mean number of vertices of the subgraph is at least a given proportion (θ 1 /(θ 1 + 1)) of the size of the original graph, and that a given proportion of the rapid points are decimated. In addition it should minimize the function q → E q (ᾱ/β), for the approximation operator norm and the intertwining error to be small (see [START_REF] Hammond | Wavelets on graphs via spectral graph theory[END_REF] and ( 12)). The Monte-Carlo estimation of this function is computationally costly, so that we propose to minimize the function q ∈ [θ 1 α; θ 2 α] → qE q X /(1 + X ) E q X / X . This is possible in practice since, from [START_REF] Avena | Two applications of random spanning forests[END_REF], we can simultaneously sample a whole continuum of forests Φ q for q ∈ [θ 1 α; θ 2 α], each of them with the correct distribution π q . Once q has been fixed, X is sampled according to π q , ᾱ and β are computed and q is chosen equal to 2ᾱ X / X , which will ensure numerical stability of the algorithm (see Section 6.4).

1.5. Related works. Previous authors have explored multiresolution analysis on graphs, and have proposed answers to questions (Q1-3). Far from being exhaustive, we describe some of these, and refer the interested reader to [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 a more complete state of the art.

Concerning the downsampling procedure, also referred to as the graph coarsening problem, many approaches have been investigated. To mention a few, one can try to decompose the graph into bipartite graphs for which the notion of "one every two points" is clear. This is the way followed in [START_REF] Narang | Perfect reconstruction two-channel wavelet filterbanks for graph structured data[END_REF][START_REF] Nguyen | Downsampling of signal on graphs via maximum spanning tree[END_REF], using either coloring-based downsampling, or maximum spanning tree. In [START_REF] Shuman | A multiscale pyramid transform for graph signals[END_REF], the partitioning of X into two subsets, is based on the sign of the eigenfunction of -L associated to the maximal eigenvalue. Finally, the authors of [START_REF] Tremblay | Subgraph-based filterbanks for graph signals[END_REF] use a community detection algorithm maximizing the modularity, to partition the original graph into many connected subgraphs of small size. The vertices of the downsampled graph are the elements of the partition, and are not properly speaking selected points of the original graph.

Turning to the weighting procedure, bipartite graphs designed methods put an edge between two selected nodes, if they share at least a neighbor in the original graph, and the weight may be proportional to the number of shared neighbors. Starting from a bipartite graph, this leads immediately to a non bipartite subgraph, thus the need to "decompose" a graph by bipartite ones. The authors of [START_REF] Shuman | A multiscale pyramid transform for graph signals[END_REF] rely on the so-called Kron reduction, which is the same as computing the Schur complement. In [START_REF] Tremblay | Subgraph-based filterbanks for graph signals[END_REF], the weight between two communities is the sum of the weight of edges linking these communities.

Finally, as far as the filtering procedure is concerned, various filter banks have been proposed, leading sometimes to orthonormal basis, or just frames. The graph wavelet filtering bank of [START_REF] Narang | Compact support biorthogonal wavelet filterbanks for arbitrary undirected graphs[END_REF] is more specifically designed for bipartite graphs, and exploit the specific features of bipartite graph Laplacian. It produces an orthonormal basis. The diffusion wavelets of [START_REF] Coifman | Diffusion wavelets[END_REF] are obtained by constructing orthonormal basis of the spaces

(V k-1 V k , k = 1, • • • , K), where V k = ((-L) 2 k
), and are thus designed to form an orthonormal basis. The spectral graph wavelets of [START_REF] Hammond | Wavelets on graphs via spectral graph theory[END_REF] are functions of type g(-tL)δ n where g is a function localized around 0 for the low-pass filters, or away from 0 for the high-pass ones. Different values of t are used to select different frequencies, thus leading to a frame. Since the computation of g(-tL) requires the knowledge of the spectral decomposition of L, polynomial approximations are performed. In [START_REF] Tremblay | Subgraph-based filterbanks for graph signals[END_REF], the filters used are the eigenfunctions of the Laplacians restricted to each small community, naturally leading to an orthonormal basis. In [START_REF] Gavish | Multiscale wavelets on trees, graphs and high dimensional data: theory and applications to semi supervised learning[END_REF], assuming that we are given a subsampling procedure encoding the geometry of the original graph, the authors construct a Haar basis adapted to this subsampling.

Compared to these works, our downsampling approach through Wilson's algorithm, is a partitioning method in the spirit of [START_REF] Tremblay | Subgraph-based filterbanks for graph signals[END_REF], a community corresponding to a tree of the forest. The mean number of selected points can be adjusted through the parameter q.

The filters we use are a special case of spectral graph wavelets of [START_REF] Hammond | Wavelets on graphs via spectral graph theory[END_REF], the function g being equal to 1/(1 + x) for the scaling function, and x/(1 + x) for the wavelet, and t corresponding to 1/q (note that K q = q (q Id -L) -1 = (Id -L/q ) -1 ). This choice is very natural in the signal processing variational approach, since K q f minimizes the function g → g, -Lg + q g -f 2 . Since we use only one value of t at each step, the filter bank we obtain is a basis, which is not orthogonal, but tends to be so when q 1 (in this case K q (x, •) δ x, while (K q -Id)(x, •) L(x, •)/q ). The weighting procedure through Kron reduction has also already be used in [START_REF] Shuman | A multiscale pyramid transform for graph signals[END_REF]. But in our approach, weighting and filtering are linked together through the intertwining relation.

To sum up, the advantages of our approach are:

• to provide a new partitioning method, allowing to tune the mean number of kept vertices;

• to link the subgraph structure and the choice of the filters through the intertwining relation;

• to produce a filter bank leading to a basis;

• to mimic the steps of a classical wavelet analysis algorithm;

• to allow the computation of various error bounds;

• to suggest a systematic method to choose the parameters q and q , and ensure numerical stability; • to have a computational complexity similar to already existing methods. Actually, apart from the sampling of the random forest which is of order α |X | /q, our approach only requires the computation of K q and of the Schur complement, already present in [START_REF] Hammond | Wavelets on graphs via spectral graph theory[END_REF][START_REF] Shuman | A multiscale pyramid transform for graph signals[END_REF]. Starting from a sparse matrix L, K q can be efficiently approximated by polynomials of small order. Due to the good repartition property of our random X , such a polynomial approximation can also be implemented for the inversion involved in the computation of L. This results in a sparse L, however the level of sparsity obtained is not enough to go on with the algorithm for large graphs. Hence, as in [START_REF] Shuman | A multiscale pyramid transform for graph signals[END_REF], a sparsification step can be added after the computation of L. Our proposition for this sparsification will be guided by the intertwining error bounds that provide a Jackson-like inequality (see Section 7).

Finally, we would like to mention that our multiresolution scheme construct a basis of the space of signals, which can be used to analyze, compress, etc, any signal on the graph. When one wants to handle just one specific signal, adaptative multiresolution schemes as in [START_REF] Elisha | Wavelet decompositions of random forests. Smoothness analysis, sparse approximation and applications[END_REF] may be more appropriate.

1.6. Organization of the paper. We begin in section 2 with notations used throughout the paper. Section 3 is devoted to the description and the properties of the random forest, and to the downsampling procedure. The weighting and the filtering procedures, X being any subset of X , are discussed in section 4. We prove in this section the bounds on operator norms. We go on with the iterative scheme and Jackson's inequality in section 5.2. The discussion on the choice of q and q and the estimates on ᾱ, β and γ, are given in section 6. We summarize the pyramidal algorithm and discuss computational issues in section 7, and we finally end with numerical experiments illustrating the method in section 8.

Notations and preliminary results.

In this section, we give the notations used throughout the paper, and we state some useful results.

Sets of functions and measures.

A function f on X will be seen as a column vector, whereas a signed measure on X will be seen as a row vector. For p ≥ 1, p (X , µ) is the space of functions on X endowed with the norm

f p = x∈X |f (x)| p µ(x) 1/p .
The scalar product of two functions f and g in 2 (X , µ) is

f, g = x∈X f (x)g(x)µ(x)
When f is a function on X , f will denote the signed measure whose density w.r.t µ is f :

f (A) = x∈A µ(x)f (x) for all subset A of X .
Similarly, when ν is a signed measure on X , ν is the density of ν w.r.t µ: ν (x) = ν(x)/µ(x) for any x ∈ X . We remind the reader the following standard results concerning the Schur complement (see for instance [START_REF] Zhang | The Schur complement and its applications[END_REF]):

Proposition 1. Assume that D is invertible. (1) M = Id B 0 D S M (D) 0 D -1 C Id . (2) det(M ) = det(D) det(S M (D)).
(3) M is invertible if and only if S M (D) is invertible. In that case,

M -1 = S M (D) -1 -S M (D) -1 BD -1 -D -1 CS M (D) -1 D -1 + D -1 CS M (D) -1 BD -1 .
(4) Assume that M defines a positive symmetric operator in 2 (X , µ). Let λ max (M ) (respectively λ min (M )) be the largest (respectively the smallest) eigenvalue of M . Then, S M (D) is also positive symmetric and

λ min (M ) ≤ λ min (S M (D)) ≤ λ max (S M (D)) ≤ λ max (M ) .
2.3. Markov process. Consider an irreducible continuous time Markov process (X(t), t ≥ 0) on X , with generator L given by ( 2) satisfying ( 1) and ( 3). We recall the definition (8) of α > 0 as the maximal rate, and that of the matrix P := L/α + Id. P is an irreducible stochastic matrix, and we denote by ( Xk , k ∈ N) a discrete time Markov chain with transition matrix P . The process (X(t), t ≥ 0) can be constructed from ( Xk , k ∈ N) and an independent Poisson process (τ i , i ≥ 0) on R + with rate α. At each event of the Poisson process, X moves according to the trajectory of X:

X(t) = +∞ i=0 Xi 1 τ i ≤t<τ i+1 .
By ( 1) and (3), µ is strictly positive. In addition, -L defines a positive symmetric operator on 2 (X , µ), and we denote by (λ i ; i = 0, • • • , n -1) the real eigenvalues of -L in increasing order. It follows from the fact that P is irreducible that

0 = λ 0 < λ 1 ≤ λ 2 • • • ≤ λ n-1 ≤ 2α . ( 14
)
The right eigenvector of -L associated to λ i is denoted by e i :

-Le i = λ i e i .
The (e i , i = 0, • • • , n -1) are normalized to have an 2 -norm equal to 1, and form an orthonormal basis of 2 (X , µ). By construction, e 0 = 1. For any subset A of X , H A is the hitting time of A by the process X:

H A := inf {t ≥ 0, X(t) ∈ A} , while H +
A is the return time to A by the process X:

H + A := inf {t ≥ τ 1 , X(t) ∈ A} . (15)
3. Random spanning forests, Wilson's algorithm and downsampling procedure.

We spend now some time on the description and the properties of the spanning random forest Φ used in the downsampling procedure . Let us call E the set of unoriented edges of G, that is the set of pairs {x, y} such that w(x, y) > 0 (and w(y, x) > 0). A spanning unoriented forest is a graph without cycles, with X as set of vertices, and a subset of E as edge set. An unoriented tree is a connected component of such a forest. By choosing in each tree one specific vertex, which we call root, we define a rooted spanning forest. Note that the number of roots is the same as the number of trees. Orienting each edge of a tree toward its root, we obtain an spanning oriented forest (s.o.f) φ. The set of roots of a spanning oriented forest φ is denoted by ρ(φ). If e = (x, y) is an oriented edge, we will use w(e) for w(x, y), and say that e ∈ φ if e is an edge of a tree of the forest φ.

3.1.

A probability measure on forests. We introduce now a real parameter q > 0, and associate to each oriented forest a weight

w q (φ) := q |ρ(φ)| e∈φ w(e) , (16)
where |ρ(φ)| is the cardinality of ρ(φ), i.e. the number of trees in the forest φ. w 1 (φ) will be denoted by w(φ) so that w q (φ) := q |ρ(φ)| w(φ) . These weights can be renormalized to define a probability measure on the set of spanning oriented forest:

π q (φ) := w q (φ) Z(q) , (17)
where the partition function Z(q) is given by

Z(q) := φ s.o.f. w q (φ) . (18)
3.2. Wilson's algorithm. A way to sample a random s.o.f. Φ from π q is given by the following iterative algorithm. Let Φ c be the current state of the forest being constructed, and let V c be the set of vertices of Φ c . At the beginning, V c is equal to ∅. While X \ V c = ∅, perform the following steps:

• Choose a point x at random in X \ V c .

• Let evolve the Markov process (X(t), t ≥ 0) from x, and stop it either when it reaches V c , or after an independent exponential time of parameter q. • Erase the loops of the trajectory drawn by X. We obtain a self-avoiding path C starting from x and oriented towards its end-point. • Add C to Φ c . Each iteration of the "while loop" stopped by the exponential time, gives birth to another tree. Wilson's algorithm is not only a way to sample π q . It provides also a powerful tool to analyse it, the reason being that it does not depend on the way we "choose a point" in the first step of the "while" loop. In addition there exists a coupling of the probability measures π q for different values of q. This means that we can construct the random forests for different values of q from the same set of random variables. This coupling is explained in [START_REF] Avena | Two applications of random spanning forests[END_REF].

3.3.

Properties of the random forest. Using Wilson's algorithm, and the explicit knowledge of the law of a loop-erased Markov process, the following statement, as all the results stated in this section, is proven in [START_REF] Avena | Two applications of random spanning forests[END_REF]: Proposition 2. Partition function. The partition function Z(q) is the characteristic polynomial of L:

Z(q) = det(q Id -L) . ( 19
)
Some other important features of the random s.o.f. Φ are listed below.

Proposition 3. Number of roots. For all k ∈ {0, • • • , n}, P q [|ρ(Φ)| = k] = I⊂{0,••• ,n-1} |I|=k i∈I q q + λ i i / ∈I λ i q + λ i .
Otherwise stated, the number of roots has the same law as n-1 i=0 B i where B 0 , • • • , B n-1 are independent, B i having Bernoulli distribution with parameter q q+λ i . Note that since λ 0 = 0, B 0 = 1 a.s. and we recover the fact that |ρ(Φ)| ≥ 1 a.s.. Proposition 4. Set of roots. ρ(Φ) is a determinantal process on X with kernel K q (see definition in (4)), i.e. for any subset A of X ,

P q (A ⊂ ρ(Φ)) = det A (K q )
, where det A is the minor defined by the rows and columns corresponding to A.

The following statement involves two independent sources of randomness, the Markov process X starting from x defined on (Ω, A, P x ), and the random forest Φ defined on (Ω f , A f , P q ). Integration on the product space (Ω × Ω f , A ⊗ A f ) w.r.t. the product measure P x,q := P x ⊗ P q is denoted by E x,q . Proposition 5. Hitting time of ρ(Φ).

(1) For any x ∈ X ,

E x,q H ρ(Φ) = P q [|ρ(Φ)| > 1] q . (20)
(2) For any x ∈ X , and m ∈ {1, • • • , n},

E x,q H ρ(Φ) | |ρ(Φ)|=m = P q [|ρ(Φ)| = m + 1] qP q [|ρ(Φ)| = m] . (21)
Note that these expressions do not depend on the starting point x, saying that in some sense, the roots of the random forest are "well spread" on X .

3.4. The downsampling procedure. In view of the results of section 3.3, X := ρ(Φ) is a natural candidate as a downsampling of X . We give now estimates on the mean of |ρ(Φ)| which unlike Proposition 3, do not depend on the knowledge of the eigenvalues λ i . To this purpose, we introduce the mean value of the eigenvalues of -L/α:

m(L) := 1 |X | |X |-1 i=0 λ i α = Trace(-L) α |X | .
Note that by [START_REF] Marchal | Loop-erased random walks, spanning trees and Hamiltonian cycles[END_REF], m(L) ∈]0, 2]. We get then

Proposition 6. Let X = ρ(Φ) and X = X \ X . Then (1) E q X ≥ |X | q q+α . (2) E q X ≥ |X | α q+2α m(L). (3) For any r ∈ [0, 1], set R r := {x ∈ X s.t. w(x) ≥ rα}. E q X ∩ R r ≥ |R r | rα q + 2α
.

Remark: Hence, in a loose sense, X contains a given proportion of points in X , and X contains a given proportion of the "rapid" points in X , i.e. points for which the rate of escape is high. Since they do not depend on the spectral decomposition of L, these estimates can be helpful concerning the choice of q. Taking for instance q ∈ [2α/3; α] ensures that the mean proportion of sampled points is greater than 2/5, and that the mean proportion of decimated points is greater than m(L)/3.

The proof of Proposition 6 relies on the following lemma.

Lemma 7. For any x ∈ X , and any q > 0,

K q (x, x) ≥ q q + w(x) . ( 22
)
1 -K q (x, x) ≥ w(x) q + λ n-1 ≥ w(x) q + 2α . ( 23
)
Proof. Remind that (e i , i = 0, • • • n -1) is an orthonormal basis of eigenfunctions of -L in 2 (X , µ), and let (1 x , x ∈ X ) be the canonical basis of R X (1 x (y) = 0 for y = x and 1 x (x) = 1). Note that

w(x) = -L(x, x) = -L1 x (x) = n-1 i=0 λ i e i , 1 x e i (x) = n-1 i=0 λ i e i (x) 2 µ(x) .
In the same way,

K q (x, x) = n-1 i=0 q q + λ i e i (x) 2 µ(x) .
Note that for any x ∈ X ,

n-1 i=0 e i (x) 2 µ(x) = n-1 i=0 1 x , e i 2 1 µ(x) = 1 x 2 µ(x) = 1 ;
Hence, for any x ∈ X , (e

2 i (x)µ(x); i = 0, • • • , n -1) is a probability distribution on {0, 1, • • • , n -1}.
The function λ ∈ R + → q/(q + λ) being convex, it follows from Jensen's inequality, that for any x ∈ X ,

K q (x, x) ≥ q q + n-1 i=0 λ i e i (x) 2 µ(x) = q q + w(x) .
Moreover,

1 -K q (x, x) = n-1 i=0 λ i q + λ i e i (x) 2 µ(x) , ≥ 1 q + λ n-1 n-1 i=0 λ i e i (x) 2 µ(x) , = w(x) q + λ n-1 .
This ends the proof of Lemma 7, since λ n-1 ≤ 2α.

We return now to the proof of Proposition 6.

Proof.

(1) It follows from Proposition 4 that E q X = x∈X P q [x ∈ ρ(Φ)] = x∈X K q (x, x). Using [START_REF] Tremblay | Subgraph-based filterbanks for graph signals[END_REF] and the definition of α, we obtain that

E q X ≥ q q+α |X |. (2) By Proposition 3, E q X = n-1 i=0 λ i q+λ i ≥ 1 q+λ n-1 Trace(-L) .
(3) By ( 23),

E q X ∩ R r = x∈Rr P q [x / ∈ ρ(Φ)] = x∈Rr (1 -K q (x, x)) ≥ x∈Rr w(x) q+λ n-1 ≥ |R r | rα q+λ n-1 .
4. The weighting and filtering procedures.

This section aims to describe the weighting and the filtering procedures we use in the multiresolution scheme, once X has been chosen. Hence, we will assume throughout this section that X is any proper subset of X . 4.1. The weighting procedure. Let us define for any x, ȳ ∈ X ,

P (x, ȳ) := P x X(H + X ) = ȳ . (24)
P is a stochastic matrix on X , and one can associate to P a weight function w and a Markov generator L on X , which are unique up to multiplicative constants. The next lemma explains how to compute L from L.

Lemma 8. Let X be any proper subset of X , and set X = X \ X . Let P be the stochastic matrix on X defined by [START_REF] Wilson | Generating random spanning trees more quickly than the cover time[END_REF]. Then, L X X is invertible. Let L be the Schur complement of L X X in L. L is an irreducible Markov generator on X , P = Id+ L/α, and for any x, ȳ ∈ X , µ(x) L(x, ȳ) = µ(ȳ) L(ȳ, x).

Proof. Let Q(x, x) := P x X(H + X ) = x , so that P = Q X X . Using Markov property at time τ 1 leads to Q(x, x) = z∈X P (x, z)P z [X(H X ) = x] = P (x, x) + z∈ X P (x, z)Q(z, x). Hence, since P = Id + L/α, Q X X = P X X + P X X Q X X Q X X = P X X + P X X Q X X ⇐⇒ P = P X X + P X X Q X X -L X X Q X X = L X X . (25)
We are now going to prove that -L X X is invertible. Let φ X be a vector of Ker(-L X X ), and let ψ := (0 X , φ X ). We get on one hand ψ, -Lψ = 0. On the other hand, ψ, -Lψ = n-1 i=0 λ i ψ, e i 2 .

Therefore, for any i ∈ {0, • • • , n -1}, λ i ψ, e i 2 . λ 0 = 0 being simple, ψ, e i = 0 for all i ≥ 1.

Hence, ψ is colinear to e 0 = 1. If X = ∅ and X = ∅, this implies that φ X ≡ 0. Hence -L X X is invertible. Going back to [START_REF] Zhang | The Schur complement and its applications[END_REF], we obtain that

Q X X = (-L X X ) -1 L X X , (26)
and

P = Id X + 1 α L X X - 1 α L X X (L X X ) -1 L X X = Id X + L α . (27)
Therefore, ∀x, ȳ ∈ X , x = ȳ, L(x, ȳ) = α P (x, ȳ) ≥ 0, and L1 X = 0. L is thus a Markov generator on X . The fact that L is reversible w.r.t. µ is a direct consequence of the fact that L is a Schur complement, and of the reversibility of L w.r.t. µ. Concerning the irreducibility statement, if x and ȳ are two distinct elements of X , the irreducibility of L implies that there exists a path in X going from x to ȳ with positive probability for P . Sampling this path on the passage times on X , we obtain a path in X going from x to ȳ with positive probability for P .

Once L is defined, w and ᾱ are defined in the same way that w and α were defined from L, i.e. ∀x ∈ X , w(x) := -L(x, x) ; ᾱ := max x∈ X w(x) .

(28)

4.2. The filtering procedure. We have now to define the "approximation" and the "detail" of a given function f on X , and to say how to reconstruct f from the approximation and the detail. This is the aim of this section.

4.2.1. The analysis operator. Let q > 0 be fixed and define the collection (ν x, x ∈ X ) of probability measures on X by ν x(x) = K q (x, x) .

Let (ρ x, x ∈ X ) be the collection of signed measures on X defined by

ρ x(x) = (K q -Id)(x, x) . ( 29
)
We associate to the collection of measures (ν x, ρ x) x∈ X ,x∈ X the corresponding collection of functions (φ x, ψ x) x∈ X ,x∈ X :

φ x(x) := ν x(x) = ν x(x) µ(x) ; ψ x(x) := ρ x(x) = ρ x(x) µ(x) . ( 30
)
The functions (ψ x, x ∈ X ) play the role of "wavelets" in our proposition of a multiresolution analysis. Note that the functions φ x are naturally normalized in 1 (X , µ), but this is not the case for the functions ψ x ( ψ x 1 = 2(1 -K q (x, x))). Note also that for any x ∈ X , ψ x, 1 = 0. In addition, we get the following result:

Lemma 9. The family (φ x, ψ x) x∈ X ,x∈ X is a basis of 2 (X , µ).
Proof. Since µ(x) > 0 for all x ∈ X , it is equivalent to prove that the matrix M whose row vectors are given by (ν x, ρ x) x∈ X ,x∈ X is invertible. M can be rewritten in terms of its block matrices according to the rows and columns indexed by X and X :

M = K q | X , X K q | X , X K q | X , X K q | X , X -Id X .
Note that the eigenvalues of K q -Id are the non positive real numbers (-

λ j q +λ j ; j = 0, • • • n -1)
. Therefore λ max (K q | X , X -Id X ) ≤ λ max (K q -Id) = 0. Equality holds iff there exists an eigenfunction of K q -Id associated to the eigenvalue 0, which is contained in Span(1 {x} ; x ∈ X ). But, 0 is a simple eigenvalue of K q -Id associated to the function 1. Since X = X , we get that λ max (K q | X , X -Id X ) < 0. Thus, K q | X , X -Id X is invertible and by Proposition 1,

det(M ) = det(K q | X , X -Id X ) det(S M (K q | X , X -Id X )) . ( 31
)
Concerning the Schur complement in (31), note that

S M (K q | X , X -Id X ) = K q | X , X -K q | X , X (K q | X , X -Id X ) -1 K q | X , X = Id X + S K q -Id (K q | X , X -Id X ) . K q -Id is a symmetric negative operator in 2 (X , µ) such that λ min (K q -Id) = -λ n-1 q +λ n-1 > -1. Using Proposition 1, we deduce that λ min (S K q -Id (K q | X , X -Id X )) ≥ λ min (K q -Id) > -1 .
Hence det(S M (K q | X , X -Id X )) > 0. We have thus proven that det(M ) = 0.

Once we have defined a basis (ξ x ) x∈X := (φ x, ψ x) x∈ X ,x∈ X of 2 (X , µ), we are able to define the analysis operator

U : 2 (X , µ) → R X f → ( ξ x , f ) x∈X = ( φ x, f , ψ x, f ) x∈ X ,x∈ X .
The hope is that when f is "regular" (a notion still to be defined), the coefficients ( f, ψ x ) x∈ X are "small". Actually, when f is a constant function, these coefficients are equal to zero. We stress the fact that these coefficients are not normalized. f can be reconstructed from the coefficients U (f ). Let ( ξx ) x∈X be the dual basis of (ξ x ) x∈X in 2 (X , µ), defined by ξx , ξ y = δ xy . Then we get the obvious reconstruction formula

f = x∈X ξ x , f ξx .
The approximation of f is then ȳ∈ X φ ȳ, f ξȳ , while the detail of f is y∈ X ψ y, f ξy . The reconstruction of f from its coefficients U (f ) can here be made explicit. This is the content of the following proposition.

The reconstruction formula.

Proposition 10 (Reconstruction formula). For any f ∈ 2 (X , µ), let f ∈ 2 ( X , µ) and f ∈ 2 ( X , µ) be defined by

∀x ∈ X , f (x) = K q f (x) = U (f )(x) ; ∀x ∈ X , f (x) = (K q -Id)f (x) = U (f )(x) . Then, f = Id X -1 q L L X X (-L X X ) -1 (-L X X ) -1 L X X q L -1 X X -Id X f f = R f + Ȓ f , where R = Id X -1 q L (-L X X ) -1 L X X , and Ȓ = L X X (-L X X ) -1 -Id X -q (-L X X ) -1 . (32)
Proof. Note that f f = K q -0 0 0 Id X f . Since K q = (Id -1 q L) -1 , we get

Id - 1 q L f f = Id -Id - 1 q L 0 0 0 Id X f = Id X L X X /q 0 L X X /q f . The inverse of Id X L X X /q 0 L X X /q is Id X L X X (-L X X ) -1 0 q (L X X ) -1 . Hence, f = Id X L X X (-L X X ) -1 0 q (L X X ) -1 Id X -1 q L X X -1 q L X X -1 q L X X Id X -1 q L X X f f , = Id X -1 q L L X X (-L X X ) -1 (-L X X ) -1 L X X q (L X X ) -1 -Id X f f .
In the sequel, we will talk about R as the approximation operator, and about Ȓ as the detail operator. In order to get stable numerical results in the analysis and the reconstruction of the signal f , it is important to control the operator norm of these operators. This is easy for the analysis operator U since for any f ∈ p (X , µ),

U (f ) p,X = K q (f ) -f 1 X p,X ≤ K q (f ) p,X + f p,X .
K q being a probability kernel, symmetric w.r.t. µ, we get by Jensen's inequality,

K q f p p,X = x∈X µ(x) K q f (x) p (33) ≤ x∈X µ(x)K q (|f | p )(x) = 1, K q (|f | p ) = K q (1), |f | p = f p p,X . (34)
Hence for any f ∈ p (X , µ),

U (f ) p,X ≤ 2 f p,X . (35)
The aim of the two following sections is to provide bounds on the norm of the reconstruction operator.

4.3. Detail operator norm. We give in this section a control on the norm of the detail operator Ȓ, in terms of β and γ defined in [START_REF] Elisha | Wavelet decompositions of random forests. Smoothness analysis, sparse approximation and applications[END_REF]. Using Markov property at time τ 1 , one can express 1/β as

1 β = max x∈ X E x H + X -τ 1 .
In the sequel, when A is a subset of X , and f is a function on A, f p,A is the norm of f in p (A, µ A ), where µ A is the conditional probability µ on A: for any B ⊂ A, µ A (B) = µ(B)/µ(A).

Proposition 11 (Detail operator norm). Let X be any proper subset of X , X = X \ X , and let Ȓ be the operator defined in [START_REF] Gavish | Multiscale wavelets on trees, graphs and high dimensional data: theory and applications to semi supervised learning[END_REF]. For any p ≥ 1, for any

f ∈ p ( X , µ X ), Ȓf p,X ≤ α β p/p * + 1 + q γ p 1/p µ( X ) 1/p f p, X , (36)
where p * is the conjugate exponent of p.

The proof is a consequence of Lemmas 12-13 given below.

Lemma 12. For any x ∈ X , and x ∈ X , we define

G X (x, x) = E x H + X 0 1 X(s)=x ds .
Then,

(-L X X ) -1 = (G X ) X X , and L X X (-L X X ) -1 = α(G X ) X X . (37)
The adjoint operator of

(-L X X ) -1 : p ( X , µ) → q ( X , µ)
is the operator

(-L X X ) -1 : q * ( X , µ) → p * ( X , µ) .
The adjoint operator of

L X X (-L X X ) -1 : p ( X , µ) → q ( X , µ)
is the operator

(-L X X ) -1 L X X : q * ( X , µ) → p * ( X , µ) , with (-L X X ) -1 L X X (x, x) = P x [X(H X ) = x] , x ∈ X , x ∈ X . (38)
Proof.

G X (x, x) = E x τ 1 0 1 X(s)=x ds + H + X τ 1 1 X(s)=x ds , = δ x (x)E x (τ 1 ) + z∈X P (x, z)E z H X 0 1 X(s)=x ds by Markov property at time τ 1 , = 1 α δ x (x) + z∈ X P (x, z)E z H + X 0 1 X(s)=x ds .
This can be rewritten as

(G X ) X X = P X X (G X ) X X (G X ) X X = 1 α Id X + P X X (G X ) X X ⇐⇒ -L X X (G X ) X X = Id X (G X ) X X = 1 α L X X (G X ) X X ⇐⇒ (-L X X ) -1 = (G X ) X X L X X (-L X X ) -1 = α(G X ) X X .
Using the symmetry of L with respect to µ, for any g ∈ q ( X , µ) and any f ∈ p ( X , µ),

(-L X X ) -1 f, g X = x,y∈ X µ(x)(-L X X ) -1 (x, y)f (y)g(x) = x,y∈ X µ(y)(-L X X ) -1 (y, x)f (y)g(x) = f, (-L X X ) -1 g X .
In the same way, for all g ∈ q ( X , µ) and all f ∈ p ( X , µ),

L X X (-L X X ) -1 f, g X = x∈ X ;y,z∈ X µ(x)L(x, y)(-L X X ) -1 (y, z)f (z)g(x) = x∈ X ;y,z∈ X µ(z)L(y, x)(-L X X ) -1 (z, y)f (z)g(x) = f, (-L X X ) -1 L X X g X .
(38) was already stated in the proof of Lemma 8 (see ( 26)).

We are now able to give bounds on the norms of the operators involved in the definition of R and Ȓ. Lemma 13. Let X be any proper subset of X , X = X \ X . β and γ being defined by (9), we get

(-L X X ) -1 L X X f p, X ≤ α β µ( X ) µ( X ) 1/p f p, X . (39)
L X X (-L X X ) -1 f p, X ≤ α β 1/p * µ( X ) µ( X ) 1/p f p, X . (40) (-L X X ) -1 f p, X ≤ 1 γ f p, X . (41)
Proof.

µ( X ) (-L X X ) -1 L X X f p p, X = x∈ X µ(x) (-L X X ) -1 L X X f (x) p , = x∈ X µ(x) |E x [f (X(H X ))]| p by (38), ≤ x∈ X µ(x)E x [|f (X(H X ))| p ] by Jensen's inequality , = x∈ X ,ȳ∈ X µ(x)((-L X X ) -1 L X X )(x, ȳ) |f | p (ȳ) by (38), = x∈ X ,ȳ∈ X µ(ȳ)(L X X (-L X X ) -1 )(ȳ, x) |f | p (ȳ) by symmetry, = α ȳ∈ X µ(ȳ) |f (ȳ)| p E ȳ H + X 0 1 X (X(s)) ds by (37), = α ȳ∈ X µ(ȳ) |f (ȳ)| p E ȳ H + X -τ 1 , ≤ α max x∈ X E x H + X -τ 1 µ( X ) f p p, X .
This gives (39). (40) follows from (39), since the operator L X X (-L X X ) -1 from p ( X , µ) to p ( X , µ) is the adjoint operator of the operator (-L X X ) -1 L X X from p * ( X , µ) to p * ( X , µ). Concerning (41), note that by Lemma 12, for any x, y ∈ X , (-L X X ) -1 (x, y) ≥ 0.

µ( X ) (-L X X ) -1 f p p, X = x∈ X µ(x) (-L X X ) -1 f (x) p , = x∈ X µ(x) y∈ X (-L X X ) -1 (x, y)f (y) p , ≤ x∈ X µ(x)   y∈ X (-L X X ) -1 (x, y)   p/p *   y∈ X (-L X X ) -1 (x, y) |f (y)| p   by Hölder's inequality , ≤   max x∈ X y∈ X (-L X X ) -1 (x, y)   p/p *   x∈ X ,y∈ X µ(x)(-L X X ) -1 (x, y) |f (y)| p   , ≤   max x∈ X y∈ X (-L X X ) -1 (x, y)   p/p *   y∈ X µ(y) |f (y)| p x∈ X (-L X X ) -1 (y, x)   by symmetry , ≤   max x∈ X y∈ X (-L X X ) -1 (x, y)   p µ( X ) f p p, X
This ends the proof of (41) since by (37

), y∈ X (-L X X ) -1 (x, y) = E x H + X 0 1 X (X(s)) ds = E x H + X = E x [H X ] .
Proof of Proposition 11. For any f ∈ p ( X , µ X ),

Ȓf p p,X = L X X (-L X X ) -1 -Id X -q (-L X X ) -1 f p p,X = µ( X ) L X X (-L X X ) -1 f p p, X + µ( X ) Id X + q (-L X X ) -1 f p p, X ≤ α β p/p * µ( X ) + 1 + q γ p µ( X ) f p p,
X by ( 40) and (41) .

4.4. Approximation operator norm. This section aims to control the operator norm of R defined in [START_REF] Gavish | Multiscale wavelets on trees, graphs and high dimensional data: theory and applications to semi supervised learning[END_REF].

Proposition 14 (Approximation operator norm). Let X be any proper subset of X , X = X \ X , and let R be the operator defined in [START_REF] Gavish | Multiscale wavelets on trees, graphs and high dimensional data: theory and applications to semi supervised learning[END_REF]. For any p ≥ 1, for any f ∈ p ( X , µ X ),

Rf p,X ≤ 1 + 2 ᾱ q p + α β 1/p µ( X ) 1/p f p, X , (42)
where ᾱ is defined by (28).

Proof. For any f ∈ p ( X , µ X ),

Rf p p,X = Id X -L q (-L X X ) -1 L X X f p p,X = µ( X ) Id X - L q f p p, X + µ( X ) (-L X X ) -1 L X X f p p, X ≤ µ( X ) Id X - L q f p p, X + µ( X ) α β f p p,
X by (39) .

It just remains to prove that for any

f ∈ p ( X , µ X ), Lf p, X ≤ 2 ᾱ f p, X .
By definition of ᾱ, the matrix P = Id X + L ᾱ is a stochastic matrix, which is symmetric w.r.t µ X . Therefore, it is contraction operator from p ( X , µ X ) to p ( X , µ X ). Hence,

Lf p, X = ᾱ (Id X -P )f p, X ≤ 2ᾱ f p, X .
4.5. Size of the detail. In this section, we provide control of the unnormalized detail coefficients f in terms of the regularity of the signal f . Proposition 15. For any p ≥ 1 and any f ∈ p (X , µ),

f p, X = (K q -Id)f p, X ≤ max x K q (x, X ) 1/p q µ( X ) 1/p
Lf p,X .

Proof. Note that by definition of K q , K q -Id = 1 q K q L. Hence,

µ( X ) (K q -Id)f p p, X = 1 q p x∈ X µ(x) E x Lf (X(T q ) p ≤ 1 q p x∈ X µ(x)E x Lf (X(T q ) p , ≤ 1 q p x∈ X ,x∈X µ(x)K q (x, x) |Lf (x)| p , ≤ 1 q p x∈ X ,x∈X µ(x)K q (x, x) |Lf (x)| p by symmetry , ≤ max x K q (x, X ) q p x∈X µ(x) |Lf (x)| p .
4.6. Link between (X , L) and ( X , L). We give bounds on the error in the interwining relation.

Proposition 16. Let X be any proper subset of X , X = X \ X . Recall the definition of β given in [START_REF] Elisha | Wavelet decompositions of random forests. Smoothness analysis, sparse approximation and applications[END_REF]. For any p ≥ 1, and any f ∈ p (X , µ),

L(K q ) X X -(K q ) X X L f p, X ≤ 2q α β 1/p * 1 µ( X ) 1/p f p,X . (43)
Proof. Writing that L is a Schur complement and using that K q L = LK q , we get

L(K q ) X X -(K q ) X X L = L X X (K q ) X X + L X X -L X X -1 L X X (K q ) X X -(LK q ) X X = -L X X (K q ) X X + L X X -L X X -1 L X X (K q ) X X = L X X -L X X -1 L X X (K q ) X X + L X X (K q ) X X = L X X -L X X -1 LK q X X . (44)
It has already been proven in Lemma 13 (see (40)), that for any f ∈ p ( X , µ X ),

L X X -L X X -1 f p, X ≤ α β 1/p * µ( X ) µ( X ) 1/p f p,

X

To get (43), it just remains to prove that for any f ∈ p (X , µ),

LK q f p, X ≤ 2q µ( X ) 1/p f p,X . ( 45
)
But this is a direct consequence of the equality LK q = q (K q -Id), and of the inequality (34).

5. The multiresolution scheme.

We assume in this section that we have at our disposal a decreasing sequence of proper subsets of X :

X 0 = X X 1 • • • X k ∅ ,
and a sequence of non negative real parameters {q i ; i = 0, • • • k -1}. We will discuss later the choice of such sequences. To stick to the previous notations, we define for i ≥ 0,

• Xi = X i+1 ; • Xi = X i \ X i+1 ;
• L 0 = L, and for i ≥ 0, L i+1 = S L i ((L i ) Xi , Xi ). L i is a Markov generator on X i which is symmetric with respect to µ.

• α i = max x i ∈X i (-L i )(x i , x i ); • 1/β i = max xi ∈X i+1 E xi H (i)+ X i+1 -τ (i) 1
, where

H (i)+ A
is the return time on the set A of a Markov process X (i) with generator L i , and τ (i) 1 is its first jump time (exponential time with parameter α i );

• 1/γ i = max xi ∈ Xi E xi H (i) X i+1 ; • for any x i , y i ∈ X i , K i (x i , y i ) = q i (q i Id X i -L i ) -1 (x i , y i ) = P x i X (i) (T q i ) = y i ,
where T q i is an exponential random variable with parameter q i independent of X (i) . Given a signal f i defined on X i , we can define its approximation coefficients at scale i + 1:

∀x i+1 ∈ X i+1 , f i+1 (x i+1 ) = fi (x i+1 ) = Ki (f i )(x i+1 ) , where Ki := (K i ) X i+1 ,X i ,
and its detail coefficients at scale i + 1:

∀x i ∈ Xi , g i+1 (x i ) = fi (x i ) = Ki (f i )(x i ) -f i (x i ) where Ki = (K i ) Xi ,X i , so that f i = Ri f i+1 + Ȓi g i+1 , with Ri = Id X i+1 -1 q i L i+1 ((-L i ) Xi Xi ) -1 (L i ) Xi X i+1
, and Ȓi =

(L i ) X i+1 Xi ((-L i ) Xi Xi ) -1 q i ((L i ) Xi Xi ) -1 -Id Xi .
Iterating the procedure, we get the usual multiresolution scheme:

f 0 = f → f 1 → f 2 • • • → f k g 1 g 2 g k 5.
1. Analysis operator norm. For k ≥ 1, let us consider the analysis operator U k :

U k : p (X ) → p (X k , µ X k ) × p ( Xk-1 , µ Xk-1 ) × • • • × p ( X0 , µ X0 ) f → [f k , g k , g k-1 , • • • , g 1 ] The space p (X k , µ X k ) × p ( Xk-1 , µ Xk-1 ) × • • • × p ( X0 , µ X0 ) is endowed with the norm [f k , g k , g k-1 , • • • , g 1 ] p := µ(X k ) f k p p,X k + k i=1 µ( Xi-1 ) g i p p, Xi-1 1/p .
Proposition 17 (Analysis operator norm). For any f ∈ p (X , µ),

U k (f ) p ≤ 2 1/p * (1 + k) 1/p f p . ( 46
)
Remark: To maintain localization properties of our wavelet basis we do not use any orthogonalization procedure. Then, working with a general graph, we do not expect a conditioning of the multilevel analysis or reconstruction operator that is independent from the graph size. However for p = ∞ our analysis operator norm is independent of the graph size, while, for p < ∞ it depends only on the number of levels, which is typically only logarithmic in the graph size. Making a careful analysis of the following proof, we will see that at least in the case p = 1 our estimation is optimal for the analysis operator we proposed. As far as the multilevel reconstruction operator is concerned, its norm is strongly related with the choice of the q i that we will discuss at the end of Section 6.

Proof.

U k (f ) p p = x k ∈X k µ(x k ) |f k (x k )| p + k i=1 xi-1 ∈ Xi-1 µ(x i-1 ) |g i (x i-1 )| p , = x k ∈X k µ(x k ) ( Kk-1 ... • • • K0 )(f )(x k ) p + k i=1 xi-1 ∈ Xi-1 µ(x i-1 ) (( Ki-1 -Id Xi-1 X i-1 ) Ki-2 • • • K0 )(f )(x i-1 ) p , ≤ x k ∈X k µ(x k ) ( Kk-1 ... • • • K0 )(f )(x k ) p + 2 p-1 k i=1 xi-1 ∈ Xi-1 µ(x i-1 ) ( Ki-1 Ki-2 • • • K0 )(f )(x i-1 ) p + ( Ki-2 • • • K0 )(f )(x i-1 ) p , ≤ 2 p-1 (a k + b k ) ,
where

a k = x k ∈X k µ(x k )( Kk-1 • • • K0 )(|f | p )(x k ) + k i=1 xi-1 ∈ Xi-1 µ(x i-1 )( Ki-1 Ki-2 • • • K0 )(|f | p )(x i-1 ) b k = k i=1 xi-1 ∈ Xi-1 µ(x i-1 )( Ki-2 • • • K0 )(|f | p )(x i-1 ) .
Note that

x k ∈X k µ(x k )( Kk-1 • • • K0 )(|f | p )(x k ) + xk-1 ∈ Xk-1 µ(x k-1 )( Kk-1 Kk-2 • • • K0 )(|f | p )(x k-1 ) = x k-1 ∈X k-1 µ(x k-1 )(K k-1 Kk-2 • • • K0 )(|f | p )(x k-1 ) = x k-1 ∈X k-1 µ(x k-1 )( Kk-2 • • • K0 )(|f | p )(x k-1 )
by symmetry. Therefore, a k = a k-1 = a 1 = f p p,X . Concerning b k , using symmetry, it holds

k i=1 xi-1 ∈ Xi-1 µ(x i-1 )( Ki-2 • • • K0 )(|f | p )(x i-1 ) = x 0 ∈X 0 µ(x 0 )|f (x 0 )| p k i=1 x 1 ∈X 1 • • • x i-2 ∈X i-2 xi-1 ∈ Xi-1 K 0 (x 0 , x 1 ) • • • K i-3 (x i-3 , x i-2 )K i-2 (x i-2 , xi-1 ) ≤ k f p p,X ,
This ends the proof of (46). We finally give an example for which one can check that in the case p = 1 our estimation ( 46) is optimal. Consider the graph of size n = l + 1 with vertex set X = {0, 1, . . . , l} and transition rates

w(x, y) =    ε 2x+1 if y = x + 1, ε 2x if y = x -1, 0
otherwise, together with the signal

f (x) = ε 2x-l , 0 ≤ x ≤ l, for ε 1. Choosing q i = q i = ε i+ 3 2 for 0 ≤ i ≤ k, the random subsets X 0 ⊃ X 1 ⊃ X 2 ⊃ • • • sampled through Wilson'
s algorithm are typically given by

X i = {i, . . . , l}, i ≤ k,
and the inequalities of the previous proof are turned into equalities when p = 1 and in the limit ε → 0.

Approximation error. Given the coefficients

[f k , g k , g k-1 , • • • , g 1 ], the reconstruction of f = f 0 is f = f 0 = R0 f 1 + Ȓ0 g 1 = R0 R1 f 2 + R0 Ȓ1 g 2 + Ȓ0 g 1 = R0 R1 • • • Rk-1 f k + k-1 j=0 ( R0 • • • Rj-1 ) Ȓj g j+1 .
The approximation of f at scale k is thus R0 R1 • • • Rk-1 f k and we have the following Jackson's type inequality:

Proposition 18 (Jackson's inequality).

For i ∈ {0, • • • , k -1} let Ri p ( Ȓi p respectively) denote the norm operator of Ri : p (X i+1 , µ X i+1 ) → p (X i , µ X i ) ( Ȓi : p ( Xi , µ Xi ) → p (X i , µ X i ) respectively). We set ri,p := µ(X i ) µ(X i+1 ) 1/p Ri p , ȓi,p := µ(X i ) µ( Xi ) 1/p Ȓi p
, and e i,p := µ(X i+1 ) µ(X i )

1/p L i+1 Ki -Ki L i p .
For any p ≥ 1 and any f ∈ p (X , µ),

f -R0 R1 • • • Rk-1 f k p,X ≤    k-1 j=0 r0,p • • • rj-1,p ȓj,p 1 q j    Lf p,X +    k-1 j=0 r0,p • • • rj-1,p ȓj,p 1 q j j-1 l=0 e l,p    f p,X . (47)
Note that it has been proven in Propositions 14, 11 and 16 that

ri,p ≤ 1 + 2 α i+1 q i p + α i β i 1/p , ȓi,p ≤ α i β i p/p * + 1 + q i γ i p 1/p
, and e i,p ≤ 2q i

α i β i 1/p * .
Moreover, we would like to stress the fact that the term involving f p,X in ( 47) is linked to the error in the intertwining relation, and would disappear if this relation was exact.

In case the intertwining relation error is small, we can interpret our inequality (47) as a Jackson type inequality. Indeed in the setting of approximation theory, Jackson type inequalities relate the smoothness of a function with the approximation properties of a multiscale basis. Let us give some more details on these classical results in the setting of wavelet analysis. Assume φ is a real valued scaling function, in the Schwartz class or compactly supported, i.e a function such that one can find a wavelet ψ such that {φ(.

-k), k ∈ Z} ∪ {x → 2 j/2 ψ(2 j . -k), j ∈ N, k ∈ Z} is an orthonormal basis of L 2 (
R), a so-called wavelet basis. The construction of such kind of basis has a deep relationship with the discrete scheme described in the introduction 1.1 as it is explained in [START_REF] Mallat | A wavelet tour of signal processing[END_REF] (see also [START_REF] Daubechies | Ten lectures on wavelets[END_REF] for details). We can define for j ≥ 0, V j := Span{φ j,k = 2 j/2 φ(2 j . -k), k ∈ Z} and P j f the orthogonal projector on V j .

The remarkable properties of wavelet basis make it possible to analyse functions in other functional spaces than L 2 (R). Suppose the function ψ of our basis has n vanishing moments with n ∈ N * (by construction ψ has always at least one vanishing moment). Let 1 ≤ p ≤ ∞ and W n,p denote the Sobolev space of functions f in L p whose derivative of order n belongs to L p . Let |f | W n,p = d n f dx n L p . We have indeed the following result in the classical setting (see [START_REF] Cohen | Numerical analysis of wavelet methods[END_REF] for a detailed proof).

Theorem 19 (Classical Jackson's inequality).

There exists a constant C > 0 such that for all f ∈ W n,p and j ≥ 0,

f -P j f L p ≤ C2 -nj |f | W n,p . (48)
As one can see, our Jackson's inequality (47) relates also the smoothness of the signal f on the graph (measured by Lf p,X ) provided that the error in the intertwining equation is small. The smallest is Lf p,X , the closer to f should be the approximation part. As in the classical inequality (48), the constant in front of the Laplacian term depends on the level of approximation and gets worse when this level gets coarser.

We can now prove Proposition 18.

Proof.

f -R0 R1 • • • Rk-1 f k p,X 0 ≤ k-1 j=0 ( R0 • • • Rj-1 ) Ȓj g j+1 p,X 0 , ≤ k-1 j=0 j-1 i=0 Ri p Ȓj p g j+1 p, Xj , ≤ k-1 j=0 j-1 i=0 Ri p Ȓj p µ(X j ) µ( Xj )
1/p L j f j p,X j q j , by Proposition 15. Since µ(X 0 ) = 1, this leads to

f -R0 R1 • • • Rk-1 f k p,X 0 ≤ k-1 j=0 r0 • • • rj-1 ȓj µ(X j ) 1/p q j L j f j p,X j .
Now,

L j f j = L j Kj-1 f j-1 = L j Kj-1 • • • K0 f 0 = Kj-1 • • • K0 L 0 f 0 + j-1 l=0 Kj-1 • • • Kl+1 L l+1 Kl -Kl L l Kl-1 • • • K0 f 0 . For any l ∈ {0, • • • , j -1}, Kl is an operator from p (X l , µ X l ) to p (X l+1 , µ X l+1 ), and (34) states that Kl p ≤ µ(X l ) µ(X l+1 ) 1/p
. Hence,

L j f j p,X j ≤ µ(X 0 ) µ(X j ) 1/p L 0 f 0 p,X 0 + j-1 l=0 µ(X l+1 ) µ(X j ) 1/p L l+1 Kl -Kl L l p µ(X 0 ) µ(X l ) 1/p f 0 p,X = µ(X 0 ) µ(X j ) 1/p L 0 f 0 p,X 0 + j-1 l=0 e l f 0 p,X ,
This gives (47).

6. About the choice of the parameters q, q .

The aim of this section is to give a guideline in the choice of the various parameters involved in the multiresolution scheme. The requirements we would like to achieve are the following ones:

(1) At each step of the downsampling, we would like to keep a fixed proportion of points in the current set. This would ensure that the number of steps k in the multiresolution scheme is of order log(|X |). In view of Proposition 6, this could be achieved by taking for all j ≥ 0,

q j ∈ [θ 1 α j , θ 2 α j ] for some θ 1 , θ 2 .
With this choice and the choice X j+1 = ρ(Φ (j) ) (where Φ (j) is the random s.o.f on X j associated to L j ) , one has for instance,

∀j ≥ 0 , |X j | q j q j + α j ≤ E q j [|X j+1 | |X j ] ≤ |X j | -|R j,r | 1 - rα j q j + 2α j ,
where R j,r = {x ∈ X j ; w j (x) ≥ rα j }. (2) To ensure numerical stability of the reconstruction operator, we would like to control at each step the norms of Rj and Ȓj . Concerning Rj (see (42)), this requires α j+1 /q j and α j /β j to be small. As for Ȓj (see (36)), one has moreover to ensure that q j /γ j is small. Note that once X j has been fixed, α j+1 , β j and γ j only depend on the choice of X j+1 , hence of q j . Since α j+1 /q j and q j /γ j should be small, the product

α j+1 q j q j γ j = α j+1
γ j should also be small, and one possible choice for q j is to minimize α j+1 /γ j .

(3) In order to get a good approximation error, one would also like the error in the intertwining problem between L j and L j+1 to be small. Referring to (47), this is achieved if at each step q j (α j /β j ) 1/p * is small with respect to the natural unit α j .

To sum up, X j being chosen, one would like to choose (q j , q j ) such that (C1): q j ∈ [θ 1 α j , θ 2 α j ];

(C2): α j /β j is small; (C3): α j+1 /q j is small; (C4): q j /γ j is small; (C5): q j α j /β j is small (consider the case p = +∞) with respect to α j .

In this respect, we need some estimates on α j+1 , γ j and β j in terms of q j . The next sections are devoted to this task.

6.1. Estimate on ᾱ.

Proposition 20. Assume that X = ρ(Φ) and let X = X \ X . Then, for any

m ∈ {1, • • • , |X | + 1}, E q ᾱ1 | X |=m ≥ E q   1 X x∈ X w(x)1 | X |=m   = q |X | -m + 1 m P q X = m -1 ; (49)
E q [ᾱ] ≥ E q   1 X x∈ X w(x)   = qE q   X X + 1   . (50)
Note that (49) is trivially true for m = |X | + 1. It is also trivially true for m = 1, since in this case L is a scalar, which is equal to 0 since L is a Markov generator.

Assuming that ᾱ = max x∈ X w(x) is not too far from the mean of the (w(x), x ∈ X ), Proposition 20 can be used to get an idea of the dependence of ᾱ with respect to q (or m). Indeed, the term | X | | X |+1 can be numerically estimated as a function of q, since there is a coupling allowing to sample the random forest for all the values of q at the same time. Getting an idea of ᾱ as a function of q from the simulations, is more time consuming since it requires the computation of a Schur complement for all values of q, and this cannot been done inductively since ρ(Φ) is not an increasing set w.r.t. q.

Proof. (50) is a direct consequence of (49) after summing over m.

E q   1 X x∈ X w(x)1 | X |=m   = 1 m x∈X E q 1 x∈ X w(x)1 | X |=m .
Hence, the result is a consequence of the following identity: for any x ∈ X , and any m ∈ {1, • • • , n + 1},

E q 1 x∈ X w(x)1 | X |=m = qP q X = m -1; x / ∈ X . (51)
The rest of the proof is devoted to proving (51). Remind that

w(x) = -L(x, x) = y∈ X ,y =x L(x, y) = y∈ X ,y =x L(x, y) + L X X (-L X X ) -1 L X X (x, y) = y∈ X ,y =x w(x, y) + y∈ X ,y =x z 1 ,z 2 / ∈ X w(x, z 1 )(-L X , X ) -1 (z 1 , z 2 )w(z 2 , y)
Remind also from Lemma 12 (see (37)) that for

z 1 , z 2 ∈ X , (-L X , X ) -1 (z 1 , z 2 ) = G X (z 1 , z 2 )
. This gives an explicit expression of w(x) as a function of X . Therefore,

E q 1 x∈ X w(x)1 | X |=m = R x,|R|=m P q X = R          y∈R,y =x w(x, y) + y∈R,y =x z 1 ,z 2 / ∈R w(x, z 1 )G R (z 1 , z 2 )w(z 2 , y)          Now, P q X = R = 1 Z(q) φ s.o.f,ρ(φ)=R q m w(φ) = q m Z(q) Z R (0) , (52)
where we set Z R (0) := φ s.o.f,ρ(φ)=R w(φ). In addition, it is proven in [START_REF] Avena | Two applications of random spanning forests[END_REF] (see Lemma 3.1, or Appendix B of the preprint) that for any subset R of X , and any

z 1 , z 2 / ∈ R, G R (z 1 , z 2 ) = E z 1 H + R 0 1 X(s)=z 2 ds = 1 Z R (0) φ s.o.f.ρ(φ)=R∪{z 2 } z 1 φ z 2 w(φ) , (53)
where z 1 φ z 2 means that z 1 is contained in the tree of φ, whose root is z 2 . Hence,

Z(q)q -m E q 1 x∈ X w(x)1 | X |=m = T 1 + T 2 , where • T 1 := R x;|R|=m Z R (0) y∈R;y =x w(x, y) = φ,ρ(φ) x |ρ(φ)|=m y∈ρ(φ);y =x w(φ)w(x, y) = (φ,y)∈S 1 (m,x) w(φ)w(x, y) where we set S 1 (m, x) := {(φ, y) such that ρ(φ) x, |ρ(φ)| = m, y ∈ ρ(φ), y = x}. • T 2 := R x;|R|=m y∈R;y =x z 1 ,z 2 / ∈R φ s.o.f.ρ(φ)=R∪{z 2 } z 1 φ z 2 w(x, z 1 )w(φ)w(z 2 , y) = (R,y,z 1 ,z 2 ,φ)∈S 2 (m,x) w(x, z 1 )w(φ)w(z 2 , y) where we set S 2 (m, x) := (R, y, z 1 , z 2 , φ) such that R x, |R| = m, y ∈ R, y = x, z 1 / ∈ R, z 2 / ∈ R, ρ(φ) = R ∪ {z 2 } , z 1 φ z 2 .
Let us compute T 1 . To any (φ, y) ∈ S 1 (m, x), we associate the forest φ 1 = φ ∪ {(x, y)}, so that w(φ 1 ) = w(φ)w(x, y). φ 1 is the forest obtained from φ by hanging the tree with root x to the root y. Hence |ρ(φ 1 )| = m -1, and the distance d(x, ρ(φ 1 )) from x to the roots of φ 1 , is equal to 1. Otherwise stated, S 1 (m, x) is sent by this operation to

F 1 (m, x) := {φ 1 such that |ρ(φ 1 )| = m -1, d(x, ρ(φ 1 )) = 1} .
Note that this correspondence is one to one. Starting from φ 1 ∈ F 1 (m, x), we recover y as the root of the tree containing x, and φ is then obtained by cutting the edge (x, y). Thus,

T 1 = φ 1 ∈F 1 (m,x)
w(φ 1 ) .

Let us now turn our attention to T 2 . To any (R, y, z 1 , z 2 , φ) ∈ S 2 (m, x), we associate the forest φ 2 = φ ∪ {(x, z 1 ), (z 2 , y)}, so that w(φ 2 ) = w(x, z 1 )w(φ)w(z 2 , y). φ 2 is the forest obtained from φ by first hanging the tree with root x to z 1 , so that x is now in the tree with root z 2 , and then z 2 is attached to y. Hence |ρ(φ 2 )| = |ρ(φ)| -2 = m -1, and d(x, ρ(φ 2 )) ≥ 2 (it could happen that z 1 = z 2 ). Otherwise stated, S 2 (m, x) is sent to

F 2 (m, x) := {φ 2 such that |ρ(φ 2 )| = m -1, d(x, ρ(φ 2 )) ≥ 2} .
Note that this correspondence is one to one. Starting from φ 2 ∈ F 2 (m, x), we recover y as the root of the tree containing x; z 1 is the point following x in the path going from x to y; z 2 is the point preceding y in the path going from x to y; R = ρ(φ 2 ) ∪ {x}, and φ is obtained from φ 2 by cutting the edges (x, z 1 ), (z 2 , y). Thus,

T 2 = φ 2 ∈F 2 (m,x) w(φ 2 ) .
It is clear that F 1 (m, x) and F 2 (m, x) are disjoint sets and that

F 1 (m, x) ∪ F 2 (m, x) = {φ such that |ρ(φ)| = m -1; x / ∈ ρ(φ)} .
This leads to

E q 1 x∈ X w(x)1 | X |=m = q m Z(q) φ,|ρ(φ)|=m-1, x / ∈ρ(φ) w(φ) = qP q [|ρ(Φ)| = m -1, x / ∈ ρ(Φ)]
This ends the proof of (51) and of Proposition 20.

Estimate on β.

Proposition 21. Assume that X = ρ(Φ) and let X = X \ X . Then, for any m ∈ {1, • • • , |X |},

E q 1 β 1 | X |=m ≥ E q   1 X x∈ X z∈ X P (x, z)E z (H X )1 | X |=m   = |X | -m αm P q X = m ; (54) E q 1 β ≥ E q   X α X   . (55)
Proof. (55) is a direct consequence of (54) after summing over m.

E q   1 X x∈ X z∈ X P (x, z)E z (H X )1 | X |=m   = 1 m x∈X E q   1 x∈ X z∈ X P (x, z)E z (H X )1 | X |=m   .
Hence, the result is a consequence of the following identity: for any x ∈ X , and any m ∈ {1, • • • , |X |},

E q   1 x∈ X z∈ X P (x, z)E z (H X )1 | X |=m   = 1 α P q X = m; x / ∈ X . ( 56
)
The rest of the proof is devoted to state (56). For x ∈ X and z ∈ X , P (x, z) = w(x, z)/α, and summing over z 2 in (53) leads to:

∀z / ∈ R , E z (H R ) = 1 Z R (0) y / ∈R φ s.o.f.ρ(φ)=R∪{y} z φ y w(φ) (57)
Therefore, using (52),

E q 1 x∈ X z∈ X P (x, z)E z (H X )1 | X |=m = q m αZ(q) R x,|R|=m z / ∈R;y / ∈R φ s.o.f.ρ(φ)=R∪{y} z φ y w(x, z)w(φ) = q m αZ(q) (R,y,z,φ)∈S 3 (m,x) w(x, z)w(φ) ,
where

S 3 (m, x) := {(R, y, z, φ) such that R x, |R| = m, y / ∈ R, z / ∈ R; ρ(φ) = R ∪ {y} , z φ y} .
To any (R, y, z, φ) ∈ S 3 (m, x), we associate the forest φ 3 = φ∪{(x, z)}, so that w(φ 3 ) = w(x, z)w(φ). φ 3 is the forest obtained from φ by hanging the tree with root x to z, so that x is now in the tree with root y. Hence |ρ(φ 3 )| = |ρ(φ)| -1 = m, and d(x, ρ(φ 3 )) ≥ 1 (it could happen that z = y). S 3 (m, x) is thus sent to F 3 (m, x) := {φ such that |ρ(φ)| = m, x / ∈ ρ(φ)}. Note that this correspondence is one to one. Starting from φ 3 ∈ F 3 (m, x), we recover y as the root of the tree containing x in φ 3 ; z is the point following x in the path going from x to y in φ 3 ; φ is obtained from φ 3 by cutting the edge (x, z) so that x is a root of φ, and R = ρ(φ) \ {y}. Thus,

E q   1 x∈ X z∈ X P (x, z)E z (H X )1 | X |=m   = q m αZ(q) φ∈F 3 (m,x) w(φ) = 1 α P q [|ρ(Φ)| = m; x / ∈ ρ(Φ)] .
This ends the proof of (56) and of Proposition 21.

6.3. Estimate on γ.

Proposition 22. Assume that X = ρ(Φ) and let X = X \ X . Then, for any m ∈ {0, • • • , |X |},

E q 1 γ 1 | X |=m ≥ E q   1 X x∈ X E x(H X )1 | X |=m   = |X | |X | -m 1 q P q X = m + 1 ;
(58)

E q 1 γ ≥ 1 q E q   |X | X + 1   . ( 59
)
Note that (58) is trivially true for m = 0, since in this case 1/γ = +∞. It is also true for m = |X | with the convention 0/0 = 0 in the right hand side, since in this case 1/γ = 0.

Proof. Here again, (59) is a direct consequence of (58) after summing over m.

E q 1 | X | x∈ X E x(H X )1 | X |=m = 1 |X |-m x∈X E q 1 x / ∈ X E x (H X )1 | X |=m = 1 |X |-m x∈X R;x / ∈R,|R|=m P q ( X = R)E x (H R ) = q m (|X |-m)Z(q) x∈X R;x / ∈R,|R|=m y / ∈R φ s.o.f.ρ(φ)=R∪{y} x φ y w(φ) = q m (|X |-m)Z(q) x∈X φ s.o.f.|ρ(φ)|=m+1 w(φ) = 1 q(|X |-m) x∈X φ s.o.f.|ρ(φ)|=m+1 P q (Φ = φ) = |X | q(|X |-m) P q (|ρ(Φ)| = m + 1) .
This ends the proof of (58) and of Proposition 22.

6.4. Discussion on the choice of X , q, q . As was said in the beginning of the section, at each step of the multiresolution scheme, one would want to choose (q, q ) so as to fulfill the constraints (C1-5). If (C3) and (C4) are satisfied, then the quantity ᾱ/γ = (ᾱ/q )(q /γ) has to be small. In the same way, if (C3) and (C5) are satisfied, then ᾱ/β is small. One has therefore to choose X at each step in such a way that ᾱ/γ and/or ᾱ/β are small. A possible choice for q is then to minimize the functions q → E q (ᾱ/γ), or q → E q ( ᾱ/β) in some interval [θ 1 α; θ 2 α] (in accordance to (C1)). These functions could be estimated using a Monte Carlo method. However, this can be quite costly, since for instance, the computation of ᾱ needs the computation of a Schur complement. If one can afford such a computation one time, doing it several times to get a Monte Carlo estimation is time consuming. Instead, in the simulations presented in section 8, we used the estimates (50) and (55). Denote

α(q) := qE q   X 1 + X   ; 1 β(q) := 1 α E q   X X   ; 1 γ(q) := 1 q E q   |X | 1 + X   .
These functions are quite easy to estimate by a Monte Carlo method. Indeed, a coalescencefragmentation process has been proposed in [START_REF] Avena | Two applications of random spanning forests[END_REF], allowing to couple the random forests sampled according to π q for various values of q. In this process, t = 1/q is seen as a time variable. Starting from t 0 , the process begins to sample a random forest whose law is π 1/t 0 , and this forest evolves in a Markovian way with fragmentations or coalescences of trees. The marginal of this forest process (Φ(t), t ≥ t 0 ) is at each time t, equal to π 1/t . We refer the reader to [START_REF] Avena | Two applications of random spanning forests[END_REF] for further details on the construction. Using this coalescence-fragmentation process, it is easy to sample the random function q → |ρ(Φ(1/q))|, and to estimate the functions α(q), β(q), γ(q) by empirical means over a sample. One can then choose q as a minimizer of q ∈ [θ 1 α; θ 2 α] → α(q)/ β(q) or q ∈ [θ 1 α; θ 2 α] → α(q)/γ(q). It can however be proven that α(q)/γ(q) is a decreasing function of q; and the reconstruction formulas show that numerical stability depends, just like our Jackson's inequality, on the product of the Rj (rather than Ȓj ), which is associated with 1/ β (rather than 1/γ). Therefore we decided in the simulations of section 8, to minimize the function q ∈ [θ 1 α; θ 2 α] → α(q)/ β(q), with θ 1 = 1/8 and θ 2 = 1.

Turning back to the choice of q , we could choose it so that the two terms of the product (ᾱ/q )(q /β) = ᾱ/β are of the same order. This leads to q = θ 1 √ ᾱβ. With this choice, one get, using (42), (50) and (55)

R ∞,∞ ≤ 1 + 2 θ 1 ᾱ β 1 + 2 θ 1 q α X 1 + X X X ≤ 1 + 2 √ θ 2 θ 1 X X ≤ |X | X if we choose θ 1 = 2 √ θ 2 .
This value of θ 1 ensures that the multiplicative term appearing in Jackson's inequality (47) for p = +∞ behaves nicely:

j-1 i=0 1 + 2 α i+1 q i ≤ j-1 i=0 |X i | |X i+1 | = |X 0 | |X j | . (60)
Another possible choice (made in the simulations results of section 8) is to ensure (60) by setting

1 + 2 α i+1 q i = |X i | |X i+1 | ⇐⇒ q i = 2α i+1 |X i+1 | Xi .
7. A summary of the procedure, and computational issues.

7.1. Summary. Fix two real numbers θ 1 , θ 2 (θ 1 = 1/8 and θ 2 = 1 in our simulations). Starting from (X , L, α, µ), one iteration of our multiresolution scheme goes through the following steps:

(1) Choice of q. Make N i.i.d. draws (N = 1 in our simulations) of the coalescence-fragmentation forest process (Φ

(i) t ; t ∈ [1/(θ 2 α); 1/(θ 1 α)], i = 1, • • • , N ).
From this, estimate the functions

α(q) q 1 N N i=1 |X | -ρ(Φ (i) 1/q ) 1 + ρ(Φ (i) 1/q ) , and
1 β(q) 1 N α N i=1 |X | -ρ(Φ (i) 1/q ) ρ(Φ (i) 1/q ) .
Choose q opt as the minimizer in [θ 1 α; θ 2 α] of the estimation of α(q) β(q) . (2) Once q opt is chosen, draw one forest Φ from π qopt . Take X = ρ(Φ), X = X \ X , μ = µ(•| X ).

(3) Compute L = L X , X -L X , X L X , X -1
L X , X , and ᾱ = max x∈ X (-L(x, x)).

(4) Choose q = 2ᾱ X / X . Compute K q = q (q Id -L) -1 . This can be done by a polynomial approximation. The scaling functions are then

x → φ x(x) = K q (x, x) µ(x) ,
while the wavelets are

x → ψ x(x) = (K q -Id)(x, x) µ(x) .
Compute also the reconstruction operators R and Ȓ. (5) Go on with ( X , L, μ, ᾱ) in place of (X , L, µ, α). 7.2. Numerical issues. The three main computational issues at each step of the algorithm are:

(1) the N sampling of the coalescence-fragmentation forest process (Φ [START_REF] Marchal | Loop-erased random walks, spanning trees and Hamiltonian cycles[END_REF], the mean time to sample Φ q is equal to n-1 i=0 1 q+λ i , corresponding to a number of random jumps for the Markov chain of order n-1 i=0 α q+λ i ≤ α|X | q . Hence, the number of operations to sample N realisations of the forest process is of order N |X | /θ 1 .

(i) t ) t∈[ 1 θ 2 α ; 1 θ 1 α ];i=1,••• ,N . As proved in
(2) the computation of the Schur complement. This requires to compute the inverse (-L X , X ) -1 .

Note that

(-L X , X ) -1 = 1 α (Id X -P X , X ) -1 = 1 α +∞ k=0 (P X , X ) k . (61)
Since (P X , X ) k (x, y) is the probability for the discrete time Markov Chain to go from x to y in k steps without entering X , and since X is "well spread", this probability is fast decaying in k, so that a good approximation of (-L X , X ) -1 is obtained by truncating the series in (61) to the first terms.

The computation of (-L X , X ) -1 is then fast if the matrix -L X , X is sparse. This may be satisfied by L at the beginning of the algorithm, but even if L is sparse, this is usually no more the case for L (see Figure 2). Therefore, after the computation of L, the authors of [START_REF] Hammond | Wavelets on graphs via spectral graph theory[END_REF][START_REF] Shuman | A multiscale pyramid transform for graph signals[END_REF] perform a sparsification method. In our frame, it would be convenient that the sparsification step does not alter too much the error bound we get on (K q ) X ,X L -L(K q ) X ,X . Denote by Ls the sparsified version of L, and let us consider the interwining error:

(K q ) X ,X L -Ls (K q ) X ,X ∞,∞ ≤ (K q ) X ,X L -L(K q ) X ,X ∞,∞ + ( L -Ls )(K q ) X ,X ∞,∞ ≤ (K q ) X ,X L -L(K q ) X ,X ∞,∞ + L -Ls ∞,∞ ≤ (K q ) X ,X L -L(K q ) X ,X ∞,∞ + max x∈ X ȳ∈ X L(x, ȳ) -Ls (x, ȳ) .
Therefore, the sparsification should for instance satisfy

max x∈ X ȳ∈ X L(x, ȳ) -Ls (x, ȳ) ≤ (K q ) X ,X L -L(K q ) X ,X ∞,∞ .
By building Ls with this constraint only, and replacing L by L s , we can however locally deteriorate a good intertwining approximation, i.e., increase a lot a small error

(x) = x∈X K q L (x, x) -L(K q ) X ,X (x, x)
by allowing this error to grow up to the order of = max x∈ X (x). Since, at least from a probabilistic point of view (see [START_REF] Avena | Approximate and exact solutions of intertwining equations through random forests[END_REF] and its appendix), the collection of these local errors (x) is more meaningful than the global error we follow a more local and restrictive approach. For θ ≥ 1 (θ = 4 in our simulation, except in Figure 3 where we show the effect of different choices of θ in our sparsification scheme), we set to 0 a maximal number of coefficients on each row x of L until the off-diagonal suppressed weight remains below a targeted error level (x)ᾱ/2θα (the coefficients of L and L are naturally on scale α and ᾱ respectively). But this has to be done in a symmetric way to ensure the symmetry of Ls . In addition one has to make Ls a Markov generator. To satisfy all these constraints, we proceed in the following way. We first sort in non-decreasing order the off-diagonal pairs (x, ȳ) (with x = ȳ) according to the quantity max( w(x, ȳ), w(ȳ, x)). We then set δ x = 0 for all x ∈ X and screen the ordered off-diagonal pairs. Each time (x, ȳ) is such that δ x + w(x, ȳ) ≤ (x) and δ ȳ + w(ȳ, x) ≤ (ȳ), we set to 0 both w(x, ȳ) and w(ȳ, x). We also increase δ x by w(x, ȳ), and δ ȳ by w(ȳ, x). We finally adjust the diagonal coefficients L(x, x) by subtracting them the associated removed weight δ x. This ensures

max x∈ X ȳ∈ X L(x, ȳ) -Ls (x, ȳ) ≤ ᾱ θα (K q ) X ,X L -L(K q ) X ,X ∞,∞ ≤ (K q ) X ,X L -L(K q ) X ,X ∞,∞ .
In doing so, we can lose the irreducibility property, i.e., we can get a disconnected reduced graph associated with ( X , Ls ). This actually does not raise any particular difficulty. The only use we made of irreducibility was to define without ambiguity our reference measure µ. When irreducibility is lost, we simply proceed by using the natural restriction of µ as reference measure. (3) the computation of K q . As in [START_REF] Shuman | A multiscale pyramid transform for graph signals[END_REF], this can be done by Chebyshev polynomial approximation. Apart from computational gain, using polynomial approximation has the advantage to produce filters K q with a good space localization as soon as L is nearly diagonal.

Numerical Results

This section is devoted to numerical illustrations of our multiresolution scheme, referred to as "the intertwining wavelets multiresolution". We show the results of some downsampling steps on Minnesota roads network (cf Figures 2 and3) containing 2642 vertices, and use the multiresolution schemes to analyse and compress the three benchmark signals of Figure 1. 8.1. Downsampling of the Minnesota roads network. Figure 2 shows the result of two levels of forest's roots sampling, combined with the weighting procedure through Schur's complement computation without sparsification. It illustrates the loss of sparsity of the weighting procedure. In Figure 3, we used the sparsification method proposed in section 7.2 with three values of the parameter θ. On these graphs, the width of one edge is proportional to its weight. 8.2. Analysis. 8.2.1. The line. We analyse the signal of Figure 1(a) using our multiresolution scheme. The results after one step are presented in Figure 4, where we can see that the big detail coefficients are located at the discontinuities of the original signal. We also compare our scheme with a classical wavelet scheme involving Daubechies12 wavelets. The results are given in Figure 5. After two steps , we end up with 370 approximation coefficients (f 2 ), instead of 256 for the classical wavelets. In both cases, the number of non vanishing detail coefficients is small. Since the intertwining wavelets ψ x are unnormalized, we normalize detail coefficients in the compression problem. More precisely, given a signal f , the unnormalized coefficients after one multiresolution step are f and f , from which we can reconstruct f by

f = R f + Ȓ f = x∈ X f (x) ξx + x∈ X f (x) ξx ,
where the dual basis ( ξx , ξx ) x∈ X ,x∈ X corresponds to the columns of matrices R and Ȓ. Unlike the classical wavelets, the basis (ξ x , x ∈ X ) and its dual basis ( ξx , x ∈ X ) are not orthonormal ones. Therefore to compress our signal, we truncate the "normalized coefficients" f (x) ξx . In a similar way, after k multiresolution steps, the unnormalized coefficients are [f k , g k , g k-1 , • • • , g 1 ], from which we can reconstruct f by:

f = R0 • • • Rk-1 f k + k-1 j=0 R0 • • • Rj-1 Ȓj g j (62) = x k ∈X k f k (x k ) ξ(k,k) x k + k-1 j=0 xj-1 ∈ Xj-1 g j (x j-1 ) ξ(k,j) xj-1 (63)
where ( ξ(k,k)

x k , x k ∈ X k ) are the columns of the matrix R0 • • • Rk-1 , while ( ξ(k,j) xj-1 , xj-1 ∈ Xj-1 ) are the colums of the matrix R0 • • • Rj-1 Ȓj (j = 0, • • • k -1). Given a threshold > 0, the compressed version of f is

f c = x k ∈X k f k (x k ) ξ(k,k) x k + k-1 j=0 xj-1 ∈ Xj-1
g j (x j-1 ) 1 |g j (x j-1 )| ξ(k,j) Another way to compress f is to keep a fixed percentage of the highest (in absolute value) normalized detail coefficients. This is the way we have done our compression experiments.

8.3.1. The line. We compared the compression results of our method with those obtained using classical Daubechies12 wavelets. We let our algorithm evolve until we get an approximation of size less than 16. In the experiment, this was achieved after 20 steps, and led to 15 approximation coefficients. For classical wavelets, we get 16 approximation coefficients after 6 steps. For both methods, a given proportion p of the details coefficients are kept to compute the compressed signals f c . Figure 8 presents the relative errors f -f c 2 / f 2 in terms of p, and shows the good behavior of the coefficients for the two methods. More precisely, we ran the spectral graph wavelets pyramidal algorithm for 7 steps, resulting in 20 approximations coefficients among approximately 2 × 2642 stored ones. We also ran our intertwining wavelets multiresolution until getting approximately the same number of approximation coefficients to get a fair comparison. This took 16 steps, resulting in 16 approximation coefficients. We kept then the same number of the biggest coefficients to construct the compressed version of the signal.

8.3.3. Sensor graph. For the signal on the sensor graph we compare once again in the same way our method with the Pyramid algorithm on 3 steps. In this case we also included the results of our procedure without sparsification, since they surprisingly show that the sparsification improve the results. In Figure 12, we compare the relative compression errors in terms of the number of kept coefficients for the Pyramid algorithm and our method with and without sparsification. In red, the error using intertwining wavelets. In blue, error using the spectral graph wavelets pyramidal algorithm. In the first graph, we do not sparsify the graph, while we perform sparsification on the second.

2. 2 .

 2 Schur complement. Let M be a matrix of size n = p + r and let M = A B C D be its block decomposition, A being a square matrix of size p and D a square matrix of size r. If D is invertible, the Schur complement of D in M is the square matrix of size p defined by S M (D) := A -BD -1 C .

Figure 1 .

 1 Figure 1. Original signals: (a) Signal on the line presented in [12]. (b) Signal on the Minnesota roads network: sign of the first Fourier mode e 1 . (c) Signal from[START_REF] Shuman | A multiscale pyramid transform for graph signals[END_REF] on the sensor graph from the GSPBox.

Figure 2 . 8 . 2 . 2 . 8 . 3 .

 282283 Figure 2. One sequence of Minnesota's subgraphs, without sparsification. (a) Original graph. (b) First level of downsampling. (c) Second level of downsampling.

Figure 3 .

 3 Figure 3. One sequence of Minnesota's subgraphs, using sparsification. On the first line, θ = 2, while θ = 4 and θ = 8 on the second and third ones. (a) First level of downsampling. (b) Second level of downsampling. (c) Third level of downsampling.

Figure 4 .

 4 Figure 4. One step of the analysis. The red crosses are the locations of the 608 chosen points in X , while the blue circles correspond to the points of X . The first figure is the original signal; the second one gives the approximations coefficients (f 1 ), and the third one the detail coefficients (g 1 ).

Figure 5 .

 5 Figure 5. Results of two steps of the multiresolution for the intertwining wavelets (first graph) and Daubechies 12 wavelets. The three parts of a graph give (g 1 , g 2 , f 2 ).

( a )

 a Original signal. (b) Approximation. Size of f 2 : 206. (c) Detail at scale 1. Size of g 1 : 188 . (d) Detail at scale 2. Size of g 2 : 106.

Figure 7 .

 7 Figure 7. Two steps of the intertwining wavelets multiresolution upsampled to the original graph: (a) f 0 ; (b) R0 R1 f 2 ; (c) Ȓ0 g 1 ; (d) R0 Ȓ1 g 2 .

Figure 8 .

 8 Figure 8. Relative compression error of signal in Figure 1(a) in terms of the percentage of kept normalized detail coefficients. In red, the error using intertwining wavelets. In blue, error using Daubechies12 wavelets.

Figure 9 .

 9 Figure 9. Original signal (blue line) and compressed one (red line) keeping 10% of normalized detail coefficients.

Figure 10 .

 10 Figure 10. Compression with 10% of kept coefficients. The first graph is the compressed signal, the second one is the error with the original one.

Figure 11 .Figure 12 .

 1112 Figure 11. Relative compression error of signal in Figure1(b), in terms of the number of kept coefficients. In red, the error using intertwining wavelets. In blue, error using the spectral graph wavelets pyramidal algorithm.

available at https://lts2.epfl.ch/gsp/

Aknowledgements

Marseille's team thanks Leiden university for its hospitality along many cumulated weeks during which we discovered our wavelets. This was supported by Frank den Hollander's ERC Advanced Grant 267356-VARIS. The four authors are especially grateful to Dominique Benielli (labex Archimède) who made possible the numerical work. L. Avena was partially supported by NWO Gravitation Grant 024.002.003-NETWORKS.

intertwining wavelets. The compressed signal computed with 10% of normalized detail coefficients is shown in Figure 9.

8.3.2. The Minnesota graph. Figure 10 gives the compressed signal computed with 10% of kept coefficients after 3 multiresolution steps.

In addition we compare the intertwining wavelets compression results with those obtained through the spectral graph wavelets pyramidal algorithm of [START_REF] Hammond | Wavelets on graphs via spectral graph theory[END_REF][START_REF] Shuman | A multiscale pyramid transform for graph signals[END_REF]. For this purpose, we used the GSPBox 1 [START_REF] Perraudin | GSPBOX: A toolbox for signal processing on graphs[END_REF]. The main features of this pyramidal algorithm are the following ones:

(1) Subsampling: X is chosen according to the sign of the highest frequency Fourier mode e n-1 .

(2) Weighting: L is computed by the Schur complement followed by a sparsification step.

(3) Approximation coefficients: f (x) = g(tL)(x), where g is a low-pass filter, and t is a positive real number. As in [START_REF] Shuman | A multiscale pyramid transform for graph signals[END_REF], we took g(x) = 1/(1 + x) and t = 2 to analyse the signal of Figure 1

The signal on X is then interpolated on the whole of X : > 0 being fixed, the interpolation is defined as

, where L = Id -L, and L is the Schur complement of L X , X in L . Using Proposition 1, one can see at once that f int | X = f . (5) The error prediction y = f -f int is stored.

At the end of one step, the signal f ∈ R X is encoded by ( f , y) ∈ R X ×X , resulting in a redundant information. Figure 11 presents the relative compression error in terms of the number of kept