Student Member, IEEE Jin Danqi

Senior Member, IEEE Jie Chen

Senior Member, IEEE Jingdong Chen

Online Proximal Learning Over Jointly Sparse Multitask Networks With ∞,1 Regularization

Keywords: Distributed optimization, diffusion strategy, joint sparsity, proximal algorithm, stochastic performance, (reweighted) ∞, 1-norm

Modeling relations between local optimum parameter vectors to estimate in multitask networks has attracted much attention over the last years. This work considers a distributed optimization problem with jointly sparse structure among nodes, that is, the local solutions have the same sparse support set. Several mixed norm have been proposed to address the jointly sparse structure in the literature. Among several candidates, the (reweighted) ∞,1-norm is element-wise separable, it is more convenient to evaluate their approximate proximal operators. Thus by introducing a (reweighted) ∞,1-norm penalty term at each node, and using a proximal gradient method to minimize the regularized cost, we devise a proximal multitask diffusion LMS algorithm which can promote joint-sparsity. Analyses are provided to characterize the algorithm behavior in the mean and mean-square sense. Simulation results are presented to show its effectiveness, as well as the accuracy of the theoretical findings.

I. INTRODUCTION

Because of their superior performance and wider stability range [START_REF] Tu | Diffusion strategies outperform consensus strategies for distributed estimation over adaptive netowrks[END_REF], diffusion strategies have been widely used in multi-agent networks to address estimation problems in a distributed and online manner. Several diffusion strategies have been introduced, and their performance analyzed in various situations, such as the diffusion LMS [START_REF] Lopes | Diffusion least-mean squares over adaptive networks: Formulation and performance analysis[END_REF], RLS [START_REF] Cattivelli | Diffusion recursive least-squares 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], as well as several of their variants [START_REF] Liu | Diffusion sparse least-mean squares over networks[END_REF], [START_REF] Lorenzo | Sparse distributed learning based on diffusion adaptation[END_REF].

