Data-driven Thresholding in Denoising with Spectral Graph Wavelet Transform

TL;DR

This paper introduces a data-driven thresholding method for denoising signals in Graph Signal Processing using Spectral Graph Wavelet Transform, optimizing parameters via SURE to improve performance over existing methods.

Contribution

It proposes a novel adaptive thresholding approach in SGWT-based denoising, accounting for correlated noise and providing a multi-scale analysis with improved results.

Findings

Better denoising performance than recent methods

Effective thresholding strategies for over-complete SGWT

Numerical validation on various graph signals

Abstract

This paper is devoted to adaptive signal denoising in the context of Graph Signal Processing (GSP) using Spectral Graph Wavelet Transform (SGWT). This issue is addressed \emph{via} a data-driven thresholding process in the transformed domain by optimizing the parameters in the sense of the Mean Square Error (MSE) using the Stein's Unbiased Risk Estimator (SURE). The SGWT considered is built upon a partition of unity making the transform semi-orthogonal so that the optimization can be performed in the transformed domain. However, since the SGWT is over-complete, the divergence term in the SURE needs to be computed in the context of correlated noise. Two thresholding strategies called coordinatewise and block thresholding process are investigated. For each of them, the SURE is derived for a whole family of elementary thresholding functions among which the soft threshold and the James-Stein threshold. This multi-scales analysis shows better performance than the most recent methods from the literature. That is illustrated numerically for a series of signals on different graphs.

Tables

Table 1. Table 1: Mean SNR performance over N = 10 𝑁 10 N=10 realizations of the low to high noise levels settings with corresponding empirical standard deviation. Left panel: f 0.01 , 2 subscript 𝑓 0.01 2 f_{0.01,2} and right panel: f 0.001 , 4 subscript 𝑓 0.001 4 f_{0.001,4} .

${\mathrm{SNR}}_{\mathrm{in}}$	16.07$\pm$0.13	10.05$\pm$0.13	4.03$\pm$0.13	16.64$\pm$0.13	10.62$\pm$0.13	4.60$\pm$0.13
${\mathrm{MSE}}^{\beta =1,G}$	19.04$\pm$0.24	14.22$\pm$0.26	9.46$\pm$0.26	24.67$\pm$0.33	19.77$\pm$0.36	14.79$\pm$0.45
${\mathrm{MSE}}^{\beta =2,G}$	20.07$\pm$0.24	15.60$\pm$0.30	10.69$\pm$0.30	26.88$\pm$0.29	22.18$\pm$0.37	17.08$\pm$0.67
${\mathrm{SURE}}_{\sigma }^{\beta =1,G}$	18.96$\pm$0.27	14.16$\pm$0.29	9.46$\pm$0.26	24.62$\pm$0.37	19.64$\pm$0.48	14.70$\pm$0.52
${\mathrm{SURE}}_{\sigma }^{\beta =2,G}$	20.04$\pm$0.37	15.49$\pm$0.43	10.64$\pm$0.32	26.73$\pm$0.32	21.91$\pm$0.52	16.88$\pm$0.59
${\mathrm{MSE}}^{\beta =1,\mathrm{LD}}$	19.10$\pm$0.24	14.28$\pm$0.27	9.58$\pm$0.27	24.68$\pm$0.34	19.79$\pm$0.36	14.83$\pm$0.46
${\mathrm{MSE}}^{\beta =2,\mathrm{LD}}$	20.08$\pm$0.24	15.61$\pm$0.29	10.72$\pm$0.30	26.90$\pm$0.26	22.20$\pm$0.36	17.13$\pm$0.67
${\mathrm{SURE}}_{\sigma }^{\beta =1,\mathrm{LD}}$	19.10$\pm$0.24	14.26$\pm$0.26	9.48$\pm$0.24	24.51$\pm$0.40	19.59$\pm$0.46	14.69$\pm$0.49
${\mathrm{SURE}}_{\sigma }^{\beta =2,\mathrm{LD}}$	20.01$\pm$0.31	15.51$\pm$0.36	10.61$\pm$0.39	26.52$\pm$0.32	21.79$\pm$0.34	16.73$\pm$0.60
Wiener	17.01$\pm$0.13	11.87$\pm$0.15	7.43$\pm$0.16	17.91$\pm$0.12	12.86$\pm$0.13	8.40$\pm$0.16
$r_{\inf }$	17.05$\pm$0.00	11.89$\pm$0.00	7.42$\pm$0.00	17.96$\pm$0.00	12.91$\pm$0.00	8.46$\pm$0.00
${\mathrm{MSE}}^{k=2}$	17.35$\pm$0.13	11.43$\pm$ 0.16	5.68$\pm$0.17	19.37$\pm$0.14	13.65$\pm$0.18	8.01$\pm$0.24
${\mathrm{MSE}}^{k=1}$	18.05$\pm$0.14	11.98$\pm$0.18	6.24$\pm$0.18	20.43$\pm$0.15	14.78$\pm$0.18	9.30$\pm$0.27
${\mathrm{MSE}}^{k=0}$	19.57$\pm$0.17	13.43$\pm$0.23	7.62$\pm$0.22	23.38$\pm$0.25	17.88$\pm$0.27	12.80$\pm$0.40

Table 2. Table 2: Mean SNR performance over N = 10 𝑁 10 N=10 realizations with f 0.01 , 5 subscript 𝑓 0.01 5 f_{0.01,5} for the Pittsburgh graph.

$\sigma$	0.004	0.005	0.01
${\hat{\sigma }}_1$	0.0065	0.0072	0.0113
${\hat{\sigma }}_2$	0.0068	0.0074	0.0114
${\mathrm{SNR}}_{\mathrm{in}}$	9.53$\pm$0.29	7.60$\pm$0.29	1.58$\pm$0.29
${\mathrm{MSE}}^{\beta =2,\mathrm{LD}}$	13.65$\pm$0.29	12.22$\pm$0.34	8.46$\pm$0.41
${\mathrm{MSE}}^{k=2}$	12.28$\pm$0.25	11.05$\pm$0.22	8.14$\pm$0.29
${\mathrm{MSE}}^{k=1}$	13.10$\pm$0.21	11.78$\pm$0.22	8.09$\pm$0.36
${\mathrm{MSE}}^{k=0}$	12.83$\pm$0.23	11.42$\pm$0.24	7.56$\pm$0.46
${\mathrm{MSE}}^{\mathrm{Soft}}$	11.22$\pm$0.23	9.62$\pm$0.27	5.22$\pm$0.22
${\mathrm{MSE}}^{\mathrm{Hard}}$	10.22$\pm$0.28	8.50$\pm$0.38	4.03$\pm$0.14
${\mathrm{SURE}}_{\sigma }^{\beta =2,\mathrm{LD}}$	13.33$\pm$0.37	12.00$\pm$0.30	7.95$\pm$0.25
${\mathrm{SURE}}_{{\hat{\sigma }}_1}^{\beta =2,\mathrm{LD}}$	12.35$\pm$0.53	11.38$\pm$0.62	8.13$\pm$0.35
${\mathrm{SURE}}_{{\hat{\sigma }}_2}^{\beta =2,\mathrm{LD}}$	12.12$\pm$0.57	11.27$\pm$0.54	8.07$\pm$0.33

