Harry Sevi
email: harry.sevi@ens-paris-saclay.fr

Gabriel Rilling
email: gabriel.rilling@cea.fr

Pierre Borgnat
email: pierre.borgnat@ens-lyon.fr

Harmonic analysis on directed graphs and applications: From Fourier analysis to wavelets ‹,‹‹

Keywords: harmonic analysis, graph signal processing, Fourier analysis, wavelets, directed graphs, random walks, semi-supervised learning

We introduce a novel harmonic analysis for functions defined on the vertices of a strongly connected directed graph (digraphs) of which the random walk operator is the cornerstone. As a first step, we consider the set of eigenvectors of the random walk operator as a non-orthogonal Fourier-type basis for functions over digraphs. We find a frequency interpretation by linking the variation of the eigenvectors of the random walk operator obtained from their Dirichlet energy to the real part of their associated eigenvalues. From this Fourier basis, we can proceed further and build multi-scale analyses on digraphs. We propose both a redundant wavelet transform and a decimated wavelet transform as an extension of spectral graph wavelets and diffusion wavelets framework respectively for digraphs. The development of our harmonic analysis on digraphs thus leads us to consider both semi-supervised learning problems and signal modeling problems on graphs applied to digraphs highlighting the efficiency of our framework.

Introduction

In a world where data available for scientific or social purposes accumulates at an exponential pace, managing, exploiting, and analyzing this torrent of data has become one of the challenges of our era. Some of them take the form of graphs, structures that arise in various fields such as neuroscience, the Internet, genomic data, road transportation, or social networks to name a few [START_REF] Barthélemy | Spatial networks[END_REF][START_REF] Barrat | Dynamical processes on complex networks[END_REF]. Thus, there is a need to develop efficient mathematical and computational approaches to process and analyze such graphs and data on graphs. Among the existing methods, a large number involve the (graph) Laplacian [START_REF] Belkin | Laplacian eigenmaps for dimensionality reduction and data representation[END_REF][START_REF] Belkin | Regularization and semi-supervised learning on large graphs[END_REF][START_REF] Giné | Empirical graph laplacian approximation of laplacebeltrami operators: Large sample results[END_REF][START_REF] Coifman | Diffusion maps[END_REF][START_REF] Luxburg | A tutorial on spectral clustering[END_REF]:

• As a fundamental subject in mathematics and physics, the Laplacian is known and used, through its spectral study to extract relevant geometric properties from a manifold [START_REF] Rosenberg | The Laplacian on a Riemannian manifold: an introduction to analysis on manifolds[END_REF] or a graph [START_REF] Chung | Spectral graph theory[END_REF] and plays a role in machine learning applications such as spectral clustering [START_REF] Luxburg | A tutorial on spectral clustering[END_REF] and semi-supervised learning [START_REF] Zhu | Semi-supervised learning with graphs[END_REF].

• In continuous space, the eigenfunctions of the Laplace-Beltrami operator represent the generalization of the Fourier basis on manifolds [START_REF] Rosenberg | The Laplacian on a Riemannian manifold: an introduction to analysis on manifolds[END_REF][START_REF] Lablée | Spectral Theory in Riemannian Geometry[END_REF].

• In discrete space, the eigenvectors of the graph Laplacian form an orthonormal Fourier-type basis for functions on undirected graphs [START_REF] Shuman | The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains[END_REF][START_REF] Ortega | Graph signal processing: Overview, challenges, and applications[END_REF], when undirected graphs are considered as a discretization or sampling of a manifold [START_REF] Giné | Empirical graph laplacian approximation of laplacebeltrami operators: Large sample results[END_REF][START_REF] Coifman | Diffusion maps[END_REF][START_REF] Belkin | Convergence of laplacian eigenmaps[END_REF][START_REF] Hein | Graph laplacians and their convergence on random neighborhood graphs[END_REF].

The eigenvalue associated with each eigenvector is then related to the notion of frequency [START_REF] Shuman | The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains[END_REF][START_REF] Ricaud | Fourier could be a data scientist: From graph fourier transform to signal processing on graphs[END_REF].

Thanks to these properties, the graph Laplacian bridges the gap between spectral graph theory and signal processing through the mathematical framework called "signal processing on graphs" [START_REF] Shuman | The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains[END_REF][START_REF] Ortega | Graph signal processing: Overview, challenges, and applications[END_REF]. This framework aims at extending the concepts of classical signal processing such as filtering or sampling, for functions defined on the vertices of a graph.

First developed in the context of undirected graphs, signal processing on graphs considers the graph Laplacian (combinatorial or normalized) [START_REF] Chung | Spectral graph theory[END_REF] as its core element. The graph Laplacian is a symmetric positive semidefinite operator. By the spectral theorem, the graph Laplacian admits an orthonormal basis of eigenvectors which can be considered as Fourier modes, the corresponding eigenvalues being associated with a notion of frequency.

The purpose of this paper is to provide some answers to the extension of the framework of signal processing on graphs to the case of digraphs. The extension of signal processing on digraphs is of interest because many networks studied in the scientific (e.g. neurosciences) or social (e.g. social networks) fields are directed. Therefore, the analysis of information on digraphs also requires the development of mathematical approaches and therefore an extension of the framework of signal processing on graphs to the case of digraphs.

Formally, the direct use of the core element of graph signal processing, i.e. the graph Laplacian, is no longer appropriate. A directed graph is naturally represented by a nonsymmetric adjacency matrix. It is always possible to naively define a graph Laplacian in the directed case, but the spectral properties associated with this undirected graph Laplacian (e.g. non-negative real eigenvalues and the existence of orthonormal eigenbases) are no longer verified in the case of digraphs. There is no simple consensus on a definition of a Laplacian for directed graphs for which the variation of eigenvectors is linked to a notion of frequency. This is the main challenge of signal processing on digraphs and this paper: which reference operator(s) should a Fourier analysis be built upon that would also generalize Laplacian-based Fourier analysis on undirected graphs?

