Roula Nassif
email: roula.nassif@aub.edu.lb

Stefan Vlaski
email: stefan.vlaski@epfl.ch

Cédric Richard
email: cedric.richard@unice.fr

Jie Chen
email: dr.jie.chen@ieee.org

Ali H Sayed
email: ali.sayed@epfl.ch

Multitask learning over graphs: An Approach for Distributed, Streaming Machine Learning

. Multitask learning is an approach to inductive transfer learning (using what is learned for one problem to assist in another problem) and helps improve generalization performance relative to learning each task separately by using the domain information contained in the training signals of related tasks as an inductive bias.

Several strategies have been derived within this community under the assumption that all data are available beforehand at a fusion center. However, recent years have witnessed an increasing ability to collect data in a distributed and streaming manner. This requires the design of new strategies for learning jointly multiple tasks from streaming data over distributed (or networked) systems. This article provides an overview of multitask strategies for learning and adaptation over networks. The working hypothesis for these strategies is that agents are allowed to cooperate with each other in order to learn distinct, though related tasks. The article shows how cooperation steers the network limiting point and how different cooperation rules allow to promote different task relatedness models. It also explains how and when cooperation over multitask networks outperforms non-cooperative strategies.

I. MULTITASK NETWORK MODELS

Consider a networked system consisting of a collection of N autonomous agents (sensors, classifiers, etc.) distributed over some geographic area and connected through a topology. The neighborhood of agent k is denoted by N k ; it consists of all agents that are connected to k by an edge-see Fig. 1 (left). A real-valued strongly convex and differentiable cost J k (w k ) is associated

with each agent k. The objective (or the task) at agent k is to estimate the parameter vector, w o k , of size M k × 1, that minimizes J k (w k ), namely,

w o k arg min wk J k (w k ). (1)
Depending on how the minimizers across the agents relate to each other, we distinguish between three categories of networks:

1) Single-task network: All costs J k (w k ) are minimized at the same location w o , namely, w o k = w o for all k -see Fig. 1 (left). 2) Clustered multitask network: The N agents are grouped into Q clusters C q (q = 1, . . . , Q) and, within each cluster C q , all the costs are minimized at the same location w o Cq , namely, w o k = w o Cq for all k ∈ C q -see Fig. Each agent k can solve (1) on its own. However, since the objectives across the network relate to each other, it is expected that by properly promoting these relationships, one may improve the network performance. In other words, it is expected that through cooperation among the agents, one may improve the network performance. One important question is how to design cooperative strategies that can lead to better performance than non-cooperative ones where each agent attempts to determine w o k on its own. This overview paper explains how multitask learning over graphs addresses this question.

Prior to multitask learning over graphs, there have been many works in the machine learning literature where learning multiple related tasks simultaneously has been considered [START_REF] Caruana | Multitask learning[END_REF]- [START_REF] Kato | Multi-task learning via conic programming[END_REF].

Multitask learning was shown, both empirically and theoretically, to improve performance relative to the traditional approach of learning each task separately. Depending on the machine learning application, several task relatedness models have been considered. For example, in [START_REF] Caruana | Multitask learning[END_REF], [START_REF] Evgeniou | Regularized multi-task learning[END_REF], the functions to be learned are assumed to share a common underlying representation. In [START_REF] Jacob | Clustered multi-task learning: A convex formulation[END_REF], it is assumed that the tasks are close to each other in some Hilbert space. Probabilistic based approaches, where a probability model capturing the relations between tasks is estimated simultaneously with functions corresponding to each task, have also been considered [START_REF] Zhang | A convex formulation for learning task relationships in multi-task learning[END_REF]. Also, graph-based approaches, where the relations between tasks are captured by an underlying graph, were also considered in the literature [START_REF] Chen | Graph-structured multitask regression and an efficient optimization method for general fused Lasso[END_REF], [START_REF] Kato | Multi-task learning via conic programming[END_REF]. All these works, however, assume that all data are available beforehand at a fusion center and propose batch-mode methods to solve multitask problems. Other existing works, such as [START_REF] Wang | Distributed multi-task learning[END_REF], consider distributed data setting. However, most of these works still require an architecture consisting of workers along with a master, where agents perform local computations followed by sending intermediate results to the master for further processing. Such solution methods are not fully distributed, which limits their range of practical applications. This paper, however, focuses on fully distributed solutions that avoid the need for central data aggregation or processing and instead rely on local computations and communication exchanges among neighborhoods. Besides providing distributed implementations, the solutions considered in this paper are able to learn continuously from streaming data. We start our exposition by describing a class of non-cooperative solutions that are able to respond in real time to streaming data.

Then, we explain how these solutions can be extended to handle multitask learning over graphs.

II. NONCOOPERATIVE LEARNING UNDER STREAMING DATA

Throughout this article, there is an explicit assumption that agents operate in the streaming data setting. That is, it is assumed that each agent k receives at each time instant i one instantaneous realization x k,i of a random data x k . The goal of agent k is to estimate the vector w o k that minimizes its risk function

J k (w k ) E xk Q k (w k ; x k )
, defined in terms of some loss function

Q k (•).
The expectation is computed over the distribution of the data x k . Agent k is particularly interested in solving the problem in the stochastic setting when the distribution of the data is generally unknown. This means that the risks J k (•) and their gradients ∇ wk J k (•) are unknown.

As such, approximate gradient vectors ∇ wk J k (•) will need to be employed. Doing so leads to the following stochastic gradient algorithm for solving (1):

w k,i = w k,i-1 -µ ∇ wk J k (w k,i-1 ), (2)
where w k,i is the estimate of w o k at iteration i and µ > 0 is a small step-size parameter. Resorting to the instantaneous realization x k,i of the random data x k , a common construction in the stochastic approximation theory is to employ the following gradient approximation at iteration i:

∇ wk J k (w k ) = ∇ wk Q k (w k ; x k,i ). (3)
We therefore focus in this paper on stochastic gradient algorithms, which are powerful iterative procedures for solving [START_REF] Caruana | Multitask learning[END_REF] in the streaming data. They enable continuous learning and adaptation in response to drifts in the location of the minimizers due to changes in the costs. We illustrate construction ( 2)-( 3) by considering scenarios from machine learning and adaptive filter theory.

Example 1. (Logistic regression network). Let γ k (i) = ±1 be a streaming sequence of (class) binary random variables and let h k,i be the corresponding streaming sequence of

M k × 1 real random (feature) vectors with R h,k = Eh k,i h k,i > 0. The processes {γ k (i), h k,i } are assumed
to be wide-sense stationary. In these problems, agent k seeks to estimate the vector w o k that minimizes the regularized logistic risk function [START_REF] Sayed | Adaptation, learning, and optimization over networks[END_REF]:

J k (w k ) = E ln 1 + e -γ k (i)h k,i wk + ρ 2 w k 2 , (4)
where ρ > 0 is a regularization parameter. Once w o k is found, γ k (i) = sign(h k,i w o k ) can then be used as a decision rule to classify new features. Using approximation (3), we obtain the following stochastic-gradient algorithm for minimizing (4):

w k,i = (1 -µρ)w k,i-1 + µγ k (i)h k,i 1 1 + e γ k (i)h k,i wk,i-1 . (5)
Example 2. (Mean-square-error (MSE) network). In such networks, each agent is subjected to streaming data {d k (i), u k,i } that are assumed to satisfy a linear regression model:

d k (i) = u k,i w o k + v k (i), (6)
for some unknown M k ×1 vector w o k to be estimated by agent k with v k (i) denoting a zero-mean measurement noise. For these networks, the risk function takes the form of an MSE cost [START_REF] Sayed | Diffusion strategies for adaptation and learning over networks[END_REF]:

J k (w k ) = 1 2 E d k (i) -u k,i w k 2 , (7)
which is minimized at w o k . The processes {u k,i , v k (i)} are zero-mean jointly wide-sense stationary with

: i) Eu k,i u ,i = R u,k > 0 if k = and zero otherwise; Ev k (i)v (i) = σ 2 v,k if k = and
zero otherwise; and iii) u k,i and v k (j) are independent of each other. Using approximation (3), we obtain the following stochastic-gradient algorithm:

w k,i = w k,i-1 + µu k,i (d k (i) -u k,i w k,i-1 ), (8)
which is the well-known least-mean-squares (LMS) algorithm [START_REF] Widrow | Adaptive Signal Processing[END_REF].

The use of the approximate gradient ∇ wk J k (•) instead of the true gradient

∇ wk J k (•) in (2)
introduces perturbations into the operation of the gradient descent iteration. This perturbation is referred to as the gradient noise defined as

s k,i (w k ) ∇ wk J k (w k ) -∇ wk J k (w k ).
The presence of this perturbation prevents the stochastic iterate w k,i from converging almost surely to the minimizer w o k when constant step-sizes are used. Some deterioration in performance will occur, and the iterate w k,i will instead fluctuate close to w o k . It is common in adaptive filtering and stochastic gradient optimization literatures to assess the size of these fluctuations by measuring their steady-state mean-square value [START_REF] Sayed | Adaptation, learning, and optimization over networks[END_REF]- [START_REF] Widrow | Adaptive Signal Processing[END_REF]. We therefore focus in this paper on highlighting the benefit of multitask learning on the network mean-square-deviation (MSD), which is defined as the steady-state average variance value:

MSD lim i→∞ 1 N N k=1 E w o k -w k,i 2 . ( 9
)
In the sequel, when discussing theoretical performance results, and to avoid excessive technicalities, it is sufficient to focus on MSE networks described in Example 2 and to assume that R u,k = R u and M k = M for all k 1 . In this way, the quality of the measurements, captured by the noise power σ 2 v,k , is allowed to vary across the network with some agents collecting noisier data than other agents. Assuming uniform regressors covariance allows us to quantify the improvement in performance that results from cooperation without biasing the results by the statistical nature of the regression data at the agents.

Performance result 1. Consider an MSE network running the non-cooperative algorithm [START_REF] Wang | Distributed multi-task learning[END_REF].