By referring to estimating an optimal parameter vector at a node as a task, and according to the relations between the optimal parameter vectors over the entire network, diffusion networks are further divided into single-task and multitask networks. In single-task networks, all nodes estimate the same parameter vector. Typical works related to single-task networks include [START_REF] Sayed | Diffusion adaptation over networks[END_REF]- [START_REF] Mao | Walkman: A communication-efficient random-walk algorithm for decentralized optimization[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 exploiting the similarities between tasks. These similarities can be promoted with appropriate regularization terms. Squared 2 -norm regularization is used in [START_REF] Chen | Multitask diffusion adaptation over networks[END_REF], and 1 -norm regularization is considered in [START_REF] Nassif | Multitask diffusion LMS with sparsity-based regularization[END_REF], [START_REF] Nassif | Proximal multitask learning over networks with sparsity-inducing coregularization[END_REF]. For the latter, a subgradient and a proximal algorithm are introduced in [START_REF] Nassif | Multitask diffusion LMS with sparsity-based regularization[END_REF] and [START_REF] Nassif | Proximal multitask learning over networks with sparsity-inducing coregularization[END_REF], respectively. In [START_REF] Chen | Multitask diffusion adaptation over networks with common latent representations[END_REF] and [START_REF] Nassif | Diffusion LMS for multitask problems with local linear equality constraints[END_REF], the authors derive solutions for other classes of multitask problems where the relations between the nodes are defined by common latent representations or local linear equality constraints, respectively. In [START_REF] Chen | Diffusion LMS over multitask networks[END_REF], the authors solve a multitask problem by estimating the combination matrix. In [START_REF] Nassif | Multitask diffusion adaptation over asynchronous networks[END_REF], the authors address multitask problems over asynchronous networks and carry out a detailed theoretical analysis. In [START_REF] Piggott | Stability of adaptive network algorithms in multitask environments[END_REF], the performance of multitask diffusion networks is analyzed for correlated noise and regressors. In [START_REF] Jin | Affine combination of diffusion stretegies over networks[END_REF], the authors propose a combination framework that aggregates several diffusion strategies. All the algorithms cited above are based on the diffusion LMS. In [START_REF] Gogineni | Diffusion affine projection algorithm for multitask networks[END_REF], the authors extend the diffusion APA to multitask framework in order to improve the robustness against correlated regressors. In [START_REF] Gogineni | Improving the performance of multitask diffusion APA via controlled inter-cluster cooperation[END_REF], the authors improve the performance of the multitask diffusion APA via controlled inter-cluster cooperation. In [START_REF] Gogineni | Partial diffusion affine projection algorithm over clustered multitask networks[END_REF], the authors propose a clustered multitask partial diffusion APA that transmits only a subset of the entries of the intermediate estimates, to provide a trade-off between the estimation performance and communication cost.

Multitask learning considerably enriches the possibilities of diffusion networks. Beyond the few examples listed above, there are also applications where the optimal parameter vectors have a jointly sparse structure, namely, local solutions have the same sparse support. Applications include, for instance, distributed spectrum sensing and channel identification in underwater and wireless communication networks with multiple sensors [START_REF] Li | Distributed jointly sparse multitask learning over networks[END_REF]. In spectrum sensing of sparse wide-band spectra in distributed wireless sensor networks, the parameter vectors at different nodes share the same intrinsic sparse structure. However, due to different channel fading effects, shadowing effects and transmission losses, the values at the nonzero entries of these parameter vectors are different [START_REF] Li | Distributed jointly sparse multitask learning over networks[END_REF]- [START_REF] Duarte | Distributed compressed sensing of jointly sparse signals[END_REF]. For underwater and wireless communication networks, it has been shown that real-world underwater channels [START_REF] Kocic | Sparse equalization for realtime digital underwater acoustic communications[END_REF] and wireless communication channels [START_REF] Masood | Efficient coordinated recovery of sparse channels in massive MIMO[END_REF] are inherently sparse with large delay spread. In addition,

The work of Jie Chen was supported in part by NSFC grant 61671382. The work of Jingdong Chen was supported in part by NSFC grant 61425005. The work of C. Richard work was supported through the UCA JEDI Investments in the Future project with the reference number ANR-15-IDEX-0001. D. Jin, J. Chen and J. Chen are with Centre of Intelligent Acoustics and Immersive Communications at School of Marine Science and Technology, Northwestern Polytechnical University, Xi'an, China (danqijin@mail.nwpu.edu.cn, dr.jie.chen@ieee.org, jingdongchen@ieee.org). C. Richard is with Université Côte d'Azur, OCA, CNRS, France (cedric.richard@unice.fr).

the channel supports for neighboring antennas or nodes are approximately the same [START_REF] Masood | Efficient coordinated recovery of sparse channels in massive MIMO[END_REF]. Indeed, times of arrival for closely spaced nodes and antennas are quite close, though the tap weights are actually different [START_REF] Li | Distributed jointly sparse multitask learning over networks[END_REF]. For more detailed information about this application, see Section VI-C. Taking jointly sparse structure into consideration can greatly improve the network performance.

Several works have been proposed to address problems with jointly sparse structure over diffusion networks. In [START_REF] Gu | Learning distributed jointly sparse systems by collaborative LMS[END_REF], the authors consider the mixed 2,0 -norm. In [START_REF] Li | Distributed jointly sparse multitask learning over networks[END_REF], the authors devise an algorithm with regularizers such as 2,0 , 2,1 and reweighted 2,1 regularizers. They also conduct theoretical analyses of the algorithm behavior in the mean and mean-square sense. However, studies in [START_REF] Li | Distributed jointly sparse multitask learning over networks[END_REF], [START_REF] Gu | Learning distributed jointly sparse systems by collaborative LMS[END_REF] are based on gradient-descent schemes, and use subgradient for non-differential regularization terms. For subgradient approaches, the subdifferential at a point may not be a singleton set, that is, it may be empty or consist of several elements. As a result, one may get stuck or have to choose one, respectively. Furthermore, even if the subdifferential is a singleton at each step, it might be highly discontinuous, so small deviations might lead to a singular behavior of the algorithm over iterations [START_REF] Banert | A general double-proximal gradient algorithm for d.c. programming[END_REF]. In this paper, we propose to use proximal operators to address the jointly sparse estimation problem and then avoid the weaknesses of subgradient-based approaches mentioned above. Proximal algorithms result in subproblems that often admit closed-form solutions or that can be solved efficiently with simple specialized methods [START_REF] Parikh | Proximal algorithms[END_REF]. In addition, they guarantee better convergence rates and stability than subgradient approaches [START_REF] Banert | A general double-proximal gradient algorithm for d.c. programming[END_REF], [START_REF] Bertsekas | Incremental proximal methods for large scale convex optimization[END_REF]. For single-task learning problems, proximal algorithms were first introduced in [START_REF] Wee | A proximal splitting approach to regularized distributed adaptive estimation in diffusion networks[END_REF] to estimate sparse parameter vectors. They were used in [START_REF] Vlaski | Proximal diffusion for stochastic costs with non-differentiable regularizers[END_REF] and [START_REF] Vlaski | Diffusion stochastic optimization with non-smooth regularizers[END_REF] to optimize general stochastic costs with non-differential regularizers, and non-smooth regularizers, respectively. For multitask learning problems, related works only include [START_REF] Nassif | Proximal multitask learning over networks with sparsity-inducing coregularization[END_REF]. In that paper, the authors consider the situation where the network is divided into several clusters. All nodes in a cluster are interested in estimating the same parameter vector, while nodes in adjacent clusters estimate parameter vectors that have a large number of similar entries. As a result the authors derive a closed-form proximal solution for an 1 -norm regularizer that promotes similarities among clusters. A theoretical analysis of the steady-state behavior is provided. The algorithm in [START_REF] Nassif | Proximal multitask learning over networks with sparsity-inducing coregularization[END_REF] is based on an extra exchange step of observations between neighboring nodes within a same cluster as well as a fusion step that averages local estimates. The proximal algorithm derived in the current work does not include these two steps.

In this paper, we introduce a proximal diffusion LMS strategy for multitask networks with jointly sparse structure. The main contributions of this work are summarized as follows:

• We derive approximate closed-form expressions for the ∞,1 -norm and reweighted ∞,1 -norm proximal operators.

• We derive a proximal multitask diffusion LMS algorithm to solve problems with jointly sparse structure.

• We conduct a theoretical analysis of the algorithm performance, including a condition for stability and a study of its transient behavior in the mean sense. Notation. Normal font x, boldface small letters x and capital letters X denote scalars, column vectors and matrices, respectively. Symbol

II. PROBLEM FORMULATION

Consider a connected network consisting of N nodes. Each node has access to streaming data {d k,n , x k,n }, where x k,n is the L × 1 real-valued regression vector at node k and time instant n, and d k,n denotes the observed real-valued signal. We assume that the data at each agent k and time instant n are related via the linear model:

d k,n = x k,n w k + z k,n , (1)
where w k ∈ IR L×1 is the unknown system vector to estimate, and z k,n is a zero-mean additive noise. We assume that z k,n is independent of any other signal. Further, we assume that vectors {w k } N k=1 are jointly sparse, namely, not only each w k is a sparse vector but, in addition, they all have the same support. The support of w k is defined as [START_REF] Huang | The benefit of group sparsity[END_REF]:

supp(w k ) {j : [w k ] j = 0}, (2)
where [w k ] j is the j-th entry of w k . By definition (2), jointly sparse structure means that:

supp(w 1 ) = • • • = supp(w k ) = • • • supp(w N ). ( 3
)
Since it is not trivial to solve problems with jointly sparse structure directly, several mixed-norms have been introduced in the literature to alleviate this problem. They include the mixed 2,1 -norm [START_REF] Meier | The group lasso for logistic regression[END_REF], ∞,1 -norm [START_REF] Chen | Recursive 1,∞ group lasso[END_REF] and their reweighted versions [START_REF] Chen | Regularized least-mean-square algorithms[END_REF]. We collect w over the entire network into an L × N matrix, and we replace the k-th column w k by the optimization variable w k . This leads us to consider:

W glob k w 1 • • • w k-1 w k w k+1 • • • w N . (4)
Since matrices defined in (4) differ from node to node, we use the subscript k in notation W glob k to distinguish them. Evaluating the mixed p,1 -norm of matrix W glob k with p = 2 or ∞, so as to promote the jointly sparse structure, results in the following two steps:

Step 1: Evaluate the p -norm of each row of W glob k , and stack the results into an intermediate vector of dimension L × 1;

Step 2: Evaluate the 1 -norm of the obtained intermediate vector to promote sparsity. Though 2,1 -norm can be more efficient in some cases [START_REF] Negahban | Simultaneous support recovery in high dimensions: Benefits and perils of block 1 / ∞-regularization[END_REF], we shall consider the ∞,1 -norm and its reweighted form to promote the joint-sparsity. It is shown in [START_REF] Esser | A convex model for nonnegative matrix factorization and dimensionality reduction on physical space[END_REF] that the ∞,1 relaxation is exact in the case of normalized nonrepeating data. [START_REF] Chen | Recursive 1,∞ group lasso[END_REF] introduced a recursive adaptive group lasso algorithm for real-time penalized least squares prediction by using the 1,∞ regularization (The anthors defined ∞,1 -norm as 1,∞ -norm in [START_REF] Chen | Recursive 1,∞ group lasso[END_REF] for a vector with group sparsity structure). Each update minimizes a convex but nondifferentiable function optimization problem. The authors developed an online homotopy method to reduce the computational complexity. [START_REF] Esser | A convex model for nonnegative matrix factorization and dimensionality reduction on physical space[END_REF] proposed a collaborative convex framework for factoring a data matrix into a nonnegative matrices product. The authors used ∞,1 regularization to select the dictionary from the data and shown that this leads to an exact convex relaxation of 0 regularization in the case of distinct noise-free data. [START_REF] Eksioglu | Group sparse RLS algorithms[END_REF] proposed a class of group sparse RLS algorithms by using different penalty term to promote the group sparsity, where the ∞,1 -norm is used as one of the penalty term and the subgradient of ∞,1 -norm is evaluated. Both ∞,1 -norm and its reweighted form are element-wise separable, which facilitates the derivation of the proximal operator in the current work. As we focus on distributed processing in the current work, only local information exchange is authorized. Thus, we restrict W glob k to the local quantity:

W k w k , w with ∈ N - k ∈ IR L×|N k | . ( 5
)
We consider that the columns of W k are sorted in increasing order according to the values of k and . In the sequel, we shall show how to evaluate the (reweighted) ∞,1 -norm of W k .

III. PROXIMAL MULTITASK DIFFUSION LMS

Before proceeding, we rewrite W k in another way in order to facilitate the presentation:

W k = w k,1 • • • w k,m • • • w k,L , (6)
where wk,m is the m-th row of matrix W k and defined by

wk,m [w k ] m , [w ] m with ∈ N - k ∈ IR 1×|N k | (7)
with [w k ] m being the m-th entry of w k . Definition [START_REF] Sayed | Diffusion adaptation over networks[END_REF] means that wk,m is a row vector. Note that (5) and ( 6) are actually the same matrix, but represented in different ways. Also, wk,m ∞ is defined as

wk,m ∞ max{|[w k ] m |, |[w ] m | with ∈ N - k } (8)
which will be used in g 1 (w k ) and g 2 (w k ) of ( 11) and [START_REF] Chen | Multitask diffusion adaptation over networks[END_REF], respectively.

To determine the unknown vectors w k with jointly sparse structure, we consider the regularized cost at node k:

J k (w k ) = J k (w k ) + λ k g i (w k ), (9)
where J k (w k ) is the mean-square error (MSE) defined as:

J k (w k ) 1 2 E |d k,n -x k,n w k | 2 . ( 10
)
The nonnegative parameter λ k is used to control the regularization strength, g 1 (w k ) and g 2 (w k ) evaluate the ∞,1 -norm and reweighted ∞,1 -norm of W k , respectively, with:

g 1 (w k ) L m=1 wk,m ∞ , (11)
g 2 (w k ) L m=1 log 1 + wk,m ∞ ε , (12)
where ε > 0 is a parameter set by the user. Interpretations about the ∞,1 -norm g 1 (w k ) and reweighted ∞,1 -norm g 2 (w k ) are provided in Appendix A. At each node k, we then consider the convex optimization problem [START_REF] Tropp | Algorithms for simultaneous sparse approximation[END_REF]:

w † k = arg min w k J k (w k ). ( 13
)
Within the context of online learning, such optimization problem is usually solved via subgradient-based methods. In this paper, we propose to devise a proximal algorithm since these algorithms are usually more stable and have a better convergence rate than subgradient iterations [START_REF] Banert | A general double-proximal gradient algorithm for d.c. programming[END_REF], [START_REF] Bertsekas | Incremental proximal methods for large scale convex optimization[END_REF], [START_REF] Wee | A proximal splitting approach to regularized distributed adaptive estimation in diffusion networks[END_REF]. Proximal gradient iteration consists of [START_REF] Parikh | Proximal algorithms[END_REF]:

w k,n+1 = prox µ k λ k ,gi(w k ) w k,n -µ k ∇J k (w k,n ) , (14)
where µ k is a positive small step-size, and:

prox λ,gi(w k ) (v) arg min w k g i (w k ) + 1 2λ w k -v 2 2 . ( 15
)
is the proximal operator. By introducing the intermediate quantity ψ k,n+1 , we further decompose [START_REF] Nassif | Proximal multitask learning over networks with sparsity-inducing coregularization[END_REF] into two steps:

ψ k,n+1 = w k,n -µ k ∇J k (w k,n ), (16)
w k,n+1 = prox µ k λ k ,gi(w k ) ψ k,n+1 . (17)
Equation ( 16) is the local update step, and ( 17) is the proximal step. Calculating the gradient of J k (w k ) at w k,n and approximating the unknown statistical quantities with instantaneous quantities, ( 16) becomes:

ψ k,n+1 = w k,n + µ k x k,n d k,n -x k,n w k,n . (18)
Combining ( 17) and ( 18) yields the proximal multitask diffusion LMS algorithm for networks with jointly sparse structure reported in Algorithm 1. This algorithm differs from the standard diffusion procedure devised in [START_REF] Sayed | Adaptive networks[END_REF], [START_REF] Sayed | Adaptation, Learning, and Optimization over Networks[END_REF] in the sense that it does not combine, in an explicit manner, the intermediate estimates ψ ,n+1 in the neighborhood of each node k. In our algorithm, the information on the support of the local estimates is shared by neighboring nodes via g i (w k ).

Algorithm 1 Proximal multitask diffusion LMS Initialize w k,0 for all k = 1, 2, • • • , N, and repeat:

ψ k,n+1 = w k,n + µ k x k,n d k,n -x k,n w k,n w k,n+1 = prox µ k λ k ,gi(w k ) ψ k,n+1 (19)
IV. PROXIMAL OPERATORS EVALUATION To apply the algorithm, we need to derive closed-form expressions for the proximal operators. Equation [START_REF] Chen | Diffusion LMS over multitask networks[END_REF] becomes:

w k,n+1 = prox µ k λ k ,gi(w k ) (ψ k,n+1 ) = argmin w k g i (w k )+ 1 2µ k λ k w k -ψ k,n+1 2 2 . (20)
As g i (w k ) is separable over its all entries, the proximal operator can be evaluated in an element-wise manner as [START_REF] Parikh | Proximal algorithms[END_REF]:

[prox µ k λ k ,gi(w k ) (ψ k,n+1 )] m = prox µ k λ k ,gi,m([w k ] m ) ([ψ k,n+1 ] m ) (21)
with:

g 1,m ([w k ] m ) wk,m ∞ (22)
g 2,m ([w k ] m ) log 1 + wk,m ∞ ε , (23)
where wk,m is the m-th row of matrix W k in [START_REF] Liu | Diffusion sparse least-mean squares over networks[END_REF].

A. Approximate proximal operator of ∞,1 -norm Substituting g 1,m into (21), we obtain:

[w k,n+1 ] m = argmin [w k ] m max{|[w k ] m |, |[w ] m | with ∈ N - k } + 1 2µ k λ k [w k ] m -[ψ k,n+1 ] m 2 . ( 24
)
For the sake of simpler notations, we shall denote [w k,n+1 ] m by ŵ as long as there is no ambiguity. We shall also denote the maximal value of

|[w ] m | for ∈ N - k as [w o k ] m . According to the relation between |[w k ] m | and [w o k ] m
, we further split the problem into the following two cases:

• Case 1: |[w k ] m | < [w o k ] m .
In this case, [START_REF] Li | Distributed jointly sparse multitask learning over networks[END_REF] becomes:

ŵ = argmin [w k ] m |[w k ] m |<[w o k ] m [w o k ] m + 1 2µ k λ k [w k ] m -[ψ k,n+1 ] m 2 . ( 25
)
Since the first term on the right-hand side (RHS) of ( 25) does not depend on [w k ] m , we conclude that:

ŵ =      [ψ k,n+1 ] m , if |[ψ k,n+1 ] m | < [w o k ] m [w o k ] m , if [ψ k,n+1 ] m ≥ [w o k ] m -[w o k ] m , if [ψ k,n+1 ] m ≤ -[w o k ] m . (26)
• Case 2: 24) becomes:

|[w k ] m | ≥ [w o k ] m . Equation (
ŵ = argmin [w k ] m |[w k ] m |≥[w o k ] m |[w k ] m | + 1 2µ k λ k [w k ] m -[ψ k,n+1 ] m 2 . ( 27
)
We shall first discard the constraint

|[w k ] m | ≥ [w o k ]
m for reasons of simplicity, and denote by ŵo the solution of the unconstrained problem. As the cost function in ( 27) is convex on IR, this constraint will be taken into account in the course of the calculation. Consider first:

ŵo = argmin [w k ] m |[w k ] m | + 1 2µ k λ k [w k ] m -[ψ k,n+1 ] m 2 = prox µ k λ k ,s([w k ] m ) [ψ k,n+1 ] m (28)
where function s

([w k ] m ) is defined as: s([w k ] m ) |[w k ] m |. (29)
The optimality condition of ( 28) says that zero belongs to the subgradient set at the minimizer ŵo , that is, [START_REF] Parikh | Proximal algorithms[END_REF] 0

∈ ∂| ŵo | + 1 µ k λ k ŵo -[ψ k,n+1 ] m (30)
This means that:

[ψ k,n+1 ] m -ŵo ∈ µ k λ k ∂| ŵo | (31)
with

∂| ŵo | =      -1, if ŵo < 0 +1, if ŵo > 0 [-1, 1], if ŵo = 0 (32)
the subdifferential of the non-differentiable function | ŵo | [START_REF] Parikh | Proximal algorithms[END_REF]. The third case in [START_REF] Bertsekas | Incremental proximal methods for large scale convex optimization[END_REF] means that, at ŵo = 0, ∂| ŵo | can be any value within [-1, 1]. Condition (31) leads to:

ŵo =      [ψ k,n+1 ] m + µ k λ k , if [ψ k,n+1 ] m < -µ k λ k [ψ k,n+1 ] m -µ k λ k , if [ψ k,n+1 ] m > µ k λ k 0, otherwise. (33)
If [START_REF] Kocic | Sparse equalization for realtime digital underwater acoustic communications[END_REF] becomes unconstrained and we have:

[w o k ] m = 0, problem (
ŵ = ŵo (34)
Otherwise, since problem ( 27) is convex on

IR, considering constraint |[w k ] m | ≥ [w o k ] m > 0 with (33) yields: ŵ = (35)                [ψ k,n+1 ] m +µ k λ k , if [ψ k,n+1 ] m ≤ -[w o k ] m -µ k λ k -[w o k ] m , if -[w o k ] m -µ k λ k < [ψ k,n+1 ] m < 0 -[w o k ] m or [w o k ] m , if [ψ k,n+1 ] m = 0 [w o k ] m , if 0 < [ψ k,n+1 ] m < [w o k ] m +µ k λ k [ψ k,n+1 ] m -µ k λ k , if [ψ k,n+1 ] m ≥ [w o k ] m +µ k λ k
To evaluate the proximal operator [START_REF] Li | Distributed jointly sparse multitask learning over networks[END_REF], several issues have to be addressed.

1. One of the main issues is that we first need to know which of ( 26), [START_REF] Vlaski | Proximal diffusion for stochastic costs with non-differentiable regularizers[END_REF] or [START_REF] Vlaski | Diffusion stochastic optimization with non-smooth regularizers[END_REF] has to be applied as the proximal operator of [START_REF] Li | Distributed jointly sparse multitask learning over networks[END_REF]. We shall now consider the following two cases:

[w o k ] m = 0 and [w o k ] m > 0. • Case A: [w o k ] m = 0. Since condition |[w k ] m | < [w o k ] m
of Case 1 cannot hold, we only consider Case 2. The proximal operator is given by ŵo in [START_REF] Wee | A proximal splitting approach to regularized distributed adaptive estimation in diffusion networks[END_REF]. Interestingly, note that [START_REF] Wee | A proximal splitting approach to regularized distributed adaptive estimation in diffusion networks[END_REF] is the proximal operator of the 1 -norm regularizer for a scalar.

• Case B: [w o k ] m > 0.
Proximal operators ( 26) and ( 35) hold simultaneously. We shall choose the solution that minimizes the cost [START_REF] Li | Distributed jointly sparse multitask learning over networks[END_REF]. As shown in Appendix B, we finally arrive at the following expression:

ŵ = (36)                [ψ k,n+1 ] m + µ k λ k , if [ψ k,n+1 ] m ≤ -[w o k ] m -µ k λ k -[w o k ] m , if -[w o k ] m -µ k λ k < [ψ k,n+1 ] m ≤ -[w o k ] m [ψ k,n+1 ] m , if [ψ k,n+1 ] m < [w o k ] m [w o k ] m , if [w o k ] m ≤ [ψ k,n+1 ] m < [w o k ] m + µ k λ k [ψ k,n+1 ] m -µ k λ k , if [w o k ] m + µ k λ k ≤ [ψ k,n+1 ] m 2.
Another issue is that ŵ cannot be evaluated with ( 33) and ( 36) since [w o k ] m is unknown. To fix this problem, we follow a local one-step approximation strategy that has already proven its effectiveness in the literature [START_REF] Chen | Diffusion LMS over multitask networks[END_REF], and which consists of using ψ ,n+1 as an approximation of w . One of the benefits of this approximation is that ψ ,n+1 can be transmitted by node to node k if the latter is in its neighborhood. An approximation of

[w o k ] m is then given by max ∈N - k { [ψ ,n+1
] m }, which allows node k to evaluate its proximal operator. This approximation will be taken into consideration in the theoretical analysis in Section V-A.

3. Before addressing the last issue, we need to point out the prominent role of Case A compared to Case B: unlike the proximal operator [START_REF] Huang | The benefit of group sparsity[END_REF] in Case B, only the proximal operator [START_REF] Wee | A proximal splitting approach to regularized distributed adaptive estimation in diffusion networks[END_REF] in Case A has the capability to drive [w k ] m to zero and promote sparsity. Observe that condition [w o k ] m = 0 has to be satisfied to trigger Case A, otherwise Case B is considered. Within the context of online learning with stochastic gradient descent algorithms, due to the existence of gradient noise, the estimates of [w o k ] m for zero-valued entries are actually nonzero-valued but vary around zero. Condition [w o k ] m = 0 of Case A is thus seldom satisfied, and [ψ k,n+1 ] m is not driven to zero. To promote the sparsity of the estimates, since the true values of the non-zero entries are usually far away from zero, we introduce a small positive threshold value τ 1 instead of zero to make a distinguish between zero-valued and nonzero-valued entries. As a consequence, we arrive at conditions [w o k ] m ≤ τ 1 to trigger Case A and [w o k ] m > τ 1 to select Case B. We summarize our method in Algorithm 2.

Algorithm 2 Approximate proximal operator of ∞,1 -norm Initialization: Choose threshold value τ 1 > 0. Proximal operator: At each instant n ≥ 0, for each node k, utilize ψ k,n+1 to evaluate w k,n+1 in an elementwise manner:

1) Calculate [w o k ] m as the maximal value of [ψ ,n+1 ] m for all ∈ N - k ; 2) If [w o k ] m ≤ τ 1 , then calculate [w k,n+1 ] m as ŵo via (33); 3) If [w o k ] m > τ 1 , then calculate [w k,n+1
] m as ŵ via [START_REF] Huang | The benefit of group sparsity[END_REF].

B. Approximate proximal operator of reweighted ∞,1 -norm Substituting the expression of g 2,m into (21), we arrive at:

