Localized Fourier Analysis for Graph Signal Processing

TL;DR

This paper introduces a localized Fourier analysis approach on graphs that partitions the spectrum to improve theoretical understanding and practical tasks like noise estimation and denoising, with efficient computation methods.

Contribution

It proposes a new spectral partitioning method for graph Fourier analysis, extending spectral graph wavelets to localized Fourier domain wavelets called LocLets, and demonstrates practical benefits.

Findings

Significant improvements in noise level estimation.

Enhanced signal denoising performance.

Efficient computation with same complexity as standard methods.

Abstract

We propose a new point of view in the study of Fourier analysis on graphs, taking advantage of localization in the Fourier domain. For a signal $f$ on vertices of a weighted graph $\mathcal{G}$ with Laplacian matrix $\mathcal{L}$, standard Fourier analysis of $f$ relies on the study of functions $g(\mathcal{L})f$ for some filters $g$ on $I_\mathcal{L}$, the smallest interval containing the Laplacian spectrum ${\mathrm sp}(\mathcal{L}) \subset I_\mathcal{L}$. We show that for carefully chosen partitions $I_\mathcal{L} = \sqcup_{1\leq k\leq K} I_k$ ($I_k \subset I_\mathcal{L}$), there are many advantages in understanding the collection $(g(\mathcal{L}_{I_k})f)_{1\leq k\leq K}$ instead of $g(\mathcal{L})f$ directly, where $\mathcal{L}_I$ is the projected matrix $P_I(\mathcal{L})\mathcal{L}$. First, the partition provides a convenient modelling for the study of theoretical properties of Fourier analysis and allows for new results in graph signal analysis (\emph{e.g.} noise level estimation, Fourier support approximation). We extend the study of spectral graph wavelets to wavelets localized in the Fourier domain, called LocLets, and we show that well-known frames can be written in terms of LocLets. From a practical perspective, we highlight the interest of the proposed localized Fourier analysis through many experiments that show significant improvements in two different tasks on large graphs, noise level estimation and signal denoising. Moreover, efficient strategies permit to compute sequence $(g(\mathcal{L}_{I_k})f)_{1\leq k\leq K}$ with the same time complexity as for the computation of $g(\mathcal{L})f$.

Tables

Table 1. Table 1: SNR performance for Swissroll ( n = 1000 𝑛 1000 n=1000 , K = 22 𝐾 22 K=22 ).

signal	$\sigma$	${\mathrm{SNR}}_{\mathrm{in}}$	$M_{\mathrm{PF}}$	$M_{\mathrm{LLet}}$	$M_{\mathrm{LLet}+\mathrm{PF}}$	${\mu }_{\mathrm{PF}}$	${\mu }_{\mathrm{LLet}}$	${\mu }_{\mathrm{LLet}+\mathrm{PF}}$
$f_{951-1000}$	0.005	16.195	17.557	20.580	20.528	17.361	20.035	19.974
$f_{501-550}$	0.005	15.859	18.298	8.244	20.821	18.044	8.140	20.245
$f_{951-1000}$	0.01	10.267	12.183	15.564	15.701	10.433	13.652	13.760
$f_{501-550}$	0.01	10.178	13.121	7.879	16.204	11.165	7.646	14.518
$f_{951-1000}$	0.015	6.763	9.430	12.661	13.129	8.961	12.127	12.388
$f_{501-550}$	0.015	6.362	9.898	7.611	14.159	9.540	7.481	13.398

Table 2. Table 2: SNR performances for denoising task.

matrix	signal	$\sigma$	${\mathrm{SNR}}_{\mathrm{in}}$	$M_{\mathrm{PF}}$	$M_{\mathrm{LLet}+\mathrm{PF}}$	${\mu }_{\mathrm{PF}}$	${\mu }_{\mathrm{LLet}+\mathrm{PF}}$
Si2	$f_{720-769}$	0.005	17.104	22.344	26.973	21.849	25.170
($n=762$, $K=22$)	$f_{370-419}$	0.005	16.820	18.034	21.778	17.813	20.673
	$f_{720-769}$	0.01	11.408	16.821	22.501	16.572	20.558
	$f_{370-419}$	0.01	11.175	12.444	19.740	12.151	17.121
	$f_{720-769}$	0.02	4.826	11.476	15.695	11.104	14.711
	$f_{370-419}$	0.02	5.047	7.354	13.542	6.925	12.677
Minnesota	$f_{2593-2642}$	0.004	13.599	17.839	20.717	17.672	20.035
($n=2642$, $K=22$)	$f_{1343-1392}$	0.004	13.741	15.999	20.388	15.822	19.234
	$f_{2593-2642}$	0.005	11.681	16.086	19.342	15.830	18.417
	$f_{1343-1392}$	0.005	11.916	14.459	18.392	14.298	18.029
	$f_{2593-2642}$	0.01	5.911	10.875	14.143	10.605	13.556
	$f_{1343-1392}$	0.01	5.843	9.660	11.762	9.409	10.952
Cage9	$f_{3485-3534}$	0.003	15.016	20.477	9.700	20.216	9.664
($n=3534$, $K=22$)	$f_{1785-1834}$	0.003	15.014	15.876	16.945	15.798	16.799
	$f_{3485-3534}$	0.005	10.503	17.290	18.410	16.772	18.185
	$f_{1785-1834}$	0.005	10.507	12.118	13.032	12.035	12.898
	$f_{3485-3534}$	0.009	5.423	13.002	13.264	12.763	12.898
	$f_{1785-1834}$	0.009	5.395	8.329	10.468	8.168	9.827
Saylr4	$f_{3515-3564}$	0.003	14.871	23.108	24.516	23.040	24.117
($n=3564$, $K=22$)	$f_{2015-2064}$	0.003	15.069	21.365	23.903	21.010	23.412
	$f_{3515-3564}$	0.005	10.478	19.135	20.966	18.943	20.268
	$f_{2015-2064}$	0.005	10.635	17.016	19.610	16.662	18.732
	$f_{3515-3564}$	0.009	5.420	15.277	17.070	14.791	16.567
	$f_{2015-2064}$	0.009	5.480	12.055	14.802	11.830	14.031
Net25	$f_{9471-9520}$	0.006	4.682	5.171	5.319	5.094	5.205
($n=9520$, $K=5$)	$f_{4971-5020}$	0.006	4.577	5.811	6.035	5.714	5.933
	$f_{9471-9520}$	0.007	3.282	5.287	5.416	5.127	5.251
	$f_{4971-5020}$	0.007	3.406	5.267	5.451	5.104	5.266
	$f_{9471-9520}$	0.008	2.208	5.171	5.268	4.928	5.034
	$f_{4971-5020}$	0.008	2.153	4.450	4.594	4.319	4.471
Net25	$f_{9471-9520}$	0.006	4.495	5.319	6.004	5.190	5.834
($n=9520$, $K=25$)	$f_{4971-5020}$	0.006	4.574	5.849	4.752	5.714	4.579
	$f_{9471-9520}$	0.007	3.307	5.264	5.765	4.951	5.447
	$f_{4971-5020}$	0.007	3.256	5.132	4.182	5.014	4.102
	$f_{9471-9520}$	0.008	2.111	4.913	5.249	4.756	5.128
	$f_{4971-5020}$	0.008	2.150	4.345	3.670	4.251	3.554

Equations

$$g(\mathcal{L})=\sum _{1\le \ell \le n}g({\lambda }_{\ell })⟨{\chi }_{\ell },\cdot ⟩{\chi }_{\ell }.$$

$$P_I(\mathcal{L})=\sum _{1\le \ell \le n}{1}_I({\lambda }_{\ell })⟨{\chi }_{\ell },\cdot ⟩{\chi }_l=\sum _{\ell :{\lambda }_{\ell }\in I}⟨{\chi }_{\ell },\cdot ⟩{\chi }_{\ell }.$$

$$g(\mathcal{L})=\sum _{1\le k\le K}g({\mathcal{L}}_{I_k}).$$

$$\mathcal{W}f=(\phi (\mathcal{L})f^T,\psi (s_1\mathcal{L})f^T,\dots ,\psi (s_J\mathcal{L})f^T{)}^T.$$

$${\mathcal{W}}^{\ast }({\eta }_0^T,{\eta }_1^T,\dots ,{\eta }_J^T{)}^T=\phi (\mathcal{L}){\eta }_0+\sum _{j=1}^J\psi (s_j\mathcal{L}){\eta }_j.$$

$$({\mathcal{W}}^{\ast }\mathcal{W}{)}^{-1}{\mathcal{W}}^{\ast }(c_{j,m}{)}_{j,m},$$

$$\psi (s\mathcal{L})f(x)=\frac{1}{2\pi }{\int }_{-\infty }^{+\infty }{\hat{\psi }}^{\ast }(s\omega )\hat{f}(\omega )e^{i\omega x}d\omega ,$$

$$\psi (s\mathcal{L})f(n)=\sum _{\ell }\psi (s{\lambda }_{\ell })\hat{f}(\ell ){\chi }_{\ell }(n).$$

$$\begin{array}{rl}& {\phi }_{m,I}=\phi ({\mathcal{L}}_I){\delta }_m=\sum _{\ell :{\lambda }_{\ell }\in I}\phi ({\lambda }_{\ell }){\hat{\delta }}_m(\ell ){\chi }_{\ell }\\ & {\psi }_{j,m,I}=\psi (s_j{\mathcal{L}}_I){\delta }_m=\sum _{\ell :{\lambda }_{\ell }\in I}\psi (s_j{\lambda }_{\ell }){\hat{\delta }}_m(\ell ){\chi }_{\ell },\end{array}$$

$${\mathcal{W}}^{I_k}f=(\phi ({\mathcal{L}}_{I_k})f^T,\psi (s_1{\mathcal{L}}_{I_k})f^T,...{)}^T,1\le k\le K.$$

$${\mathcal{W}}^{I_k\ast }({\eta }_0^T,{\eta }_1^T,\dots ,{\eta }_J^T{)}^T=\phi ({\mathcal{L}}_{I_k}){\eta }_0+\sum _{j=1}^J\psi (s_j{\mathcal{L}}_{I_k}){\eta }_j,1\le k\le K.$$

$$W_{ij}=k\left(\frac{\Vert x_i-x_j{\Vert }_2^2}{2\epsilon }\right),$$

$$\mathbb{E}({\tilde{f}}^{T}g({\mathcal{L}}_I)\tilde{f})=f^{T}g({\mathcal{L}}_I)f+{\sigma }^2\mathrm{Tr}(g({\mathcal{L}}_I)).$$

$$\Vert {\tilde{f}}_k{\Vert }_2^2=\Vert f_k{\Vert }_2^2+\Vert {\xi }_k{\Vert }_2^2+2⟨f_k,{\xi }_k⟩,$$

$${\hat{\sigma }}_{\mathrm{mean}}(c{)}^2=\frac{1}{|\{k:I_k\subset \overline{I_f}\}|}\sum _{k:I_k\subset \overline{I_f}}{c}_k\text{and}{\hat{\sigma }}_{\mathrm{med}}(c{)}^2={\mathrm{median}}_{k:I_k\subset \overline{I_f}}(c_k).$$

$$\mathbb{P}\left({\hat{\sigma }}_{\mathrm{mean}}(c{)}^2-{\sigma }^2\ge t\right)\le \exp \left(-\frac{K_f^2{t}^2}{V_f(1+B_f+\sqrt{1+\frac{2B_f{K}_{f}t}{V_f}})}\right),$$