Assume further that R u,k = R u and M k = M for all k. Under these assumptions, and for 1 Performance results under more general conditions, such as allowing for space dependent covariances and lengths and for general second-order differentiable cost functions that are not necessarily quadratic, can also be found in [START_REF] Sayed | Adaptation, learning, and optimization over networks[END_REF] for algorithm (2), in [START_REF] Nassif | Learning over multitask graphs -Part I: Stability analysis[END_REF] for the strategy introduced in Sec. IV-A, and in [START_REF] Nassif | Adaptation and learning over networks under subspace constraints -Part II: Performance analysis[END_REF] for the strategies in Sec. V.

The MSD performance expressions in these works are derived under Lipschitz gradient vectors and Hessian matrices assumptions. It should be noted that the analyses in these works allow also to recover the Excess-Risk metric at agent k, which is defined as ER k limi→∞ E(J k (w k,i ) -J k (w o k ))-see, e.g., [9, p. 388-390]. Due to space limitations, we shall only focus on presenting MSD performance expressions.

sufficiently small step-sizes, the individual steady-state variance

MSD k lim i→∞ E w o k -w k,i 2
and the network MSD defined by (9) are given by [START_REF] Sayed | Diffusion strategies for adaptation and learning over networks[END_REF]:

MSD k = µM 2 • σ 2 v,k , MSD nc = µM 2 • 1 N N k=1 σ 2 v,k , (10)
where the superscript "nc" is used to indicate that the MSD expression is for the non-cooperative solution. First, observe that the performance is on the order of µ. The smaller µ is, the better the performance will be, but the slower the convergence toward w o k will be [START_REF] Sayed | Adaptation, learning, and optimization over networks[END_REF], [START_REF] Sayed | Diffusion strategies for adaptation and learning over networks[END_REF] (the same observation is valid for future expressions [START_REF] Zhou | A regularization framework for learning from graph data[END_REF] and [START_REF] Plata-Chaves | Distributed diffusion-based LMS for node-specific adaptive parameter estimation[END_REF] with convergence to W in (17) instead).

Second, observe that agents with noisier data will perform worse than agents with cleaner data.

However, since agents are observing data arising from similar or related models w o k , it is expected that an appropriate cooperation among agents can help enhance the network performance.

III. MULTITASK LEARNING FRAMEWORK

Depending on the application, several task relatedness models can be considered. For each model, an appropriate convex optimization problem is solved in a distributed and adaptive manner. This results in different multitask strategies, and therefore different cooperation rules between agents. Rather than describing each optimization problem in isolation, we begin by introducing a general problem, which will allow us to recover various multitask strategies as special cases.

Let W col{w 1 , . . . , w N } denote the collection of parameter vectors from across the network.

We consider the following global optimization problem for the multitask formulation:

W = arg min W J glob (W) = N k=1 J k (w k ) + η 2 R(W) subject to W ∈ Ω (11)
where R(•) is a convex regularization function promoting the relationships between the tasks, Ω is a closed convex set defining the feasible region of the parameter vectors, and η > 0 is a parameter controlling the importance of the regularization. The choice of the regularizer R(•)

and the set Ω depends on the prior information on how the multitask models relate to each other.

To illustrate how problem formulation [START_REF] Widrow | Adaptive Signal Processing[END_REF] can be used, we consider the following two examples that are multitask oriented.

Example 3. (Weather forecasting [START_REF] Nassif | Learning over multitask graphs -Part I: Stability analysis[END_REF]). Consider the network in Fig. 2 (left) consisting of N = 139 weather stations located across the United States and collecting daily measurements. Let h k,i denote the feature vector consisting of collected data (temperature, wind speed, dew point, etc.) 11) takes the form (4). However, it is expected that the decision rules {w o k } at neighboring stations would be similar since they are collecting features arising from similar statistical distributions. Moreover, the strength of similarity is expected to be inversely proportional to the physical distance between the stations. This gives rise to a weighted graph (with closest nodes connected by edges) and one may expect to improve the network performance by promoting the smoothness of {w o k } with respect to the underlying graph. The simplest possible term that encourages smoothness is the graph Laplacian regularizer S(W) defined further ahead in [START_REF] Kekatos | Distributed robust power system state estimation[END_REF]. By choosing R(W) = S(W) and Ω = R M N in [START_REF] Widrow | Adaptive Signal Processing[END_REF], one arrives at a multitask formulation for the weather forecasting application that takes into account the smoothness prior over the graph. This formulation and other possible formulations are solved in Sections IV and V.

Area 1 Area 1 Area 2 Area 2 Area 3 Area 3 Area 4 Area 4 B1 B1 B5 B5 B2 B2 B4 B4 B7 B7 B3 B3 B8 B8 B10 B10 B9 B9 B14 B14 B13 B13 B12 B12 B6 B6 B11 B11
Example 4. (Power system state monitoring [START_REF] Kekatos | Distributed robust power system state estimation[END_REF]). Consider Fig. 2 (right) illustrating an IEEE 14-bus power monitoring system partitioned into 4 areas, where each area comprises a subset of buses supervised by its own control center. The local state vectors (bus voltages) to be estimated at neighboring areas may partially overlap as the areas are interconnected. This is because each control center collects measurements related to the voltages across its local buses and voltages across the interconnection between neighboring centers. For example, Area 2 supervises buses 3, 4, 7, and 8. Since it collects current readings on lines (4, 5) and [START_REF] Kato | Multi-task learning via conic programming[END_REF][START_REF] Sayed | Adaptation, learning, and optimization over networks[END_REF], its state vector extends to buses 5 (supervised by Area 1) and 9 (supervised by Area 4). In other words, if we let w n denote the state of bus n, then the cost J 2 (•) at Area 2 will depend on the extended parameter vector w 2 = col{w 3 , w 4 , w 5 , w 7 , w 8 , w 9 }. However, since the parameter vectors at Areas 1 and 4 will be w 1 = col{w 1 , w 2 , w 5 } and w 4 = col{w 9 , w 10 , w 11 , w 14 }, respectively, consensus needs to be reached on the variable w 5 between Areas 2 and 1, and on the variable w 9 between Areas 2 and 4, while minimizing the individual cost J 2 (w 2 ) penalizing deviation from data models of the form y k = H k w k + v k where H k is the measurement matrix and v k is a zero-mean noise. Thus, distributed power state estimation can be formulated as problem [START_REF] Widrow | Adaptive Signal Processing[END_REF] with R(W) = 0, whereas the constraint set Ω in this case should be selected to promote consensus over the overlapped variables. In Section V-B, we explain how such problems can be solved.

Returning to the formulation [START_REF] Widrow | Adaptive Signal Processing[END_REF], observe that even though the aggregate cost N k=1 J k (w k ) is separable in w k , the cooperation between agents is necessary due to the coupling between the tasks through the regularization and the constraint. Note that, when solving problem [START_REF] Widrow | Adaptive Signal Processing[END_REF], agent k will be responsible for estimating w k (the k-th sub-vector of W = col{w 1 , . . . , w N }), which is generally different from w o k in (1), the actual objective at agent k. However, it is expected that accurate prior information will allow the designer to choose the regularizer R(•), the set Ω, and the strength η in a way that minimizes the distance between w k and w o k . Although some existing works use primal-dual methods [START_REF] Koppel | Proximity without consensus in online multiagent optimization[END_REF] to solve multitask estimation problems, we limit our exposition to the class of primal techniques (based on propagating and estimating the primal variable) that employ stochastic-gradient iterations. Extensive studies in the literature have shown that small step-sizes enable these strategies to learn well in streaming data settings. Due to the separability property of N k=1 J k (w k ), the multitask algorithms described in the sequel will have a common structure given by:

ψ k,i = w k,i-1 -µ ∇ wk J k (w k,i-1 ) w k,i = g k {ψ ,i } ∈Nk (12a) (12b)
The first step (12a) corresponds to the stochastic gradient step on the individual cost J k (•).

We refer to this step as the self-learning step-see Fig. 3. Compared with the non-cooperative strategy [START_REF] Thrun | Learning to Learn[END_REF], observe now that the result of the gradient descent step is ψ k,i , an intermediate estimate of w o k at iteration i. This step is followed by a social learning step (12b), which uses some function g k (•) of the neighborhood iterates. As we shall see in the next sections, the form of this function depends on the regularizer ηR(•) and the set Ω in [START_REF] Widrow | Adaptive Signal Processing[END_REF], both of which allow to promote the prior information on how the tasks w o k are related. The result of this second step is w k,i , the estimate of w o k , defined by (1), at iteration i. Since we are interested in a distributed setting, agents during social learning are only allowed to collect estimators from their local neighborhood N k -see Fig. 3.

In the sequel, we show how the formulation [START_REF] Widrow | Adaptive Signal Processing[END_REF] and the social learning step (12b) specialize for regularized (Sec. IV), subspace constrained (Sec. V), and clustered (Sec. VI) multitask estimation.

IV. REGULARIZED MULTITASK ESTIMATION

In this section, we focus on the regularization term R(W) in [START_REF] Widrow | Adaptive Signal Processing[END_REF] and its implications for the learning dynamics. In multitask learning (MTL), regularization is a widely used technique to promote task relationships. In most network applications, the underlying graph structure contains information about the relatedness among neighboring tasks. As such, when considering graph-based MTL applications, incorporating the graph structure into the regularization term is a reasonable and natural step. The smoothness model (under which the tasks are similar at neighboring vertices with the strength of similarity specified by the weight between them) will play a central role in our discussion. This smoothness property is often observed in real world applications (see, e.g., Example 3) and is rich enough to convey the main ideas behind MTL, as we will see in the sequel. We will examine two main questions: 1) How to incorporate graph-based priors into the regularizer? and 2) How does the resulting MTL algorithm behave?

A. Multitask estimation under smoothness

We assume that a symmetric, weighted adjacency matrix C is associated with the connected graph illustrated in Fig. 1 (right). If there is an edge connecting agents k and , then [C] k = c k > 0 reflects the strength of the relation between k and ; otherwise, [C] k = 0. These weights are usually dictated by the physics of the problem at hand-see, e.g., [START_REF] Akyildiz | Wireless sensor networks: A survey[END_REF], [START_REF] Grady | Discrete Calculus[END_REF]Ch. 4] for graph construction methods. We introduce the graph Laplacian L, which is a differential operator defined as L = diag{C1 N } -C. Assuming that the tasks have the same length, i.e., M k = M ∀k, the smoothness of W over the graph is measured in terms of a quadratic form of the Laplacian [START_REF] Zhou | A regularization framework for learning from graph data[END_REF]:

S(W) = W LW = 1 2 N k=1 ∈Nk c k w k -w 2 , (13)
where L = L⊗I M is an extended form of the graph Laplacian (defined in terms of the Kronecker product operator ⊗). The smaller S(W) is, the smoother the signal W on the graph is. Given that the weights are nonnegative, S(W) shows that W is smooth if nodes with a large c k on the edge connecting them have similar weight values {w k , w }. Therefore, in order to enforce the prior belief that the target signal W o = col{w o 1 , . . . , w o N } is smooth with respect to the underlying weighted graph, one may choose in [START_REF] Widrow | Adaptive Signal Processing[END_REF]:

R(W) = S(W), and Ω = R M N ( 14
)
Under this choice, the stochastic gradient algorithm for solving [START_REF] Widrow | Adaptive Signal Processing[END_REF] takes the following form:

W i = ψ i -µηLW i-1 (15)
where W i is the estimate of W at instant i, and ψ i = col{ψ k,i } N k=1 is the vector collecting the intermediate estimates ψ k,i in (12a) from across all agents. Since we expect ψ i to be an improved estimate compared to W i-1 , we propose to replace W i-1 in (15) by ψ i . By doing so, we obtain algorithm [START_REF] Nassif | Learning over multitask graphs -Part I: Stability analysis[END_REF] with the social learning step given by:

w k,i = ψ k,i -µη ∈Nk c k (ψ k,i -ψ ,i ) (16)
The substitution of W i-1 by ψ i is reminiscent of incremental-type arguments in gradient descent algorithms [START_REF] Bertsekas | A new class of incremental gradient methods for least squares problems[END_REF]. Analyses in the context of adaptation over networks show that substitutions of this type lead to enhanced network stability since they allow to preserve the stability of the agents after cooperation (see, e.g., [10, p. 160] for details). Regarding algorithm [START_REF] Grady | Discrete Calculus[END_REF], it follows that, when the spectral radius of the combination matrix I -µηL is equal to one, sufficiently small step-sizes ensuring the individual agents stability will also ensure the network stability 2 [START_REF] Nassif | Learning over multitask graphs -Part I: Stability analysis[END_REF].

Proximal based approaches are also proposed in [START_REF] Wang | Distributed stochastic multi-task learning with graph regularization[END_REF] to solve multitask problems under smoothness. However, these approaches require the evaluation of the proximal operator defined by [START_REF] Wang | Trend filtering on graphs[END_REF] of the risk Q k (•) at each iteration i, which can be computationally expensive.

B. Bias-variance tradeoff

We next consider the interesting question whether multitask learning is beneficial compared to

noncooperation. The answer to this inquiry requires i) studying the performance of algorithm [START_REF] Nassif | Learning over multitask graphs -Part I: Stability analysis[END_REF] relative to the actual agents objectives {w o k } and then ii) examining when the multitask implementation [START_REF] Nassif | Learning over multitask graphs -Part I: Stability analysis[END_REF] can lead to enhanced performance in comparison to the non-cooperative solution (2). 2 In this article, a network is said to be stable if the mean-square-error 1 N W -Wi 2 converges asymptotically to a bounded region of the order of the step-size.

Algorithm [START_REF] Nassif | Learning over multitask graphs -Part I: Stability analysis[END_REF] was studied in detail in [START_REF] Nassif | Learning over multitask graphs -Part I: Stability analysis[END_REF]. It was shown that the network MSD defined by ( 9) is mainly influenced by the sum of two factors, as explained further below. The first factor is the steady-state variance of algorithm [START_REF] Nassif | Learning over multitask graphs -Part I: Stability analysis[END_REF] with respect to the regularized solution W in [START_REF] Widrow | Adaptive Signal Processing[END_REF], namely, lim i→∞

1 N E W -W i 2 .
The second one is the bias or the average distance between the regularized solution W and the unregularized one W o , namely, 1 N W o -W 2 . By increasing the regularization strength η, the variance term is more likely to decrease while the bias term is more likely to increase. Understanding this bias-variance tradeoff is critical for understanding the behavior of regularized multitask algorithms.

We therefore describe in the following the bias-variance behavior of algorithm ( 16) by considering the expressions derived in [START_REF] Nassif | Learning over multitask graphs -Part I: Stability analysis[END_REF]. These expressions are useful for illustrating the concept of multitask learning. As we will see, instead of involving the vertex domain given by the entries {c k } of the adjacency matrix, these expressions involve the graph spectral information defined by the eigendecomposition of the Laplacian L. Because the Laplacian is a real symmetric matrix, it possesses a complete set of orthonormal eigenvectors. We denote them by {v 1 , . . . , v N }. For convenience, we order the set of real, non-negative eigenvalues of

L as 0 = λ 1 < λ 2 ≤ . . . ≤ λ N ,
where, since the network is connected, there is only one zero eigenvalue with corresponding [START_REF] Chung | Spectral Graph Theory[END_REF]. Therefore, the Laplacian can be decomposed as

eigenvector v 1 = 1 √ N 1 N
L = V ΛV where Λ = diag{λ 1 , . . . , λ N } and V = [v 1 , . . . , v N ].
Performance result 2. Consider an MSE network running the multitask algorithm [START_REF] Nassif | Learning over multitask graphs -Part I: Stability analysis[END_REF] with the second step given by [START_REF] Grady | Discrete Calculus[END_REF]. Assume further that R u,k = R u ∀k and that ρ(I -µηL) ≤ 1. Under these assumptions and for sufficiently small step-sizes and smooth signal W o , it is shown that [START_REF] Nassif | Learning over multitask graphs -Part I: Stability analysis[END_REF]:

lim i→∞ 1 N E W -W i 2 ≈ N m=1 ϕ(λ m ) (17)
where

ϕ(λ m ) = µ 2N N k=1 [v m ] 2 k σ 2 v,k   M q=1 λ u,q λ u,q + ηλ m   ( 18
)
with λ u,q the q-th eigenvalue of R u and [v m ] k the k-th entry of the eigenvector v m . For the bias term, it can be shown that [START_REF] Nassif | Learning over multitask graphs -Part I: Stability analysis[END_REF]:

W o -W 2 = N m=2 ζ(λ m ), (19)
where

ζ(λ m ) = ηλ m (R u + ηλ m I M ) -1 w o m 2 , ( 20
)
with w o m = (v m ⊗ I M )W o the m-th subvector of W o = (V ⊗ I M )W o corresponding to the eigenvalue λ m .
For the steady-state variance [START_REF] Zhou | A regularization framework for learning from graph data[END_REF], observe that it consists of the summation of N terms ϕ(λ m ), each one corresponding to an eigenvalue λ m of the Laplacian. The first one ϕ(λ 1 = 0) is independent of the regularization strength η. The remaining terms ϕ(λ m = 0) decrease when η increases. Therefore, when η increases, the variance decreases. From expression [START_REF] Wang | Distributed stochastic multi-task learning with graph regularization[END_REF], we observe that the bias tends to increase by increasing the regularization strength η. However, an interesting fact arises for smooth W o . To see this, we rewrite the regularizer in (13) as:

S(W o ) = (W o ) (Λ ⊗ I M )W o = N m=2 λ m w o m 2 , ( 21
)
where we used the fact that λ 1 = 0. Intuitively, given that λ m > 0 for m = 2, . . . , N , the above expression shows that W o is considered to be smooth over the graph if w o m 2 corresponding to large eigenvalue λ m is very small. That is, a smooth W o is mainly contained in [0, λ c ], i.e.,

w o m 2 ≈ 0 if λ m > λ c
, and the smoother W o is, the smaller λ c will be. In this case, the effective sum in ( 19) is over the first c N terms (corresponding to small eigenvalues λ m ) instead of N terms. We thus conclude that as long as W o is sufficiently smooth, moderate regularization strengths η in the range (0, ∞) exist such that the decrease in variance at these values of η will dominate the increase in bias. In other words, the MSD at these values of η will be less than the MSD at η = 0, which corresponds to the noncooperative mode of operation.

Observe from ( 16) that the social learning step following from the Laplacian regularization term (13) involves a single communication step at every stochastic gradient update. When multiple steps are allowed, it is reasonable to expect that performance can be improved. In the following, we show how such solution can be designed.

C. Graph spectral regularization

The main observation behind the introduction of this regularizer is that a smooth W o over a graph exhibits a special structure in the graph spectral domain (it is mainly contained in [0, λ c ], i.e., [START_REF] Smola | Kernels and regularization on graphs[END_REF]. Graph spectral regularization is used to leverage more thoroughly the spectral information and improve the multitask network performance. In this case, the network will aim at solving problem [START_REF] Widrow | Adaptive Signal Processing[END_REF] with Ω = R M N and R(•) properly selected in order to promote the prior information available on the structure of W o in the graph spectral domain. The following class of regularization functionals on graphs can be used for this purpose [START_REF] Smola | Kernels and regularization on graphs[END_REF], [START_REF] Nassif | A regularization framework for learning over multitask graphs[END_REF]:

w o m 2 ≈ 0 if λ m > λ c )
R(W) = W r(L)W = W (r(L) ⊗ I M )W (22)
where r(•) is some well-defined non-negative function on the spectrum σ(L) = {λ 1 , . . . , λ N } of L and r(L) is the corresponding matrix function defined as [23, p. 3]:

r(L) = V r(Λ)V = N m=1 r(λ m )v m v m . (23)
Construction [START_REF] Nassif | A regularization framework for learning over multitask graphs[END_REF] uses the Laplacian as a means to design regularization operators. Requiring a positive semi-definite regularizer r(L) imposes the constraint r(λ) ≥ 0 for all λ ∈ σ(L).

Replacing ( 23) into ( 22), we obtain -compare with the regularizer in [START_REF] Smola | Kernels and regularization on graphs[END_REF] to see how an extra degree of freedom is introduced in the multitask network design:

R(W) = W (r(Λ) ⊗ I M )W = N m=1 r(λ m ) w m 2 , (24)
where W = (V ⊗I M )W and w m = (v m ⊗I M )W. The regularization in (24) promotes a particular structure in the graph spectral domain. It strongly penalizes w m 2 for which the corresponding r(λ m ) is large. Thus, one prefers r(λ m ) to be large for those w m 2 that are small and vice versa. From the discussion following [START_REF] Smola | Kernels and regularization on graphs[END_REF], it is clear that, under smoothness, the function r(λ) must be chosen to be monotonically increasing in λ. One typical choice is r(λ) = λ S with S ≥ 1.

Example 5 further ahead illustrates for instance the benefit of using λ 3 instead of λ.

Assuming the regularizer r(L) in ( 22) can be written as an S-th degree polynomial of the Laplacian L, i.e., r(L) = S s=0 β s L s for some constants {β s } (or, equivalently, r(λ) = S s=0 β s λ s ), and following similar arguments that led to [START_REF] Grady | Discrete Calculus[END_REF], one arrives at the following social step (12b) [START_REF] Nassif | A regularization framework for learning over multitask graphs[END_REF]:

       ψ s k,i = β S-s ψ k,i + ∈Nk c k (ψ s-1 k,i -ψ s-1 ,i ), s = 1, . . . , S w k,i = ψ k,i -µηψ S k,i (25)
where ψ 0 k,i = β S ψ k,i . It requires S communication steps. The resulting algorithm ( 25) is distributed since at each step, each agent is only required to exchange information locally with its neighbors. Since S communication steps are required, agent k ends up collecting information from its S-hop neighborhood.

For more general r(λ) that are not necessarily polynomial in λ, one would like to benefit from the sparsity of the graph captured by L. As long as r(L) can be approximated by some lower order polynomial in L, say r(L) ≈ S s=0 β s L s , distributed implementations of the form (25) are possible-see [START_REF] Nassif | A regularization framework for learning over multitask graphs[END_REF]. Problems of this type have already been considered in graph filters design [START_REF] Shuman | Distributed signal processing via Chebyshev polynomial approximation[END_REF], [START_REF] Sandryhaila | Discrete signal processing on graphs: Frequency analysis[END_REF]. For instance, the work [START_REF] Shuman | Distributed signal processing via Chebyshev polynomial approximation[END_REF] proposes to locally approximate r(•) by a polynomial r(•) computed by truncating a shifted Chebyshev series expansion of r(•) on [0, λ N ]. When the Since the optimization problems in Sec. IV-A and IV-C are the same with L in (13) replaced by r(L) in ( 22), the multitask strategy [START_REF] Sandryhaila | Discrete signal processing on graphs: Frequency analysis[END_REF] will behave in a similar manner as [START_REF] Grady | Discrete Calculus[END_REF]. Particularly, the bias-variance tradeoff discussion continues to apply, and expressions ( 17)-( 20) continue to hold with λ m on the RHS of ( 18) and ( 20) replaced by the function r(λ m ). By replacing λ m on the RHS of (20) by r(λ m ), one may directly observe the consequence of this regularizer on the bias term [START_REF] Wang | Distributed stochastic multi-task learning with graph regularization[END_REF], which can be made now close to zero (by choosing in the smoothness case, for example, r(λ m ) ≈ 0 if λ m ∈ [0, λ c ] and β m > 0 otherwise).

Example 5. (Graph spectral filtering). Consider the MSE network example in Fig. 4 and assume uniform data profile, i.e., R u,k = R u and σ 2 v,k = σ 2 v ∀k. In the left plot, we illustrate the entries of the tasks w o k , which are smooth over the underlying graph. In the middle plot, we illustrate the behavior of the previously described algorithms in the graph spectral domain. The top plot represents the behavior of the network output W for three different choices of regularizer r(λ) = {0, λ, λ 3 }. The bottom plot represents the behavior of the steady-state variance [START_REF] Bertsekas | A new class of incremental gradient methods for least squares problems[END_REF] with the eigenvalue λ m replaced by the function r(λ m ). Observe how the regularizer r(λ) = λ 3 penalizes low-eigenvalues less than r(λ) = λ, and consequently preserves all the signal components w m .

Observe further the graph low-pass filtering behavior [START_REF] Shuman | Distributed signal processing via Chebyshev polynomial approximation[END_REF], [START_REF] Sandryhaila | Discrete signal processing on graphs: Frequency analysis[END_REF]. Small eigenvalues λ m correspond to low frequencies, W = (V ⊗ I M )W corresponds to the graph Fourier transform, and w m = (v m ⊗ I M )W corresponds to the m-th frequency content of W. It can be shown that the m-th frequency content of the output can be bounded as w m ≤ λu,max λu,max+ηr(λm) w o m in terms of the m-th frequency content of the input W o where λ u,max is the maximum eigenvalue of R u (see [START_REF] Nassif | Learning over multitask graphs -Part I: Stability analysis[END_REF]). Since r(λ) is monotonically increasing in λ, for fixed η, as λ increases, the ratio decreases. Therefore, the network output W can be interpreted as the output of a low-pass graph filter applied to the signal W o . A similar behavior arises for the steady-state variance. For fixed η, as λ m increases, the variance at the m-th frequency, i.e., ϕ(λ m ), decreases, and for fixed λ m , as η increases, ϕ(λ m ) decreases. The regularizer r(λ) controls the shape of the filter and the strength η controls the sharpness. The non-cooperative solution (η = 0) corresponds to an all-pass graph filter. In the right plot, we illustrate the bias-variance tradeoff in the case of r(λ) = λ 3 .

Returning to the diagram in Fig. 3, observe that the self-learning step corresponds to the inference step where agent k estimates w o k from streaming data x k,i and that the social-learning step corresponds to the graph filtering step where the agents collaborate in order to perform spatial filtering and reduce the effect of the noise on the network MSD defined by [START_REF] Sayed | Adaptation, learning, and optimization over networks[END_REF]. These steps are performed simultaneously. Therefore, multitask learning over networks allows to blend real-time adaptation with graph (spatial) filtering.

D. Non-quadratic regularization

Non-quadratic regularization has been also considered in the literature [START_REF] Koppel | Proximity without consensus in online multiagent optimization[END_REF], [START_REF] Nassif | Proximal multitask learning over networks with sparsity-inducing coregularization[END_REF], [START_REF] Hallac | Network Lasso: Clustering and optimization in large graphs[END_REF]. This scenario will induce non-linearities in the social learning step (12b). In this case, multitask algorithms are derived in order to solve problem [START_REF] Widrow | Adaptive Signal Processing[END_REF] with

Ω = R M N and R(W) = N k=1 ∈Nk ρ k h k (w k , w ) (26)
where h k : R M × R M → R is a convex cost function associated with the link (k, ). In general, this function is used to enforce some constraints on the pairs of variables across an edge. Observe that [START_REF] Nassif | Proximal multitask learning over networks with sparsity-inducing coregularization[END_REF], by allowing arbitrary distance measures h k (•, •), is a generalization of the previously employed quadratic regularization. In fact, setting h k (w k , w ) = w kw 2 recovers (13).

Examples of other typical choices are the 2 -norm regularizer h k (w k , w ) = w kw and the

1 -norm regularizer h k (w k , w ) = w k -w 1 .
Instead of encouraging global smoothness, these sparsity-based regularizers can adapt to heterogeneity in the level of smoothness of the tasks w o k across nodes [START_REF] Wang | Trend filtering on graphs[END_REF]. Such heterogeneity is observed for instance in the problem of predicting housing prices [START_REF] Hallac | Network Lasso: Clustering and optimization in large graphs[END_REF]. In this problem, the objective at each node (house) in a graph (where neighboring houses are connected by edges) is to learn the weights w o k of a regression model (examples of features are number of bedrooms, square footage, etc.) to estimate the price. Due to location-based factors (such as distance to highway) that are often unknown a priori and, therefore, cannot be incorporated as features, similar houses in different, though close (neighbors), locations can have drastically different prices, i.e., drastically different w o k . The objective in this case is to encourage neighboring houses that share common models to cooperate without being influenced by the misleading information of neighbors sharing different models, i.e., perform automatic clustering. To do so, the authors in [START_REF] Hallac | Network Lasso: Clustering and optimization in large graphs[END_REF] propose to solve the network Lasso problem, i.e., problem [START_REF] Widrow | Adaptive Signal Processing[END_REF] with 2 -norm regularizer in [START_REF] Nassif | Proximal multitask learning over networks with sparsity-inducing coregularization[END_REF]. The rationale behind this choice is that 2norm encourages group sparsity, i.e., consensus across an edge w k = w . On the other hand, the 1 -norm regularizer is used in [START_REF] Nassif | Proximal multitask learning over networks with sparsity-inducing coregularization[END_REF] to promote the prior that the parameter vectors at neighboring nodes have a large number of similar entries and a small number of distinct entries. The weight ρ k ≥ 0 in (26) associated with the link (k, ) aims at locally adjusting the regularization strength.

It is usually dictated by the physics of the problem at hand. For primal adaptive techniques and due to the non-differentiability of the regularizers, proximal gradient methods can be used to solve [START_REF] Widrow | Adaptive Signal Processing[END_REF]. Assuming ρ k = ρ k and h k (w k , w ) = w kw 1 , one may arrive to a multitask algorithm of the form [START_REF] Nassif | Learning over multitask graphs -Part I: Stability analysis[END_REF] with the social learning step (12b) given by (see the derivations in [START_REF] Nassif | Proximal multitask learning over networks with sparsity-inducing coregularization[END_REF]):

w k,i = prox ηµ gk,i (ψ k,i ) (27)
where prox γg (w ) denotes the proximal operator of the function g(w):

prox γg (w ) = argmin w∈R M g(w) + 1 2γ w -w 2 , (28)
and where the function g k,i : R M → R is given by g k,i (w k ) = ∈Nk ρ k h k (w k , ψ ,i ). Notice that the proximal operator in [START_REF] Hallac | Network Lasso: Clustering and optimization in large graphs[END_REF] needs to be evaluated at each iteration. For the weighted sum of 1 -regularizer, a closed form expression can be found in [START_REF] Nassif | Proximal multitask learning over networks with sparsity-inducing coregularization[END_REF].

V. MULTITASK ESTIMATION UNDER SUBSPACE CONSTRAINTS

