Irregularity-Aware Graph Fourier Transforms Benjamin Girault, Antonio Ortega Fellow, IEEE and Shrikanth S. Narayanan Fellow, IEEE Abstract-In this paper, we present a novel generalization of the graph Fourier transform (GFT). Our approach is based on separately considering the definitions of signal energy and signal variation, leading to several possible orthonormal GFTs. Our approach includes traditional definitions of the GFT as special cases, while also leading to new GFT designs that are better at taking into account the irregular nature of the graph. As an illustration, in the context of sensor networks we use the Voronoi cell area of vertices in our GFT definition, showing that it leads to a more sensible definition of graph signal energy even when sampling is highly irregular.

Index Terms-Graph signal processing, graph fourier transform.

I. INTRODUCTION

With the ever growing deluge of data, graph signal processing has been proposed for numerous applications in recent years, thanks to its ability to study signals lying on irregular discrete structures. Examples include weather stations [START_REF] Girault | Stationary Graph Signals using an Isometric Graph Translation[END_REF], taxi pickups [START_REF] Chen | Monitoring Manhattan's Traffic at 5 Intersections?[END_REF] or bicycle rides [START_REF] Hamon | Networks as signals, with an application to a bike sharing system[END_REF], [START_REF] Hamon | Analysis of temporal networks using signal processing methods : Application to the bike-sharing system in Lyon, Theses[END_REF], people in a social network [START_REF] Girault | Graphe de contacts et ondelettes : étude d'une diffusion bactérienne[END_REF], or motion capture data [START_REF] Kao | Graphbased approach for motion capture data representation and analysis[END_REF]. While classical signal processing is typically built on top of regular sampling in time (1D Euclidean space), or space (2D Euclidean space), graphs can be applied in irregular domains, as well as for irregular sampling of Euclidean spaces [START_REF] Belkin | Towards a theoretical foundation for Laplacian-based manifold methods[END_REF]. Thus, graph signal processing can be used to process datasets while taking into consideration irregular relationships between the data points.

Successful use of graph signal processing methods for a given application requires identifying: i) the right graph structure, and ii) the right frequency representation of graph signals.

The choice of graph structure has been studied in recent work on graph learning from data [START_REF] Samuel I Daitch | Fitting a graph to vector data[END_REF]- [START_REF] Mei | Signal processing on graphs: Causal modeling of unstructured data[END_REF]. In this paper we focus on the second question, namely, given a graph structure, how to extend the classical Fourier transform definition, which relies on a regular structure, to a graph Fourier transform, which is linked to a discrete irregular structure.

State of the art methods derive the definition of the graph Fourier transform (GFT) from algebraic representations of the graph such as the adjacency matrix A, whose entries are the weights of the edges connecting vertices. If i and j are two vertices connected through an edge weighted by

The authors are with the Signal and Image Processing Institude, University of Southern California, Los Angeles, CA 90089, USA (firstname.lastname@usc.edu).

Work partially funded by NSF under grants CCF-1410009, CCF-1527874, CCF-1029373.

The research is based upon work supported by the Office of the Director of National Intelligence (ODNI), Intelligence Advanced Research Projects Activity (IARPA), via IARPA Contract No 2017-17042800005. The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of the ODNI, IARPA, or the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for Governmental purposes notwithstanding any copyright annotation thereon. w(ij), then A ij = w(ij). In the context of sensor networks, edges are often defined by selecting the K-nearest neighbors to each vertex, with weights given by a Gaussian kernel of the Euclidean distance. This is a setting shown to have good properties in the context of manifold sampling [START_REF] Hein | Graph Laplacians and their Convergence on Random Neighborhood Graphs[END_REF] when the number of samples is large. The adjacency matrix is then used to compute the Laplacian matrix L = D -A with D the diagonal degree matrix verifying D ii = d i = j w(ij). Conventional methods in graph signal processing use the eigenvectors of either A [START_REF] Sandryhaila | Discrete Signal Processing on Graphs[END_REF] or L [START_REF] Shuman | The Emerging Field of Signal Processing on Graphs: Extending High-Dimensional Data Analysis to Networks and Other Irregular Domains[END_REF] as graph Fourier modes.

One motivation of this paper comes from physical sensing applications, where data collection is performed by recording on a discrete set of points, such as weather readings in a weather station network. In these applications the distribution of the sensors is irregular, and so it is not easy to map these measurements back onto a grid in order to use the classical signal processing toolbox. Importantly, in many of these applications, the specific arrangement of sensors is unrelated to the quantity measured. As an obvious example, the exact location and distribution of weather stations does not have any influence over the weather patterns in the region. Thus, it would be desirable to develop graph signal processing representations that i) have a meaningful interpretation in the context of the sensed data, and ii) are not sensitive to changes in graph structure, and in particular its degree of regularity. 1In order to motivate our problem more precisely on this sensor network example, consider a popular definition of GFT, based on the eigenvalues of L, with variation for a graph signal x defined as:

x H Lx = 1 2 ij w(ij) |x i -x j | 2 .
A signal x with high variation will have very different values (i.e., large |x ix j | 2 ) at nodes connected by edges with large weights w(ij). This GFT definition is such that variation for constant signals is zero, since L1 = 0. This is a nice property, as it matches definitions of frequency in continuous domain, i.e., a constant signal corresponds to minimal variation and thus frequency zero. Note that this property is valid independently of the number and position of sensors in the environment, thus achieving our goal of limiting the impact of graph choice on signal analysis.

In contrast, consider the frequency representation associated to an impulse localized to one vertex in the graph. The impulse signal δ i has variation δ H i Lδ i = d i , where d i is the degree of vertex i. Impulses do not have equal variation, so that a highly localized phenomenon in continuous space (or a sensor anomaly affecting just one sensor) would produce different spectral signatures depending on the degree of the vertex where the measurement is localized. As a consequence, if the set of sensors changes, e.g., because only a subset of sensors are active at any given time, the same impulse at the same vertex may have a significantly different spectral representation.

As an alternative, choosing the symmetric normalized Laplacian L = D -1/2 LD -1/2 to define the GFT would lead to the opposite result: all impulses would have the same variation, but a constant signal would no longer correspond to the lowest frequency of the GFT. Note the importance of graph regularity in this trade-off: the more irregular the degree distribution, the more we deviate from desirable behavior for the impulses (L) or for the lowest frequency (L).

As further motivation, consider the definition of signal energy, which is typically its 2 -norm in the literature:

E x = x 2
2 [14]- [START_REF] Girault | Translation on Graphs: An Isometric Shift Operator[END_REF]. Assume that there is an area within the region being sensed where the energy of the continuous signal is higher. For this given continuous space signal, energy estimates through graph signal energy will depend significantly on the position of the sensors (e.g., these estimates will be higher if there are more sensors where the continuous signal energy is higher). Thus, estimated energy will depend on the choice of sensor locations, with more significant differences the more irregular the sensor distribution is.

In this paper, we propose a novel approach to address the challenges associated with irregular graph structures by replacing the dot product and the 2 -norm, with a different inner product and its induced norm. Note that here irregularity refers to how irregular the graph structure is, with vertices having local graph structures that can vary quite significantly. In particular, a graph can be regular in the sense of all vertices having equal degree, yet vertices not being homogeneous w.r.t. to some other respect. Different choices of inner product can be made for different applications. For any of these choices we show that we can compute a set of graph signals of increasing variation, and orthonormal with respect to the chosen inner product. These graph signals form the basis vectors for novel irregularity-aware graph Fourier transforms that are both theoretically and computationally efficient. This framework applies not only to our motivating example of a sensor network, but also to a wide range of applications where vertex irregularity needs to be taken into consideration.

In the rest of this paper, we first introduce our main contribution of an irregularity-aware graph Fourier transform using a given graph signal inner product and graph signal energy (Sec. II). We then explore the definition of graph filters in the context of this novel graph Fourier transform and show that they share many similarities with the classical definition of graph filters, but with more richness allowed by the choice of inner product (Sec. III). We then discuss specific choices of inner product, including some that correspond to GFTs known in the literature, as well novel some novel ones (Sec. IV). Finally, we present two applications of these irregularity-aware graph Fourier transforms: vertex clustering and analysis of sensor network data (Sec. V).

II. GRAPH SIGNAL ENERGY: LEAVING THE DOT PRODUCT

One of the cornerstones of classical signal processing is the Fourier transform, and one of its essential properties is orthogonality. Indeed, this leads to the Generalized Parseval Theorem: The inner product in the time x, y and spectral x, y domains of two signals are equal. In this section, we propose a generalization to the definition of the GFT. Our key observation is to note that the choice of an inner product is a parameter in the definition of the GFT, and the usual dot product is not the only choice available for a sound definition of a GFT verifying Parseval's Theorem.

A. Graph Signal Processing Notations

Let G = (V, E, w) be a graph with V = {1, . . . , N } its vertex set, 2 E ⊆ V × V its edge set, and w : E → R + the weight function of the edges (as studied in [START_REF] Girault | Signal Processing on Graphs -Contributions to an Emerging Field[END_REF], these weights should describe the similarity between vertices). For simplicity, we will denote the edge from i to j as ij. A graph is said undirected if for any edge ij, ji is also an edge, and their weights are equal: w(ij) = w(ji). A self-loop is an edge connecting a vertex to itself.

