Fei Hua
email: fei.hua@oca.eu

Member, IEEE Roula Nassif
email: roula.nassif@epfl.ch

Haiyan Wang
email: hywang@sust.edu.cn

Fellow, IEEE Ali H Sayed
email: ali.sayed@epfl.ch

C Richard
email: cedric.richard@unice.fr

Online Distributed Learning over Graphs with Multitask Graph-Filter Models

Keywords: Graph signal processing, graph filter, diffusion LMS, node-varying graph filter, clustering. I

In this work, we are interested in adaptive and distributed estimation of graph filters from streaming data. We formulate this problem as a consensus estimation problem over graphs, which can be addressed with diffusion LMS strategies. Most popular graph-shift operators such as those based on the graph Laplacian matrix, or the adjacency matrix, are not energy preserving. This may result in an ill-conditioned estimation problem, and reduce the convergence speed of the distributed algorithms. To address this issue and improve the transient performance, we introduce a preconditioned graph diffusion LMS algorithm. We also propose a computationally efficient version of this algorithm by approximating the Hessian matrix with local information. Performance analyses in the mean and mean-square sense are provided. Finally, we consider a more general problem where the filter coefficients to estimate may vary over the graph. To avoid a large estimation bias, we introduce an unsupervised clustering method for splitting the global estimation problem into local ones. Numerical results show the effectiveness of the proposed algorithms and validate the theoretical results.

In order to cope with graph signals, GSP relies on two ingredients: the graph shift operator (GSO) on one hand, which accounts for the topology of the graph, and the graph Fourier transform (GFT) on the other hand, which allows to represent graph signals in the graph frequency domain. Built upon the definition of the GFT, graph filters play a central role in processing graph signal spectra by selectively amplifying or attenuating frequency components. Various architectures of graph filters have been proposed in the literature, including finite impulse response (FIR) [START_REF] Shuman | The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains[END_REF], [START_REF] Sandryhaila | Discrete signal processing on graphs[END_REF] and infinite impulse response (IIR) [START_REF] Shi | Infinite impulse response graph filters in wireless sensor networks[END_REF], [START_REF] Liu | Filter design for autoregressive moving average graph filters[END_REF] filters. From a perspective of scalability, and considering energy constraints and band-limited communications that may be encountered in large networks of distributed nodes such as sensor networks, significant efforts have been made recently to derive distributed graph filters. These filtering procedures allow each node to exchange only local information with its neighboring nodes [START_REF] Loukas | Distributed autoregressive moving average graph filters[END_REF]- [START_REF] Shuman | Distributed signal processing via chebyshev polynomial approximation[END_REF].

Much of the GSP literature has focused on static graph signals, that is, signals that need not evolve with time. However, a wide spectrum of network-structured problems requires adaptation to time-varying dynamics. Sensor networks, social networks, vehicular networks, communication networks, and power grids are some typical examples. Prior to the more recent GSP literature, many earlier works on adaptive networks have addressed problems dealing with this challenge by developing processing strategies that are well-suited to data streaming into graphs; see, e.g., [START_REF] Sayed | Diffusion strategies for adaptation and learning over networks: an examination of distributed strategies and network behavior[END_REF]- [START_REF]Diffusion adaptation over networks[END_REF]. Several diffusion strategies have been introduced, and their performance studied in various situations, such as diffusion LMS [START_REF] Lopes | Diffusion least-mean squares over adaptive networks: Formulation and performance analysis[END_REF], RLS [START_REF] Cattivelli | Diffusion recursive leastsquares for distributed estimation over adaptive networks[END_REF], and APA [START_REF] Li | Distributed adaptive estimation based on the APA algorithm over diffusion networks with changing topology[END_REF]. By referring to the problem of estimating an optimal parameter vector at a node as a "task", and depending on the relations between the parameter vectors across the entire network, adaptive networks can be divided into single or multitask networks. In single-task networks, all nodes estimate the same parameter vector. Typical works include [START_REF] Sayed | Diffusion strategies for adaptation and learning over networks: an examination of distributed strategies and network behavior[END_REF]- [START_REF]Diffusion adaptation over networks[END_REF]. With multitask networks, multiple but related parameter vectors are inferred simultaneously in a cooperative manner, so as to improve the estimation accuracy by using the similarities between tasks [START_REF] Chen | Diffusion LMS over multitask networks[END_REF]- [START_REF]Proximal multitask learning over networks with sparsity-inducing coregularization[END_REF].

Related works: In this work, we are interested in online learning of linear graph models for representing streaming graph signals in a distributed manner. We focus on diffusion strategies because they are scalable, robust, and enable network adaptation and learning. Recently, some research works have considered time-varying graph signals. An adaptive graph signal reconstruction algorithm based on the LMS is proposed in [START_REF] Lorenzo | Adaptive least mean squares estimation of graph signals[END_REF] but it operates in a centralized manner. In [START_REF] Mei | Signal processing on graphs: Causal modeling of unstructured data[END_REF], the authors focus on estimating a network structure capturing the dependencies among time series in the form of a possibly directed, weighted adjacency matrix. A causal autoregressive process is introduced in the time series to capture the intuition that information travels over the network at some fixed speed. In [START_REF] Isufi | Forecasting time series with VARMA recursions on graphs[END_REF], vector autoregressive (VAR) and vector autoregressive moving average (VARMA) models are proposed for predicting time-varying processes on graphs. Joint time-vertex stationarity is introduced for time-varying graph signals in [START_REF] Loukas | Stationary time-vertex signal processing[END_REF], [START_REF] Loukas | Predicting the evolution of stationary graph signals[END_REF], and a joint time-vertex harmonic analysis for graph signals is proposed in [START_REF] Grassi | A time-vertex signal processing framework: Scalable processing and meaningful representations for time-series on graphs[END_REF]. It is shown that joint stationarity facilitates estimation or recovery tasks when compared to purely time or graph methods.

All the aforementioned works focus on centralized solutions whereas distributed algorithms may be more appropriate within the context of big data applications. In [START_REF] Isufi | 2-dimensional finite response graph-temporal filters[END_REF], [START_REF] Isufi | Filtering random graph processes over random time-varying graphs[END_REF], the behavior of some distributed graph filters on time-varying graph signals is studied. Considering graph signal sampling and reconstruction, several distributed algorithms have been proposed to track time-varying band-limited graph signals, e.g., LMS-based algorithms in [START_REF] Lorenzo | Distributed adaptive learning of graph signals[END_REF], RLS-based methods in [START_REF] Lorenzo | Adaptive graph signal processing: Algorithms and optimal sampling strategies[END_REF], Kalman-based methods in [START_REF] Isufi | Observing and tracking bandlimited graph processes[END_REF], and kernel-based algorithms in [START_REF] Romero | Kernel-based reconstruction of space-time functions on dynamic graphs[END_REF]. In [START_REF] Qiu | Time-varying graph signal reconstruction[END_REF], the authors are interested in timevarying graph signals with temporal smoothness prior. They devise distributed gradient descent algorithms to reconstruct the signals. Most of these works assume that the graph signals are band-limited. Another limitation is that the graph Fourier decomposition (eigenvectors) is needed beforehand, which is impractical for large networks.

