Decentralized Sparse Multitask RLS over Networks

# Decentralized Sparse Multitask RLS over Networks

Xuanyu Cao and K. J. Ray Liu, Fellow, IEEE
Email: {apogne, kjrliu}@umd.edu
Department of Electrical and Computer Engineering, University of Maryland, College Park, MD
###### Abstract

Distributed estimation, decentralized optimization, adaptive networks, recursive least squares, sparsity, multitask, alternating direction method of multipliers, subgradient method

## I Introduction

In the last decade, distributed adaptive signal processing has emerged as a vital topic because of the vast applications in need of decentralized real-time data processing over networked systems. For multi-agent networks, distributed adaptive algorithms only rely on local information exchange, i.e., information exchange among neighbor nodes, to estimate the unknowns. This trait endows distributed adaptive algorithms with low communication overhead, robustness to node/link failures and scalability to large networks. In the literature, the centralized least mean squares (LMS) and recursive least squares (RLS) [1] have been extended to their decentralized counterparts [2, 3] to deal with estimation problems over networks. Furthermore, many natural signals are inherently sparse with most entries equal to zero such as the image signals and audio signals in [4, 5, 6, 7]. Sparsity of signals are particularly conspicuous in the era of big data: for many applications, redundant input features (e.g., a person’s salary, education, height, gender, etc.) are collected to be fed into a learning system to predict a desired output (e.g., whether a person will resign his/her job). Most input features are unrelated to the output so that the weight vector between the input vector and the output is highly sparse. As such, several sparse adaptive algorithms have been proposed such as the sparse LMS in [8, 9], the sparse RLS in [10] and the distributed sparse RLS in [11].

Most of the decentralized sparse adaptive algorithms are focused on the single task estimation problem, in which all nodes receive data associated with the same unknown vector and collaborate to estimate it. On the contrary, many applications are inherently multitask-oriented, i.e., each node has its own unknown vector different from others’. For instance, in a sensor network, each node may want to estimate an unknown vector related to its specific location and thus different nodes have different unknown vectors to be estimated. In fact, several decentralized multitask adaptive algorithms have been proposed in the literature including the multitask diffusion LMS in [12], its asynchronous version in [13] and its application in the study of tremor in Parkinson’s disease [14]. In particular, a sparse multitask LMS algorithm is proposed in [15] to promote sparsity of the estimated multitask vectors.

To the best of our knowledge, all the existing distributed adaptive algorithms for multitask estimation problems are based on various decentralized versions of LMS. The RLS based sparse multitask estimation problems have not aroused much attention. It is well known that the RLS possesses much faster convergence speed than the LMS. Hence, the RLS is more suitable for applications in need of fast and accurate tracking of the unknowns than the LMS, especially when the devices are capable of dealing with computations of moderately high complexity (which is the case as the computational capability of devices is increasing drastically). This motivates us to study the decentralized sparse multitask RLS problem over networks. The main contributions of this paper are summarized as follows.

• A global networked RLS minimization problem is formulated. In accordance with the multitask nature of the estimation problem, each node has its own weight vector. Since neighbor nodes often share analogous properties and thus similar weight vectors, we add regularization term to penalize deviations of neighbors’ weight vectors. To enforce sparsity of the weight vectors, we further introduce regularization.

• A decentralized online alternating direction method of multipliers (ADMM) algorithm is proposed for the formulated sparse multitask RLS problem. The proposed ADMM algorithm is simplified so that each iteration consists of simple closed-form computations and each node only needs to store and update one matrix and six dimensional vectors, where is the dimension of the weight vectors. We show that the gaps between the outputs of the proposed ADMM algorithm and the optimal points of the formulated RLS problems converge to zero.

• To overcome the relatively high computational cost of the proposed ADMM algorithm, we further present a decentralized online subgradient method, which enjoys lower computational complexity. We theoretically analyze its convergence behaviors and show that the tracking error of the weight vectors is upper bounded by some constant related to the network topology and algorithm parameters.

• Numerical simulations are conducted to corroborate the effectiveness of the proposed algorithms. Their advantages over the single task sparse RLS algorithm in [11] are highlighted. We also observe an accuracy-complexity tradeoff between the proposed two algorithms.