A graph is algebraically represented by its adjacency matrix 3 A verifying A ij = 0 if ij is not an edge, and A ij = w(ij) otherwise. 4 We use also a normalization of the adjacency matrix by its eigenvalue of largest magnitude µ max such that A (norm) = 1 µmax A [START_REF] Sandryhaila | Discrete Signal Processing on Graphs: Frequency Analysis[END_REF]. We denote d i = j A ij the degree of vertex i, and D the diagonal degree matrix having d as diagonal. When the graph is undirected, the Laplacian matrix is the matrix L = D -A. Two normalizations of this matrix are frequently used in the literature: the normalized Laplacian matrix L = D -1/2 LD -1/2 = I -D -1/2 AD -1/2 and the random walk Laplacian

L RW = D -1 L = I -D -1 A.
In this graph signal processing literature, the goal is to study signals lying on the vertices of a graph. More precisely, a graph signal is defined as a function of the vertices x : V → R or C. The vertices being indexed by {1, . . . , N }, we represent this graph signal using the column vector

x = [x 1 , . . . , x N ]
T . Following the literature, the GFT is defined as the mapping from a signal x to its spectral content x by orthogonal projection onto a set of graph Fourier modes (which are themselves graph signals) {u 0 , . . . , u N-1 } and such that

x i = l x(l) [u l ] i = [U x] i , where U = [u 0 • • • u N-1 ]. We denote F -1 = U the inverse GFT matrix and F = U -1 the GFT matrix.
Furthermore, we assume that there is an operator ∆ : C V → R + quantifying how much variation a signal shows on the Here GQV is the Graph Quadratic Variation (see Sec. IV-A3).

Name (M) ∆(x) [15] Comb. Lapl. (L) 1 2 ij w(ij)|x i -x j | 2 [15] Norm. Lapl. (L) 1 2 ij w(ij) x i √ d i - x j √ d j 2 [18] GQV ((I -A (norm) ) 2 ) i |x i -[A (norm) x] i | 2
graph. ∆ typically depends on the graph chosen. Several examples are given on Tabs. I and II. The graph variation ∆ is related to a frequency as variation increases when frequency increases. The graph frequency of the graph Fourier mode u l is then defined as λ l = ∆(u l ). 5In summary, the two ingredients we need to perform graph signal processing are i) the graph Fourier matrix F defining the projection on a basis of elementary signals (i.e., the graph Fourier modes), and ii) the variation operator ∆ defining the graph frequency of those elementary signals. Our goal in this paper is to propose a new method to choose the graph Fourier matrix F, one driven by the choice of graph signal inner product. For clarity and conciseness, we focus on undirected graphs, but our definitions extend naturally to directed graphs. An in-depth study of directed graphs will be the subject of a future communication.

B. Norms and Inner Products

In the literature on graph signal processing, it is desirable for the GFT to orthogonal, i.e., the graph Fourier transform is an orthogonal projection on an orthonormal basis of graph Fourier modes. Up to now, orthogonality has been defined using the dot product: x, y = y H x, with . H the transpose conjugate operator. We propose here to relax this condition and explore the benefits of choosing an alternative inner product on graph signals to define orthogonality.

First, note that ., . Q is an inner product on graph signals, if and only if there exists a Hermitian positive definite (HPD) matrix Q such that:

x, y Q = y H Qx, for any graph signals x and y. We refer to this inner product as the Q-inner product. With this notation, the standard dot product is the I-inner product. Moreover, the Q-inner product induces the Q-norm:

x Q = x, x Q = x H Qx.
Therefore, an orthonormal set for the Q-inner product is a set {u l } l verifying:

u k H Qu l = 0 if k = l 1 otherwise, i.e., U H QU = I with U = [• • • u l • • • ].
Although any HPD matrix Q defines a proper inner product, we will be mostly focusing on diagonal matrices Q = diag(q 1 , . . . , q N ). In that case, the squared Q-norm of a graph

Name

∆(x) [18] GTV i x i -[A (norm) x] i [19] GDV ij w(ij)[x i -x j ] +
signal x, i.e., its energy, is a weighted sum of its squared components:

x 2 Q = i∈V q i |x i | 2 .
Such a norm is simple but will be shown to yield interesting results in our motivating example of sensor networks (see Sec. V-B). Essentially, if q i quantifies how important vertex q i is on the graph, the energy above can be used to account for irregularity of the graph structure by correctly balancing vertex importance. Examples of diagonal matrix Q include I and D, but more involved definitions can be chosen such as with the illustrating example of a ring in Sec. II-D or our novel approach based on Voronoi cells in Sec. IV-D. Moreover, endowing the space of graph signals with the Qinner product yields the Hilbert space H G (Q) = (C V , ., . Q ) of graph signals of G. This space is important for the next section where we generalize the graph Fourier transform.

Remark 1. To simplify the proofs, we observe that many results from matrix algebra that rely on diagonalization in an orthonormal basis w.r.t. the dot product can actually be extended to the Q-inner product using the following isometric operator on Hilbert spaces:

ϕ : H G (Q) -→ H G (I) x -→ Q 1/2 x .
Since Q is invertible, ϕ is invertible, and an orthonormal basis in one space is mapped to an orthonormal basis in the other space, i.e., {u l } l is an orthonormal basis of H G (Q) if and only if {ϕ(u l )} l is an orthonormal basis of H G (I).

For example, if h is a linear operator of H G (Q) such that h(x) = Hx, then h : y → ϕ(h(ϕ -1 (y))) = Q 1/2 HQ -1/2 y is an operator of H G (I). Also, the eigenvectors of H are orthonormal w.r.t. the Q-inner product if and only if the eigenvectors of Q 1/2 HQ -1/2 are orthonormal w.r.t. the dot product. This applies for example to L RW = D -1 L the random walk Laplacian. Since L = D 1/2 (D -1 L)D -1/2 has orthonormal eigenvectors w.r.t. the dot product, then the eigenvectors of L RW are orthonormal w.r.t. the D-inner product (see Sec. IV-B).

C. Contribution 1: Generalized Graph Fourier Modes

Definition 1 (Generalized Graph Fourier Modes). Given the Hilbert space H G (Q) of graph signals and the variation operator ∆, the set of (∆, Q)-graph Fourier modes is defined as an orthonormal basis of graph signals {u l } l solution to the following sequence of minimization problems, for increasing L ∈ {0, . . . , N -1}:

min u L ∆(u L ) subj. to U H L QU L = I, (1)
where

U L = [u 0 . . . u L ].
This definition is mathematically sound since the search space of each minimization problem is non empty. Note also that this definition relies on only two assumptions: i) Q is an inner product matrix, and ii) ∆ maps graph signals to real values. 6Among the examples of graph variation examples given in Tabs. I and II, those of Tab. I share the important property of being quadratic forms: Definition 2 (Hermitian Positive Semi-Definite Form). The variation operator ∆ is a Hermitian positive semi-definite (HPSD) form if and only if ∆(x) = x H Mx and M is a Hermitian positive semi-definite matrix.

When ∆ is an HPSD form it is then algebraically characterized by the matrix M. In what follows we denote the graph variation operator as M whenever ∆(x) = x H Mx is verified. Examples from the literature of graph variation operators that are HPSD forms are shown on Tab. I, and non-HPSD ones are shown on Tab. II.

The following theorem is of utmost importance as it relates the solution of (1) to a generalized eigenvalue problem when ∆ is an HPSD form: Theorem 1. If ∆ is an HPSD form with HPSD matrix M, then {u l } l is a set of (∆, Q)-graph Fourier modes if and only if the graph signals {u l } l solve the following generalized eigenvalue problems for increasing eigenvalues λ l = ∆(u l ):

Mu l = λ l Qu l , with u l 2 Q = 1. Proof. Observing that x H Mx = Q -1 Mx, x Q , we can show that M Hermitian is equivalent to Q -1 M
being a selfadjoint operator on the Hilbert space H G (Q). Therefore, using the spectral theorem, there exists an orthonormal basis

{u l } l of H G (Q) of eigenvectors of Q -1 M, such that: ∆(u l ) = u l H Mu l = u l H QQ -1 Mu l = λ l u l H Qu l = λ l .
The relation above between the variation of the graph Fourier modes and the eigenvalues of M also shows that, just as approaches of the literature based on the Laplacian [START_REF] Shuman | The Emerging Field of Signal Processing on Graphs: Extending High-Dimensional Data Analysis to Networks and Other Irregular Domains[END_REF] or on the adjacency matrix [START_REF] Sandryhaila | Discrete Signal Processing on Graphs[END_REF], the set of graph Fourier modes is not unique when there are multiple eigenvalues:

If ∆(u l ) = ∆(u k ), then if v l = αu l + (1 -α)u k and v k = (1 -α)u l + αu k , we have ∆(u l ) = ∆(v l ) = ∆(v k ) and v l , v k Q = 0, for any α ∈ [0, 1]
. This remark will be important in Sec. III, as it is desirable for the processing of graph signals to be independent of the particular choice of (∆, Q)-graph Fourier modes.

D. Example with the Ring Graph: Influence of Q