Contributions: Recently, several works successfully applied adaptive algorithms to graph signals. In [START_REF] Lorenzo | Adaptive least mean squares estimation of graph signals[END_REF], [START_REF] Lorenzo | Distributed adaptive learning of graph signals[END_REF], [START_REF] Lorenzo | Adaptive graph signal processing: Algorithms and optimal sampling strategies[END_REF] for instance, graph Fourier coefficients are learned from streaming graph signals under band-limited assumption to perform adaptive reconstruction and tracking of time-varying graph signals from partial observations. In this work, we are interested in online distributed learning of linear graph models without assumption of band-limited processes. We use graph filter models in the time-vertex domain where there is no need to decompose the graph shift operator. The formulated optimization problem relies on minimizing a global cost consisting of the aggregate sum of individual costs at each vertex. To address this problem, we blend concepts from adaptive networks [START_REF] Sayed | Adaptation, learning, and optimization over networks[END_REF] and graph signal processing to devise graph diffusion LMS strategies. Considering that most popular shift operators are not energy preserving and may result in a slow convergence speed, we introduce a preconditioned optimization strategy to improve the transient performance. As this may lead to an increased computational complexity, we further propose a computationally efficient algorithm. Explicit theoretical performance analyses in the mean and mean-square-error sense are provided. We also give alternative theoretical results that are tractable for large networks. Simulation results show the efficiency of the proposed algorithms and validate the theoretical models. Finally, we extend these node-invariant filter models to more flexible ones where each node in the graph seeks to estimate a local node-varying graph filter. This allows us to exploit more degrees of freedom in the filter coefficients to better model graph signals. We introduce an unsupervised clustering strategy to determine which nodes in the graph share the same graph filter and may collaborate to estimate its parameters. Numerical results on a real-word dataset illustrate the efficiency of the proposed methods.

The rest of the paper is organized as follows. Section II formulates the problem and provides the centralized solution. Section III introduces the distributed algorithms, namely, the graph diffusion LMS strategy and its preconditioned counterparts. Section IV provides their theoretical analyses in the mean and the mean-square sense. A clustering method is devised to estimate local node-varying graph filters in Section V. Numerical results in Section VI show the effectiveness of these algorithms and validate the theoretical models.

Notations: We use normal font letters to denote scalars, boldface lowercase letters to denote column vectors, and boldface uppercase letters to denote matrices. The m-th entry of a vector x is denoted by x m or [x] m when needed, the (m, n)th entry of a matrix X is denoted by x mn or [X] m,n when needed, the k-th row of a matrix X is denoted by [X] k,• . We use the symbol ⊗ to denote the Kronecker product and the symbol Tr(•) to denote the trace operator. The operator col{•} stacks the column vector entries on top of each other. The symbol vec(X) refers to the vectorization operator that stacks the columns of a matrix on top of each other. Operator x = diag(X) stores the diagonal entries of X into vector x, and X = diag(x) is a diagonal matrix containing the vector x along its main diagonal. The symbol bdiag{•} forms a matrix from block arguments by placing each block immediately below and to the right of its predecessor. The symbol • b,∞ denotes the block maximum norm of a matrix. The symbols 1 and I are the vector of all ones and the identify matrix of appropriate size, respectively. The symbols ρ(•) and λ max (•) denote the spectral radius and the maximum eigenvalue of its matrix argument, respectively.

II. PROBLEM FORMULATION AND CENTRALIZED SOLUTION

We consider an undirected, weighted and connected graph G = (N , E, W ) of N nodes, where N = {1, 2, . . . , N } is the set of nodes, and E is the set of edges such that if node k is connected to node , then (k, ) ∈ E. We denote by N k the neighborhood of node k including itself, namely,

N k = { : = k ∨ ( , k) ∈ E}. Matrix W ∈ R N ×N
is the adjacency matrix whose (k, )-th entry w k assigns a weight to the relation between vertices k and . Since the graph is undirected, W is a symmetric matrix. The degree matrix D diag(W 1) is a diagonal matrix whose i-th diagonal entry is the degree of node i which is equal to the sum of all the weights of edges incident at node i. The combinatorial Laplacian matrix is defined as L D -W which is a real, symmetric, positive semi-definite matrix. We further assume that the graph is endowed with a graph shift operator defined as an N × N shift matrix S whose entry s k can be non-zero only if k = or (k, ) ∈ E. Although the shift matrix can be any matrix that captures the graph topology for the problem at hand [START_REF] Ortega | Graph signal processing: Overview, challenges, and applications[END_REF], popular choices are the graph 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], the adjacency matrix [START_REF] Sandryhaila | Discrete signal processing on graphs[END_REF], and their normalized counterparts.

A graph signal is defined as x = [x 1 , . . . , x N ] ∈ R N where x k is the signal sample associated with node k. Let x(i) denote the graph signal at time i. Operation Sx is called graph shift. It can be performed locally at each node k by linearly combining the samples x from its neighboring nodes ∈ N k .

A. Graph filter and data model

In this paper, we focus on linear shift-invariant FIR graph filters H : R N ×N → R N ×N of order M , which are polynomials of the graph-shift operator [START_REF] Sandryhaila | Discrete signal processing on graphs[END_REF]:

H M -1 m=0 h o m S m , (1)
where

h o = {h o m } M -1
m=0 are the scalar filter coefficients. With the definitions of graph signal and graph shift operator, one common filtering model assumes that the filtered graph signal y(i) is generated from the input graph signal x(i) as follows [START_REF] Sandryhaila | Discrete signal processing on graphs[END_REF], [START_REF] Nassif | A graph diffusion LMS strategy for adaptive graph signal processing[END_REF]:

y(i) = Hx(i) + v(i) = M -1 m=0 h o m S m x(i) + v(i), (2)
where v(i) = [v 1 (i), . . . , v N (i)] ∈ R N denotes an i.i.d. zeromean noise independent of any other signal and with covari-

ance matrix R v = diag{σ 2 v,k } N k=1 .
For each node k, the filtered signal y k (i) can be computed by linearly combining the input signals at nodes located in an (M -1)-hop neighborhood [START_REF] Shuman | The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains[END_REF]. This model however assumes the instantaneous diffusion of information over the graph since S m x(i) translates x(i) without time delay. As this assumption may appear as a serious limitation, we consider the more general model embedding the temporal dimension as follows [START_REF] Isufi | 2-dimensional finite response graph-temporal filters[END_REF], [START_REF]Distributed diffusion adaptation over graph signals[END_REF]:

y(i) = M -1 m=0 h o m S m x(i -m) + v(i). (3)
Observe that the input signal x(i) in (2) has been replaced by x(i -m) in (3), i.e., the m-hop spatial shift S m is now carried out in m time slots. This model implements an FIR filter in both graph domain and temporal domain. By retaining the following shifted signals that form the N × M -1 matrix:

X r = x(i -1), Sx(i -2), . . . , S M -2 x(i -M + 1) , (4)
note that only one shift is required at time instant i to produce the filtered signal y(i). This means that y k (i) can be computed using only local information available within the one-hop neighborhood of node k. Let Z(i) denote the N × M matrix given by:

Z(i) x(i), Sx(i -1), . . . , S M -1 x(i -M + 1) , (5)
then model (3) can be written alternatively as:

y(i) = Z(i)h o + v(i) (6)
From model ( 6), sample y k (i) at node k can be written as:

y k (i) = z k (i)h o + v k (i), (7)
where z k (i) is the k-th row of Z(i) given by:

z k (i) col [x(i)] k , [Sx(i-1)] k , . . . , [S M -1 x(i-M +1)] k . (8
) Observe in [START_REF] Chen | Discrete signal processing on graphs: Sampling theory[END_REF] that each node shares the same filter coefficient vector h o . The corresponding graph filter (1) is referred to as node-invariant graph filter. A more flexible model was introduced in [START_REF] Segarra | Optimal graph-filter design and applications to distributed linear network operators[END_REF], and called a node-variant graph filter. It allows the filter coefficients to vary across nodes as follows:

H M -1 m=0 diag(h (m) )S m , (9)
with

h (m) ∈ R N . If h (m) = h m 1
for all m, model ( 9) reduces to the node-invariant model [START_REF] Hua | A preconditioned graph diffusion LMS for adaptive graph signal processing[END_REF]. Otherwise, each node applies different weights to the shifted signal S m x. Then y k (i) in ( 7) can be re-written as:

y k (i) = z k (i)h o k + v k (i), (10)
where