The roadmap of the remaining part of this paper is as follows. In Section \@slowromancapii@, the sparse multitask RLS problem is formally formulated. In Section \@slowromancapiii@, we propose and simplify a decentralized online ADMM algorithm for the formulated RLS problem. In Section \@slowromancapiv@, we propose a decentralized online subgradient method for the formulated problem in order to reduce computational complexity. In Section \@slowromancapv@, numerical simulations are conducted. In Section \@slowromancapvi@, we conclude this work.

## Ii The Statement of the Problem

We consider a network of nodes and some edges between these nodes. We assume that the network is a simple graph, i.e., the network is undirected with no self-loop and there is at most one edge between any pair of nodes. Denote the set of neighbors of node (those who are linked with node by an edge) as . The network can be either connected or disconnected (there does not necessarily exist a path connecting every pair of nodes). Time is divided into discrete slots denoted as . Each node has an unknown (slowly) time-variant dimensional weight vector to be estimated. The formulated network is therefore a multitask learning network since different nodes have different weight vectors, as opposed to the traditional single task learning network [2], which is usually transformed into a consensus optimization problem framework [16, 17, 18]. Each node has access to a sequence of private measurements , where is the input regressor at time and is the output observation at time . The measurement data are private in the sense that node has access only to its own measurement sequence but not others’. The data at node are assumed to conform to a linear regression model with (slowly) time-variant weight vector :

 dn(t)=uTn(t)˜wn(t)+en(t), (1)

where is the output measurement noise at time . In multitask learning networks, the benefit of cooperation between nodes comes from the fact that neighboring nodes have similar weight vectors [12], where similarity is embodied by some specific distance measures. By incorporating terms promoting similarity between neighbors and enforcing cooperation in the network, an estimator may achieve potentially higher performance than its non-cooperative counterpart.

Moreover, many signals in practice are highly sparse, i.e., most entries in the signal are equal to zero, with examples encompassing image signals, audio signals, etc. [4, 5, 6, 7]. The sparsity of signals is especially conspicuous in today’s big data era because redundant data are collected as input features among which most are unrelated to the targeted output, leading to sparsity of the corresponding weight vectors. Furthermore, as per convention in adaptive algorithms [1], we assume that the weight vectors varies with time very slowly. This suggests that past data are of great merit to estimate the current weight vector, which justifies the advantage of the RLS (studied in this paper) over the LMS (studied in all existing works on multitask estimation [12, 13, 15, 14]) in terms of convergence speed.

In all, we propose an RLS based estimator to track the unknown weight vectors while enforcing similarity between neighbors’ weight vectors and sparsity of all weight vectors. The estimator at time is the optimal solution of the following optimization problem:

 Minimizew1,...,wN  N∑n=1T∑t=1λT−t(dn(t)−uTn(t)wn)2+βN∑n=1∑m∈Ωn∥wn−wm∥22+γN∑n=1∥wn∥1, (2)

where are the forgetting factor of the RLS algorithm, regularization coefficient for similarity between neighbors’ weight vectors and regularization coefficient for sparsity, respectively. If , then problem (2) enforces consensus of weight vectors across nodes, and thus degenerates to the sparse RLS problem in [11]. Note that the measurement data arrives in a sequential manner, which necessitates an online (real time) algorithm to solve (2) due to the prohibitive computation and storage cost of offline methods. Further note that the private measurement data are distributed among network nodes. Thus, a distributed algorithm for (2) is imperative as centralized algorithms are vulnerable to link failures and can incur large communication costs, not to mention the privacy concerns of the private data. Therefore, we are aimed at finding distributed online algorithm for solving (2). In the following two sections, we propose two different distributed online algorithms with complementary merits in accuracy and computational complexity.

## Iii The Decentralized Online ADMM

In this section, we propose an alternating direction method of multipliers (ADMM) based decentralized online algorithm for solving (2). It is further simplified so that its iteration consists of simple closed-form computations and each node only needs to store and update one matrix and six dimensional vectors. We show that the gaps between the outputs of the proposed ADMM algorithm and the optimal points of (2) converge to zero. Before the development of the algorithm, we first present some rudimentary knowledge of ADMM in the following subsection.