To give intuitions on the influence of Q on the (∆, Q)graph Fourier modes, we study the example of a ring graph with the Laplacian matrix as graph variation operator matrix: M = L. In Fig. 1, we show the (L, Q)-graph Fourier modes of a ring with 8 vertices for three choices of Q:

• Q = diag(0.1, 1, . . . , 1) with a less important vertex 1, • Q = I, i.e., the combinatorial Laplacian GFT, • Q = diag(10, 1, . . . , 1) with a more important vertex 1. Several observations can be made. First of all, any (L, I)graph Fourier mode which is zero on vertex 1 (i.e., u 1 , u 3 , u 5 ) is also an (L, Q)-graph Fourier mode with the same graph variation. Intuitively, this vertex having zero value means that it has no influence on the mode (here, Qu l = u l ), hence its importance does not affect the mode.

For the remaining modes, we make several observations. First of all, we consider the spectral representation of a highly localized signal on vertex 1, such as δ 1 pictured in the last column of Fig. 1. While Q = I involves many spectral components of roughly the same power, the other two cases q 1 = 0.1 and q 1 = 10 are distinctive. Indeed, in the first case, we have δ 1 [START_REF] Belkin | Towards a theoretical foundation for Laplacian-based manifold methods[END_REF] that is large (almost an order of magnitude larger), interpreted by a graph Fourier mode [u 7 ] that is close to our localized signal. On the other hand, for q 1 = 10 the two largest Fourier components of δ 1 are δ 1 (0) and δ 1 (1) with [u 1 ] being close to our localized signal. In other words, q i shifted the impulse δ 1 to the higher spectrum when small (unimportant vertex i) or the lower spectrum when large (important vertex i).

These cases give intuitions on the impact of Q, and ultimately on how to choose it.

E. Discussion on M

When ∆ is and HPSD from, we can rewrite the minimization of (1) as a generalized Rayleigh quotient minimization. Indeed, the L th minimization problem is also given by:

min x:U H L Qx=0 ∆ x x Q ,
which is exactly a Rayleigh quotient minimization when ∆ is an HPSD form of HPSD matrix M since:

∆ x x Q = x H Mx x H Qx .
For example, this allows the study of bandpass signals using spectral proxies as in [START_REF] Anis | Efficient Sampling Set Selection for Bandlimited Graph Signals Using Graph Spectral Proxies[END_REF].

Having an HPSD form for ∆ also has two advantages, the first one being a simple solution to (1) as stated in Thm. 1. But more importantly, this solution is efficient 7 to compute through the generalized eigenvalue problem since both Q and M are Hermitian [START_REF] Golub | Matrix Computations[END_REF]Section 8.7]. Therefore, instead of inverting Q and computing eigenvalues and eigenvectors of the non-Hermitian matrix Q -1 M, the two matrices M and Q can be directly used to compute the (M, Q)-graph Fourier modes using their Hermitian properties.

λ 0 = 0 λ 1 = 0.59 λ 2 = 0.74 λ 3 = 2 λ 4 = 2.42 λ 5 = 3.41 λ 6 = 3.8 λ 7 = 21.05 q 1 = 0.1 1 1 1 1 1 1 1 1 10 20 -0.30 0.05 λ l δ 1 (l) λ 0 = 0 λ 1 = 0.59 λ 2 = 0.59 λ 3 = 2 λ 4 = 2 λ 5 = 3.41 λ 6 = 3.41 λ 7 = 4 Q = I 1 1 1 1 1 1 1 1 2 4 -0.5 0.5 λ l δ 1 (l) λ 0 = 0 λ 1 = 0.24 λ 2 = 0.59 λ 3 = 1.31 λ 4 = 2 λ 5 = 2.8 λ 6 = 3.41 λ 7 = 3.85 q 1 = 10 1 1 1 1 1 1 1 1 2 4 -0.5 2 λ l δ 1 (l)
Fig. 1. (L, Q)-graph Fourier modes of the ring graph with N = 8 vertices. Q = diag(q 1 , 1, . . . , 1) with one choice of q 1 per row of the table, for increasing values of q 1 . Modes are not aligned on their index, but when they correspond either exactly (3 rd , 5 th , and 8 th column), or approximately, to highlight the impact of q 1 on the graph Fourier modes. Colors scale from deep blue for -1 to deep red for 1 and color being lighter as values get closer to 0 until white for 0. Note that [u 7 ] 1 = -3.08 for q 1 = 0.1 is out of this range. The last column shows the graph Fourier transform δ 1 of the impulse δ 1 .

In some applications, it may also be desirable to have Q depend on M. In our motivating example the key variations of the constant signal 1 and of the impulses δ i are given by:

∆(1) = ij M ij ∆(δ i ) = M ii .
In the case of 1, variation should be zero, hence ij M ij = 0. However, for δ i , the variations above cannot be directly compared from one impulse to another as those impulses have different energy. Variations of the energy normalized signals yields:

∆ δ i δ i Q = M ii Q ii .
We advocated for a constant (energy normalized) variation. In effect, this leads to a relation between Q and M given by constant M ii /Q ii , i.e., the diagonals of Q and M should be equal, up to a multiplicative constant. For instance, choosing M = L, which leads to ∆(1) = 0, and ∆(δ i ) = d i , our requirement on Q leads to q i = d i , hence Q = D. This is the random walk Laplacian approach described in Sec. IV-B.

F. Relation to [19]

Finally, we look at the recent advances in defining the graph Fourier transform from the literature, and in particular to [START_REF] Sardellitti | On the Graph Fourier Transform for Directed Graphs[END_REF] which aims at defining an orthonormal set of graph Fourier modes for a directed graph with non-negative weights. After defining the graph directed variation as:

GDV(x) := ij w(ij) [x i -x j ] + ,
where [x] + = max(0, x), the graph Fourier modes are computed as a set of orthonormal vectors w.r.t. the dot product that solve:

min {u l } l l GDV(u l ) subj. to U H U = I. (2)
The straightforward generalization of this optimization problem to any graph variation operator ∆ and Q-inner product, is then:

min {u l } l l ∆(u l ) subj. to U H QU = I. (3)
Unfortunately, we cannot use this generalization when ∆ is an HPSD form. Indeed, in that case, the sum in ( 3) is exactly tr(U H MU) = tr(Q -1 M) for any matrix U verifying U H QU = I. For example, U = Q -1/2 solves (3), but under the assumption that Q is diagonal, this leads to a trivial diagonal GFT, modes localized on vertices of the graph. Note that, for any graph variation, a solution to our proposed minimization problem in (1) is also a solution to the generalization in (3): Property 1. If {u l } l is a set of (∆, Q)-graph Fourier modes, then it is a solution to:

min {u l } l l ∆(u L ) subj. to U H QU = I, with U = [u 0 . . . u N-1 ].
For example, the set of (GDV, I)-graph Fourier modes is a solution to (2), hence a set of graph Fourier modes according to [START_REF] Sardellitti | On the Graph Fourier Transform for Directed Graphs[END_REF]. This property is important as it allows, when direct computation through a closed form solution of the graph Fourier modes is not possible, to approximate the (∆, Q)graph Fourier modes by first using the techniques of [START_REF] Sardellitti | On the Graph Fourier Transform for Directed Graphs[END_REF] to obtain the (Q 1/2 ∆ Q -1/2 , I)-graph Fourier modes and then using Remark 1.

Finally, another recent work uses a related optimization function to obtain graph Fourier modes [START_REF] Shafipour | A Digraph Fourier Transform With Spread Frequency Components[END_REF]:

min {u l } l l (f l -GQDV(u l )) 2 subj. to U H U = I, with f l = l-1
N -1 f max and the graph quadratic directed variation defined as:

GQDV(x) := ij w(ji) [x i -x j ] 2 + .
Beyond the use of a squared directed difference

[x i -x j ] 2 + ,
the goal of this optimization is to obtain evenly distributed graph frequencies, and not orthogonal graph Fourier modes of minimally increasing variation of [START_REF] Sardellitti | On the Graph Fourier Transform for Directed Graphs[END_REF] and of our contribution. Note that this alternative approach can implement the constraint U H QU = I to use an alternative definition of graph signal energy. This is however out of the scope of this paper.

G. Contribution 2: Generalized GFT

Given the (∆, Q)-graph Fourier modes in the previous section, and the definition of the inverse graph Fourier transform F -1 = U found in the literature and recalled in Sec. II-A, we can now define the generalized graph Fourier transform. Note that ∆ is not assumed to be an HPSD form anymore.

Definition 3 (Generalized Graph Fourier Transform). Let U be the matrix of (∆, Q)-graph Fourier modes. The (∆, Q)-graph Fourier transform ((∆, Q)-GFT) is then:

F = U H Q,
and its inverse is:

F -1 = U.
The inverse in Def. 3 above is a proper inverse since FF -1 = U H QU = I. In Sec. II-A, we introduced the graph Fourier transform as an orthogonal projection on the graph Fourier modes. Def. 3 is indeed such a projection since:

x(l) = [Fx] l = U H Qx l = u l H Qx = x, u l Q .
Another important remark on Def. 3 concerns the complexity of computing the graph Fourier matrix. Indeed, this matrix is the inverse of the inverse graph Fourier matrix F -1 = U which is directly obtained using the graph Fourier modes. Computation of a matrix inverse is in general costly and subject to approximation errors. However, just as classical GFTs for undirected graphs, we can obtain this inverse without preforming an actual matrix inverse. Indeed, our GFT matrix is obtained through a simple matrix multiplication U H Q = U -1 that uses the orthonormal property of the graph Fourier modes basis. Additionally, when Q is diagonal, this matrix multiplication is extremely simple and easy to compute. One property that is essential in the context of an orthonormal set of graph Fourier modes is Parseval's Identity where the energy in the vertex and spectral domains are equal:

Property 2 (Generalized Parseval's Theorem). The (∆, Q)- GFT is an isometric operator from H G (Q) to H G (I). Mathe- matically: x, y Q = x, y I . Proof. x, y Q = y H Qx = y H F -1 H QF -1 x = y H x.
Finally, Prop. 2 bears similarities with Remark 1. Indeed, in both cases, there is an isometric map from the Hilbert space H G (Q) to either H G (I) or H G (I). However, in Prop. 2, the graph signal x is mapped to the spectral components Fx of x instead of Q 1/2 x such that both cases are distinct. Intuitively, ϕ and Q 1/2 reshape the space of graph signals to account for irregularity, while the graph Fourier matrix F decomposes the graph signal into elementary graph Fourier modes of distinct graph variation.

III. GRAPH FILTERS

In this section we explore the concept of operator on graph signals, and more specifically, operators whose definition is intertwined with the graph Fourier transform. Such a relation enforces a spectral interpretation of the operator and ultimately ensures that the graph structure plays an important role in the output of the operator.

A. Fundamental Matrix of the GFT

Before defining graph filters, we need to introduce what we call the fundamental matrix of the (∆, Q)-GFT given its graph Fourier matrix F and the diagonal matrix of graph frequencies Λ:

Z := F -1 ΛF = UΛU H Q.
Although ∆ does not appear in the definition above, Z does depend on ∆ through the graph Fourier matrix F. Some authors use the term shift for this matrix when Q = I. However, the literature uses very often the Laplacian matrix [START_REF] Shuman | The Emerging Field of Signal Processing on Graphs: Extending High-Dimensional Data Analysis to Networks and Other Irregular Domains[END_REF] as the fundamental matrix, and as a close equivalent to a second order differential operator [START_REF] Hein | Graph Laplacians and their Convergence on Random Neighborhood Graphs[END_REF], it does not qualify as the equivalent of a shift operator. Therefore, we choose the more neutral term of fundamental matrix. Yet, our definition is a generalization of the literature where such a matrix is always diagonalizable in the graph Fourier basis, with graph frequencies as eigenvalues. Tab. III shows several classical choices of fundamental matrices depending on Q and ∆.

Further assuming that ∆ is an HPSD form, we have Z = Q -1 M. As we will see in Sec. IV, this is consistent with the graph signal processing literature. Noticeably, this matrix is also uniquely defined under these conditions, even though the graph Fourier matrix may not be. Moreover, complexity of computing this matrix is negligible when Q is diagonal. A diagonal matrix Q also leads to Z having the same sparsity as M since M and Z differ only for non-zero elements of M. Algorithms whose efficiency is driven by the sparsity of M are generalized without loss of efficiency to Q -1 M = Z. This will be the case for the examples of this paper.

B. Definitions of Graph Filters

We recall here the classical definitions of graph filters and describe how we straightforwardly generalize them using the fundamental matrix Z of the GFT.

First of all, given a graph filter H, we denote by H = FHF -1 the same graph filter in the spectral domain such that Hx = H x. This notation allows to properly study the spectral response of a given graph filter. We can now state the three definitions of graph filters.

The first one directly extends the invariance through time shifting. This is the approach followed by [START_REF] Sandryhaila | Discrete Signal Processing on Graphs[END_REF], replacing the adjacency matrix by the fundamental matrix of the graph: Definition 4 (Graph Filters by Invariance). H is a graph filter if and only if it commutes with the fundamental matrix Z of the GFT:

HZ = ZH.
The second definition extends the convolution theorem for temporal signals, where convolution in the time domain is equivalent to pointwise multiplication in the spectral domain [START_REF] Girault | Signal Processing on Graphs -Contributions to an Emerging Field[END_REF]Sec. 3.2]8 . The following definition is identical to those in the literature, simply replacing existing choices of GFT with one of our proposed GFTs: Definition 5 (Graph Filters by Convolution Theorem). H is a graph filter if and only if there exists a graph signal h such that:

Hx(l) = h(l) x(l).
The third definition is also a consequence of the convolution theorem, but with a different interpretation. Indeed, in the time domain, the Fourier transform s of a signal s is a function of the (signed) frequency. Given that there is a finite number of graph frequencies, any function of the graph frequency can be written as a polynomial of the graph frequency (through polynomial interpolation). We obtain the following definition using the fundamental matrix instead [START_REF] Shuman | Vertexfrequency analysis on graphs[END_REF], [START_REF] Hammond | Wavelets on graphs via spectral graph theory[END_REF]: Definition 6 (Graph Filters by Polynomials). H is a graph filter if and only if it is a polynomial in the fundamental matrix Z of the GFT:

H = k h k Z k .
Interestingly, these definitions are equivalent for a majority of graphs where no two graph frequencies are equal. However, in the converse case of a graph with two equal graph frequencies λ k = λ l , these definitions differ. Indeed, Def. 6 implies that h(k) = h(l), while h(k) and h(l) are arbitrary according to Def. 5. Also, a graph filter according to Def. 4 does not necessarily verify H diagonal, since Hu l = H k,l u k + H l,l u l is a valid graph filter, even with H k,l = 0. Choosing one definition over another is an application-dependent choice driven by the meaning of two graph Fourier modes of equal graph frequencies. For example, if these modes are highly related, then Def. 6 is a natural choice with equal spectral response of the filter, whereas in the opposite case of unrelated modes, Def. 5 allows to account for more flexibility in the design of the filter.

C. Mean Square Error

Mean Square Error (MSE) is classically used to study the error made by a filter when attempting to recover a signal from a noisy input. More precisely, given an observation y = x + n of a signal x with additive random noise n (with zero mean), MSE is defined as:

MSE H (x) := E Hy -x 2 2 = i E ([Hy] i -x i ) 2 . ( 4
)
In other words, this is the mean energy of the error made by the filter. However, the definition of energy used above is the dot product. As stated before, this energy does not account for irregularity of the structure.

In the general case of an irregular graph, we defined in Sec. II-B graph signal energy using the Q-norm. We define in this section the energy of an error by the Q-norm of that error, thus generalizing mean square error into the Q-MSE:

MSE (Q) H (x) := E Hy -x 2 Q .
Using the zero mean assumption, this yields the classical bias/variance formula:

MSE (Q) H (x) := Hx -x 2 Q bias term + E Hn 2 Q noise variance term (5)
The definition of graph Fourier transform introduced in Def. 3 is then a natural fit to study this MSE thanks to Parseval's identity (Prop. 2):

MSE (Q) H (x) = H -I x 2 2 + E H n 2 2 .
Since H is diagonal (or block diagonal if using Def. 4), studying the bias and the variance from a spectral point of view is much simpler. Recalling the meaning of the graph Fourier modes in terms of graph variation, this also allows to quantify the bias and variance for slowly varying components of the signal (lower spectrum) to quickly varying components (higher spectrum) giving an intuitive interpretation to MSE and the bias/variance trade off across graph frequencies.

To have a better intuition on how using the Q-norm allows to account for irreguarity, we now assume that Q is diagonal. Let = Hn be the filtered noise. The noise variance term above becomes:

E 2 Q = i q i E | i | 2
In other words, those vertices with higher q i have a higher impact on the overall noise energy. Using E

2 Q = tr
H Q , we also have:

E Hn 2 Q = tr H H HFΣ n QF -1
, where Σ n = E nn H is the covariance matrix of the noise. We introduce a definition of Q-white noise with a tractable power spectrum:

Definition 7 (Q-White Noise). The graph signal noise n is said Q-White Noise (Q-WN) if and only if it is centered E [n] = 0 and its covariance matrix verifies Σ n = σ 2 Q -1 , for some σ ≥ 0.
There are two important observations to make on this definition. In the vertex domain, if we assume a diagonal Q, this noise has higher power on vertices with smaller q i :

E |n i | 2 = σ 2 qi .
Assuming Q-WN is equivalent to assuming less noise on vertices with higher q i . Therefore, this definition of noise can account for the irregularity of the graph structure through Q.

Second, the importance of this definition is best seen in the spectral domain. Indeed, the spectral covariance matrix verifies:

Σ n = FΣ n F H = σ 2 I.
In other words, a Q-WN has a flat power spectrum. Note that this is true independently of the variation operator ∆ chosen to define the (∆, Q)-GFT matrix F, hence the name of Q-WN (and not (∆, Q)-WN).

The noise variance term in [START_REF] Girault | Graphe de contacts et ondelettes : étude d'une diffusion bactérienne[END_REF] under the assumption of a Q-WN n becomes:

E Hn 2 Q = σ 2 tr H H H .
In other words, it is completely characterized by the spectral response of the filter H. Note that using Def. 4, H is not necessarily Hermitian. Using Def. 5 or Def. 6, H is diagonal and the noise variance term simplifies to:

E Hn 2 Q = σ 2 l h(l) 2 .