Besides regularized-based algorithms, projection-based algorithms have received considerable attention in the literature of deterministic [START_REF] Bertsekas | A new class of incremental gradient methods for least squares problems[END_REF], [START_REF] Lorenzo | Distributed signal recovery based on in-network subspace projections[END_REF]- [START_REF] Mota | Distributed optimization with local domains: Applications in MPC and network flows[END_REF] and stochastic [START_REF] Sayed | Adaptation, learning, and optimization over networks[END_REF], [START_REF] Sayed | Diffusion strategies for adaptation and learning over networks[END_REF], [START_REF] Nassif | Adaptation and learning over networks under subspace constraints -Part II: Performance analysis[END_REF]- [START_REF] Sahu | CIRF E: A distributed random fields estimator[END_REF] optimization. The objective in this case is to design distributed networks that are able to project onto low-dimensional subspaces while minimizing the individual costs, i.e., solve problems of the form [START_REF] Widrow | Adaptive Signal Processing[END_REF] with [START_REF] Lorenzo | Distributed signal recovery based on in-network subspace projections[END_REF], [START_REF] Nassif | Adaptation and learning over networks under subspace constraints -Part II: Performance analysis[END_REF]: R(W) = 0, and Ω = Range(U) [START_REF] Lorenzo | Distributed signal recovery based on in-network subspace projections[END_REF] where Range(•) denotes the range space operator and U is an M t × P full-column rank matrix with M t = N k=1 M k and P M t . The reader will soon realize that consensus-type problems are instances of this formulation. Also, multitask estimation under smoothness can benefit from this formulation: as explained earlier in Sec. IV-C, when the first c eigenvectors of the Laplacian are available, the designer can project onto Range(U) with U = [v 1 , . . . , v c ] ⊗ I M instead of using regularization.

Let P U = U(U U) -1 U denote the projection onto the range space of U. Assuming that the network topology and the signal subspace U are such that the following feasibility problem:

find A such that A U = U, U A = U , ρ(A -P U ) < 1, [A] k = 0, if / ∈ N k and = k (30)
admits at least one solution, one may arrive to a multitask strategy of the form [START_REF] Nassif | Learning over multitask graphs -Part I: Stability analysis[END_REF] with the social learning step (12b) given by [START_REF] Nassif | Adaptation and learning over networks under subspace constraints -Part II: Performance analysis[END_REF]:

w k,i = ∈Nk A k ψ ,i (31)
where A k = [A] k is the (k, )-th block (of size M k × M ) of the N × N block matrix A. A matrix A satisfying the constraints in ( 30) is semi-convergent [START_REF] Lorenzo | Distributed signal recovery based on in-network subspace projections[END_REF], [START_REF] Nassif | Adaptation and learning over networks under subspace constraints -Part II: Performance analysis[END_REF]. Particularly, it holds that:

lim i→∞ A i = P U . ( 32
)
The first two constraints in [START_REF] Nedic | Distributed subgradient methods for multi-agent optimization[END_REF] state that the P columns of U are right and left eigenvectors of A associated with the eigenvalue 1. Together with these two constraints, the third constraint in [START_REF] Nedic | Distributed subgradient methods for multi-agent optimization[END_REF] ensures that A has P eigenvalues at one, and that all other eigenvalues are strictly less than one in magnitude. The last constraint in (30) corresponds to the sparsity constraint which characterizes the network topology and ensures local exchange of information at each instant i.

Before explaining how some typical choices of U lead to well-studied distributed inference problems, we note that the distributed algorithm [START_REF] Mota | Distributed optimization with local domains: Applications in MPC and network flows[END_REF] has an attractive property: in the small stepsize regime, the iterates generated by [START_REF] Mota | Distributed optimization with local domains: Applications in MPC and network flows[END_REF] achieve the steady-state performance of the following gradient projection algorithm [START_REF] Nassif | Adaptation and learning over networks under subspace constraints -Part II: Performance analysis[END_REF]:

W i = P U W i -µcol ∇ wk J k (w k,i-1 ) N k=1 (33)
which is centralized since, at each instant i, agent k needs to send its estimate ψ k,i in (12a) to a fusion center, which performs the projection, and then sends the result w k,i back to the agent.

Performance result 3. Consider an MSE network running algorithm [START_REF] Nassif | Learning over multitask graphs -Part I: Stability analysis[END_REF] with the social step (12b)

given by [START_REF] Mota | Distributed optimization with local domains: Applications in MPC and network flows[END_REF] with A = [A k ] satisfying the constraints in [START_REF] Nedic | Distributed subgradient methods for multi-agent optimization[END_REF]. Assume that the network is seeking W o ∈ Range(U). Assume further that U = U ⊗ I M where U = [u 1 , . . . , u P ] is semi-orthogonal, and that R u,k = R u and M k = M for all k. Under these assumptions, and for sufficiently small step-sizes, the network MSD defined by ( 9) is given by [START_REF] Nassif | Adaptation and learning over networks under subspace constraints -Part II: Performance analysis[END_REF]:

MSD = µM 2N P m=1 N k=1 [u m ] 2 k σ 2 v,k . (34)
Notice that the projection framework will not induce bias in the estimation. This is because W o ∈ Range(U), and, therefore, the vector W in ( 11) is equal to W o , the network objective.

Moreover, the benefit of cooperation can be readily seen by assuming uniform variances σ 2 v,k = σ 2 v for all k. In this case, comparing [START_REF] Plata-Chaves | Distributed diffusion-based LMS for node-specific adaptive parameter estimation[END_REF] with [START_REF] Sayed | Diffusion strategies for adaptation and learning over networks[END_REF] in the non-cooperative case, we conclude that MSD = ( P /N )MSD nc where P /N 1. Therefore, the cooperative strategy outperforms the non-cooperative one by a factor of N/ P .

A. Single-task estimation

In single-task estimation, the agents are seeking a common minimizer w o -see Fig. 1 (left). This problem is encountered in many applications. Examples include target localization and distributed sensing (see, e.g., [START_REF] Sayed | Diffusion strategies for adaptation and learning over networks[END_REF]). Single-task estimation can be recast in the form [START_REF] Widrow | Adaptive Signal Processing[END_REF] where R(•) and Ω are chosen according to [START_REF] Lorenzo | Distributed signal recovery based on in-network subspace projections[END_REF] with U = 1 √ N (1 N ⊗I M ). Several algorithms for solving such consensustype problems have been proposed in the literature, including incremental [START_REF] Bertsekas | A new class of incremental gradient methods for least squares problems[END_REF], consensus [START_REF] Nedic | Distributed subgradient methods for multi-agent optimization[END_REF],

and diffusion [START_REF] Sayed | Adaptation, learning, and optimization over networks[END_REF], [START_REF] Sayed | Diffusion strategies for adaptation and learning over networks[END_REF] strategies. Due to lack of space, we will describe only the class of diffusion strategies, which can be written in the form [START_REF] Nassif | Learning over multitask graphs -Part I: Stability analysis[END_REF] with the social step (12b) given by:

w k,i = ∈Nk a k ψ ,i (35)
where a k corresponds to the (k, )-th entry of an N × N doubly-stochastic matrix A satisfying:

a k ≥ 0, N =1 a k = 1, N k=1 a k = 1, and a k = 0 if / ∈ N k . (36)
Several rules for selecting locally these combination coefficients have been proposed in the literature, such as the Metropolis rule and Laplacian rule; see, e.g., [START_REF] Sayed | Adaptation, learning, and optimization over networks[END_REF]. Observe that step [START_REF] Sahu | CIRF E: A distributed random fields estimator[END_REF] can be written in the form of [START_REF] Mota | Distributed optimization with local domains: Applications in MPC and network flows[END_REF] with A k = a k I M and A = A ⊗ I M , and that the resulting matrix A will satisfy the constraints in (30) over a strongly connected network.

B. Multitask estimation with overlapping parameter vectors

It is assumed that the individual costs J k (•) depend only on a subset of the components of a global parameter vector w = [w 1 , . . . , w M ] ∈ R M ×1 [START_REF] Mota | Distributed optimization with local domains: Applications in MPC and network flows[END_REF], [START_REF] Alghunaim | Distributed coupled multi-agent stochastic optimization[END_REF]- [START_REF] Sahu | CIRF E: A distributed random fields estimator[END_REF]. This situation is observed in Example 4 where the network global parameter vector is w = col{w n } 14 n=1 and where the states w k to be estimated at neighboring areas partially overlap. It can be verified that this problem can also be recast in the form [START_REF] Lorenzo | Distributed signal recovery based on in-network subspace projections[END_REF] with U properly selected. To solve this consensus-type problem, and motivated by the single-task diffusion strategies, the works [START_REF] Alghunaim | Distributed coupled multi-agent stochastic optimization[END_REF], [START_REF] Plata-Chaves | Distributed diffusion-based LMS for node-specific adaptive parameter estimation[END_REF] propose the following algorithm. Assume agent k is interested in estimating the entry w n of w and let N n k denote the set of neighbors of k that are interested also in estimating w n . In order to reach consensus on w n , agent k assigns to its entry w n a set of non-negative coefficients {a n k } satisfying

a n k = 0 if / ∈ N n k , ∈N n k a n k = 1, ∈N n k a n k = 1, (37)
and performs the following convex combination:

w n k,i = ∈N n k a n k ψ n ,i (38)
where ψ n ,i is the entry of the M × 1 intermediate estimate ψ ,i (obtained from (12a)) corresponding to the variable w n and w n k,i is the estimate of w n at node k and instant i. It can also be verified that solution [START_REF] Zhao | Distributed clustering and learning over networks[END_REF] can be written in the form [START_REF] Mota | Distributed optimization with local domains: Applications in MPC and network flows[END_REF] with the block A k properly selected.

For MSE networks, a recursive least-squares (RLS) approach is proposed in [START_REF] Sahu | CIRF E: A distributed random fields estimator[END_REF] to solve overlapping multitask estimation. In general, second-order gradient methods enjoy faster convergence rates than first-order methods at the expense of increasing the computational complexity.

VI. CLUSTERED MULTITASK ESTIMATION