ADMM is an optimization framework widely applied to various signal processing applications, including wireless communications [19], power systems [20] and multi-agent coordination [21]. It enjoys fast convergence speed under mild technical conditions [22] and is especially suitable for the development of distributed algorithms [23, 24]. ADMM solves problems of the following form:

 Minimizex,zf(x)+g(z)  s.t.  Ax+Bz=c, (3)

where are constants and are optimization variables. and are two convex functions. The augmented Lagrangian can be formed as:

 Lρ(x,z,y)=f(x)+g(z)+yT(Ax+Bz−c)+ρ2∥Ax+Bz−c∥22, (4)

where is the Lagrange multiplier and is some constant. The ADMM then iterates over the following three steps for (the iteration index):

 xk+1=argminxLρ(x,zk,yk), (5) zk+1=argminzLρ(xk+1,z,yk), (6) yk+1=yk+ρ(Axk+1+Bzk+1−c). (7)

The ADMM is guaranteed to converge to the optimal point of (3) as long as and are convex [23, 24]. It is recently shown that global linear convergence can be ensured provided additional assumptions on problem (3) holds [22].

### Iii-B Development of Decentralized Online ADMM for (2)

To apply ADMM to (2), we first transform it to the form of (3). We introduce auxiliary variables and , where denotes the cardinality of a set. Denote the index of -th neighbor of node as . Thus, problem (2) can be equivalently transformed into the following problem:

 Minimize  N∑n=1T∑t=1λT−t(dn(t)−un(t)Txn)2            +βN∑n=1⎡⎢⎣|Ωn|∥xn∥22−2⎛⎝|Ωn|∑i=1vn,i⎞⎠Txn+|Ωn|∑i=1∥vn,i∥22⎤⎥⎦+γN∑n=1∥wn∥1s.t.  xn=wn,n=1,...,N,      vn,i=wg(n,i),n=1...,N,i=1...,|Ωn|, (8)

where the optimization variables are . Note that optimization problem (8) is in the form of (3) (regarding ’s and ’s as the variable in (3) and ’s as the variable in (3)). Thus, we can apply ADMM to problem (8). Introducing Lagrange multiplier , we can form the augmented Lagrangian of (8) as follows:

 (9)

In the following, for ease of notation, we use to represent all the and similarly for . We apply the ADMM updates (5), (6) and (7) to problem (8) as follows:

 {xk+1,vk+1}=argminx,vLρ(x,v,wk,yk,zk), (10) wk+1=argminwLρ(xk+1,vk+1,w,yk,zk), (11) yk+1n=ykn+ρ(xk+1n−wk+1n), (12) zk+1n,i=zkn,i+ρ(vk+1n,i−wk+1g(n,i)). (13)

In the following, we detail how to implement the updates of the primal variables, i.e., (10) and (11), in a distributed and online fashion.

#### Iii-B1 Updating x and v

The update of and in (10) can be decomposed across nodes. For each node , the subproblem is:

 (14)

Define the data dependent input correlation matrix and input-output cross correlation vector of node at time to be:

 Rn(T)=T∑t=1λT−tun(T)un(T)T, (15) pn(T)=T∑t=1λT−tdn(t)un(t). (16)

Note that the objective function in (14) is a convex quadratic function. Hence, the necessary and sufficient condition for optimality of problem (14) is that the gradient of the objective function vanishes. The gradient of the objective function, which is denoted as , with respect to and can be computed as follows:

 ∇xnJkn(T)=(2Rn(T)+2β|Ωn|I+ρI)xn−2β|Ωn|∑i=1vn,i−2pn(T)+ykn−ρwkn, (17) ∇vn,iJkn(T)=−2βxn+(2β+ρ)vn,i+zkn,i−ρwkg(n,i). (18)