Equations

$$A\Vert f{\Vert }_2^2\le \sum _{i\in I}|⟨f,r_i⟩{|}^2\le B\Vert f{\Vert }_2^2.$$

$$\mathfrak{F}=\left\{\sqrt{{\psi }_j}(\mathcal{L}){\delta }_i,j=0,\dots ,J,i\in V\right\}.$$

$${\psi }_0(x)=\omega (x)\text{and}{\psi }_j(x)=\omega (b^{-j}x)-\omega (b^{-j+1}x)\text{for}j=1,\dots ,J,\text{where}J=\lfloor \frac{\log {\lambda }_1}{\log b}\rfloor +2.$$

$$\mathcal{W}f={\left(\sqrt{{\psi }_0}(\mathcal{L})f^T,\dots ,\sqrt{{\psi }_J}(\mathcal{L})f^T\right)}^T\in {\mathbb{R}}^{n(J+1)}.$$

$${\mathcal{W}}^{\ast }{\left({\eta }_0^T,{\eta }_1^T,\dots ,{\eta }_J^T\right)}^T=\sum _{j\ge 0}\sqrt{{\psi }_j}(\mathcal{L}){\eta }_j.$$

$$\mathbf{E}\Vert h(\tilde{F})-F{\Vert }^2=\mathbf{E}\left[-n{\sigma }^2+\Vert h(\tilde{F})-\tilde{F}{\Vert }^2+2\sum _{i,j=1}^{n(J+1)}\mathbf{Cov}({\Xi }_i,{\Xi }_j){\partial }_j{h}_i(\tilde{F})\right],$$

$$\mathbf{SURE}(h)=-n{\sigma }^2+\Vert h(\tilde{F})-\tilde{F}{\Vert }^2+2\sum _{i,j=1}^{n(J+1)}\mathbf{Cov}({\Xi }_i,{\Xi }_j){\partial }_j{h}_i(\tilde{F}).$$

$$\mathbf{E}\Vert h(\tilde{F})-F{\Vert }^2=\mathbf{E}\Vert h(\tilde{F})-\tilde{F}{\Vert }^2+\mathbf{E}\Vert \Xi {\Vert }^2+2\mathbf{E}⟨h(\tilde{F})-\tilde{F},\Xi ⟩.$$

$$\mathbf{E}⟨g(F+\Xi ),\Xi ⟩=\sum _{i=1}^{n(J+1)}\mathbf{E}[g_i(F+\Xi ){\Xi }_i]=\sum _{i=1}^{n(J+1)}\mathbf{Cov}(g_i(F+\Xi ),{\Xi }_i).$$

$$\mathbf{Cov}(g_i(F+\Xi ),{\Xi }_i)=-n{\sigma }^2+\sum _{j=1}^{n(J+1)}\mathbf{Cov}({\Xi }_i,{\Xi }_j)\mathbf{E}[({\partial }_j{h}_i)(F+\Xi )],$$

$${\partial }_i{h}_i(F+\Xi )=1_{[t_i,\infty )}(|{\tilde{F}}_i|)\left[1+(\beta -1)\frac{t^{\beta }}{|{\tilde{F}}_i{|}^{\beta }}\right].$$

$$\mathbf{SURE}(h)=-n{\sigma }^2+\sum _{i=1}^{n(J+1)}{\tilde{F}}_i^2{\left(1\wedge \frac{t_i^{\beta }}{|{\tilde{F}}_i{|}^{\beta }}\right)}^2+2\sum _{i=1}^{n(J+1)}\mathbf{V}({\Xi }_i)1_{[t_i,\infty )}(|{\tilde{F}}_i|)\left[1+\frac{(\beta -1)t_i^{\beta }}{|{\tilde{F}}_i{|}^{\beta }}\right].$$

$$h_i(x)=x_i\max \left\{1-\frac{t_{\ell }^{\beta }}{\Vert x{\Vert }_{B_{\ell }}^{\beta }},0\right\},x\in {\mathbb{R}}^{n(J+1)},\text{and}\ell \in L:i\in B_{\ell }.$$

$${\partial }_j{h}_i(\tilde{F})={\tilde{F}}_i{1}_{[t_{\ell },\infty )}(\Vert \tilde{F}{\Vert }_{B_{\ell }})\beta t_{\ell }^{\beta }{\tilde{F}}_j\Vert \tilde{F}{\Vert }_{B_{\ell }}^{-\beta -2},$$

$${\partial }_i{h}_i(\tilde{F})=1_{[t_{\ell },\infty )}(\Vert \tilde{F}{\Vert }_{B_{\ell }})\left(1-t_{\ell }^{\beta }\Vert \tilde{F}{\Vert }_{B_{\ell }}^{-\beta }+\beta t_{\ell }^{\beta }{\tilde{F}}_i^2\Vert \tilde{F}{\Vert }_{B_{\ell }}^{-\beta -2}\right).$$

$$\begin{gathered} \mathbf{SURE}(h)=-n{\sigma }^2+\sum _{\ell \in L}{\left(1\wedge \frac{t_{\ell }^{\beta }}{\Vert \tilde{F}{\Vert }_{B_{\ell }}^{\beta }}\right)}^2\Vert \tilde{F}{\Vert }_{B_{\ell }}^2 \\ +2\sum _{\ell \in L}{1}_{[t_{\ell },\infty )}(\Vert \tilde{F}{\Vert }_{B_{\ell }})\left[\left(1-\frac{t_{\ell }^{\beta }}{\Vert \tilde{F}{\Vert }_{B_{\ell }}^{\beta }}\right)\sum _{i\in B_{\ell }}\mathbf{V}({\Xi }_i)+\frac{\beta t_{\ell }^{\beta }}{\Vert \tilde{F}{\Vert }_{B_{\ell }}^{\beta +2}}\sum _{i,j\in B_{\ell }}\mathbf{Cov}({\Xi }_i,{\Xi }_j){\tilde{F}}_i{\tilde{F}}_j\right] \end{gathered}$$

$$\mathbf{E}[\|h(\widetilde{F})-F\|^{2}]=\mathbf{E}\Bigg{[}-\mathrm{Tr}(\mathcal{W}\Gamma\mathcal{W}^{\ast})+\|h(\widetilde{F})-\widetilde{F}\|^{2}+2\sum_{i,j=1}^{n(J+1)}\mathbf{Cov}(\Xi_{i},\Xi_{j})\partial_{j}h_{i}(\widetilde{F})\Bigg{]},$$