Hereafter, the random walk operator on graphs [START_REF] Lovász | Random walks on graphs[END_REF][START_REF] Coulhon | Random walks on graphs with regular volume growth[END_REF][START_REF] Aldous | Reversible markov chains and random walks on graphs[END_REF] is proposed as being a suitable reference operator for extending the framework of signal processing on digraphs. Like the graph Laplacian, the random walk operator is associated with the notion of diffusion. However, unlike the graph Laplacian, its definition is as straightforward on directed graphs and undirected ones. Assuming that the random walk operator is diagonalizable, it potentially admits complex conjugate eigenvectors as well as associated complex conjugate eigenvalues. Our framework is based on the fact that the variational analysis of the eigenvectors of the random walk operator through their associated Dirichlet energy [START_REF] Coulhon | Random walks on graphs with regular volume growth[END_REF] is related to the real part of the eigenvalue [START_REF] Sevi | Multiresolution analysis of functions on directed networks[END_REF]. Therefore, the Dirichlet energy of a given eigenvector of the random walk operator on a digraph can be considered as its frequency, using the same analogy as in [START_REF] Shuman | The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains[END_REF], now for digraphs. Provided with this notion of frequency, the eigenvectors of the random walk operator form a suitable Fourier basis on digraphs. It is furthermore consistent with the undirected setting when the Fourier basis is based on the random-walk Laplacian operator (see section 4. [START_REF] Barthélemy | Spatial networks[END_REF].

By now considering the random walk operator as a reference operator and its associated eigenvectors as a Fourier basis for functions defined on digraphs, the last part of our harmonic analysis on digraphs consists in building a multi-resolution analysis for functions defined on the vertices of a digraph. Complex graphs, both directed and undirected, have structures at different scales. Motivated by the efficiency of multi-resolution wavelet analysis in traditional signal processing [START_REF] Mallat | A theory for multiresolution signal decomposition: the wavelet representation[END_REF] and the multi-scale dimension of graph data, a number of multi-resolution graph constructions have emerged in recent years [START_REF] Coifman | Diffusion wavelets[END_REF][START_REF] Maggioni | Diffusion polynomial frames on metric measure spaces[END_REF][START_REF] Hammond | Wavelets on graphs via spectral graph theory[END_REF][START_REF] Rustamov | Wavelets on graphs via deep learning[END_REF][START_REF] Dong | Sparse representation on graphs by tight wavelet frames and applications[END_REF].

Here, we propose a new harmonic analysis built around the random walk operator and whose multi-resolution constructions generalize spectral graph wavelets [START_REF] Maggioni | Diffusion polynomial frames on metric measure spaces[END_REF][START_REF] Hammond | Wavelets on graphs via spectral graph theory[END_REF] and diffusion wavelets [START_REF] Coifman | Diffusion wavelets[END_REF][START_REF] Bremer | Diffusion wavelet packets[END_REF] on digraphs.

Literature review and related work

Wavelets [START_REF] Mallat | A theory for multiresolution signal decomposition: the wavelet representation[END_REF][START_REF] Cohen | Biorthogonal bases of compactly supported wavelets[END_REF][START_REF] Sweldens | The lifting scheme: A custom-design construction of biorthogonal wavelets[END_REF][START_REF] Kingsbury | Complex wavelets for shift invariant analysis and filtering of signals[END_REF][START_REF] Candès | New tight frames of curvelets and optimal representations of objects with piecewise c2 singularities[END_REF][START_REF] Peyré | Orthogonal bandelet bases for geometric images approximation[END_REF] have been thoroughly investigated for over two decades. Their efficiency and success in analyzing functions defined on the real line have led to their generalization for functions defined on higher dimensional spaces such as the sphere [START_REF] Schröder | Spherical wavelets: Efficiently representing functions on the sphere[END_REF][START_REF] Antoine | Wavelets on the 2-sphere: A group-theoretical approach[END_REF][START_REF] Wiaux | Exact reconstruction with directional wavelets on the sphere[END_REF] or other manifolds [START_REF] Coifman | Diffusion wavelets[END_REF][START_REF] Maggioni | Diffusion polynomial frames on metric measure spaces[END_REF][START_REF] Dong | Sparse representation on graphs by tight wavelet frames and applications[END_REF][START_REF] Maggioni | Biorthogonal diffusion wavelets for multiscale representations on manifolds and graphs[END_REF][START_REF] Coulhon | Heat kernel generated frames in the setting of dirichlet spaces[END_REF][START_REF] Geller | Continuous wavelets on compact manifolds[END_REF][START_REF] Rahman | Multiscale representations for manifold-valued data[END_REF]. The development of wavelets on manifolds and the multi-scale features of graphs and data on graphs recently led to the extension of wavelets constructions to the discrete space setting. These wavelet-type constructions include wavelets on unweighted graphs [START_REF] Crovella | Graph wavelets for spatial traffic analysis[END_REF], lifting wavelets [START_REF] Jansen | Multiscale methods for data on graphs and irregular multidimensional situations[END_REF], diffusion wavelets [START_REF] Coifman | Diffusion wavelets[END_REF][START_REF] Bremer | Diffusion wavelet packets[END_REF][START_REF] Maggioni | Biorthogonal diffusion wavelets for multiscale representations on manifolds and graphs[END_REF], diffusion polynomial frames [START_REF] Maggioni | Diffusion polynomial frames on metric measure spaces[END_REF], spectral graph wavelets [START_REF] Hammond | Wavelets on graphs via spectral graph theory[END_REF][START_REF] Tremblay | Graph wavelets for multiscale community mining[END_REF][START_REF] Narang | Perfect reconstruction two-channel wavelet filter banks for graph structured data[END_REF][START_REF] De Loynes | Localized fourier analysis for graph signal processing[END_REF], Haar-like wavelets [START_REF] Szlam | Diffusion-driven multiscale analysis on manifolds and graphs: top-down and bottom-up constructions[END_REF][START_REF] Gavish | Multiscale wavelets on trees, graphs and high dimensional data: Theory and applications to semi supervised learning[END_REF][START_REF] Irion | Hierarchical graph laplacian eigen transforms[END_REF][START_REF] Tremblay | Subgraph-based filterbanks for graph signals[END_REF], average interpolating wavelets [START_REF] Rustamov | Average interpolating wavelets on point clouds and graphs[END_REF], graph wavelets via deep learning [START_REF] Rustamov | Wavelets on graphs via deep learning[END_REF], multi-scale pyramid transform [START_REF] Shuman | A multiscale pyramid transform for graph signals[END_REF], tight wavelet frames on graphs [START_REF] Dong | Sparse representation on graphs by tight wavelet frames and applications[END_REF][START_REF] Göbel | Construction of tight frames on graphs and application to denoising[END_REF][START_REF] Behjat | Signal-adapted tight frames on graphs[END_REF] and intertwining wavelets on graphs via random forests [START_REF] Avena | Intertwining wavelets or multiresolution analysis on graphs through random forests[END_REF]. The wavelet constructions mentioned above were all developed on undirected graphs. In the present work, we intend to construct a harmonic analysis on digraphs. The first stage of the proposed harmonic analysis is the development of a Fourier analysis.

In recent years, Sandryhaila and Moura generalized some fundamental concepts of traditional signal processing such as filtering to digraphs using the adjacency matrix as the central component of their framework [START_REF] Sandryhaila | Discrete signal processing on graphs[END_REF]. They also proposed the generalized eigenvectors of the adjacency matrix obtained by Jordan decomposition as a Fourier-type basis on digraphs [START_REF] Sandryhaila | Discrete signal processing on graphs: Frequency analysis[END_REF]. To the best of our knowledge, no wavelet design has been proposed in this framework, except the development of the filterbank approach generalized in [START_REF] Tay | Bipartite graph filter banks: Polyphase analysis and generalization[END_REF] which can be applied to bipartite directed graphs and uses the adjacency matrix for a polyphase representation of the filters.

More recently, new Fourier-type bases have been proposed which are based on a novel measure of graph directed variation (GDV) for functions defined on directed graphs, which is defined as the Lovàsz extension of the graph cut size [START_REF] Sardellitti | On the Graph Fourier Transform for Directed Graphs[END_REF][START_REF] Shafipour | A Directed Graph Fourier Transform With Spread Frequency Components[END_REF]. In both propositions, Fourier modes are defined as a set of orthonormal functions that approximately solve a non-convex optimization problem involving the GDVs of the Fourier modes. The two approaches differ insofar as the minimized objective in [START_REF] Sardellitti | On the Graph Fourier Transform for Directed Graphs[END_REF] is simply the sum of the GDVs of the Fourier modes whereas in [START_REF] Shafipour | A Directed Graph Fourier Transform With Spread Frequency Components[END_REF] the objective is designed to spread the individual GDVs. It is worth stressing that because of the non-convexity of the objective, the Fourier modes are not properly defined by the optimization problems alone and depend on the specific algorithms used to approximately solve the problems. In both cases, the Fourier modes could also be defined on undirected graphs but they would not coincide with standard approaches based on Laplacian operators. The reasons are that the Fourier modes are not eigenvectors of any given operator and that the GDV on undirected graphs differs from the common quadratic variation measures used on undirected graphs. Unlike these approaches, the Fourier basis that we propose is generally not orthogonal but it does generalize the Laplacian-based Fourier basis for undirected graphs. Furthermore, we associate a notion of frequency to the Fourier modes, eigenvectors of the random walk operator, via their Dirichlet energy which is conveniently related to the corresponding eigenvalues. Most of the Fourier-like bases approaches on directed graphs have been also discussed in [START_REF] Marques | Signal processing on directed graphs: The role of edge directionality when processing and learning from network data[END_REF].

Mhaskar has proposed the first wavelet-type frame construction for functions defined on digraphs [START_REF] Mhaskar | A unified framework for harmonic analysis of functions on directed graphs and changing data[END_REF] as an extension of the diffusion polynomial frames [START_REF] Maggioni | Diffusion polynomial frames on metric measure spaces[END_REF]. Although our redundant wavelet construction on digraphs is inspired in part by that of polynomial diffusion frames, it differs mainly in the choice of the Fourier basis and the considered operator. Mhaskar used in [START_REF] Mhaskar | A unified framework for harmonic analysis of functions on directed graphs and changing data[END_REF] the set of left-singular vectors of the weighted adjacency matrix of a digraph as a Fourier basis while we propose the set of eigenvectors of the random walk operator. Recently, adaptive directional Haar tight framelets have been proposed and can deal with digraph signal representations [START_REF] Xiao | Adaptive directional haar tight framelets on bounded domains for digraph signal representations[END_REF] inspired by a previous harmonic analysis on digraphs [START_REF] Chui | Representation of functions on big data associated with directed graphs[END_REF].

Finally, Furutani et al. [START_REF] Furutani | Graph signal processing for directed graphs based on the hermitian laplacian[END_REF] proposed an extension of graph signal processing to digraphs using an operator called the Hermitian Laplacian or Magnetic Laplacian [START_REF] De Verdière | Magnetic interpretation of the nodal defect on graphs[END_REF][START_REF] Berkolaiko | Nodal count of graph eigenfunctions via magnetic perturbation[END_REF] as their reference operator. The main properties of this latter operator are to be hermitian and preserve the directionality of the digraph. The magnetic Laplacian has been highlighted through machine learning applications such that community detection on digraphs [START_REF] Fanuel | Magnetic eigenmaps for community detection in directed networks[END_REF] and graph neural networks on digraphs [START_REF] Zhang | Magnet: A magnetic neural network for directed graphs[END_REF]. Other operators based on digraphs such that the Hermitian adjacency matrix [START_REF] Guo | Hermitian adjacency matrix of digraphs and mixed graphs[END_REF] have been investigated through clustering on digraphs [START_REF] Cucuringu | Hermitian matrices for clustering directed graphs: insights and applications[END_REF][START_REF] Laenen | Higher-order spectral clustering of directed graphs[END_REF], but no Fourier analyses of these operators have been carried out.

Contributions

The contributions of this article are the following:

1. Construction of Fourier bases on digraphs. We propose the random walk operator associated with a random walk on a digraph as a reference operator. We propose the eigenvectors of the random walk operator as a Fourier basis on digraphs and determine a variational characterization of the eigenvectors of the random walk operator -see Proposition 4.1.

Construction of multi-resolution analyses on digraphs.

We propose an overcomplete spectral wavelet transform on digraphs in Section 5.1. This construction extends the framework of spectral wavelets on undirected graphs [START_REF] Hammond | Wavelets on graphs via spectral graph theory[END_REF] and diffusion polynomial frames [START_REF] Maggioni | Diffusion polynomial frames on metric measure spaces[END_REF]. We also propose a critically sampled wavelet transform in Section 5.2 generalizing the framework of diffusion wavelets [START_REF] Coifman | Diffusion wavelets[END_REF][START_REF] Bremer | Diffusion wavelet packets[END_REF] to the directed case.

3. Efficiency of the theoretical framework through applications. The development of our harmonic analysis on digraphs leads us to consider semi-supervised learning problems with ℓ 2 regularization in Section 6.1 and ℓ 1 regularization in section 6.3.3 and signal modeling by filtering on digraphs in Section 6.2 highlighting the effectiveness of our theoretical framework.

Outline of the paper

This paper is structured as follows. We introduce the essential aspects of graph theory and review the foundations of graph signal processing in Section 2. We present operators defined on digraphs built from the random walk operator and their properties in Section 3 and propose a new class of operators on digraphs expressed as the convex combination between the random walk operator and its time-reversed version. Section 4.1 is dedicated to the presentation of a Fourier transform on digraphs as a set of eigenvectors of the random walk operator, and we propose a Fourier-type analysis for functions defined on digraphs by studying the variation of the random walk's eigenvectors. Section 5 is dedicated to multiresolution analyses on digraphs. We present a redundant wavelets construction, as well as a critically sampled wavelets construction with the random walk operator as a reference operator. To illustrate our Fourier-type analysis, we study the machine learning problems in Section 6 such as semi-supervised learning and signal modeling on digraphs, and show through numerical experiments the efficiency of our framework compared to the existing approaches and we illustrate these multiresolution analyses through some examples. We conclude in Section 7.

Preliminaries

Fundamentals of signal processing on graphs

Signal processing on graphs [START_REF] Shuman | The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains[END_REF][START_REF] Ortega | Graph signal processing: Overview, challenges, and applications[END_REF][START_REF] Sandryhaila | Discrete signal processing on graphs[END_REF] is a mathematical framework in which the core concepts of harmonic analysis are generalized to functions defined on the vertices of a graph. In this section, we introduce the essential aspects of the signal processing on graphs framework and give some additional remarks.

Graph theory setup

Let G " pV, E, wq be a weighted directed graph (digraph) where V " tv 1 , . . . , v N u is a finite set of vertices, E Ď V ˆV is a set of directed edges. Each edge e ij " pv i , v j q is directed and represents a link from vertex v i to vertex v j . The edge weight function w : V ˆV Ñ R `associates a nonnegative real value to every vertex pair: wpv i , v j q ą 0 if e ij P E, otherwise wpv i , v j q " 0. We denote by |V| " N , the cardinality of the vertex set V, that is the total number of vertices in V. A weighted digraph G can be entirely represented by its weighted adjacency matrix W " tw ij u 1ďi,jďN P R N ˆN `where w ij " wpv i , v j q are the weights associated to the respective edges e ij . We define respectively the out-degree and the in-degree of a vertex v i P V by d ì " ř N j"1 w ij and d í "

ř N j"1 w ji . For the sake of simplicity, we refer to the out-degree d ì of a vertex v i P V as its degree that we denote by d i . A digraph is called strongly connected if there is a path in each direction between each pair of vertices of the graph.

A limitation of our framework is that it only applies to the case of weighted digraphs with nonnegative weights, that is W " tw ij u 1ďi,jďN P R N ˆN `. Other approaches exist for the case of digraph with positive and negative weights but they are not compatible with the one presented here [START_REF] Sandryhaila | Discrete signal processing on graphs[END_REF][START_REF] Sandryhaila | Discrete signal processing on graphs: Frequency analysis[END_REF][START_REF] Bauer | Normalized graph laplacians for directed graphs[END_REF].

We assume through the theoretical sections of this paper that the digraph G is strongly connected and the random walk operators are diagonalizable. These are two necessary conditions for our theoretical framework. If the digraph is not strongly connected, our framework may still apply by modifying the digraph to make it strongly connected. This can be achieved using e.g. the Pagerank method [START_REF] Langville | Google's PageRank and beyond: The science of search engine rankings[END_REF], which will be used in Section 6.2. Since diagonalizable matrices form a dense set, the case of a non-diagonalizable random walk operator is a priori very rare in practice and typically never occurs in numerical eigendecompositions.

Graph signals

Let f : V Þ Ñ C be a function defined on the vertex set V of a given digraph G. We define a graph signal f as a column vector representation of the function f applied at each node v i P V, i.e.

f " rf pv 1 q, . . . , f pv N qs J P C N .

We assume throughout this paper that the graph signals are defined in ℓ 2 pV, νq which is the Hilbert space of functions defined over the vertex set V of G endowed with the inner product associated with ν : V Þ Ñ R `, a positive arbitrary vertex measure xf , gy ν "

ÿ xPV f pxqgpxqνpxq, (1)
for all functions f , g and where f pxq denotes the complex conjugate of f pxq. Note that we can easily generalize to any other ℓ p norm.

Linear filters on graphs

A linear operator on graph is represented by a matrix H P C N ˆN that acts on a graph signal input s P C N and produces a graph signal output s P C N according to the matrix vector product s " Hs.

Within the framework of graph signal processing and inspired by conventional signal processing, we are mainly interested in graph filters. As in [START_REF] Sandryhaila | Discrete signal processing on graphs[END_REF], a graph filter is supposed linear and "shift-invariant" (i.e., LSI) in the sense that it commutes with a chosen reference operator R. In [START_REF] Sandryhaila | Discrete signal processing on graphs[END_REF], the adjacency matrix is chosen as reference operator, and it has been discussed since (e.g. [START_REF] Shuman | The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains[END_REF][START_REF] Tremblay | Design of graph filters and filterbanks[END_REF]) that other choices are possible. Hence, we consider any reference operator R on the graph. Definition 2.1 ( [START_REF] Sandryhaila | Discrete signal processing on graphs[END_REF][START_REF] Tremblay | Design of graph filters and filterbanks[END_REF]). A graph filter is a linear shift-invariant (LSI) operator with matrix representation H so that it commutes with the chosen reference operator R, i.e.

HR " RH.

(

) 2
Specifically, we can characterize graph filters as being expressed as a polynomial of the reference operator R, thanks to the following Theorem 1 from [START_REF] Sandryhaila | Discrete signal processing on graphs[END_REF]. Note that there are some specific exceptions (when the characteristic and minimum polynomials of R are not equal).

Theorem 2.1. [START_REF] Sandryhaila | Discrete signal processing on graphs[END_REF] Let R be a reference operator on a digraph G. Let us assume that the characteristic and minimum polynomials of R are equal, i.e. p R pxq " m R pxq. Then a graph filter H is linear and shift-invariant if and only if it is expressed as a polynomial of R H " T ÿ t"0 h t R t , h t P C, @t " 0, . . . , T.

(

As an example, considering a diagonalizable R whose eigenspaces are of dimension 1, then any graph filter H (linear and commuting with R) has an expression as a finite polynomial of R (as there is only N degrees of freedom to define its effect in the eigenspaces of R). Now, it is often useful to consider graph filters without reducing them to their polynomial form, and a filter can be also seen as any function applied in the eigenspace of R to the eigenvalues of R, 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]. More general situations (with R non diagonalizable) are discussed in the Section 2 of [START_REF] Tremblay | Design of graph filters and filterbanks[END_REF].

Remark 2.1. This notion of invariance with respect to a reference operator on graphs may seem rather restrictive insofar as many diagonalizable reference operators on graphs may have eigenspaces of dimension larger than 1. However, it could also be argued that the set of all commutating filters may include unexpected elements and that the set of polynomial filters may be a more relevant generalization of the LSI filters in classical signal processing. As an example, let us consider the undirected cycle graph. Its (combinatorial) Laplacian is the circulant matrix with three main diagonals

‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ '
, for which the discrete Fourier basis is an eigenbasis. However, most eigenspaces of the Laplacian matrix have dimension equal to 2 and are spanned by pairs of complex conjugate exponentials v j rks " e ikj{N and v ´j rks " e ´ikj{N . The operator that turns a complex exponential into its conjugate, swapping v j with v ´j (i.e. swapping positive and negative frequencies), is usually not considered a LSI filter in classical signal processing. It commutes with the Laplacian but cannot be expressed as a polynomial of the Laplacian. In contrast, filters defined as polynomials of the Laplacian correspond to symmetric filters in classical signal processing, which seems more appropriate as a generic filter class on the undirected cycle graph.

A graph filter can also be characterized by its eigenvalues on each eigenspace of R. Let tE j u m j"1 denote the eigenspaces of R and tE j u m j"1 the corresponding spectral projectors, characterized by E 2

i " E i and E i E j " δ ij E j . A graph filter with respect to a diagonalizable operator R is a linear combination of the spectral projectors

H " m ÿ j"1 γ j E j , γ j P C, E j P C N ˆN , @j " 1, . . . , m, (4)
where γ j is the eigenvalue of H associated to the eigenspace E j .

Canonical operators on digraphs

In this section, we introduce the fundamental linear operators and their properties in order to develop a harmonic analysis on digraphs.

Random walk operators on graphs

We define a random walk on a directed graph as follows.

Definition 3.1 ( [START_REF] Lovász | Random walks on graphs[END_REF][START_REF] Aldous | Reversible markov chains and random walks on graphs[END_REF]). A random walk on a strongly connected digraph G " pV, E, wq is a homogeneous Markov chain X " pX n q ně0 with a finite state space V and whose transition probabilities are proportional to the edges' weights. The entries of the transition matrix P, written as ppx, yq, @x, y P V are defined as ppx, yq " PpX n`1 " y|X n " xq " wpx, yq ř y wpx, yq .

A transition matrix must have non-negative entries and rows that sum up to one. This is guaranteed by the fact that the wpx, yq are non-negative and that the denominators are always non-zero thanks to the graph being strongly connected. From a graph theory perspective, the transition matrix P P R N ˆN is equal to

P " D ´1W,
where D " diagtd 1 , . . . , d n u is the diagonal matrix of the out-degrees of the vertices and W is the weighted adjacency matrix. We call P n , the n-step transition matrix whose entries written as p n px, yq, @x, y P V are defined as p n px, yq " PpX n " y|X 0 " xq.

For any graph signal f , measure µ : V Þ Ñ R `and states x, y P V, we adopt the following conventions [START_REF] Brémaud | Markov chains: Gibbs fields[END_REF]:

Pf pxq " ÿ yPV ppx, yqf pyq, (5)
µPpyq "

ÿ xPV µpxqppx, yq. (6)
where f and µ are respectively the column and row vector representations of f and µ.

Let us consider the transition matrix P P R N ˆN and a Kronecker delta function δ k at vertex k P V, as a graph signal. If P acts on δ k , i.e. Pδ k , the mass at vertex v k in δ k propagates towards vertices where ppv i , v k q ‰ 0, that is vertices with edges pointing to v k . In other words, the mass propagates towards parent vertices and not children, which may seem counter-intuitive. Moreover, when δ k is applied to the left of P, i.e. δ k P, the mass at vertex v k in δ k propagates from the vertex v k towards vertices where ppv k , v i q ‰ 0, that is edges edges coming from v k . In other words, the mass propagates from the parent vertex towards the children's vertices. When applied to an arbitrary signal x, the random walk operator P thus behaves as a local averaging of the children's vertices. Given the importance of local averaging in signal processing, this suggests that the random walk operator P may be natural as a core element of graph signal processing. We define the irreducibility of a random walk on a graph as follows.

Definition 3.2. A random walk X is said to be irreducible if for any x, y P V, the probability from x to reach y is strictly positive, in other words: @x, y P V, Dm ă 8 : PpX n`m " y|X n " xq ą 0.

Irreducibility is equivalent to the strong connectivity of G, that is there is a (directed) path from any vertex v i P V to any vertex v j P V. Remark 3.1. In the case of connected undirected graphs, the graph is always strongly connected and the associated random walk is always irreducible.

We now recall the notions of recurrent and transient state. Definition 3.3 ([75]). Let X " pX n q ně0 be a random walk on G. We define the return time to state x P V as T i " inftn ě 1, X n " iu.

A vertex x P V is called recurrent if PpT i |X 0 " iq " 1, and otherwise is called transient.

Proposition 3.1 ([75]

). An irreducible random walk X with a finite state space is always recurrent: all states are recurrent.

Remark 3.2. An irreducible random walk X with finite state space is positive recurrent. Consequently, the random walk X admits a unique stationary distribution [START_REF] Brémaud | Markov chains: Gibbs fields[END_REF].

We define the periodicity of a random walk on a graph as follows.

Definition 3.4. Let X be a random walk on G. The period of a vertex x P V is:

ϱpxq " gcdtn P N `: p pnq px, xq ą 0u.

Thus starting in x, the Markov chain can return to x only at multiples of the period ϱ. The state x is aperiodic if ϱpxq " 1 and periodic is ϱpxq ą 1. Having defined the notions of irreducibility and periodicity, we are able to define the notion of ergodicity of a random walk. Definition 3.5. Let X be a random walk on G. The random walk X is ergodic if it is irreducible and aperiodic.

We set out the proposition for the stationary distribution of an ergodic random walk on G.

Proposition 3.2 ([76]

). Let G be a digraph with f inite state space V. If a random walk X with its transition matrix P is ergodic, i.e. irreducible and aperiodic, the measures P n px, .q converge towards the row vector π " rπpv 1 q, ¨¨¨, πpv N qs P R N `as n Ñ 8, i.e. the unique stationary distribution. In particular, πP " π, with :

N ÿ i"1
πpv i q " 1, πpv i q ě 0.

Reversibility

Given the discrete-time Markov chain X " pX n q ně0 and M ą 0 a finite time horizon, we define the time reversed Markov chain as X ˚" pX n q " pX M ´nq for n " 0, ¨¨¨, M . We denote by P ˚" tp ˚px, yqu x,yPV the transition matrix associated with the Markov chain X ˚. It verifies p ˚px, yq " PpX n`1 " y|X n " xq " PpX M ´n´1 " yq PpX M ´n " xq ppy, xq, @x, y P V.

If we assume that X is ergodic with stationary distribution π, the time reversed Markov chain X ˚is also ergodic with stationary distribution π and the entries p ˚px, yq P P ˚are p ˚px, yq " πpyq πpxq ppy, xq, @x, y P V or in its transition matrix version

P ˚" Π ´1P J Π,
where Π " diagtπpv 1 q, ¨¨¨, πpv N qu is the diagonal matrix of the stationary distribution.

A random walk X with stationary distribution π and transition matrix P is called reversible if πpxqppx, yq " πpyqppy, xq, @x, y P V.

Let us use the function space ℓ 2 pV, πq (endowed with its inner product as in Eq. ( 1)) for the stationary distribution π.

In this space, P ˚is the adjoint of P, that is xf , Pgy π " xP ˚f , gy π , for all f , g. The ergodic random walk X with its transition matrix P is reversible if and only if we have the following relation P " P ˚.

Remark 3.3. Ergodic random walks on finite undirected graphs are reversible. That means that the transition matrix associated to the time reversed random walk X ˚, P ˚is equal to the transition matrix P associated to the original ergodic random walk X namely P ˚" P.

Remark 3.4. In the undirected setting, the stationary distribution π admits a closed form expression. Indeed, on a given weighted undirected graph G " pV, E, wq represented by its symmetric adjacency matrix W " tw xy u x,yPV P R N ˆN `and the degree of a vertex x P V is cpxq " ř yPV w xy . As a result, the associated random walk is reversible with stationary distribution π defined by πpxq " cpxq{c G where c G " ř xPV cpxq. However, in the directed setting, the stationary distribution π does not admit an analytical form. In order to calculate it, we use iterative methods such as the power iteration method [START_REF] Golub | Matrix computations[END_REF] or Markov Chain Monte Carlo(MCMC) methods [START_REF] Atchadé | Adaptive markov chain monte carlo: theory and methods[END_REF].

Remark 3.5. Ergodic random walks defined on digraphs are typically non-reversible. Nevertheless, reversible ergodic random walks may be constructed on digraphs by modifying the original non-reversible ergodic random walk. Such a modification is discussed in Section 3.2.2.

Eigenvalue distribution

By the Perron-Frobenius theorem [START_REF] Brémaud | Markov chains: Gibbs fields[END_REF], if the random walk X is ergodic with a diagonalizable transition matrix P, the diagonalization of P admits a simple dominant eigenvalue λ max " 1. The other eigenvalues tλ ‰ 1u satisfy |λ| ă 1, which means that all eigenvalues different of λ max lie within the unit circle.

Random walk generalizations

Given a random walk X on a digraph G with transition matrix P, we are able to build new types of random walks based on P with various purposes.

Lazy random walks

The periodicity of a random walk X can be overcome by considering the lazy random walk version of X . The transition matrix P associated to the lazy random walk X , based upon P is expressed as P " I `P 2 .

The lazy random walk can be seen as the random walk on a modified graph G where an edge from each vertex of G to itself is added with a weight equal to the vertex's degree in G. The lazy random walk X is always aperiodic, hence when G is strongly connected, X is also ergodic with a stationary measure π verifying π " π P " πP.

More generally, let us define P, the class of transition matrices associated with generalized lazy random walks by P "

"

Pγ : Pγ " p1 ´γqP `γI ˇˇˇγ P r0, 1q

* .

We note that the elements Pγ P P share the same eigenspaces as P, hence the graph filters with respect to P are also graph filters with respect to P.

Additive reversibilization

From a non-reversible ergodic random walk X with transition matrix P and unique stationary measure π, an additive reversibilization of X can be constructed [START_REF] Fill | Eigenvalue bounds on convergence to stationarity for nonreversible markov chains, with an application to the exclusion process[END_REF], denoted by X , whose transition matrix is the average between P and its time reversed P ˚:

P " P `P2 . (7)
X is a reversible random walk with the same unique stationary distribution π. More generally, let us define P the class of convex combinations between the random walk matrix P and its time reversed version P ˚as P "

"

Pα : Pα " p1 ´αqP `αP ˚ˇα P r0, 1s

* . (8)
We note that the elements Pα P P share the same stationary distribution π but do not have the same eigenspaces. We also note that all the random walks associated to transition matrices Pα P P are non-reversible except for α " 1{2, namely P1{2 " P.

Remark 3.6. Given an ergodic random walk pX , P, πq, we note that all lazy and reversibilized versions, respectively Pγ P P and Pα P P share the same unique stationary distribution π.

Remark 3.7. As a generalization of random walks, we may also consider the multiplicative reversibilization of the ergodic non-reversible random walk P, that is PP ˚which is also a reversible Markov chain with stationary distribution π. It is useful e.g. for defining the convergence bounds of non-reversible Markov chains [START_REF] Fill | Eigenvalue bounds on convergence to stationarity for nonreversible markov chains, with an application to the exclusion process[END_REF].

Laplacians on digraphs

In this section, we introduce several definitions of Laplacians on digraphs, which are extended from the undirected case [START_REF] Marques | Signal processing on directed graphs: The role of edge directionality when processing and learning from network data[END_REF][START_REF] Anis | Efficient sampling set selection for bandlimited graph signals using graph spectral proxies[END_REF][START_REF] Girault | Irregularity-aware graph fourier transforms[END_REF].

Normalized directed graph Laplacian

The normalized Laplacian on a digraph G is expressed in terms of the transition matrix P of an ergodic random walk X on G and is defined as follows. Definition 3.6 ([82,[START_REF] Butler | Interlacing for weighted graphs using the normalized laplacian[END_REF]). Let G " pV, Eq be a digraph with |V| " N . Let X be an ergodic random walk on G with transition matrix P and unique stationary distribution π. The normalized Laplacian on G is defined by

L " I ´Π1{2 PΠ ´1{2 `Π´1{2 P J Π 1{2 2 , ( 9
)
where I is the identity matrix and Π " diagtπpv 1 q, . . . , πpv N qu is the diagonal matrix of the stationary distribution.

Random walk Laplacian

Another definition is the random walk Laplacian L RW , defined as

L RW " I ´P. (10)
The normalized Laplacian on digraphs L is connected to the random walk Laplacian through the following relation

L " Π 1{2 L RW Π ´1{2 .
More generally, we define the random walk Laplacian on digraphs associated to a transition matrix Pα P P by

L RW,α " I ´P α `P α 2 .
In the following proposition, we show that the random walk Laplacian L RW,α is equal to L RW .

Proposition 3.3. For any Pα P P, we have

L RW,α " L RW , @α P r0, 1s. (11)
Proof.

L RW,α " I ´P α `P α 2 " I ´αP `p1 ´αqP ˚`αP ˚`p1 ´αqP 2 " I ´αpP `P˚q `p1 ´αqpP `P˚q 2 " I ´P `P2 " L RW .

Combinatorial Laplacian

Finally, we also introduce the combinatorial Laplacian L, defined as

L " Π ´ΠP `PJ Π 2 . ( 12
)
The latter is related to the random walk Laplacian through

L " Π L RW .
Remark 3.8. As mentioned in remark 3.3, the transition matrix of a random walk on an undirected graph is equal to the transition matrix of its time reversed Markov chain. Furthermore, the stationary distribution π admits a closed form as we mention in remark 3.4. As a result, the definitions eqs. ( 9) and (10), generalize the usual definitions for undirected graphs [START_REF] Chung | Spectral graph theory[END_REF]. The combinatorial Laplacian on digraphs (12) (thereafter denoted L directed ) generalizes the one on undirected graphs up to a multiplicative factor:

L undirected " D ´W " p ř i d i qL directed .

Canonical operators on digraphs and Hilbert spaces

Let G " pV, Eq be a directed graph, ℓ 2 pVq and ℓ 2 pV, πq be the Hilbert spaces associated respectively with the counting measure and with the stationary measure π of the ergodic random walk X on G. Let φ : ℓ 2 pVq Ñ ℓ 2 pV, πq be the linear mapping defined as follows

φ : f Þ Ñ Π ´1{2 f , @f P ℓ 2 pVq. ( 13
)
Definition 3.7 ( [START_REF] Kubrusly | Hilbert space operators[END_REF]). Let H, K be two Hilbert spaces. A linear transformation V : H Þ Ñ K is an isometry if and only if }Vpx ´yq} K " }x ´y} H , @x, y P H.