$$\mathbb{P}\left({\hat{\sigma }}_{\mathrm{mean}}(c{)}^2-{\sigma }^2\le -t\right)\le \exp \left(-\frac{K_f^2{t}^2}{2V_f}\right);$$

$$\mathbb{P}\left(\widehat{\sigma}^{2}_{\rm med}\geq\beta^{-1}\sigma^{2}+2\sigma^{2}\beta^{-1}t\right)\leq\exp\left(\frac{K_{f}}{2}\ln\Big{[}4p^{+}(t)(1-p^{+}(t))\Big{]}\right),$$

$$\mathbb{P}\left(\widehat{\sigma}^{2}_{\rm med}\leq\beta\sigma^{2}-\sigma^{2}\beta t\right)\leq\exp\left(\frac{K_{f}}{2}\ln\Big{[}4p^{-}(t)(1-p^{-}(t))\Big{]}\right),$$

$$p^{+}(t)=\mathbb{P}({\Gamma }_{n_{\infty }}\ge n_{\infty }+2n_{\infty }t)\text{and}{p}^{-}(t)=\mathbb{P}({\Gamma }_{n_0}\le n_0-n_{0}t).$$

$$\mathbb{P}(|⟨f_k,{\xi }_k⟩|\ge t_{\alpha ,\sigma }\Vert f_k{\Vert }_2)\le \alpha \text{where}{t}_{\alpha ,\sigma }=\sigma \times \sqrt{-2\ln \left(\frac{\alpha }{4}\right)}.$$

$$\frac{\Vert f_k{\Vert }_2^2+2t_{\alpha ,\sigma }\Vert f_k{\Vert }_2}{{\sigma }^2}\ge {\Phi }_{\frac{n_k}{n_{\ell }}{\Gamma }_{n_{\ell }}-{\Gamma }_{n_k}}^{-1}\left(1-\frac{3\alpha }{2}\right).$$

$$b_k=\frac{\Vert {\xi }_k{\Vert }_2^2+\Vert f_k{\Vert }_2^2+2⟨{\xi }_k,f_k⟩}{n_k}\text{and}{b}_{\ell }=\frac{\Vert {\xi }_{\ell }{\Vert }_2^2}{n_{\ell }}$$

$${\hat{\sigma }}_{\mathrm{mean}}^r(c{)}^2=\frac{1}{|\{k:I_k\subset \overline{J_f}(r)\}|}\sum _{k:I_k\subset \overline{J_f}(r)}{c}_k\text{and}{\hat{\sigma }}_{\mathrm{med}}^r(c{)}^2={\mathrm{med}}_{k:I_k\subset \overline{J_f}(r)}(c_k).$$

$$\mathbb{E}\left[\big{\|}f-\widetilde{f}_{I}\big{\|}_{2}^{2}\right]=\mathbb{E}\left[\big{\|}f-\widetilde{f}\big{\|}_{2}^{2}\right]-\sigma^{2}\big{|}\overline{I}\cap{\rm sp}(\mathcal{L})\big{|}.$$

$$p_k=\mathbb{P}({\sigma }^2{\Gamma }_{n_k}>\Vert {\tilde{f}}_k{\Vert }_2^2)\text{and}{\Gamma }_{n_k}\sim {\chi }^2(n_k);$$

$$\frac{\Vert f_k{\Vert }_2}{\sigma }\left(\frac{\Vert f_k{\Vert }_2}{\sigma }-2\frac{t_{\alpha /2,\sigma }}{\sigma }\right)\ge {\Phi }_{{\Gamma }_{n_k}-{\Gamma }_{n_k}^{\prime }}^{-1}\left(1-\frac{\alpha }{2}\right),$$

$${\left(\frac{\Vert f_k{\Vert }_2+\sqrt{t_{\beta ,k}}}{\sigma }\right)}^2<{\Phi }_{{\Gamma }_{n_k}}^{-1}\left(1-\frac{\alpha }{1-\beta }\right).$$

$$\Vert f-f_I{\Vert }_2^2\le |\{k,I_k\subset I_f,p_k>\alpha \}|{\left(t_{\alpha /2,\sigma }+\sqrt{t_{\alpha /2,\sigma }^2+{\left(\sigma {\Phi }_{{\Gamma }_{n_1}-{\Gamma }_{n_1}^{\prime }}^{-1}\left(1-\frac{\alpha }{2}\right)\right)}^2}\right)}^2.$$

$$\mathbb{E}\Vert f_I-{\tilde{f}}_I{\Vert }_2^2=|\{k,p_k\le \alpha \}|n_1{\sigma }^2.$$

Full text

Keywords: Nonparametric regression; Multiscale statistics; Variance estimation; Concentration inequalities; Graph signal processing; Spectral graph theory; Graph Laplacian; Harmonic analysis on graphs

1 Introduction

Graphs provide a generic representation for modelling and processing data that reside on complex domains such as transportation or social networks. Numerous works combining both concepts from algebraic and spectral graphs with those from harmonic analysis (see for example [11, 13, 3] and references therein) have allowed to generalize fundamental notions from signal processing to the context of graphs thus giving rise to Graph Signal Processing (GSP). For an introduction to this emerging field and a review of recent developments and results see [48] and [40]. In general, two types of problems can be distinguished according to whether the underlying graph is known or unknown. The first case corresponds to the setup of a sampled signal at certain irregularly spaced points (intersections of a transportation network, nodes in a computer network, …). In the second case, a graph is constructed from the data itself, it is generally interpreted as a noisy realization of one or several distributions supported by a submanifold of the Euclidean space. In this latter context, the theoretical submanifold is somehow approximated using standard methods such as $k$-NN, $\varepsilon$-graph and their Gaussian weighted versions. In any of these cases, the framework is actually similar: it consists of a graph (given by the application or by the data) and signals are real-valued functions defined on the vertices of the graph.

Notions of graph Fourier analysis for signals on graphs were introduced and studied over the past several years [48, 45, 49, 46]. The graph Fourier basis is given by the eigenbasis $(\chi_{\ell})_{\ell}$ of the Laplacian matrix $\mathcal{L}$. The Graph Fourier Transform (GFT) consists in representing a signal $f$ in the Fourier basis $(\langle f,\chi_{\ell}\rangle)_{\ell}$, and by analogy with the standard case, the eigenvalues of $\mathcal{L}$ play the role of frequencies. From this definition, it follows that many filtering techniques are written in terms of vectors $g(\mathcal{L})f$, for some filter functions $g$ which act on the spectrum of $\mathcal{L}$ (scaling, selecting, …). Fourier analysis on graphs has been successfully applied to many different fields such as stationary signals on graphs [41], graph signal energy study [26], convolutional neural networks on graphs [18].

Graph wavelets are an important application of graph Fourier analysis, and several definitions of graph wavelets were proposed [15, 13, 24, 34, 57, 27]. When performing Fourier analysis of a signal, there is no guarantee that localization of a signal in the frequency domain (a.k.a Fourier domain) implies localization in the graph domain. This phenomenon is illustrated by the fact that the eigenvectors corresponding to the upper part of Laplacian spectrum tend to be more oscillating than those from the bottom of the spectrum (see for example [58, Fig. 1.6, p. 28] for an illustration). To overcome this problem, [28] developed a fairly general construction of a frame enjoying the usual properties of standard wavelets: each vector of the frame is defined as a function $g(s\mathcal{L})\delta_{m}$ (where $\delta_{m}$ is a Kronecker signal, having zero values at every vertex except $m$) and is localized both in the graph domain and the spectral domain at fine scale $s$. The transform associated with this frame is named Spectral Graph Wavelet Transform (SGWT), and it was used in numerous subsequent works [54, 2, 27].

Signals which are sparse in the Fourier domain form an important class of graph signals. Indeed, there is a tight relationship between sparsity in the Fourier domain and the notion of regularity of a signal $f$ on the vertices of a graph $\mathcal{G}$ which comes from the Laplacian matrix $\mathcal{L}$ of $\mathcal{G}$. Intuitively, a smooth signal will not vary much between two vertices that are close in the graph. This regularity property can be read in the Fourier domain: a very smooth signal will be correctly represented in the Fourier domain with a small number of eigenvectors associated with the lower spectral values; on the contrary, non-smooth signals (i.e. highly oscillating) are represented with eigenvectors corresponding to the upper part of the spectrum. Both the types of signal are said frequency sparse.

In this paper, we propose to exploit localization in the Fourier domain to improve graph Fourier analysis. More precisely, we consider vectors of the form $g(\mathcal{L}_{I_{k}})f$ instead of vectors $g(\mathcal{L})f$ in graph Fourier analysis, where $\mathcal{L}_{I_{k}}$ is defined as the matrix $\mathcal{L}P_{I_{k}}(\mathcal{L})$ and $P_{I_{k}}(\mathcal{L})$ denotes the projection onto the eigenspaces whose eigenvalue is contained in subset $I_{k}$. Localized Fourier analysis is motivated by problems and properties defined on strict subsets of the spectrum ${\rm sp}(\mathcal{L})$ (e.g. any problem defined in terms of frequency sparse graph signals). As a central application of Fourier localization, we introduce the Fourier localized counterpart of SGWT, that we call LocLets for Localized graph wavelets. We prove that various frame constructions can be written in terms of LocLets, hence benefiting from all the advantages of localization discussed in this paper.

Defining $I_{\mathcal{L}}$ as the smallest interval containing the entire spectrum ${\rm sp}(\mathcal{L})$, the local Fourier analysis consists in choosing a suitable partition $I_{\mathcal{L}}=\sqcup_{k}I_{k}$ into subintervals on which standard Fourier analysis is performed. Such an analysis on disjoint intervals naturally benefits from several interesting properties. In particular, when $f$ is modeled by a Gaussian random vector with independent entries, the disjointness of subintervals preserves these properties in the sense that random variables $(g(\mathcal{L}_{I_{k}})f)_{k}$ are still Gaussian and independent. This simple observation has important consequences to study the graph problem at stake. In this work, it allows us to propose some noise level estimator from the random variables sequence $(g(\mathcal{L}_{I_{k}})f)_{k}$, and to provide a theoretical analysis of the denoising problem. Disjointness of subsets $(I_{k})_{k}$ also provides simple strategies to parallelize Fourier analysis computations.

We also consider the general problem given by a noisy signal on a graph $\widetilde{f}=f+\xi$, where $\xi$ is some random Gaussian vector with noise level $\sigma$. We provide results for two important tasks: the estimation of $\sigma$ when the latter is unknown, and the denoising of noisy signal $\widetilde{f}$ in order to recover signal $f$. We show that for frequency sparse signals, localization allows to adapt to the unknown Fourier support of signal $f$. Theoretical guarantees and practical experiments show that localized Fourier analysis can improve state-of-the-art denoising techniques, not only in precision of the estimator $\widehat{f}$ of $f$, but also in time computations.

We provide an efficient method to choose a partition $I_{\mathcal{L}}=\sqcup_{k}I_{k}$ for the Fourier localized vectors $g(\mathcal{L}_{I_{k}})f$ to be sufficiently informative. Using well-known techniques for efficient graph Fourier analysis (a.k.a Chebyshev filter approximations), we propose scalable methods to perform localized Fourier analysis with no computational overhead over standard fast Fourier graph analysis. In particular, all methods introduced in the paper avoid the computation of the entire eigendecomposition of $\mathcal{L}$ which is a major computational advantage when considering large graphs in applications.