Letting the gradients with respect to and be zero, we rewrite the update in (14) as:

 ⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣2Rn(T)+2β|Ωn|I+ρI\tabularcell@hbox−2βI\tabularcell@hbox−2βI⋯−2βI[2pt/2pt]−2βI\tabularcell@hbox(2β+ρ)I\tabularcell@hbox−2βI\tabularcell@hbox\tabularcell@hbox(2β+ρ)I\huge 0⋮\tabularcell@hbox\tabularcell@hbox⋱−2βI\tabularcell@hbox\tabularcell@hbox\huge 0(2β+ρ)I⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣xk+1n[2pt/2pt]vk+1n,1vk+1n,2⋮vk+1n,|Ωn|⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣2pn(T)−ykn+ρwkn[2pt/2pt]−zkn,1+ρwkg(n,1)−zkn,2+ρwkg(n,2)⋮−zkn,|Ωn|+ρwkg(n,|Ωn|)⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦ (19)

To inverse the matrix in (19), we need to use the following matrix inversion lemma.

###### Lemma 1.

For arbitrary matrices such that all the matrix inversions at the R.H.S. of (20) exist, we have:

 (20)

Define a new matrix:

 Fn(T)=[2Rn(T)+(ρ+2βρ|Ωn|2β+ρ)I]−1. (21)

By invoking the matrix inversion lemma (20), we can solve for the update (19) in closed form:

 xk+1n=Fn(T)(2pn(T)−ykn+ρwkn)+2β2β+ρFn(T)|Ωn|∑i=1(−zkn,i+ρwkg(n,i)), (22) vk+1n,i=2β2β+ρFn(T)(2pn(T)−ykn+ρwkn)+12β+ρ(−zkn,i+ρwkg(n,i)) +(2β2β+ρ)2Fn(T)|Ωn|∑j=1(−zkn,j+ρwkg(n,j)). (23)

#### Iii-B2 Updating w

We note that the update for in (11) can be decomposed not only across nodes but also across each entry of the vector . For each node , the -th entry of can be updated as follows:

 wk+1n(l)=argminwn(l){γ|wn(l)|−ykn(l)wn(l)−⎛⎝∑g(m,i)=nzkm,i(l)⎞⎠wn(l)+ρ2[wn(l)−xk+1n(l)]2 +ρ2∑g(m,i)=n[wn(l)−vk+1m,i(l)]2} (24) =argminwn(l){γ|wn(l)|+ρ2(1+|Ωn|)[wn(l)−1ρ(1+|Ωn|)(ykn(l)+ρxk+1n(l) +∑g(m,i)=n(zkm,i(l)+vk+1m,i(l)))]2} (25) =Sγρ(1+|Ωn|)⎛⎝1ρ(1+|Ωn|)(ykn(l)+ρxk+1n(l)+∑g(m,i)=n(zkm,i(l)+vk+1m,i(l)))⎞⎠, (26)

where the soft-threshold function is defined for as follows:

 Sκ(a)=⎧⎨⎩a−κ,~{}~{}if~{}~{% }a>κ,0,~{}~{}if~{}~{}|a|≤κ,a+κ,~{}~{}if~{}~{}a<κ. (27)

In (26), we have made use of the following fact.

###### Lemma 2.

For any , we have:

 Sλρ(v)=argminx (λ|x|+ρ2(x−v)2). (28)

Once we extend the definition of to vectors in a entrywise way, we can write the update for compactly as:

 wk+1n=Sγρ(1+|Ωn|)⎛⎝1ρ(1+|Ωn|)(ykn+ρxk+1n+∑g(m,i)=n(zkm,i+vk+1m,i))⎞⎠. (29)

#### Iii-B3 Online Algorithm with Varying T

So far, the derived ADMM algorithm is only suitable for one particular time shot . Since it takes iterations for ADMM to converge to the optimal point, for each time , we ought to run multiple rounds of ADMM iterations for some sufficiently large . After the ADMM has converged for this particular time , we update the data related quantities ( and ) and move to the next time slot. However, since the underlying weight vectors are varying across time (i.e., the underlying linear system is non-stationary), it is meaningless to estimate the weight vectors very accurately for every time slot. Thus, in the following, we choose , i.e., only one iteration of ADMM update is executed in each time slot. This is inspired by many existing adaptive algorithms such as the LMS algorithm, where only one step of gradient descent is performed at each time slot [1]. As such, we replace with in the previously derived updates (22), (23), (29) and get updates that are suitable for varying time :

 xn(T)=Fn(T)(2pn(T)−yn(T−1)+ρwn(T−1)) +2β2β+ρFn(T)|Ωn|∑i=1(−zn,i(T−1)+ρwg(n,i)(T−1)), (30) vn,i(T)=2β2β+ρFn(T)(2pn(T)−yn(T−1)+ρwn(T−1)) +12β+ρ(−zn,i(T−1)+ρwg(n,i)(T−1)) +(2β2β+ρ)2Fn(T)|Ωn|∑j=1(−zn,j(T−1)+ρwg(n,j)(T−1)), (31) (32)