Proposition 3.4. The linear transformation φ : ℓ 2 pVq Þ Ñ ℓ 2 pV, πq is an isometry.

Proof. Given g, h P ℓ 2 pVq, and the linear transformation φ, we have

xφpgq, φphqy π " xΠ ´1{2 g, Π ´1{2 hy π " xg, hy.
Let us introduce the operator T P ℓ 2 pVq defined as

T " Π 1{2 PΠ ´1{2 . ( 14
)
Definition 3.8 ( [START_REF] Kubrusly | Hilbert space operators[END_REF]). Let H, K bet two Hilbert spaces. An invertible bounded linear transformation

V : H Þ Ñ K intertwines an operator M : H Þ Ñ H to an operator S : K Þ Ñ K if VM " SV.
The operators M and S are called similar.

Proposition 3.5. The linear operator φ intertwines the operator T on ℓ 2 pVq to the operator P on ℓ 2 pV, πq.

Proof. φpTq " Π ´1{2 T " Π ´1{2 Π 1{2 PΠ ´1{2 " PΠ ´1{2 .
Hence the operators T and P are similar with respect to φ. Identically, the random walk Laplacian L RW and the normalized Laplacian L are also similar with respect to φ. This suggests that P and L RW are best viewed as operators on ℓ 2 pV, πq while T and L are essentially the same operators on ℓ 2 pVq.

Fourier analysis on digraphs

In this section, we propose a new analysis on digraphs that is different from the existing approaches on digraphs [START_REF] Sandryhaila | Discrete signal processing on graphs: Frequency analysis[END_REF][START_REF] Sardellitti | On the Graph Fourier Transform for Directed Graphs[END_REF][START_REF] Shafipour | A Directed Graph Fourier Transform With Spread Frequency Components[END_REF][START_REF] Mhaskar | A unified framework for harmonic analysis of functions on directed graphs and changing data[END_REF]. It is based on the study of variations in the eigenvectors of the random walk operator. Before discussing it in greater detail, we introduce the elements that allow us to study the variation of signals on digraphs.

Digraph Fourier transform

Let P be the transition matrix of the ergodic random walk X . We assume P is diagonalizable, that is P admits an eigenvalue decomposition

P " ΞΘΞ ´1,
where Ξ " rξ 1 , ¨¨¨, ξ N s is an eigenbasis with elements ξ j , j " 1, . . . , N and Θ " diagtϑ 1 , ¨¨¨, ϑ N u is the corresponding diagonal eigenvalue matrix, not necessarily distinct. Given a graph signal s, its digraph Fourier transform denoted by ŝ " rŝ 1 , . . . , ŝN s is ŝ " Ξ ´1s.

The values tŝ j u N j"1 characterize the content of the graph signal s in the graph Fourier domain.

The graph Fourier transform is not uniquely defined as it depends on the choice of eigenbasis Ξ. First all eigenvectors ξ j are defined up to a multiplicative scalar. This can be easily resolved by normalizing the vectors and e.g. making the first non-zero coefficient of each eigenvector real and positive. Second, whenever there are eigenspaces with dimension greater than 1, there are infinite choices for the orientations of the basis vectors within such eigenspaces. There is generally no solution to this second issue. This implies that the values of individual Fourier coefficients ŝj within such eigenspaces are rather meaningless. However, it only impacts cases with eigenspaces of dimension greater than 1, which are not that common in applications.

The inverse digraph Fourier transform corresponding to an eigenbasis Ξ is given by s " Ξŝ.

It reconstructs the original signal from its frequency contents by forming a linear combination of eigenvectors weighted by the Fourier coefficients. More generally, any Pα P P admits an eigenvalue decomposition

Pα " Ξ α Θ α Ξ ´1 α .
Consequently, one can build an infinity of Fourier-type bases tΞ α u αPr0,1s for functions defined on digraphs. Among these, the case α " 1{2 is particularly interesting insofar as P1{2 " P is self-adjoint in ℓ 2 pV, πq. This implies that there is an orthonormal eigenbasis Ξ 1{2 of ℓ 2 pV, πq which yields the following theorem.

Theorem 4.1 (Generalized Parseval's Theorem). Given an eigenbasis Ξ 1{2 of P that is orthonormal in ℓ 2 pV, πq, the corresponding Fourier transform is an isometric operator from ℓ 2 pV, πq to ℓ 2 pt1, . . . , N uq xΞ ´1 1{2 x, Ξ ´1 1{2 yy " xx, yy π , @x, y P ℓ 2 pV, πq.

Proof.

Ξ J 1{2 ΠΞ 1{2 " I ùñ xΞ ´1 1{2 x, Ξ ´1 1{2 yy " x J Ξ ´J 1{2 Ξ ´1 1{2 y " x J Πy " xx, yy π
Another interesting property of P is that its eigenspaces are the same as those of the random walk Laplacian L RW [START_REF] Rosenberg | The Laplacian on a Riemannian manifold: an introduction to analysis on manifolds[END_REF].

Remark 4.1. In the undirected setting, a suitable Fourier-type basis is usually defined as an orthonormal eigenbasis of the combinatorial Laplacian L or its normalized counterpart L [START_REF] Chung | Spectral graph theory[END_REF]. These Fourier-type bases are both orthonormal for ℓ 2 pVq. Another less common choice is a Fourier-type basis based upon the random walk Laplacian L RW " I ´P. As in the directed case with P, it leads to a Fourier basis that is orthonormal in ℓ 2 pV, πq, where π is proportional to the degrees of the vertices (see remark 3.4). Some properties of the different definitions are investigated in [START_REF] Girault | Irregularity-aware graph fourier transforms[END_REF].

.

Regularity of signals on graphs

The behavior of a graph signal over a graph can be analyzed by measuring the variation or regularity. We first define the length of the graph gradient at a given vertex as follows.

Definition 4.1 ([20]

). Let G " pV, Eq be a digraph and q : E Ñ R `be a positive measure defined on the edge set E. The length of the graph gradient of a graph signal f at vertex x P V on an arbitrary graph G under the measure ν is

|∇f pxq| " ˆ1 2 ÿ yPV,px,yqPE |f pxq ´f pyq| 2 qpx, yq ˙1{2 .
Intuitively, the length of the graph gradient measures the variation of a graph signal around a given vertex. We now introduce the Dirichlet energy as a measure of the variation of a signal over a strongly connected graph.

Definition 4.2 ([76]

). The Dirichlet energy of a graph signal f associated to the ergodic random walk X with transition matrix P and stationary distribution π on a strongly connected graph G is

D 2 π,P pf q " 1 2 ÿ px,yqPE πpxqppx, yq|f pxq ´f pyq| 2 , (16)
" xf , L RW f y π .

For all Pα P P, the associated Dirichlet energy is the same:

D 2 π, Pα pf q " D 2 π,P pf q " xf , L RW f y π .
We also introduce the Rayleigh quotient of a graph signal f associated to the Dirichlet energy D 2 π,P pf q as R π,P pf q " D 2 π,P pf q }f } 2 π .

For any Pα P P, we have R π,P pf q " R π, Pα pf q.

Frequency analysis on digraphs

The following proposition settles a key connection between the regularity of the eigenvectors of the transition matrix P of the ergodic random walk X and their associated eigenvalues.

Proposition 4.1. Let ξ P C N be an eigenvector of a transition matrix P of an ergodic random walk X , with stationary distribution π, associated to the eigenvalue ϑ P C. The Rayleigh quotient of ξ is given by R π,P pξq " 1 ´Repϑq, where Repϑq denotes the real part of ϑ P C.

Proof. R π,P pξq " 1 }ξ} 2 π ˆxξ, ξy π ´1 2 xξ, Pξy π ´1 2 xξ, P ˚ξy π " 1 }ξ} 2 π ˆ}ξ} 2 π ´ϑ `θ 2 }ξ} 2 π Ṙπ,P pξq " " 1 ´Repϑq ‰ .
The latter proposition thus indicates that the variation of any eigenvector ξ of P, as described by its Rayleigh quotient, is associated to the real part of its respective eigenvalue ϑ. More generally, the variation of any eigenvector ξ α of Pα P P is characterized exactly in the same manner as in Proposition 4.1, i.e. R π,P pξ α q " 1 ´Repϑ α q.

Therefore, we are now able to associate to each eigenvector ξ of P and generally to each eigenvector ξ α of a given Pα P P, a value ω characterizing its variation that we call frequency expressed intuitively as ω " 1 ´Repϑq, ω P r0, 2s.

Remark 4.2. In the undirected setting, the random walk is reversible such that the Dirichlet energy associated to the random walk operator P is

D π,P pf q " xf , L RW f y π ,
where L RW " I ´P. The random walk operator P is self-adjoint in ℓ 2 pV, πq. As a result, given an eigenvector ξ of P associated to an eigenvalue ϑ, the associated Rayleigh quotient is R π,P pξq " 1 ´ϑ.

Therefore, the variation of the eigenvectors of P is directly related to their respective eigenvalues. A similar result holds when the Fourier basis is based on the combinatorial Laplacian L " D ´W thanks to the following identity

D π,P pf q " xf , L RW f y π " xf , Lf y.
As L is a symmetric semi-definite operator, its eigenvalues are non-negative and real. Hence the Rayleigh quotient of an eigenvector ϕ of L associated to an eigenvalue λ is R π,P pϕq " λ, which is the definition of frequency considered in [START_REF] Shuman | The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains[END_REF].

On the content of subspaces associated to the random walk

The eigenvalues of an ergodic transition matrix P associated to a random walk on a strongly connected digraph G are either real or come as complex-conjugate pairs pϑ, θq P C 2 , ϑ " α `iβ. In the following, we examine how these individual eigenspaces can be understood with respect to each other.

Let us first consider pairs of eigenspaces related to complex-conjugate eigenvalues α ˘iβ. Assuming that the eigenspaces have dimension equal to 1, the eigenvectors can be chosen as complex-conjugate ξ, ξ P C N . These eigenvectors share the same frequency 1 ´α and can be seen as a generalization of the Fourier modes at frequencies ω and ´ω in one-dimensional classical Fourier analysis. One can also define the corresponding real-valued Fourier modes

ξ cos " ξ `ξ 2 , ξ sin " ξ ´ξ 2i ,
which generalize the cosine and sine functions to graph signals. As with the discrete Fourier transform, there are also frequencies for which the ω and ´ω frequency subspaces are the same (the zero frequency and possibly the 1{2 frequency in discrete signal processing). The case of real eigenvalues in graph signal processing can be seen as a generalization of these.

When the eigenspaces have dimension greater than 1, the situation is similar to the multidimensional classical Fourier analysis setting where the classical frequency is a vector ω " rω 1 , . . . , ω n s J . In graph signal processing, the ω frequency would then be a generalization of the norm of the classical frequency vector }ω}.

Continuing the analogy with classical Fourier analysis, in more than one dimensions Fourier modes with a given frequency ω are all oscillations at the same frequency but with different orientations. We argue that the case in graph signal processing of eigenvalues with the same real part but different non-conjugate imaginary parts may be a generalization of different orientations. This is illustrated in Section 4.5.2 on the directed toroidal graph. For arbitrary graphs, low frequency Fourier modes are typically oscillatory but some higher frequency modes may be highly localized which does not fit the common interpretation of an oscillation [START_REF] Nakatsukasa | Mysteries around the graph laplacian eigenvalue 4[END_REF][START_REF] Cloninger | Natural graph wavelet packet dictionaries[END_REF]. The interpretation that the imaginary parts are related to different orientations mostly holds for oscillatory modes. On digraphs, we can expect the same non trivial localization properties as in [START_REF] Nakatsukasa | Mysteries around the graph laplacian eigenvalue 4[END_REF][START_REF] Cloninger | Natural graph wavelet packet dictionaries[END_REF] to hold, hence the interpretation for orientations for localized modes should be considered differently.

Random walk graph filters

Having now tools for the Fourier analysis of digraphs, we define, thanks to that, random walk graph filter as being a LSI filter with respect to the reference operator P:

H " T ÿ t"0 h t P t .
Assuming that P is diagonalizable, we can represent a graph filter H as a linear combination of the spectral projectors tE ϑ k u m k"1 associated with P, as defined in (4):

H " m ÿ k"1 γ k E ϑ k , γ k P C, k " 1, . . . , m. ( 18
)
The filter is characterized in the Fourier domain by its spectral response coefficients, or Fourier coefficients, γ k , k " 1, . . . , m. As in classical signal processing, the application of a filter to a signal in the Fourier domain is obtained by multiplying the filter's Fourier coefficients with the signal's corresponding Fourier coefficients. A particularity here is that a filter may only have one Fourier coefficient for a Fourier subspace of dimension larger than 1 while the signal is generally described by as many coefficients as the dimension of the subspace.

From a filter design point of view, we may also be interested in filters defined by a frequency response. This is a special case of [START_REF] Ricaud | Fourier could be a data scientist: From graph fourier transform to signal processing on graphs[END_REF] where the spectral coefficients only depend on the frequency: γ k " hpω k q " hp1 ´Repϑ k qq, where h : R `Þ Ñ C is the frequency response of the filter.

Let us define the projector associated with the frequency ω as

S ω " ÿ ϑ:1´Repϑq"ω E ϑ . (19)
It is a projector that projects on the sum of the eigenspaces corresponding to the frequency ω. We call S ω the mono-frequency projector associated with the frequency ω. A filter based on a frequency response h can be expressed as

H ω " ÿ ωPω hpωqS ω , ω P r0, 2s. (20)
A major advantage of such a formulation is that it allows one to shift, contract or dilate a filter along the frequency axis by simply applying these operations to its frequency response. Unfortunately, this is not possible for all graph filters because a filter may have different Fourier coefficients γ k for different eigenspaces even when they share the same frequency ω.

Fourier analysis on finite groups: a graph signal perspective

In order to show the consistency of our Fourier analysis with respect to traditional signal processing, we depict the directed counterparts of two well-known objects: the cycle graph associated with the cyclic group pZ{nZq and the toroidal graph that is the direct product of cyclic groups, i.e. T " pZ{n 1 Zq ' ¨¨¨' pZ{n r Zq [START_REF] Terras | Fourier analysis on finite groups and applications[END_REF].

Fourier analysis on the directed cycle graph

A directed cycle graph C N is a graph with N vertices containing a single cycle through all the vertices and were all the edges are directed in the same direction. The directed cycle graph C N is represented by its adjacency matrix C N which is circulant As the out-degree matrix of C N is D " I N , the transition matrix for the random walk random walk X on C N is P " C N . P is diagonalizable, i.e. P " ΞΘΞ ´1 with Θ " diagtϑ 1 , ¨¨¨, ϑ N u P C N ˆN the eigenvalue matrix where ϑ k " e 2πipk´1q{N , k " 1, . . . , N. and Ξ " rξ 1 , ¨¨¨, ξ N s P C N ˆN is the discrete Fourier basis where

C N " » - - - - - - 0 0 . . .
ξ k " 1 ? N r1, ϑ 1 k , . . . , ϑ N ´1 k s J k " 1, ¨¨¨, N.
The random walk X is irreducible and periodic. As our frequency analysis only deals with ergodic random walks, we need X to be aperiodic. In order to overcome the periodicity issue, we consider a γ-lazy random walk X characterized by its transition matrix Pγ P P defined at Section 3.2.1. For such a lazy random walk, the stationary distribution π exists and is constant. As a consequence, the Rayleigh quotient associated with a given eigenvector ξ k is given by R pγq π,P pξ k q " p1 ´γq

" 1 ´cos ˆ2πpk ´1q N ˙ȷ " ωpγq k , k " 1, ¨¨¨, N.
Since the γ-lazy random walks X have the same eigenspaces as the non-lazy random walk, taking the limit γ Ñ 0 of these Rayleigh quotients allows us to define a frequency for the eigenspaces of the non-lazy random walk too.

ω k " " 1 ´cos ˆ2πpk ´1q N ˙ȷ, k " 1, ¨¨¨, N.
Thus, the ordering of the frequencies tω k u N k"1 associated with the eigenvectors

tξ k u N k"1
by the Rayleigh quotient coincides exactly with the classical signal processing approach. However, the actual values of the frequencies are warped by the function x Þ Ñ 1 ´cospxq compared to their classical counterparts. As the warping could in theory be reversed by modifying our definition of frequency, the notion of frequency that we propose is consistent with the classical notion of frequency.

Finally, the trick we used here with taking the limit γ Ñ 0 of γ-lazy random walks could be used for any graph where the random walk is irreducible but periodic. This effectively allows us to define a meaningful notion of frequency for any graph that is strongly connected.

Fourier analysis on the directed toroidal graph

A directed toroidal graph T m,n is the Cartesian product of the directed cycle graphs C m and C n , namely T m,n " C m ˝Cn . We introduce necessary definitions for the understanding of the section.

Definition 4.3 ([88]

). Let G " pU, Eq and H " pV, Fq be two graphs with respective vertex sets U " pu 1 , . . . , u n q and V " pv 1 , . . . , v m q. The Cartesian product of G and H is the graph G ˝H with vertex set U ˆV in which two vertices x " pu i , v j q and y " pu p , v q q are adjacent if and only if either u i " u p and pv j , v q q P F or v j " v q and pu i , u p q P E.

Definition 4.4 ([88]

). Let G " pU, Eq and H " pV, Fq be two graphs. The adjacency matrix of the Cartesian product G ˝H, denoted A G˝H is

A G˝H " A G b I |V| `I|U| b A H .
where b is the Kronecker product symbol. Lemma 4.1. Suppose λ 1 , ¨¨¨, λ n are eigenvalues of A G and µ 1 , ¨¨¨, µ m are eigenvalues of A H . Then the eigenvalues of A G˝H are all λ i `µj for 1 ď i ď n and 1 ď j ď m. Moreover, if u and v are eigenvectors for A G and A H with eigenvalues λ and µ respectively, then the vector w " u b v is an eigenvector of A G˝H with eigenvalue λ `µ.

Directed toroidal graph

Let T m,n " C m ˝Cn be the directed toroidal graph characterized by its adjacency matrix A Tm,n " C m b I n `Im b C n . The eigenvalues of C m are λ k " e 2iπpk´1q{m for k P v1, mw and the eigenvectors are denoted by rϕ 1 , . . . , ϕ m s. Identically, the eigenvalues of C n are µ ℓ " e 2iπpℓ´1q{n for ℓ P v1, nw and the eigenvectors are denoted by rψ 1 , . . . , ψ n s. The directed toroidal graph T m,n is a directed 2-regular graph. Consequently, the out-degree matrix of T m,n is D Tm,n " 2I mˆn . We thus define a random walk X on T m,n with a diagonalizable transition matrix P " D ´1 Tm,n A Tm,n . We specify its spectral properties in the following lemma. Lemma 4.2. Let T m,n be a directed toroidal graph and X be a random walk associated with the diagonalizable transition matrix P. The directed toroidal graph T m,n is a directed 2-regular graph. Consequently, the spectrum of P is σpPq "

" ζ i,j " λ i `µj 2
, i P v1, mw, j P v1, nw * , with associated eigenvectors

" ϕ i b ψ j , i P v1, mw, j P v1, nw * .
For the cycle graph, we observed that our notion of frequency was consistent with the classical notion of frequency used in signal processing up to some monotonous transform. This is only approximately true for the toroidal graph. Indeed, it is usually possible to find pairs of eigenvectors of the toroidal graph for which the frequency order differs compared to the classical case, although it does not seem to occur when the frequency difference is large enough. An example where classical frequencies and graph frequencies behave slightly differently is presented below.

As in Section 4.5.1, the associated random walk X is irreducible and periodic. In order to overcome this periodicity issue, we consider a γ-lazy random walk X with transition matrix Pγ P P.

Figure 1 illustrates the eigenvalue distribution of the transition matrix of a directed toroidal graph. As we can appreciate, the distribution of eigenvalues is remarkable. We denote by Λ g , the set of eigenvalues located on the green dotted line and Λ o the set of eigenvalues located on the orange dotted line. Sets Λ g , and Λ o , include respectively eigenvalues having the same real part and different imaginary parts i.e. RepΛ g q " ´0. [START_REF] Hammond | Wavelets on graphs via spectral graph theory[END_REF] and RepΛ o q " 0.25. In order to illustrate the insights discussed in Section 4.3.1, we exhibit in figure 2 the two-dimensional representation of some eigenvectors from the transition matrix of a directed toroidal graph T 36,60 corresponding to a graph frequency ω " p7´?5q{8. The eigenvalues associated with these eigenvectors have identical real parts but different non-conjugate imaginary parts. Analytically, the Rayleigh quotients of the selected eigenvectors are equal, as explained in Section 4.3.1. As we can see in figure 2, the selected eigenvectors are all (discretized) complex exponential waveforms with what look like identical or very similar frequencies. The main difference between them appears to be their orientation, which seems related to the different imaginary parts of the eigenvalues. On a side note, we also observe that the graph frequencies of these 4 eigenvectors are identical but this is not exactly true for their classical frequencies, which are equal to 0.260, 0.269, 0.260 and 0.269 from top to bottom respectively. . Eigenvalues on the dotted green line: Repλ g q " ´0.25, @λ g P Λ g . Eigenvalues on the dotted orange line: Repλ o q " 0.25, @λ o P Λ o .

Multiresolution analyses on digraphs

In the previous Sections, we proposed a Fourier-like basis on digraphs, as the set of random walk eigenvectors and determined a frequency analysis by studying their variation. This allows the definition of a low-pass or band-pass filter, hence we are now able to construct multiresolution analyses of functions over digraphs. In the first instance, we propose wavelet frames made of analysis and synthesis graph filter banks. This multi-scale construction is closely related to the diffusion polynomial frames [START_REF] Maggioni | Diffusion polynomial frames on metric measure spaces[END_REF] and spectral graph wavelets [START_REF] Hammond | Wavelets on graphs via spectral graph theory[END_REF]. In a second instance, we propose a critically sampled wavelet construction on digraphs generalizing the diffusion wavelets framework [START_REF] Coifman | Diffusion wavelets[END_REF]. For the sake of simplicity, we will present our theoretical constructions in the ℓ 2 pVq space. The linear transformation φ defined in eq. ( 13) and its inverse may be used to express the constructed wavelets and scaling functions in ℓ 2 pV, πq or to map measured signals from ℓ 2 pV, πq into ℓ 2 pVq.

Spectral wavelet transform on digraphs

In this Section, we propose a redundant wavelet transform on digraphs. We follow the construction of spectral wavelets on undirected graphs [START_REF] Hammond | Wavelets on graphs via spectral graph theory[END_REF] and polynomial diffusion frames [START_REF] Maggioni | Diffusion polynomial frames on metric measure spaces[END_REF]. The novelty is the construction of filters designed in the frequency domain via a frequency response function as a linear combination of mono-frequency projectors defined in eq. [START_REF] Lovász | Random walks on graphs[END_REF]. We introduce the necessary elements for the construction of redundant wavelets on digraphs.

Theoretical framework

Let G " pV, Eq be a strongly connected digraph with cardinality |V| " N . G is characterized by its adjacency matrix W P R N ˆN `. On G, we define a random walk X imaginary part real part characterized by its transition matrix P " D ´1W P R N ˆN `. X is assumed ergodic with stationary distribution π. Let T " Π 1{2 PΠ ´1{2 P ℓ 2 pVq be the operator similar to P in ℓ 2 pVq. We assume T diagonalizable. In the graph signal processing framework, the concept of low-pass, high-pass, and band-pass frequency response is intrinsically related to the frequency interpretation of the eigenvectors (or equivalently the eigenspaces) of the reference operator. In Section 4.3, we have determined a frequency interpretation of the eigenvectors of the random walk operator on digraphs. Let ω 1 ď ¨¨¨ď ω N be the random walk's graph frequencies ordered from the lowest frequency ω 1 to the highest frequency ω N . Let H ω be the graph filter defined in [START_REF] Coulhon | Random walks on graphs with regular volume growth[END_REF]. Let us define the cut-off random walk's graph frequency ω cut . The ideal low-pass and high-pass frequency responses, respectively h and g are defined by

hpω m q " 1 ´gpω m q " " 1 if ω m ď ω cut 0 if ω m ą ω cut .
In addition, given two random walk's graph frequencies ω ℓ and ω k , ω k ă ω ℓ , an ideal band-pass frequency response u is defined by

upω k q " " 1 if ω m ď ω k ď ω ℓ 0 if ω k ą ω ℓ and ω k ă ω m .
Hence, an ideal low-pass (resp. high-pass) operator is associated with a low-pass (resp. high-pass) frequency response. As a practical example, it is convenient to use low pass and high pass filters having cut-off frequency ω cut equal to the median of the bandwidth i.e. such that exactly half of the frequencies ω m are lower frequencies that ω cut . More generally, the term low-pass and bass-pass filters can be assumed in a relaxed way, using the same intuitive meaning as in classical signal processing. We call filters low-pass if they do not significantly affect the frequency content of low-frequency signals but attenuate, i.e., reduce, the magnitude of high-frequency signals. Analogously, highpass filters pass high-frequency signals while attenuating low-frequency ones, and bandpass filters pass signals with frequency content within a specified frequency band while attenuating all others. Definition 5.1. A low pass operator at dilation t, H t on a strongly connected digraph G is defined by

H t " p ÿ ℓ"1 hptω ℓ qS ℓ , t P N.
with h : R `Ñ R `, a function giving a low-pass frequency response and tS ℓ u p ℓ"1 , the set of projectors (19) associated to the frequencies ω ℓ of the eigenvectors of T. The scaling function at dilation t and translation y is denoted by h t,y " H t δ y where δ y is the Kronecker delta function at vertex y P V. Definition 5.2. A band pass operator at dilation t, G t on a strongly connected digraph G is defined as

G t " p ÿ ℓ"1 gptω ℓ qS ℓ , t P N.
with g : R `Ñ R `, a function giving a band-pass frequency response and tS ℓ u p ℓ"1 , the set of projectors (19) associated to the frequencies ω ℓ of the eigenvectors of T. The wavelet function at dilation t and translation y is denoted by

g t,y " G t δ y ,
where δ y is the Kronecker delta function at vertex y P V.

Having low-pass and band-pass operators on digraphs, we are able to build analysis and synthesis filter banks on digraphs. Definition 5.3. We define a bank of synthesis filters K as a set comprising a low-pass filter at dilation t J , H t J and a series of band-pass filters at increasing dilations tt j u J j"1 , tG t j u J j"1 :

K " tH t J , G t 1 , . . . , G t J u.
Definition 5.4. We define a bank of analysis filters K as a set comprising a filter at dilation t J , HG t J and a series of filters at increasing dilations tt j u J j"1 , t Gt j u J j"1 :

K " t Ht J , Gt 1 , . . . , Gt J u,
where

HG t J " p ÿ ℓ"1 hptω ℓ qS ℓ P R N ˆN ,
and Gt j "

p ÿ ℓ"1
gptω ℓ qS ℓ P R N ˆN , @j " 1, . . . , J.

We also define ht,k and gt,k as row vectors as follows:

ht,k " δ J k Ht , gt,k " δ J k Gt . Proposition 5.1. Given a fixed set of increasing dilations tt j u J j"1 , the synthesis and analysis filters sets respectively K and K, the perfect reconstruction condition

H t J Ht J `J ÿ j"1 G t j Gt j " I (21)
is guaranteed if and only if

hpt J ω ℓ q hpt J ω ℓ q `J ÿ j"1
gpt j ω ℓ qgpt j ω ℓ q " 1, @ℓ " 1, . . . , p.

Proof.

H t J Ht J `J ÿ j"1 G t j Gt j " p ÿ ℓ,ℓ 1 "1 hpt J ω ℓ q hpt J ω ℓ 1 qS ℓ S ℓ 1 `J ÿ j"1 p ÿ ℓ,ℓ 1 "1 gpt j ω ℓ qgpt j ω ℓ 1 qS ℓ S ℓ 1 " p ÿ ℓ"1 hpt J ω ℓ q hpt J ω ℓ qS ℓ `J ÿ j"1 p ÿ ℓ"1 gpt j ω ℓ qgpt j ω ℓ qS ℓ " p ÿ ℓ"1 " hpt J ω ℓ q hpt J ω ℓ q `J ÿ j"1 gpt j ω ℓ qgpt j ω ℓ q ‰ S ℓ
Consequently, ( 21) is guaranteed if and only if hpt J ω ℓ q hpt J ω ℓ q `J ÿ j"1 gpt j ω ℓ qgpt j ω ℓ q " 1, @ℓ " 1, . . . , p.

We now define a frame as follows.

Definition 5.5. Let G " pV, Eq be a digraph and µ a measure on V. A countable family of elements tf k u n k"1 P ℓ 2 pVq is said to be a frame if for any graph signal f P ℓ 2 pVq we have

A}f } 2 ď n ÿ k"1 |xf , f k y| 2 ď B}f } 2 ,
for some constants 0 ă A ď B ă 8 which are called frame bounds.

We introduce the rectangular matrix K " pH t J , G t 1 , . . . , G t J q P R N ˆN pJ`1q where : K " `ht J ,1 , . . . , h t J ,N , g t 1 ,1 , . . . , g t J ,1 , . . . , g t J ,N ˘P R N ˆN pJ`1q .

We also introduce the following rectangular matrix K " p Ht J , Gt 1 , . . . , Gt J q J P R N pJ`1qˆN K "

¨h t J ,1 . . . ht J ,N gt 1 ,1 . . . gt J ,1 . . . gt J ,N ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' P R N pJ`1qˆN . (22)
Proposition 5.2. Let K and K be respectively the synthesis and analysis filter banks. We assume that the filter banks K and K verify the perfect reconstruction condition stated in proposition 5.1. Consequently, K is a frame for ℓ 2 pVq with lower frame bound 1{}K} 2 and upper frame bound } K} 2 .

Proof. The perfect reconstruction condition implies that }f }