$$\Gamma ={\sigma }^2(\mathrm{Id}+2\alpha W+{\alpha }^{2}W{W}^{\ast }).$$

$$\mathbf{E}[{\tilde{f}}^{T}g(\mathcal{L})\tilde{f}]=f^{T}g(\mathcal{L})f+\mathbf{E}[{\xi }^{T}g(\mathcal{L})\xi ]=f^{T}g(\mathcal{L})f+{\sigma }^2\mathrm{Tr}g(\mathcal{L}),$$

$${\hat{\sigma }}_1^2=\frac{{\tilde{f}}^{T}g(\mathcal{L})\tilde{f}}{\mathrm{Tr}g(\mathcal{L})}.$$

$${\hat{\sigma }}_1^2=\frac{{\tilde{f}}^T\mathcal{L}\tilde{f}}{\mathrm{Tr}\mathcal{L}}=\frac{{\sum }_{i,j\in V}{w}_{ij}|\tilde{f}(i)-\tilde{f}(j){|}^2}{2\mathrm{Tr}\mathcal{L}}.$$

$${\hat{\sigma }}_2^2=\frac{{\sum }_{i=nJ+1}^{n(J+1)}(\mathcal{W}\tilde{f}{)}_i^2}{\mathrm{Tr}{\psi }_J(\mathcal{L})}.$$

Full text

Keywords: Spectral Graph Theory . Denoising . Stein Unbiased Risk Estimation . Spectral Graph Wavelet Transform . Tight Frame . Variance Estimation

1 Introduction

The emerging field of Graph Signal Processing (GSP) aims to bridge the gap between signal processing and spectral graph theory (see for instance Chung, (1997); Belkin and Niyogi, (2008) and references therein). One objective is to generalize fundamental analysis operations from regular grid signals to irregular structures as graphs. There is an extensive literature on GSP, in particular we refer the reader to Shuman et al., (2013) for an introduction to this field and Ortega et al., (2018) for an overview of recent developments, challenges and applications. As a matter of fact, GSP have already been applied in machine/deep learning: convolutional neural networks (CNN) on graphs Bruna et al., (2014); Henaff et al., (2015); Defferrard et al., (2016), semi-supervised classification with graph CNN Kipf and Welling, (2017); Hamilton et al., (2017), community detection Tremblay and Borgnat, (2014), to name just a few. In the context of GSP, the authors of Coifman and Lafon, (2006); Gavish et al., (2010); Hammond et al., (2011) have developed wavelet transforms on graphs. More specifically, in Hammond et al., (2011) a fairly general construction of a frame enjoying the usual properties of standard wavelets is developed: each vector of the frame is localized both in the graph domain and the spectral domain. The transform associated with this frame is named Spectral Graph Wavelet Transform (SGWT). Many studies based on SGWT (or some variants) explore the denoising performance of this approach using different strategies Leonardi and Van De Ville, (2013); Onuki et al., (2016); Wang et al., (2016); Deutsch et al., (2016); Irion and Saito, (2017); Dong et al., (2016); Göbel et al., (2018) from signal adapted tight frames to regularization method.

The denoising approach chosen in this paper involves several thresholding processes in the transformed domain of the wavelet coefficients. Actually, this can be seen as an extension to SGWT of the methodology of Donoho and Johnstone, (1995); Cai, (1999). With this approach, the main challenge is the efficient calibration of the parameters minimizing the MSE risk in a complete data-driven way. Recently, in the setting of discrete wavelets transform on a regular grid—the so-called regular case—the Stein’s unbiased risk estimate (SURE) has proven to be a powerful tool for signal/image restoration Luisier et al., (2007); Pesquet et al., (2009); Vaiter et al., (2013). Based on the Stein’s lemma, this estimator acts as a proxy for the MSE which cannot be computed in practice since the original signal is unknown. In this paper, the SURE is explicitly computed for an arbitrary thresholding process in Theorem 1 and for correlated noise in the graph domain in Corollary 1. Also, let us point out that contrary to the regular wavelet transform, the SGWT is no longer orthogonal so that a white Gaussian noise in the graph domain is transformed in a correlated noise. Consequently, the divergence term of the resulting SURE involves the covariance of the transformed noise making the numerical evaluation less simple than in the regular case. Afterward, the SURE is specified to the case of coordinatewise and block thresholding. The latter is inspired by image denoising problems for which a Stein risk estimator has been proposed in Peyré et al., (2011) to tune both the block-sparsity structure and the threshold. A similar selection strategy has been developed by Navarro et al., (2013) in the context of deconvolution. The R package gasper which implements the method introduced in this paper is available on github111https://github.com/fabnavarro/gasper (de Loynes et al.,, 2020) as well as the scripts to reproduce the results presented222https://github.com/fabnavarro/SGWT-SURE.

Many denoising methods consist in shrinking the coefficients of a sparse representation of the considered signals in a certain transformed domain. In Talebi and Milanfar, (2014), the aim is to attenuate the coefficients in the Fourier domain corresponding to deletion of the high frequencies contained in the noisy signal. The selection of the parameters is done by minimizing the SURE estimating the MSE in a similar way to our methodology. Due to the limited ability of ideal lowpass graph filter to separate the low-frequency noise, an adaptive weighted graph filter is proposed in Chen et al., (2020). These two graphs filtering methods can be approximately mimicked with the SWGT based one by removing the finest scales and an appropriate choice of the thresholding. However, the localization properties of SGWT both in space and frequency allow to better process signals with heterogeneous regularity in the initial domain. Intuitively, the thresholding process in the SGWT transformed domain eliminates the high frequency part of the noise where the signal is locally regular and tends to keep the peaks of the original signal whereas the Fourier transform tends to smooth them. This comes at the cost of a redundant frame in which the transformed noise is no longer stationary neither decorrelated.

Regarding the regularization method implemented in Onuki et al., (2016), the regularization parameter is also selected optimizing an MSE proxy based on a similar argument. Nonetheless, beyond the fact that the philosophy is different (regularization versus thresholding), one stress that the empirical risk bias is explicitly determined while the MSE estimation in Onuki et al., (2016) is only validated numerically. Another penalization method is given in Wang et al., (2016) that extends the approach from Tibshirani and Taylor, (2011) within the framework of graphs. For this method, the divergence term is computed explicitly; this gives rise to a data-driven parameter selection method so that this approach is an interesting concurrent to our methodology.