h o k ∈ R M is the filter coefficient vector at node k collected into h (m) , i.e., [h o k ] m = [h (m) ] k .
In this work, we seek to estimate h o k from the filtered graph signal y k (i) and inputs z k (i), in a distributed, collaborative and adaptive manner. Distributed algorithms such as the diffusion LMS exist in the literature to address single-task and multitask inference problems with similar data models as [START_REF] Chen | Discrete signal processing on graphs: Sampling theory[END_REF] or [START_REF] Sandryhaila | Discrete signal processing on graphs[END_REF]. In this work, however, regressors z k (i) in ( 8) are raised from graph shifted signals. This paper aims to exploit the graph shift structure in the regression data and incorporate it into the formulation of the distributed algorithm -see expression [START_REF] Sayed | Diffusion strategies for adaptation and learning over networks: an examination of distributed strategies and network behavior[END_REF] further ahead. In the sequel, first, we shall study the case where the filter coefficients are common for all nodes, i.e. h o k = h o , ∀k ∈ N . We shall show how to estimate h o from streaming data {y(i), x(i)} in a centralized way and then, in a distributed way. Next, we shall assume that there are clusters of nodes within the graph, and each node in the same cluster uses the same filter. This model is called a hybrid node-varying graph filter [START_REF] Gama | Convolutional neural networks via node-varying graph filters[END_REF]. We shall introduce an unsupervised clustering method to allow each node to identify which neighboring nodes it should collaborate with.

B. Centralized solution

Before introducing the distributed method, we first introduce the centralized solution. Consider the data model ( 6) and assume that {y(i), x(i), v(i)} are zero-mean jointly wide-sense stationary random processes. Estimating h o from {y(i), Z(i)} can be performed by solving the following problem:

h o = arg min h J(h), (11)
where J(h) denotes the mean-square-error criterion:

J(h) = E y(i) -Z(i)h 2 = E{y (i)y(i)} -2h r Zy + h R Z h, (12)
and the M × M matrix R Z and the M × 1 vector r Zy are given by:

R Z E{Z (i)Z(i)}, r Zy E{Z (i)y(i)}. ( 13
)
By setting the gradient vector of J(h) to zero, the optimal parameter vector h o can be found by solving:

R Z h o = r Zy . (14)
It can be verified that the (m, n)-th entry of R Z is given by:

[R Z ] m,n = Tr (S m-1 ) S n-1 R x (m -n) (15)
where R x (m) E{x(i)x (i -m)}. The m-th entry of the vector r Zy is given by:

[r Zy ] m = Tr (S m-1 ) R xy (m) , (16)
with R xy (m) E{y(i)x (i-m)} denoting cross correlation function, which is assumed independent of time i. Instead of solving ( 14), h o can be sought iteratively by using the gradient-descent method:

h(i + 1) = h(i) + µ r Zy -R Z h(i) , (17)
with µ > 0 a small step-size. Since the statistical moments are usually unavailable beforehand, one way is to replace them by the instantaneous approximations R Z ≈ Z (i)Z(i) and r Zy ≈ Z (i)y(i). This yields the LMS graph filter:

h(i + 1) = h(i) + µZ (i) y(i) -Z(i)h(i) . (18)
This stochastic-gradient algorithm is referred to as the centralized graph-LMS algorithm. In this centralized setting, each node k at each time instant i sends its data {x k (i), y k (i)} to a fusion center which will update h(i) according to [START_REF] Isufi | Autoregressive moving average graph filtering[END_REF]. Note that the step-size µ in ( 18) must satisfy 0 < µ <

2 λmax(R Z )
in order to guarantee stability in the mean under certain independence conditions on the data [START_REF] Sayed | Adaptive Filters[END_REF].

III. DIFFUSION LMS STRATEGIES OVER GRAPH SIGNALS

In this section, we seek to estimate the graph filter coefficients in a distributed fashion. First, we review the graph diffusion LMS strategy [START_REF]Distributed diffusion adaptation over graph signals[END_REF]. Then, a preconditioned algorithm is proposed to improve the transient performance. We also devise a computationally efficient counterpart of this algorithm.

A. Graph diffusion LMS

Consider the local data model ( 7) at node k. It is worth noting that, by retaining the past shifted signals {[S m-1 x(im)] : m = 1, . . . , M -1} at each node in the network from previous iterations, z k (i) can be computed locally at node k from its one-hop neighbors at each iteration i. Let R z,k E{z k (i)z k (i)} denote the M × M covariance matrix with (m, n)-th entry given by [START_REF]Distributed diffusion adaptation over graph signals[END_REF]:

[R z,k ] m,n = Tr [S m-1 )] k,• [S n-1 ] k,• R x (m -n) . (19) Considering the local cost J k (h) at node k: J k (h) = E|y k (i) -z k (i)h| 2 , (20)
the global cost ( 12) is now the aggregation of the local costs over the graph:

J(h) = N k=1 J k (h). (21)
In order to minimize [START_REF] Defferrard | Convolutional neural networks on graphs with fast localized spectral in[END_REF] in a decentralized fashion, there are several useful techniques, e.g., incremental strategy [START_REF] Bertsekas | A new class of incremental gradient methods for least squares problems[END_REF], consensus strategy [START_REF] Xiao | A scheme for robust distributed sensor fusion based on average consensus[END_REF] and diffusion strategy [START_REF]Diffusion adaptation over networks[END_REF]. Diffusion strategies are attractive since they are scalable, robust, and enable continuous learning and adaptation. In particular, the adapt-then-combine (ATC) diffusion LMS takes the following form at node k [START_REF]Distributed diffusion adaptation over graph signals[END_REF]:

z k (i) = x k (i), ∈N k s k [z (i -1)] 1 , . . . , ∈N k s k [z (i -1)] M -1 , (22a)
ψ k (i + 1) = h k (i) + µ k z k (i) y k (i) -z k (i)h k (i) , (22b)
h k (i + 1) = ∈N k a k ψ (i + 1), (22c)
where µ k > 0 is a local step-size parameter and {a k } are non-negative combination coefficients chosen to satisfy:

a k > 0, N =1 a k = 1, and a k = 0 if / ∈ N k . ( 23
)
This implies that the matrix A with ( , k)-th entry a k is a left-stochastic matrix, which means that the sum of each of its columns is equal to 1. In the first step (22a), each node k uses the first M -1 entries of z (i -1) from its one-hop neighbors and its own input sample x k (i) to compute z k (i). Note that the first M -1 entries of z k (i) then need to be retained for the next iteration. In the adaptation step (22b), each node k updates its local estimate h k (i) to an intermediate estimate ψ k (i + 1). In the combination step (22c), node k aggregates all the intermediate estimates ψ (i + 1) from its neighbors to obtain the updated estimate h k (i + 1).

B. Graph diffusion preconditioned LMS

The regressor z k (i) used in the adaptation step (22b) results from shifted graph signals while the shift matrix S is not energy preserving in general [START_REF] Gavili | On the shift operator, graph frequency, and optimal filtering in graph signal processing[END_REF]. This is due to the fact that the magnitude of the eigenvalues of the shift operator S are not uniformly equal to 1; the energy of the shifted signal S m x changes exponentially with m. Thus, the eigenvalue spread of R z,k may be large and the LMS update may suffer from slow convergence speed in this case [START_REF] Sayed | Adaptive Filters[END_REF]. To address this issue, albeit at an increased computational cost, we resort to a form of Newton's method. Focusing on the adaptation step, we have:

ψ k (i + 1) = h k (i) -µ k [∇ 2 h J k (h k (i))] -1 [∇ h J k (h k (i))], (24)
where ∇ 2 h J k (•) denotes the Hessian matrix for J k (•) and ∇ h J k (•) is its gradient vector, if available. For the quadratic cost function [START_REF] Coutino Minguez | Advances in distributed graph filtering[END_REF], expression [START_REF]Diffusion adaptation over networks[END_REF] would lead to:

ψ k (i + 1) = h k (i) + µ k R -1 z,k r zy,k -R z,k h k (i) , (25
) where r zy,k = E{z k (i)y k (i)}. Note that the second term on the RHS of (25) requires second-order moments. Since they are rarely available beforehand, we can replace r zy,k -R z,k h k (i) by the instantaneous approximation:

r zy,k -R z,k h k (i) ≈ z k (i)e k (i) (26) with e k (i) = y k (i) -z k (i)h k (i).
The adaptation step (25) becomes:

ψ k (i + 1) = h k (i) + µ k R -1 z,k (i)z k (i)e k (i), (27)
where R z,k (i) is an estimate for R z,k (i) which can possibly be obtained recursively:

R z,k (i) = (1-µ) R z,k (i-1)+µ z k (i)z k (i) , i ≥ 1, ( 28
)
where µ is a small factor that can be chosen in (0, 0.1] in practice. It can be verified that E{ R z,k (i)} = R z,k is an unbiased estimate when i → ∞. As discussed before, S may not be energy preserving and results in a large eigenvalue spread of R z,k , which may even be close to singular. The inverse R -1 z,k would then be ill-conditioned and lead to undesirable effects. To address this problem, it is common to use regularization [START_REF] Sayed | Adaptive Filters[END_REF]. We obtain the diffusion LMS-Newton algorithm:

ψ k (i + 1) = h k (i) + µ k I + R z,k (i) -1 z k (i)e k (i), (29a)
h k (i + 1) = ∈N k a k ψ (i + 1), (29b)
with ≥ 0 a small regularization parameter. Compared with the diffusion LMS algorithm [START_REF] Sayed | Diffusion strategies for adaptation and learning over networks: an examination of distributed strategies and network behavior[END_REF], algorithm (29) requires first to recursively estimate the Hessian matrix according to [START_REF] Chen | Diffusion LMS over multitask networks[END_REF] and then calculate

I + R z,k (i) -1
. This algorithm can lead to improved performance as shown in the sequel, but at the expense of additional computation cost.

In order to reduce the computational complexity of the LMS-Newton algorithm, we propose to use a preconditioning matrix P k does not depend on the graph signal x(i) in the adaptation step, instead of the Hessian matrix R z,k or its estimate R z,k . Since the large eigenvalue spread of the input covariance matrix R z,k results mainly from the shift matrix S and the filter order M , we construct an M ×M preconditioning matrix P k as follows:

P k diag{ [S (m-1) ] k,• 2 } M m=1 . (30)
The rationale behind [START_REF] Nassif | Multitask diffusion adaptation over asynchronous networks[END_REF] is that, in the case where

x(i) is i.i.d. with variance σ 2 , it follows from (19) that R z,k = σ 2 P k .
According to [START_REF] Nassif | Multitask diffusion adaptation over asynchronous networks[END_REF], matrix P k does not depend on x(i) and can be evaluated beforehand at each node k during an initial step. Each node k only requires to know the edge weights in its M -hop neighborhood, which can be performed in a decentralized manner. Interestingly, P k is a diagonal matrix, which means that the matrix product in the adaptation step does not require expensive matrix inversion. Following the same line of reasoning as for the Newton algorithm (29), a regularization term I M can be added to P k . This leads to:

D k = ( I M + P k ) -1 . ( 31
)
We arrive at the following preconditioned graph diffusion LMS strategy:

ψ k (i + 1) = h k (i) + µ k D k z k (i)e k (i), (32a)
h k (i + 1) = ∈N k a k ψ (i + 1). (32b)
At each iteration i, node k uses the local information to update the intermediate estimate ψ k (i + 1) in the adaptation step (32a). Then, in the combination step (32b), the intermediate estimates ψ (i + 1) from the neighborhood of node k are combined to get h k (i + 1). Although the preconditioning matrix D k is not the true Hessian matrix, we prove in Section IV that the algorithm converges to the optimal solution h o provided that it is stable.

C. Comparison with the graph diffusion LMS

We explain how preconditioning with (30) improves performance. For comparison purposes, let us first focus on the adaptation step of diffusion LMS. At each node k, the m-th entry of h k (i) is updated as follows:

[ψ k (i + 1)] m = [h k (i)] m + µ k z k (i)e k (i) m . ( 33
)
During the transient phase, the m-th entry [h k (i)] m exponentially converges to its optimal value with a time constant [START_REF] Sayed | Adaptive Filters[END_REF]:

τm ≈ 1 2µ k λ m ( 34
)
where λ m denotes the m-th eigenvalue of R z,k . Given µ k , the convergence rate of each entry of h k (i) then depends on the corresponding eigenvalue. Disparity between entries increases as the eigenvalue spread defined as λ max /λ min increases. The preconditioning matrix D k is diagonal at each node k, which means that the m-th entry of h k (i) in (32a) converges to its optimal value with a time constant:

τm ≈ 1 2µ k d k,m λ m ( 35
)
where d k,m denotes the m-th diagonal entry of D k . The convergence speed now depends on d k,m λ m . Considering the case where P k is proportional to R z,k , then d k,m is inversely proportional to λ m , which mitigates the effects of the eigenvalues spread. We shall analyze and illustrate in Section IV and Section VI, respectively, how this preconditioning improves convergence speed in more general cases.

IV. PERFORMANCE ANALYSIS

We shall now analyze the stochastic behavior of the diffusion preconditioned LMS (PLMS) algorithm [START_REF] Lorenzo | Adaptive least mean squares estimation of graph signals[END_REF] in the sense of mean and mean-square error. We introduce the following weight error vectors at each node k:

hk (i) = h o -h k (i), ψk (i) = h o -ψ k (i), (36)
and we collect them across the nodes into the network weight error vectors :

h(i) col{ h1 (i), h2 (i), . . . , hN (i)}, (37)
ψ(i) col{ ψ1 (i), ψ2 (i), . . . , ψN (i)}. ( 38
)
We refer to mean stability of the error vector h(i) if the limit superior of E h(i) is bounded. Furthermore, we will claim that the algorithm converges in the mean to the optimum if E h(i) converges to zero as i tends to +∞ regardless of the starting point. Mean-square stability refers to the case where the superior limit of E h(i) 2 is bounded.

Let us introduce the following N × N block matrices with individual entries of size M × M :

A A ⊗ I M , (39)
M bdiag{µ k I M } N k=1 , (40)
D bdiag{D k } N k=1 . (41)
The estimation error in (32a) can be written as:

e k (i) = y k (i) -z k (i)h k (i) = z k (i) hk (i) + v k (i). ( 42
)
Subtracting h o from both sides of (32a) and using the above relation, then stacking ψk (i) across the nodes, we obtain

ψ(i + 1) = (I N M -MDR z (i)) h(i) -MDp zv (i), (43)
where R z (i) is an N × N block matrix with entries of size M × M define as:

R z (i) bdiag{z k (i)z k (i)} N k=1 , (44)
and p zv (i) is an N × 1 block column vector with entries of size M × 1 given by:

p zv (i) col{z k (i)v k (i)} N k=1 . (45)
Subtracting h o from both sides of (32b), we obtain the block weight error vector:

h(i + 1) = A ψ(i + 1). (46)
Finally, combing [START_REF] Romero | Kernel-based reconstruction of space-time functions on dynamic graphs[END_REF] and [START_REF] Nassif | A graph diffusion LMS strategy for adaptive graph signal processing[END_REF], the network weight error vector h(i) of algorithm [START_REF] Lorenzo | Adaptive least mean squares estimation of graph signals[END_REF] evolves according to:

h(i + 1) = B(i) h(i) -A MDp zv (i), (47)
with

B(i) = A I N M -MDR z (i) . (48)
To proceed with the analysis, we introduce the following assumption.

Assumption 1 (independent inputs). The inputs z k (i) arise from a zero-mean random process that is temporally white with R z,k > 0.

