In the framework of "Discrete Signal Processing on Graphs" [1]}, [2]}, the semi-supervised learning problem for directed graphs is formulated as follows
\(\underset{f}{\operatorname{argmin}} \big \lbrace c\Vert \mathbf {M}_{l}(f-y)\Vert ^2+\varrho _2\Vert f-\mathbf {W}^{norm}f\Vert _2^2\big \rbrace , \quad c,\varrho _2 > 0,\) 
with \(\mathbf {W}^{norm}\) , a normalized version of the the adjacency matrix \(\mathbf {W}\)  so that its spectral norm is equal to one.
Semi-supervised learning on \(\mathcal {G}_{\operatorname{sym}}\)
We use the same methods for the semi supervised learning on \(\mathcal {G}\) .
If the graph signal \(y\in \ell ^2(\mathcal {V})\) , the optimization problem is
\(\underset{f}{\operatorname{argmin}}\big \lbrace c\Vert \mathbf {M}_{l}(f-y)\Vert ^2+c\Vert \mathbf {M}_{u}f\Vert ^2+\varrho _2\langle f,\mathcal {L}_{\operatorname{sym}} f\rangle \big \rbrace , \quad c,\varrho _2 > 0.\) 
with \(\mathcal {L}_{\operatorname{sym}}\)  the normalized Laplacian defined as
\(\mathcal {L}_{\operatorname{sym}}=\mathbf {I}-\Pi _{\operatorname{sym}}^{1/2}\mathbf {P}_{sym}\Pi _{\operatorname{sym}}^{-1/2}, \quad \Pi _{\operatorname{sym}}=\operatorname{diag}\lbrace \pi _{sym}(v_1),\ldots ,\pi _{sym}(v_N)\rbrace .\) 

The optimization problem is
\(\underset{\tilde{f}}{\operatorname{argmin}}\big \lbrace c\Vert \mathbf {M}_{l}(\tilde{f}-\tilde{y})\Vert ^2+c\Vert \mathbf {M}_{u}\tilde{f}\Vert ^2+\varrho _2\langle \tilde{f},\mathcal {L}_{\operatorname{sym}} \tilde{f}\rangle \big \rbrace , \quad c,\varrho _2 > 0.\) 
We can rewrite the optimization problem (REF ) as follows
\(\underset{f}{\operatorname{argmin}}\big \lbrace c\Vert \mathbf {M}_{l}(f-y)\Vert _{\pi }^2+c\Vert \mathbf {M}_{u}f\Vert _{\pi }^2+\varrho _2\langle f,\mathbf {L}_{\operatorname{RW,\operatorname{sym}}} f\rangle _{\pi }\big \rbrace , \quad c,\varrho _2 > 0.\) 
with \(\mathbf {L}_{\operatorname{RW,\operatorname{sym}}}\)  the random walk Laplacian defined as
\(\mathbf {L}_{\operatorname{RW,\operatorname{sym}}}=\mathbf {I}-\mathbf {P}_{\operatorname{sym}}.\) 
Solution of the problem
The problem (REF ) is quadratic and convex and therefore admits a closed form solution.
Semi-supervised learning on \(\mathcal {G}\)
The solution of the optimization problem (REF ) is
\(f^{*}=\operatorname{sign}\big [(\mathbf {I}+\gamma \mathcal {L})^{-1}\mathbf {M}_{l}y\big ],\quad \gamma =\varrho _2/c, \quad f^{*}\in \ell ^2(\mathcal {V}).\) 

The solution of the optimization problem (REF ) is
\(f^{*}=\operatorname{sign}\big [(\mathbf {I}+\gamma \mathbf {L}_{\operatorname{RW}})^{-1}\mathbf {M}_{l}y\big ],\quad \gamma =\varrho _2/c, \quad f^{*}\in \ell ^2(\mathcal {V},\pi ).\) 

The solution of the optimization problem (REF ) is
\(f^{*}=\operatorname{sign}\big [\big (\mathbf {M}_{l}+\gamma \mathbf {R}_{M}^{-1}\big )y\big ],\quad \gamma =\varrho _2/c, \quad f^{*}\in \ell ^2(\mathcal {V}).\) 
with \(\mathbf {R}_{M}=(\mathbf {I}-\mathbf {W}^{norm})^{\top }(\mathbf {I}-\mathbf {W}^{norm})\) .
Semi-supervised learning on \(\mathcal {G}_{\operatorname{sym}}\)
The solution of the optimization problem (REF ) is
\(f^{*}=\operatorname{sign}\big [(\mathbf {I}+\gamma \mathcal {L}_{\operatorname{sym}})^{-1}\mathbf {M}_{l}y\big ],\quad \gamma =\varrho _2/c, \quad f^{*}\in \ell ^2(\mathcal {V}).\) 

The solution of the optimization problem (REF ) is
\(f^{*}=\operatorname{sign}\big [(\mathbf {I}+\gamma \mathbf {L}_{\operatorname{RW,sym}})^{-1}\mathbf {M}_{l}y\big ],\quad \gamma =\varrho _2/c, \quad f^{*}\in \ell ^2(\mathcal {V},\pi _{\operatorname{sym}}).\) 
We summarize the different methods in Table REF .
<TABLE>Experiments
Let us consider the dataset of the political blogs of the 2004 US presidential campaign [3]}. The dataset consists of 1224 political blogs where each political blog is associated with a political orientation, either republican or democrat. This dataset is modeled by a graph \(\mathcal {G}=(\mathcal {V},\mathcal {E})\)  where each vertex \(v\in \mathcal {V}\)  is associated to a blog and an edge between two vertices \(\lbrace v_i,v_j\rbrace \)  indicates the presence of hyperlinks from the blog associated to the vertex \(v_i\)  to the vertex \(v_j\) . The political orientations of blogs are modeled by a graph signal \(f_0=\lbrace f_0(v_1),\cdots ,f_0(v_N)\rbrace \)  where \(f_0(v_i)\in \mathcal {O}\)  where each vertex is associated with a label belongs to the set \(\mathcal {O}=\lbrace -1,1\rbrace \)  characterizing its political orientation, the political orientation associated with Democrats being “+1" and the political orientation associated with Republicans being “-1". The directed graph \(\mathcal {G}\)  associated to the political blogs is not strongly connected. Our theoretical framework leads us to only consider strongly connected graphs. As a result, we consider the largest strongly connected subgraph of \(\mathcal {G}\) , denoted by \(\mathcal {G}^{\prime }=(\mathcal {V}^{\prime },\mathcal {E}^{\prime })\)  which is made up of \(|\mathcal {V}^{\prime }|=N^{\prime }=793\)  vertices (hence roughly 65% of the vertex set \(\mathcal {V}\) ) and its associated graph signal \(f_0^{\prime }\) . The respective performances of the approaches semisupdirzhu,semisupdirre,semisupmoura and their symmetrized counterparts semisupdirzhusym,semisupdirresym are compared in fig:ssllscc.

Numerical simulations
<FIGURE>Figure REF  displays the performance of the semi-supervised approaches semisupdirzhu,semisupdirre,semisupmoura,semisupdirzhusym,semisupdirresym. The performance rate is an average of 500 realizations calculated by determining the parameter \(\gamma \in [0,10]\)  on the training set giving the best performance rate.
Firstly, we remark that for all proportions of known labels, the performance rates of semi supervised approaches based on \(\mathcal {L}, \mathbf {L}_{\operatorname{RW}}\)  and its symmetrized counterparts \(\mathcal {L}_{\operatorname{sym}}, \mathbf {L}_{\operatorname{RW,sym}}\)  are better than the performance rates of the semi-supervised approach based on \(\mathbf {R}_{M}\) . On the top left of fig:ssllscc, we compare the semi-supervised approaches involving the operators \(\mathcal {L}, \mathbf {L}_{\operatorname{RW}}\)  and \(\mathbf {R}_M\) . The performance rate of the semi-supervised approach associated with \(\mathcal {L}\)  is the best compared to the approaches involving \(\mathbf {L}_{\operatorname{RW}},\mathbf {R}_M\)  for low rates of known labels. Furthermore, the performance rate approach associated with \(\mathcal {L}\)  is slightly better than the approach involving \(\mathbf {L}_{\operatorname{RW}}\) . On the top right of  fig:ssllscc, we compare the semi-supervised approaches involving the operators \(\mathcal {L}_{\operatorname{sym}}, \mathbf {L}_{\operatorname{RW,sym}}\)  and \(\mathbf {R}_M\) . Similarly the top left figure, the performance rate of the semi-supervised approach associated with \(\mathcal {L}_{\operatorname{sym}}\)  is slightly better than the approach involving \(\mathbf {L}_{\operatorname{RW}}\) . On the bottom left of fig:ssllscc, we compare the semi-supervised approaches involving the operators \(\mathbf {L}_{\operatorname{RW}},\mathbf {L}_{\operatorname{RW,sym}}\)  and \(\mathbf {R}_{M}\) . The performance rates between the approach involving \(\mathbf {L}_{\operatorname{RW}}\)  and \(\mathbf {L}_{\operatorname{RW,sym}}\)  are quite identical for any proportion of known labels as well as to the bottom right of fig:ssllscc comparing the semi-supervised approaches involving \(\mathbf {L}_{\operatorname{RW}},\mathbf {L}_{\operatorname{RW,sym}}\)  and \(\mathbf {R}_{M}\) . The comparison of performance rates of the semi-supervised approaches involving the different operators listed in the sslmethods associated with this example leads us to the following conclusions :