The paper is structured as follows. Section 2 presents the relevant notions and techniques necessary to perform localized Fourier analysis. In Section 2.2, we introduce LocLets, the Fourier localized extension of SGWT. Section 3 is devoted to the study of the denoising problem for signals on graphs. The section provides results about noise level estimation, and signal denoising. Additional properties of LocLets, such as computational aspects and relationships with known wavelet transforms, are further developed in Section 4. In Section 5, we analyze the experiments made to support the interesting properties of localized Fourier analysis highlighted in this paper. Finally, the proofs are gathered in Section 7.

2 Localized Fourier analysis for graph signals

In this section, we introduce the central notion studied in this paper: localization of graph Fourier analysis. First, we recall the relevant notions of graph Fourier analysis. Then we provide examples from previous works that motivates localization in the Fourier domain. We also introduce LocLets, an important application of Fourier localization to SGWT. Finally, we discuss briefly the particular case of graphs sampled from a manifold.

2.1 Functional calculus and Fourier analysis for graph signals

Let $\mathcal{G}=(\mathcal{V},\mathcal{E})$ be an undirected weighted graph with $\mathcal{V}$ the set of vertices, $\mathcal{E}$ the set of edges, $n=|\mathcal{V}|$ the number of nodes (the size of $\mathcal{G}$), and $(W_{ij})_{i,j\leq n}$ the weights on edges. Let us introduce the diagonal degree matrix whose diagonal coefficients are given by $D_{ii}=\sum_{1\leq j\leq n}W_{ij}$ for $1\leq i\leq n$. The resulting non-normalized Laplacian matrix $\mathcal{L}$ of graph $\mathcal{G}$ is defined as $\mathcal{L}=D-W$. The $n$ non-negative eigenvalues of $\mathcal{L}$, counted without multiplicity, are denoted by $\lambda_{1},\ldots,\lambda_{n}$ in the decreasing order. In the sequel, $\mathrm{sp}(\mathcal{L})$ stands for the spectrum of $\mathcal{L}$. The corresponding eigenvectors are denoted $\chi_{1},\ldots,\chi_{n}$.

Given a graph $\mathcal{G}$, the GFT of a real-valued function $f$ defined on the vertices of $\mathcal{G}$ is nothing but the representation of $f$ in the orthonormal basis of eigenvectors of $\mathcal{L}$. Namely, for a signal $f:\mathcal{G}\rightarrow\mathbb{R}$, the $\ell$-th Fourier coefficient of $f$, denoted $\widehat{f}(\ell)$, is given by $\widehat{f}(\ell)=\langle f,\chi_{\ell}\rangle$. The Fourier support ${\rm supp}(\widehat{f})$ of signal $f$ is the set of indices $\ell$ such that $\widehat{f}(\ell)\neq 0$. We will see in Section 2.2 that graph wavelets can be defined in a similar manner.

Functional calculus is a powerful technique to study matrices, and constitutes the heart of GSP. For a function $g$ defined on some domain $D_{g}$, $\mathrm{sp}(\mathcal{L})\subset D_{g}$, functional calculus reads as

[TABLE]

Interpreting the eigenvalues $\lambda_{\ell}$, $\ell=1,\ldots,n$, as the fundamental frequencies associated with a graph, the linear map $g(\mathcal{L})$ is generally seen as a filter operator in terms of signal analysis.

Also, spectral projections of matrix $\mathcal{L}$ can be made explicit with the help of functional calculus, by setting $g=\mathrm{1}_{I}$. More precisely, for any subset $I\subset I_{\mathcal{L}}$, consider the map $P_{I}(\mathcal{L})$ given by:

[TABLE]

Then, $P_{I}(\mathcal{L})$ is nothing but the spectral projection on the linear subspace spanned by the eigevectors associated with the eigenvalues belonging to $I$. In the sequel, $n_{I}=|I\cap\mathrm{sp}(\mathcal{L})|$ will stand for the number of eigenvalues contained in subset $I\cap\mathrm{sp}(\mathcal{L})$.

Spectral projections are a practical tool to focus on some part of the spectrum $\mathrm{sp}(\mathcal{L})$. More precisely, let $I_{\mathcal{L}}=\sqcup_{1\leq k\leq K}I_{k}$ be a partition of interval $I_{\mathcal{L}}$ into disjoint subsets $(I_{k})_{k}$. Since intervals $I_{k}$ are disjoints, functional analysis of $\mathcal{L}$ reduces to that of its projections $\mathcal{L}_{I_{k}}=\mathcal{L}P_{I_{k}}(\mathcal{L})$ in the sense of the identity:

[TABLE]

In this paper, one will study the extent to which Fourier analysis on large graphs is improved when considering local Fourier analysis on each subset $I_{k}$ instead of global Fourier analysis on $I_{\mathcal{L}}$.

2.2 LocLets: a localized version of SGWT

This section introduces an important application of the localized graph Fourier analysis, namely the notion of localized SGWT.

2.2.1 Construction of a SGWT

Let $f:\mathcal{G}\rightarrow\mathbb{R}$ be a signal on the graph $\mathcal{G}$. Let $\varphi,\psi:\mathbb{R}\rightarrow\mathbb{R}$ be respectively the scaling and kernel functions (a.k.a. father and mother wavelet functions), and let $s_{j}>0$, $1\leq j\leq J$, be some scale values. The discrete SGWT is defined in [28] as follows:

[TABLE]

The adjoint matrix $\mathcal{W}^{*}$ of $\mathcal{W}$ is:

[TABLE]

We also recall from [28] that a discrete transform reconstruction formula using SGWT coefficients $(c_{j,m})_{\begin{subarray}{c}0\leq j\leq J\\ 1\leq m\leq n\end{subarray}}$ is obtained by the formula

[TABLE]

where $(\mathcal{W}^{*}\mathcal{W})^{-1}$ stands for a pseudo-inverse of the matrix $\mathcal{W}^{*}\mathcal{W}$.

Remark 1**.**

From a theoretical point of view, no more assumptions on the scale values $s_{j}$ are required. However, the choice in practice of the scale values depends simultaneously on the mother wavelet $\psi$ and the maximal spectral value of the Laplacian $\mathcal{L}$. In [28], these values $s_{j}$ are suitably chosen accordingly with the graphs considered in the experiments. Besides, for our experiments, we follow the construction of [27] in which $s_{j}=b^{j}$ for some real $b>1$ (see also Section 4.2.1 for the details).

2.2.2 Motivation for considering Fourier localized graph signals

The definition of SGWT as given in [28] is closely related to its counterpart from traditional wavelet transform by the action of the transform on the Fourier domain. Equation (9) in [28] highlights the following decomposition of wavelet coefficients in the frequency domain:

[TABLE]

where $\hat{\psi}^{*}(s\omega)=\hat{\overline{\psi}_{s}}(\omega)$ denotes the Fourier coefficient of a scaled version $\overline{\psi}_{s}$ of the mother wavelet $\psi$, and $x\mapsto e^{i\omega x}$ are the eigenfunctions of the one-dimensional Laplacian $\mathcal{L}=\frac{\partial^{2}}{\partial x^{2}}$ on the real line. On the other hand, a similar formula holds for the case of graph wavelets where integration over frequencies is replaced by summation over Laplacian eigenvalues. It is given by Equation (22) in [28]:

[TABLE]

Equations 2 and 3 suggest strong similarities between the role of frequencies for wavelets on the real line and the role of Laplacian eigenvalues for graph wavelets. For example, an explicit mapping between frequencies and Laplacian eigenvalues is exhibited in part III.C of [34] for the case of the cycle graph. As an important consequence, filtering in the frequency domain for traditional time signals is similar to filtering Laplacian spectrum in the graph setting. This similarity has guided authors of [33] who propose an analog of Meyer wavelets [39] for graphs by applying traditional Meyer filters on the graph Fourier domain. As for the traditional case, smoothness for graph Meyer wavelets can be guaranteed by calibrating Meyer filters so that the Fourier support of wavelet functions is concentrated in the bottom of the Laplacian spectrum. A larger list of tight frames candidates obtained from traditional wavelets defined on the Fourier domain is provided in Table 1 of [34].

Other examples of signals with limited frequency support have inspired graph counterparts with restricted Fourier supports. In [1, 8, 42, 43], $k$-bandlimited graph signals are supported on the eigenspaces associated with the $k$ smallest eigenvalues and define an equivalent of $\omega$-bandlimited signals. Compressive sampling is an application of this type of signals both in the traditional [7, 20, 6] and in the graph [42] settings. While low-pass frequency filters are usually used to capture the smooth part of a signal, high-pass filters have the ability to point out anomalies in smooth signals. Such applications of high-pass filters can be found in [9] for optical laser measurements and time signals, but also in [45] for anomalies detection in temperature measurements on a graph of weather stations. The latter application was also considered in Section V.D of [8] to illustrate the use of filter banks that split the graph signal into two bandlimited signals. Authors of [55, 56] introduce graph sampling methods in the Fourier domain taking advantage of periodic patterns and extending to the graph setting sampling techniques designed for shift-invariant signals [21]. The sampling task for time signals was addressed in [30, 36] in the case of multiband signals supported on a disjoint union of frequency subsets. This suggests providing methods adapted to disjoint subsets in the Fourier domain also for the case of multiband graph signals.

All the examples discussed above share a common feature: their Fourier support is localized in the Fourier domain in the sense that such signals $f$ satisfy $f=P_{I}(\mathcal{L})f$ for some subset $I\subset I_{\mathcal{L}}$, $I\cap\rm{sp}(\mathcal{L})\neq\rm{sp}(\mathcal{L})$. For $k$-bandlimited signals, $I$ is a set containing only the $k$ smallest eigenvalues of $\mathcal{L}$; for high-pass filters, $I$ contains only the largest eigenvalues; for periodic signals, $I$ is a set such that $I\cap\rm{sp}(\mathcal{L})=\{\lambda_{i~{}\rm{mod}~{}m}~{}|~{}\lambda_{i}\in\rm{sp}(\mathcal{L})\}$ for some period $m$; and for multiband signals, $I$ is a disjoint union $I=\sqcup_{k}I_{k}$ of intervals $(I_{k})_{k}$. The current paper proposes methods to study these graph signals that we call Fourier localized signals.

2.2.3 Definition of LocLets

Spectral graph wavelet functions are given by $(\varphi(\mathcal{L})\delta_{m},\psi(s_{j}\mathcal{L})\delta_{m})_{1\leq j\leq J,1\leq m\leq n}$. We define a LocLet function to be the projection of a graph wavelet function onto a subset of eigenspaces of $\mathcal{L}$.

Definition 2**.**

Let $(\varphi(\mathcal{L})\delta_{m},\psi(s_{j}\mathcal{L})\delta_{m})_{1\leq j\leq J,1\leq m\leq n}$ be the family functions induced by a SGWT. Then, for any subset $I\subset I_{\mathcal{L}}$, $1\leq j\leq J$ and $1\leq m\leq n$, set:

[TABLE]

where $\widehat{\delta_{m}}(\cdot)=\langle\delta_{m},\chi_{\cdot}\rangle$ is the graph Fourier transform of $\delta_{m}$. The functions $(\varphi_{m,I},\psi_{j,m,I})_{1\leq j\leq J,1\leq m\leq n}$ are called Localized waveLets functions (LocLets). The functions $\varphi_{m,I},\psi_{j,m,I}$ are said to be localized at $I$.