The paper is organized as follows. Section 2 introduces the notation and briefly reviews the notions of tight frame and SGWT of Hammond et al., (2011). Section 3 is devoted to denoising and the SURE estimator for generic thresholding process in the transformed domain. Then, the SURE is specified in the cases of coordinatewise, block thresholding processes and for correlated noise in the graph domain. In Section 4 numerical comparisons with the classical Wiener filter (oracle version) and the trend filtering introduced in Wang et al., (2016) for denoising are discussed. Several signals and graphs, including examples from real datasets, are considered. For these experimental results, the construction of the frame follows Göbel et al., (2018). In terms of denoising performance, other tight frames such as spectrum adapted and/or signal adapted tight frames from Shuman et al., (2015) and Behjat et al., (2016) might give better results. Still, Theorem 1 actually applies to any tight frame and the question of exhibiting the most efficient one is beyond the scope of the paper (see Shuman, (2020) for a comprehensive survey).

2 Spectral Graph Wavelet Transform

2.1 Graphs, Frames and Tight Frames

Let $G$ be an undirected weighted graph, with set of vertices $V$, and weights $(w_{ij})_{i,j\in V}$ satisfying $w_{ij}=w_{ji}$ for $i,j\in V$. The size of the graph is the number of nodes $n=|V|$. The (unnormalized) graph Laplacian matrix $\mathcal{L}\in\mathbb{R}^{V\times V}$ associated with $G$ is the symmetric matrix defined as $\mathcal{L}=D-W$, where $W$ is the matrix of weights with coefficients $(w_{ij})_{i,j\in V}$, and $D$ the diagonal matrix with diagonal coefficients $D_{ii}=\sum_{j\in V}w_{ij}$. A signal $f$ on the graph $G$ is a function $f:V\rightarrow\mathbb{R}$.

Let $\mathfrak{F}=\{r_{i}\}_{i\in I}$ be a frame of vectors of $\mathbb{R}^{V}$, that is a family of vectors in $\mathbb{R}^{V}$ such that there exist $A,B>0$ satisfying for all $f\in\mathbb{R}^{V}$

[TABLE]

The linear map $T_{\mathfrak{F}}:\mathbb{R}^{V}\rightarrow\mathbb{R}^{I}$ defined for $f\in\mathbb{R}^{V}$ by $T_{\mathfrak{F}}f=(\langle f,r_{i}\rangle)_{i\in I}$ is called the analysis operator. The synthesis operator is the adjoint of $T_{\mathfrak{F}}$: namely, it is the linear map $T_{\mathfrak{F}}^{\ast}:\mathbb{R}^{I}\rightarrow\mathbb{R}^{V}$ defined for a vector of coefficients $(c_{i})_{i\in I}$ by $T^{\ast}_{\mathfrak{F}}(c_{i})_{i\in I}=\sum_{i\in I}c_{i}r_{i}$. As a frame is in particular a generating family of $\mathbb{R}^{V}$, a signal $f\in\mathbb{R}^{V}$ can be recovered from its coefficients $T_{\mathfrak{F}}f$ with the help of the synthesis operator.

2.2 Construction of Tight Frames

A frame $\mathfrak{F}$ is said to be tight if $A=B=1$ in Equation (1)—the latter is then termed the Parseval identity. From now on, the frames considered are supposed to be tight. Let us recall the generic construction of such a frame (c.f. Kereta et al., (2019) for instance).

Since $\mathcal{L}$ is self-adjoint, it admits the spectral decomposition $\mathcal{L}=\sum_{\ell}\lambda_{\ell}\langle\chi_{\ell},\cdot\rangle\chi_{\ell},$ where $\lambda_{1}\geq\lambda_{2}\geq\cdots\geq\lambda_{n}=0$ denote the (ordered) eigenvalues of the matrix $\mathcal{L}$, and $(\chi_{\ell})_{1\leq\ell\leq n}$ are the associated normalized and pairwise orthogonal eigenvectors. Then, for any function $\rho:\mathrm{sp}(\mathcal{L})\rightarrow\mathbb{R}$ defined on the spectrum $\mathrm{sp}(\mathcal{L})$ of matrix $\mathcal{L}$, the functional calculus formula reads $\rho(\mathcal{L})=\sum_{\ell}\rho(\lambda_{\ell})\langle\chi_{\ell},\cdot\rangle\chi_{\ell}$. A finite collection $(\psi_{j})_{j=0,\ldots,J}$ is a finite partition of unity on the compact $[0,\lambda_{1}]$ if $\psi_{j}:[0,\lambda_{1}]\rightarrow[0,1]~{}\textrm{for all}~{}j\in\mathcal{J}~{}\textrm{and}~{}\forall\lambda\in[0,\lambda_{1}],~{}\sum_{j=0}^{J}\psi_{j}(\lambda)=1$. Given a finite partition of unity $(\psi_{j})_{j=0,\ldots,J}$, the Parseval identity implies that the following set of vectors is a tight frame:

[TABLE]

Also, following Leonardi and Van De Ville, (2013); Göbel et al., (2018), a partition of unity can be easily defined as follows: let $\omega:\mathbb{R}^{+}\rightarrow[0,1]$ be some function with support in $[0,1]$, satisfying $\omega\equiv 1$ on $[0,b^{-1}]$, for some $b>1$, and set

[TABLE]

2.3 Discrete SGWT Associated with a Partition of Unity

Let $(\psi_{j})_{j=0,\ldots,J}$ be a partition of unity of $[0,\lambda_{1}]$. The SGWT of a signal $f\in\mathbb{R}^{V}$ is given by

[TABLE]

The adjoint linear transformation $\mathcal{W}^{\ast}$ of $\mathcal{W}$ is:

[TABLE]

The tightness of the underlying frame implies that $\mathcal{W}^{\ast}\mathcal{W}=\mathrm{Id}_{\mathbb{R}^{V}}$ so that a signal $f\in\mathbb{R}^{V}$ can be recovered by applying $\mathcal{W}^{\ast}$ to its wavelet coefficients $((\mathcal{W}f)_{i})_{i=1,\ldots,n(J+1)}\in\mathbb{R}^{n(J+1)}$ (see Hammond et al., (2011)).

3 Adaptive Denoising with SGWT