The semi supervised learning approaches involving Dirichlet energies based upon \(\mathcal {L}\)  and \(\mathbf {L}_{\operatorname{RW}}\)  (and their symmetrized counterparts) yields better performance rates compared to the approach proposed by Sandryhaila and Moura [1]} based upon \(\mathbf {R}_{M}\) .
Consequently, semi-supervised learning approaches on directed graphs with Dirichlet energies based on the random walk operator seem more appropriate.

The Hilbert space where we consider the graph signal seems to have an influence on the performance rates of the semi-supervised problem. Here, considering the graph signal \(y\)  belongs to \(\ell ^2(\mathcal {V})\) , or equivalently considering the operator \(\mathcal {L}\)  in the regularization term yields to better performance than considering the graph signal \(y\)  belonging to \(\ell ^2(\mathcal {V},\pi )\)  or equivalently considering the operator \(\mathcal {L}_{\operatorname{RW}}\)  in the regularization term.


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 \(\mathcal {G}=(\mathcal {V},\mathcal {E})\)  be a directed graph with cardinality \(|\mathcal {V}|=N\)  and \(\mu :\mathcal {V}\rightarrow \mathbb {R}_{+}\)  a positive measure. The directed graph \(\mathcal {G}\)  is characterized by its adjacency matrix \(\mathbf {W}\in \mathbb {R}^{N\times N}\) . Given a random graph signal \(y\) , we wish to determine the missing values of the signal by using a graph filter \(\mathbf {H}\)  expressed as a polynomial finite sum of a reference operator \(\mathbf {R}\) , that is
\(\mathbf {H}=\sum _{k=0}^K\theta _k\mathbf {R}^{k},\quad \theta _k\in \mathbb {R},\quad k=0,\cdots ,K.\) 
The whole graph signal \(f_0\in \ell ^2(\mathcal {V},\mu )\)  is defined as
\(f_0=[f_0(v_1),\ldots ,f_0(v_N)]^{\top }.\) 
The random signal \(y=[y_1,\ldots ,y_N]^{\top }\)  is a multivariate random variable indexed on the vertex set \(\mathcal {V}\) . The random variables \(y_j, j=1,\ldots ,N\)  can be expressed as \(y_j=\varepsilon _j f_0(v_j)\) , where the \(\varepsilon _j\)  are assumed independent and Bernoulli distributed with parameters \(\delta _j\) : \(\varepsilon _j\sim \operatorname{Ber}(\delta _j)\) . The random variable \(\bar{\varepsilon }\)  defined as
\(p=\frac{1}{N}\sum _{i=1}^N \varepsilon _i,\) 
is the proportion of known values of \(f_0\) , all collected in \(y\) . This signal \(y\)  can be viewed as a random sampling of the graph signal \(f_0\) . The optimization problem is the following
\(\underset{\theta =\lbrace \theta _k\rbrace _{k=0}^K\in \mathbb {R}^{K+1}}{\operatorname{argmin}}\mathbb {E}\bigg [\Vert f_0-\sum _{k=0}^K\theta _k\mathbf {R}^ky\Vert _{\mu }^2 \bigg ],\) 