IV. CHOOSING THE GRAPH SIGNAL ENERGY MATRIX Q

The goal of this section is twofold. First, we look at the literature on the graph Fourier transform to show how closely it relates to our definition, including the classical Fourier transform for temporal signals. This is summarized in Tab. III. The second goal is to give examples of graph Fourier transforms that use the newly introduced degree of freedom allowed by the introduction of the Q-inner product.

A. The Dot Product: Q = I 1) Temporal Signals: This case derives from classical digital signal processing [START_REF] Oppenheim | Discrete-time Signal Processing[END_REF], with x[k] a periodic temporal signal of period T and sampled with sampling frequency N/T (N samples per period). This sampling corresponds to a ring graph with N vertices. In this case, energy is classically defined as a scaled 2 -norm of the vector

[x[0], . . . , x[N -1]] T , i.e., E x = T N x H Ix. DFT modes are eigenvectors of the continuous Laplacian operator ∆ T x = d 2 x
dt 2 , thus corresponding to the variation ∆(x) = ∆ T x, x = x H ∆ T x. Finally, DFT modes are orthogonal w.r.t. the dot product.

2) Combinatorial and Normalized Laplacians: This classical case relies on the computations of the eigenvectors of the combinatorial Laplacian L = D -A or of the normalized Laplacian L = D -1/2 LD -1/2 to define the graph Fourier modes [START_REF] Shuman | The Emerging Field of Signal Processing on Graphs: Extending High-Dimensional Data Analysis to Networks and Other Irregular Domains[END_REF]. These graph Fourier transforms are exactly the (L, I)-GFT and the (L, I)-GFT.

3) Graph Shift: This is a more complex case where the (generalized) eigenvectors of the adjacency matrix are used as graph Fourier modes [START_REF] Sandryhaila | Discrete Signal Processing on Graphs[END_REF]. When the graph is undirected, it can be shown that this case corresponds to the ((I -A (norm) ) 2 , I)-GFT [START_REF] Sandryhaila | Discrete Signal Processing on Graphs: Frequency Analysis[END_REF], with ∆(x) = x -A (norm) x 2 2 = GQV(x), the Graph Quadratic Variation. An alternative definition of smoothness based on the 1 -norm of x -A (norm) x leading to the Graph Total Variation ∆(x) = x -A (norm) x 1 = GTV(x) is also investigated in [START_REF] Sandryhaila | Discrete Signal Processing on Graphs: Frequency Analysis[END_REF]. However, this norm leads to a different graph Fourier transform when used to solve (1), since an 1 -norm promotes sparsity (smaller support in graph Fourier modes) while an 2 -norm promotes smoothness. Finally, when the graph is directed, the resulting GFT of [START_REF] Sandryhaila | Discrete Signal Processing on Graphs[END_REF] has no guarantee that its graph Fourier modes are orthogonal. Note that the same ∆(x) = x -A (norm) x 2 2 can be used in the directed case, leading to M = (I-A (norm) ) H (I-A (norm) ). Solving ( 1) is then equivalent to computing the SVD of (I -A (norm) ) and using its right singular vectors as (M, I)graph Fourier modes. However, these differ from the eigenvectors of (I -A (norm) ) used in [START_REF] Sandryhaila | Discrete Signal Processing on Graphs[END_REF] since those are orthogonal.

4) Directed Variation Approach (∆ = GDV): We presented in Sec. II-F a recently proposed GFT [START_REF] Sardellitti | On the Graph Fourier Transform for Directed Graphs[END_REF]. This approach minimizes the sum of directed variation of an orthogonal set of graph signals (see ( 3)). As explained in Sec. II-F, the (∆, I)-GFT is a directed GFT as described in [START_REF] Sardellitti | On the Graph Fourier Transform for Directed Graphs[END_REF].

B. The Degree-Norm (Random Walk Laplacian

): Q = D
The random walk Laplacian, L RW = D -1 L, is not widely used in the graph signal processing literature, given its lack of symmetry, so that its eigenvectors are not necessarily orthogonal. Therefore, the graph Fourier matrix F is not unitary, hence naively computing it through the matrix inverse F = U -1 is not efficient. Yet, this normalization is successfully used in clustering in a similar manner to the combinatorial and normalized Laplacians [START_REF] Von | A tutorial on spectral clustering[END_REF]. Noticeably, it can be leveraged to compute an optimal embedding of the vertices into a lower dimensionnal space, using the first few eigenvectors [START_REF] Belkin | Laplacian eigenmaps for dimensionality reduction and data representation[END_REF]. In graph signal processing we can cite [START_REF] Sunil | Compact Support Biorthogonal Wavelet Filterbanks for Arbitrary Undirected Graphs[END_REF], [START_REF] Gadde | Bilateral filter: Graph spectral interpretation and extensions[END_REF] as examples of applications of the random walk Laplacian. Another example, in image processing, is [START_REF] Liu | Random Walk Graph Laplacian-Based Smoothness Prior for Soft Decoding of JPEG Images[END_REF] where the authors use the random walk Laplacian as smoothness prior for soft decoding of JPEG images through x H LD -1 Lx = L RW x 2 D . Our framework actually allows for better insights on this case from a graph signal processing perspective 9 . Indeed, considering the inner product with Q = D, we obtain the (L, D)-GFT, whose fundamental matrix is Z = D -1 L = L RW . In [START_REF] Liu | Random Walk Graph Laplacian-Based Smoothness Prior for Soft Decoding of JPEG Images[END_REF], this leads to minimizing L RW x 2 D = Zx 2 D which is equivalent to minimizing the energy in the higher spectrum since Z is a high pass filter. This GFT is orthonormal w.r.t. the D-inner product, leading to a graph signal energy definition based on the degree-norm:

E x = x 2 D = x H Dx.
As stated in Remark 1, this normalization is related to the normalized Laplacian through the relation L RW = D -1/2 LD 1/2 , such that L RW and L share the same eigenvalues and their eigenvectors are related: If x is an eigenvector of L RW , ϕ(x) is an eigenvector of L.

By combining the variation of the combinatorial Laplacian with the normalization of the normalized Laplacian, the random walk Laplacian achieves properties that are desirable for

Q Orthon. ∆(1) = 0 ∆(δ i ) = cst Z [15] Comb. Lapl. Non-neg. L = D -A I L [15] Norm. Lapl. Non-neg. L = D -1 2 LD -1 2 I L [14] Graph Shift Complex (I -A norm ) H × I A (I -A norm ) [19] Graph Cut Non-neg. GDV I Approx. n.a. RW Lapl. Non-neg. L RW = D -1 L I L RW L D
a sensor network: ∆ (δ i / δ i D ) = di di = 1 since δ i D = d i , and L1 = 0 = 0 • D1 i.e., the graph frequencies are normalized, the constant signals have zero variation, and all impulses have the same variation, while having different energies 10 .

This case is actually justified in the context of manifold sampling in the extension [START_REF] Hein | Graph Laplacians and their Convergence on Random Neighborhood Graphs[END_REF] of [START_REF] Belkin | Towards a theoretical foundation for Laplacian-based manifold methods[END_REF]. The authors show that under some conditions on the weights of the edges of the similarity graph, the random walk Laplacian is essentially the continuous Laplacian, without additive or multiplicative term, even if the probability density function of the samples is not uniform.

C. Bilateral Filters

: Q = I + D
Bilateral filters are used in image processing to denoise images while retaining clarity on edges in the image [START_REF] Tomasi | Bilateral Filtering for Gray and Color Images[END_REF]:

y i = 1 1 + d i   x i + j w(ij)x j   with w(ij) = exp - p i -p j 2 2σ 2 d exp - I(i) -I(j) 2 2σ 2 i
where p i is the position of pixel i, I(i) is its intensity, and σ d and σ i are parameters. Intuitively, weights are smaller when pixels are either far (first Gaussian kernel) or of different intensities (second Gaussian kernel). This second case corresponds to image edges. We can rewrite this filtering operation in matrix form as:

y = (I + D) -1 (I + A)x = I -(I + D) -1 Lx.
In other words, using Z = (I + D) -1 L, we obtain that bilateral filtering is the polynomial graph filter I -Z. This corresponds exactly to a graph filter with the (L, I + D)-GFT.

Moreover, given a noisy observation y = x + n of the noiseless image x, noise on the output of this filter is given by (I-Z)y-x = Zx+(I-Z)n. This noise can be studied using the (I + D)-MSE introduced in Sec. III-C. Indeed, we can experimentally observe that pixels do not contribute equally to the overall error, with pixels on the edges (lower degree) being less filtered than pixels in smooth regions (higher degree). experimentally, we observe that (I -Z)n is an (I + D)-WN, thus validating the use of the (I + D)-MSE. This approach of quantifying noise is coherent with human perception of noise as the human eye is more sensitive to small changes in smooth 10 Impulse energies show here the irregularity of the graph structure, and as such are naturally not constant.

areas. We will develop the study of bilateral filters with respect to our setting in a future communication.

Finally, note that the expression of bilateral filtering shown here is different than in [START_REF] Gadde | Bilateral filter: Graph spectral interpretation and extensions[END_REF]:

y i = 1 d i j w(ij)x j .
This expression is equivalent to considering the filter I -Z with the (L, D)-GFT, and is not the original approach of [START_REF] Tomasi | Bilateral Filtering for Gray and Color Images[END_REF].