Moreover, the updates (12) and (13) for dual variables can be rewritten as:

 yn(T)=yn(T−1)+ρ(xn(T)−wn(T)), (33) zn,i(T)=zn,i(T−1)+ρ(vn,i(T)−wg(n,i)(T)). (34)

The correlation matrices and cross-correlation vectors can be updated as follows:

 Rn(T+1)=λRn(T)+un(T+1)un(T+1)T, (35) pn(T+1)=λpn(T)+dn(T+1)un(T+1). (36)

And is computed according to (21).

###### Remark 1.

The computation of in (21) necessitates inversion of an matrix, which incurs a computational complexity of unless special structure is present. For the special case of (which is suitable for time-invariant weight vectors), this burden can be alleviated as follows. According to (21), (35) and the condition that , we have:

 Fn(T) =[2(Rn(T−1)+un(T)un(T)T)+(ρ+2βρ|Ωn|2β+ρ)I]−1 (37) =[F−1n(T−1)+2un(T)un(T)T]−1 (38) =Fn(T−1)−Fn(T−1)un(T)un(T)TFn(T−1)12+un(T)TFn(T−1)un(T). (39)

However, in the general case where , the matrix inversion incurred by the computation of is inevitable, which is the most computationally intensive part of the proposed ADMM algorithm.

So far, the ADMM updates involve primal variables and dual variables . For each node , and include -dimensional vectors, which is costly to sustain in terms of communication and storage overhead, especially when the numbers of neighbors (degrees) are large. This motivates us to simplify the ADMM updates (30)-(36) so that the number of vectors at each node is independent of its degree. To this end, we first define the following auxiliary variables:

 z––n(T)=|Ωn|∑i=1zn,i(T), (40) ¯¯¯zn(T)=∑g(m,i)=nzm,i(T), (41) v––n(T)=|Ωn|∑i=1vn,i(T), (42) ¯¯¯vn(T)=∑g(m,i)=nvm,i(T), (43) ¯¯¯¯¯wn(T)=∑m∈Ωnwm(T), (44) ηn(T)=Fn(T)(2pn(T)−yn(T−1)+ρwn(T−1)), (45) θn(T)=Fn(T)(−z––n(T−1)+ρ¯¯¯¯¯wn(T−1)), (46) ¯¯¯ηn(T)=∑m∈Ωnηm(T), (47) ¯¯¯θn(T)=∑m∈Ωnθm(T). (48)

Thus, the update for in (30) can be rewritten as:

 xn(T)=ηn(T)+2β2β+ρθn(T). (49)

Using (31) yields the update for and :

 v––n(T) =2β|Ωn|2β+ρηn(T)+(2β2β+ρ)2|Ωn|θn(T)+12β+ρ(−z––n(T−1)+ρ¯¯¯¯¯wn(T−1)). (50)
 ¯¯¯vn(T) =2β2β+ρ¯¯¯ηn(T)+(2β2β+ρ)2¯¯¯θn(T)+12β+ρ(−¯¯¯zn(T−1)+ρ|Ωn|wn(T−1)). (51)

The update for can be rewritten as:

 wn(T)=Sγρ(1+|Ωn|)(1ρ(1+|Ωn|)(yn(T−1)+ρxn(T)+¯¯¯zn(T−1)+ρ¯¯¯vn(T))). (52)

Similarly, from (34), we can spell out the updates for and :

 z––n(T)=z––n(T−1)+ρ(v––n(T)−¯¯¯¯¯wn(T)), (53) ¯¯¯zn(T)=¯¯¯zn(T−1)+ρ(¯¯¯vn(T)−|Ωn|wn(T)). (54)