[w k,n+1 ] m = argmin [w k ] m log 1 + max |[w k ] m |, |[w ] m | : ∀ ∈ N - k ε + 1 2µ k λ k [w k ] m -ψ k,n+1 m 2 . ( 37
)
Considering the same definition for [w o k ] m as before, and denoting [w k,n+1 ] m by ŵ for the sake of conciseness, the problem can be split into two cases depending on the order relation between

[w o k ] m and |[w k ] m |: • Case 1: |[w k ] m | < [w o k ] m . Equation (37) becomes: ŵ = argmin [w k ] m |[w k ] m |<[w o k ] m log 1 + [w o k ] m ε + 1 2µ k λ k [w k ] m -ψ k,n+1 m 2 . ( 38
)
This problem has a closed-form solution given by:

ŵ =      [ψ k,n+1 ] m , if |[ψ k,n+1 ] m | < [w o k ] m [w o k ] m , if [ψ k,n+1 ] m ≥ [w o k ] m -[w o k ] m , if [ψ k,n+1 ] m ≤ -[w o k ] m . (39)
• Case 2:

|[w k ] m | ≥ [w o k ] m .
In this case, (37) becomes:

ŵ = argmin [w k ] m |[w k ] m |≥[w o k ] m log 1 + |[w k ] m | ε + 1 2µ k λ k [w k ] m -ψ k,n+1 m 2 . ( 40
)
By first discarding the constraint

|[w k ] m | ≥ [w o k ]
m as we did it before, we obtain the following unconstrained optimization problem:

ŵo = argmin [w k ] m log 1 + |[w k ] m | ε + 1 2µ k λ k [w k ] m -ψ k,n+1 m 2 . ( 41
)
Deriving the solution ŵo directly from the optimality condition is a tough problem. We shall now alleviate this issue with a Majorization-Minimization (MM) approach [START_REF] Hunter | Departments Of Biomathematics, and Human Genetics[END_REF], [START_REF] Candès | Enhancing sparsity by reweighted 1 minimization[END_REF]. Conceptually, MM algorithms work by iteratively minimizing a simple surrogate function majorizing a given objective function. By introducing the auxiliary variable [u k ] m , problem (41) can be written as:

ŵo = argmin [w k ] m ,[u k ] m |[w k ] m |≤[u k ] m log 1 + [u k ] m ε + 1 2µ k λ k [w k ] m -ψ k,n+1 m 2 . (42)
Define:

f ([u k ] m ) log 1 + [u k ] m ε . ( 43
)
It can be observed that function f ([u k ] m ) is concave and below its tangent. So we have:

f ([u k ] m ) ≤ f ([u (t) k ] m ) + f ([u (t) k ] m ) • [u k ] m -[u (t) k ] m (44)
with f (•) the first-order derivative of f (•), and where the superscript (t) denotes the iteration index. Using [START_REF] Hunter | Departments Of Biomathematics, and Human Genetics[END_REF], we can construct a surrogate function for [START_REF] Eksioglu | Group sparse RLS algorithms[END_REF] as follows:

h([u k ] m , [w k ] m ) log 1+ [u (t) k ] m ε + [u k ] m -[u (t) k ] m ε + [u (t) k ] m + 1 2µ k λ k [w k ] m -ψ k,n+1 m 2 (45)
Using the driving principle of MM algorithms, this yields the following optimization problem:

[w (t+1) k ] m ,[u (t+1) k ] m = argmin [w k ] m ,[u k ] m |[w k ] m |≤[u k ] m [u k ] m ε + [u (t) k ] m + 1 2µ k λ k [w k ] m -ψ k,n+1 m 2 (46)
which is equivalent to:

[w (t+1) k ] m = argmin [w k ] m |[w k ] m | ε + |[w (t) k ] m | + 1 2µ k λ k [w k ] m -ψ k,n+1 m 2 . ( 47
)
Note that, in the construction of ( 46), we have dropped several constant terms unrelated to [u k ] m and [w k ] m . Using the optimality condition of (47) at the minimizer [w

(t+1) k
] m , we have:

ψ k,n+1 m -[w (t+1) k ] m ∈ µ k λ k ε + |[w (t) k ] m | • ∂ [w (t+1) k ] m . (48)
Considering ( 32) with (48) yields:

[w (t+1) k ] m = (49)                [ψ k,n+1 ] m + µ k λ k ε+|[w (t) k ] m | , if [ψ k,n+1 ] m < -µ k λ k ε+|[w (t) k ] m | [ψ k,n+1 ] m - µ k λ k ε+|[w (t) k ] m | , if [ψ k,n+1 ] m > µ k λ k ε+|[w (t) k ] m | 0 otherwise.
Performing one iteration of ( 49) with [w

(t) k ] m = [w k,n
] m at each instant n can be sufficient for approximating ŵo with enough precision, as illustrated in the sequel. This value is then used to obtain an approximate solution ŵ of problem [START_REF] Negahban | Simultaneous support recovery in high dimensions: Benefits and perils of block 1 / ∞-regularization[END_REF] by taking the constraint [START_REF] Negahban | Simultaneous support recovery in high dimensions: Benefits and perils of block 1 / ∞-regularization[END_REF] becomes unconstrained and we have:

|[w k ] m | ≥ [w o k ] m into account. If [w o k ] m = 0, problem (
ŵ = ŵo ( 50
)
where ŵo is obtained by [START_REF] Arenas-Garcia | Mean-square performance of a convex combination of two adaptive filters[END_REF].

Otherwise, since problem ( 40) is convex on

IR, considering constraint |[w k ] m | ≥ [w o k ] m > 0 with (49) yields: ŵ = (51)                          [ψ k,n+1 ] m + µ k λ k ε + |[w (t) k ] m | , if [ψ k,n+1 ] m ≤ -[c k ] m -[w o k ] m , if -[c k ] m < [ψ k,n+1 ] m < 0 -[w o k ] m or [w o k ] m , if [ψ k,n+1 ] m = 0 [w o k ] m , if 0 < [ψ k,n+1 ] m < [c k ] m [ψ k,n+1 ] m - µ k λ k ε + |[w (t) k ] m | , if [ψ k,n+1 ] m ≥ [c k ] m
with [c k ] m a constant defined as:

[c k ] m [w o k ] m + µ k λ k ε+|[w (t) k ] m | . ( 52
)
As with the ∞,1 -norm, we have to fix several issues. 1. We first need to check which of ( 39), ( 49) and [START_REF] Silva | A transient analysis for the convex combination of two adaptive filters with transfer of coefficients[END_REF] has to be applied to calculate the proximal operator. We consider the following two cases:

[w o k ] m = 0 and [w o k ] m > 0. • Case A: [w o k ] m = 0.
We focus on Case 2 since condition for Case 1 cannot hold. In Case 2, observe that any

[w k ] m calculated as [w (t+1) k ] m in (49) satisfies |[w k ] m | ≥ [w o k ] m = 0.
Thus the proximal operator is given by [START_REF] Arenas-Garcia | Mean-square performance of a convex combination of two adaptive filters[END_REF].

• Case B: [w o k ] m > 0.
Case 1 and Case 2 can hold simultaneously. We must choose, among the resulting proximal operators, the one that minimizes cost [START_REF] Meier | The group lasso for logistic regression[END_REF]. The derivation is provided in Appendix C. We obtain the following solution: We summarize our method in Algorithm 3.

ŵ = (53)                          [ψ k,n+1 ] m + µ k λ k ε + |[w (t) k ] m | , if [ψ k,n+1 ] m ≤ -[c k ] m -[w o k ] m , if -[c k ] m < [ψ k,n+1 ] m ≤ -[w o k ] m [ψ k,n+1 ] m , if [ψ k,n+1 ] m < [w o k ] m [w o k ] m , if [w o k ] m ≤ [ψ k,n+1 ] m < [c k ] m [ψ k,n+1 ] m - µ k λ k ε + |[w (t) k ] m | , if [ψ k,n+1 ] m ≥ [c k ] m 2
Algorithm 3 Approximate proximal operator of reweighted ∞,1 -norm Initialization: Choose threshold value τ 2 > 0 and ε > 0.

Proximal operator: At each instant n ≥ 0, for each node k, utilize ψ k,n+1 to evaluate w k,n+1 in an element-wise manner:

1) Calculate [w o k ] m as the maximal value of [ψ ,n+1 ] m for all ∈ N - k ; 2) If [w o k ] m ≤ τ 2 , then calculate [w k,n+1 ] m as [w (t+1) k ] m via (49) with [w (t) k ] m = [w k,n ] m ; 3) If [w o k ] m > τ 2 , then calculate [w k,n+1
] m as ŵ via (53).

V. PERFORMANCE AND CONVERGENCE ANALYSES

In this section, we analyse the performance and convergence property of Algorithm 2 and Algorithm 3 in an unified framework. Quantities specifically related to ∞,1 -norm or reweighted ∞,1 -norm are distinguished by the superscripts (1) and (2) , respectively. Observe that the closed-forms expressions of both approximate proximal operators can be written compactly as:

prox µ k λ k ,gi(w k ) ψ k,n+1 = ψ k,n+1 -γ (i) k,n+1 , (54) with γ (i)
k,n+1 a (L × 1)-dimensional vector. For the ∞,1 -norm, the m-th entry of γ

(1)
k,n+1 is given by:

γ (1) k,n+1 m =      -µ k λ k , if ψ k,n+1 m < -µ k λ k ψ k,n+1 m , if ψ k,n+1 m ≤ µ k λ k µ k λ k , if ψ k,n+1 m > µ k λ k (55) if [w o k ] m ≤ τ 1 , and: γ (1)
k,n+1 m = (56)                -µ k λ k , if ψ k,n+1 m ≤ -[w o k ] m -µ k λ k ψ k,n+1 m +[w o k ] m , if -[w o k ] m -µ k λ k < [ψ k,n+1 ] m ≤-[w o k ] m 0, if [ψ k,n+1 ] m < [w o k ] m ψ k,n+1 m -[w o k ] m , if [w o k ] m ≤ [ψ k,n+1 ] m < [w o k ] m +µ k λ k µ k λ k , if ψ k,n+1 m ≥ [w o k ] m +µ k λ k if [w o k ] m > τ 1 .
For the reweighted ∞,1 -norm, the m-th entry of γ

k,n+1 is given by:

γ (2) k,n+1 m =                    -µ k λ k ε+|[w (t) k ] m | , if ψ k,n+1 m < -µ k λ k ε+|[w (t) k ] m | ψ k,n+1 m , if ψ k,n+1 m ≤ µ k λ k ε+|[w (t) k ] m | µ k λ k ε+|[w (t) k ] m | , if ψ k,n+1 m > µ k λ k ε+|[w (t) k ] m | (57) if [w o k ] m ≤ τ 2
, and:

γ (2) k,n+1 m = (58)                          - µ k λ k ε+|[w (t) k ] m | , if [ψ k,n+1 ] m ≤ -[c k ] m [ψ k,n+1 ] m + [w o k ] m , if -[c k ] m < [ψ k,n+1 ] m ≤ -[w o k ] m 0, if [ψ k,n+1 ] m < [w o k ] m [ψ k,n+1 ] m -[w o k ] m , if [w o k ] m ≤ [ψ k,n+1 ] m < [c k ] m µ k λ k ε+|[w (t) k ] m | , if [ψ k,n+1 ] m ≥ [c k ] m if [w o k ] m > τ 2 .
We observe in ( 55)-( 56) that: γ

(1)

k,n+1 m ≤ µ k λ k (59) Thus, γ (1)
k,n+1 is absolutely bounded: |γ 57)-(58), we conclude that: |γ

(1) k,n+1 | µ k λ k L . From (
(2) k,n+1 | b k L (60)
where b k is a constant defined as:

b k µ k λ k ε . ( 61
)
Substituting ( 54) into [START_REF] Piggott | Stability of adaptive network algorithms in multitask environments[END_REF], we obtain:

ψ k,n+1 = w k,n + µ k x k,n d k,n -x k,n w k,n w k,n+1 = ψ k,n+1 -γ (i) k,n+1 . (62)
We shall now analyse the performance of proximal multitask diffusion LMS algorithm based on expression (62). Define

wk,n+1 w k,n+1 -w k , (63) ψk,n+1 ψ k,n+1 -w k . ( 64
)
By collecting w k , w k,n+1 , ψ k,n+1 , wk,n+1 , ψk,n+1 , γ

k,n+1 over the entire network into block column vectors, we obtain quantities w , w n+1 , ψ n+1 , wn+1 , ψn+1 , γ (i) n+1 , respectively. To facilitate the theoretical analysis, we introduce the following assumptions on the regression data and step-size. These assumptions are widely used in the analysis of adaptive filters [START_REF] Sayed | Adaptive Filters[END_REF]- [START_REF] Chen | Nonnegative least-mean-square algorithm[END_REF] and diffusion networks [START_REF] Sayed | Diffusion adaptation over networks[END_REF], [START_REF] Sayed | Adaptation, Learning, and Optimization over Networks[END_REF], [START_REF] Nassif | Proximal multitask learning over networks with sparsity-inducing coregularization[END_REF].

A1 (Independent Regressors): The regression vector x k,n , generated from a zero-mean random process, is temporally stationary, white (over n) and spatially independent (over k) with covariance matrix R x,k = E{x k,n x k,n } > 0.

A2 (Small step-sizes): The step-sizes µ k of the network are small enough, so that terms on the higher-order powers of the step-sizes can be ignored.

A. Mean behavior analysis

Subtracting w k from both sides of the first equation of (62), and using signal model [START_REF] Tu | Diffusion strategies outperform consensus strategies for distributed estimation over adaptive netowrks[END_REF] and block notations, we obtain:

ψn+1 = B n wn + U h n , (65)
where

B n I -U M n (66) with U diag{µ 1 , µ 2 , • • • , µ N } ⊗ I L , (67)
M n diag x 1,n x 1,n , x 2,n x 2,n , • • • , x N,n x N,n , (68)
h n col{x 1,n z 1,n , x 2,n z 2,n , • • • , x N,n z N,n }. ( 69
)
Subtracting w k from both sides of the second equation of (62) and using block notations, we obtain:

wn+1 = ψn+1 -γ (i) n+1 . (70)
Combining ( 65) and (70) leads to:

wn+1 = B n wn + U h n -γ (i) n+1 . (71)
Using A1 and taking the expectation of (71), we obtain:

E{ wn+1 } = BE{ wn } -E{γ (i) n+1 }, (72)
where

B E{B n } = I -U M , (73)
E{γ (1) k,n+1 } m =                -µ k λ k , if E{ψ k,n+1 } m ≤-E{[w o k ] m }-µ k λ k E{ψ k,n+1 } m +E{[w o k ] m }, if -E{[w o k ] m }-µ k λ k < E{ψ k,n+1 } m ≤ -E{[w o k ] m } 0, if E{ψ k,n+1 } m < E{[w o k ] m } E{ψ k,n+1 } m -E{[w o k ] m }, if E{[w o k ] m }≤ E{ψ k,n+1 } m < E{[w o k ] m }+µ k λ k µ k λ k , if E{ψ k,n+1 } m ≥ E{[w o k ] m }+µ k λ k (77)
with

M E{M n } = diag R x,1 , R x,2 , • • • , R x,N , (74)
E{γ (i) n+1 } = col E{γ (i) 1,n+1 }, • • • , E{γ (i) N,n+1 } . ( 75
)
To analyze iteration (72), we need to derive the explicit expression of E{γ (i) n+1 } firstly. Since it has different expressions for ∞,1 -norm and reweighted ∞,1 -norm, we shall now evaluate it separately.

1) ∞,1 -norm: Taking the expectation of ( 55), [START_REF] Porter | The bellhop manual and user's guide: Preliminary draft[END_REF], and approximating the inequality conditions by conditions on the expectation terms, the m-th entry of E{γ [START_REF] Tu | Diffusion strategies outperform consensus strategies for distributed estimation over adaptive netowrks[END_REF] k,n+1 } is given by:

E{γ (1) k,n+1 } m =      -µ k λ k , if E{ψ k,n+1 } m < -µ k λ k E{ψ k,n+1 } m , if E{ψ k,n+1 } m ≤ µ k λ k µ k λ k , if E{ψ k,n+1 } m > µ k λ k (76) if E{[w o k ] m } ≤ τ 1 , and it is defined by (77) if E{[w o k ] m } > τ 1 .
To facilitate the derivation of the transient mean behavior in (76) and (77), we approximate the quantities appeared in "if" conditions by their corresponding expectations. It is difficult, if not impossible, to conduct theoretical analysis without this approximation. Note that this approximation is only used in the transient mean behavior analysis and have not been used in any other places. It is observed that the theoretical transient mean behavior obtained with this approximation match well with the results obtained via the Monte-Carlo simulation in Fig. 2 of Section VI-A. Vector E{ψ k,n+1 } is the k-th block vector of E{ψ n+1 }, which is evaluated as follows:

E{ψ n+1 } = E{ ψn+1 } + w , (78)
with:

E{ ψn+1 } = BE{ wn }. (79)
On the other hand, E{[w o k ] m } is approximated by the maximal value of E{ψ ,n+1 } m for all ∈ N - k . 2) Reweighted ∞,1 -norm: As the distribution of [w

(t)
k ] m is unknown, we cannot evaluate exactly the expectations involving [w

(t)

k ] m in (57) and (58). We therefore consider the first-order Taylor series expansion as in [START_REF] Arenas-Garcia | Mean-square performance of a convex combination of two adaptive filters[END_REF]- [START_REF] Silva | A transient analysis for the convex combination of two adaptive filters with transfer of coefficients[END_REF]; the performance results obtained with the first-order Taylor series expansion match well with the Monte-Carlo results in Fig. 3 of Section VI-A. By expanding f ([w

(t) k ] m ) µ k λ k /[ε + |[w (t) k ] m |] around E{[w (t)
k ] m } and taking expectation, we obtain the following approximation:

E µ k λ k ε + |[w (t) k ] m | ≈ µ k λ k ε + |E{[w (t) k ] m }| . ( 80
)
Taking the expectation of (57), (58) and using (80), as well as approximating the inequality conditions by conditions on the expectation terms, we obtain the explicit expression of E{γ

k,n+1 }, with its m-th entry given by (83) or (84) when

E{[w o k ] m } ≤ τ 2 and E{[w o k ] m } > τ 2 , respectively, with quantities: E{[w (t) k ] m } = E{w k,n } m (81) E{[c k ] m } = E{[w o k ] m } + µ k λ k ε+ E{w k,n } m . ( 82
)
Note that several approximations, such as the first-order Taylor series expansion [START_REF] Arenas-Garcia | Mean-square performance of a convex combination of two adaptive filters[END_REF]- [START_REF] Silva | A transient analysis for the convex combination of two adaptive filters with transfer of coefficients[END_REF] and approximated the inequality conditions, are only used in the transient mean behavior of the reweighted ∞,1 -norm, not in any other places. It is observed that the theoretical transient mean behavior obtained with these approximations match well with the results obtained via the Monte-Carlo simulation in Fig. 3 of Section VI-A.

From iteration (72), we obtain the following Theorem 1 for the mean stability of proximal multitask diffusion LMS algorithm [START_REF] Piggott | Stability of adaptive network algorithms in multitask environments[END_REF].