Let $I_{\mathcal{L}}=\sqcup_{1\leq k\leq K}I_{k}$ be some partition. Then, the localized SGWT transfom of $f$ with respect to partition $(I_{k})_{1\leq k\leq K}$, denoted by $\mathcal{W}^{(I_{k})_{k}}f$, is defined as the family $\mathcal{W}^{(I_{k})_{k}}f=(\mathcal{W}^{I_{k}}f)_{k}$ where

[TABLE]

Similarly to Equation (1), the adjoint transform is given by

[TABLE]

As already observed, localized SGWT of a signal $f$ contains more precise information about signal $f$ than its standard SGWT. The latter can easily be obtained from the former since subsets $(I_{k})_{k}$ are pairwise disjoint and formula $g(s\mathcal{L})=\sum_{1\leq k\leq K}g(s\mathcal{L}_{I_{k}})$ holds for all filter $g$, and in particular for $g=\varphi$ or $g=\psi$. When the partition $I_{\mathcal{L}}=\sqcup_{1\leq k\leq K}I_{k}$ is carefully chosen, we show that the SGWT localization provides interesting features such as independence of random variables in denoising modelling, or considerable improvements in denoising tasks.

Remark 3**.**

The functions $\varphi_{m,I},\psi_{j,m,I}$ are localized in the Fourier domain in the sense that the support of their Fourier transforms are contained in a subset $I$ of $I_{\mathcal{L}}$. A different localization property, observable in the graph domain, is considered in [28]. We refer to the latter property as graph domain localization in the current paper.

A property about the graph domain localization at fine scales is stated in [28, Theorem 5.5] but this result appears to be not informative in general for the case of Fourier localized functions $\varphi_{m,I},\psi_{j,m,I}$. For instance, any function of the form $\psi 1_{I}$ vanishes in a neighborhood of [math] as soon as $0\notin I$, as observed in Section 2.3 of [27]. However, other graph domain localization results were obtained by the authors of [14] for the case of frames considered in [27] and discussed in our Section 4.2.1. In addition, the weaker graph localization property [28, Lemma 5.2] for powers of the Laplacian still holds in our Fourier localized setting since in practice we approximate any function $x\to\psi(s_{j}x)1_{I}(x)$ by a Chebyshev polynomial of order $N$ (see Section 4).

2.3 Fourier localization for weighted graphs sampled from manifolds

In some applications, the underlying graph is unknown and is built from the data. In this case, the resulting graph is thought as a random sampling of a low-dimensional sub-manifold embedded in a higher dimensional Euclidean space.

More precisely, let $\mathcal{M}$ be a Riemannian manifold of dimension $d$ embedded in $\mathbb{R}^{m}$ with $m>d$. A popular way to define a graph from a finite set of points $\{x_{1},\ldots,x_{n}\}\subset\mathcal{M}$ consists in defining a weighted adjacency matrix $W=(W_{ij})_{i,j\leq n}$ as follows:

[TABLE]

where $\|\cdot\|_{2}^{2}$ stands for the Euclidean distance in $\mathbb{R}^{m}$ and $\varepsilon>0$ is some parameter called the bandwidth of the kernel $k$. A typical choice for the kernel $k$ is the exponential function $k(x)=\exp(-x)$, $x\in\mathbb{R}$. As an example, the swissroll graph of Section 5 is built following this idea.

A whole part of the literature is dedicated to the question of the convergence of the discrete (normalized or non-normalized) Laplacian matrices $L_{n,\varepsilon}$ toward the Laplace-Beltrami operator $\Delta_{\mathcal{M}}$ (see [44] for a detailed exposition of this classical object from differentiable geometry). The discretized operator $L_{n,\varepsilon}$ depending on two parameters, the convergence as $n\to\infty$ and/or $\varepsilon\to 0$ have been considered in [29, 53, 12, 25, 3, 38]. Loosely speaking, theses results are devoted to the approximation of $\Delta_{\mathcal{M}}f$ by $L_{n,\varepsilon}f$ at the sample points.

Furthermore, it is shown in [59, 23] that, under mild conditions, eigenvalues and eigenfunctions of $\Delta_{\mathcal{M}}$ are well approximated by those of $L_{n,\varepsilon}$. As a consequence of particular interest, bandlimited (or even multiband) signals on the manifold $\mathcal{M}$, when sampled at the points $\{x_{1},\ldots,x_{n}\}\subset\mathcal{M}$, are expected to be bandlimited (multiband) with respect to the graph Laplacian $L_{n,\varepsilon}$ (with slight differences when considering the normalized or the non-normalized Laplacian). It is worth noting that the sample points $\{x_{1},\ldots,x_{n}\}$ do not have to lie exactly in the manifold $\mathcal{M}$ but can be disrupted by a noise. The spectral properties are preserved by standard spectral perturbation arguments (see [32]) providing the noise level is sufficiently small. Such a perturbation argument is discussed at some point in [12] and remains valid in our context.

To conclude this discussion, let us point out that the choice of a Gaussian kernel in (4), while popular, is quite arbitrary. The results in [23] are stated for a rather general kernel including non-smooth kernels. In addition, a variable bandwidth kernel is also considered in [5].

3 Local Fourier analysis and graph functions denoising

The denoising problem is stated as follows: given an observed noisy signal $\widetilde{f}$ of the form $\widetilde{f}=f+\xi$ where $\xi$ is a $n$-dimensional Gaussian vector distributed as $\mathcal{N}(0,\sigma^{2}\mathrm{Id})$, provide an estimator of the a priori unknown signal $f$.

This section shows how localized Fourier analysis helps in estimating the noise level $\sigma$ when it is unknown, and in recovering the original signal $f$ when the latter is sparse in the Fourier domain. In what follows, we will focus on random variables of the form $\|P_{I_{k}}\widetilde{f}_{I_{k}}\|_{2}$ where $\widetilde{f}$ is the noisy signal and $I_{k}$ is a subset in the partition $I_{\mathcal{L}}=\sqcup_{k}I_{k}$. To keep the notations light, $n_{k}$, $f_{k}$, $\xi_{k}$ and $\widetilde{f}_{k}$ will stand for $n_{I_{k}}$, $P_{I_{k}}f$, $P_{I_{k}}\xi$ and $P_{I_{k}}\widetilde{f}$ respectively. In addition, the cumulative distribution function of a random variable $X$ will be denoted by $\Phi_{X}$.

3.1 Noise level estimation for frequency sparse signals

Since in real application the noise level $\sigma$ remains unknown in general, new estimators $\widehat{\sigma}$ based on localization properties in the spectrum are introduced in the sequel.

3.1.1 Noise level estimation from projections along ${\rm sp}(\mathcal{L})$

For any filter $g$ defined on $I_{\mathcal{L}}$ and any subset $I\subset I_{\mathcal{L}}$, simple computations give rise to

[TABLE]

Since both $\widetilde{f}^{T}g(\mathcal{L}_{I})\widetilde{f}$ and ${\rm Tr}(g(\mathcal{L}_{I}))$ are known, Equation (5) suggests building estimators from the expression $\frac{\widetilde{f}^{T}g(\mathcal{L}_{I})\widetilde{f}}{{\rm Tr}(g(\mathcal{L}_{I}))}$. In [17], the noise level is estimated by $\frac{\widetilde{f}^{T}\mathcal{L}\widetilde{f}}{{\rm Tr}(\mathcal{L})}$ which can be seen as the graph analog of the Von Neumann estimator from [60]. The main drawback of this estimator is its bias.

Theoretically, without any assumption on the signal $f$, the bias term $\frac{f^{T}g(\mathcal{L}_{I})f}{{\rm Tr}(g(\mathcal{L}_{I}))}$ is minimized when $g=\mathrm{1}_{\{\lambda_{\ell^{\ast}}\}}$ where $\ell^{\ast}=\mathrm{argmin}\{|\widehat{f}(\ell)|:\lambda_{\ell}\in\mathrm{sp}(\mathcal{L})\}$. The computation of such filters would require the complete reduction of $\mathcal{L}$ which does not scale well with the size of the graph. Instead, these ideal filters will be approximated by filters of the form $g=\mathrm{1}_{I_{k}}$, for $I_{k}$ a subset in the partition $I_{\mathcal{L}}=\sqcup_{k}I_{k}$. It is worth noting that with $k^{\ast}={\rm argmin}_{k}\|f_{k}\|_{2}$, the function $g^{\ast}=\mathrm{1}_{I_{k^{\ast}}}$ achieves the minimal bias of the estimator among all filters of the form $g=\sum_{k}\alpha_{k}\mathrm{1}_{I_{k}}$.

Discarding some intervals $I_{k}$ with $n_{k}=0$, it can be assumed without loss of generality that $n_{k}\neq 0$ for all $1\leq k\leq K$. Also, observe that the random variable $\|\widetilde{f}_{k}\|_{2}^{2}$ can be decomposed as follows

[TABLE]

where $\frac{\|\xi_{k}\|_{2}^{2}}{\sigma^{2}}$ and $\frac{\langle f_{k},\xi_{k}\rangle}{\sigma}$ are random variables distributed as $\chi^{2}(n_{k})$ and $\mathcal{N}(0,\|f_{k}\|_{2}^{2})$ respectively.

Proposition 4**.**

Let $(c_{k})_{1\leq k\leq K}$ be the sequence of non-negative random variables defined, for all $k=1,\ldots,K$, by $c_{k}=\|\widetilde{f}_{k}\|_{2}^{2}/n_{k}$. Then,

the random variables $c_{1},\ldots,c_{K}$ are independent; 2. 2.

for all $k,k^{\prime}$ such that $f_{k}=f_{k^{\prime}}$, $c_{k}$ and $c_{k^{\prime}}$ are identically distributed if and only if $n_{k}=n_{k^{\prime}}$; 3. 3.

for $k$ such that $f_{k}=0$, $c_{k}$ is distributed as $\frac{\sigma^{2}}{n_{k}}\Gamma_{n_{k}}$ where $\Gamma_{n_{k}}\sim\chi^{2}(n_{k})$.

3.1.2 The case of frequency-sparse signals

When the signal $f$ is sparse in the Fourier domain, the condition $f_{k}=0$ is met for most of the intervals $I_{k}\subset I_{\mathcal{L}}$. Let us define $I_{f}=\sqcup_{k:I_{k}\cap{\rm supp}\widehat{f}\neq\emptyset}I_{k}$ to be the union of subsets $I_{k}$ intersecting the Fourier support ${\rm supp}(\widehat{f})$ of $f$. Also, denote by $\overline{I_{f}}=I_{L}\backslash I_{f}$ its complement set. In order to take advantage of Fourier sparsity, let us introduce the quantities $\widehat{\sigma}_{\rm mean}$ and $\widehat{\sigma}_{\rm med}$ as follows:

[TABLE]

The following concentration inequalities show that $\widehat{\sigma}_{\rm mean}$ and $\widehat{\sigma}_{\rm med}$ are natural estimators of the noise level $\sigma$.

Proposition 5**.**

Let $K_{f}=|\{k:I_{k}\subset\overline{I_{f}}\}|$, $n_{0}=\min\{n_{k}:k,I_{k}\subset\overline{I_{f}}\}$, $n_{\infty}=\max\{n_{k}:k,I_{k}\subset\overline{I_{f}}\}$, $V_{f}=2\sigma^{4}\sum_{k:I_{k}\subset\overline{I_{f}}}{1}/{n_{k}}$ and $B_{f}={2\sigma^{2}}/n_{0}$. Then the following concentration inequalities hold:

for all $t\geq 0$,

[TABLE]

and for all $0\leq t\leq\sigma^{2}$,

[TABLE] 2. 2.

for all $t\geq 0$, with $\beta=n_{0}/n_{\infty}$,

[TABLE]

and for all $0\leq t\leq 1$ such that $p^{-}(t)\leq 1/2$,

[TABLE]

where

[TABLE]

Obviously, the Fourier support ${\rm supp}(\widehat{f})$ and the subset $\overline{I_{f}}$ remain generally unknown in applications and have to be approximated. Let us recall that the main issue for estimating $\sigma$ comes from the bias term $\frac{\|f_{k}\|_{2}^{2}}{n_{k}}$ in Equation (5), and in particular when the value $\sigma^{2}$ is negligible compared to $\frac{\|f_{k}\|_{2}^{2}}{n_{k}}$. Therefore, a suitable candidate to approximate $\overline{I_{f}}$ will be some subset $\overline{J_{f}}\subset I_{\mathcal{L}}$ for which the impact of larger values $\frac{\|f_{k}\|_{2}^{2}}{n_{k}}$ is minimized. This is made clear by Proposition 6 below. The latter involves the following concentration bounds for Gaussian random variables: for all $0<\alpha<1$

[TABLE]

Proposition 6**.**

Let $0<\alpha<1$. Let $t_{\alpha,\sigma}$ be defined by Equation (8). Assume that $f_{\ell}=0$ and that the following inequality holds:

[TABLE]

Then, the quantities

[TABLE]

satisfy $\mathbb{P}(b_{k}\geq b_{\ell})\geq 1-\alpha$.

By invariance under permutations, one may assume without loss of generality that the values $c_{k}$ are ordered in the decreasing order. Proposition 6 quantifies the fact that the highest values of $c_{k}$ correspond most likely to the indices $k$ for which $f_{k}\neq 0$. Consequently, setting $\overline{J_{f}}(r)=\sqcup_{k\in\{r,r+1,...K-r\}}I_{k}$ for all $1\leq r\leq\frac{K}{2}$, the estimators introduced in Equation (7) may be rewritten replacing the unknown subset $\overline{I_{f}}$ by its known approximation $\overline{J_{f}}(r)$. So we define the estimators

[TABLE]

It is worth noting that from the symmetry of the subset $\overline{J}_{f}(r)$, it follows that the value $\widehat{\sigma}_{\rm med}^{r}$ actually does not depend on parameter $r$, and one will write $\widehat{\sigma}_{\rm med}$ in place of $\widehat{\sigma}_{\rm med}^{r}$.

3.2 Denoising Frequency Sparse Signals

Let us begin with a result illustrating that localized Fourier analysis in $I_{\mathcal{L}}$ provides strong benefits in noise reduction tasks when the underlying signal is frequency sparse.

Proposition 7**.**

Assume $f=f_{I}$ for some subset $I\subset I_{\mathcal{L}}$. Then

[TABLE]

In particular, denoising of $\widetilde{f}$ boils down to denoising of $\widetilde{f}_{I}=f_{I}+\xi_{I}$.

While Proposition 7 asserts a trivial denoising solution in the Fourier domain, i.e. simply destroying the projection $\widetilde{f}_{\overline{I}}=\xi_{\overline{I}}$, this approach is no longer that immediate when considering the graph domain observations since the Fourier support of $f$ is unknown in practice and needs to be estimated. Based on the $\chi^{2}$-statistics, Algorithm 1 is designed for this purpose. To the best of our knowledge, previous works that proposed method for Fourier support recovery for graph noisy signals [47] involve the complete eigendecomposition of matrix $\mathcal{L}$. The methodology suggested below makes use of projectors on eigenspaces which can be approximating with Chebyshev polynomials as detailed in the next Section 4.

Heuristically, if $I$ contains the support of the Fourier transform of $f$, on the complementary subset $\overline{I}$ we only observe pure white Gaussian noise so that $\|P_{\overline{I}}\widetilde{f}\|_{2}^{2}=\|\widetilde{f}_{\overline{I}}\|_{2}^{2}$ is distributed as $\sigma^{2}\chi^{2}(n_{I})$ with $n_{I}=|\overline{I}\cap\mathrm{sp}(L)|$. On the other hand, on $I$ the square of the Euclidean norm of a non-centered Gaussian vector is observed. Consequently, the quantity $\mathbb{P}\left(\chi^{2}(n_{I})>\sigma^{-2}\|P_{I}\widetilde{f}\|^{2}_{2}\right)$ is typically very close to zero whereas $\mathbb{P}\left(\chi^{2}(n-n_{I})>\sigma^{-2}\|P_{\overline{I}}\widetilde{f}\|^{2}_{2}\right)$ remains away from [math]. To put it in a nutshell, sliding a window along the spectrum of $\mathcal{L}$, Algorithm 1 performs a series of $\chi^{2}$-test.

With the objective to provide theoretical guarantees that $\chi^{2}$-tests approach ${\rm supp}(\widehat{f})$ correctly, it is important to turn the condition on the $p_{k}$-value into a condition involving only the values $\|f_{k}\|_{2}$ and $\sigma$. The next lemma shows that for sufficiently large values of the ratio $\frac{\|f_{k}\|_{2}}{\sigma}$, the inequality $p_{k}\leq\alpha$ holds so that the corresponding components ${\rm supp}(\widehat{f}_{k})$ of the Fourier domain are legitimately included in the support estimate $I$.

Lemma 8**.**

Let $0<\alpha<1$ and let $\Gamma_{n_{k}},\Gamma_{n_{k}}^{\prime}$ be two i.i.d $\chi^{2}(n_{k})$ random variables. Assume that:

[TABLE]

where $t_{\alpha,\sigma}$ is defined by Equation (8). Then $p_{k}\leq\alpha$.

In contrast to Lemma 8, the lemma below states that condition $p_{k}>\alpha$ holds for sufficiently small values of ratio $\sigma^{-1}\|f_{k}\|_{2}$.

Lemma 9**.**

Let $0<\alpha<1$ and let $\Gamma_{n_{k}}$ be a $\chi^{2}(n_{k})$ random variable. For $0<\beta<1$, set $t_{\beta,k}=\sigma^{2}\Phi_{\Gamma_{n_{k}}}^{-1}(1-\beta)$. Assume that

[TABLE]

Then $p_{k}>\alpha$.

Compared to Proposition 7, the result below quantifies the error resulting by approximating the support running Algorithm 1. Note that the requirement to have a constant sequence $(n_{k})_{k}$ is used for statement clarity but similar assertions hold for the case $n_{k}\neq n_{k^{\prime}}$.

Proposition 10**.**

Set $f_{I}=\sum_{k:p_{k}\leq\alpha}P_{I_{k}}f$. Assume that $n_{k}=n_{1}$ for all $1\leq k\leq K$. Then,

the Fourier support approximation $\ell_{2}$-error satisfies

[TABLE] 2. 2.

the Noise $\ell_{2}$-error on Fourier support:

[TABLE]

Lemma 9 asserts that the set $\{k,p_{k}>\alpha\}$ is small when most of the values $\|f_{k}\|_{2}$ are large enough compared to noise level $\sigma$ for $I_{k}\cap{\rm supp}\widehat{f}\neq\emptyset$. In such a case, Fourier support approximation $\ell_{2}$-error is small. Regarding the noise $\ell_{2}$-error, the inclusion $\{k,p_{k}\leq\alpha\}\subset\{k,I_{k}\cap{\rm supp}\widehat{f}\neq\emptyset\}$ holds by Lemma 9. Moreover, Lemma 8 asserts that the set $\{k,p_{k}\leq\alpha\}$ contains the entire set $\{k,I_{k}\cap{\rm supp}\widehat{f}\neq\emptyset\}$ for sufficiently large values of $\sigma^{-1}\|f_{k}\|_{2}$ when $I_{k}\cap{\rm supp}\widehat{f}\neq\emptyset$. For such favorable situations, the noise $\ell_{2}$-error is exactly $n_{1}\sigma|\{k,I_{k}\subset I_{f}\}|$, the amount of noise on the extended support $I_{f}$.

The second step gives an estimate of the original signal using a thresholding procedure on each element $\widetilde{f}_{I}$ and $\widetilde{f}_{\overline{I}}$. On the one hand, the methodology developed in [27] is prohibitive in terms of time and space complexity as soon as the underlying graphs become moderately large. On the other hand, the fast SGWT remains an approximating procedure. If a signal happens to be very frequency-sparse, then an even more optimal strategy is possible: first, the support $I$ in the frequency domain is approximated with the help of Algorithm 1; then, the procedure of [27] is applied to $P_{I}f$ (the low-rank part) and LocLets on $P_{\overline{I}}(\mathcal{L})f$. This idea is made precisely in Algorithm 3.

Estimator $\widehat{f}$ produced in Algorithm 3 satisfies a tighter oracle bound inequality than the one given in [27, Theorem 3]. This theoretical guarantee is widely supported by our experiments described in Section 5. Following notations from [27, Equation (21)], we denote by $OB(f_{I})$ the oracle bound obtained from an oracle estimator of $f_{I}$ from a noisy $\widetilde{f}_{I}$ exploiting some knowledge about the unknown signal $f_{I}$. We refer to [27] for precise details.

Theorem 11**.**

Let $I,\widehat{f}$ be respectively the support approximation and the estimator of $f$ obtained from Algorithms 1 and 3 with threshold value $t_{2}=0$. Then we have

[TABLE]

The right-hand side in the inequality of Theorem 11 has a more explicit expression in terms of $\alpha,\sigma$ using Proposition 10. Up to the error made by approximating the support with Algorithm 1, the $\ell_{2}$-risk is essentially bounded by the $\ell_{2}$-risk of the Parseval frame procedure from [27] on the low-rank projection $f_{I}$ of $f$, that is

[TABLE]

To conclude, Theorem 11 provides a theoretical guarantee that the support approximation improves the denoising performances obtained from [27].

4 Properties of LocLets

In this section, we highlight important properties for the application of Fourier localization in practice. First we discuss computational analysis, and methods to apply our techniques to large graphs. Then we study the relationships of LocLets with well-known graph wavelet constructions.

4.1 Fast LocLet Transform and Computational Analysis

In the case of large graphs, GSP requires a special care for being efficient since functional calculus relies a priori on the complete reduction of the Laplacian. Actually, several efficient methods were designed to retrieve only partial information from the eigendecomposition as matrix reduction techniques (see for instance [37, 54]) or polynomial approximations [28, 50, 19]. In this paper, the widely adopted latter approach with Chebyshev polynomials approximation is preferred and briefly recalled below (we refer the reader to [50, Section III.C.] for a brief but more detailed description of Chebyshev approximation).

4.1.1 Chebyshev approximations

Roughly speaking, the idea is to approximate the function $g$ with its Chebyshev expansion $g_{N}$ at order $N$. More precisely, the Chebyshev polynomials of the first kind $(T_{i})_{i\geq 0}$ are defined from the second order recursion

[TABLE]