D. Voronoi Cell Areas: Q = C
In the motivating example of a sensor network, we wished to find a graph Fourier transform that is not biased by the particular sampling being performed. Considering the energy of the signal being measured, this energy should not vary with the sampling. In the continuous domain, we observe that energy of a continuous signal s is defined as:

E s := | s(x)| 2 dx.
In the definition above, dx acts as an elementary volume, and the integral can be interpreted as a sum over elementary volumes of the typical value | s(x)| within that volume times the size of the volume dx. The discrete version of this integral is therefore:

E s ≈ i | s(x i )| 2 vol(i) = i |s i | 2 vol(i),
with s the sampled signal on the points {x i } i . Given a particular sampling, the question is then what is a good value of vol(i)? A simple and intuitive approximation is to approximate this volume with the volume of the subspace of points whose closest sample is x i . This subspace is exactly the Voronoi cell of i, and the volume is therefore the area of the cell in 2D. Let c i be this area, and C = diag(c 1 , . . . , c N ). We obtain:

E s ≈ i |s i | 2 c i = s H Cs = s 2 C .
Intuitively, this corresponds to interpolating the sampled signal into a piecewise constant signal for which the signal values are equal within each Voronoi cell, and then computing the continuous energy of this interpolated signal. Other interpolation schemes could be considered. For example, if we assume a weighted linear interpolation scheme, then we obtain the interpolated signal s(x) = i f (x, x i )s i and its energy is:

E s = | s(x)| 2 dx = i,j f (x, x i ) H f (x, x j )dx s H i s j ,
and we have As soon as there is one location x whose interpolation involves two samples i,j, then this interpolation scheme corresponds a non-diagonal matrix Q. However, as we will show in Sec. V-B, the approach based on the Voronoi cell areas gives already good results compared to the state of the art of graph signal processing for sensor networks.

E s = s 2 Q with q ij = f (x, x i ) H f (x, x j )dx.

V. EXPERIMENTS

The choice of a graph Fourier transform, by selecting a graph variation operator ∆ and a matrix Q such as those discussed in Sec. IV, is application dependent. In particular, meaningful properties for the graph Fourier transform can be drastically different. For example, in a clustering context, the practitioner is interested in extracting the structure, i.e., in identifying tight groups of vertices. This is in constrast to sensor networks where the particular arrangement of sensors should have as little influence over the results as possible, i.e., tight groups of stations should not be a source of bias (e.g., weather stations).

For this reason, there is no unique way of studying how good a particular graph Fourier transform is. In this section, we use these two applications, i.e., clustering and sensor networks, to show how our framework can be leveraged to achieve their application-specific goals.

All experiments were performed using the GraSP toolbox for Matlab [START_REF] Girault | Grasp: A matlab toolbox for graph signal processing[END_REF].

A. Clustering

In this section, we study the problem of clustering defined as grouping similar objects into groups of object that are similar while being dissimilar to objects in other groups. In a graph setting, we are interested in grouping vertices of a graph such that there are many edges within each group, while groups are connected by few edges. Our goal is not to provide a new approach to perform clustering, but rather use this problem as a showcase for how using a well defined matrix Q can help achieve the goals of a target application.

In the context of clustering, spectral clustering extracts groups using the spectral properties of the graph [START_REF] Von | A tutorial on spectral clustering[END_REF]. More precisely, using C-means (k-means with k = C) on the first C graph Fourier modes yields interesting partitions. In [START_REF] Von | A tutorial on spectral clustering[END_REF], the author interprets this approach using graph cuts, for (L, I)-GFT (combinatorial Laplacian), (L, I)-GFT (normalized Laplacian) and (L, D)-GFT (random walk Laplacian). We extend this interpretation here to any variation operator ∆ and any diagonal innner product matrix Q.

For each cluster c, let V c ⊂ V be the subset of its vertices. Then the set {V c } c of all these subsets is a partition of the set of vertices V. Let the normalized cut of {V c } c associated to the (∆, Q)-GFT be:

(∆, Q) -Ncut ({V c } c ) := c ∆ h (Q) c with h (Q) c
the normalized indicator function of cluster V c verifying:

h (Q) c := 1 1 Vc Q 1 Vc
and 1 Vc the indicator function of cluster c. This extends the normalized cut interpretation of [START_REF] Von | A tutorial on spectral clustering[END_REF]. Using the notation

H = [h (Q) 1 • • • h (Q) C ] T
, and since Q is diagonal, we obtain the orthonormality property H H QH = I. Finding a partition of vertices minimizing the normalized cut (∆, Q) -Ncut ({V c } c ) is therefore equivalent to finding orthonormal normalized indicator functions of minimal variation.

Spectral clustering is performed by first relaxing the constraint that the graph signals h

(Q) c
are normalized indicator functions. Using Prop. 1, the first C (∆, Q)-graph Fourier modes are solutions to this relaxed optimization problem. The final partition is then obtained through C-means on the spectral features, where the feature vector of vertex i is

[u 0 ] i , . . . , [u C-1 ] i
T . The accuracy of spectral clustering is therefore driven by the choice of variation operator ∆ and inner product matrix Q. To illustrate this, and the statement that choosing Q and ∆ is application-dependent, we look at several choices in the context of a 2-class clustering problem with skewed clusters. Combined with an irregular distribution of inputs within each cluster, we are in the context of irregularity where our framework thrives through its flexibility.

The resulting graph is shown on Fig. 2, with a sparse cluster on the left (30 inputs), and a dense one on the right (300 inputs). Each cluster is defined by a 2D Gaussian distribution, with an overlapping support, and samples are drawn from these two distributions. We already observe that inputs are irregularly distributed, and some inputs are almost isolated. These almost isolated inputs are important as they represent outliers. We will see that correctly choosing Q alleviates the influence of those outliers. The graph of inputs is built in a K-nearest neighbor fashion (K = 10) with weights chosen using a Gaussian kernel of the Euclidean distance (σ = 0.4).