Let $f\in\mathbb{R}^{V}$ be some signal on a graph $G$ and $\xi$ be an $n$-dimensional Gaussian vector distributed as $\mathcal{N}(0,\sigma^{2}\mathrm{Id})$. The aim of denoising is to recover the unknown signal $f$ from the observed noisy version $\tilde{f}=f+\xi$. Basically, our denoising procedure will consist of three steps: (1) compute the SGWT transform $\mathcal{W}\tilde{f}\in\mathbb{R}^{n(J+1)}$; (2) apply a given thresholding operator $h(\cdot)$ to the coefficients $\mathcal{W}\tilde{f}$; (3) apply the inverse SGWT transform to obtain an estimation $\hat{f}$ of the original signal. Here, the main challenge in denoising consists in choosing a suitable thresholding operator with respect to the noisy signal $\tilde{f}$ and the underlying graph. The performance measure in the sequel will be the MSE between the original signal $f$ and the denoised signal $\hat{f}$: $\|f-\hat{f}\|_{2}^{2}$. First, it is worth noting that the Parseval identity allows direct optimization in the transformed domain of wavelet coefficients. Secondly, in practice, obviously the original signal remains unknown. To overcome this difficulty, the MSE is generally substituted with the Stein’s Unbiased Risk Estimator which no longer depends on the original signals (see Donoho and Johnstone, (1995) for instance). Nonetheless, contrary to the usual wavelet transform, the white noise $\xi$ is mapped onto a correlated Gaussian noise. In the next section, the SURE is derived taking into account these correlations.

3.1 The SURE Estimator in the Transformed Domain

By linearity, the denoising problem $\widetilde{f}=f+\xi$ is transferred to the denoising problem $\widetilde{F}=F+\Xi$ with $\Xi\sim\mathcal{N}(0,\sigma^{2}\mathcal{W}\mathcal{W}^{\ast})$, $\widetilde{F}=\mathcal{W}\widetilde{f}$ and $F=\mathcal{W}f$. It is worth noting that $\mathcal{W}\mathcal{W}^{\ast}$ is an orthogonal projector of rank $n$ in $\mathbb{R}^{n(J+1)}$ (see Göbel et al., (2018) for instance). Consequently, its spectral decomposition reads $\mathcal{W}\mathcal{W}^{\ast}=U\Sigma U^{\ast}$ with $U$ a unitary matrix of $\mathbb{R}^{n(J+1)}$ and $\Sigma=\left(\begin{smallmatrix}\mathrm{Id}_{\mathbb{R}^{n}}&0\\ 0&0\\ \end{smallmatrix}\right)$.

A thresholding process is a map $h:\mathbb{R}^{n(J+1)}\rightarrow\mathbb{R}^{n(J+1)}$. Typically, the map $h$ is a coordinatewise or a block shrinkage in applications. The following result extending the SURE’s expression to correlated noise is based on the Stein’s lemma in Stein, (1981) in which $h$ is assumed to be weakly differentiable. One refers the reader to Stein, (1981) for the precise definition.

Theorem 1** ($h$-SURE).**

Let $h$ be a weakly differentiable thresholding process for the denoising problem $\widetilde{F}=F+\Xi$. Then the theoretical MSE is given by

[TABLE]

where $h_{i}$ is the $i$-th component of $h$.

It is worth noting that $\mathbf{Cov}(\Xi_{i},\Xi_{j})=\sigma^{2}(\mathcal{W}\mathcal{W}^{\ast})_{i,j}$ so that, as soon as the thresholding process $h$ is specified and the noise variance $\sigma^{2}$ estimated, the SURE of $h$ defined below can be completely computed from the noisy observations as in the regular case:

[TABLE]

Proof.

The theoretical MSE can be rewritten as follows

[TABLE]

The second term is equal to $n\sigma^{2}$ since almost surely $\|\Xi\|^{2}=\|U^{\ast}\Xi\|^{2}=\|P_{K}U^{\ast}\Xi\|^{2}$ where $K=\mathrm{ker}(\mathcal{W}\mathcal{W}^{\ast})^{\perp}$ and $P_{K}$ the orthogonal projection onto $K$. Finally, setting $g(x)=h(x)-x$, $x\in\mathbb{R}^{n(J+1)}$, it remains to compute the last term $\mathbf{E}\langle g(\widetilde{F}),\Xi\rangle=\mathbf{E}\langle g(F+\Xi),\Xi\rangle$ where $F$ is deterministic and $\Xi\sim\mathcal{N}(0,\sigma^{2}\mathcal{W}\mathcal{W}^{\ast})$. A simple computation gives

[TABLE]

Then, following Liu, (1994), each term in the sum above is given by

[TABLE]

since $\mathrm{Tr}(\sigma^{2}\mathcal{W}\mathcal{W}^{\ast})=n\sigma^{2}$. This ends the proof of Theorem 1. ∎

3.2 Coordinatewise Thresholding Process

For a coordinatewise thresholding process, the map $h$ is of the form $h(x)=(\tau(x_{i},t_{i}))_{i=1,\ldots,n(J+1)}$ where $(t_{i})_{i=1,\ldots,n(J+1)}$ are the thresholds. In practice, we may choose $\tau(x,t)=x\max\{1-t^{\beta}|x|^{-\beta},0\}$ with $\beta\geq 1$. The most popular choices are the soft thresholding ($\beta=1$), the James-Stein thresholding ($\beta=2$) and the hard thresholding ($\beta=\infty$). The latter will not be considered here since it does not lead to a sufficiently regular thresholding process for Theorem 1 to be applied.

For any $\beta\in[1,\infty)$, the derivative $\partial_{j}h_{i}$ vanishes whereas

[TABLE]

Consequently, the SURE associated with $h$ is given by

[TABLE]

The usual expression of the SURE is recovered from the identity above remarking that $\mathbf{V}[\Xi_{i}]$ are identically equal to $\sigma^{2}$ when the transformed noise is uncorrelated.

Let us notice that the coordinatewise soft thresholding ($\beta=1$) satisfies an oracle inequality as shown in Göbel et al., (2018). Similarly to the regular case, it states that up to a log factor, the soft thresholding estimator can mimic an oracle projection.

3.3 Optimization: Donoho and Johnstone’s Trick

The SURE can be optimized in the same way as in the standard case using the Donoho and Johnstone’s trick of Donoho and Johnstone, (1995) whose the justification is recalled below.

For the sake of simplicity, we first consider the case of the coordinatewise thresholding process with a uniform threshold: $t_{i}=t$ for all $i=1,\ldots,n(J+1)$. Denote by $a_{1},\ldots,a_{n(J+1)}$ the absolute values of the noisy wavelet coefficients $|\widetilde{F}_{i}|$ in the increasing order. The trick comes from the observation that, on each interval $(a_{k},a_{k+1})$, the last term of Equation (2) is non-decreasing whereas the second term is an increasing function of $t$. Consequently, the SURE hits its minimum at some value $a_{k^{\ast}}$, $k^{\ast}=1,\ldots,n(J+1)$.

If the thresholds $t_{i}$ are no longer uniform but merely tied inside blocks with values $t_{1},\ldots,t_{L}$, the same trick is still valid: group the terms in the sums along the different parameters $t_{1},\ldots,t_{L}$ and optimize each partial sum with respect to $t_{k}$, $k=1,\ldots,L$.

3.4 Block Thresholding Process

In order to take advantage of the localization properties of SGWT and the regularity of the original signal, we may introduce block thesholding processes similar to Cai, (1999).