for $i\geq 2$. Then, the matrix $\mathcal{L}$ is normalized as $\widetilde{\mathcal{L}}=\frac{2}{\lambda_{1}}\mathcal{L}-I_{n}$ so that $\mathrm{sp}(\widetilde{\mathcal{L}})\subset[-1,1]$. This gives rise to some function $\widetilde{g}:[-1,1]\rightarrow\mathbb{R}$ with the property $g(\mathcal{L})=\widetilde{g}(\widetilde{\mathcal{L}})$. In fact, $\widetilde{g}(x)=g(\frac{\lambda_{1}}{2}(x+1))$ for all $x\in[-1,1]$. Then $g(\mathcal{L})$ has the following truncated Chebyshev expansion $g(\mathcal{L})\approx g_{N}(\mathcal{L})$:

[TABLE]

where $N$ is the maximal degree of polynomials $T_{i}$ used in the expansion, and $a_{i}(\widetilde{g})$ is the $i$-th coefficient in the $N$-th order Chebyshev expansion of function $\widetilde{g}$. Following [28], for any filter $g$ on $\mathrm{sp}(\mathcal{L})$ and any signal $f$ on graph $\mathcal{G}$, the approximation $g_{N}(\mathcal{L})$ provides a vector value close to $g(\mathcal{L})f$ with time complexity $O(|\mathcal{E}|N)$.

The object presented in the sequel involves in particular the spectral projection $P_{I}(\mathcal{L})f$ of a signal $f$ for any subset $I\subset I_{\mathcal{L}}$ which can be derived from the Chebyshev expansion of the indicator function $g=\mathrm{1}_{I}$. This observation actually appears in several recent works [19, 22]. More importantly for our study, this efficient estimation is part of the Hutchinson stochastic trace estimator technique [31], providing us with an effective method to estimate $n_{I}=\rm{Tr}(\mathcal{L}_{I})$. Finally, the present paper focuses on the computation of a sequence $g(\mathcal{L}_{I_{k}})_{1\leq k\leq K}$ (or its vector counterpart $g(\mathcal{L}_{I_{k}})f$) instead of a single $g(\mathcal{L})$ (resp. $g(\mathcal{L})f$). While a naive estimation would suggest that the computational complexity is then multiplied by a factor $K$ compared to the complexity of the computation of $g(\mathcal{L})$, we argue in the following that there is in fact no significant computational overhead.

4.1.2 Sharing Chebyshev polynomial matrices among filters

Let us assume that it is needed to compute the estimated values of $g_{k}(\mathcal{L})f$ for a given signal $f$ for several filters $g_{k}$, $k=0,\ldots,K$. Then the following two-step strategy can be adopted: (1) pre-compute Chebyshev expansions $\widetilde{g}_{k}(x)\approx\widetilde{g}_{k,N}(x)=\sum_{0\leq i\leq N}a_{i}(\widetilde{g_{k}})T_{i}(x)$ for all $k=0,\ldots,K$; independently, compute Chebyshev approximation vectors $T_{i}(\widetilde{\mathcal{L}})f$ for all $0\leq i\leq N$; (2) combine the previous results to compute the Chebyshev approximation $g_{k,N}(\mathcal{L})f$ of $g_{k}(\mathcal{L})f$:

[TABLE]

The complexity of the first step is dominated by the $N$ matrix-vector multiplications required to obtain $T_{i}(\widetilde{\mathcal{L}})f$. So the first step has complexity $O(|\mathcal{E}|N)$. The second step adds $N$ weighted matrices $a_{i}(\widetilde{g}_{k})T_{i}(\widetilde{\mathcal{L}})$ together, which is an operation of complexity $O(Nn^{2})$ at most. As an important matter of fact, the overall complexity for this procedure is bounded by $O(|\mathcal{E}|N+Nn^{2})$, which is independent of the number of filters $g_{k}$, and the same as for the computation of $g(\mathcal{L})$.

Sharing matrices among filters has several examples of applications in the current paper:

Computation of $g(\mathcal{L}_{I_{k}})f$ for all $1\leq k\leq K$: the equation $g(\mathcal{L}_{I_{k}})f=g(\mathrm{1}_{I_{k}}(\mathcal{L})\mathcal{L})f$ holds so that we can consider filters $g_{k}(x)=g(\mathrm{1}_{I_{k}}(x)x)$. 2. 2.

Computation of $g(s\mathcal{L})f$ for several scale values $s$: consider filters of the form $g_{s}(x)=g(sx)$. 3. 3.

Computation of $n_{I_{k}}$ for all $1\leq k\leq K$: Hutchinson’s stochastic estimation computes averages of $f_{i}^{T}P_{I_{k}}(\mathcal{L})f_{i}$ for some random vectors $f_{i}$ ($i\leq n_{H}$) whose computational complexity is dominated by the approximation of vectors $P_{I_{k}}f_{i}$. Considering filters $g_{k}(x)=\mathrm{1}_{I_{k}}(x)$, and sharing random vectors $(f_{i})_{i}$ among all approximations of $n_{k}$, we end up with a complexity of $O(n_{H}N|\mathcal{E}|)$, independent of value $K$.

In particular, Algorithm 1 has complexity $O(n_{H}N|\mathcal{E}|+Nn^{2})$. Indeed, its efficiency is calibrated on the computations of sequences $(\|\widetilde{f}_{k}\|_{2})_{1\leq k\leq K}$ and $(n_{k})_{1\leq k\leq K}$ whose computational analysis was discussed previously. It is worth observing that values $n_{k}$ do not depend on signal $f$ and should be estimated only once in the case where several signals $\widetilde{f}_{1},\widetilde{f}_{2},\ldots$ are to be denoised.

4.1.3 Optimizing storage of LocLets coefficients

The storage of wavelet coefficients $(\mathcal{W}^{I_{k}}f)_{1\leq k\leq K}$ requires a priori $K$ times the storage cost associated with the original transform $\mathcal{W}f$. When matrix reduction techniques are used to compute wavelets transform [54], one may reduce the storage consumption of the localized SGWT by suitably choosing the impulse functions $(\delta_{m})_{m}$. For instance, assume that for each subset $I_{k}$ a Lanczos basis $(v_{m}^{k})_{m}$ of the subspace spanned by $\{\chi_{\ell},\ell\in I_{k}\}$ is given. Then the size of sequences $(v_{m}^{k})_{m}$ and $(v_{m}^{k})_{m,k}$ are respectively of order $O(|I_{k}\cap{\rm sp}(\mathcal{L})|)=O(n_{k})$ and $O(n)$. Thus, with impulse functions $(v_{m}^{k})_{n}$ in place of $\delta_{m}$ for transform $\mathcal{W}^{I_{k}}$, the storage requirements of localized transform $(\mathcal{W}^{I_{k}}f)_{1\leq k\leq K}$ and the original one $\mathcal{W}f$ are of the same order $O(Jn)$.

4.2 Connections with Well-know Frames

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

[TABLE]

A frame is said to be tight if $A=B$. This section gives two examples of frames introduced in the literature which can be realized as a LocLets representation and thus benefit from the advantages given by localization in the spectrum.

4.2.1 Parseval frames

Parseval frames are powerful representations to design wavelets with nice reconstruction properties [34, 27]. In this section, we investigate the extent to which Parseval frames can be obtained from some LocLet representation. We show that for a particular choice of partition $I_{\mathcal{L}}=\sqcup_{k}I_{k}$, there exist frames which are Parseval frames and composed only of LocLets functions.

A finite collection $(\psi_{j})_{j=0,\ldots,J}$ is a finite partition of unity on the compact $[0,\lambda_{1}]$ if

[TABLE]

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]

Some constructions of partition of unity involve functions $(\psi_{j})_{j}$ that have almost pairwise disjoint supports i.e. ${\rm supp}(\psi_{j})\cap{\rm supp}(\psi_{j^{\prime}})=\emptyset$ as soon as $|j-j^{\prime}|>1$. For such partition of unity, set $I_{0}={\rm supp}(\psi_{0})$, $I_{J}=I_{0}={\rm supp}(\psi_{J})$ and $I_{j}={\rm supp}(\psi_{j})\cap{\rm supp}(\psi_{j+1})$ for all $1\leq j\leq J-1$. Then, the sequence $(I_{j})_{0\leq j\leq J}$ defines a finite partition of $[0,\lambda_{1}]$, $[0,\lambda_{1}]=\sqcup_{0\leq j\leq J}I_{j}$, such that:

[TABLE]

An alternative tight frame can be constructed using a LocLet representation as shown in the following proposition.

Proposition 12**.**

Assume Equations (11) and (12) hold and set, for all $1\leq k\leq J$, $\varphi_{n,k}=\sqrt{\psi_{0}}(\mathcal{L}_{I_{k}})\delta_{n}$, $\psi_{1,n,k}=\sqrt{\psi_{k}}(\mathcal{L}_{I_{k}})\delta_{n}$ and $\psi_{2,n,k}=\sqrt{\psi_{k}}(\mathcal{L}_{I_{k+1}})\delta_{n}$ for all $1\leq k\leq J$. Then $(\varphi_{n,k},\psi_{j,n,k})_{1\leq j\leq 2,~{}1\leq m\leq n,~{}1\leq k\leq J}$ is a tight frame.

The resulting tight frame of Proposition 12 is actually frame of LocLets if additionally the functions $\psi_{j}$ is of the form $\psi_{j}=\psi_{1}(s_{j}.)$ for some scale parameter $s_{j}$, $1\leq j\leq J$. This is typically the case for the frames introduced in [34, 27]. In these papers, the partition of unity is 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}]$ and set $\psi_{0}(\cdot)=\omega(\cdot)$ and for $j=1,\ldots,J$

[TABLE]

In particular, the functions $\psi_{k}$ have supports in intervals $J_{k}=[b^{k-2},b^{k}]$. Thus, one may define disjoint intervals $(I_{k})_{k}$ as follows: $I_{k}=[b^{k-1},b^{k}]$. We have $J_{k}=I_{k}\cup I_{k+1}$, so that Equations (12) hold whereas the scaling property $\psi_{j}=\psi_{1}(b^{-1}.)$ is straightforward. By Proposition 12, the set of vectors

[TABLE]

is a tight frame of LocLets. Observe that the transform $(\mathcal{W}^{I_{k}})_{I_{k}}$each component $\mathcal{W}^{I_{k}}$ of the LocLet transform $(\mathcal{W}^{I_{k}})_{I_{k}}$ only admit two scale parameters $s_{k},s_{k-1}$.

4.2.2 Spectrum-adapted tight frames

Let us consider another family of tight frames tailored to the distribution of the Laplacian $\mathcal{L}$ eigenvalues proposed in [51]. As shown below, these frames can be written in terms of a warped version of LocLets, and up to some approximation, in terms of (non-warped) LocLets. First, let us briefly recall the construction from [51].

The notion of warped SGWT is introduced in [51] to adapt the kernel to the spectral distribution. Given a warping function $\omega:I_{\mathcal{L}}\rightarrow\mathbb{R}$, the warped SGWT is defined as:

[TABLE]

As for our spectral localization, the objective of warping is to take benefits from the distribution of ${\rm sp}(\mathcal{L})$ along interval $I_{\mathcal{L}}$. While the two techniques show similarities (e.g. estimation of ${\rm sp}(\mathcal{L})$ distribution), they are meant to answer different problems: warped SGWT is a technique to adapt the whole spectrum to some task (e.g. producing a tight frame), whereas localized SGWT is designed to answer problems related to localized subsets in the spectrum (e.g. denoising a frequency sparse signal). Here we show that the advantages of both LocLets and warped SGWT are obtained when the two methods are combined in a warped LocLet representation.