Solution of the problem
The solution of the problem (REF ) is the following
\(\widehat{\theta }=\mathbf {Z}^{-1}\mathbf {M}\mathbf {Q}^{\top }f_0,\) 
where \(\mathbf {M}=\operatorname{diag}\lbrace \mu (v_1),\cdots ,\mu (v_N)\rbrace \)  is the diagonal matrix of the associated measure \(\mu \) , \(\mathbf {Z}\)  the matrix where each entry \(Z_{k\ell }\in \mathbf {Z}\)  is expressed as
\(Z_{k\ell }=\operatorname{Tr}\big ([\mathbf {R}^{k}]^{\top }\mathbf {M}\mathbf {R}^{\ell }\mathbb {E}(y y^{\top })\big ),\quad \forall \lbrace k,\ell \rbrace  \in 0,K ^{2},\) 
and \(\mathbf {Q}=[q_{0},\cdots ,q_{K}]\in \mathbb {R}^{N\times (K+1)}\)  where each vector \(q_{j}\)  is
\(q_{j}=\mathbf {R}^{j}\mathbb {E}(y),\quad \forall j\in 0,K .\) 
Consequently, the resulting graph filter \(\widehat{\mathbf {H}}\)  is
\(\widehat{\mathbf {H}}=\sum _{k=0}^K\hat{\theta }_k\mathbf {R}^{k},\quad \hat{\theta }_k\in \mathbb {R},\quad k=0,\cdots ,K.\) 

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 \(\mathbf {R}\) , the order of the filter \(K\)  and the random sampling strategy defined by the Bernoulli parameters \(\delta _k\) . Here, we evaluate primarily the influence of \(\mathbf {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 [3]}, described in section REF .
We consider the directed graph \(\mathcal {G}^{\prime }\) . The directed graph \(\mathcal {G}^{\prime }\)  is represented by its adjacency matrix \(\mathbf {W}\) . On \(\mathcal {G}^{\prime }\) , we consider the following operators

\(\mathbf {W}_{norm}\)  the normalized version of the adjacency matrix \(\mathbf {W}\)  whose norm is equal to one.

The random walk operator \(\mathbf {P}\)  associated with an ergodic random walk \(\mathcal {X}\)  with unique stationary distribution \(\pi \) .

\(\bar{\mathbf {P}}\)  the additive reversibilization of \(\mathbf {P}\) .

\(\mathbf {T}=\Pi ^{1/2}\mathbf {P}\Pi ^{-1/2}\)  the operator similar to \(\mathbf {P}\) .

\(\bar{\mathbf {T}}=\Pi ^{1/2}\bar{\mathbf {P}}\Pi ^{-1/2}\)  the operator similar to \(\bar{\mathbf {P}}\) .

\(\bar{\mathbf {P}}_{\alpha } \in \bar{\mathcal {P}}\)  with \(\bar{\mathcal {P}}=\big \lbrace \bar{\mathbf {P}}_{\alpha }: \bar{\mathbf {P}}_{\alpha }=(1-\alpha )\mathbf {P}+\alpha \mathbf {P}^{*}\big | \alpha \in [0,1]\big \rbrace \)

\(\bar{\mathbf {T}}_{\alpha } \in \bar{\mathcal {T}}\)  with \(\bar{\mathcal {T}}=\big \lbrace \bar{\mathbf {T}}_{\alpha }: \bar{\mathbf {T}}_{\alpha }=(1-\alpha )\Pi ^{1/2}\mathbf {P}\Pi ^{-1/2}+\alpha \Pi ^{1/2}\mathbf {P}^{*}\Pi ^{-1/2}\big | \alpha \in [0,1]\big \rbrace \)

We also consider \(\mathcal {G}_{\operatorname{sym}}^{\prime }\) , the symmetrized version of \(\mathcal {G}^{\prime }\) . The undirected graph \(\mathcal {G}_{\operatorname{sym}}^{\prime }\)  is represented by its adjacency matrix \(\mathbf {W}_{\operatorname{sym}}=(\mathbf {W}+\mathbf {W}^{\top })/2\) . On \(\mathcal {G}_{\operatorname{sym}}^{\prime }\) , we define the following operators

\(\mathbf {W}_{\operatorname{sym}}^{norm}\)  the normalized version of the adjacency matrix \(\mathbf {W}_{\operatorname{sym}}\)  whose the norm is equal to one.

The random walk operator \(\mathbf {P}_{\operatorname{sym}}\)  corresponding to the ergodic random walk \(\mathcal {X}_{\operatorname{sym}}\)  with stationary distribution \(\pi _{\operatorname{sym}}\) .

\(\mathbf {T}_{\operatorname{sym}}=\Pi _{\operatorname{sym}}^{1/2}\mathbf {P}_{\operatorname{sym}}\Pi _{\operatorname{sym}}^{-1/2}\)  the operator similar to \(\mathbf {P}_{\operatorname{sym}}\) .

Finally, we consider \(\mathcal {G}\) , the full graph. The directed graph \(\mathcal {G}\)  is represented by its adjacency matrix \(\mathbf {W}\) . The graph \(\mathcal {G}\)  is not strongly connected. Consequently, we can not build directly a random walk operator on \(\mathcal {G}\) . We need to deal with its property of not being strongly connected, as required by our framework.
Indeed, some vertices in \(\mathcal {G}\)  have an out-degree equal to zero and we cannot create a transition matrix directly from the adjacency matrix. To overcome this problem, we propose two approaches that make the graph strongly connected and the associated random walk ergodic:
Rank-one perturbation: from the original adjacency matrix \(\mathbf {W}\) , we construct a new adjacency matrix \(\mathbf {W}_{\epsilon }\)  as follows
\(\mathbf {W}_{\epsilon }=\mathbf {W}+\epsilon \mathbf {J},\) 
where \(\mathbf {J}=\mathbf {11}^{\top }/N\)  is a rank-one matrix and \(\epsilon \)  is small. The weak perturbation of the adjacency matrix \(\mathbf {W}\)  by the matrix \(\mathbf {J}\)  ensures that the random walk on \(\mathcal {G}\)  is ergodic with stationary measure \(\pi _{\epsilon }\)  and its associated transition matrix \(\mathbf {P}_{\epsilon }\)  is well-defined. For our experiments, we choose \(\epsilon =10^{-4}\) .

Construction of the Google matrix of \(\mathcal {G}\)  [6]}. This is achieved in two steps. Firstly, we construct an adjacency matrix \(\tilde{\mathbf {W}}\)  from \(\mathbf {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 \(\tilde{\mathbf {W}}\)  we can construct the transition matrix \(\mathbf {S}\) . Secondly,
we define the Google matrix \(\mathbf {P}_{G}\)  as
\(\mathbf {P}_{G}=(1-\gamma )\mathbf {S}+\gamma \mathbf {J}.\) 
where \(\gamma =0.85\)  [6]}.
We will compare the reconstruction accuracy using a proportion of correctly reconstructed labels obtained via
\(\hat{f}=\operatorname{sign}(\widehat{\mathbf {H}} y).\) 
Subproblem 1
We solve the problem (REF ) by learning a polynomial graph filter with \(K=10\)  on \(\mathcal {G}^{\prime }\)  and \(\mathcal {G}_{\operatorname{sym}}^{\prime }\)  We consider as reference operators the ones listed on \(\mathcal {G}^{\prime }\)  and \(\mathcal {G}_{\operatorname{sym}}^{\prime }\) . In this subproblem, the random variables \(y_j,j=1,\ldots ,N^{\prime }\)  are distributed according to one of the two following cases

Random variables \(y_j=\varepsilon _jf_{0}^{\prime }(v_j)\)  where \(\varepsilon _j\sim \operatorname{Ber}(p)\)  with \(p\)  the proportion of known labels.

Random variables \(y_j=\varepsilon _j f_{0}^{\prime }(v_j)\)  where \(\varepsilon _j\sim \operatorname{Ber}(\alpha \pi _j)\) , such that \(\sum _j\mathbb {E}(\varepsilon _j)=pN^{\prime }\)  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  (REF ).
We summarize the different cases in table REF .
<TABLE><FIGURE>
Numerical simulations: Case 1
We evaluate the reconstruction performance of the graph signal \(f_0^{\prime }\)  in figure REF .
For all proportions \(p\)  of known labels, the average rate of correctly reconstructed labels in the graph signal \(f_0^{\prime }\)  from a learned polynomial graph filter of \(\mathbf {P}\)  or \(\bar{\mathbf {P}}\)  is significantly better than the one based on a learned polynomial graph filter of \(\mathbf {W}_{norm}\)  in the case where we consider either the subgraph \(\mathcal {G}^{\prime }\)  or its symmetrized version \(\mathcal {G}_{\operatorname{sym}}^{\prime }\) . Furthermore, the reconstruction performance using a filter based on \(\mathbf {P}\)  is slightly better than using a filter based on \(\bar{\mathbf {P}}\) . The reconstruction performance using a filter based on \(\mathbf {T}\)  is the best. For small proportions \(p\)  of known labels, the average rate of correctly reconstructed labels from a learned polynomial graph filter of \(\bar{\mathbf {T}}\)  is slightly better. We also notice that the reconstruction performance using a filter based on \(\mathbf {T}_{\operatorname{sym}}\)  is identical to using a filter based on \(\mathbf {P}_{\operatorname{sym}}\)  for all proportions \(p\)  of known labels.
<FIGURE>
Numerical simulations: Case 2
The reconstruction performance is shown in figure REF . We observe a significantly better performance using a filter based on \(\mathbf {P}\)  or \(\bar{\mathbf {P}}\)  than based on \(\mathbf {W}_{norm}\)  in the case where we consider either the subgraph \(\mathcal {G}^{\prime }\)  or its symmetrized version \(\mathcal {G}_{\operatorname{sym}}^{\prime }\)  and a slightly better performance using \(\mathbf {P}\)  rather than \(\bar{\mathbf {P}}\) . For all proportions \(p\)  of known labels, the average rates of correctly reconstructed labels are identical using either a filter based on \(\mathbf {T}\)  or \(\bar{\mathbf {T}}\) 
are better than using a filter based on \(\mathbf {P}\)  or \(\bar{\mathbf {P}}\) . We notice the reconstruction performance based on \(\mathbf {W}_{norm}\)  has a large variability and even decreases for larger \(p\) .

Case 1 vs. Case 2
The reconstruction performance using a filter based on \(\mathbf {P}\)  or \(\bar{\mathbf {P}}\)  is slightly better when the samples in \(y\)  are selected according to a distribution proportional to the stationary distribution \(\pi \)  (case 2) than when the distribution is uniform (case 1). It is also the case by using a filter based on \(\mathbf {T}\)  or \(\bar{\mathbf {T}}\) . This is expected as vertices with larger \(\pi _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 \(\mathbf {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 \(\mathbf {T}\)  or \(\bar{\mathbf {T}}\)  with respect to using a filter based on \(\mathbf {P}\)  or \(\bar{\mathbf {P}}\)  suggest that the graph signal seems more suitable to live in \(\ell ^2(\mathcal {V})\)  than \(\ell ^2(\mathcal {V},\pi )\) .

Subproblem 2
Here we consider the same problem (REF ) as in the previous section but now compare the performance using reference operators \(\bar{\mathbf {P}}_{\alpha }\in \bar{\mathcal {P}}\)  and their equivalents \(\bar{\mathbf {T}}_{\alpha }=\Pi ^{1/2}\bar{\mathbf {P}}_{\alpha }\Pi ^{-1/2}\)  The random variables \(y_j, j=1,\ldots ,N^{\prime }\)  are distributed according to one of the two following cases:

Random variables \(y_j=\varepsilon _jf_{0}^{\prime }(v_j)\)  where \(\varepsilon _j\sim \operatorname{Ber}(p)\)  with \(p\)  the proportion of known labels.

Random variables \(y_j=\varepsilon _j f_{0}^{\prime }(v_j)\)  where \(\varepsilon _j\sim \operatorname{Ber}(\alpha \pi _j)\) , such that \(\sum _j\mathbb {E}(\varepsilon _j)=pN^{\prime }\)  with \(p\)  the proportion of known labels.

We summarize the different cases in subproblemtwo.
<TABLE><FIGURE>
Numerical simulations: Case 1
We evaluate the reconstruction performance of \(f_0^{\prime }\)  using filters either based on \(\bar{\mathbf {P}}_{\alpha }, \alpha \in [0,1]\)  at the top of fig:recsignalcaseonsubptwo or either based on \(\bar{\mathbf {T}}_{\alpha }, \alpha \in [0,1]\)  at the bottom of fig:recsignalcaseonsubptwo. We first consider the top of fig:recsignalcaseonsubptwo. For small proportions of \(p\)  of known labels, we notice that the reconstruction performance using a filter based on \(\bar{\mathbf {P}}_{\alpha }, \alpha \simeq 0.2\)  is the best while for higher proportions \(p\) , the reconstruction performance using a filter based on \(\bar{\mathbf {P}}_{\alpha },\alpha \simeq 0.3\)  is the best. We notice poorer performance overall when \(\alpha \)  increases from \(\alpha =0.5\)  to \(\alpha =1\) . At the bottom of fig:recsignalcaseonsubptwo, the reconstruction performances using a filter based on \(\bar{\mathbf {T}}_{\alpha }, \alpha \simeq 0.2\)  are the best except for the rate \(p=0.01\) . We notice poorer performance overall when \(\alpha \)  increases from \(\alpha =0.5\)  to \(\alpha =1\) .
For all proportions \(p\)  of known labels, the reconstruction performance are slightly better overall if we consider filters based on \(\bar{\mathbf {T}}_{\alpha },\alpha \in [0,1]\)  than filters based on \(\bar{\mathbf {P}}_{\alpha },\alpha \in [0,1]\) . That means the problem is more adapted when the graph signal \(f_0^{\prime }\)  is in \(\ell ^2(\mathcal {V})\)  than \(f_0^{\prime }\)  is in \(\ell ^2(\mathcal {V},\pi )\) . Furthermore, the study of this case suggests that, for this signal, graph filters \(\bar{\mathbf {P}}_{\alpha }\)  perform better when the convex combination of \(\mathbf {P}\)  and \(\mathbf {P}^*\)  involves more the random walk \(\mathbf {P}\) , i.e. \(\alpha \in (0,0.5)\) , than the random \(\mathbf {P}^{*}\) , i.e. \(\alpha \in (0.5,1)\) .
<FIGURE>
Numerical simulations: Case 2
We evaluate the reconstruction performance of \(f_0^{\prime }\)  using filters either based on \(\bar{\mathbf {P}}_{\alpha }, \alpha \in [0,1]\)  at the top of fig:recsignalcasetwosubptwo or either based on \(\bar{\mathbf {T}}_{\alpha }, \alpha \in [0,1]\)  at the bottom of fig:recsignalcasetwosubptwo. Let us consider the top of fig:recsignalcasetwosubptwo. We notice that the reconstruction performances using a filter based on \(\bar{\mathbf {P}}_{\alpha },\alpha \simeq 0.3\)  are the best except for the rates \(p=\lbrace 0.01,0.02\rbrace \) . We notice poorer performance overall when \(\alpha \)  increases from \(\alpha =0.5\)  to \(\alpha =1\) . This suggests that, for this signal, graph filters based on \(\bar{\mathbf {P}}_{\alpha }\)  perform better when the convex combination of \(\mathbf {P}\)  and \(\mathbf {P}^*\)  involves more the forward random walk \(\mathbf {P}\) , i.e. \(\alpha \in (0,0.5)\) , than the backward random walk \(\mathbf {P}^*\) , i.e. \(\alpha \in (0.5,1)\) . At the bottom of fig:recsignalcasetwosubptwo, the reconstruction performances using a filter based on \(\bar{\mathbf {T}}_{\alpha }, \alpha \simeq 0.3\)  are the best for all proportions \(p\)  of known labels. Similarly to the previous cases, we notice poorer performance overall when \(\alpha \)  increases from \(\alpha =0.5\)  to \(\alpha =1\) .
Firstly, the reconstruction performance is better using filters based on \(\bar{\mathbf {P}}_{\alpha }\)  than \(\bar{\mathbf {T}}_{\alpha }\) . That means the problem is more suitable for graph signal \(f_0^{\prime } \in \ell ^2(\mathcal {V},\pi )\)  is slightly better than performance with respect to the graph signal belonging to \(\ell ^2(\mathcal {V})\)  when the samples in \(y\)  are selected from a distribution proportional to the stationary distribution \(\pi \)  than when the distribution is uniform. Conversely, the problem is more suitable for graph signal \(f_0^{\prime }\)  is \(\ell ^2(\mathcal {V})\)  than \(\ell ^2(\mathcal {V},\pi )\)  when the distribution in \(y\)  is uniform while the samples in \(y\)  are selected according to a distribution proportional to the stationary distribution \(\pi \) . As a result, the distribution of the graph signal has an influence in the reconstruction performance through the use of learned filters based either on \(\bar{\mathbf {P}}_{\alpha }\)  or \(\bar{\mathbf {T}}_{\alpha }\) .

Subproblem 3
We now consider the same problem (REF ) on the whole graph \(\mathcal {G}\) . In this subproblem, the random variables \(y_j,j=1\ldots ,N\)  are distributed as \(y_j=\varepsilon _jf_{0}(v_j), \varepsilon _j\sim \operatorname{Ber}(p)\) . For each proportion of known labels \(p\) , the good reconstruction rates are derived by averaging 500 graph signals realizations of \(y\) .
We summarize the following case in the subproblemthree.
<TABLE>
Numerical simulations
The reconstruction performance is shown in fig:recsignalcaseonesubpthree. For all \(p\)  values, the reconstruction performance using filters based on \(\mathbf {W}_{norm}\)  is significantly worse than with all other reference operators. Compared to the other reference operators, \(\mathbf {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 clearly better performance generally when using the reversibilizations \(\bar{\mathbf {P}}_{\epsilon }\)  and \(\bar{\mathbf {P}}_{G}\)  compared to the non-reversible random walks \({\mathbf {P}}_{\epsilon }\)  and \({\mathbf {P}}_{G}\) . Furthermore, the best reconstruction performance is obtained using a filter based on \(\bar{\mathbf {T}}_{G}\)  and the reconstruction performance using a filter based on \(\bar{\mathbf {T}}_{\epsilon }\)  is identical to \(\bar{\mathbf {T}}_{G}\)  for higher rate of known labels.
This differs from the results on \(\mathcal {G}^{\prime }\)  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.
<FIGURE>
Subproblem 4
Here we consider the same problem (REF ) as in the previous section but now compare the performance using reference operators \(\bar{\mathbf {P}}_{\alpha }\in \bar{\mathcal {P}}\)  and their respective equivalents \(\bar{\mathbf {T}}_{\alpha }=\Pi ^{1/2}\bar{\mathbf {P}}_{\alpha }\Pi ^{-1/2}\) . The random walk operator is built from the Google matrix approach (Approach 1). The random variables \(y_j, j=1\ldots ,N\)  are distributed as \(y_j=\varepsilon _jf_{0}(v_j), \varepsilon _j\sim \operatorname{Ber}(p)\) . We summarize the case in the subproblemfour.
<TABLE><FIGURE>
Numerical simulations: Case 1
We evaluate the reconstruction performance of \(f_0^{\prime }\)  using filters either based on \(\bar{\mathbf {P}}_{\alpha }, \alpha \in [0,1]\)  at the top of the fig:recsignalcaseonesubpfour or either based on \(\bar{\mathbf {T}}_{\alpha }, \alpha \in [0,1]\)  at the bottom of the fig:recsignalcaseonesubpfour. For all proportions \(p\)  of known labels, the reconstruction performances using a filter based on \(\bar{\mathbf {P}}_{\alpha },\alpha \simeq 0.6\)  are the best. We also notice that the reconstruction performance using a filter with \(\alpha \in \lbrace 0.7,0.8\rbrace \)  are better than \(\alpha \simeq 0.5\) . This suggests that, for this signal, graph filters based on \(\bar{\mathbf {P}}_{\alpha }\)  perform better when the convex combination of \(\mathbf {P}\)  and \(\mathbf {P}^*\)  involves more the random walk \(\mathbf {P}^*\) , i.e. \(\alpha \in [ 0.6,0.8]\) , than the random walk \(\mathbf {P}\) , i.e. \(\alpha \in [0.1,0.5]\) . At the bottom of fig:recsignalcaseonesubpfour, the reconstruction performances using a filter based on \(\bar{\mathbf {T}}_{\alpha }, \alpha \simeq 0.5\)  are the best, for all proportions \(p\)  of known labels. We notice poorer performance overall when \(\alpha \)  increases from \(\alpha =0.5\)  to \(\alpha =1\) . For all proportions \(p\)  of known labels, the reconstruction performance are better overall if the graph signal \(f_0\)  is in \(\ell ^2(\mathcal {V})\)  than \(f_0\)  is in \(\ell ^2(\mathcal {V},\pi )\) . In other words, the reconstruction performance using filters based on \(\bar{\mathbf {T}}_{\alpha }\)  is better than using filters based on \(\bar{\mathbf {P}}_{\alpha }\) .
To summarize :

The Hilbert space where we consider the graph signal or the distribution of the sampled graph signal has an influence on the performance rates of the reconstruction of the graph signal. In the example discussed in the current section, we obtain better performances when we consider the Hilbert space of the graph signals associated to the measure corresponding to the sampling probability.

The construction of graph filters based on the random walk operator yields better performance rates than filters based on the adjacency matrix.

The choice of the parameter \(\alpha \)  in the convex combination between \(\mathbf {P}\)  and \(\mathbf {P}^{*}\)  has an influence on the performance rates. In the example discussed in the current section, the optimal convex combination \(\mathbf {P}\)  and \(\mathbf {P}^{*}\)  (resp. \(\mathbf {T}\)  and \(\mathbf {T}^{*}\) ) is rather forward than backward with one exception.


Multi-resolution analyses on directed graphs
In the previous sections, we proposed a Fourier like basis on directed graphs, as the set of random walk's eigenvectors and determine a frequency analysis by studying the smoothness of the random walk's eigenvectors. We are thus able to construct multiresolution analyses of functions over directed graphs. In a 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 construction [8]} and spectral graph wavelets [9]}. In a second instance, we propose a critically sampled wavelet construction on directed graphs generalizing the diffusion wavelets framework [10]}. In our theoretical constructions, we propose to work in the \(\ell ^2(\mathcal {V})\)  space. We use the linear transformation \(\varphi \)  defined in isometricop to go into the \(\ell ^2(\mathcal {V},\pi )\)  space.
Redundant wavelet transform on directed graphs
In this section, we propose a redundant wavelet transform on directed graphs. We follow the construction of spectral wavelets on undirected graphs [9]} and polynomial diffusion frames [8]}. The novelty is the construction of filters designed in the frequency domain via a frequency response function as a linear combination of projectors onto the associated mono-frequency random walk subspaces. We introduce the necessary elements for the construction of redundant wavelets on directed graphs.
Theoretical framework
Let \(\mathcal {G}=(\mathcal {V},\mathcal {E})\)  be a strongly connected directed graph with cardinality \(|\mathcal {V}|=N\) . The directed graph \(\mathcal {G}\)  is characterized by its adjacency matrix \(\mathbf {W}\in \mathbb {R}^{N\times N}_{+}\) .
On \(\mathcal {G}\) , we defined a random walk \(\mathcal {X}\)  characterized by its transition matrix \(\mathbf {P}=\mathbf {D}^{-1}\mathbf {W}\in \mathbb {R}^{N\times N}_{+}\) . The \(\mathcal {X}\)  is ergodic with stationary distribution \(\pi \) . We introduce \(\mathbf {T}=\Pi ^{1/2}\mathbf {P}\Pi ^{-1/2}\in \ell ^2(\mathcal {V})\) , the operator similar to \(\mathbf {P}\) . We assume \(\mathbf {T}\)  diagonalizable. We define a low-pass operator as follows.
Definition 8.1 A low pass operator at dilation \(t\) , \(\mathbf {H}_{t}^{\mathcal {G}}\)  on a strongly connected directed graph \(\mathcal {G}\)  is defined by
\(\mathbf {H}_{t}^{\mathcal {G}}=\sum _{\ell =1}^p h(t\omega _{\ell })\mathbf {S}_{\ell }, \quad t\in \mathbb {N}.\) 
with \(h:[0,2]\rightarrow \mathbb {R}\) , a function giving a low-pass frequency response and \(\lbrace \mathbf {S}_{\ell }\rbrace _{\ell =1}^p\) , the set of mono-frequency random walk projectors associated to the mono-frequency random walk subspaces \(\lbrace \mathfrak {S}_{\ell }\rbrace _{\ell =1}^p\)  of \(\mathbf {T}\) . The scaling function at dilation \(t\)  and translation \(y\)  is denoted by
\(\mathfrak {h}_{t,y}^{\mathcal {G}}= \mathbf {H}_{t}^{\mathcal {G}}\delta _{y}\) 
where \(\delta _{y}\)  is the Kronecker delta function at the vertex \(y\in \mathcal {V}\) .

We define a bandpass operator as follows.
Definition 8.2 A band pass operator at dilation \(t\) , \(\mathbf {G}_{t}^{\mathcal {G}}\)  on a strongly connected directed graph \(\mathcal {G}\)  is defined as
\(\mathbf {G}_{t}^{\mathcal {G}}=\sum _{\ell =1}^p g(t\omega _{\ell })\mathbf {S}_{\ell }, \quad t\in \mathbb {N}.\) 
with \(g:[0,2]\mapsto \mathbb {R}\) , a function giving a band-pass frequency response and \(\lbrace \mathbf {S}_{\ell }\rbrace _{\ell =1}^p\) , the set of mono-frequency random walk projectors associated to the mono-frequency random walk subspaces \(\lbrace \mathfrak {S}_{\ell }\rbrace _{\ell =1}^p\)  of \(\mathbf {T}\) . The wavelet function at dilation \(t\)  and translation \(y\)  is denoted by
\(\mathfrak {g}_{t,y}^{\mathcal {G}}=\mathbf {G}_{t}^{\mathcal {G}}\delta _{y},\) 
where \(\delta _{y}\)  is the Kronecker delta function at the vertex \(y\in \mathcal {V}\) .

Having defined the low-pass and band-pass operators on directed graphs, we are able to build analysis and synthesis filter banks on directed graphs.
Definition 8.3 We define a bank of synthesis filters \(\mathcal {K}\)  as the set of a low-pass filter at dilation \(t_J\) , \(\mathbf {H}_{t_J}^{\mathcal {G}}\)  and a series of band-pass filters at increasing dilations \(\lbrace t_j\rbrace _{j=1}^J\) , \(\lbrace \mathbf {G}_{t_j}\rbrace _{j=1}^J\)  :
\(\mathcal {K}=\lbrace \mathbf {H}_{t_J}^{\mathcal {G}},\mathbf {G}_{t_1}^{\mathcal {G}},\ldots ,\mathbf {G}_{t_J}^{\mathcal {G}}\rbrace .\) 

Definition 8.4 We define a bank of analysis filters \(\tilde{\mathcal {K}}\)  as the set of a
filter at dilation \(t_J\) , \(\tilde{\mathbf {H}}_{t_J}^{\mathcal {G}}\)  and a series of filters at increasing dilations \(\lbrace t_j\rbrace _{j=1}^J\) , \(\lbrace \tilde{\mathbf {G}_{t_j}}\rbrace _{j=1}^J\)  :
\(\tilde{\mathcal {K}}=\lbrace \tilde{\mathbf {H}}_{t_J}^{\mathcal {G}},\tilde{\mathbf {G}}_{t_1}^{\mathcal {G}},\ldots ,\tilde{\mathbf {G}}_{t_J}^{\mathcal {G}}\rbrace ,\) 
where
\(\tilde{\mathbf {H}}_{t_J}^{\mathcal {G}}=\sum _{\ell =1}^p \tilde{h}(t\omega _{\ell })\mathbf {S}_{\ell }\in \mathbb {R}^{N\times N},\) 
and
\(\tilde{\mathbf {G}}_{t_j}^{\mathcal {G}}=\sum _{\ell =1}^p \tilde{g}(t\omega _{\ell })\mathbf {S}_{\ell }\in \mathbb {R}^{N\times N}, \quad \forall j=1,\ldots ,J.\) 
We also define \(\tilde{\mathfrak {h}}_{t,k}\)  and \(\tilde{\mathfrak {g}}_{t,k}\)  as row vectors as follows:
\(\tilde{\mathfrak {h}}_{t,k}=\delta _k^{\top }\tilde{\mathbf {H}}_{t}^{\mathcal {G}}, \quad \tilde{\mathfrak {g}}_{t,k}=\delta _k^{\top }\tilde{\mathbf {G}}_{t}^{\mathcal {G}}.\) 

Proposition 8.1 
Given a fixed set of increasing dilations \(\lbrace t_j\rbrace _{j=1}^J\) , the synthesis and analysis filters sets respectively \(\mathcal {K}\)  and \(\tilde{\mathcal {K}}\) , the perfect reconstruction condition
\(\mathbf {H}_{t_J}^{\mathcal {G}}\tilde{\mathbf {H}}_{t_J}^{\mathcal {G}}+\sum _{j=1}^J\mathbf {G}_{t_j}^{\mathcal {G}}\tilde{\mathbf {G}}_{t_j}^{\mathcal {G}}=\mathbf {I}\) 
is guaranteed if and only if
\({h}(t_J\omega _{\ell }) \tilde{h}(t_J\omega _\ell )+\sum _{j=1}^J g(t_j\omega _\ell )\tilde{g}(t_j\omega _\ell )=1, \quad \forall \omega _\ell \in [0,2].\) 

\(\mathbf {H}_{t_J}^{\mathcal {G}}\tilde{\mathbf {H}}_{t_J}^{\mathcal {G}}+\sum _{j=1}^J\mathbf {G}_{t_j}^{\mathcal {G}}\tilde{\mathbf {G}}_{t_j}^{\mathcal {G}}&=\sum _{\ell ,\ell ^{\prime }=1}^p h(t_J\omega _{\ell })\tilde{h}(t_J\omega _{\ell ^{\prime }})\mathbf {S}_{\ell } \mathbf {S}_{\ell ^{\prime }}+\sum _{j=1}^{J} \sum _{\ell ,\ell ^{\prime }=1}^p g(t_j\omega _{\ell })\tilde{g}(t_j\omega _{\ell ^{\prime }})\mathbf {S}_{\ell } \mathbf {S}_{\ell ^{\prime }}\\&=\sum _{\ell =1}^p h(t_J\omega _{\ell })\tilde{h}(t_J\omega _{\ell })\mathbf {S}_{\ell } +\sum _{j=1}^{J} \sum _{\ell =1}^p g(t_j\omega _{\ell })\tilde{g}(t_j\omega _{\ell })\mathbf {S}_{\ell }\\&=\sum _{\ell =1}^p \big [h(t_J\omega _{\ell })\tilde{h}(t_J\omega _{\ell }) +\sum _{j=1}^{J} g(t_j\omega _{\ell })\tilde{g}(t_j\omega _{\ell })\big ]\mathbf {S}_{\ell }\) 
Consequently, (REF ) is guaranteed if and only if
\(h(t_J\omega _{\ell })\tilde{h}(t_J\omega _{\ell }) +\sum _{j=1}^{J} g(t_j\omega _{\ell })\tilde{g}(t_j\omega _{\ell })=1 \quad \forall \ell =1\cdots ,p.\) 
We now define a frame as follows.
Definition 8.5 Let \(\mathcal {G}=(\mathcal {V},\mathcal {E})\)  be a directed graph and \(\mu \)  a measure on \(\mathcal {V}\) . A countable family of elements \(\lbrace f_k\rbrace _{k=1}^n \in \ell ^2(\mathcal {V})\)  is said to be a frame if for any graph signal \(f \in \ell ^2(\mathcal {V})\)  we have
\(A\Vert f\Vert ^2\le \sum _{k=1}^n|\langle f,f_k\rangle |^2\le B\Vert f\Vert ^2,\) 
for some constants \(0<A\le B<\infty \)  which are called frame bounds.

We introduce the rectangular matrix \(\mathbf {K}=(\mathbf {H}_{t_J}^{\mathcal {G}},\mathbf {G}_{t_1}^{\mathcal {G}},\ldots ,\mathbf {G}_{t_J}^{\mathcal {G}})\in \mathbb {R}^{N\times N(J+1)}\)  where :
\(\mathbf {K}=\begin{pmatrix}\mathfrak {h}_{t_J,1}^{\mathcal {G}},\ldots ,\mathfrak {h}_{t_J,N}^{\mathcal {G}},\mathfrak {g}_{t_1,1}^{\mathcal {G}},\ldots ,\mathfrak {g}_{t_J,1}^{\mathcal {G}},\ldots ,\mathfrak {g}_{t_J,N}^{\mathcal {G}}\end{pmatrix}\in \mathbb {R}^{N\times N(J+1)}.\) 
We also introduce the following rectangular matrix \(\tilde{\mathbf {K}}=(\tilde{\mathbf {H}}_{t_J}^{\mathcal {G}},\tilde{\mathbf {G}}_{t_1}^{\mathcal {G}},\ldots ,\tilde{\mathbf {G}}_{t_J}^{\mathcal {G}})\in \mathbb {R}^{N(J+1)\times N}\)
\(\tilde{\mathbf {K}}=\begin{pmatrix}\tilde{\mathfrak {h}}_{t_J,1}^{\mathcal {G}}\\\vdots \\\tilde{\mathfrak {h}}_{t_J,N}^{\mathcal {G}}\\\tilde{\mathfrak {g}}_{t_1,1}^{\mathcal {G}}\\\vdots \\\tilde{\mathfrak {g}}_{t_J,1}^{\mathcal {G}}\\\vdots \\\tilde{\mathfrak {g}}_{t_J,N}^{\mathcal {G}}\end{pmatrix}\in \mathbb {R}^{N(J+1)\times N}.\) 
Proposition 8.2 Let \(\mathbf {K}\)  and \(\tilde{\mathbf {K}}\)  be respectively the synthesis and analysis filter banks. We assume that the filter banks \(\mathbf {K}\)  and \(\tilde{\mathbf {K}}\)  verifies the perfect reconstruction condition enunciated at perfectreccond. Consequently, \(\tilde{\mathbf {K}}\)  is a frame for \(\ell ^2(\mathcal {V})\)  with lower frame bound \(1/\Vert \mathbf {K}\Vert ^2\)  and upper frame bound \(\Vert \tilde{\mathbf {K}}\Vert ^2\) .

We have
\(\Vert \tilde{\mathbf {K}}f\Vert ^2=\sum _{k=1}^N|\langle \tilde{\mathfrak {h}}_{t_J,k}^{\mathcal {G}},f\rangle |^2+\sum _{j=1}^J\sum _{k=1}^N \langle \tilde{\mathfrak {g}}_{t_j,k}^{\mathcal {G}},f\rangle |^2.\) 
Consequently, we can write
\(& A\Vert f\Vert ^2\le \Vert \tilde{\mathbf {K}}f\Vert ^2\le B\Vert f\Vert ^2,\\&A\Vert f\Vert ^2\le \sum _{k=1}^N|\langle \tilde{\mathfrak {h}}_{t_J,k}^{\mathcal {G}},f\rangle |^2+\sum _{j=1}^J\sum _{k=1}^N \langle \tilde{\mathfrak {g}}_{t_j,k}^{\mathcal {G}},f\rangle |^2\le B\Vert f\Vert ^2,\) 

where \(A=1/\Vert \mathbf {K}\Vert ^2\)  et \(B=\Vert \tilde{\mathbf {K}}\Vert ^2\) .
The right side is obtained by Cauchy-Schwarz inequality applied on \(\Vert \tilde{\mathbf {K}}f\Vert ^2\) . The left side is obtained by the fact that \(\Vert f\Vert ^2=\Vert \mathbf {K}\tilde{\mathbf {K}}f\Vert ^2\le \Vert \mathbf {K}\Vert ^2\Vert \tilde{\mathbf {K}}f\Vert ^2\)  and application of the Cauchy-Schwarz inequality. As a result, \(\tilde{\mathbf {K}}\)  is a frame.
Definition 8.6 (Wavelet decomposition of a signal on a directed graph) 
Any graph signal \(f\in \ell ^2(\mathcal {V})\)  can be expressed as follows :
\(f=\sum _{k=1}^N \langle f,\tilde{\mathfrak {h}}_{t_J,k}\rangle \,\mathfrak {h}_{t_J,k}+\sum _{j=1}^J\sum _{k=1}^N\langle f,\tilde{\mathfrak {g}}_{t_j,k}\rangle \,\mathfrak {g}_{t_j,k}.\) 


Critically sampled wavelet transform on directed graphs
Notations
We first recall relevant and useful notations for a clear understanding of the construction of diffusion wavelets [10]}, [14]}. 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 [10]}, [14]} except that we adopt a column vector convention whereas a row vector convention was used in [10]}, [14]}.
Let \(\mathsf {V}_0=\ell ^2(\mathcal {V})\)  be the space of functions defined over the vertices of a directed graph \(\mathcal {G}=(\mathcal {V},\mathcal {E})\) . If \(\mathbf {L}\)  is a linear transformation of \(\mathsf {V}_0\)  into \(\mathsf {V}_0\) , \([\mathbf {L}]_{B_1}^{B_2}\)  indicates the matrix representing the linear transformation \(\mathbf {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 \([X]_{B}\)  where the columns of \([X]_{B}\)  are the coordinates of the vectors \(X\)  in the coordinate system defined by \(B\) . Generally, \([B_1]_{B_2}=[\mathbf {I}]_{B_1}^{B_2}\)  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, \([\mathbf {I}]_{B}^{B}=\mathbf {I}\)  the identity matrix. We will also abuse this notation if \(B_2\)  spans a subspace of the space spanned by \(B_1\) .
If \(B_1\)  and \(B_2\)  are two bases, the representations of \(X\)  in \(B_1\)  and \(B_2\)  are related as follows: \([X]_{B_2} = [\mathbf {I}]_{B_1}^{B_2}[X]_{B_1}\) . If \(B_1,B_2,B_3,B_4\)  are arbitrary bases, the matrix representation of \(\mathbf {L}\)  with respect to the basis \(B_1\)  in the domain and \(B_2\)  in the range in the different bases is expressed as
\([\mathbf {L}]_{B_1}^{B_2}=[\mathbf {I}]_{B_4}^{B_2}[\mathbf {L}]_{B_3}^{B_4}[\mathbf {I}]_{B_2}^{B_1}.\) 
Furthermore, if \(B_1\)  and \(B_2\)  are linearly independent sets of vectors that do not span the whole space \(\mathsf {V}_0\)  then we will still use the notation \([\mathbf {L}]_{B_1}^{B_2}\) , but in that case it will represent the restriction of the linear transformation \(\mathbf {L}\)  to the domain and range subspaces spanned by \(B_1\)  and \(B_2\) .

Diffusion wavelets
The construction of the diffusion wavelets [10]}, [14]} enables a multi-resolution analysis of functions on graphs generalizing the concept of classical multiresolution analysis [21]}. The starting point of this construction is a diffusion operator \(\mathbf {T}\) . Similarly to the classical multiresolution analysis, diffusion wavelets is characterized by a family of nested subspaces \(\mathsf {V}_0\supseteq \mathsf {V}_1,\cdots \mathsf {V}_j\supseteq \cdots \) , where each subspace \(\mathsf {V}_j \)  is spanned by a basis of scaling functions \(\Phi _j\) . The complement of \(\mathsf {V}_{j+1}\)  into \(\mathsf {V}_{j}\)  is called \(\mathsf {W}_{j}\)  and is spanned by a set of diffusion wavelets \(\Psi _j\) .