Consider a partition $(B_{\ell})_{\ell\in L}$ of $\{1,\ldots,n(J+1)\}$ and set $\|x\|_{B_{\ell}}^{2}=\sum_{i\in B_{\ell}}(x_{i})^{2}$. In this case, the thesholding process $h=(h_{i})_{i=1,\ldots,n(J+1)}$ reads

[TABLE]

If $i,j$ are in different blocks, then $\partial_{j}h_{i}$ vanishes. Additionally, if $i,j$ are in $B_{\ell}$ but $i\neq j$ then

[TABLE]

whereas

[TABLE]

Consequently, a straightforward computation leads to

[TABLE]

Once again, for uncorrelated transformed noise, the usual expression easily follows from the identity above. Note also that the optimization of the SURE in this case requires more sophisticated techniques as the divergence term is no longer monotone.

3.5 Correlated Noise in the Graph Domain

The SURE can also be stated in the context of correlated noise at the cost of some prior information on the covariance structure. More precisely, in the denoising problem $\widetilde{f}=f+\xi$ with correlated noise, it is supposed that $\xi\sim\mathcal{N}(0,\Gamma)$ for some covariance matrix $\Gamma$. The denoising problem reads in the transformed problem as $\widetilde{F}=F+\Xi$ with $\Xi\sim\mathcal{N}(0,\mathcal{W}\Gamma\mathcal{W}^{\ast})$.

Corollary 1**.**

Under the assumption of Theorem 1, the theoretical $\mathrm{MSE}$ is given by

[TABLE]

where $h_{i}$ is the $i$-th component of $h$.

Let us point out that the parameters selection can be made without computing explicitly $\mathrm{Tr}(\mathcal{W}\Gamma\mathcal{W}^{\ast})$ since it does not depend on $h$—even though, the $\mathrm{MSE}$ estimate is obviously shifted by this quantity. Besides, the correlation structure $\Gamma$ is actually hidden in the quantities $\mathbf{Cov}(\Xi_{i},\Xi_{j})$, namely, for $1\leq i,j\leq n(J+1)$: \mathbf{Cov}(\Xi_{i},\Xi_{j})=\big{(}\mathcal{W}\Gamma\mathcal{W}^{\ast}\big{)}_{i,j}. Consequently, computationally speaking, there is no additional burden compared to the white noise case.

Proof.

The proof follows the lines of Theorem 1 with \mathbf{E}[\Xi_{i}^{2}]=\big{(}\mathcal{W}\Gamma\mathcal{W}^{\ast}\big{)}_{i,i},1\leq i\leq n(J+1). ∎

In applications, it is usually reasonable to assume some structure on the covariance matrix $\Gamma$ reflecting the topology of the underlying graph. Typically, the noise on two given vertices may be correlated if those vertices are close enough in the graph. For example, let $\xi_{0}\sim\mathcal{N}(0,\sigma^{2}\mathrm{Id})$, we set $\xi=\xi_{0}+\alpha W\xi_{0}$ where $W$ is the graph matrix of weights and $\alpha>0$ some tuning parameter describing the global intensity of the correlation. Then, it follows,

[TABLE]

Other choices are obviously possible.

3.6 Complexity

Regarding the space complexity, we need to store the frame and the weights appearing in the SURE for a cost of $O(n^{2}(J+1)^{2})$. With given Laplacian eigendecomposition, the time complexity of the optimization of the SURE is (in average) of order $O(n(J+1)\log(n(J+1)))$ for the coordinate-wise estimator following [8]. For the block estimator, the use of a grid search is a limitation.

4 Numerical Results

This section presents the empirical performance of the proposed automatic threshold selection for different signals defined on different graphs: the Minnesota roads graph (seen as a reference in many recent studies, see Behjat et al., (2016) and references therein) with synthetic signals and the Facebook graph with signals from Wang et al., (2016), the Pittsburgh Census Tract graph, a graph built from a dataset on New York City taxis with a real signal as well as numerical experiments in the correlated and block cases. All the experiments are conducted with the R package gasper.

4.1 The Minnesota Roads Graph

The Minnesota roads graph is a planar graph consisting of 2642 vertices and 6606 edges. Each vertex is described by its $(x,y)$-coordinates. The function $\omega$ chosen in the experiments is a piecewise linear function with support in $[0,1]$ and constant equal to 1 on $[0,b^{-1}]$ with $b=2$. From $\lambda_{1}\approx 6.89$, we deduce that the number of scales is $J+1=5$—see Section 2.

On this graph, two classes of synthetic signals are generated inspired by the methodology introduced in Behjat et al., (2016). Let us briefly recall the construction: let $\eta\in(0,1)$ and $k\in\mathbb{N}$ be two parameters; a signal $f_{\eta,k}$ is obtained by letting the adjacency matrix $W$ acts on an i.i.d. realization $x_{\eta}$ of Bernoulli random variables of parameter $\eta$, in symbols $f_{\eta,k}=W^{k}x_{\eta}/\lambda_{1}^{k}$. This method generates signals with different regularities. This has to be understood in the sense of the graph topology and not that one given by the embedding space $\mathbb{R}^{2}$. In the experiment, two signals are generated with parameters $\eta=0.01,k=2$ and $\eta=0.001,k=4$ respectively (see Figure 1 left and right respectively).

We compare the performance in terms of SNR (computed on the functions after reconstruction) for different denoising strategies, different noise levels and each synthetic signal $f_{\eta,k}$. For each noise level $\sigma=0.005,0.01,0.02$ and $\sigma=0.001,0.002,0.004$, a sample of $N=10$ white Gaussian noise is simulated and a global (G) versus level-dependent (LD) coordinatewise thresholding are performed with respect to the soft ($\beta=1$) and James-Stein ($\beta=2$) thresholding rules. For each strategy, we compare the average behavior of the SNR for parameter selected with the oracle ($\mathrm{MSE}^{\beta=1,2}$ obtained by minimizing the MSE using the original signal $f$) and the SURE with known $\sigma$ ($\mathrm{SURE}_{\sigma}^{\beta=1,2})$. Also, the standard deviation on the sample is provided.

These results are first compared to the classical Wiener filter. More precisely, the Wiener filter consists of attenuating the Fourier coefficient $\mathcal{F}(\widetilde{f})$ of $\widetilde{f}$. Below, we only consider the oracle linear attenuation $\mathcal{F}(\widetilde{f})[i]{\mathcal{F}(f)[i]}/({\mathcal{F}(f)[i]^{2}+\sigma^{2}})$. While this estimator is unrealistic since it depends on $f$, any Wiener filter has worse performance than this oracle. Table 1 is completed by the performance of the Wiener filter on each signal. Also, the theoretical value $r_{\inf}$ of the oracle risk given in Mallat, (2009) is recalled for comparison purpose.