Let $\omega$ be some warping function on $I_{\mathcal{L}}$ chosen in the form $\omega(\cdot)=\log(C\omega_{0}(\cdot))$ where $\omega_{0}$ stands for the cumulative spectral distribution of $\mathcal{L}$ and $C$ is some normalization constant as shown in [51]. Then, let $\gamma>0$ be an upper bound on $\rm{sp}(\mathcal{L})$ and let $R,J$ be two integers such that $2\leq R\leq J$. Setting $\omega_{\gamma,J,R}=\frac{\gamma}{J+1+R}$, Corollary 2 in [51] asserts that the family $(g_{m,j})_{m,j}$ of functions defined below is a tight frame

[TABLE]

where functions $\widehat{g}_{j}$ arise from some kernel $\widehat{g}$ as

[TABLE]

Typically in [51], the kernel $\widehat{g}$ takes the form

[TABLE]

for some sequence $(a_{j})_{j}$ satisfying $\sum_{j}(-1)^{j}a_{j}=0$.

The following proposition states that Equation (13) admits an alternative form involving only (warped) LocLets functions.

Proposition 13**.**

Setting $\psi(\lambda)=\widehat{g}(\log(C\lambda))$ for $\lambda>0$, consider the family of warped LocLets defined for all $0\leq k\leq R-1$, $1\leq m\leq n$ and $1\leq j\leq J$ by

[TABLE]

Then, the following identity holds for all $j=1,\ldots,J$ and all $m=1,\ldots,n$

[TABLE]

5 Experiments on suites of large matrices

This section details experiments made on large graphs to validate the Fourier localization techniques introduced in that paper. After describing the experimental settings, we describe the outcomes of several experiments showing strong advantages in the use of Fourier localization in practice.

5.1 Choice of spectral partition $I_{\mathcal{L}}=\sqcup_{k}I_{k}$

In order to keep the problem combinatorially tractable, it is necessary to reduce the choice of possible partitions of $I_{\mathcal{L}}$ into subintervals $I_{k}$. That is why, the partitions considered in the sequel are regular in the sense that all intervals have the same length $\lambda_{1}/K$ for some integer $K\geq 1$. Thereafter, the parameter $K$ is chosen so that the eigenvalues are distributed as evenly as possible in each interval $I_{k}$. Without prior information, it is indeed natural not to favor one part of the spectrum over another. Most importantly, in the view of the concentration property of the median around the noise level $\sigma^{2}$ of Proposition 5, it is essential to keep the parameter $\beta$ as close to one as possible.

In order to implement the ideas above, it is necessary to estimate the spectral measure of $\mathcal{L}$ which can be described by the so-called spectral density function:

[TABLE]

There are several techniques for such an approximation among which the Kernel Polynomial Method (see, e.g. [52, 61]). The latter approximates the spectral density $\varphi_{\mathcal{L}}$ with the help of matrix Chebyshev expansion $(\varphi_{\mathcal{L}}^{N})_{N}$ (see [35] for a detailed presentation).

Now, let $(I_{k})_{1\leq k\leq K}$ be some regular partition of $I_{\mathcal{L}}$ and $(n_{k})_{1\leq k\leq K}$ be the corresponding numbers of eigenvalues in each $I_{k}$. Choosing the parameter $K\geq 1$ so that the entropy defined by

[TABLE]

is maximal ensures that the eigenvalues are as equally distributed in each interval as possible. In application, the Kernel Polynomial Method provides an approximation $n_{k}^{N}$ of $n_{k}$ and the corresponding empirical entropy $E_{N}(K)$ is used as a proxy for the theoretical one.

Empirically, the entropy increases logarithmically and then stabilizes from a certain elbow value $K_{\rm elbow}$ as illustrated in Figure 1. This elbow value is displayed in dashed lines in Figure 1. In the experiments, we choose this value $K_{\rm elbow}$ motivated by two reasons. First, as the intervals become shorter it is more difficult to obtain a uniform distribution of the eigenvalues into those intervals. The second reason is related to the quality of the estimate $n_{k}^{N}$ of $n_{k}$ as the sample size decreases. To illustrate this fact, we consider the Mean Relative Error (MRE) defined by

[TABLE]

As highlighted by Figure 1, the empirical entropy actually stabilizes when the Chebyshev approximation, in terms of MRE, is no longer sharp enough.

5.2 The experimental settings

Following [22], we propose to validate our techniques on an extended suite of large matrices extracted from the Sparse Matrix Collection in [16]. Most of these matrices have an interpretation as the Laplacian matrix of a large graph. We define matrix $\mathcal{L}$ from the following matrices of the suite:si2 ($n=769$), minnesota ($n=2642$), cage9 ($n=3534$), saylr4 ($n=3564$) and net25 ($n=9520$). We extend this graph collection with the well-studied swissroll graph Laplacian matrix ($n=1000$).

We sample randomly signals whose supports are sparse in the Fourier domain. We will use the notation $f_{i-j}$ for normalized signals supported on a sub-interval of $I_{\mathcal{L}}$ containing exactly the eigenvalues $\lambda_{i},\lambda_{i+1},\ldots,\lambda_{j}$. As an example, $f_{n-n}$ is a constant signal while $f_{1-2}$ is a highly non-smooth signal supported on the eigenspaces of large eigenvalues $\lambda_{1},\lambda_{2}$. For experiments, the signals were calculated from the knowledge of $\mathrm{sp}(\mathcal{L})$, and relevant projections of random functions on the graph.

We have compared the performances of Algorithms 2 and 3 against the thresholding procedure described in [27]. As the denoising method in [27] requires the computation of the whole spectral decomposition of the Laplacian, it does not scale to large graphs. We stress here that we provide a fair comparison with [27], only in terms of denoising performance, and with no computational considerations. Moreover, we choose for LocLets to use the most naive thresholding procedure by considering a global and scale independent threshold level.

For all the experiments below, the SGWT and LocLets are built upon the scale and kernel functions giving rise to the Parseval frame of  [27], whose construction is recalled in Section 4.2.1. More precisely, set respectively $\varphi=\sqrt{\zeta_{0}}$ and $\psi=\sqrt{\zeta_{1}}$ for the scale and kernel functions with $\zeta_{0}(x)=\omega(x)$, and $\zeta_{1}(x)=\omega(b^{-1}x)-\omega(x)$, where we choose $b=2$ and $\omega$ is piecewise linear, vanishes on $[1,\infty)$ and is constant equal to one on $(-\infty,b]$. The scales are of the form $s_{j}=b^{-j+1}$ for $j=1,\ldots,J$ where $J$ is chosen similarly to [27].

In what follows, ‘PF’ stands for Parseval Frame and refers to the estimator of [27]; the estimators implemented by Algorithm 2 and Algorithm 3 are referred to as ‘LLet’ and ‘LLet+PF’ respectively. The notation ‘$\mathrm{SNR}_{\mathrm{in}}$’ refers to the trivial model releasing the noisy signal $\widetilde{f}$, corresponding to the classical input noise level measurement, and serves as a worst-case baseline for other models. Below, the latter methodology is shown to outperform all the others for very frequency-sparse signals. It is also worth recalling that ‘LLet+PF’ benefits from the dimension reduction property of LocLets. More precisely, whereas the whole eigendecomposition of $\mathcal{L}$ is required to apply ‘PF’, for Parseval frame denoising in the context of ‘LLet+PF’, only a low-rank spectral decomposition is needed, namely the decomposition of $\mathcal{L}_{I}$ for $I$ the estimate of ${\rm supp}\widehat{f}$.