A consequence of Assumption 1 is that z k (i) is independent of h (j) for all and j ≤ i. This independence assumption is not true in the current work. Two successive regressors z k involve common entries that cannot be statistically independent as in a conventional FIR implementation. However, when the step-size is sufficiently small, conclusions derived under this assumption tend to be realistic. For more details and discussions, see [START_REF] Sayed | Adaptive Filters[END_REF]Section 16.4]. Since this assumption helps to simplify the derivations without constraining the conclusions, it is widely used in the literature of adaptive filters and adaptive networks [START_REF]Diffusion adaptation over networks[END_REF], [START_REF] Sayed | Adaptive Filters[END_REF]. We shall see in Section VI that the resulting expressions match well the simulation results for sufficiently small step-sizes.

A. Mean-error behavior analysis

Taking expectations of both sides of (47), using the fact that Ep zv (i) = 0, and applying Assumption 1, we find that the network mean error vector evolves according to:

E h(i + 1) = B E h(i), (49)
where:

B EB(i) = A (I N M -MDR z ), (50)
R z ER z (i) = bdiag{R z,k } N k=1 . ( 51
)
Theorem 1 (Convergence in the mean). Assume that data model [START_REF] Chen | Discrete signal processing on graphs: Sampling theory[END_REF] and Assumption 1 hold. Then, for any initial condition, algorithm (32) converges asymptotically in the mean toward the optimal vector h o if, and only if, the step-sizes in M are chosen to satisfy:

ρ A (I N M -MDR z ) < 1, (52)
where ρ(•) denotes the spectral radius of its matrix argument.

In the case where the signal x(i) is i.i.d, a sufficient condition for (52) to hold is to choose µ k such that:

0 < µ k < 2 λ max (D k R z,k ) , k = 1, . . . , N. (53)
Proof. The weight error vector h(i) converges to zero if, and only if, the coefficient matrix B in ( 49) is a stable matrix, namely, ρ(B) < 1. Since any induced matrix norm is lower bounded by the spectral radius, we have the following relation in terms of block maximum norm [START_REF]Diffusion adaptation over networks[END_REF]:

ρ(B) ≤ A (I N M -MDR z ) b,∞ ≤ A b,∞ • I N M -MDR z b,∞ = I N M -MDR z b,∞ , (54)
where the last equality follows from the fact that A is left stochastic, which implies that A b,∞ = 1 from Lemma D.4 of [START_REF]Diffusion adaptation over networks[END_REF]. Matrix R z is block diagonal if x(i) is i.i.d. Since D is also diagonal, their product is symmetric, and MDR z is a block diagonal symmetric matrix. Then, following Lemma D.5 of [START_REF]Diffusion adaptation over networks[END_REF], its block maximum norm agrees with its spectral radius:

I N M -MDR z b,∞ = ρ(I N M -MDR z ). (55)
Combining ( 54) and ( 55), we verify that condition (53) ensures the stability of B.

B. Mean-square-error behavior analysis

We shall now study the mean-square-error behavior of algorithm [START_REF] Lorenzo | Adaptive least mean squares estimation of graph signals[END_REF]. Let Σ be any N M × N M positive semidefinite matrix that we are free to choose. The freedom in selecting Σ will allow us to derive different performance measures about the network and the nodes. We consider the weighted mean-square error vector, i.e., E h(i) 2 Σ , where

h(i) 2 Σ h (i)Σ h(i).
From Assumption 1 and Ep zv (i) = 0, using [START_REF]Distributed diffusion adaptation over graph signals[END_REF], we obtain the following variance relation:

E h(i + 1) 2 Σ = E h(i) 2 Σ + E A MDp zv (i) 2 Σ , ( 56
)
where Σ E{B (i)ΣB(i)}. Let σ denote the (N M ) 2 × 1 vector obtained by vectorizing matrix Σ, namely, σ = vec(Σ). With some abuse of notation, we shall use • 2 σ to also refer to the quantity • 2 Σ when it is more convenient. Let σ = vec(Σ ). Considering that vec(U ΣW ) = (W ⊗U )σ, it can be verified that:

σ = F σ, ( 57
)
where F is the (N M ) 2 × (N M ) 2 matrix given by

F E{B (i) ⊗ B (i)} = I (N M ) 2 -I N M ⊗ R z DM -R z DM ⊗ I N M + O(M 2 ) (A ⊗ A), (58)
where O(M 2 ) denotes E{R z (i)DM⊗R z (i)DM}, which depends on the square of the step-sizes, {µ 2 k }. While we can continue the analysis by taking this factor into account as was done in other studies [START_REF] Sayed | Fundamentals of adaptive filtering[END_REF], it is sufficient for the exposition to focus on the case of sufficiently small step-sizes where terms involving higher powers of the step-sizes {µ k } can be ignored. Following the same line of reasoning, for sufficiently small step-sizes {µ k }, F can be approximated by:

F ≈ B ⊗ B . (59)
The second term on the RHS of ( 56) can be written as:

E A MDp zv (i) 2 Σ = Tr(ΣG), (60)
where

G A MDSDMA, (61)
S E{p zv (i)p zv (i)} = bdiag{σ 2 v,k R z,k } N k=1 . (62)
Using the property Tr(ΣW ) = [vec(W )] σ, combining ( 57) and ( 60), the variance relation ( 56) can be re-written as:

E h(i + 1) 2 σ = E h(i) 2 F σ + [vec(G )] σ. ( 63
)
Theorem 2 (Stability in the mean-square). Assume that data model [START_REF] Chen | Discrete signal processing on graphs: Sampling theory[END_REF] and Assumption 1 hold. Algorithm (32) converges in the mean-square sense if the matrix F in (58) is stable. Assume further that the step-sizes are sufficiently small such that (59) is a reasonable approximation. In that case, the stability of F is ensured if B is stable.

Proof. Iterating (63) starting from i = 0, we obtain

E h(i+1) 2 σ = E h(0) 2 F i+1 σ +[vec(G )] i j=0 F j σ, (64)
with initial condition h(0) = h o -h(0). Provided F is stable, F i → 0 as i → ∞, then the first item on the RHS of (64) converges to zero and the second item converges to a finite value. The weighted mean-square error converges to a finite value as i → ∞ which implies that the algorithm (32) will converge in the mean-square sense if F is stable. Under the sufficiently small step-sizes assumption where the higher-order terms of F in (58) can be neglected, approximation ( 59) is reasonable. The eigenvalues of F are all the products of the eigenvalues of B, which means that ρ(F ) = [ρ(B)] 2 . It follows that F is stable if B is stable. Therefore, when the graph signal x(i) is i.i.d, according to Theorem 1, condition (53) ensures mean-square stability of the algorithm under the assumed approximation [START_REF]1981-2010 U.S. climate normals[END_REF].

Theorem 3 (Network transient MSD). Assume sufficiently small step-sizes that ensure mean and mean-square stability.

The network transient mean-square deviation (MSD) defined as ζ(i) = 1 N E h(i) 2 evolves according to the following recursion for i ≥ 0:

ζ(i + 1) = ζ(i) + 1 N [vec(E{ h(0) h (0)})] (F -I (N M ) 2 ) + [vec(G )] F i vec(I N M ). (65)
Proof. Comparing (64) at time i + 1 and i, E h(i + 1) 2 σ is related to E h(i) 2 σ as follows:

E h(i + 1) 2 σ = E h(i) 2 σ + [vec(G )] F i σ+ [vec(E{ h(0) h (0)})] (F -I (N M ) 2 )F i σ. ( 66
)
Substituting σ by 1 N vec(I N M ) leads to (65). Although expression (65) gives a compact form of the transient MSD model, it may not be practical to use since F is of size (N M ) 2 × (N M ) 2 and may become huge for large networks or high order filters. For example, in the simulation Section VI, we considered a network consisting of N = 60 nodes and a graph filter of degree M = 5. Matrix F is of size 90000 × 90000 and requires a prohibitive amount of computational time and memory space (about 60GB in that case). To tackle this issue, we make use of the following properties of the Kronecker product:

vec(XY Z) = (Z ⊗ X)vec(Y ) (67)
Tr(XY ) = vec(Y ) vec(X). (68)
This leads to:

Corollary 1 (Alternative network transient MSD expression).

ζ(i + 1) = ζ(i) + 1 N Tr B i G(B i ) + h(0) h (0) (B i+1 ) B i+1 -(B i ) B i . ( 69
)
While the update with F has a computation complexity of order O((N M ) 2 ), using B only requires matrix manipulations of order O(N M ).

Corollary 2 (Network steady-state MSD). Consider sufficiently small step-sizes to ensure mean and mean-square convergence. The network steady-state MSD is given by

ζ = 1 N [vec(G )] (I (N M ) 2 -F ) -1 vec(I N M ). (70)
Proof. The network steady-state MSD is defined as:

ζ = lim i→∞ 1 N E{ h(i) 2 }. ( 71
)
If F is stable, we obtain from (63) as i → ∞:

lim i→∞ E h(i) 2 (I (N M ) 2 -F )σ = [vec(G )] σ. ( 72
)
We obtain (70) by substituting σ in (72) by 1

N (I (N M ) 2 - F ) -1 vec(I N M ).
Following the same line of reasoning as for the transient MSD model, the steady-state MSD given by (70) can be equivalently expressed as:

Corollary 3 (Alternative steady-state MSD expression). ζ = 1 N ∞ i=0 Tr B i G(B i ) ( 73
)
This expression is obtained by a series expansion of (70). In practice, a limited number of iterations can be used instead of the upper limit index i → ∞ to obtain an accurate result.

V. UNSUPERVISED CLUSTERING FOR HYBRID NODE-VARYING GRAPH FILTER

In Section III, we investigated the scenario where the nodes in a graph share a common filter coefficient vector. Now we extend this model to the more flexible case [START_REF] Sandryhaila | Discrete signal processing on graphs[END_REF], which allows the filter coefficients to vary across the graph. We further assume that the graph is decomposed into Q clusters of nodes C q and, within each cluster C q , there is a common filter coefficient vector h o q to estimate, namely,

h o k = h o q , if k ∈ C q . ( 74
)
We assume that there is no prior information on the clusters composition and that the nodes do not know which other nodes share the same estimation task. Applying the algorithm [START_REF] Lorenzo | Adaptive least mean squares estimation of graph signals[END_REF] within this context may result a bias due to aggregate intermediate estimates from different data models. To address this issue, automatic network clustering strategies may be used [START_REF] Chen | Diffusion LMS over multitask networks[END_REF], [START_REF] Zhao | Distributed clustering and learning over networks[END_REF]- [START_REF] Khawatmi | Decentralized clustering and linking by networked agents[END_REF] in order to inhibit cooperation between nodes from different clusters. These methods are based on local stand-alone estimation strategies that may not be efficient for the current context. Basically, the polynomial form (1) of graph filters does not make the estimation of the filter coefficients reliable enough for the higher degrees. In the following, we tackle this problem by devising a clustering strategy based on the PLMS.

First, we introduce the N × N instantaneous clustering matrix E i , whose ( , k)-th entry shows if node k believes at time i that its neighboring node belongs to the same cluster or not, namely,

[E i ] k = 1, if ∈ N k and k believes that h o k = h o , 0, otherwise. ( 75
)
At each time instant i, node k infers which neighbors belong to its cluster based on the non-zeros entries of the k-th column of E i . We collect these entries into a set N k,i , so that node k only combines the intermediate estimates from its neighbors in N k,i . Condition [START_REF] Sayed | Adaptive networks[END_REF] on the combination coefficients becomes:

a k > 0, N =1 a k = 1, and a k = 0 if / ∈ N k,i , (76)
Since the clustering information is unknown beforehand, we propose to learn E i in an online way by computing a Boolean variable b k (i) defined as follows:

b k (i) = 1, if ψ (i + 1) -h k (i) 2 ≤ β, 0, otherwise, (77)
with β > 0 a preset threshold. Variable b k (i) is defined from the 2 -norm distance between the estimates at two neighboring nodes. If this distance is smaller than the threshold β, the two nodes are then assigned to the same cluster. Note that the distance between ψ (i + 1) and h k (i) is used in (77), instead of the distance between h (i + 1) and h k (i + 1), in order to merge the learning and the clustering processes. Consider the learning process at any node k defined in [START_REF] Lorenzo | Adaptive least mean squares estimation of graph signals[END_REF]. Information about clusters is used in the combination step (32b), where N k now denotes the neighboring nodes of k that share the same estimation task as node k. This information should be available as soon as possible in order to avoid estimation bias.

Considering the distance between h (i + 1) and h k (i + 1) to decide if nodes k and are in the same cluster would allow to update the composition of sets N k and N in the combination step (32b) used to calculate parameter vectors h (i + 2) and h k (i + 2). This latency time can be shortened by considering the distance between ψ (i + 1) and h k (i) right after the adaptation step (32a), and using this information to define N k in the combination step. This strategy usually fails if left as is, because the estimation of higher-order coefficients is not reliable enough and results in bad clustering performance. We propose to estimate this distance from M k principal components of the estimates, to be defined, of the estimates. Because R z,k cannot be reasonably used to perform a Principal Component Analysis (PCA) of the input data, as it is rarely available beforehand and would involve significant additional computational effort, we suggest instead using matrix P k . As for the PLMS, the rationale behind this is that R z,k = σ 2 P k when x(i) is i.i.d. with variance σ 2 . Another interest lies in that P k is a diagonal matrix, which greatly simplifies calculations. Without loss of generality, consider that the diagonal entries of P k are in decreasing order. Projecting data onto the first M k principal axes then reduces to selecting their first M k entries and set the other entries to zero. Dimension M k can be determined by setting the ratio of explained variance to total variance to some desired level τ as follows:

minimize M k subject to

M k m=1 πk,m ≥ τ (78)
with πk,m = [p k ] m /Tr(P k ) an approximation of the proportion of total variance [START_REF] Jolliffe | Principal component analysis: a review and recent developments[END_REF]. In practice, we can use a predefined threshold τ ∈ [0.9, 1) to decide how many entries should be retained. The Boolean variable (77) becomes:

b k (i) = 1, if ψ (i+1)-h k (i) 2 h k (i) 2 ≤ β, 0, otherwise, (79)
with ψ (i + 1) and h k (i) the first M k entries of ψ (i + 1) and h k (i) respectively. Compared to (77), note that the distance in (79) that accounts for the similarity between ψ (i + 1) and h k (i) has been normalized. We suggest to choose β ∈ (0, 0.01] to lower the false detection rate. To reduce noise effects, we further introduce a smoothing step:

t k (i) = νt k (i -1) + (1 -ν)b k (i), (80)
where t k (i) is a trust level, and ν is a forgetting factor in (0, 1) to balance the past and present cluster assignments. Once the trust level t k (i) exceeds a preset threshold θ, which can be chosen in [0.5, 1), node k concludes that node belongs to its cluster, that is,

[E i ] k = 1, if t k (i) ≥ θ, 0, otherwise. (81)
Based on [E i ] k , each node k determines at each time instant i those nodes that it believes they belong to the same cluster, updates the combination coefficients according to (76), and finally combines the estimates from its neighbors with (32b).

VI. NUMERICAL RESULTS

A. Experiment with i.i.d. input data

We first considered a zero-mean i.i.d. Gaussian graph signal x(i)