Our methodology is also compared to the so-called graph trend filtering (i.e. for $k=0,1,2$) introduced in Wang et al., (2016). The graph trend filtering is a regularization method with a penalty term involving the graph difference operator at a given order (see Wang et al., (2016)). In the experiments, we make use of the matlab toolbox gtf333Available here: https://sites.cs.ucsb.edu/ yuxiangw/resources.html provided by the authors of Wang et al., (2016).

Generally speaking, we observe from Table 1 that our method performs better than the trend filtering motivating the use of multiscale analysis. This idea is confirmed by the comparison with the Wiener filter, in particular in the lower $\mathrm{SNR}$ regime that is for higher noise levels. Also, similarly to the regular case, numerical experiments shows that the James-Stein threshold ($\beta=2$) is slightly more efficient than the soft threshold in particular for the global thresholding process.

Let us point out that there is no fundamental difference in terms of performance between the global and level dependent thresholding in this experiment. In fact, the level dependent thresholding always performs at least as good as the global one. Since the additional computational cost is acceptable, the level dependent thresholding appears to be a good choice without any further a priori knowledge.

4.2 The Facebook Graph

Here we examined and compare the denoising performance of level dependent SGWT thresholding (LD) against the trend filtering and Laplacian smoothing Smola and Kondor, (2003) on a nonplanar graph considered in Wang et al., (2016): the Facebook graph from the Stanford Network Analysis Project444http://snap.stanford.edu/data/ego-Facebook.html. This undirected graph, collected from survey participants using this Facebook app, is composed of 4039 nodes representing Facebook users, and 88,234 edges representing friendships (see Leskovec and Mcauley, (2012) for more details). For signal $f$, we consider the different regularities used in Wang et al., (2016) as well as the same noise levels, for 5 realizations (see (Wang et al.,, 2016, Section 5.1) for more details). More precisely, we simply run the scripts provided in the matlab toolbox provided by the authors of Wang et al., (2016). The results are shown in the Figure 2 (to be compared with (Wang et al.,, 2016, Figure 5 p.14 and Figure 9 p.27)).

With the exception of the dense Poisson case, where all methods provide comparable results, LD globally provides better performance than trend filtering (whatever the value of $k$), especially at high noise levels where the maximal gain in terms of SNR is greater than 15dB and exceeds 5dB and 10dB respectively for the 5 highest noise levels for inhomogeneous and homogeneous random walk cases. For this graph, the CPU times associated with trend filtering (for $k=1$ and $k=2$) for one type of signal and 5 realizations (and a 51-point grid search) are of the order of 3 to 4 days (depending on the case), for LD around 25 minutes (including diagonalization and frame calculation which only need to be calculated once).

The calibration of certain parameters has not been studied. However, the latter can considerably influence the performance of the SGWT. For the homogeneous random walk case, we examine here the influence of the number of scales retained for the construction of the frame and controller by $b$. For $b=5,4,3$, the frame contains respectively $7,8,9$ scales. The results are shown in Figure 2, where it can be seen that the frame composed of 7 scales produces the best results. Note that these performances might be improved, for example by making the SURE depend on the $\beta$ parameter characterizing the threshold rule. The frame considered here has only one parameter, other more flexible constructions, based on a partition of the unit or other types of tight frames such as spectrum adapted and/or signal adapted tight frames from Shuman et al., (2015) and Behjat et al., (2016) could also lead to an improvement.

4.3 Pittsburgh Census Tract Graph

For the sake of completeness, our methodology is also compared to the trend filtering on the Pittsburgh Census Tract graph considered in Wang et al., (2016) which consists of 402 vertices and 2382 edges. The very same piecewise linear function $\omega$ with $b=2$ is still used for this graph. The number of scales is then $J+1=7$. For this experiment, only the level dependent threshold procedure is considered.

We consider the signal and 10 realizations of the noisy signal generated in Wang et al., (2016) (corresponding to an average noise level of $4.84\pm 0.37\mathrm{dB}$). The resulting SNR for the oracles of SGWT and fused lasso (i.e. $k=0$) are respectively $11.51\pm 0.52\mathrm{dB}$ and $9.85\pm 0.54\mathrm{dB}$.

Additionally, we run a comparison with the trend filtering (i.e. for $k=0,1,2$) for the signal $f_{\eta,k}$ with $\eta=0.01$ and $k=5$ with the different noise levels $\sigma=0.004$, $\sigma=0.005$ and $\sigma=0.01$. A comparison with another wavelet estimator proposed in Sharpnack et al., (2013) is also provided, considering two thresholding rules (i.e. “soft” and “hard”). For these 5 competitors we only report the oracles results.

Even though the SURE no longer depends on the original signal, it does depend on $\sigma^{2}$ in both methodologies. Since in real applications, the noise level remains unknown in general, we introduce two naive estimators of $\sigma$. In fact, a straightforward computation shows that for any function $g:\mathbb{R}_{+}\rightarrow\mathbb{R}_{+}$:

[TABLE]

so that a biased estimator of $\sigma^{2}$ is given by

[TABLE]

As soon as the original signal is reasonably smooth so that $f^{T}g(\mathcal{L})f$ is negligible compared to $\mathrm{Tr}~{}g(\mathcal{L})$, then $\hat{\sigma}^{2}$ is an accurate enough estimation of $\sigma^{2}$. As a first choice, we choose $g(x)=x$. Thanks to Dirichlet’s formula, it follows:

[TABLE]

This is nothing but the graph analogue of the Von Neumann estimator of von Neumann, (1941) explaining the terminology Graph Von Neumann estimator (GVN).

A second natural choice is given by $g(x)=\psi_{J}(x)$ corresponding to the filter at the finest scale. The resulting estimator is called High Pass Filter Von Neumann (HPFVN). The value of the estimator is easily computed from the coefficients as follows:

[TABLE]

Table 2 summarizes the findings with the nomenclature of Table 1 (where LD stands for level-dependent (LD) coordinatewise thresholding) of the main document. Lines $\mathrm{SURE}_{\hat{\sigma}_{1,2}}^{\beta=1,2}$ stand for the SURE procedure in which the noise level is estimated by the GVN ($\hat{\sigma}_{1}$) and the HPFVN ($\hat{\sigma}_{2}$). HPFVN and GVN provide a very comparable sigma estimate in this setting. Additionally, a visual comparison of our methodology with the fused lasso is illustrated in Figure 4 (corresponding to one realization in the context of the third column of the Table 2). We can see that our approach provides a gain of about 2.5dB compared to the fused lasso.

Again, in these experiments, the multiscale analysis shows better performances than the trend filtering. Besides, the estimation of $\sigma$ is sufficiently accurate to have a fully data-driven procedure.