E{γ (2) k,n+1 } m =                    -µ k λ k ε+ E{[w (t) k ] m } , if E{ψ k,n+1 } m < -µ k λ k ε+ E{[w (t) k ] m } E{ψ k,n+1 } m , if E{ψ k,n+1 } m ≤ µ k λ k ε+ E{[w (t) k ] m } µ k λ k ε+ E{[w (t) k ] m } , if E{ψ k,n+1 } m > µ k λ k ε+ E{[w (t) k ] m } (83) E{γ (2) k,n+1 } m =                          - µ k λ k ε+ E{[w (t) k ] m } , if E{ψ k,n+1 } m ≤ -E{[c k ] m } E{ψ k,n+1 } m +E{[w o k ] m }, if -E{[c k ] m } < E{ψ k,n+1 } m ≤ -E{[w o k ] m } 0, if E{ψ k,n+1 } m < E{[w o k ] m } E{ψ k,n+1 } m -E{[w o k ] m }, if E{[w o k ] m }≤ E{ψ k,n+1 } m < E{[c k ] m } µ k λ k ε+ E{[w (t) k ] m } , if E{ψ k,n+1 } m ≥ E{[c k ] m } (84)
Theorem 1. (Mean stability) Assume data model [START_REF] Tu | Diffusion strategies outperform consensus strategies for distributed estimation over adaptive netowrks[END_REF] and assumption A1 hold. Then for any initial conditions, distributed networks endowed with the proximal multitask diffusion LMS algorithm [START_REF] Piggott | Stability of adaptive network algorithms in multitask environments[END_REF] are stable in the mean, if step-sizes µ k satisfy:

0 < µ k < 2 λ max {R x,k } , k = 1, • • • , N, (85)
where λ max {•} denotes the maximal eigenvalue of its matrix argument. The block maximum norm of the bias can be upper bounded as:

lim n→∞ E{ wn+1 } b,∞ ≤ √ L • max k µ k λ k 1 -B b,∞ (86)
lim n→∞ E{ wn+1 } b,∞ ≤ 1 ε √ L • max k µ k λ k 1 -B b,∞ (87)
for the ∞,1 -norm and reweighted ∞,1 -norm, respectively. Proof: By iterating the RHS of (72) from time instant n = 0, and proving the convergence of the obtained series, we arrive at condition (85) for step-size to ensure the mean stability of [START_REF] Piggott | Stability of adaptive network algorithms in multitask environments[END_REF]. For more details, see Appendix D. Remark 1: Equation (85) provides an upper bound for step-size µ k to ensure the mean stability of the distributed networks with the proximal multitask diffusion LMS algorithm [START_REF] Piggott | Stability of adaptive network algorithms in multitask environments[END_REF]. The upper bound is closely related to the second-order statistics R x,k of the input signals. Equations ( 86) and (87) indicate that the proximal multitask diffusion LMS algorithm ( 19) is biased. The upper bound of the biases are proportional to the length of the system vector L, the step-size µ k and the regularization parameter λ k . The bias can be reduced by using a sufficiently small step-size µ k or regularization parameter λ k . Besides, for the reweighted ∞,1 -norm, the bias is also inversely proportional to parameter ε. In addition, the upper bound of bias for the reweighted ∞,1 -norm has an improvement with a factor 1 ε than that of the ∞,1 -norm.

B. Mean-square behavior analysis

Under A1 and using (71), then for any semi-positive definite matrix Σ of compatible dimension, the weighted mean-square behavior of wn+1 evaluates as:

E{ wn+1 2 Σ } = E{ wn 2 Σ } + E{ U h n 2 Σ } + E{ γ (i) n+1 2 Σ } -2E{ w n B n Σγ (i) n+1 } -2E{h n U Σγ (i) n+1 }, (88)
where x 2

Σ

x Σx, and

Σ E{B n ΣB n }. ( 89
)
Let σ vec{Σ} and σ vec{Σ }, where vec{•} operator stacks the columns of its matrix argument on top of each other. Using the property of vec{•} operator, (89) becomes:

σ = E{B n ⊗ B n }σ. ( 90
)
Under A2 and ignoring terms on the second-order of the maximal step-size, we have the approximation for (90):

σ ≈ F σ (91)
with F B ⊗ B . Define:

H U diag σ 2 z,1 R x,1 , • • • , σ 2 z,N R x,N U . ( 92
)
We then have:

E{ U h n 2 Σ } = [vec{H}] σ. ( 93
)
To make the analysis tractable, we adopt approximation:

E{ w n B n Σγ (i) n+1 } ≈ E{ w n B Σγ (i) n+1 }. ( 94
)
Since |γ

(1)

k,n+1 | µ k λ k L and |γ (2) k,n+1 | µ k λ k ε L , we conclude that γ (i)
n+1 is at most of the same order as the step-size. This implies that the last term on the RHS of (88) contains higher-order powers of the step-size, and can be ignored according to assumption A2. Finally, by using ( 91), ( 93) and (94), expression (88) becomes:

E{ wn+1 2 σ } = E{ wn 2 F σ }+[vec{H}] σ + E{ γ (i) n+1 2 σ } -2E{ w n B Σγ (i) n+1 }, ( 95
)
where we use the notations E{ wn+1 Theorem 2. (Mean-square stability) Assume data model ( 1) and assumptions A1, A2 hold. Further assume that approximation (91) is reasonable for sufficiently small step-sizes. Then for any initial conditions, distributed networks endowed with proximal multitask diffusion LMS algorithm [START_REF] Piggott | Stability of adaptive network algorithms in multitask environments[END_REF] is stable in the mean-square sense, if the step-sizes µ k are sufficiently small and satisfy (85).

Proof: By proving that the last two terms on the RHS of ( 95) is bounded, and iterating (95) from time instant n = 0, we obtain the condition for step-size to ensure the mean-square stability of [START_REF] Piggott | Stability of adaptive network algorithms in multitask environments[END_REF]. For more details, see Appendix E. Remark 2: The weighted mean-square behavior of wn+1 , that is E{ wn+1 2 Σ }, can be decomposed as:

E wn+1 2 Σ = E w n+1 -w -E{w n+1 }+E{w n+1 } 2 Σ = E wn+1 -E{ wn+1 } 2 Σ Variance term + E{ wn+1 } Bias term 2 Σ (96)
By substituting wn+1 of (71) and E{ wn+1 } of (72) into (96), we obtain an equivalent form of (95). Relation (96) is called the bias-variance decomposition. From (96), we observe that the stability of E{ wn+1 2 Σ } in Theorem 2 ensures the stabilities of both the Variance term and the Bias term E{ wn+1 } of Theorem 1.

VI. SIMULATION RESULTS

In this section, we present simulation results to validate the effectiveness of the algorithm. With the exception of the simulation results presented in Section VI-A, used to validate theoretical results in the mean behavior analysis and obtained by averaging over 500 independent Monte-Carlo runs, all other simulated curves were obtained by averaging over 100 independent Monte-Carlo runs.

A. Theoretical validation

We considered a connected network consisting of 16 nodes and 36 edges. The number of edges at each node was between 2 and 7, without taking into account the possible self-loops connecting each node to itself. Other characteristics of the network are listed in Table I. Each regressor x k,n was generated from a zero-mean Gaussian distribution with covariance matrix R x,k = σ 2

x,k I 30 . Each additive noise z k,n was generated from a zero-mean Gaussian distribution with variance σ 2 z,k . Variances σ 2

x,k and σ 2 z,k at each node were generated randomly from a Gaussian distribution as shown in Fig. 1. Note that these variances settings are the same as in [START_REF] Chen | Multitask diffusion adaptation over networks[END_REF], [START_REF] Chen | Diffusion LMS over multitask networks[END_REF].

To validate the theoretical results reported in the mean behavior analysis, we considered a stationary system identification scenario. The unknown system coefficients w k were generated such that the entire network has a jointly sparse structure with sparsity degree of 10/30. Each nonzero element of w k was generated independently from a standard Gaussian distribution. The simulation results presented in this section were obtained by averaging over 500 independent runs. The results are illustrated in Fig. 2 and Fig. 3 for the proximal multitask LMS with ∞,1 -regularization and reweighted ∞,1 -regularization, respectively. The good match between the theoretical results and the Monte-Carlo curves illustrates the accuracy of theoretical results in the mean behavior analysis.

We also considered the average error over the entire network defined by: Though there is a small bias between theoretical and Monte-Carlo curves, the results reported in Fig. 4 and Fig. 5 confirm the accuracy of the theoretical models.

E av (n) 1 N N k=1 w k,n -w k (97)

B. Numerical Simulations 1) Comparison with existing algorithms:

We firstly considered a non-stationary jointly sparse system identification scenario with w k varying over time. Each nonzero entry of w k was generated independently from a standard Gaussian distribution. The evolution of w k was divided into four stationary stages and three transient stages. During stationary stages, sparse vectors w k were set to sparsity degree of 3/30, 5/30, 8/30 and 10/30, respectively. The transient stages were designed by using linear interpolation over 500 time instants. The regressors x k,n were generated as those in Section VI-A for white inputs, while they were generated according to a zero-mean Gaussian distribution with covariance matrix R x,k = σ 2