with covariance R x = diag{σ 2 x,k } N k=1 . Variances σ 2
x,k were randomly generated from the uniform distribution U(1, 1.5). In this setting, the graph signal sample x k (i) was independent of x (j) for all and j ≤ i. We assumed the linear data model [START_REF] Chen | Discrete signal processing on graphs: Sampling theory[END_REF]. The graph filter order was set to M = 5 and the coefficients h o m were randomly generated from the uniform distribution U(0, 1). Noise v(i) was zero-mean Gaussian with covariance

R v = diag{σ 2 v,k } N k=1 . Variances σ 2
v,k were randomly generated from the uniform distribution U(0.1, 0.15). We considered this data model with an Erdős- Rényi random graph and on a sensor network graph. Both consisted of N = 60 nodes. The Erdős-Rényi random graph was generated in a similar construction as in [START_REF] Mei | Signal processing on graphs: Causal modeling of unstructured data[END_REF]. Namely, it was obtained by generating an N × N symmetric matrix S whose entries were governed by the Gaussian distribution N (0, 1), and then thresholding edges to be between 1.2 and 1.8 in absolute value. Then, the edges were soft thresholded by 1.1 to be between 0.1 and 0.7 in magnitude. The shift matrix S was normalized by 1.1 times its largest eigenvalue. The sensor network was generated by using GSPBOX [START_REF] Perraudin | GSPBOX: A toolbox for signal processing on graphs[END_REF]. Each node was connected to its 5 nearest neighbors. The shift matrix was the normalized adjacency matrix, that is, S = W 1.1λmax(W ) . In this case, all the eigenvalues of S are smaller than 1 and the energy of the shifted signal S m x diminishes for large m. The smallest eigenvalue λ min (R z,k ) was very small, and, for some node, it was close to 0.

With this setting, we compared the diffusion LMS algorithm [START_REF] Sayed | Diffusion strategies for adaptation and learning over networks: an examination of distributed strategies and network behavior[END_REF], the diffusion LMS-Newton (LMSN) algorithm [START_REF] Chen | Diffusion LMS for multitask problems with overlapping hypothesis subspaces[END_REF], and the diffusion preconditioned LMS (PLMS) algorithm [START_REF] Lorenzo | Adaptive least mean squares estimation of graph signals[END_REF]. Simulated results were averaged over 500 Monte-Carlo runs. For the LMSN and PLMS algorithms, we set the regularization parameter as = 0.01. We ran algorithms ( 22), ( 29) and ( 32) by setting a ,k = 1

|N k | for ∈ N k . We used a uniform step-size for all nodes, i.e., µ k = µ for all k. We also considered the -normalized LMS ( -NLMS) method for comparison purposes. In this case, the adaptation step (32a) is substituted by:

ψ k (i + 1) = h k (i) + µ k z k (i) 2 + z k (i)e k (i). (82)
With the Erdős-Rényi graph, we compared the performance of the LMS, PLMS, LMSN and -NLMS algorithms. We set µ = {0.08, 0.008, 0.01, 0.05}, respectively. The network MSD performance of each algorithm is reported in Fig. 1.

The theoretical transient and steady-state MSD are also reported. With the sensor network graph, we compared the performance of the LMS, PLMS and LMSN algorithms. We set µ = {0.08, 0.005, 0.0055}, respectively. The performance of each algorithm is reported in Fig. 2(a). In Fig. 1, we observe that the diffusion -NLMS converged slower than all other algorithms. Observe that the diffusion LMSN and PLMS algorithms converged faster than the LMS algorithm for both graphs, and the diffusion PLMS performed similarly compared with the LMSN in terms of convergence rate. Also, note that the theoretical results match well the simulated curves.

In a second experiment, we considered the normalized graph Laplacian matrix S = D -1 2 LD -1 2 , and the adjacency matrix W , as shift operators. For the normalized graph Laplacian, λ max (R z,k ) was large for all nodes. Therefore, for the diffusion LMS algorithm, the step-size was chosen relatively small to guarantee convergence. We set µ = {0.004, 0.01, 0.008} for the LMS, LMSN and PLMS. The results are reported in Fig. 2(b). For the adjacency matrix, we used uniform step-sizes µ = {0.02, 0.018} for the LMSN and PLMS, respectively. The step-size was set to µ k = 0.05

• 2 λmax(R z,k )
for each node k for the LMS update in order to achieve the same steady-state MSD. The results are reported in Fig. 2(c). We observe in Fig. 2 that the diffusion LMSN and PLMS algorithms converged faster than the LMS algorithm with the three graph shift operators. The PLMS algorithm achieved the same performance as the LMSN algorithm with a lower computational complexity.

B. Experiment with correlated input data

We tested the algorithms over the sensor network graph with correlated graph signals. We first considered a zeromean i.i.d. Gaussian graph signal driven by a non-diagonal covariance matrix R x . This means that the input data were correlated over the vertex domain, but uncorrelated over time. Matrix R x was generated as R x = V diag{σ 2

x,k } N k=1 V , with σ 2

x,k randomly chosen from the uniform distribution U(1, 1.5) and V is the graph Fourier transform matrix. The graph shift operator was defined by the normalized adjacency matrix. The filter degree was set as M = 3. We observe in Fig. 3 that the proposed diffusion PLMS algorithm performed as well as the LMSN algorithm, and converged faster than the diffusion LMS algorithm. Next, we considered a graph signal x(i) that was correlated over both vertex and time domains. We assumed that x(i) is a Gaussian process with zero mean and covariance matrix R x satisfying the discrete Lyapunov equation:

SR x S -R x + I = 0. ( 83
)
Graph signal sample x(i) was related to x(i -1) as follows:

x(i) = Sx(i -1) + w(i) (84)
with S the normalized adjacency matrix and w(i) a zero-mean i.i.d. Gaussian noise with covariance I N . It can be checked that x(i) is wide-sense stationary with E{x(i)} = 0, and R x (τ ) = S τ R x (0) for all τ > 0, where R x (0) satisfies the Lyapunov equation ( 83). The graph filter order was set as M = 3. The step-sizes were set to µ = {0.1, 0.038, 0.03} for the LMS, LMSN and PLMS, respectively. The regularization parameter was set to = 0.1. Fig. 4 depicts the simulated and theoretical MSD performance. We observe that, due to correlation over time, the diffusion LMSN method converged faster than the PLMS. However, the proposed PLMS algorithm still performed better than the diffusion LMS method.

C. Clustering method for node-varying graph filter

Finally, we considered a scenario where nodes do not share the same filter coefficients. We assumed the linear data model [START_REF] Sandryhaila | Discrete signal processing on graphs[END_REF]. The graph shift operator was defined by the normalized adjacency matrix, and the graph filter order was set as M = 3. The nodes were grouped into three clusters: , with basic clustering mechanism (77) and M k = M for all k, the oracle PLMS algorithm where the clusters are assumed to be known a priori, the PLMS algorithm without clustering mechanism, and the non-cooperative algorithm where a k = 1 if k = and zero otherwise. All algorithms used the adaptation step (32a) with the same step-size µ k = 0.01 for all k. Parameters {τ, β, θ, ν} were set to {0.9, 0.01, 0.5, 0.98}, respectively. As shown in Fig. 5, the non-cooperative method did not achieve acceptable MSD level. The main reason is that, with the normalized adjacency matrix as graph shift operator S, the entries of z k (i) in ( 8) corresponding to higher powers of S are significantly diminished, resulting in poor estimation performance of filter coefficients when nodes cannot cooperate. The PLMS algorithm without clustering mechanism did not achieve good performance too because it has been designed to converge toward a consensual solution h o , which does not make sense for this scenario. The PMLS with clustering mechanism (71) and M k = M achieved slightly improved performance because the estimation of higher-order coefficients in h k was not reliable enough, leading to incorrect clustering. The proposed PLMS algorithm with clustering mechanism (78)-(81) performed as well as the oracle algorithm. Fig. 6 (a) shows the topology of the graph given by the adjacency matrix A (and the shift matrix S).