2 " }K Kf } 2 ď }K} 2 } Kf } 2 .
Hence, we have the frame inequalities

1 }K} 2 }f } 2 ď } Kf } 2 ď } K} 2 }f } 2 , ( 23
)
where the second inequality is just the definition of the spectral norm. }K} 2 and } K} 2 are finite because of finite dimensions, hence the frame bounds verify 0 ă }K} 2 ď } K} 2 ă 8.

Definition 5.6 (Wavelet decomposition of a signal on a digraph). Any graph signal f P ℓ 2 pVq can be expressed as follows :

f " N ÿ k"1 xf , ht J ,k y h t J ,k `J ÿ j"1 N ÿ k"1
xf , gt j ,k y g t j ,k .

Diffusion wavelet transform on digraphs

We first recall relevant and useful notations for a clear understanding of the construction of diffusion wavelets [START_REF] Coifman | Diffusion wavelets[END_REF][START_REF] Maggioni | Biorthogonal diffusion wavelets for multiscale representations on manifolds and graphs[END_REF]. Two notations are introduced: one for the representation of linear transformations as matrices and a second for the representation of sets of vectors as matrices where columns correspond to a vector in the set. These notations were used in [START_REF] Coifman | Diffusion wavelets[END_REF][START_REF] Maggioni | Biorthogonal diffusion wavelets for multiscale representations on manifolds and graphs[END_REF] except that we adopt a column vector convention whereas a row vector convention was used in [START_REF] Coifman | Diffusion wavelets[END_REF][START_REF] Maggioni | Biorthogonal diffusion wavelets for multiscale representations on manifolds and graphs[END_REF].

Let V 0 " ℓ 2 pVq be the space of functions defined over the vertices of a digraph G " pV, Eq. If L is a linear transformation of V 0 into V 0 , rLs B 2 B 1 indicates the matrix representing the linear transformation L with respect to the basis B 1 in the domain and B 2 in the range. A set of vectors X represented in a given basis B will be written in the matrix form rXs B where the columns of rXs B are the coordinates of the vectors X in the coordinate system defined by B. Generally, rB 1 s B 2 " rIs B 2 B 1 represent the basis vectors B 1 in terms of the basis B 2 . We will note that, for a given basis B if the input and output bases are the same, rIs B B " I the identity matrix. We will also abuse this notation if B 1 spans a subspace of the space spanned by B 2 . If B 2 spans a subspace of the space spanned by B 1 instead, then we use the convention rIs B 2 B 1 " prIs B 1 B 2 q `, where X `denotes the Moore-Penrose pseudoinverse of X. If B 1 and B 2 are two bases, the representations of X in B 1 and B 2 are related as follows:

rXs B 2 " rIs B 2 B 1 rXs B 1 . If B 1 , B 2 , B 3 , B 4
are arbitrary bases, the matrix representation of L with respect to the basis B 1 in the domain and B 2 in the range in the different bases is expressed as

rLs B 2 B 1 " rIs B 2 B 4 rLs B 4 B 3 rIs B 3 B 1 .
Furthermore, if B 1 and B 2 are linearly independent sets of vectors that do not span the whole space V 0 then we will still use the notation rLs B 2 B 1 , but in that case it will represent the restriction of the linear transformation L to the domain and range subspaces spanned by B 1 and B 2 .

Diffusion wavelets

The construction of the diffusion wavelets [START_REF] Coifman | Diffusion wavelets[END_REF][START_REF] Maggioni | Biorthogonal diffusion wavelets for multiscale representations on manifolds and graphs[END_REF] enables a multi-resolution analysis of functions on graphs generalizing the concept of classical multiresolution analysis [START_REF] Mallat | A theory for multiresolution signal decomposition: the wavelet representation[END_REF]. The starting point of this construction is a diffusion operator T. Similarly to the classical multiresolution analysis, diffusion wavelets is characterized by a family of nested subspaces

V 0 Ě V 1 Ě ¨¨¨Ě V j Ě ¨¨¨,
where each subspace V j is spanned by a basis of scaling functions Φ j . The complement of V j`1 into V j is called W j and is spanned by a set of diffusion wavelets Ψ j .