x,k R † for colored inputs, where R † is an 30×30 Hermite matrix with eigenvalue spread λ max {R † }/λ min {R † } = 21, and symbol λ min {•} denotes the minimal eigenvalue of its matrix argument. For comparison purpose, non-cooperative diffusion LMS algorithm, non-cooperative sparse diffusion LMS [START_REF] Lorenzo | Sparse distributed learning based on diffusion adaptation[END_REF] with zero-attracting (ZA) regularizer and reweighted zero-attracting (RZA) regularizer, multitask diffusion LMS with adaptive combiner [START_REF] Chen | Diffusion LMS over multitask networks[END_REF] and jointly sparse multitask diffusion LMS [START_REF] Li | Distributed jointly sparse multitask learning over networks[END_REF] with 2,1 -regularization, reweighted 2,1regularization (RW 2,1 ) and 2,0 -regularization were taken into consideration. We adopted a uniform step-size 0.01 for all algorithms. For these algorithms, we set their parameters so that they reach their best performance. To enable reproducible research, details about parameters used by these algorithms can be found in Table III of Appendix F.

The results are illustrated in Fig. 6 for white inputs. We observe that multitask LMS with adaptive combiner is the worst one among all competing algorithms, since it utilizes similarities between neighboring nodes to improve estimation accuracy. This does not necessarily exist in jointly sparse scenarios and may deteriorate the MSD performance. Since jointly sparse system can be regarded as a special case of general sparse systems, by using additional information about system sparsity, sparse diffusion LMS with ZA regularizer and RZA regularizer have better performance than the non-cooperative LMS. Similarly, all jointly sparse multitask algorithms considered in this comparative experiment perform better than the non-cooperative LMS. Observe that sparse diffusion LMS algorithms can perform slightly better than jointly sparse diffusion LMS algorithms as the parameter vectors to estimate become sparser, while they do not perform as well when these parameter vectors are less sparse. These findings illustrate the interest of exploiting joint sparsity as prior information. On the one hand, the proposed proximal multitask LMS with reweighted ∞,1 -regularization performs better than all other algorithms as evidenced by its lowest steadystate MSD. On the other hand, the proposed proximal multitask LMS with ∞,1 -regularization has similar performance to the sparse diffusion LMS with RZA regularizer when there are more zeros in the jointly sparse system to estimate, but the former performs better than the latter as the number of nonzero entries increases.

The results are illustrated in Fig. 7 for colored input signals. It can be observed that the convergence rate of all algorithms is slower than in the case of white inputs, and some algorithms have poorer performance than the non-cooperative algorithm when the parameter vectors to estimate are less sparse. Besides some conclusions that have been drawn for white inputs, we observe that the proximal multitask LMS algorithm with reweighted ∞,1 -regularization still has the best performance.

Then we examined the performance of all algorithms as a function of the sparsity degree of the jointly sparse system to estimate. We considered the steady-state MSD as the measure of performance in this experiment. The results are illustrated in Fig. 8. Besides some conclusions that have been drawn before, we observe that some of the multitask algorithms with sparse and jointly sparse regularizers may have poorer performance than the non-cooperative LMS when the number of nonzero entries increases. Three algorithms, including the proximal multitask LMS with reweighted ∞,1 -regularization, uniformly show better performance than all other algorithms for all sparsity degrees, and have similar performance to the non-cooperative LMS for totally non-sparse systems.

2) Effects of parameters setting: To examine the effects of parameters setting, including the step-size µ k , the regularization parameter λ k , the threshold values τ 1 , τ 2 and the parameter ε, we considered a stationary system identification problem with sparsity degree of 3/30. Each parameter was set to a same value for all nodes in the network. We examined the influence of one selected parameter at a time, setting all other parameters to fixed values, in order to facilitate comparison.

The effects of the step-size µ k are illustrated in Fig. 9 and Fig. 10 for the proximal LMS with ∞,1 -regularization and reweighted ∞,1 -regularization, respectively. Observe that µ k allows to control the trade-off between convergence rate and steady-state performance. A larger step-size results in a faster convergence rate at the cost of a larger steady-state MSD. A small step-size results in a more accurate estimation at the cost of a slower convergence speed. The effects of the regularization parameter λ k are illustrated in Fig. 11 and Fig. 12. We observe that increasing λ k improves the convergence speed and the MSD at steady-state at first, and then degrades them. In this experiment, the critical values were 0.08 and 0.01 for the ∞,1 regularizer and the reweighted ∞,1 regularizer, respectively.

The effects of parameters τ 1 and τ 2 are shown in Fig. 13. Similar behaviors can be observed on the MSD at steady-state of both algorithms with respect to τ 1 and τ 2 . The best values were obtained over interval [0.05, 0.8] for both regularizers. First, these results show the need for introducing these two parameters. Second, they show that these two parameters can be appropriately selected over a large interval. Consider the ∞,1 -norm regularizer. When [w o k ] m ≤ τ 1 , observe that the approximate proximal operator is given by [START_REF] Wee | A proximal splitting approach to regularized distributed adaptive estimation in diffusion networks[END_REF] and corresponds to the 1 -norm regularizer, which enjoys the following properties: On the one hand, it shrinks the estimates of zero-valued parameters to zero when [w o k ] m ≤ τ 1 ; On the other hand, for small nonzero-valued entries, though introducing a bias when shrunk to zero, it lowers the estimates variance and results in a satisfying performance. Since jointly sparse systems are a special case of general sparse ones, for which the 1 -norm regularizer is prescribed, it is expected that our algorithm works well for a large range of τ 1 values. The same reasoning applies to the reweighted ∞,1 -norm regularizer and τ 2 .

The effects of parameter ε in the reweighted ∞,1 regularizer are illustrated in Fig. 14. We set both parameters µ k and λ k to 0.01. We observe that the best value of ε resulting in the lowest MSD at steady-state was 0.002, but its effect on the steady-state MSD is weaker than the other parameters.

C. Practical application

We shall now validate the proposed algorithm in a practical application. This application follows the experimental setup described in [START_REF] Li | Distributed jointly sparse multitask learning over networks[END_REF]. In demodulation and decoding in underwater and wireless communication networks, a standalone low-cost sensor may not be able to decode or demodulate the source signal reliably due to the extra low signal-to-noise ratio (SNR) condition. It is better to adopt multiple sensors to collaboratively recover the information coming from a same source [START_REF] Zhu | Distributed in-network channel decoding[END_REF], [START_REF] Zhu | Distributed consensus-based demodulation: algorithms and error analysis[END_REF] by combining and exchanging information within the sensor network. Practically, it is necessary to estimate the channel impulse responses between the source and each sensor before decoding. On the one hand, it has been shown that real-world underwater channels [START_REF] Kocic | Sparse equalization for realtime digital underwater acoustic communications[END_REF], and wireless communication channel [START_REF] Masood | Efficient coordinated recovery of sparse channels in massive MIMO[END_REF], are sparse with large delay spread. On the other hand, the channel supports for neighboring antennas or nodes are approximately the same [START_REF] Masood | Efficient coordinated recovery of sparse channels in massive MIMO[END_REF]. Indeed the times of arrival for closely spaced nodes and antennas are quite close, though the tap weights are different [START_REF] Li | Distributed jointly sparse multitask learning over networks[END_REF]. We shall now check that taking the jointly sparse property into consideration can improve the network performance.

Consider the problem of identifying underwater acoustic channels. We set the scenario presented in Fig. 15, where we have one source (the submarine) and two receiving sensor arrays. Both arrays are linear arrays and consist of 10 with the same distance (3.75 m) between each of them. We assume that each node has the ability to communicate and process data. These two arrays are positioned in parallel with interspace 50 m, which us the distributed network with 20 nodes. Within each array, the nodes are connected one by one in a chain, and the corresponding pairs of nodes between two arrays are also connected. The underwater acoustic channels to estimate were generated via the BELLHOP model [START_REF] Porter | The bellhop manual and user's guide: Preliminary draft[END_REF], which has been developed for predicting acoustic pressure fields in ocean environments. We used a white signal as input, and the additive noise at each node was white Gaussian. SNR conditions are listed in Table II. We compared our algorithms with the algorithms considered in Section VI-B1. Parameters used by all them are listed in Table III. The results are depicted in Fig. 16. As can be observed, the proposed proximal multitask diffusion LMS algorithm with reweighted ∞,1 regularizer achieved the best performance in terms of the steady-state MSD.

VII. CONCLUSION

We considered the problem of estimating a set of parameter vectors in a distributed manner, where the local solutions have the same sparse support. We devised a proximal diffusion algorithm with (reweighted) ∞,1 -norm regularization, with closed-form expressions for the regularizers. We conducted theoretical analyses of the algorithms behavior in the mean and mean-square sense. Simulation results illustrated the effectiveness of the proposed algorithms, as well as the accuracy of theoretical results. Integrating weighted network connection information to enhance this jointly-sparse estimation will be considered in future work.

Using (122)-(125), we have:

φ wn , γ (i) n+1 ≤ E{ γ (i) n+1 2 σ }+ 2E{ w n B Σγ (i) n+1 } ≤ κ (i) 1 + τ (i) max • κ (i) 2
(126) for all n. Given a weighting matrix Σ, the positive constant κ where E{ w0 2 } is the initial condition. The stability of E{ wn+1 2 σ } requires the convergence of terms on the RHS of (129), which is ensured by condition ρ(F ) < 1. Since F = B ⊗ B , it is enough to select sufficiently small step-sizes µ k satisfying (85) to ensure ρ(F ) < 1.

APPENDIX F PARAMETERS USED BY ALGORITHMS IN SECTION VI-B1