λ l h (Q) blue (λ l ) 0 0.5 1 1.5 2 0 0.1 λ l MSE (Q) H l (b) (L, I) -1 0 1 -1 0 1 u 0 u 1 s d ⊥ d ⊥ s (c) (L, I), feat. norm.
λ l h (Q) blue (λ l ) 0 0.5 1 1.5 2 0 0.1 λ l MSE (Q) H l (d) (L, D) 0.1 -0.3 0 0.1 u 0 u 1 s d ⊥ d 0 1k 2k 3k 0 0.35 λ l h (Q) blue (λ l ) 0 1k 2k 3k 0 0.1 λ l MSE (Q) H l (e) (L, C) 0 0.2 0.4 0 0.2 0.4 u 0 u 1 s d ⊥ d ⊥ s 0 0.4 0.8 0 0.5 λ l h (Q) blue (λ l ) 0 0.4 0.8 0 1 λ l MSE (Q) H l (f) ( x-A (norm) x 1 , I)
Fig. 3. Several spectral clustering results on the graph of Fig. 2. The smallest cluster is in dark red, and the largest cluster is in light blue. The second row of plots shows the input in the feature space, i.e., to each vertex i is associated a point [u 0 ] i , [u 1 ] i T in the spectral feature space. The set s (resp. d ) corresponds to the vertices in the sparse (resp. dense) cluster and correctly clustered. The set ⊥s (resp. ⊥ d ) corresponds to vertices incorrectly clustered in the sparse (resp. dense) cluster. Colors in both rows of figures match for consistency. The third row shows the (M, Q)-GFT of the corresponding normalized indicator function h (Q) blue (the first Fourier component is not shown as it is quite large, except for (f)). The fourth row shows the Q-MSE of the ideal low-pass filter H l with varying cut-off frequency λ l applied to this normalized indicator function. Generalized GFT: Fig. 3 shows the results of spectral clustering on this graph, using several GFTs (first two rows), the analysis of h (Q) blue in the spectral domain (third row), and an example of Q-MSE of several ideal low pass filters on the same indicator function. First of all, the Random Walk Laplacian case based on the (L, D)-GFT gives the best results, consistent with [START_REF] Von | A tutorial on spectral clustering[END_REF] where the author advocates this approach. [START_REF] Von | A tutorial on spectral clustering[END_REF]Prop. 5] gives an intuition on why this case works best using a random walk interpretation: If a starting vertex is chosen at random from the stationary distribution of the random walk on the graph, then the (L, D) -Ncut is exactly the probability to jump from one cluster to another. Minimizing this cut, is then minimizing transitions between clusters. In the spectral domain, we see that the indicator function has a lot of energy in the first few spectral components, with a rapid decay afterwards. We notice also, a slight increase of energy around graph frequency 1 due to the red vertices of the sparse (red) cluster which are in the middle of the dense (blue) cluster (see Fig. 2(a)). As expected, the Q-MSE decreases when the cut-off frequency increases. But interestingly, if we use the classical I-MSE, it is not decreasing with the cut-off frequency because of the use of the Q-norm instead of the 2 -norm (see Fig. 4(a). The combinatorial Laplacian case, i.e., using the (L, I)-GFT, suffers from the presence of outliers. Indeed, if i is such an outlier, then having it in a separate cluster, i.e., V 1 = {i} and V 2 = V\{i} for C = 2, leads to ∆(h

(I) 1 ) = ∆(h (I) 2 ) = d i .
Since i is isolated, it has extremely low degree, and the resulting cut is small, making it a good clustering according to the (L, I) -Ncut. Not accounting for irregularity of the degree in the definition of the normalized indicator function h (I) c leads therefore to poor clustering in the presence of outliers. A similar behavior is observed for the (L, C)-GFT with a normalization by the Voronoi cell area that is large for isolated vertices, hence an even smaller normalized cut. In the spectral domain, more Fourier components of the normalized indicator function are large, especially around the graph frequency 0. This results in a Q-MSE that is large for our ideal low-pass filters. In other words, more lowpass components are required to correctly approximate the indicator function using a lowpass signal.

Next, is the normalized Laplacian case, i.e., using the (L, I)-GFT. Just as the Random Walk Laplacian case, singleton clusters with outliers are not associated to small cuts since (L, I) -Ncut({i}, V\{i}) = (L, D) -Ncut({i}, V\{i}). However, a careful study of the graph Fourier modes reveals the weakness of this case. Indeed, using the well known relation between the (L, I)-graph Fourier modes {u l } l and the (L, D)graph Fourier modes{v l } l with ϕ(v l ) = u l = D 1/2 v l , this case deviates from the Random Walk Laplacian by introducing a scaling of graph Fourier modes by the degree matrix. Differences are then noticeable in the spectral feature space (bottom plots in Fig. 3). Indeed, the graph Fourier modes verify

[u l ] i = √ d i [v l ] i
, such that low degree vertices, e.g., outliers and vertices in the boundary of a cluster, have spectral features of smaller magnitude:

[u 0 ] i , [u 1 ] i T [v 0 ] i , [v 1 ] i T .
In other words, spectral features of low degree vertices are moved closer to the origin of the feature space, and closer to each other. C-means cluster those vertices together, resulting in the poor separation of clusters we see on Fig. 3(b). In the spectral domain, we clearly see the interest of using the random walk Laplacian approach: the normalized indicator functions are not lowpass. This is further seen on the Q-MSE that decreases continuously with the cutoff frequency without the fast decay we wish for.

The spectral clustering literature is actually not using the raw spectral features to run C-means, but normalize those features by their 2 -norm prior to running C-means [START_REF] Andrew | On spectral clustering: Analysis and an algorithm[END_REF]. This is actually equivalent to the projection of the (L, D)-graph Fourier modes {v l } l on the unit circle:

[ u 0 ] i = cos (θ i ) [ u 1 ] i = sin (θ i ) θ i = atan 1 D [v 1 ] i .
We can use this relation to characterize how spectral features of vertices are modified compared to the (L, D)-GFT approach, and how C-means behaves. Indeed, separation of features (i.e., how far the spectral features of two vertices are) is modified by the transformation above. Looking at the norm of the gradient of these spectral features yields:

--→ grad [ u 0 ] i [ u 1 ] i 2 = 1 D 1 + ( 1 D [v 1 ] i ) 2
Therefore, the gradient is smaller for higher values of Finally, we performed the same specral clustering approach using the (∆, I)-GFT with ∆(x) = x-A (norm) x 1 the graph total variation. Results on Fig. 3(f) show that this approach actually identifies the stronger small clusters rather than larger clusters. More importantly, in this case, the spectral features of both clusters are not separated, such that C-means cannot separate clusters based on the spectral features. Furthermore, it turns out that the clusters given by these spectral features are not the best with an (GTV, I) -Ncut of 8.21, while those given by the combinatorial Laplacian yields a (GTV, I) -Ncut of 6.61. In other words, the relaxation we made that consider arbitrary h (I) c instead of normalized indicator functions leads to large errors. The spectral domain analysis suffers from the difficulty to compute efficiently the graph Fourier modes, since the graph variation is not HPSD. However, the first few components of the normalized indicator function show that this signal is now low pass, with a lot of energy in these components. Q-MSE is not sharply decreasing either and stays very high.

B. Sensor Networks

Here we explore our motivating example of a sensor network, where the goal is the opposite of that in the one above: clusters of vertices should not bias the results. More precisely, we are interested in studying signals measured by sensor networks where the location of the sensors is unrelated to the quantity measured, e.g., temperature readings in a weather station network. In this section, we show that we can obtain a good definition of inner product achieving this goal.

Let G (U) be a graph whose vertices are uniformly sampled on a 2D square plane, and G (NU) a similar graph obtained with a non-uniform sampling distribution (see Fig. 5(e) where the areas of higher density are indicated). Edges of both graphs are again selected using K-nearest neighbors, with K = 10. Weights are chosen as a Gaussian kernel of the Euclidean distance w ij = expd(i, j) 2 /(2σ 2 ) , with σ = 0.3 (empirically chosen to have a good spread of weights). Example graphs with 500 sampled vertices are shown in Figs. 5(a) and 5(e). In this experiment, we are interested in studying how we can remove the influence of the particular arrangement of sensors such as to obtain the same results from one graph realization to another. To that end, we generate 500 graph realizations, each with 1000 vertices, to compute empirical statistics.

The signals we consider are pure sine waves, which allows us to control the energies and frequencies of the signals we use. We experiment with several signals for each frequency and report the worst case statistics over those signals of equal (continuous) frequency. This is done by choosing several values of phase shift for each frequency.

Let s ν (x, y; ϕ) = cos(2πνx i +ϕ) be a continuous horizontal Fourier mode of frequency ν phase shifted by ϕ. For a given graph G generated with any of the two schemes above, we define s ν (i; ϕ, G) = s ν (x i , y i ; ϕ) the sampled signal on the graph. Its energy w.r.t. the Q-norm is then:

E (Q) ν (ϕ; G) := s ν (_; ϕ, G) 2
Q . Note that this energy depends on the graph if Q does, i.e., we have Q(G) which may vary from one graph to another such as Q(G) = D(G). We use the shorter notation Q to keep notations shorter.

To get better insights on the influence of the sensor sampling, we are interested in statistical quantities of the graph signal energy. First, we study the empirical mean:

µ (Q) (ν; ϕ) := E (Q) ν (ϕ; G) G = 1 N G N G n=1 E (Q) ν (ϕ; G n ),
with N G the number of sampling realizations. It is interesting to study how the mean varies depending on ϕ and ν in order to check whether the graph signal energy remains constant, given that the energy of the continuous signals being sampled is constant. The first thing we observe is that the averaging over ϕ, i.e., over many signals of equal continuous frequency, yields the same average mean for all continuous frequencies. However, µ (Q) (ν; ϕ) does depend on ϕ. To show this, we use the maximum absolute deviation of the normalized mean μ

(Q) (ν; ϕ) := µ (Q) (ν; ϕ)/ µ (Q) (ν; ϕ) ϕ : m (Q) (ν) := max ϕ 1 -μ(Q) (ν; ϕ) .
Use of the normalized mean is necessary to be able to compare the various choices of Q since they can yield quite different average mean µ (Q) (ν; ϕ) ϕ .

However, the mean is only characterizing the bias of the graph signal energy approximation to the continuous signal energy. We also need to characterize the variance of this estimator. To that end, we consider the empirical standard deviation of the graph signal energy:

σ (Q) (ν; ϕ) := 1 N G -1 N G n=1 E (Q) ν (ϕ; G n ) -µ (Q) (ν; ϕ) 2 .
This quantity shows how much the actual sampling being performed influences the signal energy estimator. Since the mean energy is influenced by the choice of Q, we need to normalize the standard deviation to be able to compare results between various choices of Q. This yields the coefficient of variation, and we report its maximum over all signals of equal (continuous) frequency:

CV (Q) (ν) := max ϕ σ (Q) (ν; ϕ) µ (Q) (ν; ϕ) ,
from which we can study the variance of the graph signal energy estimator depending on the continuous frequency. We experiment here with three choices for Q. The first one a standard GFT from the literature, corresponding to the 2norm: Q = I. The second one corresponds to the random-walk Laplacian with Q = D. The third one is our novel approach based on the Voronoi cell area inner product Q = C (see Sec. IV-D).

Results for G (U) and G (NU) are given in Fig. 5. The coefficient of variation shows here the strong advantage of using the Voronoi cell area inner product: CV (C) (ν) is very small in the lower continuous spectrum, and almost zero (see Figs. 5(c) and 5(g)). In other words, for signals that do not vary too quickly compared to the sampling density, the Cnorm gives an estimated energy with small variance, removing the influence of the particular arrangement of sensors. This is true for both the uniform distribution and the non-uniform distribution. Both the dot product Q = I and the degree norm Q = D have larger variance here.

Finally, the maximum absolute deviation on Figs. 5(d) and 5(h) shows again the advantage of the C-norm. Indeed, considering several signals of equal continuous frequencies, these should have equal average energy. While this is the case for both the dot product and the C-norm in the lower spectrum for a uniform sampling, using a non-uniform sampling yields a strong deviation for the dot product. The C-norm appears again as a good approximation to the continuous energy with small deviation between signals of equal frequencies, and between samplings.

Remark 2. In [START_REF] Hein | Graph Laplacians and their Convergence on Random Neighborhood Graphs[END_REF], [START_REF] Coifman | Diffusion maps[END_REF], the authors advocate for the use of the random walk Laplacian as a good approximation for the continuous (Laplace-Beltrami) operator of a manifold. We showed here that this case is not ideal when it comes to variance of the signal energy estimator, but in those communications, the authors actually normalize the weights of the graph using the degree matrix prior to computing the random walk Laplacian, thus working on a different graph. If A is the adjacency prior its normalization, and D the associated degree matrix, then A = D -1 A D -1 . This normalization is important when we look at the degree. Indeed, without it, we see on Fig. 5(a)-(b) and 5(e)-(f) that degree and Voronoi cell area evolve in opposed directions: large degrees correspond to small areas. Pre-normalizing by the degree corrects this. Using this normalization leads to better results with respect to the energy (omitted here to save space), however, Voronoi cells still yield the best results.

These results show that it is possible to better account for irregularity in the definition of the energy of a graph signal: The Voronoi cell area norm yields very good results when it comes to analysing a signal lying in a Euclidean space (or on a manifold) independently of the sampling performed.

VI. CONCLUSION

We showed that it is possible to define an orthonormal graph Fourier transform with respect to any inner product. This allows in particular to finely account for an irregular graph structure. Defining the associated graph filters is then straightforward using the fundamental matrix of the graph. Under some conditions, this led to a fundamental matrix that is easy to compute and efficient to use. We also showed that the literature on graph signal processing can be interpreted with this graph Fourier transform, often with the dot product as inner product on graph signals. Finally, we showed that we are able to obtain promising results for sensor networks using the Voronoi cell areas inner product.

This work calls for many extensions, and giving a complete list of them would be too long for this conclusion. At this time, we are working on the sensor network case to obtain graph signal energies with even less variance, especially once the signals are filtered. We are also working on studying bilateral filters, and image smoothers in general [START_REF] Milanfar | A Tour of Modern Image Filtering: New Insights and Methods, Both Practical and Theoretical[END_REF], and extending them to graph signals with this new setting. We explore also an extension of the study of random graph signals and stationarity with respect to this new graph Fourier transform. Finally, many communications on graph signal processing can be reinterpreted using alternative inner product that will give more intuitions on the impact of irregularity.

ACKNOWLEDGEMENT

The authors would like to thank the anonymous reviewers for their careful reading and constructive comments that helped improve the quality of this communication.

APPENDIX

A. Low Pass Filtering of an Indicator Function

We now show an example to compare the quality of different low pass approximations to the normalized cluster indicator function, under different GFT definitions. This is motivated by the fact that spectral clustering techniques make use of the low frequencies of the graph spectrum. Fig. 6 shows these approximations in the vertex domain. Note that these correspond to the Q-MSE plots of Fig. 3.

Comparing results for l = 1 and l = 5 we note that for both L and L, the low pass approximation is good, and significantly better than that achieved based on the GFT corresponding to graph total variation. However, the (L, I)-GFT suffers from isolated vertices (l = 1 case on Fig. 6), and the approximation needs more spectral components for an output closer to the indicator function (l = 5 case), while the (L, I)-GFT is clearly biased with the degree (smaller amplitudes on the vertices of the cluster boundaries).

This confirms that the indicator function approximation (and thus spectral clustering performance) is better for the (L, D)-GFT, which does not have the bias towards isolated vertices shown by the (L, I)-GFT (impulses on those isolated vertices are associated to lowpass signal, while they are associated to the graph frequency 1 for the (L, D)-GFT), and without the bias of the first Fourier mode towards vertex degrees of the (L, I)-GFT. blue . H l keeps only the first l + 1 spectral components (i.e., up to the graph frequency λ l ). Note that H l depends on the (∆, Q)-GFT.

Fig. 2 .

 2 Fig. 2. Two cluster dataset input with a sparse cluster (left) and a dense cluster (right). (a) Ground truth clusters. (b) Vertex degrees from the 10nearest neighbor graph Gaussian weights.

  Fig.3. Several spectral clustering results on the graph of Fig.2. The smallest cluster is in dark red, and the largest cluster is in light blue. The second row of plots shows the input in the feature space, i.e., to each vertex i is associated a point[u 0 ] i , [u 1 ] iT in the spectral feature space. The set s (resp. d ) corresponds to the vertices in the sparse (resp. dense) cluster and correctly clustered. The set ⊥s (resp. ⊥ d ) corresponds to vertices incorrectly clustered in the sparse (resp. dense) cluster. Colors in both rows of figures match for consistency. The third row shows the (M, Q)-GFT of the corresponding normalized indicator function h(Q) blue (the first Fourier component is not shown as it is quite large, except for (f)). The fourth row shows the Q-MSE of the ideal low-pass filter H l with varying cut-off frequency λ l applied to this normalized indicator function. Generalized GFT: (a): (L, I), (b) (L, I), (c) (L, I) with normalized feature vectors, (d) (L, D), (e) (L, C), and (f) (∆, I) with ∆(x) = x -A (norm) x 1 .

Fig. 4 .

 4 Fig. 4. I-MSE of the ideal low-pass filter H l applied to the normalized indicator function h (Q) blue corresponding to (a) Fig. 3(d) and (b) to Fig. 3(e).

  |[v 1 ] i |: This spectral feature normalization brings features closer for higher values of |[v 1 ] i |. As features of the dense and sparse clusters are identified by the magnitude of |[v 1 ] i |, this normalization actually brings features of the two clusters closer. This yields less separation in the feature space and poor clustering output from C-means (see Fig. 3(c)).

Fig. 5 .

 5 Fig. 5. Study of the graph energy of single cosine continuous signal sν (x, y; ϕ) = cos(2πνx i + ϕ) sampled uniformly (a)-(d) and non-uniformly (e)-(h) as sν (i; ϕ) = sν (x i , y i ; ϕ). Example sampling with 500 samples (vertices) showing vertex degree ((a) and (e)), and Voronoi cell area ((b) and (f)) with colors. (c)-(d) and (g)-(h) 1000 samples (vertices), with results averaged over 500 sampling (graph) realizations. (c) and (g) coefficient of variation CV (Q) (ν) of the signal energy depending on the continuous frequency ν. (d) and (h) maximum absolute deviation m (Q) (ν) of the normalized mean signal energy depending on the continuous frequency ν.

Fig. 6 .

 6 Fig.6. Output of the ideal low pass filter H l applied to the normalized indicator function h(Q) blue . H l keeps only the first l + 1 spectral components (i.e., up to the graph frequency λ l ). Note that H l depends on the (∆, Q)-GFT.

TABLE I .

 I Hermitian positive semidefinite (HPSD) graph variation operators.

TABLE II .

 II Non-HPSD graph variation operators. Here, GTV is the Graph Total Variation (see Sec. IV-A3), and GDV is the graph directed variation (see Sec. II-F).

TABLE III .

 III State of the art of the graph Fourier transform.

	Ref	Name	Directed	Weights	∆

TABLE IV .

 IV F 1 score of the sparse (left) cluster on Fig.2for the clustering results of Fig.3.

	Approach	Accuracy	Sparse F 1 score
	(L, I)	91.21%	6.45%
	(L, I) (L, I), feat. norm.	93.64% 46.67%	46.15% 4.35%
	(L, D)	96.67%	77.55%
	(L, C)	91.21%	6.45%
	( x -A (norm) x 1 , I)	86.06%	4.17%

Note that studies of stationarity of graph signals, such as[START_REF] Girault | Stationary Graph Signals using an Isometric Graph Translation[END_REF], are concerned with the variability across multiple observations of signals on the same graph. Instead here we focus on how the choice of different graphs (i.e., different vertex/sensor positions placed in space) affects the spectral representation of the corresponding graph signals (e.g., different discrete samples of the same continuous domain signal).

Although vertices are indexed, the analysis should be independent of the specific indexing chosen, and any other indexing shall lead to the same outcome for any vertex for a sound GSP setting.

Throughout this paper, bold uppercase characters represent matrices and bold lowercase characters represent column vectors.

Note that A ij corresponds here to the edge from i to j, opposite of the convention used in[START_REF] Sandryhaila | Discrete Signal Processing on Graphs[END_REF]. This allows for better clarity at the expense of the interpretation of the filtering operation Ax. Here, (Ax) i = j A ij x j is the diffusion along the edges in reverse direction, from j to i.

We also used before λ l = ∆(u l ) when ∆ is a quadratic variation[START_REF] Girault | Translation on Graphs: An Isometric Shift Operator[END_REF]. In this paper, however, the difference between the two does not affect the results.

Q may depend on ∆, but no such relation is required in Def. 1. An example is given at the end of this section where q i = ∆(δ i ) is desirable.

Efficiency comes here from the cost of computing generalized eigenvalues and eigenvectors compared to the when M is not HPSD or Q is not HPD.

In[START_REF] Shuman | Vertexfrequency analysis on graphs[END_REF], the authors also define convolutions as pointwise multiplications in the spectral domain. However, the notation used for the GFT is unclear since x(λ l ) (instead of x(l)) can be interpreted as the graph signal x having equal spectral components for equal graph frequencies. Such a requirement actually falls within the setting of Def. 6.

The property that the eigenvectors of the random walk Laplacian are orthogonal w.r.t. the Q-norm is well known in the literature. We are however the first to make the connection with a properly defined graph Fourier transform.