Construction
The construction of diffusion wavelets proceeds with a diffusion operator \(\mathbf {T}\)  defined on a directed strongly connected graph \(\mathcal {G}=(\mathcal {V},\mathcal {E})\) . In the original diffusion wavelets framework by Coifman and Maggioni [10]}, the graph is assumed undirected and they suggest using the reversible random walk operator \(\mathbf {P}\) .
More generally, any operator \(\mathbf {T}\)  can be used as long as it is a low-pass filter when it is applied to localized graph signals. Given a directed graph \(\mathcal {G}\) , an appropriate choice is to use the operator \(\mathbf {T}=\Pi ^{1/2}\mathbf {P}\Pi ^{-1/2}\)  where \(\mathbf {P}\)  is the random walk operator with low-pass frequency response with respect to the definition proposed in the section .
Assuming an ordering of vertices in \(\mathcal {V}\) , \(\mathbf {T}\)  is originally represented on the canonical basis \(\Phi _0=\lbrace \delta _k\rbrace _{k\in \mathcal {V}}\)  of \(\mathsf {V}_0\)  where \(\delta _k\)  is the Kronecker delta function corresponding to the vertex \(k\in \mathcal {V}\) . Using the notations from the section REF , \([\mathbf {T}]_{\Phi _0}^{\Phi _0}\)  is the matrix representation of the linear operator \(\mathbf {T}\)  with respect to the basis \(\Phi _0\)  in the domain and \(\Phi _0\)  in the range.
The columns of the matrix \([\mathbf {T}]_{\Phi _0}^{\Phi _0}\)  can be described as a set of functions \(\widetilde{\Phi }_1=\lbrace \widetilde{\phi }_{1,k}\rbrace _{k\in \mathcal {V}}\)  represented in the basis \(\Phi _0\) . Each element \(\widetilde{\phi }_{1,k},k\in \mathcal {V}\)  of \(\widetilde{\Phi }_{1}\)  corresponds to the action of \(\mathbf {T}\)  on the Kronecker \(\delta \) -function at vertex \(k\)  and is represented as \([\widetilde{\phi }_{1,k}]_{\Phi _0}=[\mathbf {T}]_{\Phi _0}^{\Phi _0}\delta _k\) .
The low-pass property of \(\mathbf {T}\)  leads to the fact that each element \(\widetilde{\phi }_{1,k}\)  of \(\widetilde{\Phi }_{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 \(\widetilde{\phi }_{1,k},k\in \mathcal {V}\)  are therefore interpreted as scaling functions. Due to the fact that the support of each element \(\widetilde{\phi }_{1,k}\)  covers a small neighborhood around their respective vertices, these elements \(\widetilde{\phi }_{1,k}\)  can generally be well approximated by a linear combination of the other functions \(\widetilde{\phi }_{1,l}\)  with \(l\ne k\) . The next stage of the construction is a column subset selection stage [23]}. The aim is to find a small number of columns of \(\tilde{\Phi }_1\) , forming a set \(\mathcal {C}\)  such that the residual \(\Vert [\tilde{\Phi }_1]_{\Phi _0}-\Pi _\mathcal {C}[\tilde{\Phi }_1]_{\Phi _0}\Vert _{\beta }\)  is as minimal as possible where \(\Pi _{\mathcal {C}}\)  is the projection matrix onto the space spanned by the columns of \(\mathcal {C}\)  and \(\beta \)  denotes the spectral norm or the Frobenius norm.
More precisely, this step involves the selection of a subset \(\lbrace \widetilde{\phi }_{1,k},k\in \mathcal {I}_1\rbrace \)  from \(\widetilde{\Phi }_{1}\) , \(|\mathcal {I}_1|\le |\mathcal {V}|\) . We select a subset \(\lbrace \widetilde{\phi }_{1,k},k\in \mathcal {I}_1\rbrace \)  such that all \(\widetilde{\phi }_{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 \(\mathcal {I}_1\)  iteratively by using a modified Gram-Schmidt orthogonalization procedure [10]}.
The subset \(\lbrace \widetilde{\phi }_{1,k},k\in \mathcal {I}_1\rbrace \)  spans a subspace \(\mathsf {V}_1\)  which corresponds to the first approximation subspace of the multi-resolution analysis. We will thus denote \(\lbrace \widetilde{\phi }_{1,k},k\in \mathcal {I}_1\rbrace \)  as \(\Phi _1\) , where \(\Phi _1\)  is by definition a basis of \(\mathsf {V}_1\)  that is generally not orthogonal. The vectors in \(\Phi _1\)  are the scaling functions at scale 1.
In order to define the next scales of scaling functions, we consider the compression step of the diffusion operator \(\mathbf {T}\)  on the subspace \(\mathsf {V}_1\) , that is its restriction on \(\mathsf {V}_1\) . The latter can be represented in \(\Phi _1\)  as
\([\mathbf {T}]_{\Phi _1}^{\Phi _1}=[\mathbf {I}]_{\Phi _0}^{\Phi _1}[\mathbf {T}]_{\Phi _0}^{\Phi _0} [\mathbf {I}]_{\Phi _1}^{\Phi _0}.\) 
where \([\mathbf {I}]_{\Phi _0}^{\Phi _1}\)  represents the restriction of the signal space \(\mathsf {V}_0\)  to \(\mathsf {V}_1\)  (with respective bases \(\Phi _0\)  and \(\Phi _1\) ) and \([\mathbf {I}]_{\Phi _1}^{\Phi _0}\)  represents the embedding of \(\mathsf {V}_1\)  in \(\mathsf {V}_0\) .
The next approximation space \(\mathsf {V}_2\subset \mathsf {V}_1\)  and its associated basis \(\Phi _2\)  can be obtained in the same way as the definition of \(\mathsf {V}_1\)  and its basis \(\Phi _1\)  except that we now consider the operator \(\mathbf {T}^2\)  restricted to \(\mathsf {V}_1\)  instead of \(\mathbf {T}\)  in \(\mathsf {V}_0\) .
The columns of \([\mathbf {T}^2]_{\Phi _1}^{\Phi _1} = ([\mathbf {T}]_{\Phi _1}^{\Phi _1})^2\)  can be interpreted as scaling functions at scale 2, \(\widetilde{\Phi }_2=\lbrace \widetilde{\phi }_{2,k}\rbrace \) , represented in the basis \(\Phi _1\) . From these functions we extract a subset \(\Phi _2 = \lbrace \widetilde{\phi }_{2,k},k\in \mathcal {I}_2\rbrace \)  such that any function in \(\widetilde{\Phi }_2\)  is well-approximated by a linear combination of functions of \(\Phi _2\) .
After \(j\)  iterations of this procedure we have defined \(j\)  approximation subspaces \(\mathsf {V}_j\subset \mathsf {V}_{j-1}\subset \cdots \subset \mathsf {V}_1\)  with corresponding bases \(\Phi _j, \Phi _{j-1},\ldots , \Phi _1\) . At each step the basis \(\Phi _j\)  is defined by its representation in the basis \(\Phi _{j-1}\)  based on the restriction of the operator \(\mathbf {T}^{2^j}\)  to \(\mathsf {V}_{j-1}\) . In order to represent these functions in the original basis \(\Phi _0\)  of \(\mathsf {V}_0\)  we can use
\([\Phi _j]_{\Phi _0} = [\mathbf {I}]_{\Phi _j}^{\Phi _0}=[\mathbf {I}]_{\Phi _{0}}^{\Phi _0} [\mathbf {I}]_{\Phi _{1}}^{\Phi _0}\cdots [\mathbf {I}]_{\Phi _{j-1}}^{\Phi _{j-2}}[\mathbf {I}]_{\Phi _{j}}^{\Phi _{j-1}}.\) 
Since every function in \(\Phi _0\)  is defined on \(\mathsf {V}_0\) , so is every function on \(\Phi _j\) . Hence any function in the approximation space \(\mathsf {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 \(\Psi _j\)  for the subspaces \(\mathsf {W}_j\)  by selecting a subset of the columns of the band pass operator \([\mathbf {I}_{\mathsf {V}_j}-\Phi _{j+1}\Phi _{j+1}^{\dagger }]_{\Phi _j}^{\Phi _j}\)  which is the orthogonal projector on the complement of \(\mathsf {V}_{j+1}\)  into \(\mathsf {V}_j\) . The wavelets capture the detail lost from going from \(\mathsf {V}_j\)  to \(\mathsf {V}_{j+1}\) . As our framework falls into bi-orthogonal scope, we need to build the dual wavelet bases. For a scale \(j\) , we have a wavelet base \(\Psi _j\)  and we need to construct the associated dual base \(\widehat{\Psi }_j\)  obtained such that we have the following relation
\(\widehat{\Psi }_j^{\dagger }\Psi _j=\mathbf {I}.\) 
where \(\widehat{\Psi }_j^{\dagger }\)  is the pseudo-inverse of \(\Psi _j\) . We mention that in [14]}, biorthogonal 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 REF . The idea is merely to replace the powers of \(\mathbf {T}\)  as multiresolution scaling operators by more arbitrary low-pass filters. More precisely, where the operator \(\mathbf {T}^{2^j}\)  is used to define the approximation space \(\mathsf {V}_{j+1}\)  in section REF , we propose to use instead a low-pass graph filter \(\mathbf {H}_j\) . If the graph filters \(\mathbf {H}_j\)  correspond to the scaling operators defined in spectralwaveletdg 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 REF .

Applications
This section is devoted to illustrate our multiresolution analysis framework through the visualization of wavelets and scaling functions on regular-type structures and semi-supervised learning on directed graphs. The application framework will be the same as in section .
Multiresolution analysis on the directed cycle graph
We show an example of multi-resolution analysis on the directed cycle graph \(\mathcal {C}_N\)  with \(N=256\) . We use the same assumptions made in section REF . We construct both orthogonal and biorthogonal multi-resolution analyses on \(\mathcal {C}_N\)  through the framework of the diffusion wavelets applied on the dyadic powers of \(\mathbf {T}\) , i.e. \(\lbrace \mathbf {T}^{2^j}\rbrace _{j=1}^J\) . We set the number of scales at \(J=6\) .
<FIGURE><FIGURE>fig:ortbiortscf shows the orthogonal and biorthogonal scaling functions at scales 1,3 and 5. At each scale, we represent 3 or 4 scaling functions.
fig:ortbiortwavf 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 at scale 3 and 5 that biorthogonal wavelet functions have support slightly smaller than the orthogonal wavelets functions.
Furthermore, we note that orthogonal scaling functions as well orthogonal wavelet functions have more oscillations than their biorthogonal counterparts. That means that orthogonal scaling and wavelet functions have a poorer frequency localization. The conclusions about the spatial localization of the scaling functions are identical to the conclusions in [26]}. 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 propagate towards the parents' node of the selected vertex. We can also note the orthogonal wavelet transform is more robust with respect to perturbations than the biorthogonal wavelet transform. Let \(\mathcal {W}_{ortho}=[\Psi _{1},\ldots ,\Psi _{5},\Phi _{5}]\)  be the inverse wavelet transform and \(\mathcal {W}_{biorth}\)  its biorthogonal counterpart. We have that the condition numbers \(\kappa (\mathcal {W}_{ortho})\approx 33\)  and \(\kappa (\mathcal {W}_{biortho})\approx 2\times 10^{4}\) . Consequently, the biorthogonal wavelet transform is much more sensitive to perturbations than its orthogonal counterpart.

Multiresolution analysis 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. The Watts-Strogatz model [27]} is a undirected random graph model exhibiting small world properties including short average path lengths and high clustering [28]}. The construction of a directed graph 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 \(\mathcal {G}\sim \operatorname{DWS}(N,K,\beta )\)  a graph constructed following a directed Watts-Strogatz model with \(N\)  vertices, \(K\)  nearest neighbors and rewiring probability \(\beta \) .
In a first instance, we analyze the scaling functions built from low-pass filters based on \(\bar{\mathbf {T}}_{\alpha }\)  at a given scale of \(\mathcal {G}_1\sim \operatorname{DWS}(64,2,0)\) . 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 \(\bar{\mathbf {T}}_{\alpha }\)  at a given scale of \(\mathcal {G}_2\sim \operatorname{DWS}(64,2,0.02)\) . Lastly, we show both orthogonal and biorthogonal scaling and wavelets functions built following the construction Section REF . We summarize the following cases in table REF .
<TABLE>Case 1
In this case, we analyze some scaling functions built from low-pass filters based on \(\bar{\mathbf {P}}_{\alpha }\)  at scale \(t=2^4\)  of a directed Watts-Strogatz graph \(\mathcal {G}_2\sim \operatorname{DWS}(64,2,0)\) . We consider the following low-pass filters
\(\bar{\mathbf {T}}_{\alpha }=\Pi ^{1/2}\bar{\mathbf {P}}_{\alpha }\Pi ^{-1/2},\quad \forall \alpha \in [0,1].\) 
and we construct a collection of low-pass graph filters at a given scale following the construction presented in spectralwaveletdg. More precisely, we build the filters \(\mathbf {H}_{\alpha }\)  as follows
\(\mathbf {H}_{\alpha }=\sum _{\omega \in \omega } h(t\omega )\mathbf {S}_{\omega ,\alpha },\) 
with \(t=2^{4}\)  and \(h(x)=\exp (-x)\) .
<FIGURE>mracaseone represents the 50th scaling function \(\mathfrak {h}_{\alpha ,50}=\mathbf {H}_{\alpha }\delta _{50}\)  at scale \(t=2^4\)  built from low-pass filters \(\mathbf {T}_{\alpha }\) , for all \(\alpha \in [0,1]\)  according the construction spectralwaveletdg on the directed Watts-Strogatz graph \(\mathcal {G}\sim \operatorname{DWS}(64,2,0)\) . 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 \(\mathbf {T}_{\alpha }\)  admit the same discrete Fourier basis and the associated scaling functions \(\mathfrak {h}_{\alpha ,50}\)  are exactly the same. We observe that the scaling function is centered around the node 50 and has a symmetric shape.
In this case, for non symmetric graph filters \(\mathbf {H}_{\alpha }\) , the scaling functions have symmetric shape and are centered around its given node.