Table III provides all parameters used by the algorithms in Section VI-B1. Notations used for these parameters are the same as those in the original references. Some pairs of columns, standing for different parameters, are merged into a single column for compactness. The corresponding symbols can be distinguished by the symbol " | ". We used a uniform step-size 0.01 for all algorithms.

  [ • ] m denotes the m-th entry of its vector argument. The superscript (•) denotes the transpose operator. The mathematical expectation is denoted by E{•}. The Gaussian distribution with mean µ and variance σ 2 is denoted by N (µ, σ 2 ). Operator | • | takes the absolute value of its scalar or vector argument. Operator max{•,•} extracts the maximum value of its two arguments. Operator diag{•} generates a diagonal matrix from its argument. Symbol denotes a component-wise inequality. The symbol ⊗ denotes the Kronecker product. The set N k denotes the neighbors of node k, including k itself, and |N k | denotes its cardinality. The N - k denotes the neighbors of node k, excluding node k. Vector L is the all-one vector of dimension L × 1.

.

  Unlike the proximal operator (53) in Case B, only the proximal operator (49) in Case A has the capability to drive [w k ] m to zero and promote sparsity. Thus, as in Section IV-A, we relax the condition [w o k ] m = 0 by introducing a small positive tolerance τ 2 to distinguish between zero and nonzero entries. This leads to the condition [w o k ] m ≤ τ 2 for Case A, and [w o k ] m > τ 2 for Case B.

2 Σ } and E{ wn+1 2 σ

 22 } interchangeably. From iteration (95), we obtain the following Theorem 2 to ensure the mean-square stability of the proximal multitask diffusion LMS algorithm[START_REF] Piggott | Stability of adaptive network algorithms in multitask environments[END_REF].

Fig. 1 .Fig. 2 .Fig. 3 .Fig. 4 .

 1234 Fig. 1. Agent input and noise variances.

Fig. 5 .Fig. 6 .

 56 Fig. 5. The average error Eav(n) of the reweighted ∞,1 -regularization; Theory: red solid line; Monte-Carlo simulations: blue dashed line.

Fig. 7 .

 7 Fig. 7. Comparison of the proposed algorithms with several state-of-the-art algorithms for colored inputs.

Fig. 8 .

 8 Fig. 8. Steady-state MSDs of all algorithms as a function of the sparsity degree.

Fig. 9 .

 9 Fig. 9. Effect of the step-size µ k on the MSD of the proximal LMS algorithm with ∞,1 -regularization.

Fig. 10 .

 10 Fig. 10. Effect of the step-size µ k on the MSD of the proximal LMS algorithm with reweighted ∞,1 -regularization.

Fig. 11 .

 11 Fig. 11. Effect of the regularization parameter λ k on the MSD of the proximal LMS algorithm with ∞,1 -regularization.

Fig. 12 .

 12 Fig. 12. Effect of the regularization parameter λ k on the MSD of the proximal LMS algorithm with reweighted ∞,1 -regularization.

Fig. 13 .

 13 Fig.[START_REF] Nassif | Multitask diffusion LMS with sparsity-based regularization[END_REF]. Effect of parameters τ 1 and τ 2 on the steady-state MSD of the proximal LMS algorithm with ∞,1 -regularization (black) and reweighted ∞,1regularization (red), respectively.

Fig. 14 .

 14 Fig. 14. Effect of parameter ε on the steady-state MSD of the proximal LMS algorithm with reweighted ∞,1 -regularization.

Fig. 15 .Fig. 16 .

 1516 Fig. 15. Experimental setup for the practical application.

2 σ } as: E{ wn+1 2 σ }≤ E{ wn 2 F 2 F 2 σ }≤ E{ w0 2 F

 222222 written as a scaled multiple of [vec{H}] σ as:κ (i) 3 = p • [vec{H}] σ(127)with p ≥ 0[START_REF] Lorenzo | Sparse distributed learning based on diffusion adaptation[END_REF]. Using (95) and (125)-(127), we obtain an upper bound of E{ wn+1 σ }+[vec{H}] σ+ φ wn , γ(i) n+1 ≤ E{ wn σ }+(1 + p)• [vec{H}] σ.(128)Iterating (128) from n = 0, we obtain:E{ wn+1

TABLE I CHARACTERISTICS

 I OF THE NETWORK USED FOR MODEL VALIDATION. L IS THE LAPLACIAN MATRIX ASSOCIATED WITH THE GRAPH,

TABLE II SNR

 II LEVEL IN DECIBEL (DB) FOR THE PRACTICAL APPLICATION. SINCE IT VARIES ACCORDING TO NODES, WE ENUMERATE THE MAXIMUM, MINIMUM AND MEAN VALUES.

	SNR Level Maximum Minimum Mean
	SNR	7.89	6.19	7.07

TABLE III PARAMETERS

 III USED BY ALL ALGORITHMS FOR SIMULATION.

	Algorithms proposed ∞,1 proximal LMS proposed RW ∞,1 proximal LMS ZA-Sparse LMS [6] RZA-Sparse LMS [6] 2,0 jointly sparse [24] 2,1 jointly sparse [24] RW 2,1 jointly sparse [24]	λ k | γ 0.08 0.01 | 0.03 | 0.03	ε 0.01 0.45 8 100	η 0.1 0.03 1.5	τ 1 | τ 2 0.1 | 0.05

APPENDIX A INTERPRETATIONS ABOUT THE (REWEIGHTED) ∞,1 -NORM

The ∞,1 -norm of matrix W k is defined as g 1 (w k ) of [START_REF] Mao | Walkman: A communication-efficient random-walk algorithm for decentralized optimization[END_REF]. We focus on the interpretation of wk,m ∞ . From [START_REF] Sayed | Adaptive networks[END_REF]

• Case c: The m-th row wk,m corresponds to the non-zero valued entries, and k ] m for all m. Similar interpretations can be obtained for the reweighted ∞,1 -norm of [START_REF] Chen | Multitask diffusion adaptation over networks[END_REF]. APPENDIX B DERIVATION OF [START_REF] Huang | The benefit of group sparsity[END_REF] The initial problem [START_REF] Li | Distributed jointly sparse multitask learning over networks[END_REF] leads to [START_REF] Bazerque | Distributed spectrum sensing for cognitive radio networks by exploiting sparsity[END_REF] and [START_REF] Kocic | Sparse equalization for realtime digital underwater acoustic communications[END_REF] in Case 1 and Case respectively. For ease of presentation, we define:

We select solution [START_REF] Duarte | Distributed compressed sensing of jointly sparse signals[END_REF] or [START_REF] Vlaski | Diffusion stochastic optimization with non-smooth regularizers[END_REF] depending on the value taken by costs (98) and (99). To save space, we partially present the derivation. The rest of the derivation can be obtained by following the same routine.

, we obtain:

] m -µ k λ k of Case 2 into (99), we obtain:

Thus the proximal operator is given by:

we arrive at proximal operator (103) directly.

Substituting ŵ = [w o k ] m of Case 2 into (99), we obtain:

Thus the proximal operator is given by:

By following the same routine, we obtain the proximal operator for [ψ k,n+1 ] m ≤ 0. Finally, by combining all these results, we arrive at the expression of proximal operator [START_REF] Huang | The benefit of group sparsity[END_REF] when

APPENDIX C DERIVATION OF [START_REF] Simões | FADE: Fast and asymptotically efficient distributed estimator for dynamic networks[END_REF] Define:

The initial problem (37) leads to J3 ([w k ] m ) and J4 ([w k ] m ) in Case 1 and Case 2, respectively. We select solution ( 39), ( 49) or ( 51) by comparing costs (107) and (108) as that in Section IV-A. We partially present the derivation to save space.

When

] m of Case 1 into (107), we obtain:

Substituting ŵ = [w o k ] m of Case 2 (108), we obtain:

This means that the proximal operator is given by:

, since ŵ of both Case 1 and Case 2 is given by: ŵ

we arrive at proximal operator (112) directly.

of Case 1 into (107), we obtain:

of Case 2 is an approximation of the minimizer ŵo of (41), and according to the relation J4 ( ŵo ) ≤ J4 ([w o k ] m ), we define the proximal operator as:

By following the same routine, we obtain the proximal operator for [ψ k,n+1 ] m ≤ 0. Finally, by combining all these results, we arrive at expression [START_REF] Simões | FADE: Fast and asymptotically efficient distributed estimator for dynamic networks[END_REF] when

APPENDIX D PROOF OF THEOREM 1

Iterating (72) from n = 0, we obtain:

where E{ w0 } is the initial condition. The convergence of (115) requires that both terms on the RHS to be convergent. For the first term, it requires that spectral radius ρ(B) < 1 to ensure the convergence. For the second term, it is sufficient to prove that

A series is absolutely convergent if each term is bounded by a term of an absolutely convergent series [START_REF] Lorenzo | Sparse distributed learning based on diffusion adaptation[END_REF], [START_REF] Nassif | Proximal multitask learning over networks with sparsity-inducing coregularization[END_REF]. Define s m B j E{γ (i) n+1-j } m . Since the block maximum norm of a block vector is larger than or equal to the largest absolute value of its entry, we have:

where • b,∞ is the block maximum norm [START_REF] Sayed | Diffusion adaptation over networks[END_REF]. The quantity E{γ (i) n+1-j } b,∞ is finite for all j and n, and bounded by some constant γ (i) max . Actually, from (76), ( 77), ( 83), (84) and following the routine for the boundness of γ (i) k,n+1 , we obtain:

where the quantity b k has been defined in (61). Thus, we conclude that E{γ

Consequently, convergence of ( 116) is ensured by condition ρ(B) < 1.

Step-sizes µ k satisfying (85) ensure the mean stability of the network. With step-sizes µ k satisfying (85) to ensure ρ(B) < 1, the block maximum norm of the bias can be bounded as n → ∞:

Substituting ( 119) and (120) into (121), we arrive at (86) and (87), respectively.

APPENDIX E PROOF OF THEOREM 2

Since Σ is a positive semi-definite matrix, and vector γ (i) n+1 is uniformly bounded for all time instant n and i = 1 or 2, we have:

for all n, where κ