Now we move into explaining how clustered multitask estimation can be solved. Clustered multitask learning was first considered in [START_REF] Jacob | Clustered multi-task learning: A convex formulation[END_REF] within the machine learning community. Then, it was extended to adaptation and learning over networks in the work [START_REF] Chen | Multitask diffusion adaptation over networks[END_REF]. As we shall see, clustered multitask estimation merges subspace constraints with regularization. Let M k = M for all k. In clustered multitask networks, agents within a cluster C q are interested in estimating the same vector w o Cq -see Fig. 1 (middle). Without loss of generality, we index agents according to their cluster indexes such that agents from the same cluster will have consecutive indexes. Let N q denote the number of agents in cluster C q . Since agents within C q need to reach a consensus on w o Cq , clustered multitask estimation problems can be recast in the form [START_REF] Widrow | Adaptive Signal Processing[END_REF] with:

Ω = Range(U), U = diag 1 N q (1 Nq ⊗ I M ) Q q=1 (39)
Therefore, the cluster consensus step takes the form [START_REF] Mota | Distributed optimization with local domains: Applications in MPC and network flows[END_REF] with A = A⊗I M and A = diag{A q } Q q=1 where the N q × N q blocks A q are chosen according to the constraints in [START_REF] Nedic | Distributed subgradient methods for multi-agent optimization[END_REF]; one typical choice is doubly-stochastic blocks. The resulting N × N matrix A = [a k ] will satisfy:

a k ≥ 0, A1 N = 1 N , 1 N A = 1 N , and
a k = 0 if / ∈ N k ∩ C(k), (40)
where N k ∩ C(k) denotes the neighboring nodes of k that are inside its cluster. The choice of the regularizer R(W) in ( 11) depends on the prior information on how the models across the clusters relate to each other. One typical choice is [START_REF] Nassif | Proximal multitask learning over networks with sparsity-inducing coregularization[END_REF], [START_REF] Chen | Multitask diffusion adaptation over networks[END_REF]:

R(W) = N k=1 ∈Nk\C(k) ρ k h k (w k , w ) (41)
where N k \ C(k) denotes neighboring nodes of k that are outside its cluster and h k (w k , w ) is a cost associated with the inter-cluster link (k, ). This function is used to enforce some constraints on the pairs of variables across an inter-cluster edge. Examples are h k (w k , w ) = w kw 2

to enforce graph smoothness [START_REF] Chen | Multitask diffusion adaptation over networks[END_REF] and h k (w k , w ) = w kw 1 to enforce sparsity priors [START_REF] Nassif | Proximal multitask learning over networks with sparsity-inducing coregularization[END_REF].

Clustered multitask algorithms have in general the structure [START_REF] Nassif | Learning over multitask graphs -Part I: Stability analysis[END_REF] with step (12b) given by:

φ k,i = ∈Nk∩C(k) a k ψ ,i w k,i = g k (φ k,i , {φ ,i } ∈Nk\C(k) ) (42a) (42b)
In this algorithm, the self-learning step (12a) is followed by an intra-cluster social learning step (42a) where node k receives the intermediate estimates ψ ,i from its intra-cluster neighbors N k ∩C(k) and combines them in a convex manner through the coefficients {a k } in [START_REF] Plata-Chaves | Heterogeneous and multitask wireless sensor networks[END_REF] to obtain the intermediate value φ k,i . The second step (42b) is an inter-cluster social learning step where agent k receives the intermediate estimates {φ ,i } from its neighbors that are outside its cluster N k \ C(k) and combines them properly using the function g k (•) to obtain w k,i . This step helps to incorporate the available prior information on how the models across the clusters are related into the adaptation mechanism. The function g k (•) depends on the regularizer R(•). For example, for 1 -norm co-regularizers h k (with ρ k = ρ k ), one may arrive to an inter-cluster learning step (42b) given by w k,i = prox ηµ gk,i (φ k,i ) with g k,i (w k ) = ∈Nk\C(k) ρ k h k (w k , φ ,i ) [START_REF] Nassif | Proximal multitask learning over networks with sparsity-inducing coregularization[END_REF].

VII. CONCLUSION

In this article, we explained how prior knowledge about tasks relationships can be incorporated into the adaptation mechanism and how different priors yield different multitask strategies. It then follows that choosing the optimal strategy for a given problem is equivalent to choosing the task relatedness model which best fits the underlying problem. Choosing the most practically viable solution then balances further this model fit against computational and communication constraints.

There are several other aspects and strategies for multitask learning over graphs that were not covered in this article due to space limitations. For instance, we only focused on multitask networks endowed with parameter estimation tasks. However, distributed detection was also considered from a multitask perspective (see, e.g., [START_REF] Teklehaymanot | Robust diffusion-based unsupervised object labelling in distributed camera networks[END_REF]). Online network clustering was also considered. The objective in this case is to design diffusion networks that are able to adapt their combination coefficients in [START_REF] Sahu | CIRF E: A distributed random fields estimator[END_REF] in order to exclude harmful neighbors sharing distinct tasks [START_REF] Zhao | Distributed clustering and learning over networks[END_REF], [START_REF] Chen | Diffusion LMS over multitask networks[END_REF]. Readers can refer to [START_REF] Plata-Chaves | Heterogeneous and multitask wireless sensor networks[END_REF] to have a list of other literature works that are multitask oriented.

Multitask learning over graphs is worth exploring further, as there are many potential ideas to build on. For instance, the expressions show the sensitivity of the results to the underlying graph structure. It would be useful to infer the entries c k of the adjacency matrix in [START_REF] Grady | Discrete Calculus[END_REF] simultaneously with the self-learning step (12a). This leads to learning the relations between the tasks simultaneously with the tasks. Automatically determining the optimal regularization strength η and allowing edge regularizers h k (w k , w ) beyond just the 1 -norm constitute also clear extensions. Finally, we believe that the number of multitask learning applications in "distributed, streaming machine learning" is vast, and hope to witness increased utilization of the algorithms and theoretical results established in the domain of "learning and adaptation over networks".

7 Fig. 1 .

 71 Fig. 1. Network models. (Left) Single-task network. (Middle) Clustered multitask network. (Right) Multitask network.

3 )

 3 1 (middle). Similarities or relationships may exist among the distinct minimizers {w o Cq }. Multitask network: The individual costs are minimized at distinct, though related, locations {w o k } -see Fig. 1 (right).

Fig. 2 .

 2 Fig. 2. Examples of multitask applications. (Left) Weather forecasting. (Right) Distributed power system monitoring.

Fig. 3 .

 3 Fig. 3. A common diagram for the multitask strategies described in this work. The structure involves two main steps: i) a self-learning step (12a), and ii) a social learning step (12b).

Fig. 4 .

 4 Fig. 4. Illustrative example for spectral regularization. (Left) Estimation under smoothness. (Middle) Behavior of algorithm (12) in the graph spectral domain with w m = (v m ⊗ IM )W . (Right) Bias-variance tradeoff for r(λ) = λ 3 .


R A Borsoi

C Richard

A Ferrari

J Chen

J C M Bermudez

Online Graph-Based Change Point Detection in Multiband Image Sequences

Keywords: Hyperspectral images, change detection, graphs, multitemporal, superpixels

The automatic detection of changes or anomalies between multispectral and hyperspectral images collected at different time instants is an active and challenging research topic. To effectively perform change-point detection in multitemporal images, it is important to devise techniques that are computationally efficient for processing large datasets, and that do not require knowledge about the nature of the changes. In this paper, we introduce a novel online framework for detecting changes in multitemporal remote sensing images. Acting on neighboring spectra as adjacent vertices in a graph, this algorithm focuses on anomalies concurrently activating groups of vertices corresponding to compact, well-connected and spectrally homogeneous image regions. It fully benefits from recent advances in graph signal processing to exploit the characteristics of the data that lie on irregular supports. Moreover, the graph is estimated directly from the images using superpixel decomposition algorithms. The learning algorithm is scalable in the sense that it is efficient and spatially distributed. Experiments illustrate the detection and localization performance of the method.

I. INTRODUCTION

The increasing availability of multitemporal multi-and hyperspectral devices allows for a detailed analysis of the evolution of a scene over time. This can potentially benefit different applications, ranging from agricultural and forestry monitoring, natural disaster and urban landscape analysis [START_REF] Liu | A review of change detection in multitemporal hyperspectral images: Current techniques, applications, and challenges[END_REF], to surveillance and other industry problems [START_REF] Banerjee | Hyperspectral video for illumination-invariant tracking[END_REF], [START_REF] Stefansson | Hyperspectral time series analysis: hyperspectral image data streams interpreted by modeling known and unknown variations[END_REF]. Among these applications, the detection of changes or anomalies between images acquired at different time instants is an active and challenging research topic [START_REF] Liu | A review of change detection in multitemporal hyperspectral images: Current techniques, applications, and challenges[END_REF], [START_REF] Jianya | A review of multi-temporal remote sensing data change detection algorithms[END_REF]. The earliest change detection (CD) works consisted of post-classification algorithms, which compare supervised classification results of the two images [START_REF] Liu | A review of change detection in multitemporal hyperspectral images: Current techniques, applications, and challenges[END_REF], [START_REF] Jianya | A review of multi-temporal remote sensing data change detection algorithms[END_REF]. Although conceptually simple, those methods depend on training data, whose availability is often limited. This motivated the development of unsupervised techniques, which usually consists in analysing some kind of distance between feature-based representations of the pixels in both images [START_REF] Liu | A review of change detection in multitemporal hyperspectral images: Current techniques, applications, and challenges[END_REF]. Other popular approaches search for a transformation of the data to a lower-dimensional feature space that highlights the changes (e.g., multivariate alteration detection (MAD), temporal PCA, or unmixing) [START_REF] Liu | A review of change detection in multitemporal hyperspectral images: Current techniques, applications, and challenges[END_REF].

The work of R. Borsoi Traditionally, most CD algorithms were focused on processing a single pair of images (i.e., the bitemporal problem). However, those methods fail to explore the information contained in longer data streams, which has recently become the source of attention. For instance, [START_REF] Bertoluzza | A novel framework for bi-temporal change detection in image time series[END_REF] proposed to improve the performance of binary CD by constraining detection results to be consistent across all closed loops of subsequences sampled from the observed time series. CD algorithms known a priori have also been applied to new pairs of images using domain adaptation techniques [START_REF] Zanotta | An adaptive semisupervised approach to the detection of user-defined recurrent changes in image time series[END_REF]. In [START_REF] Vaduva | A latent analysis of earth surface dynamic evolution using change map time series[END_REF], different features were computed in local spatial windows centered at each pixel, for the whole time series. These temporal vectors of features were then classified using methods such as, e.g., K-means. More applied work often considers precise season-trend parametric models over, e.g., vegetation indices and detect changes by examining the model parameters [START_REF] Verbesselt | Near real-time disturbance detection using satellite image time series[END_REF].