Case 2
In this case, we analyze some scaling functions built from low-pass filters based on \(\bar{\mathbf {P}}_{\alpha }\)  at scale 10 of \(\mathcal {G}_2\sim \operatorname{DWS}(64,2,0.02)\) . We also consider the following low-pass filters
\(\bar{\mathbf {T}}_{\alpha }=\Pi ^{1/2}\bar{\mathbf {P}}_{\alpha }\Pi ^{-1/2},\quad \forall \alpha \in \lbrace 0,0.5,1\rbrace .\) 
and we construct a collection of low-pass graph filters at a given scale following the construction presented in spectralwaveletdg. More precisely, we build the filters \(\mathbf {H}_{\alpha }\)  as follows
\(\mathbf {H}_{\alpha }=\sum _{\omega \in \omega } h(t\omega )\mathbf {S}_{\omega ,\alpha }\) 
with \(t=2^{4}\)  and \(h(x)=\exp (-x)\) . We observe the 50th scaling functions at the scale \(t=2^{4}\) , that is \(\mathfrak {h}_{\alpha ,50}= \mathbf {H}_{\alpha }\delta _{50}\) .
<FIGURE>mracasetwopartone represent the 50th scaling function at the scale \(t=2^4\) , for different \(\alpha \) , \(\mathfrak {h}_{\alpha ,50}\) . The graph \(\mathcal {G}_2\)  admits a directed edge from the node 50 to the node 11 and a directed edge from the node 30 to the node 50. As \(\alpha \)  increases towards 1, we observe that the scaling function \(\mathfrak {h}_{\alpha ,50}\)  propagates around the child node 11. This means the more \(\alpha \)  increases, the more the influence of the backward random walk \(\mathbf {P}^{*}\)  is important. Wee note for \(\alpha =0.5\) , the scaling function \(\mathfrak {h}_{\alpha ,50}\)  is centered around the child node 11 an the parent node 30. This case seems to be the equilibrium.
Finally, as \(\alpha \)  goes to 0, we observe that \(\mathfrak {h}_{\alpha ,50}\)  propagates around the parent node 30.