To conclude, let us stress that the SGWT is computationally more efficient than the GTF ($k=1,2$). For the latter, on a standard laptop (Intel Core i7@2.7GHz-16Go LP-DDR3@2133MHz), each 10 realizations consumed about 4h16m cpu time in mean (min:1h45m, max:6h19m) with the same grid search as for the Pittsburg in Wang et al., (2016). The CPU consumption for $k=0$ is more decent with a mean of 3m for each 10 realizations exploiting the idea and the c++ code of Chambolle and Darbon, (2009). Incidentally, this main drawback of GTF was noticed in Padilla et al., (2018) forcing a preprocessing of the graph using Depth-first search algorithm. The SGWT on its side consumed 3min for the diagonalization and 42s for each 10 realizations.

4.4 Real Dataset: New York City Taxis

Our methodology has been also tested on a real data fetched from NYC taxis555https://s3.amazonaws.com/nyc-tlc/trip+data/yellow_tripdata_2018-01.csv databases. We build a graph with 265 vertices consisting of the LocationID (Pick-Up and Drop-Off) and define Gaussian weights $w_{ij}=\exp(-\tau d_{i,j}^{2})$ where $d_{i,j}$ is the mean distance taken on all the trips between $i$ and $j$ or $j$ and $i$. The signal $f$ considered is defined upon the variable “total amount” on which an artificial noise is added. For an average input $\mathrm{SNR}_{\mathrm{in}}=5.23\pm 0.38\mathrm{dB}$ on $N=25$ observations, we obtain for the level dependent SURE with $\beta=2$ and $\sigma$ known, an output $\mathrm{SNR}$ of $10.41\pm 0.62\mathrm{dB}$ compared to the performance of the oracle Wiener filter, $\mathrm{SNR}=7.50\pm 0.40\mathrm{dB}$ and the oracle fused lasso, $\mathrm{SNR}=5.69\pm 0.43\mathrm{dB}$.

4.5 Correlated Noise

On the Minnesota roads graph and for the signal $f_{\eta,k}$ of Section 4 with $\eta=0.01$ and $k=2$, we add one realization of a correlated noise with covariance matrix given by (3) where $\alpha=0.5$ is added leading to a noisy signal with $\mathrm{SNR}_{\mathrm{in}}=2.07\mathrm{dB}$. We run the level-dependent coordinate-wise thresholding process. The $\mathrm{SNR}$ given by the oracle involving the unknown signal $f_{\eta,k}$ and the one given by the SURE estimator adapted to correlated noise are similar: $8.54\mathrm{dB}$ versus $8.51\mathrm{dB}$. In this case, the SURE for uncorrelated noise shows very bad performances: we found $4.33$ for the corresponding SNR.

The quality of the SURE for uncorrelated noise is closely related to the intensity of the correlation tuned by the parameter $\alpha$. As an example, if we choose $\alpha=0.1$, the SURE adapted to correlated noise still performs very well with an estimated SNR of $10.44\mathrm{dB}$ compared to the oracle $10.46\mathrm{dB}$. The SURE for uncorrelated noise is nonetheless not that bad since it estimates the SNR at $9.78\mathrm{dB}$.

Consequently, the SURE estimate for uncorrelated noise is robust to small correlations which is particularly interesting in applications since it can be difficult to estimate the correlation structure.

4.6 Further Experiments with Block Thresholding

Finally, we report some experimental results in the context of block thresholding in the setting of Section 4.3. For each scale $j=0,1,\ldots,7$, the $n$ wavelet coefficients are split into $L$ blocks of uniform length (except for the last block that can be shorter). The best performance of block thresholding is achieved for blocks of size $|L|=47$. The $\mathrm{SNR}_{\mathrm{in}}=9.64\pm 0.33\mathrm{dB}$ for 25 realizations. For the oracle global coordinate-wise threshold with $\beta=2$, we obtain a SNR of $11.07\pm 0.33\mathrm{dB}$ compared to the block procedure with a uniform threshold $11.82\pm 0.40\mathrm{dB}$. The block procedure performs better than the coordinate-wise one for a uniform threshold but is actually worse compared to the level dependent coordinate-wise thresholding process. The level-dependent method is expected to give better performance, but would require a more sophisticated optimization algorithm than grid search to be computationally acceptable.

5 Conclusion and Perspectives

In this paper, we have introduced a version of the SURE designed for SGWT, allowing automatic parameter selection in denoising tasks of signals on graphs. Closed-form expressions for coordinatewise and block SUREs have been provided for a wide range of threshold rules. Finally, the case of a correlated noise in the graph has also been considered. Many experiments on the Minnesota graph, the Facebook graph built, the Pittsburgh graph and the NYC taxis graph from real data are conducted.

For signals of different regularity on those graphs, it has been shown that the SURE provides an efficient estimate of the theoretical MSE. These experiments also show that multi-scale analysis is a serious competitor to existing methods. Indeed, the SGWT shows performances equivalent and sometimes even much better than the GTF, especially at high noise levels.

To be complete, the variance parameter $\sigma^{2}$ should be estimated. This has been (partially) addressed in the supplementary material by introducing the two estimators HPFVN and GVN. The performance of these estimators highly depends on the underlying signal. Further investigations seem to be necessary.

Theoretically, the coordinatewise soft thresholding in the transformed domain satisfies an oracle inequality as shown in Göbel et al., (2018). Similarly to the regular case, this oracle inequality states that the estimator mimics the oracle projection up to a log factor. The proof relies on the fact that the multivariate risk is expressed as a sum of univariate risks so that the Donoho’s machinery applies. Using this fact, a maximal inequality for the SURE might be stated as well.

Regarding, the numerical complexity, the main limitation is the need of a complete reduction of the Laplacian. In the same vein as Hammond et al., (2011), many of the involved steps might be numerically optimized using Chebyshev polynomials. Actually, the only problematic step in the method is the computation of the weights $(\mathcal{W}\mathcal{W}^{\ast})_{i,j}$ appearing in the SURE. However, their expression in terms of covariance suggests that Monte-Carlo estimation could work. Besides, the space-time complexity might be reduced taking advantage of the low-rank property of $\mathcal{W}\mathcal{W}^{\ast}$ implying several linear constraints on the weights (precisely $nJ$). Finally, for the block procedure to be completely useful, an adapted optimization algorithm should be implemented.

Some questions not addressed in the paper remains open. As already announced in introduction of the paper, the choice of a suitable frame for different graphs and different families of signals is still an open problem in spite of advances in recent years. Most likely, good choices of frame should involve a notion of graph limit such as the one introduced for graphons (see Lovász, (2012)) or the more probabilistic Benjamin-Schramm limit introduced in Benjamini and Schramm, (2001). This formal study should also give rise to less naïve estimators of the noise level and above all give recommendations according to the class of signals considered.