Other algorithms were devised to operate online, such as [START_REF] Salmon | Land cover change detection using the internal covariance matrix of the extended Kalman filter over multiple spectral bands[END_REF], which detects changes on vegetation index data by monitoring the state covariance matrix of the extended Kalman filter. Another approach used density ratio estimation (learned from training data) to detect changes [START_REF] Anees | A relative density ratio-based framework for detection of land cover changes in MODIS NDVI time series[END_REF]. In [START_REF] Grobler | Using Page's cumulative sum test on MODIS time series to detect land-cover changes[END_REF] a cumulative sum test is used to detect change points in time series data. However, the pre-and post-change densities have to be estimated a priori. Recent work also considered multitemporal unmixing to separate the observed images into time series of the material signatures (endmembers) and their fractional abundances for each pixel [START_REF] Thouvenin | Online unmixing of multitemporal hyperspectral images accounting for spectral variability[END_REF], [START_REF] Borsoi | Kalman filtering and expectation maximization for multitemporal spectral unmixing[END_REF]. This can be used to perform CD in the lower dimensional space [START_REF] Liu | A review of change detection in multitemporal hyperspectral images: Current techniques, applications, and challenges[END_REF]. Other methods based on, e.g., total-variation [START_REF] Frecon | On-the-fly approximation of multivariate total variation minimization[END_REF] and low-rank and sparse decomposition [START_REF] Candès | Robust principal component analysis?[END_REF] can also be considered to perform CD. However, reconciling a low computational complexity with flexibility to incorporate a priori knowledge about the nature of the changes can be difficult in these techniques.

Despite the significant interest in this problem, there is still a need for algorithms that operate online and explore the temporal and the intrinsic structure of the data, while also being unsupervised and computationally efficient. In this work, we propose an online framework for detecting changes in multitemporal multiband images. We assume that changes occur concurrently in groups of connected, spectrally homogeneous pixels, which are then represented as adjacent vertices in a graph constructed from the superpixel representation of the images. Instead of considering graph regularizations [START_REF] Yokoya | Multisensor coupled spectral unmixing for time-series analysis[END_REF] or graph-based representations of the images learned from the data [START_REF] Héas | Modeling trajectory of dynamic clusters in image time-series for spatio-temporal reasoning[END_REF], which lead to computationally inefficient and offline solutions, the proposed method fully benefits from recent advances in graph signal processing to better exploit the characteristics of the data that lie on irregular supports [START_REF] Shuman | The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains[END_REF]. Moreover, differently from previous bitemporal CD approaches that use superpixels to compute features to classification algorithms [START_REF] Gong | Superpixel-based difference representation learning for change detection in multispectral remote sensing images[END_REF] or as an ad hoc approach to introduce spatial homogeneity in CD algorithms [START_REF] Wu | Superpixel-based unsupervised change detection using multi-dimensional change vector analysis and SVMbased classification[END_REF], we use the superpixel decomposition to generate a graph directly from the observed images. This allows us to exploit theoretically-grounded tools from graph signal processing. The learning algorithm is computationally efficient and spatially distributed, and benefits from temporal information by using a simple strategy that propagates two statistics with different learning rates. Simulations illustrate both the detection and localization performance of the proposed method.

II. PROBLEM FORMULATION

We consider a sequence of multi-or hyperspectral images Y t ∈ R L×N with N pixels and L bands, for t = 1, . . . , T . The unsupervised binary CD problem can be formulated as follows: given the image sequence Y t , find c t ∈ R N , for t = 2, . . . , T . where c t,j = 1 if there occurred a change in pixel j between instants t -1 and t, and c t,j = 0 otherwise. More precisely, we consider the measurements to be given by:

Y t = Q t + E t , (1)
where Q t denotes contextually relevant image content, and E t represents additive noise and perturbations that are noninformative changes, which include, e.g., noise, illumination and other acquisition variations [START_REF] Borsoi | Spectral variability in hyperspectral data unmixing: A comprehensive review[END_REF], [START_REF] Borsoi | Super-resolution for hyperspectral and multispectral image fusion accounting for seasonal spectral variability[END_REF]. This allows us to define the CD problem equivalently as:

c t,j = 1 ⇐⇒ q t,j = q t-1,j , t = 2, . . . , T, (2)
where q t,j is the j th column of Q t , for t = 1, . . . , T .

In this work, we make some simplifying assumptions that facilitate the definition of the change detection problem:

A1: Changes occur concurrently in sets of pixels that share semantically similar information (such as buildings, roads, crop fields, etc.). Mathematically, we assume that every pixel i for which c t,i = 1 is part of one among a set of compact, connected, spectrally homogeneous groups of pixels C (k), k = 1, . . . , N C for which c t,j = 1 for all pixels j in the same group as i.

A2: Only a single change occurs in a given streaming signal y t,j of length T (1 ≤ t ≤ T ).

A1 is reasonable as changes usually occur in compact regions of the image composed of several pixels. A2 allows us to explore temporal information contained in longer sequences [START_REF] Liu | A review of change detection in multitemporal hyperspectral images: Current techniques, applications, and challenges[END_REF], differently from multi-and hyperspectral image CD algorithms that focus in the bitemporal problem (i.e., when T = 2). Moreover, very long sequences where multiple changes occur can be divided into subsequences of smaller length T in which A2 is satisfied, as long as multiple abrupt changes do not happen in close succession.

A1 and A2 motivate a reformulation of the CD problem as the detection of abrupt variations in the data sequence. This simplified model is more accurate as the ratio of the sampling period of the discrete image sequence to the time required for the occurrence of a significant change increases. We use this simplified model for ease of exposition.

The change point detection (CPD) problem consists of finding t c ∈ N such that q t,j = q j , t < t c ,

q t,j = q j + δ j , t ≥ t c . (3)
where vector q j ∈ R L is the j th column of matrix Q, which captures the static underlying image content, and δ j ∈ R L is the j th column of matrix ∆, which captures the changes. This formulation is adequate for different applications such as surveillance or other industry applications [START_REF] Banerjee | Hyperspectral video for illumination-invariant tracking[END_REF], [START_REF] Stefansson | Hyperspectral time series analysis: hyperspectral image data streams interpreted by modeling known and unknown variations[END_REF], which can benefit from computationally efficient online algorithms.

A. Superpixel-based graph construction

An important aspect of this problem is how groups of pixels C (k) should be defined. Fundamentally, C (k) should contain pixels that share semantically similar information (e.g., belonging to the same objects or regions that may experience changes concurrently). We focus at linking pixels that are spatially adjacent and spectrally similar, attempting to respect sharp image borders that might indicate different objects or image regions. Recently, different works proposed to divide multiband images into disjoint, compact regions by means of image (over)-segmentation methods such as superpixels, ultrametric contour maps, or binary partition trees [START_REF] Achanta | SLIC superpixels compared to state-of-the-art superpixel methods[END_REF]- [START_REF] Veganzones | Hyperspectral image segmentation using a new spectral unmixing-based binary partition tree representation[END_REF]. In particular, superpixels are able to group spectrally similar pixels in compact spatial neighborhoods of average size S 2 with excellent preservation of image borders [START_REF] Achanta | SLIC superpixels compared to state-of-the-art superpixel methods[END_REF]. The superpixel decomposition has been sucesfully applied to construct multiscale transformations used in spectral unmixing applications [START_REF] Borsoi | A fast multiscale spatial regularization for sparse hyperspectral unmixing[END_REF]- [START_REF] Borsoi | A blind multiscale spatial regularization framework for kernel-based spectral unmixing[END_REF]. In this paper, we consider the superpixel decomposition in order to divide the pixels n = 1, . . . , N into sets C (k) that group the pixels belonging to each superpixel of the image.

III. PROPOSED STRATEGY

Traditional CD techniques operate individually on each pixel of an image. In this work we propose a framework that is able to explore, in a theoretically-grounded manner, the information that changes occur in groups of semantically similar pixels in possibly long image sequences, while at the same time benefiting from temporal information to operate online and computationally efficiently. We base our work in [START_REF] Ferrari | Distributed change detection in streaming graph signals[END_REF], which we extend to consider multidimensional data, and use a graph constructed from the superpixel decomposition.

We assume that the images are defined over an undirected graph G = {V, E, W } with N vertices V = {1, . . . , N }, M edges (i, j) ∈ E ⊆ V × V and a (symmetric) adjacency matrix W ∈ R N ×N whose elements are nonnegative and specify the similarity between the different pixels in the image. We also denote by L ∈ R N ×N the normalized graph Laplacian of G.

According to the discussion in Section II-A, G is constructed such that (i, j) ∈ E if and only if pixels i and j belong to the same superpixel, in which case [W ] (i,j) = 1.

This definition makes C (k) equivalent to clusters in G, and the CD problem becomes closely related to testing for changes in clusters of pixels in the graph. A change will be characterized by a δ j = 0 for j in a well-connected cluster of vertices of G. A detection strategy that is closely related to this problem is the graph Fourier scan statistic (GFSS) [START_REF] Sharpnack | Detecting anomalous activity on networks with the graph Fourier scan statistic[END_REF]. The GFSS addresses a counterpart of this problem for scalar (onedimensional) graph signals x ∈ R N measured over G by testing whether the average value of x is constant on all vertices, or if there is a cluster C of vertices with a different average value [START_REF] Sharpnack | Detecting anomalous activity on networks with the graph Fourier scan statistic[END_REF]. The GFSS considers the following test statistic:

r centrGFSS = g x 2 , (4)
where g(x) is the graph-filtered version of x, where the graph filter g is defined as

g x = N n=2 h(µ n ) u n x u n , (5)
with u n ∈ R N being the n th eigenvector of L, µ n the eigenvalue associated with u n , and