For all our experiments, we set $\alpha=0.001$ for Algorithm 1. For the denoising experiment, we compute the best SNR result $r_{D}$ over a large grid of values $(t_{1},t_{2})$, and for each denoising method $D$ with $D\in\{`\mathrm{SNR}_{\mathrm{in}}\textrm{',`PF',`LLet',`LLet+PF'}\}$. Then, we calculate two metrics: the maximum $\mathrm{M}_{D}$ and average value $\mu_{D}$ of the values $r_{D}$ over $10$ random perturbations of the signal $f$. We recall that a good quality in denoising is reflected by a large value of the SNR metric.

5.3 Analysis of our experiments

5.3.1 Noise level estimation

We have evaluated the performances of estimators $\widehat{\sigma}_{\rm mean}^{r}$ and $\widehat{\sigma}_{\rm med}$ in the estimation of the unknown noise level $\sigma$ from $10$ realizations of the noisy signal $\widetilde{f}=f+\xi$ for a given noise level $\sigma$. Figure 2 (resp. Figure 3 ) shows the best performances of each estimator on the minnesota (resp. net25) graph for the non-regular but frequency sparse signal $f_{1392-1343}$ (resp. $f=f_{4971-5020}$), when parameter $K$ ranges in $\{5,10,20,30,40,50\}$ and for level of noise $\sigma=0.01$ (resp. $\sigma=0.001$).

Figure 2 illustrates that both estimators $\widehat{\sigma}_{\rm mean}^{r}$ and $\widehat{\sigma}_{\rm med}$ can provide good estimates of $\sigma$. Best performances are obtained for values of parameter $K$ below the elbow value $K_{elbow}(minnesota)=22$ introduced in Section 5.1. We observe that performances drop considerably if almost no localization is used (for instance, for parameter values $K=1$ or $K=2$, $\widehat{\sigma}\sim 0.021$ in the experiment of Figure 2, far from the performances for $K\geq 5$ for estimating $\sigma=0.01$).

Figure 3 shows that localization is necessary, namely $K\geq 10$ or even $K\geq 20$, in order to reach the best performances for the large net25 graph. Contrary to experiments for the minnesota graph, estimators $\widehat{\sigma}_{\rm mean}^{r}$ and $\widehat{\sigma}_{\rm med}$ underestimate the value of $\sigma$. Also, best values of $K$ range between $10$ and $30$ for net25 graph, compared to best values $K=5$ and $K=10$ for minnesota graph (see Figure 2). This illustrates the idea that noise level estimation strongly depends on the underlying graph structure. As a consequence, a parameter $K$ selection has to take graph and signal information into account to be relevant. Interestingly, the elbow values $K_{elbow}(minnesota)=22$ and $K_{elbow}(net25)=22$ provide performances which are not optimal, but close to the best possible ones.

In Figure 4, performances for various values of parameter $r$ are displayed for a fixed parameter $K=K_{elbow}(minnesota)$. While it is true that $\widehat{\sigma}_{\rm mean}^{r}$ can perform better than $\widehat{\sigma}_{\rm med}^{r}$, it happens only for very specific values of $r$, which a priori depend on the signal regularity. Without any further parameter selection, these observations suggest using the most robust estimator $\widehat{\sigma}_{\rm med}$ in practice.

5.3.2 Sparse signal denoising

As a first denoising experiment, we have compared the performances of ‘LLet’, ‘PF’ and ‘LLet+PF’ for a fixed value $K=K_{elbow}$ given by the rule of thumb described in Section 5.1. For each matrix in the Extended Matrices Suite, we have experimented the denoising task on two frequency-sparse signals, one regular and the other non-regular. Several values of noise level $\sigma$ were used, corresponding to values of $\mathrm{SNR}_{\mathrm{in}}$ ranging in $[4,18]$. Results from experiments are displayed in Tables 1 and 2. The first obvious observation is that ‘LLet+PF’ performs better than its competitors in almost all situations. The gain is sometimes considerable since we observed a gap of $5$dB in $\mu_{D}$-metric between ‘LLet+PF’ and its closest concurrent ‘PF’ in some cases, and up to $7$dB in $\mathrm{M}_{D}$-metric. These experiments confirm the theoretical guarantees obtained in Theorem 11. The benefits of localization are reduced for graph net25: Table 2 shows that the more conservative choice $K=5$ is better than $K=25$. It appears that for net25, the spectrum $\mathrm{sp}(\mathcal{L})$ is localized at a small number of distinct eigenvalues, hence diminishing the advantages of localizing with our methods.

Another interesting observation is that ‘LLet’ may outperform ‘PF’ in some specific signal and noise level configurations, as shown in Table 1. This is a very favorable result for localized Fourier analysis, since ‘LLet’ appears to be a technique which is more accurate and more efficient as well compared to ‘PF’ in some situations. However in many cases, ‘LLet’ performances drop down compared to the more stable thresholding techniques ‘PF’ and ‘LLet+PF’, which use thresholds adapted to the wavelet basis.

We also provide experimental results to understand the extent to which our results depend on the partition size parameter $K$. A few remarks are suggested by Figure 5:

• •

The best performances are not obtained for the elbow value $K_{elbow}$, suggesting searching for a more task-adapted size of partition $K$.

• •

Good performances persist for values of $K$ much larger than $K_{elbow}$, and in particular for regular signals.

• •

For large values of $K$, there is a severe drop in performances. As explained before, the error generated by Chebyshev’s approximation grows with the number of intervals in the partition, which makes the approximation of the support more difficult.

6 Conclusion and future works

We have introduced a novel technique to efficiently perform graph Fourier analysis. This technique uses functional calculus to perform Fourier analysis on different subsets of the graph Laplacian spectrum. In this paper, we have demonstrated that localization in the spectrum provides interesting improvements in theoretical results for some graph signal analysis tasks. New estimators of the noise level were introduced, taking advantage of the convenient modelling of the denoising problem given by localization, and for which concentration results were proved. Localization allows also to study theoretically the denoising procedure with wavelets, and fits with the design of many well-known techniques (e.g. tight frames for graph analysis). Through many experiments, we have validated that localization techniques introduced in this paper improve on state-of-the-art methods for several standard tasks.

Although we provide a rule of thumb to choose a partition $I_{\mathcal{L}}=\sqcup_{1\leq k\leq K}I_{k}$ for which denoising results show good performances, experiments suggest that our elbow rule is not optimal in most cases. There is certainly an interesting topic in searching for a suitable partition $I_{\mathcal{L}}=\sqcup_{1\leq k\leq K}I_{k}$ that would be more adapted to a specific task (e.g. denoising). To extend the current work, it would also be interesting to consider other common tasks in GSP, such as de-convolution or in-painting.

7 Proofs

Proof of Proposition 4.

We have $\tilde{f}_{k}=f_{k}+\xi_{k}$, where $\xi_{k}=\sum_{\ell:\lambda_{\ell}\in I_{k}}\widehat{\xi}(\ell)\chi_{\ell}$. Random variables $(\widehat{\xi}(\ell))_{\ell}$ are all distributed as $\mathcal{N}(0,\sigma^{2})$ and independent by orthogonality of the eigenbasis $(\chi_{\ell})_{\ell}$. In particular for $k\neq k^{\prime}$, vectors $\xi_{k}$ and $\xi_{k^{\prime}}$ are independent as expressions involving variables $\widehat{\xi}(\ell)$ over disjoint subsets $I_{k}$ and $I_{k^{\prime}}$. Thus the random variables $(c_{k})_{1\leq k\leq K}$ are also independent. 2. 2.

When $n_{k}=n_{k^{\prime}}$, $\xi_{k}$ and $\xi_{k^{\prime}}$ are identically distributed and the result follows from Equality (6).

When $n_{k}\neq n_{k^{\prime}}$, we have $\mathbb{E}(c_{k})\neq\mathbb{E}(c_{k^{\prime}})$ as the following equality holds for all $1\leq k\leq K$:

[TABLE] 3. 3.

Since $(\widehat{\xi}(\ell))_{\ell}$ are independent normal variables $\mathcal{N}(0,\sigma^{2})$, the statement is clear from the expression $c_{k}=\frac{1}{n_{k}}\sum_{\ell:\lambda_{\ell}\in I_{k}}|\widehat{\xi}(\ell)|^{2}$.

∎

The following lemma is useful for the proof of Proposition 5.

Lemma 14**.**

Let $Z\sim\mathcal{B}(n,p)$ for some parameters $n\geq 1$ and $p\leq 1/2$. Then

[TABLE]

Proof.

A simple consequence of [10, Theorem 1] implies that for all $n\geq 1$ and all $a\geq p$

[TABLE]

Now the result follows since $na=\lceil n/2\rceil$ implies $\frac{1}{2}\leq a\leq\frac{1}{2}+\frac{1}{n}$ so that

[TABLE]

∎

Proof of Proposition 5.

For $k$ such that $I_{k}\subset\overline{I_{f}}$, we have $c_{k}=\frac{\sigma^{2}}{n_{k}}\Gamma_{n_{k}}$, which follows the $\Gamma(\frac{n_{k}}{2},\frac{2\sigma^{2}}{n_{k}})$ distribution. Then concentration inequalities for $\widehat{\sigma}_{\rm mean}(c)^{2}$ are a direct consequence of Theorem 2.57 in [4], applied with $a_{k}=\frac{n_{k}}{2}$ and $b_{k}=\frac{2\sigma^{2}}{n_{k}}$. 2. 2.

For all $k=1,\ldots,K_{f}$, we define

[TABLE]

As a matter of fact, $(\gamma^{-})_{k=1,\ldots,K_{f}}$ and $(\gamma_{k}^{+})_{k=1,\ldots,K_{f}}$ are two sequences of i.i.d. random variables with $\gamma_{1}^{-}\sim\chi^{2}(n_{0})$ and $\gamma^{+}_{1}\sim\chi^{2}(n_{\infty})$ such that

[TABLE]

Then, for all $t>0$,

[TABLE]

Similarly, for all $t\in(0,1)$,

[TABLE]

To conclude, apply Lemma 14 to Inequalities (14) and (15) to obtain our result.

∎

Proof of Proposition 6.

The concentration bound of Equation (8) implies that

[TABLE]

Since subsets $I_{k}$ and $I_{\ell}$ are disjoint, random variables $\|\xi_{k}\|_{2}^{2}$ and $\|\xi_{\ell}\|_{2}^{2}$ are independent. Thus, $\frac{n_{k}}{n_{\ell}}\|\xi_{\ell}\|_{2}^{2}-\|\xi_{k}\|_{2}^{2}$ is distributed as $\frac{n_{k}}{n_{\ell}}\Gamma_{n_{\ell}}-\Gamma_{n_{k}}$ where $\Gamma_{n_{k}}$ and $\Gamma_{n_{\ell}}$ are independent random variables with $\Gamma_{n_{k}}\sim\chi^{2}(n_{k})$ and $\Gamma_{n_{\ell}}\sim\chi^{2}(n_{\ell})$. Therefore, the statement of Proposition 6 follows. ∎

Proof of Proposition 7.

First, the following equalities hold:

[TABLE]

As $(f-\widetilde{f})_{I}$ and $\widetilde{f}_{\overline{I}}$ are orthogonal vectors, it follows that

[TABLE]

It remains to notice that $\mathbb{E}(||\widetilde{f}_{\overline{I}}||^{2})=\sigma^{2}|\overline{I}\cap\mathrm{sp}(\mathcal{L})|$. ∎

Proof of Lemma 8.

By Equation (6) and the concentration bound of Equation (8), it follows that

[TABLE]

where $\theta(f_{k},\alpha,\sigma)=\sigma^{-2}(\|f_{k}\|_{2}-2t_{\alpha/2,\sigma})\|f_{k}\|_{2}$. Consequently, $1-\Phi_{\Gamma_{n_{k}}-\Gamma_{n_{k}}^{\prime}}(\theta(f_{k},\alpha,\sigma))\leq\alpha/2$ and $p_{k}\leq\alpha$. ∎

Proof of Lemma 9.

Using an estimate on the $\chi^{2}(n_{k})$ tail distribution and independence of $\Gamma_{n_{k}}$ and $\Gamma_{n_{k}}^{\prime}=\sigma^{-2}\|\xi_{k}\|_{2}^{2}$, it follows

[TABLE]

∎

Proof of Proposition 10.

To prove Inequality (9), first observe that $f=\sum_{k:I_{k}\subset I_{f}}f_{k}$ so that

[TABLE]

The summands which are not present in both terms are exactly those satisfying either $I_{k}\subset I_{f}$ and $p_{k}>\alpha$ or $I_{k}\cap I_{f}=\emptyset$ and $p_{k}\leq\alpha$. Noting that $f_{k}=0$ when $I_{k}\cap I_{f}=\emptyset$, it comes

[TABLE]

Applying Lemma 8 for all indices $1\leq k\leq K$ satisfying $p_{k}>\alpha$, one deduce

[TABLE]

from which, since $n_{k}=n_{1}$ for all $k$, Inequality (9) follows. 2. 2.

Since $\sigma^{-2}||\xi_{k}||_{2}^{2}$ is distributed as a $\chi^{2}(n_{1})$ random variable, the second Inequality (10) follows

[TABLE]

∎

Proof of Theorem 11.

Since threshold value is $t_{2}=0$ on $\overline{I}$, $\widehat{f}=\widehat{f}_{I}$. Then, clearly $f_{\overline{I}}(f_{I}-\widehat{f}_{I})=0$ almost surely so that

[TABLE]

Applying Theorem 3 from [27] to $\mathbb{E}\|f_{I}-\widehat{f_{I}}\|_{2}^{2}$ yields our statement. ∎

Proof of Proposition 12.

Recalling that, for any function $g$ defined on $\mathrm{sp}(\mathcal{L})$ and any subset $I\subset I_{\mathcal{L}}$,

[TABLE]

it follows by Equations (11) and (12).

[TABLE]

∎

Proof of Proposition 13.

Remarking that $\widehat{g}_{j}(\lambda)=\psi(s_{j}\omega_{0}(\lambda))$ with $s_{j}=e^{-j\omega_{\gamma,J,R}}$, Equation (13) implies that

[TABLE]

Setting $J_{j}=[C^{-1}e^{(j-R)\omega_{\gamma,J,R}},C^{-1}e^{j\omega_{\gamma,J,R}}]$ and recalling that ${\rm supp}(\widehat{g})=[-R\omega_{\gamma,J,R},0]$, it follows that $\lambda\in\mathrm{supp}(\widehat{g}_{j})$ if and only if $\omega_{0}(\lambda)\in J_{j}$ if and only if $s_{j}\omega_{0}(\lambda)\in J_{0}$. Moreover, $J_{0}=\sqcup_{k}I_{k}$ and $J_{j}=s_{j}^{-1}J_{0}$ yield $J_{j}=\sqcup_{k}s_{j}^{-1}I_{k}$ with $s_{j}^{-1}I_{k}\cap s^{-1}_{j^{\prime}}I_{k^{\prime}}$ excepted when $j=j^{\prime}$ ad $k=k^{\prime}$. Consequently, Equation (16) can be reformulated as

[TABLE]

∎