Case 3
Here, we consider orthogonal and biorthogonal scaling functions on a graph \(\mathcal {G}\sim \operatorname{DWS}(64,2,0.02)\)  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 \(\mathbf {H}=\mathbf {T}^{2}\)  and we look at the dyadic powers of \(\mathbf {H}\)  i.e. \(\lbrace \mathbf {H}^{2^j}\rbrace _{j=1}^J\) . We set the number of scales at \(J=5\) . We observe the orthogonal an biorthogonal scaling functions at node 49 at scale 3 obtained by the diffusion wavelets procedure and we compare to the scaling function at node 49 built from a graph filter based on \(\bar{\mathbf {T}}\)  at the scale \(t=16\) .
<FIGURE>fig:scwscasethree shows the orthogonal and bi-orthogonal scaling functions at node 49 built w.r.t. the diffusion wavelets framework and the scaling function at node 49 built w.r.t. 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 a large size support. By contrast, the orthogonal scaling function has a poor spatial localization around the nodes 11 and 30.

Semi-supervised learning on directed graphs with \(\ell _1\) -regularization
We discuss the semi-supervised learning approach in the case of functions over directed graphs with \(\ell _1\)  regularization on the wavelet coefficients. Our aim is to show that the performance of the semi-supervised learning problem with \(\ell _1\)  regularization is competitive to the existing approaches, e.g. the semi-supervised learning problem studied in section . We discuss the following semi-supervised learning problem, again on the example of the political blogs dataset [3]}. We use the same notations described in sect:applicationsone and work on the subgraph \(\mathcal {G}^{\prime }\) . The \(\ell _1\) -regularized semi-supervised learning problem is
\(w^{*}=\underset{w}{\operatorname{argmin}}\:\Vert \tilde{y}-\mathbf {M}\mathbf {K}w\Vert _2^2+\lambda \Vert w\Vert _1,\quad \lambda \in \mathbb {R}^{+}.\) 
The graph signal \(\tilde{y}=\mathbf {M}y\)  is the partially labeled graph signal with \(\mathbf {M}=\lbrace m_{ij}\rbrace _{\lbrace 1\le i,j\le N\rbrace }\)  is the mask operator, i.e. the diagonal matrix where \(m_{ii}=\mathbb {1}_{v_i\in \mathcal {O}}\)  where \(\mathcal {O}\subset \mathcal {V}\)  is the subset of known labels.
The matrix \(\mathbf {K}=(\mathbf {H}_J,\mathbf {G}_1,\cdots ,\mathbf {G}_J)\)  is the synthesis filter bank. If we set \(\mathbf {X}=\mathbf {MK}\) , the equation (REF ) can be rewritten as
\(w^{*}=\underset{w}{\operatorname{argmin}}\:\Vert \tilde{y}-\mathbf {X}w\Vert _2^2+\lambda \Vert w\Vert _1,\quad \lambda \in \mathbb {R}.\) 
The formulation (REF ) is identical to the formulation of the problem of signal restoration with redundant wavelet transforms in [30]}, except that \(\mathbf {K}\)  is the synthesis wavelet transform for functions defined over directed graphs. Furthermore, previous approaches of semi-supervised learning on undirected graphs have been investigated using overcomplete graph wavelets [31]} or critically sampled spline graph wavelets [32]}. The \(\ell _1\) -regularized synthesis semi-supervised learning problem is convex and can be solved efficiently using e.g. proximal splitting methods [33]}, [34]}, [35]}. From the solution \(w^{*}\)  of semisuplone, we define the restored signal \(f^{*}\)  as
\(f^{*}=\operatorname{sign}(\mathbf {K}w^{*}).\) 
where \(\mathcal {L}\)  is the directed normalized Laplacian defined at (REF ) and the graph bandpass filters \(\mathbf {G}_1,\mathbf {G}_2\)  are
We compare the performances of formulations (REF ) and (REF ). The formulation (REF ) requires a synthesis filter bank \(\mathbf {K}\) . We choose a two scales synthesis graph filter
\(\mathbf {K}\) . We choose a two scales synthesis graph filter bank \(\mathbf {K}=(\mathbf {H}_2,\mathbf {G}_1,\mathbf {G}_2)\)  based on a heat kernel construction [36]}
\(\mathbf {H}_j=e^{-\mathcal {L}t_j},\quad t_j=2^j, \quad j=1,2.\) 
with \(\mathcal {L}\)  is the directed Laplacian defined at (REF )
and the graph bandpass filters \(\mathbf {G}_1,\mathbf {G}_2\)  as
\(\mathbf {G}_2=\mathbf {H}_{1}-\mathbf {H}_2,\quad \mathbf {G}_1=\mathbf {I}-\mathbf {H}_1.\) 
<FIGURE>sslwavvsdir shows the performance of semi-supervised learning between the formulation (REF ) and (REF ) on the largest strongly connected graph obtained from US political blogs.
The performance is obtained by averaging 200 realizations and determining the value of the parameters \(\lambda \)  and \(\gamma \)  associated respectively to the approaches (REF ) and (REF ) giving the best performance rates. On sslwavvsdir, the performance based on the \(\ell _1\)  penalization of redundant graph wavelets coefficients on directed graphs is competitive against the problem with \(\ell _2\)  Dirichlet regularization term for all percentage of known labels. Normally, one should expect better performances in the \(\ell _1\)  case but it seems not the case because our graph wavelets are not good in terms of spatial localization. Still, in order to improve the performance of the method (REF ), it could be wise to build new wavelets with a better space-frequency characterization, and choose the best suitable number of scales with respect to the characteristics of the data set.

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 to existing approaches.

Acknowledgements
The authors would like to thank Romain Cosentino, Randall Balestriero, Yanis Bahroun, Benjamin Girault and Antonio Ortega for the useful discussions.