h(µ) = min 1, γ/µ , µ > 0 . (6)
Parameter γ > 0 in (6) controls the behavior of the graph filter that is applied to the vertex-wise test statistic, which attenuates high-frequency variations in x while enhancing signals that are consistent across well connected clusters of G (i.e., a low-pass graph filter) [START_REF] Shuman | The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains[END_REF]. Although the graph filter in (5) requires the eigendecomposition of L, which may be computationally costly, several strategies have been proposed to address this issue. Those typically involve approximating g by a finite impulse response [START_REF] Segarra | Optimal graph-filter design and applications to distributed linear network operators[END_REF], or by autoregressive moving average [START_REF] Isufi | Autoregressive moving average graph filtering[END_REF] graph filters that can be implemented efficiently and in a distributed manner (even for large N ). The GFSS, however, cannot be directly applied to the present problem since the underlying hypotheses are not satisfied (i.e., the graph signal is not constant over all vertices under the null hypothesis, and the nodal observations are vector valued pixels). If the underlying trend Q were known, it could be subtracted from Y t , and each band of the resulting signals could be tested using the GFSS. Since Q is not known, we adopt the adaptive strategy first considered in [START_REF] Ferrari | Distributed change detection in streaming graph signals[END_REF], [START_REF] Keriven | NEWMA: a new method for scalable model-free online change-point detection[END_REF], where two adaptive filters with different time constants λ and Λ, 0 < λ < Λ < 1, are used according to:

V t = (1 -λ)V t-1 + λY t , V t = (1 -Λ)V t-1 + ΛY t . (7)
Both V t and V t are weighted averages of Y t , which can mitigate the influence of additive noise E t in (1). Moreover, V t captures the short-term variations in Y t , from which the long-term average V t can be subtracted [START_REF] Keriven | NEWMA: a new method for scalable model-free online change-point detection[END_REF] as:

D t = V t -V t , (8)
where

D t = [d (1)
t , . . . , d

t ]. Assuming model ( 3), we have that D t is asymptotically zero-mean in the absence of a change-point (i.e., t < t c ). On the other hand, after a change-point occurs (i.e., t ≥ t c ), D t will be non-zero-mean at the clusters affected by the change in ∆, since V t will adapt more quickly than V t to represent the changes in Y t due to the difference between the time constants of both filters in (7) (i.e., a change can be detected by verifying how the short-term average differs from the long-term trend).

The graph signal D t is, however, vector valued (at each pixel), what still hinders the use of the GFSS test statistic as defined in [START_REF] Jianya | A review of multi-temporal remote sensing data change detection algorithms[END_REF]. Moreover, as noted in [START_REF] Ferrari | Distributed change detection in streaming graph signals[END_REF], the global test statistic r centrGFSS does not allow one to know in which cluster of G the change occurred. An alternative strategy proposed in [START_REF] Ferrari | Distributed change detection in streaming graph signals[END_REF] is to consider the coherent sum of the graph-filtered difference signal at the neighborhood of each vertex of G, which localizes the changes in the graph. Since we are testing for changes affecting all bands of pixels in clusters C (k) (i.e., we do not localize the changes at the band level), we propose a distributed test statistic for each vertex of G by summing the squared norms of the coherent sums of g(d

( )
t ) over all bands = 1, . . . , L. This leads to the following test statistic for each pixel:

r daGFSS (n) = L =1 m∈N (n) g d ( ) t 2 , (9)
where N (n) denotes the set of neighbors of vertex n in the graph G. Assuming that vec(E t ) ∼ N (0, σ 2 I), where vec(•) is the vectorization operator, it can be shown using [29, Proposition 2] that d

( ) t is asymptotically distributed as d ( ) t a ∼ N (0, ησ 2 I), with η = λ 2 -λ + Λ 2 -Λ - 2λΛ λ + Λ -λΛ . ( 10
)
This means that the graph-filtered d

( ) t is distributed as g(d ( ) t ) a ∼ N (0, R), where R = ησ 2 N n=2 h(µ n ) 2 u n u n , (11)
for = 1, . . . , L.

Similarly, m∈N (n) g(d

( ) t ) a ∼ N (0, σ 2 R (n))
, where

σ 2 R (n) = m,p∈N (n) ησ 2 N k=2 h(µ k ) 2 u k u k m,p , (12)
and [ • ] m,p denotes the (m, p) th element of a matrix. Note that the summands in (9) are statistically independent for different values of , which means that under the null hypothesis, r daGFSS (n) follows a Chi-squared distribution with L degrees of freedom. Then, a p-value of p n can be obtained by setting the test threshold ξ n according to

ξ n = σ R (n)Γ -1 inc L 2 ; 1 -(1 -p n )Γ L 2 , (13)
where Γ(•) is the common gamma function and Γ inc (• ; •) is the upper incomplete gamma function. Determining whether t is a change point of Y t is now equivalent to testing if a change occurred on at least one of the vertices, n = 1, . . . , N , which consists of a multiple testing problem. In order to limit the probability of type I errors in this context, we employ a False Discovery Rate (FDR) controlling procedure. Since the test statistics r daGFSS (n) are statistically dependent, we employ the Bonferroni correction as in [START_REF] Ferrari | Distributed change detection in streaming graph signals[END_REF], since it can be applied efficiently and distributedly in this case. Thus, for an FDR of α we set p n = α/N .

Algorithm 1: Remote sensing Graph CD Algorithm

Input : Images {Y t}, parameters 0 < λ < Λ < 1, γ > 0.

1 Compute the graph G based on the first observed HI using the procedure of Section II-A ;

2 Initialize V 1 = Y 1 and V 1 = Y 1 ; 3 for t = 2, 3, . . . do 4
Compute V t, V t and Dt using eqs. ( 7) and ( 8) ;

5

Compute the test statistic r daGFSS (n) using eq. ( 9);

IV. EXPERIMENTAL RESULTS

In this section, we evaluate the performance of the proposed method using two examples. We compare the proposed method with the traditional change vector analysis (CVA), with MAD [START_REF] Nielsen | Multivariate alteration detection (MAD) and MAF postprocessing in multispectral, bitemporal image data: New approaches to change detection studies[END_REF] and with the IRMAD [START_REF] Nielsen | The regularized iteratively reweighted MAD method for change detection in multi-and hyperspectral data[END_REF] in the first example, and also with the DSFA [START_REF] Du | Unsupervised deep slow feature analysis for change detection in multi-temporal remote sensing images[END_REF] in the second one. The first example illustrates the change-point detection performance for an abrupt change, while the second considers a more general scenario where varied amounts of changes are present in the scene. The superpixel decomposition is computed using the SLIC algorithm [START_REF] Achanta | SLIC superpixels compared to state-of-the-art superpixel methods[END_REF]. Snapshots of the sequences used in both examples are shown in Fig. 1.

Example 1: In this example, we consider a video sequence with T = 70 frames of a concrete floor, which is static at the beginning up until a water bottle is emptied in front of the camera. This generates an abrupt change point at t c = 16, after which frequent and intermittent changes are observed until the end of the sequence. The images contained L = 9 bands and were resized to N = 100 pixels to accelerate processing. Additive white Gaussian noise with a signal to noise ratio of 10 dB was added to each frame. To compare the performance of the algorithms, we ran a Monte Carlo simulation with 1000 realizations, and evaluated the mean detection delay (E{ t c -t c }, where t c is the detected change point) as a function of the probability of false alarm (evaluated as the probability of having a detection when t < t c ). The parameters of the proposed algorithm were set as S = 6, γ = 0.1, λ = 0.01, and Λ = 0.8. The results are presented in Fig. 2-(a), and show that the proposed method is able to provide a detection delay smaller than the other approaches.

Example 2: In this experiment, we consider a real multispectral video sequence of an outdoor scene with intermittent object motion (e.g., passing cars and pedestrians) as well as more subtle, non-informative changes (e.g., tree leaves moving). The video contained 50 × 50 pixels, L = 7 spectral bands and T = 1000 frames. Ground truth segmentation results of the moving targets of interest (e.g., cars and pedestrians) were available for this sequence [START_REF] Benezeth | Background subtraction with multispectral video sequences[END_REF], which we use to compute the true detections c t,n and evaluate the detection and localization performance of the algorithms. Denoting the estimated changes by c t,n , we compute the probability of detection and false alarm by averaging them across all pixels and across all frames:

Pd = ∀t,n : ct,n=1 c t,n ∀t,n c t,n , (14) Pfa
The parameters of the proposed algorithm were set as S = 5, γ = 0.1, λ = 0.15, and Λ = 0.5.

The results are shown in Fig. 2-(b) and 3. It can be seen that even though this sequence is challenging since the changes happen during short and intermittent periods, the proposed algorithm is able to obtain a better performance if the probability of false alarm is not very small.

V. CONCLUSIONS

In this paper, we proposed an online graph-based change detection algorithm for multitemporal multiband images. Acting on neighboring spectra as adjacent vertices in a graph, this algorithm focuses on anomalies concurrently activating groups of vertices corresponding to compact, well-connected and spectrally homogeneous image regions. It fully benefits from recent advances in graph signal processing and learns a graph to characterize the support of the data directly from the input image using superpixel decomposition methods. This allows us to exploit both spatial and spectral information. The method is unsupervised as it requires minimal knowledge about the nature of the changes, and is computationally efficient and scalable as it is spatially distributed. Preliminary experiments indicate superior performance of the proposed algorithm both in the detection and localization of changes.

  and J. Bermudez was funded in part by CNPq under grants 304250/2017-1, 409044/2018-0, 141271/2017-5 and 204991/2018-8. The work of C. Richard was funded in part by the ANR under grant ANR-19-CE48-0002, and by the 3IA Côte d'Azur Senior Chair program.

Figure 1 .

 1 Figure 1. Snapshot of the multispectral video sequences used in example 1 (left) and in example 2 (right).

6 for n = 1 ,

 61 . . . , N do 7 if r daGFSS (n) > ξn then 8 Flag t as a change point for vertex n and set ct,n ← 1 ;

Figure 2 .Figure

 2 Figure 2. Experimental results. (a) Mean detection delay as a function of the probability of false alarm (Example 1). (b) Probability of detection as a function of the probability of false alarm (Example 2).