Now, we are ready to formally present the proposed decentralized online ADMM algorithm for solving (2), which is summarized in Algorithm 1. Notice that the algorithm is completely distributed: each node only needs to communicate with its neighbors. It is also online (real-time): each node only needs to store and update one matrix and six dimensional vectors . All other involved quantities in Algorithm 1 are intermediate and can be derived from these stored matrices and vectors.

### Iii-C Convergence of the Algorithm 1

In this subsection, we briefly discuss about the convergence of Algorithm 1 and show that the gap between its output, , and the optimal point of problem (2), which we denote as , converges to zero. We make the following assumptions.

###### Assumption 1.

The true weight vector is time-invariant, i.e., the linear regression data model is .

###### Assumption 2.

For each node , the input process is independent across time with time-invariant correlation matrix .

###### Assumption 3.

For each node , the noise process has zero mean, i.e., and is independent across time and independent from the input process .

Note that all of these assumptions are standard when analyzing the performance of adaptive algorithms in the literature [1]. From the definition of and , we know that they are weighted sum of i.i.d. terms. According to the strong law of large numbers for weighted sums [25, 10], as , converges to . Similarly, converges to . When , the optimization problem at time , i.e., problem (2), is to minimize (w.r.t. ):

 N∑n=1T∑t=1λT−t(wTnun(t)un(t)Twn−2dn(t)un(t)Twn)+βN∑n=1∑m∈Ωn∥wn−wm∥22+γN∑n=1∥wn∥1 ≈11−λN∑n=1(wTnRnwn−2˜wT% nRnwn)+βN∑n=1∑m∈Ωn∥wn−wm∥22+γN∑n=1∥wn∥1. (57)

Note that the R.H.S. of (57) does not depend on . Note that ADMM is guaranteed to converge to the optimal point for static convex optimization problem of the form (3) and the R.H.S. of (57) can be transformed into the form of (3) as we do in Subsection III-B. So, the output of Algorithm 1, , converges to the minimum point of the R.H.S. of (57). Due to (57) and the definition of , we know that also converges to the minimum point of the R.H.S. of (57). Hence, the difference between the output of Algorithm 1, i.e., , and the optimal point of (2), i.e., , converges to zero.

## Iv The Decentralized Online Subgradient Method

The implementation of the proposed Algorithm 1 necessitates an inversion of an matrix at each time and each node, which may not be suitable for nodes with low computational capability. In fact, a relatively high computational overhead is a general drawback of dual domain methods (e.g., ADMM) in optimization theory [17]. On the contrary, primal domain methods such as gradient descent method, though having relatively slow convergence speed, enjoys low computational complexity [26]. As such, in this section, we present a distributed online subgradient method for problem (2) to trade off convergence speed and accuracy for low computational complexity.

### Iv-a Development of the Decentralized Online Subgradient Method

Recall the optimization problem at time , i.e., problem (2). Denote the objective function of (2) as . We derive the subdifferential (the set of subgradients [27]) of at to be:

 ∂HT(w)=⎡⎢ ⎢ ⎢ ⎢ ⎢⎣2R1(T)w1−2p1(T)+2β(2|Ω1|w1−2∑m∈Ω1wm)+γsgn(w1)⋮2RN(T)wN−2pN(T)+2β(2|ΩN|wN−2∑m∈ΩNwm)+γsgn(wN)⎤⎥ ⎥ ⎥ ⎥ ⎥⎦, (58)

where the sign (set) function is defined as:

 sgn(x)=⎧⎨⎩1,  if  x>0,−1,  if  x<0,,  if  x=0. (59)

The extension of the function to vectors is entrywise. The subgradient method is to simply use the iteration , where is any subgradient of at and is the step size [27]. This naturally leads to the following decentralized online update:

 wn(T)=wn(T−1)−α[2Rn(T)wn(T−1)−2pn(T)+4β∑m∈Ωn(wn(T−1)−wm(T−1)) +γsgn(wn(T−1))], (60)

where is any number within the interval 111There is a standard abuse of notation for the function: in (58) and (59), is defined to be the interval while in (IV-A), is defined to be any arbitrary number within . In the following, the latter definition will be used.. By introducing an auxiliary variable