Construction

The construction of diffusion wavelets proceeds with a diffusion operator T defined on a directed strongly connected graph G " pV, Eq. In the original diffusion wavelets framework by Coifman and Maggioni [START_REF] Coifman | Diffusion wavelets[END_REF], the graph is assumed undirected and they suggest using the reversible random walk operator P.

More generally, any operator T can be used as long as it acts like a low-pass filter. Given a digraph G, an appropriate choice seems to be the operator T " Π 1{2 PΠ ´1{2 where P is the random walk operator. However, the random-walk operator is not guaranteed to have a low-pass frequency response as the highest frequency eigenvectors may have eigenvalues with modulus arbitrarily close to one (as e.g. for the cycle graph). In such a case, a low-pass filter can be defined using a lazy random walk instead: T " Π 1{2 pγI p1 ´γqPqΠ ´1{2 . Assuming an ordering of vertices in V, T is originally represented on the canonical basis Φ 0 " tδ k u kPV of V 0 where δ k is the Kronecker delta function corresponding to the vertex k P V. Using the notations from Section 5.2, rTs Φ 0 Φ 0 is the matrix representation of the linear operator T with respect to the basis Φ 0 in the domain and Φ 0 in the range.

The columns of the matrix rTs Φ 0 Φ 0 can be described as a set of functions r Φ 1 " t r ϕ 1,k u kPV represented in the basis Φ 0 . Each element r ϕ 1,k , k P V of r Φ 1 corresponds to the action of T on the Kronecker δ-function at vertex k and is represented as r r ϕ 1,k s Φ 0 " rTs Φ 0 Φ 0 δ k . The low-pass property of T leads to the fact that each element r ϕ 1,k of r Φ 1 is a function localized around its vertex k whose support extends to its close neighbors. In terms of the diffusion wavelets construction, the elements r ϕ 1,k , k P V are therefore interpreted as scaling functions. Due to the fact that the support of each element r ϕ 1,k covers a small neighborhood around their respective vertices, these elements r ϕ 1,k can generally be well approximated by a linear combination of the other functions r ϕ 1,l with l ‰ k. The next stage of the construction is a column subset selection stage [START_REF] Civril | Column subset selection for approximating data matrices[END_REF]. The aim is to find a small number of columns of Φ1 , forming a set C such that the residual }r Φ1 s Φ 0 ´ΠC r Φ1 s Φ 0 } β is as minimal as possible where Π C is the projection matrix onto the space spanned by the columns of C and β denotes the spectral norm or the Frobenius norm.

More precisely, this step involves the selection of a subset t r

ϕ 1,k , k P I 1 u from r Φ 1 , |I 1 | ď |V|.
We select a subset t r ϕ 1,k , k P I 1 u such that all r ϕ 1,k are well-enough approximated by linear combinations of the functions in the subset. In classical signal processing terms, this step is analogous to a subsampling of a set of scaling functions in the classical discrete wavelet transform. Coifman and Maggioni used a greedy approach to build the set I 1 iteratively by using a modified Gram-Schmidt orthogonalization procedure [START_REF] Coifman | Diffusion wavelets[END_REF].

The subset t r ϕ 1,k , k P I 1 u spans a subspace V 1 which corresponds to the first approximation subspace of the multi-resolution analysis. We will thus denote t r ϕ 1,k , k P I 1 u as Φ 1 , where Φ 1 is by definition a basis of V 1 that is generally not orthogonal. The vectors in Φ 1 are the scaling functions at scale 1.

To define the next scales of scaling functions, we consider the compression step of the diffusion operator T on the subspace V 1 , that is its restriction on V 1 . The latter can be represented in Φ 1 as

rTs Φ 1 Φ 1 " rIs Φ 1 Φ 0 rTs Φ 0 Φ 0 rIs Φ 0 Φ 1 . where rIs Φ 1
Φ 0 represents the restriction of the signal space V 0 to V 1 (with respective bases Φ 0 and Φ 1 ) and rIs Φ 0 Φ 1 represents the embedding of V 1 in V 0 . The next approximation space V 2 Ă V 1 and its associated basis Φ 2 can be obtained in the same way as the definition of V 1 and its basis Φ 1 except that we now consider the operator

T 2 restricted to V 1 instead of T in V 0 .
The columns of rT 2 s Φ 1 Φ 1 « prTs Φ 1 Φ 1 q 2 can be interpreted as scaling functions at scale 2, r Φ 2 " t r ϕ 2,k u, represented in the basis Φ 1 . From these functions we extract a subset Φ 2 " t r ϕ 2,k , k P I 2 u such that any function in r Φ 2 is well-approximated by a linear combination of functions of Φ 2 .

After j iterations of this procedure we have defined j approximation subspaces V j Ă V j´1 Ă ¨¨¨Ă V 1 with corresponding bases Φ j , Φ j´1 , . . . , Φ 1 . At each step the basis Φ j is defined by its representation in the basis Φ j´1 based on the restriction of the operator T 2 j to V j´1 . In order to represent these functions in the original basis Φ 0 of V 0 we can use

rΦ j s Φ 0 " rIs Φ 0 Φ j " rIs Φ 0 Φ 0 rIs Φ 0 Φ 1 ¨¨¨rIs Φ j´2 Φ j´1 rIs Φ j´1
Φ j . Since every function in Φ 0 is defined on V 0 , so is every function on Φ j . Hence any function in the approximation space V j can be extended naturally to the whole space V 0 .

Regarding the construction of the wavelets, we propose to construct the wavelet bases Ψ j for the subspaces W j by selecting a subset of the columns of the band pass operator rI V j ´Φj`1 Φ : j`1 s Φ j Φ j which is the orthogonal projector on the complement of V j`1 into V j . The wavelets capture the detail lost from going from V j to V j`1 . As our framework falls into a bi-orthogonal scope, we need to build the dual wavelet bases. For a scale j, we have a wavelet base Ψ j and we need to construct the associated dual base p Ψ j obtained such that we have the following relation

Ψ 1: j Ψ j " I. where Ψ 1:
j is the pseudo-inverse of Ψ j . We mention that in [START_REF] Maggioni | Biorthogonal diffusion wavelets for multiscale representations on manifolds and graphs[END_REF], bi-orthogonal wavelet transform is proposed but only the scaling functions are actually defined. Our construction is identical but we also propose a definition for wavelets.

A generalization using more arbitrary scaling operators

We propose a generalization of the diffusion wavelet framework that enables it to be combined for instance with the wavelet construction presented in Section 5.2. The idea is merely to replace the powers of T as multiresolution scaling operators by more arbitrary low-pass filters. More precisely, where the operator T 2 j is used to define the approximation space V j`1 in section Section 5.2.2, we propose to use instead a low-pass graph filter H j . If the graph filters H j correspond to the scaling operators defined in Section 5.1 with appropriately increasing scales, this approach provides a way to reduce the redundancy of the sets of scaling functions. Similarly, the same approach could be used to reduce the redundancy of the wavelet functions defined in Section 5.2.

Applications

This section is devoted to illustrating our theoretical framework through examples of semi-supervised learning and signal modeling on digraphs. It turns out that in the following applications, we are able to work either in the Hilbert space ℓ 2 pVq or ℓ 2 pV, πq where π is the unique stationary distribution of an ergodic random walk. Furthermore, by the results obtained in Section 3.4, we define the isometries φ and φ ´1 going from ℓ 2 pVq to ℓ 2 pV, πq.

φ

: f Þ Ñ Π ´1{2 f , φ ´1 : f Þ Ñ Π 1{2 f , @f P ℓ 2 pV, πq. (24)
For the sake of simplicity, we will always map signals living in ℓ 2 pV, πq to ℓ 2 pVq, so that we can work in ℓ 2 pVq where the inner product coincides with the one between vectors in C N .

In practice, given a graph signal f measured in the real world, we must first decide in which space the signal lives. This is a somewhat arbitrary choice between the signal living in a space associated with the counting measure or in a space associated with the stationary measure π of an ergodic random walk. If we decide that the signal lives in ℓ 2 pV, πq, we then convert the signal into ℓ 2 pVq using the application φ ´1, we perform all the necessary operations in ℓ 2 pVq and then go back to the ℓ 2 pV, πq space using φ. In ℓ 2 pVq, the operator T defined in eq. ( 14) corresponds to the random walk operator and L is the operator such that xf , Lf y is the Dirichlet energy.

Semi-supervised learning on digraphs via ℓ 2 -regularization

We discuss the semi-supervised learning problem on digraphs with a regularization term of type ℓ 2 . The following problem aims to show the efficiency and relevance of the Dirichlet energy ( 16) as a regularization term for signal processing on digraphs.

Definition of the problem

Let G " pV, Eq be a strongly connected digraph where V " pv 1 , . . . , v N q is the vertex set, E is the edge set and with cardinality |V| " N . G is represented by its adjacency matrix W P R N ˆN `. Let X be the associated random walk on G. The random walk X is represented by its random walk operator P P R N ˆN `with unique stationary distribution π.

We also introduce G sym the symmetrized version of G. The undirected graph G sym is represented by its adjacency matrix W sym " pW `WJ q{2 P R N ˆN `. Let X sym be the associated random walk on G sym . The random walk X sym is represented by its random walk operator P sym P R N ˆN `with unique stationary distribution π sym . Let y : V Ñ t´1, 1u be a "label" function defined on the vertex set V. From the function y, we obtain the graph signal y defined as follows y " rypv 1 q, . . . , ypv N qs J .

Given the values ypv k q t1ďkďN u on only a subset of labeled vertices U Ă V, the aim is to estimate the labels of the remaining unlabeled vertices. This is the semi-supervised learning problem on graphs [START_REF] Belkin | Regularization and semi-supervised learning on large graphs[END_REF][START_REF] Zhu | Semi-supervised learning with graphs[END_REF][START_REF] Zhou | Learning from labeled and unlabeled data on a directed graph[END_REF]. We introduce the formulation of the problem in the next section.

Considered approaches

The standard approach with ℓ 2 regularization is based on the following optimization problem:

argmin f ␣ c}M l pf ´yq} 2 `ϱ1 }M u f } 2 `ϱ2 Spf q ( , c, ϱ 1 , ϱ 2 ą 0. ( 25
)
The term }M l pf ´yq} 2 is the data fidelity term with M l " diagtm 1 , . . . , m N u a diagonal matrix where m i " 1 v i PU , for all i " 1 . . . , N (hence 0 on vertices with unknown labels).