D. Reconstruction on U.S. temperature dataset

We considered a dataset that collects hourly temperature measurements at N = 109 stations for T = 8759 hours across the United States in 2010 [START_REF]1981-2010 U.S. climate normals[END_REF]. An undirected graph, illustrated in Fig. 10, was constructed according to the nodes coordinates by using the k-NN approach (k = 7) of GSPBOX. In the first experiment, the dataset was divided into a training set containing T train = 6570 hours data (about 75% of total). The remaining data were assigned to the test set. The goal of this experiment was to learn a graph filter that minimizes the reconstruction error over the training set, i.e.,

min Ttrain i=1 N k=1 |y k (i) - M m=1 h m,k [S m x(i -m + 1)] k | 2 , (85)
where y(i) is the ground truth temperature at time i, and x(i) is the partial observation given by x(i) = diag(1 S )y(i). Here 1 S denotes the set indicator vector, whose k-th entry is equal to one if node k is sampled, and zero otherwise. The sampling set, illustrated in Fig. 10 (a) was fixed over time in the first experiment. The normalized adjacency matrix was set as graph shift operator. Graph filter degree was set to M = 4. Note that if h m,k = h m for all k, problem (85) refers to the single-task problem where all the nodes seek to find common graph filter coefficients; see model [START_REF] Hua | A preconditioned graph diffusion LMS for adaptive graph signal processing[END_REF]. Otherwise, problem (85) refers to the multitask problem; see model [START_REF] Tsitsvero | Signals on graphs: Uncertainty principle and sampling[END_REF]. We ran different models and algorithms on the training set to learn graph filter coefficients. In Fig. 11, we provide the true temperature and the reconstructed ones obtained by the different algorithms at an unobserved node, black circled in Fig. 10, over the last 120 hours samples of the test set. For comparison purposes, reconstruction results of the Kernel Kalman Filter (KKF) and the Kernel Ridge Regression (KRR) in [START_REF] Romero | Kernel-based reconstruction of space-time functions on dynamic graphs[END_REF] are also reported in Fig. 11. We observe that the single-task diffusion LMSN and multitask diffusion LMS were not able to reconstruct the true temperature, whereas the multitask diffusion PLMS and the multitask diffusion LMSN showed a good reconstruction performance, and performed better than the KKF and KRR at the selected unobserved node. To evaluate the performance over all unobserved nodes on the test set, we considered the normalized mean square error (NMSE) defined as: where ŷ(i) denotes the reconstructed estimate at time i, 1 S is the set indicator vector whose k-th entry is equal to one if node k has not been sampled, and zero otherwise. The results are reported in Table I. We observe that the multitask diffusion PLMS performed as well as the LMSN at a lower computational cost, and both performed better than the KFF and KRR. Finally, Figure 12 reports the original topology and the clusters learned by the multitask diffusion PLMS. In the second experiment, we divided the dataset into two parts. The first part contained the first 4200 hours sampled at the nodes showed in Fig. 10 (a), and the second part contained the remaining hours sampled at the nodes showed in Fig. 10 (b). This means that the sampling set abruptly changed at time t = 4201 (the black circled node was unobserved in both case). We applied the multitask diffusion PLMS method over the entire dataset. In Fig. 13 (a), the reconstructed temperature at the unobserved black circled node is reported from time t = 4100 to t = 4300 by using the filter coefficients learned up to time t = 4000. As expected, we can notice that the reconstruction performance was successful from t = 4100 to t = 4200, and dramatically deteriorated after t = 4201. This is due to the fact that the sampling set changed, which led to a drift in the filter coefficients to estimate. Figure 13 (b) depicts the reconstruction behavior over last 120 hours using the filter coefficients learned up to time t = 8600. It can be observed that the proposed method was able to track the drift in the filter coefficients.

VII. CONCLUSION

Diffusion LMS strategies were considered to estimate graph filter coefficients in an adaptive and distributed manner. A diffusion LMS with Newton-like descent procedure was first proposed to achieve improved convergence rate, since usual algorithms may suffer from ill conditioning effects due to the use of non-energy preserving graph shift operators. A preconditioned diffusion LMS strategy, which does not require computationally intensive matrix inversion and only uses local information, was then devised to reduce the computational burden. Its convergence behavior was analyzed in the mean and mean-square-error sense. Finally, for hybrid node-varying graph filters, a clustering mechanism to be used with the preconditioned diffusion LMS was proposed. Simulation results validated the theoretical models and showed the efficiency of the proposed algorithms. In this work, we assumed that the graph shift operator is known and time invariant. In future works, we will consider ways to estimate the graph shift operator and the filter coefficients simultaneously. Because of their improved flexibility, which allows to address a variety of non-linear identification problems, we will also focus on non-parametric methods such as [START_REF] Koppel | Decentralized online learning with kernels[END_REF], [START_REF] Pradhan | Exact nonparametric decentralized online optimization[END_REF] and see how they can be incorporated into GSP framework.

Fig. 1 :

 1 Fig. 1: Network MSD performance with the Erdős-Rényi graph.

Fig. 2 :

 2 Fig. 2: Network MSD performance for different types of shift operators with the sensor network.

Fig. 3 :

 3 Fig. 3: Network MSD performance with a vertex domain correlated input signal.

Fig. 4 :

 4 Fig. 4: Network MSD performance with input graph signal correlated over both vertex and time domains.

C 1 =

 1 {1, . . . , 20}, C 2 = {21, . . . , 40}, and C 3 = {41, . . . , 60}. The optimal graph filter coefficients h o k were set according to [0.5 0.4 0.9] if k ∈ C 1 , [0.3 0.1 0.4] if k ∈ C 2 , and [0.9 0.3 0.7] if k ∈ C 3 . We considered for comparison purpose the PLMS algorithm with clustering mechanism (78)-(81)

Fig. 6 (Fig. 5 :

 65 Fig. 5: Network MSD performance for different clustering algorithms.

Fig. 6 :Fig. 7 :

 67 Fig. 6: Graph topology and clusters. Next, we considered that optimal parameter vectors h o k change over time while clusters remain unchanged. Nodes were grouped into two clusters C 1 = {1, . . . , 30} and C 2 = {31, . . . , 60}. The optimal parameter vectors changed for

Fig. 8 :

 8 Fig. 8: Network MSD performance with model and clusters change.

Fig. 9 :Fig. 10 :

 910 Fig. 9: (Top) Ground truth cluster. (Bottom) Inferred clusters at steady-state of a single Monte Carlo run. From left to right: Stage 1, Stage 2, Stage 3.

  1 S ) (y(i) -ŷ(i))

Fig. 11 :

 11 Fig. 11: True temperatures and reconstructed ones at an unobserved node. µ LMS = 10 -5 , µ PLMS = µ LMSN = 10 -4 .

Fig. 12 :

 12 Fig. 12: U.S. temperature graph topology and learned clusters.

Fig. 13 :

 13 Fig. 13: True temperatures and reconstructed ones at an unobserved node. For clarity purposes, focus on the intervals (a) [4100, 4300] and (b) [8640, 8759].

TABLE I :

 I NMSE of different algorithms.

	Algorithm	NMSE
	KKF	0.1093
	KRR	0.0479
	Multitask diffusion LMS	0.1152
	Multitask diffusion PLMS	0.0090
	Multitask diffusion LMSN 0.0031

The work of C. Richard was funded in part by the PIA program under its IDEX UCAJEDI project (ANR-15-IDEX-0001) and by ANR under grant ANR-19-CE48-0002. The work of F. Hua was supported in part by China Scholarship Council and NSFC grant 61771396. The work of H. Wang was partially supported by NSFC under grants 61571365, 61671386 and National Key R&D Program of China (2016YFC1400200). The work of A. H. Sayed was supported in part by NSF grant CCF-1524250. Part of the results in this paper were presented at EUSIPCO [1] and ASILOMAR [2] conferences.