In the second term, }M u f } 2 , M u selects the unlabeled vertices and is defined as M u " I N ´Ml . This second term is a Tikhonov-like regularization term. The third term, Spf q, is the variational regularization term derived from the graph Spf q " xf , Xf y, X P R N ˆN .

The semi-supervised learning on graphs problem is formulated in the same manner for the digraph G or its symmetrized version G sym , except that the variational regularization term Spf q differs. In all these situations, for learning two classes, the solution is obtained by constraining the previous continuous problem of having labels taking values only in t`1, ´1u; a relaxed version of the problem is to solve the continuous optimization problem and take the sign of its solution. We now compare three methods for the directed case and the symmetrized case. Note that we restrict the discussion to c " ρ 1 (or ρ 1 " 0 as it is usually done in the literature [START_REF] Belkin | Regularization and semi-supervised learning on large graphs[END_REF][START_REF] Zhou | Learning from labeled and unlabeled data on a directed graph[END_REF], and because it leaves only one trade-off, between the supervised loss (possibly with a classical Tikhonov term) and the graph regularization term. It simplifies also the expression of the closed form solutions.

(a) Semi-supervised learning on G, the directed graph:. Method 1: The first method is a version of eq. ( 25), that was previously investigated by [START_REF] Zhou | Learning from labeled and unlabeled data on a directed graph[END_REF].

Assuming the graph signal lives in ℓ 2 pVq, the optimization problem is

argmin f ␣ c}M l pf ´yq} 2 `c}M u f } 2 `ϱ2 xf , Lf y ( , c, ϱ 2 ą 0. ( 26
)
with L the directed normalized Laplacian [START_REF] Luxburg | A tutorial on spectral clustering[END_REF]. The term xf , Lf y is the Dirichlet energy in ℓ 2 pVq.

Method 2: For this method, we also use the formulation [START_REF] Rustamov | Wavelets on graphs via deep learning[END_REF], except that we now assume that the graphs signal lives in ℓ 2 pV, πq. Hence, we define ỹ P ℓ 2 pVq as follows ỹ " φ ´1pyq " Π 1{2 y, y P ℓ 2 pV, πq.

The optimization problem is therefore

argmin f ␣ c}M l p f ´ỹq} 2 `c}M u f } 2 `ϱ2 x f , L f y ( , c, ϱ 2 ą 0. ( 27
)
Equation ( 27) is the same as eq. ( 26) except for y which is now ỹ. The optimization problem ( 27) can be rewritten as

argmin f ␣ c}M l pf ´yq} 2 π `c}M u f } 2 π `ϱ2 xf , L RW f y π ( , c, ϱ 2 ą 0. ( 28
)
with L RW the random walk Laplacian on directed graphs [START_REF] Rosenberg | The Laplacian on a Riemannian manifold: an introduction to analysis on manifolds[END_REF]. The term xf , L RW f y π is then the Dirichlet energy [START_REF] Belkin | Convergence of laplacian eigenmaps[END_REF].

Method 3: In the framework of "Discrete Signal Processing on Graphs" [START_REF] Sandryhaila | Discrete signal processing on graphs[END_REF][START_REF] Sandryhaila | Discrete signal processing on graphs: Frequency analysis[END_REF], the semisupervised learning problem for digraphs is formulated as follows

argmin f ␣ c}M l pf ´yq} 2 `ϱ2 }f ´Wnorm f } 2 2 ( , c, ϱ 2 ą 0, (29)
with W norm , a normalized version of the the adjacency matrix W so that its spectral norm is equal to one.

Semi-supervised learning on G sym , the symmetrized graph:. We consider the same methods applied now on G sym .

Method 1: Given the graph signal y P ℓ 2 pVq, the optimization problem is

argmin f ␣ c}M l pf ´yq} 2 `c}M u f } 2 `ϱ2 xf , L sym f y ( , c, ϱ 2 ą 0. ( 30
)
with L sym the normalized Laplacian defined as

L sym " I ´Π1{2 sym P sym Π ´1{2 sym , Π sym " diagtπ sym pv 1 q, . . . , π sym pv N qu.
Method 2: The optimization problem is

argmin f ␣ c}M l p f ´ỹq} 2 `c}M u f } 2 `ϱ2 x f , L sym f y ( , c, ϱ 2 ą 0. ( 31
)
This can be rewritten as

argmin f ␣ c}M l pf ´yq} 2 π `c}M u f } 2 π `ϱ2 xf , L RW,sym f y π ( , c, ϱ 2 ą 0. ( 32
)
with L RW,sym the random walk Laplacian defined as L RW,sym " I ´Psym .

Solutions in closed form

The problem ( 25) is quadratic and convex and therefore admits a closed form solution.

Semi-supervised learning on G.

Method 1: The solution of the optimization problem ( 26) is

f ˚" pI `γLq ´1M l y, γ " ϱ 2 {c, f ˚P ℓ 2 pVq. ( 33
)
Method 2: The solution of the optimization problem ( 27) is

f ˚" pI `γL RW q ´1M l y, γ " ϱ 2 {c, f ˚P ℓ 2 pV, πq. ( 34
)
Method 3: The solution of the optimization problem ( 29) is

f ˚" `Ml `γR ´1 M ˘y, γ " ϱ 2 {c, f ˚P ℓ 2 pVq.
with R M " pI ´Wnorm q J pI ´Wnorm q.

Semi-supervised learning on G sym .

Method 1: The solution of the optimization problem ( 30) is

f ˚" pI `γL sym q ´1M l y, γ " ϱ 2 {c, f ˚P ℓ 2 pVq. ( 35
)
Method 2: The solution of the optimization problem ( 31) is

f ˚" pI `γL RW,sym q ´1M l y, γ " ϱ 2 {c, f ˚P ℓ 2 pV, π sym q. ( 36
)
Remember that in each case, we take the approximated solution for 2 classes which keeps only the sign of the solution as the label of the inferred class: y ˚" sgnpf ˚q. Moreover, in order to optimize the parameter γ, we use the same approach as many authors [START_REF] Belkin | Regularization and semi-supervised learning on large graphs[END_REF][START_REF] Zhou | Learning from labeled and unlabeled data on a directed graph[END_REF]: a grid search γ P r0, 100s on a controlled example (with ground-truth) to have a numerical validation that the method can work well using some optimal γ. We summarize the different methods in Table 1. Application of semi-supervised learning on directed stochastic block model In this section, we apply a semi-supervised learning framework on a generalization of the stochastic block model [START_REF] Holland | Stochastic blockmodels: First steps[END_REF] (SBM) for digraphs. In contrast to the classical SBM, the additional parameters of the model we propose are used to assign different probabilities to the directions of the edges across different clusters. Let G " pV, Eq be a directed graph with cardinality |V| " N . Our directed stochastic block model GpS, p in , p out , p, qq is defined by S " tn 1 , . . . , n k u the set of the number of points associated with a given cluster with k ě 2, p in the probability that there is an edge between two vertices within the same clusters, p out the probability that there is an edge between two vertices belonging to two different clusters, p the probability that there is an edge orientation between two vertices belonging to the same cluster and q the probability that there is an edge orientation from cluster ℓ to j, ℓ, j ď k given that there is an edge between two vertices belonging to clusters ℓ and j respectively. Our purpose is to highlight the effectiveness of the Dirichlet energy based on the random walk operator compared to the Dirichlet energy based on the symmetrization of the adjacency matrix for semi-supervised learning problems on digraphs.

Method G G sym Method 1 L L sym Method 2 L RW L RW,sym Method 3 [55] R M -
Numerical simulations. We generate one instance from G 1 pS, p in , p out , p, qq where S " tn 1 " 150, n 2 " 150u, p in " 0.1, p out " 0.01, p " 0.5, q " 0.5 and one instance from G 2 pS, p in , p out , p, qq where S " tn 1 " 150, n 2 " 150u, p in " 0.1, p out " 0.01, p " 0.5, q " 0.9. The parameters of G 1 and G 2 are identical except for the q parameter. The edges of G 2 are much more oriented from cluster 1 to cluster 2 pq " 0.9q while the edges of G 1 are uniformly oriented pq " 0.5q between cluster 1 and cluster 2. Figure 3 shows the performance of the semi-supervised learning on the instances generated by G 1 and G 2 . Regarding the instance generated by G 1 , we remark that for low proportions of known labels, the semi-supervised approach based on the normalized directed Laplacian L performs better than its unnormalized counterpart L (L dir in Figure 3) and the symmetrized counterparts based on L sym and the unnormalized one L sym (L in Figure 3) while for higher proportions of known labels, the semi-supervised approach based on symmetrized unnormalized L sym gives the best performance. Regarding the instance generated by G 2 , the semi-supervised approach based on L strongly outperforms the symmetrized counterparts for all proportions of known labels. We show that when the orientations of the edges are important (q " 0.9), the semi-supervised approach based on the normalized directed Laplacian works is much more effective than the semi-supervised approach based on the symmetrized counterpart. Consequently, the directionality plays a strong role and the Dirichlet energy based on the random walk operator is meaningful for strong edges orientation.

Application of semi-supervised learning on political blogs

Let us consider the dataset of the political blogs of the 2004 US presidential campaign [START_REF] Adamic | The political blogosphere and the 2004 us election: divided they blog[END_REF]. The dataset consists of 1224 blogs, each associated with a political orientation, either republican or democrat. The hyperlinks between blogs provide the dataset with a graph structure G " pV, Eq, where each vertex v P V is associated to a blog and an edge between two vertices tv i , v j u indicates the presence of hyperlinks from blog i to blog j. The political orientations of the blogs are modeled by a graph signal f 0 " tf 0 pv 1 q, ¨¨¨, f 0 pv N qu where f 0 pv i q " `1 if blog i has a Democrat orientation and f 0 pv i q " ´1 if it is Republican. The directed graph G is not strongly connected, which is a problem with our framework that requires strongly connected graphs. As a result, we consider the largest strongly connected subgraph of G, denoted by G 1 " pV 1 , E 1 q which is made up of |V 1 | " N 1 " 793 vertices (hence roughly 65% of the vertex set V) and its associated graph signal f 1 0 . The respective performances of the approaches eqs. ( 26), ( 28) and ( 29) and their symmetrized counterparts eqs. ( 30) and ( 32) are compared in fig. 4. Numerical simulations. Figure 4 shows the performances of the semi-supervised approaches eqs. ( 26), ( 28) to [START_REF] Sweldens | The lifting scheme: A custom-design construction of biorthogonal wavelets[END_REF] and [START_REF] Candès | New tight frames of curvelets and optimal representations of objects with piecewise c2 singularities[END_REF], in terms of the proportion of correct matching between y ˚and f 0 (or f 1 0 ). For each method, several values of γ P r0, 10s were tested and we kept the one giving the best performance over 500 realizations. The lines displayed in fig. 4 are averaged estimation accuracy over 500 realizations and the error bars show the corresponding standard deviations.

Firstly, we remark that for all proportions of known labels, the semi supervised approaches based on L, L RW and its symmetrized counterparts L sym , L RW,sym perform better than the semi-supervised approaches based on R M . On the top left of fig. 4, we compare the semi-supervised approaches involving the operators L, L RW and R M . The semisupervised approach associated with L performs the best compared to the approaches involving L RW , R M for low proportions of known labels. Furthermore, L performs slightly better than the approach involving L RW . On the top right of fig. 4, we compare the semi-supervised approaches involving the operators L sym , L RW,sym and R M . Similarly the top left figure, the semi-supervised approach associated with L sym performs slightly better than the approach involving L RW . On the bottom left of fig. 4, we compare the semi-supervised approaches involving the operators L RW , L RW,sym and R M . The approach involving L RW and L RW,sym perform quite similarly for any proportion of known labels. The same can be observed in the bottom right of fig. 4 comparing the semi-supervised approaches involving L RW , L RW,sym and R M . This performance comparison leads us to the following conclusions:

• The semi supervised learning approaches involving Dirichlet energies based upon L and L RW (and their symmetrized counterparts) yield better performance compared to the approach proposed by Sandryhaila and Moura [START_REF] Sandryhaila | Discrete signal processing on graphs[END_REF] based upon R M . Consequently, semi-supervised learning approaches on digraphs with Dirichlet energies based on the random walk operator seem more appropriate.

• The Hilbert space where we assume the graph signal lives seems to have an influence on the performance. Here, we obtain better performance when assuming that the graph signal y belongs to ℓ 2 pVq rather than ℓ 2 pV, πq, or equivalently when using L in the regularization term rather than L RW .

Signal modeling on directed graphs via filtering

In this section, we consider a way to model the relationships between the values of a graph signal using a graph filter. The model is expressed as a graph filter that takes some values from the graph signal and reconstructs the other values. A possible application of such a model could be the lossy compression of the signal, where knowing only the graph, the coefficients of the graph filter and a few signal values enable the reconstruction of the whole signal. The major difference with the previous application is that we assume the whole graph signal known to learn the filter.

Problem formulation

Let G " pV, Eq be a digraph with cardinality |V| " N and µ : V Ñ R `a positive measure. The digraph G is characterized by its adjacency matrix W P R N ˆN `. Let f 0 P ℓ 2 pV, µq be a known graph function and f 0 " rf 0 pv 1 q, . . . , f 0 pv N qs J , the corresponding graph signal. Let y be a random sampling of f 0 according to some known sampling strategy, where for a given sampling realization y j " f 0 pv j q or y j " 0 depending on whether vertex j is selected in that sampling realization. The goal of the problem is to find a graph filter H that optimally reconstructs f 0 from the samples in y. More specifically, we consider graph filters expressed as finite polynomials of a reference operator R with order K.

H " K ÿ k"0 θ k R k , θ k P R, k " 0, ¨¨¨, K. (37)
The reconstruction of f 0 from y is then f0 " Hy "

K ÿ k"0 θ k R k y. (38)
In order to find an optimal filter p H " ř K k"0 θk R k , we propose to minimize the expected quadratic reconstruction error:

θ " argmin θ"tθ k u K k"0 PR K`1 E › › ›f 0 ´K ÿ k"0 θ k R k y › › › 2 µ , (39)
Finally, we only consider sampling strategies where samples in y are selected independently according to a given distribution δ, that is y j " f 0 pv j q with probability δ j and y j " 0 otherwise. The random variables y j , j " 1, . . . , N can be expressed as y j " ε j f 0 pv j q, where the ε j are independent and Bernoulli distributed with parameters δ j : ε j " Berpδ j q. The random variable ε defined as

p " 1 N N ÿ i"1 ε i , (40)
is the proportion of known values of f 0 collected in y.

Solution in closed form

The solution of the problem (39) is the following

p θ " Z ´1MQ J f 0 , (41)
where M " diagtµpv 1 q, ¨¨¨, µpv N qu is the diagonal matrix of the associated measure µ, Z the matrix where each entry Z kℓ P Z is expressed as

Z kℓ " Tr `rR k s J MR ℓ Etyy J u ˘, @tk, ℓu P v0, Kw 2 ,
and Q " rq 0 , ¨¨¨, q K s P R N ˆpK`1q where each vector q j is q j " R j Etyu, @j P v0, Kw.

Finally, according to the definition of y, we have Ety j u " δ j f 0 pv j q, Ety 2 j u " δ j f 2 0 pv j q and Ety i y j u " δ i δ j f 0 pv i qf 0 pv j q for i ‰ j, which would allow us to replace the above terms Etyu and Etyy J u with their values. For the experiments below, we will assume that the sampling strategy is unknown and estimate Epyq, Epyy J q empirically.

An exploration of modeling accuracy with different graph operators

The reconstruction quality using the previous approach depends on several factors: the choice of reference operator R, the order of the filter K and the random sampling strategy defined by the Bernoulli parameters δ k . Here, we evaluate primarily the influence of R and compare two possible sampling strategies. As an example, we also consider the dataset of the political blogs of the 2004 US presidential campaign [START_REF] Adamic | The political blogosphere and the 2004 us election: divided they blog[END_REF] described in Section 6.1.5.

As before, we will use the largest strongly connected directed subgraph G 1 , represented by its adjacency matrix W. On G 1 , we consider the following operators • W norm the normalized version of the adjacency matrix W whose spectral norm is equal to one.

• The random walk operator P associated with an ergodic random walk X with unique stationary distribution π.

• P the additive reversibilization of P.

• T " Π 1{2 PΠ ´1{2 the operator similar to P.

• T " Π 1{2 PΠ ´1{2 the operator similar to P.

• Pα P P with P " ␣ Pα : Pα " p1 ´αqP `αP ˚ˇα P r0, 1s (

• Tα P T with T " ␣ Tα : Tα " p1 ´αqΠ 1{2 PΠ ´1{2 `αΠ 1{2 P ˚Π´1{2 ˇˇα P r0, 1s (

We also consider G 1 sym , the symmetrized version of G 1 . The undirected graph G 1 sym is represented by its adjacency matrix W sym " pW `WJ q{2. On G 1 sym , we define the following operators • W norm sym the normalized version of the adjacency matrix W sym whose spectral norm is equal to one.

• The random walk operator P sym corresponding to the ergodic random walk X sym with stationary distribution π sym .

• T sym " Π 1{2 sym P sym Π ´1{2 sym the operator similar to P sym . Finally, we consider G, the full graph that is not strongly connected. Consequently, we can not build directly a random walk operator on G unless we transform G into a strongly connected graph, as required by our framework. We thus consider two classic approaches that make the graph strongly connected and the associated random walk ergodic: Approach 1: Rank-one perturbation: from the original adjacency matrix W, we construct a new adjacency matrix W ϵ as follows

W ϵ " W `ϵJ,
where J " 11 J {N is a rank-one matrix and ϵ is small. The weak perturbation of the adjacency matrix W by the matrix J ensures that the random walk on G is ergodic with stationary measure π ϵ and its associated transition matrix P ϵ is well-defined.

For our experiments, we choose ϵ " 10 ´4 to have }W ´Wϵ } F small enough.

Approach 2: Construction of the Google matrix of G [START_REF] Langville | Google's PageRank and beyond: The science of search engine rankings[END_REF][START_REF] Ermann | Google matrix analysis of directed networks[END_REF]. This is achieved in two steps. Firstly, we construct an adjacency matrix W from W by adding a weight one from all dangling nodes, that is nodes with no out-edges towards all the nodes in the graph. From W we can construct the transition matrix S. Secondly, we define the Google matrix P G as P G " γS `p1 ´γqJ.

where γ " 0.85 [START_REF] Langville | Google's PageRank and beyond: The science of search engine rankings[END_REF][START_REF] Ermann | Google matrix analysis of directed networks[END_REF].

We will compare the reconstruction accuracy using a proportion of correctly reconstructed labels obtained via f " sgnp p Hyq.

Subproblem 1. We solve the problem (39) by learning a polynomial graph filter with K " 10 on G 1 and G 1 sym with each of the reference operators listed above. In this subproblem, the random variables y j , j " 1, . . . , N 1 are distributed according to one of the two following cases • uniform sampling: y j " ε j f 1 0 pv j q where ε j " Berppq with p the proportion of known labels.

• π-weighted sampling: y j " ε j f 1 0 pv j q where ε j " Berpαπ j q, such that ř j Etε j u " pN 1 with p the proportion of known labels.

The proportions of correctly reconstructed labels are measured for various p values using 500 realizations of y to estimate its mean and covariance used in the solution [START_REF] Crovella | Graph wavelets for spatial traffic analysis[END_REF].

We summarize the different cases in table 2.

Case

Distribution y Reference operators

Case 1 y j " ε j f 1 0 pv j q, ε j " Berppq W norm , P, P, T, T, W norm sym , P sym , T sym Case 2 y j " ε j f 1 0 pv j q, ε j " Berpαπ j q, ř j Etε j u " pN 1 W norm , P, P, T, T, W norm sym , P sym , T sym Numerical simulations: Case 1. We evaluate the reconstruction performance of the graph signal f 1 0 in figure 5. For all proportions p of known labels, the average rate of correctly reconstructed labels using a graph filter baseed on P or P is significantly better than when using a graph filter based on W norm . This holds on the subgraph G 1 ands its symmetrized version G 1 sym . Furthermore, the reconstruction performance using a filter based on P is slightly better than using a filter based on P. The reconstruction performance using a filter based on T is the best overall. For small proportions p of known labels, the average rate of correctly reconstructed labels from filter based on T is slightly better. We also notice that the reconstruction performance using a filter based on T sym is identical to using a filter based on P sym for all proportions p of known labels. Numerical simulations: Case 2. The reconstruction performance is shown in figure 6. Filters based on W norm perform significantly worse than any other choice both on the subgraph G 1 or its symmetrized version G 1 sym . Also, the reconstruction accuracy when using W norm has a much larger variability and the average performance even decreases for larger p. Among other methods, the differences are not significant except for filters based on T which perform slightly worse than others for large p.

Case 1 vs. Case 2. The reconstruction performance using a filter based on P or P is slightly better when the samples in y are selected according to a distribution proportional to the stationary distribution π (case 2) than when the distribution is uniform (case 1). It is also the case when using a filter based on T or T. This is expected as vertices with larger π k values correspond to better-connected blogs, which are likely to have an influence on a larger number of other blogs. Although the reconstruction performance using a filter based on W norm is better in case 2 than in case 1, the large variability and the poorer performance at large p in case 2 suggest that these filters provide generally poorer models of the signal. The better reconstruction performance using a filter based on T or T with respect to using a filter based on P or P suggests that for this application it seems more suitable for the graph signal to live in ℓ 2 pVq than ℓ 2 pV, πq.

Subproblem 2. Here we consider the same problem [START_REF] Geller | Continuous wavelets on compact manifolds[END_REF] as in the previous section but now compare the performance using reference operators Pα P P and their equivalents Tα " Π 1{2 Pα Π ´1{2 . The random variables y j , j " 1, . . . , N 1 are distributed according to one of the two following cases:

• uniform sampling: y j " ε j f 1 0 pv j q where ε j " Berppq with p the proportion of known labels.

• π-weighted sampling: y j " ε j f 1 0 pv j q where ε j " Berpαπ j q, such that ř j Etε j u " pN 1 with p the proportion of known labels.

We summarize the different cases in table 3.

Case

Distribution y Reference operators Case 1 y j " ε j f 1 0 pv j q, ε j " Berppq Pα P P, Tα P T Case 2 y j " ε j f 1 0 pv j q, ε j " Berpαπ j q, ř j Etε j u " pN 1 Pα P P, Tα P T Numerical simulations: Case 1. We evaluate the reconstruction performance of f 1 0 using filters either based on Pα , α P r0, 1s at the top of fig. 7 or either based on Tα , α P r0, 1s at the bottom of fig. 7. We first consider the top of fig. 7. For small proportions p of known labels, we notice that the best reconstruction performance using a filter based on Pα is obtained for α » 0.2 while for higher proportions p, the best performance is obtained for α » 0.3. We notice poorer performance overall when α increases from α " 0.5 to α " 1. At the bottom of fig. 7, the reconstruction performances using a filter based on Tα , α » 0.2 are the best except for the rate p " 0.01. We also notice poorer performance overall when α increases from α " 0.5 to α " 1. For all proportions p of known labels, filters based on Tα appear to perform slightly better than filters based on Pα . This again suggests that this graph signal f 1 0 is better represented as living in ℓ 2 pVq than ℓ 2 pV, πq. Furthermore, the better performance when α P p0, 0.5q suggests that, for this signal, graph filters Pα perform better when the convex combination of P and P ˚involves more the random walk P, i.e. α P p0, 0.5q, than the time-reversed random walk P ˚.

Numerical simulations: Case 2. We evaluate the reconstruction performance of f 1 0 using filters either based on Pα , α P r0, 1s at the top of fig. 8 or based on Tα , α P r0, 1s at the bottom of fig. 8. Let us consider the top of fig. 8. We notice that the reconstruction performances using a filter based on Pα , α » 0.3 are the best except for the rates p " t0.01, 0.02u. We notice poorer performance overall when α increases from α " 0.5 to α " 1. This suggests that, for this signal, graph filters based on Pα perform better when the convex combination of P and P ˚involves more the forward random walk P, i.e. α P p0, 0.5q, than the backward random walk P ˚, i.e. α P p0.5, 1q. At the bottom of fig. 8, the reconstruction performances using a filter based on Tα , α » 0.3 are the best for all proportions p of known labels. Similarly to the previous cases, we notice poorer performance overall when α increases from α " 0.5 to α " 1.

On this example, the reconstruction performance is better using filters based on Pα than Tα . This suggests that the graph signal f 1 0 is best viewed as living in ℓ 2 pV, πq compared to ℓ 2 pVq when the samples in y are selected from a distribution proportional to π. 3. The choice of the parameter α in the convex combination between P and P ˚has an influence on the performance. For the example presented here, the optimal convex combination of P and P ˚(resp. T and T ˚) always involves more the forward random walk (P) than the backward random walk P ˚, except in one case. The message is that varying α should be considered and useful. A limitation of the present work is that we don't provide a way to choose it and one will have to rely classically on grid search or cross-validation. The optimal parameter is expected to depend upon the type of data studied (whether the information, hence the structure, flows more in the direction of the arrow, or the reverse). Let us remark that, in a specific case, it can be possible to include the estimation of α in the optimisation process, as it has been done in [START_REF] Hua | Learning combination of graph filters for graph signal modeling[END_REF].

Multiresolution analyses applications

This section is devoted to illustrate our multiresolution analysis framework through the visualization of wavelets and scaling functions on regular-type structures and semisupervised learning on digraphs. The application framework will be the same as in section 6.

Diffusion wavelets on the directed cycle graph

We show an example of multi-resolution analysis on the directed cycle graph C N with N " 256. We use the same assumptions made in Section 4.5.1. We construct both orthogonal and biorthogonal multi-resolution analyses on C N through the framework of the diffusion wavelets applied on the dyadic powers of T, i.e. tT 2 j u J j"1 . We set the number of scales at J " 6.

Figure 11 shows the orthogonal and biorthogonal scaling functions at scales 1,3 and 5. At each scale, we represent 3 or 4 scaling functions. Figure 12 shows the orthogonal and biorthogonal wavelet functions also at scales 1,3 and 5. We note at each scale that the support of the orthogonal scaling functions is larger than the support of the biorthogonal scaling functions. We also note on scales 3 and 5 that biorthogonal wavelet functions have support slightly smaller than the orthogonal wavelet functions. Furthermore, we note that orthogonal scaling functions, as well as orthogonal wavelet functions, have more oscillations than their biorthogonal counterparts. That means that orthogonal scaling and wavelet functions have poorer frequency localization. The conclusions about the spatial localization of the scaling functions are identical to the conclusions in [START_REF] Bremer | Diffusion wavelet packets[END_REF]. The novelty here is we work on the directed circle graph and we show the orthogonal and biorthogonal wavelets functions in addition to the orthogonal and biorthogonal scaling functions. Finally, the scaling and wavelets functions are not centered around the vertex where they have been selected. They instead propagate along with the graph towards the parent vertices of the selected vertex. We can also note that the orthogonal wavelet transform is more robust with respect to perturbations than the biorthogonal wavelet transform. Let W ortho " rΨ 1 , . . . , Ψ 5 , Φ 5 s be the inverse wavelet transform and W biorth its biorthogonal counterpart. The condition numbers on this example are κpW ortho q « 33 and κpW biortho q « 2 ˆ10 4 . Consequently, the biorthogonal wavelet transform is much more sensitive to perturbations than its orthogonal counterpart.

Spectral wavelets on the directed Watts-Strogatz graph

In this section, we show examples of multi-resolution analyses on the directed version of a graph from the Watts-Strogatz model based on the spectral wavelets construction describes at Section 5.1. The Watts-Strogatz model [START_REF] Watts | Collective dynamics of "small-world" networks[END_REF] is an undirected random graph model exhibiting small world properties including short average path lengths and high clustering [START_REF] Newman | The structure and function of complex networks[END_REF]. The construction of a digraph from the Watts-Strogatz model starts with a directed cycle graph with N vertices. Each node is connected to its k next nodes following the direction of the directed cycle graph. For the sake of simplicity, we consider the k next nodes connected to a given node i following the direction of the directed cycle graph as its "closest" neighbors. Starting from an arbitrary vertex, we apply the following procedure to each vertex in a clockwise manner. At vertex i, the edge that connects i to each of its next nodes is randomly rewired with probability p or remains untouched with probability 1 ´p. This procedure is repeated cyclically for each successive vertex until the vertex i is selected again. We denote by G " DWSpN, K, βq a graph constructed following a directed Watts-Strogatz model with N vertices, K nearest neighbors and rewiring probability β.

In a first instance, we analyze the scaling functions built from low-pass filters based on Tα at a given scale of G 1 " DWSp64, 2, 0q. That is a special case of the directed Watts-Strogatz model with no rewiring. In a second instance, we analyze some scaling functions built from low-pass filters based on Tα at a given scale of G 2 " DWSp64, 2, 0.02q. Lastly, we show both orthogonal and biorthogonal scaling and wavelets functions built following the construction in Section 5.2. We summarize the following cases in table 6. Case 1. In this case, we analyze some scaling functions built from low-pass filters based on Pα at scale t " 2 4 of a directed Watts-Strogatz graph G 2 " DWSp64, 2, 0q. We consider the following low-pass filters

Case

Graph

Tα " Π 1{2 Pα Π ´1{2 , @α P r0, 1s.

and we construct a collection of low-pass graph filters at a given scale following the construction presented in Section 5.1. More precisely, we build the filters H α as follows

H α " ÿ ωPω hptωqS ω,α ,
with t " 2 4 and hpxq " expp´xq. Figure 13 represents the 50 th scaling function h α,50 " H α δ 50 at scale t " 2 4 built from low-pass filters T α , for all α P r0, 1s according to the construction Section 5.1 on the directed Watts-Strogatz graph G " DWSp64, 2, 0q. The particularity of this case is that the following graph has no rewiring. Consequently, the associated adjacency matrix is circulant. All the low-pass filters T α admit the same discrete Fourier basis, and all corresponding frequencies are the same. Hence the associated scaling functions h α,50 are the same. We observe that the scaling function is centered around node 50 and has a symmetric shape. In this case, for nonsymmetric graph filters H α , the scaling functions have symmetric shapes and are centered around their given node. Case 2. In this case, we analyze some scaling functions built from low-pass filters based on Pα at scale 10 of G 2 " DWSp64, 2, 0.02q. We also consider the following low-pass filters Tα " Π 1{2 Pα Π ´1{2 , @α P t0, 0.5, 1u.

and we construct a collection of low-pass graph filters at a given scale following the construction presented in section 5.1. More precisely, we build the filters H α as follows

H α " ÿ ωPω hptωqS ω,α
with t " 2 4 and hpxq " expp´xq. We observe the 50 th scaling functions at scale t " 2 4 , that is h α,50 " H α δ 50 . The 50 th scaling function at the scale t " 2 4 , for different α, h α,50 , are represented in fig. 14. The graph G 2 admits a directed edge from node 50 to node 11 and a directed edge from node 30 to node 50. As α increases towards 1, we observe that the scaling function h α,50 propagates around the child node 11. This means that the more α increases, the more the influence of the backward random walk P ˚is important. For α " 0.5, the scaling function h α,50 diffuses both ways and is centered around the child node 11 and the parent node 30. Finally, as α goes to 0, we observe that h α,50 mostly propagates towards the parent node 30.

Case 3. Here, we consider orthogonal and biorthogonal scaling functions on a graph G " DWSp64, 2, 0.02q used on the Case 2. These orthogonal and biorthogonal scaling functions are built from the diffusion wavelets procedure. We start the procedure with the low-pass filter H " T 2 and we look at the dyadic powers of H i.e. tH 2 j u J j"1 . We set the number of scales at J " 5. We observe the orthogonal and biorthogonal scaling functions at node 49 at scale 3 obtained by the diffusion wavelets procedure and we compare them to the scaling function at node 49 built from a graph filter based on T at the scale t " 16.

Figure 15 shows the orthogonal and bi-orthogonal scaling functions at node 49 built using the diffusion wavelets framework and the scaling function at node 49 built using the spectral graph wavelets framework. We note that the biorthogonal scaling function and the scaling function built by the spectral graph wavelet framework are similar, i.e. well localized around nodes 11 and 30 and localized around 49 with large size support. By contrast, the orthogonal scaling function has poor spatial localization around the nodes 11 and 30.

Semi-supervised learning on directed graphs with ℓ 1 -regularization

We discuss the semi-supervised learning approach in the case of functions over digraphs with ℓ 1 regularization on the wavelet coefficients. We aim to show that the performance of the semi-supervised learning problem with ℓ 1 regularization is competitive compared to the existing approaches, e.g. the semi-supervised learning problem studied in section 6. We use the same notations as in section 6 and work on the subgraph G 1 . The ℓ 1 -regularized solution for the semi-supervised learning problem is w ˚" argmin w }ỹ ´MKw} 2 2 `λ}w} 1 , λ P R `.

The graph signal ỹ " My is the partially labeled graph signal with M " tm ij u t1ďi,jďN u is the mask operator, i.e. the diagonal matrix where m ii " 1 v i PO where O Ă V is the subset of known labels.

The matrix K " pH J , G 1 , ¨¨¨, G J q is the synthesis filter bank. If we set X " MK, the equation ( 44) can be rewritten as w ˚" argmin w }ỹ ´Xw} 2 2 `λ}w} 1 , λ P R.

The formulation ( 44) is identical to the formulation of the problem of signal restoration with redundant wavelet transforms in [START_REF] Selesnick | Signal restoration with overcomplete wavelet transforms: Comparison of analysis and synthesis priors[END_REF], except that K is the synthesis wavelet transform for functions defined over digraphs. Furthermore, previous approaches of semisupervised learning on undirected graphs have been investigated using overcomplete graph wavelets [START_REF] Shuman | Semi-supervised learning with spectral graph wavelets[END_REF] or critically sampled spline graph wavelets [START_REF] Ekambaram | Wavelet-regularized graph semi-supervised learning[END_REF]. The ℓ 1 -regularized synthesis semi-supervised learning problem is convex and can be solved efficiently using e.g. proximal splitting methods [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF][START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF][START_REF] Becker | Nesta: A fast and accurate first-order method for sparse recovery[END_REF]. From the solution w ˚of eq. ( 44), we define the restored signal f ˚as f ˚" signpKw ˚q.

We compare the performance of formulations [START_REF] Narang | Perfect reconstruction two-channel wavelet filter banks for graph structured data[END_REF] and [START_REF] Rustamov | Wavelets on graphs via deep learning[END_REF]. Formulation (44) requires a synthesis filter bank K. We choose a two scales synthesis graph filter bank K " pH 2 , G 1 , G 2 q based on a heat kernel construction [START_REF] Barlow | Random walks and heat kernels on graphs[END_REF] H j " e ´Lt j , t j " 2 j , j " 1, 2.

with L the directed Laplacian defined at [START_REF] Luxburg | A tutorial on spectral clustering[END_REF] and the graph bandpass filters G 1 , G 2 as

G 2 " H 1 ´H2 , G 1 " I ´H1 .
Semi-supervised learning on directed stochastic block model. We use the same parametrization setting used in Section 6.1.4 for the generation of the directed SBM graphs G 1 and G 2 . Figure 16 shows the performance of the semi-supervised learning formulation (44) on instances generated by directed SBM graphs G 1 and G 2 . The performance is obtained by averaging 200 realizations and the parameters λ (the ℓ 1 penalization) and γ (the ℓ 2 regularization) are selected following a grid search on a controlled example (with ground-truth) that give the best performances with λ P r0, 100s and γ P r0, 100s. In figure Section 6.1.4, the performance based on the ℓ 1 penalization of the spectral graph wavelets coefficients on digraphs is better than the ℓ 1 penalization based on the spectral graph wavelets coefficients build from the symmetrized normalized Laplacian and for both instances of graphs and for all proportions of known labels. Consequently, taking into account the direction of edges for the semi-supervised learning problem based on ℓ 1 penalization of the spectral graphs coefficients is useful, and the proposed construction allows it. Semi-supervised learning on political blogs. Figure 17 shows the performance of semisupervised learning formulation (44) on the largest strongly connected graph obtained from US political blogs. The performance is obtained by averaging 200 realizations and choosing the parameters λ and γ that maximize the performance. In Figure 17, the performance based on the ℓ 1 penalization of spectral graph wavelets coefficients on digraphs is better than the ℓ 1 penalization of the spectral graph wavelets coefficients built from the symmetrized normalized Laplacian from all percentage of known labels. As the result, the semi-supervised approach based on the ℓ 1 penalization of the spectral wavelets coefficients on digraphs is more effective than the ℓ 1 penalization of the spectral graph wavelets coefficients built from the symmetrized normalized Laplacian on this dataset.

Conclusion

We introduced a novel harmonic analysis on directed graphs. First, we proposed a frequency analysis for functions defined on directed graphs based on the eigenvectors of the random walk operator on a directed graph. From this Fourier-type frequency interpretation, we showed how to construct redundant wavelets on directed graphs as well as critically sampled wavelets by generalizing the diffusion wavelets framework. Finally, we illustrated our harmonic analysis through examples of semi-supervised learning and graph signal modeling on directed graphs and showed the relevance of our framework compared to existing approaches.

Figure 1 :

 1 Figure1: Eigenvalue distribution of the random walk matrix of a directed toroidal graph T 6,4 . Eigenvalues on the dotted green line: Repλ g q " ´0.25, @λ g P Λ g . Eigenvalues on the dotted orange line: Repλ o q " 0.25, @λ o P Λ o .

Figure 2 :

 2 Figure 2: Two-dimensional representation of the eigenvectors for the toroidal graph T 36,60 with a graph frequency equal to p7 ´?5q{8. Each row shows one eigenvector, with its real part on the left hand side and its imaginary part on the right hand side. Among the 8 eigenvectors with this frequency, we only show those corresponding to eigenvalues with a non-negative imaginary part (shown on the left hand side).

Figure 3 :

 3 Figure 3: Semi-supervised learning on G 1 (left) and G 2 (right). The notations L dir,norm and L norm correspond respectively to L and L sym .

Figure 4 :

 4 Figure 4: Semi-supervised learning on G 1 of the political blogs of the 2004 US presidential campaign. The notations L norm and L norm,sym correspond respectively to L and L sym .

Figure 5 :

 5 Figure 5: (Subproblem 1, Case 1) Left: Reconstruction of the graph signal f 1 0 on G 1 . Right: Reconstruction of the graph signal f 1 0 on G 1 sym .The notation P bar (resp. T bar ) correspond to P (resp. T).

  0

Figure 6 :

 6 Figure 6: (Subproblem 1, Case 2) Left: Reconstruction of the graph signal f 1 0 on the subgraph G 1 . Right: Reconstruction of the graph signal f 1 0 on G 1 sym . The notation P bar (resp. T bar ) correspond to P (resp. T).

Figure 8 :

 8 Figure 8: (Subproblem 2, case 2) Up: Reconstruction of the graph signal f 1 0 on G 1 , Pα P P. Down: Reconstruction of the graph signal f 1 0 on G 1 , Tα P T .

Figure 10 :Bi-orthogonal scaling functions Scale 5 Figure 11 :

 10511 Figure 10: (Subproblem 4, case 1) Up: Reconstruction of the graph signal f 0 on G, Pα P P. Down: Reconstruction of the graph signal f 0 on G, Tα P T .

Bi-orthogonal wavelet functions Scale 1 Figure 12 :

 112 Figure 12: Orthogonal and biorthogonal wavelet functions on the directed cycle graph C 256 .

Figure 13 :

 13 Figure 13: (Case 1) 50 th scaling function at scale t=16 on the graph G " DWSp64, 2, 0q, for all α P r0, 1s.

Figure 16 :Figure 17 :

 1617 Figure 16: Semi-supervised learning on G 1 (left) and G 2 (right).

Table 1 :

 1 Table of the different operators associated with the semi-supervised learning methods according to G and G sym .

	1							1
								0.95
	0.95							0.9
								0.85
	0.9							0.8				L
	Accuracy						Accuracy	0.75				L norm L dir L dir,norm
	0.85							0.7
								0.65
						L
	0.8					L norm		0.6
						L dir L dir,norm		0.55
	0.75	0	0.01 0.02	0.05	0.07 Rate of known labels	0.1	0.15	0.5	0.01 0.02	0.05	Rate of known labels 0.07	0.1	0.15

Table 2 :

 2 Summary table of the different cases of the subproblem 1.

				Graph signal reconstruction
		1			Subproblem 1 Case 1
		0.98
	Reconstruction accuracy	0.92 0.94 0.96							W sym norm
		0.9							P sym
									T sym
		0.88
		0	0.02	0.04	0.06	0.08	0.1	0.12	0.14	0.16
					Rate known labels

Table 3 :

 3 Summary table of the different cases of the subproblem 2.

  Reconstruction of the graph signal f 1 0 on G 1 , Pα P P. Down: Reconstruction of the graph signal f 1 0 on G 1 , Tα P T .

		Graph signal reconstruction					Graph signal reconstruction
		Subproblem 2, Case 2					Subproblem 2, Case 2
		P	α		1			P	α
		0.99
	Reconstruction accuracy	0.94 0.95 0.96 0.97 0.98	α=0 α=0.1 α=0.2 α=0.3 α=0.4	Reconstruction accuracy	0.92 0.94 0.96 0.98					α=0.5 α=0.6 α=0.7 α=0.8 α=0.9
			α=0.5							α=1
					0.9
	0.02 0.02 0.02 0.02 Figure 7: (Subproblem 2, Case 1) Up: 0 0.04 0.06 0.08 0.1 0.12 Rate known labels 0.91 0.92 0.93 0.94 0.95 0.96 0.97 0.98 0.99 Reconstruction accuracy Reconstruction graph signal Subproblem 2, Case 1 P α 0.04 0.06 0.08 0.1 0.12 0.14 0.14 α=0 α=0.1 α=0.2 α=0.3 α=0.4 α=0.5 0.16 0.16 Rate known labels 0.9 0.92 0.94 0.96 0.98 Reconstruction accuracy Reconstruction graph signal Subproblem 2, Case 1 T α α=0 α=0.1 α=0.2 α=0.3 α=0.4 α=0.5 0.04 0.06 0.08 0.1 0.12 0.14 0.16 Rate known labels 0 0.04 0.06 0.08 0.1 0.12 0.14 0.16 Rate known labels 0.92 0.99 1 Graph signal reconstruction Subproblem 2, Case 2 T α 0.94 0.95 0.96 0.97 0.98 Reconstruction accuracy α=0 α=0.1 α=0.2 α=0.3 α=0.4 0.93 α=0.5	Reconstruction accuracy Reconstruction accuracy Reconstruction accuracy	0.86 0.88 0.9 0.92 0.94 0.96 0.98 0.88 0.9 0.92 0.94 0.96 0.98 0.85 1 0.9 0.95	0 0	0.02 0.02 0.02 0.02	0.04 Reconstruction graph signal 0.06 0.08 0.1 Rate known labels Subproblem 2, Case 1 P α 0.04 0.06 0.08 0.1 0.12 Rate known labels Reconstruction graph signal Subproblem 2, Case 1 T α 0.04 0.06 0.08 0.1 Rate known labels 0.04 0.06 0.08 0.1 Rate known labels Graph signal reconstruction Subproblem 2, Case 2 T α	0.12 0.14 0.12 0.12	0.14 α=0.5 α=0.6 α=0.7 α=0.8 α=0.9 α=1 0.16 α=0.5 α=0.6 α=0.7 α=0.8 α=0.9 α=1 0.14 0.14 α=0.5 α=0.6 α=0.7 α=0.8 α=0.9 α=1	0.16 0.16 0.16

Table 6 :

 6 Case 2: Scaling functions built wrt. sect. 5.1 G 2 " DWSp64, 2, 0.02q Tα P T Case 3: Scaling functions built wrt. sect. 5.2 G 2 " DWSp64, 2, 0.02q T Table of the different cases for the directed Watts-Strogatz graph.

	Reference operators

Acknowledgements

We thank Yanis Bahroun, Randall Balestriero, Romain Cosentino, Benjamin Girault and Antonio Ortega for the useful discussions regarding this work. Last, we are grateful to the reviewers for their helpful suggestions.

‹ This work was supported by the ANR-14-CE27-0001 GRAPHSIP grant, the ACADEMICS Grant given by the IDEXLYON project of the Université de Lyon, as part of the "Programme Investissements d'Avenir" ANR-16-IDEX-0005, the ANR-19-CE48-0002 DARLING grant, and the GraphNEx CHIST-ERA ANR-21-CHR4-0009 grant.

‹‹ Preliminary results were presented in [1,2]. ˚Corresponding author.

Case

Distribution y

Reference operators Case 1 y j " ε j f 1 0 pv j q, ε j " Berppq W norm , P ϵ , P G , Pϵ , PG , T ϵ , T G , Tϵ , TG Previously, we observed that it is more suitable for the graph signal f 1 0 to belong to ℓ 2 pVq than ℓ 2 pV, πq when the distribution in y is uniform. This suggests that it may be a good idea in practice to use the measure of the Hilbert space to the sampling distribution if it is known. Subproblem 3. We now consider the same problem [START_REF] Geller | Continuous wavelets on compact manifolds[END_REF] on the whole graph G. In this subproblem, the random variables y j , j " 1 . . . , N are distributed as y j " ε j f 0 pv j q, ε j " Berppq. For each proportion of known labels p, the correct reconstruction rates are obtained by averaging over 500 realizations of y. We summarize the following case in table 4.

Numerical simulations. The reconstruction performance is shown in fig. 9. For all p values, the reconstruction performance using filters based on W norm is significantly worse than with all other reference operators. Compared to the other reference operators, W norm is the only one that can be used without modifying the graph to make it strongly connected. Still, its performance is always worse.

Among the other reference operators, we notice a better performance generally when using the reversibilizations Pϵ and PG compared to the non-reversible random walks P ϵ and P G . Furthermore, the best reconstruction performance is obtained using a filter based on TG and the reconstruction performance using a filter based on Tϵ is identical to TG for higher proportions of known labels. This differs from the results on G 1 where both would perform similarly.

Finally, the reconstruction performance does not depend as much on the approach we use to make the graph strongly connected. It seems though that the Google approach slightly outperforms the rank-one approach, both for the non-reversible random walks and their reversibilizations. Subproblem 4. Here we consider the same problem [START_REF] Geller | Continuous wavelets on compact manifolds[END_REF] as in the previous section but now compare the performance using reference operators Pα P P and their respective equivalents Tα " Π 1{2 Pα Π ´1{2 . The random walk operator is built from the Google matrix approach (Approach 1). The random variables y j , j " 1 . . . , N are distributed as y j " ε j f 0 pv j q, ε j " Berppq. We summarize the case in the table 5

Case

Distribution y Reference operators

Case 1 y j " ε j f 0 pv j q, ε j " Berppq Pα P P, Tα P T Numerical simulations: Case 1. We evaluate the reconstruction performance of f 1 0 using filters either based on Pα , α P r0, 1s at the top of fig. 10 or either based on Tα , α P r0, 1s at the bottom of fig. 10. For all proportions p of known labels, the best reconstruction performance using a filter based on Pα is obtained for α » 0.6. We also notice that the reconstruction performance using a filter with α P t0.7, 0.8u is better than with α » 0.5. This suggests that, for this signal, graph filters based on Pα perform better when the convex combination of P and P ˚involves more the time-reversed random walk P ˚than the random walk P (i.e. when α ą 0.5). At the bottom of fig. 10, the best performance is obtained for a filter Tα with α » 0.5 overall. The performance is worse when α increases from α " 0.5 to α " 1. For all proportions p of known labels, the reconstruction performance are better overall if the graph signal f 0 is in ℓ 2 pVq than in ℓ 2 pV, πq. In other words, the reconstruction performance using filters based on Tα is better than using filters based on Pα .

Key takeaways

The subproblems investigated and the numerical simulations allow us to note the following important points:

1. Both the Hilbert space where we choose for the graph signal and the sampling distribution impact the performance. For the example presented here, we obtain better performance when we consider the Hilbert space with inner product associated with a given measure corresponds to the sampling probability.

2. Graph filters based on the random walk operator yield better performance than those based on the adjacency matrix.