Decentralized Sparse Multitask RLS over Networks

TL;DR

This paper introduces decentralized algorithms for sparse multitask recursive least squares over networks, enabling efficient, real-time collaborative estimation of multiple unknown vectors with theoretical convergence analysis.

Contribution

It proposes novel decentralized online ADMM and subgradient algorithms for sparse multitask RLS, with simplified implementation and convergence guarantees.

Findings

Algorithms effectively estimate multiple unknown vectors in networked systems.

The subgradient method has proven convergence with an explicit error bound.

Numerical simulations demonstrate the algorithms' accuracy and low complexity.

Abstract

Distributed adaptive signal processing has attracted much attention in the recent decade owing to its effectiveness in many decentralized real-time applications in networked systems. Because many natural signals are highly sparse with most entries equal to zero, several decentralized sparse adaptive algorithms have been proposed recently. Most of them is focused on the single task estimation problems, in which all nodes receive data associated with the same unknown vector and collaborate to estimate it. However, many applications are inherently multitask oriented and each node has its own unknown vector different from others'. The related multitask estimation problem benefits from collaborations among the nodes as neighbor nodes usually share analogous properties and thus similar unknown vectors. In this work, we study the distributed sparse multitask recursive least squares (RLS) problem over networks. We first propose a decentralized online alternating direction method of multipliers (ADMM) algorithm for the formulated RLS problem. The algorithm is simplified for easy implementation with closed-form computations in each iteration and low storage requirements. Moreover, to further reduce the complexity, we present a decentralized online subgradient method with low computational overhead. We theoretically analyze the convergence behavior of the proposed subgradient method and derive an error bound related to the network topology and algorithm parameters. The effectiveness of the proposed algorithms is corroborated by numerical simulations and an accuracy-complexity tradeoff between the proposed two algorithms is highlighted.

Equations

$$d_n(t)={\mathbf{u}}_n^{\mathsf{T}}(t){\tilde{\mathbf{w}}}_n(t)+e_n(t),$$

$${\text{Minimize}}_{{\mathbf{w}}_1,...,{\mathbf{w}}_N}\sum _{n=1}^N\sum _{t=1}^T{\lambda }^{T-t}{\left(d_n(t)-{\mathbf{u}}_n^{\mathsf{T}}(t){\mathbf{w}}_n\right)}^2+\beta \sum _{n=1}^N\sum _{m\in {\Omega }_n}\Vert {\mathbf{w}}_n-{\mathbf{w}}_m{\Vert }_2^2+\gamma \sum _{n=1}^N\Vert {\mathbf{w}}_n{\Vert }_1,$$

$${\text{Minimize}}_{\mathbf{x},\mathbf{z}}f(\mathbf{x})+g(\mathbf{z})\text{s.t.}\mathbf{Ax}+\mathbf{Bz}=\mathbf{c},$$

$${\mathfrak{L}}_{\rho }(\mathbf{x},\mathbf{z},\mathbf{y})=f(\mathbf{x})+g(\mathbf{z})+{\mathbf{y}}^{\mathsf{T}}(\mathbf{Ax}+\mathbf{Bz}-\mathbf{c})+\frac{\rho }{2}\Vert \mathbf{Ax}+\mathbf{Bz}-\mathbf{c}{\Vert }_2^2,$$

$${\mathbf{x}}^{k+1}=\arg \underset{\mathbf{x}}{\min }{\mathfrak{L}}_{\rho }\left(\mathbf{x},{\mathbf{z}}^k,{\mathbf{y}}^k\right),$$

$${\mathbf{z}}^{k+1}=\arg \underset{\mathbf{z}}{\min }{\mathfrak{L}}_{\rho }\left({\mathbf{x}}^{k+1},\mathbf{z},{\mathbf{y}}^k\right),$$

$${\mathbf{y}}^{k+1}={\mathbf{y}}^k+\rho \left({\mathbf{Ax}}^{k+1}+{\mathbf{Bz}}^{k+1}-\mathbf{c}\right).$$

$$\begin{array}{rl}& \text{Minimize}\sum _{n=1}^N\sum _{t=1}^T{\lambda }^{T-t}{\left(d_n(t)-{\mathbf{u}}_n(t{)}^{\mathsf{T}}{\mathbf{x}}_n\right)}^2\\ & +\beta \sum _{n=1}^N\left[|{\Omega }_n|\Vert {\mathbf{x}}_n{\Vert }_2^2-2{\left(\sum _{i=1}^{|{\Omega }_n|}{\mathbf{v}}_{n,i}\right)}^{\mathsf{T}}{\mathbf{x}}_n+\sum _{i=1}^{|{\Omega }_n|}\Vert {\mathbf{v}}_{n,i}{\Vert }_2^2\right]+\gamma \sum _{n=1}^N\Vert {\mathbf{w}}_n{\Vert }_1\\ & \text{s.t.}{\mathbf{x}}_n={\mathbf{w}}_n,n=1,...,N,\\ & {\mathbf{v}}_{n,i}={\mathbf{w}}_{g(n,i)},n=\mathrm{1...},N,i=\mathrm{1...},|{\Omega }_n|,\end{array}$$

$$\begin{array}{rl}& {\mathfrak{L}}_{\rho }(\{{\mathbf{x}}_n,{\mathbf{v}}_{n,i},{\mathbf{w}}_n,{\mathbf{y}}_n,{\mathbf{z}}_{n,i}{\}}_{n=1,...,N,i=1,...,|{\Omega }_n|})\\ & =\sum _{n=1}^N\sum _{t=1}^T{\lambda }^{T-t}{\left(d_n(t)-{\mathbf{u}}_n(t{)}^{\mathsf{T}}{\mathbf{x}}_n\right)}^2+\beta \sum _{n=1}^N\left[|{\mathcal{N}}_n|\Vert {\mathbf{x}}_n{\Vert }_2^2-2{\left(\sum _{i=1}^{|{\Omega }_n|}{\mathbf{v}}_{n,i}\right)}^{\mathsf{T}}{\mathbf{x}}_n+\sum _{i=1}^{|{\Omega }_n|}\Vert {\mathbf{v}}_{n,i}{\Vert }_2^2\right]\\ & +\gamma \sum _{n=1}^N\Vert {\mathbf{w}}_n{\Vert }_1+\sum _{n=1}^N{\mathbf{y}}_n^{\mathsf{T}}({\mathbf{x}}_n-{\mathbf{w}}_n)+\sum _{n=1}^N\sum _{i=1}^{|{\Omega }_n|}{\mathbf{z}}_{n,i}^{\mathsf{T}}({\mathbf{v}}_{n,i}-{\mathbf{w}}_{g_{n,i}})\\ & +\frac{\rho }{2}\sum _{n=1}^N\Vert {\mathbf{x}}_n-{\mathbf{w}}_n{\Vert }_2^2+\frac{\rho }{2}\sum _{n=1}^N\sum _{i=1}^{|{\Omega }_n|}\Vert {\mathbf{v}}_{n,i}-{\mathbf{w}}_{g(n,i)}{\Vert }_2^2\end{array}$$

$$\left\{{\mathbf{x}}^{k+1},{\mathbf{v}}^{k+1}\right\}=\arg \underset{\mathbf{x},\mathbf{v}}{\min }{\mathfrak{L}}_{\rho }\left(\mathbf{x},\mathbf{v},{\mathbf{w}}^k,{\mathbf{y}}^k,{\mathbf{z}}^k\right),$$

$${\mathbf{w}}^{k+1}=\arg \underset{\mathbf{w}}{\min }{\mathfrak{L}}_{\rho }\left({\mathbf{x}}^{k+1},{\mathbf{v}}^{k+1},\mathbf{w},{\mathbf{y}}^k,{\mathbf{z}}^k\right),$$

$${\mathbf{y}}_n^{k+1}={\mathbf{y}}_n^k+\rho \left({\mathbf{x}}_n^{k+1}-{\mathbf{w}}_n^{k+1}\right),$$

$${\mathbf{z}}_{n,i}^{k+1}={\mathbf{z}}_{n,i}^k+\rho \left({\mathbf{v}}_{n,i}^{k+1}-{\mathbf{w}}_{g(n,i)}^{k+1}\right).$$

$$\displaystyle\begin{split}&\left\{\mathbf{x}_{n}^{k+1},\left\{\mathbf{v}_{n,i}^{k+1}\right\}_{i=1,...,|\Omega_{n}|}\right\}\\ &=\arg\min_{\mathbf{x}_{n},\left\{\mathbf{v}_{n,i}\right\}_{i=1,...,|\Omega_{n}|}}~{}\Bigg{\{}\sum_{t=1}^{T}\lambda^{T-t}\left(d_{n}(t)-\mathbf{u}_{n}(t)^{\mathsf{T}}\mathbf{x}_{n}\right)^{2}\\ &+\beta\left[|\Omega_{n}|\|\mathbf{x}_{n}\|_{2}^{2}-2\left(\sum_{i=1}^{|\Omega_{n}|}\mathbf{v}_{n,i}\right)^{\mathsf{T}}\mathbf{x}_{n}+\sum_{i=1}^{|\Omega_{n}|}\|\mathbf{v}_{n,i}\|_{2}^{2}\right]\\ &+\mathbf{y}_{n}^{k\mathsf{T}}\mathbf{x}_{n}+\sum_{i=1}^{|\Omega_{n}|}\mathbf{z}_{n,i}^{k\mathsf{T}}\mathbf{v}_{n,i}+\frac{\rho}{2}\left\|\mathbf{x}_{n}-\mathbf{w}_{n}^{k}\right\|_{2}^{2}+\frac{\rho}{2}\sum_{i=1}^{|\Omega_{n}|}\left\|\mathbf{v}_{n,i}-\mathbf{w}_{g(n,i)}^{k}\right\|_{2}^{2}\Bigg{\}}.\end{split}$$

$${\mathbf{R}}_n(T)=\sum _{t=1}^T{\lambda }^{T-t}{\mathbf{u}}_n(T){\mathbf{u}}_n(T{)}^{\mathsf{T}},$$

$${\mathbf{p}}_n(T)=\sum _{t=1}^T{\lambda }^{T-t}{d}_n(t){\mathbf{u}}_n(t).$$

$${\nabla }_{{\mathbf{x}}_n}{J}_n^k(T)=(2{\mathbf{R}}_{\mathbf{n}}(\mathbf{T})+2\beta |{\Omega }_n|\mathbf{I}+\rho \mathbf{I}){\mathbf{x}}_n-2\beta \sum _{i=1}^{|{\Omega }_n|}{\mathbf{v}}_{n,i}-2{\mathbf{p}}_n(T)+{\mathbf{y}}_n^k-\rho {\mathbf{w}}_n^k,$$

$${\nabla }_{{\mathbf{v}}_{n,i}}{J}_n^k(T)=-2\beta {\mathbf{x}}_n+(2\beta +\rho ){\mathbf{v}}_{n,i}+{\mathbf{z}}_{n,i}^k-\rho {\mathbf{w}}_{g(n,i)}^k.$$

$$\displaystyle\begin{split}&\left[\begin{array}[]{c;{2pt/2pt}cccc}2\mathbf{R}_{n}(T)+2\beta|\Omega_{n}|\mathbf{I}+\rho\mathbf{I}&-2\beta\mathbf{I}&-2\beta\mathbf{I}&\cdots&-2\beta\mathbf{I}\\ \hdashline[2pt/2pt]-2\beta\mathbf{I}&(2\beta+\rho)\mathbf{I}&&&\\ -2\beta\mathbf{I}&&(2\beta+\rho)\mathbf{I}&&\text{\huge 0}\\ \vdots&&&\ddots&\\ -2\beta\mathbf{I}&&\text{\huge 0}&&(2\beta+\rho)\mathbf{I}\end{array}\right]\left[\begin{array}[]{c}\mathbf{x}_{n}^{k+1}\\ \hdashline[2pt/2pt]\mathbf{v}_{n,1}^{k+1}\\ \mathbf{v}_{n,2}^{k+1}\\ \vdots\\ \mathbf{v}_{n,|\Omega_{n}|}^{k+1}\end{array}\right]\\ &=\left[\begin{array}[]{c}2\mathbf{p}_{n}(T)-\mathbf{y}_{n}^{k}+\rho\mathbf{w}_{n}^{k}\\ \hdashline[2pt/2pt]-\mathbf{z}_{n,1}^{k}+\rho\mathbf{w}_{g(n,1)}^{k}\\ -\mathbf{z}_{n,2}^{k}+\rho\mathbf{w}_{g(n,2)}^{k}\\ \vdots\\ -\mathbf{z}_{n,|\Omega_{n}|}^{k}+\rho\mathbf{w}_{g(n,|\Omega_{n}|)}^{k}\\ \end{array}\right]\end{split}$$

$$\displaystyle\left[\begin{array}[]{cc}\mathbf{A}&\mathbf{B}\\ \mathbf{C}&\mathbf{D}\end{array}\right]^{-1}=\left[\begin{array}[]{cc}\left(\mathbf{A}-\mathbf{B}\mathbf{D}^{-1}\mathbf{C}\right)^{-1}&-\left(\mathbf{A}-\mathbf{B}\mathbf{D}^{-1}\mathbf{C}\right)^{-1}\mathbf{BD}^{-1}\\ -\mathbf{D}^{-1}\mathbf{C}\left(\mathbf{A}-\mathbf{B}\mathbf{D}^{-1}\mathbf{C}\right)^{-1}&\mathbf{D}^{-1}+\mathbf{D}^{-1}\mathbf{C}\left(\mathbf{A}-\mathbf{B}\mathbf{D}^{-1}\mathbf{C}\right)^{-1}\mathbf{BD}^{-1}.\end{array}\right]$$

$${\mathbf{F}}_n(T)={\left[2{\mathbf{R}}_n(T)+\left(\rho +\frac{2\beta \rho |{\Omega }_n|}{2\beta +\rho }\right)\mathbf{I}\right]}^{-1}.$$

$${\mathbf{x}}_n^{k+1}={\mathbf{F}}_n(T)\left(2{\mathbf{p}}_n(T)-{\mathbf{y}}_n^k+\rho {\mathbf{w}}_n^k\right)+\frac{2\beta }{2\beta +\rho }{\mathbf{F}}_n(T)\sum _{i=1}^{|{\Omega }_n|}\left(-{\mathbf{z}}_{n,i}^k+\rho {\mathbf{w}}_{g(n,i)}^k\right),$$

$${\mathbf{v}}_{n,i}^{k+1}=\frac{2\beta }{2\beta +\rho }{\mathbf{F}}_n(T)\left(2{\mathbf{p}}_n(T)-{\mathbf{y}}_n^k+\rho {\mathbf{w}}_n^k\right)+\frac{1}{2\beta +\rho }\left(-{\mathbf{z}}_{n,i}^k+\rho {\mathbf{w}}_{g(n,i)}^k\right)$$

$$+{\left(\frac{2\beta }{2\beta +\rho }\right)}^2{\mathbf{F}}_n(T)\sum _{j=1}^{|{\Omega }_n|}\left(-{\mathbf{z}}_{n,j}^k+\rho {\mathbf{w}}_{g(n,j)}^k\right).$$

$$\displaystyle w_{n}^{k+1}(l)=\arg\min_{w_{n}(l)}\Bigg{\{}\gamma|w_{n}(l)|-y_{n}^{k}(l)w_{n}(l)-\left(\sum_{g(m,i)=n}z_{m,i}^{k}(l)\right)w_{n}(l)+\frac{\rho}{2}\left[w_{n}(l)-x_{n}^{k+1}(l)\right]^{2}$$

$$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\frac{\rho}{2}\sum_{g(m,i)=n}\left[w_{n}(l)-v_{m,i}^{k+1}(l)\right]^{2}\Bigg{\}}$$

$$\displaystyle=\arg\min_{w_{n}(l)}\Bigg{\{}\gamma|w_{n}(l)|+\frac{\rho}{2}(1+|\Omega_{n}|)\Bigg{[}w_{n}(l)-\frac{1}{\rho(1+|\Omega_{n}|)}\Bigg{(}y_{n}^{k}(l)+\rho x_{n}^{k+1}(l)$$

$$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\sum_{g(m,i)=n}\left(z_{m,i}^{k}(l)+v_{m,i}^{k+1}(l)\right)\Bigg{)}\Bigg{]}^{2}\Bigg{\}}$$

$$\displaystyle=\mathcal{S}_{\frac{\gamma}{\rho(1+|\Omega_{n}|)}}\left(\frac{1}{\rho(1+|\Omega_{n}|)}\Bigg{(}y_{n}^{k}(l)+\rho x_{n}^{k+1}(l)+\sum_{g(m,i)=n}\left(z_{m,i}^{k}(l)+v_{m,i}^{k+1}(l)\right)\Bigg{)}\right),$$

$${\mathcal{S}}_{\kappa }(a)=\{\begin{array}{l}a-\kappa ,\text{ if }a>\kappa ,\\ 0,\text{ if }|a|\le \kappa ,\\ a+\kappa ,\text{ if }a<\kappa .\end{array}$$

$${\mathcal{S}}_{\frac{\lambda }{\rho }}(v)=\arg \underset{x}{\min }\left(\lambda |x|+\frac{\rho }{2}(x-v{)}^2\right).$$

$$\displaystyle\mathbf{w}_{n}^{k+1}=\mathcal{S}_{\frac{\gamma}{\rho(1+|\Omega_{n}|)}}\left(\frac{1}{\rho(1+|\Omega_{n}|)}\Bigg{(}\mathbf{y}_{n}^{k}+\rho\mathbf{x}_{n}^{k+1}+\sum_{g(m,i)=n}\left(\mathbf{z}_{m,i}^{k}+\mathbf{v}_{m,i}^{k+1}\right)\Bigg{)}\right).$$

$${\mathbf{x}}_n(T)={\mathbf{F}}_n(T)\left(2{\mathbf{p}}_n(T)-{\mathbf{y}}_n(T-1)+\rho {\mathbf{w}}_n(T-1)\right)$$

$$+\frac{2\beta }{2\beta +\rho }{\mathbf{F}}_n(T)\sum _{i=1}^{|{\Omega }_n|}\left(-{\mathbf{z}}_{n,i}(T-1)+\rho {\mathbf{w}}_{g(n,i)}(T-1)\right),$$

$${\mathbf{v}}_{n,i}(T)=\frac{2\beta }{2\beta +\rho }{\mathbf{F}}_n(T)\left(2{\mathbf{p}}_n(T)-{\mathbf{y}}_n(T-1)+\rho {\mathbf{w}}_n(T-1)\right)$$

$$+\frac{1}{2\beta +\rho }\left(-{\mathbf{z}}_{n,i}(T-1)+\rho {\mathbf{w}}_{g(n,i)}(T-1)\right)$$

$$+{\left(\frac{2\beta }{2\beta +\rho }\right)}^2{\mathbf{F}}_n(T)\sum _{j=1}^{|{\Omega }_n|}\left(-{\mathbf{z}}_{n,j}(T-1)+\rho {\mathbf{w}}_{g(n,j)}(T-1)\right),$$

$$\displaystyle\mathbf{w}_{n}(T)=\mathcal{S}_{\frac{\gamma}{\rho(1+|\Omega_{n}|)}}\left(\frac{1}{\rho(1+|\Omega_{n}|)}\Bigg{(}\mathbf{y}_{n}(T-1)+\rho\mathbf{x}_{n}(T)+\sum_{g(m,i)=n}\left(\mathbf{z}_{m,i}(T-1)+\mathbf{v}_{m,i}(T)\right)\Bigg{)}\right).$$

$${\mathbf{y}}_n(T)={\mathbf{y}}_n(T-1)+\rho \left({\mathbf{x}}_n(T)-{\mathbf{w}}_n(T)\right),$$

Full text

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 $l_{1}$ 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 $M\times M$ matrix and six $M$ dimensional vectors, where $M$ 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 $N$ 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 $n$ (those who are linked with node $n$ by an edge) as $\Omega_{n}$. 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 $t=1,2,...$. Each node $n$ has an unknown (slowly) time-variant $M$ dimensional weight vector $\widetilde{\mathbf{w}}_{n}(t)\in\mathbb{R}^{M}$ 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 $n$ has access to a sequence of private measurements $\{d_{n}(t),\mathbf{u}_{n}(t)\}_{t=1,2,...}$, where $\mathbf{u}_{n}(t)\in\mathbb{R}^{M}$ is the input regressor at time $t$ and $d_{n}(t)\in\mathbb{R}$ is the output observation at time $t$. The measurement data are private in the sense that node $n$ has access only to its own measurement sequence but not others’. The data at node $n$ are assumed to conform to a linear regression model with (slowly) time-variant weight vector $\widetilde{\mathbf{w}}_{n}(t)$:

[TABLE]

where $e_{n}(t)$ is the output measurement noise at time $t$. 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 $\widetilde{\mathbf{w}}_{n}(t)$ 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 $\{\widetilde{\mathbf{w}}_{n}(t)\}_{n=1,2,...,N}$ while enforcing similarity between neighbors’ weight vectors and sparsity of all weight vectors. The estimator at time $T$ is the optimal solution of the following optimization problem:

[TABLE]

where $0<\lambda<1,\beta>0,\gamma>0$ are the forgetting factor of the RLS algorithm, regularization coefficient for similarity between neighbors’ weight vectors and regularization coefficient for sparsity, respectively. If $\beta=\infty$, 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 $\{\mathbf{u}_{n}(t),d_{n}(t)\}$ 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 $M\times M$ matrix and six $M$ 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.

III-A Preliminaries of ADMM

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:

[TABLE]

where $\mathbf{A}\in\mathbb{R}^{p\times n},B\in\mathbb{R}^{p\times m},c\in\mathbb{R}^{p}$ are constants and $\mathbf{x}\in\mathbb{R}^{n},\mathbf{z}\in\mathbb{R}^{m}$ are optimization variables. $f:\mathbb{R}^{n}\mapsto\mathbb{R}$ and $g:\mathbb{R}^{m}\mapsto\mathbb{R}$ are two convex functions. The augmented Lagrangian can be formed as:

[TABLE]

where $\mathbf{y}\in\mathbb{R}^{p}$ is the Lagrange multiplier and $\rho>0$ is some constant. The ADMM then iterates over the following three steps for $k\geq 0$ (the iteration index):

[TABLE]

The ADMM is guaranteed to converge to the optimal point of (3) as long as $f$ and $g$ 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 $\mathbf{x}_{n}\in\mathbb{R}^{M},n=1,...,N,$ and $\mathbf{v}_{n,i}\in\mathbb{R}^{M},n=1,...N,i=1,...,|\Omega_{n}|$, where $|\cdot|$ denotes the cardinality of a set. Denote the index of $i$-th neighbor of node $n$ as $g(n,i)$. Thus, problem (2) can be equivalently transformed into the following problem:

[TABLE]

where the optimization variables are $\mathbf{w}_{n},\mathbf{x}_{n},\mathbf{v}_{n,i},n=1,...,N,i=1,...,|\Omega_{n}|$. Note that optimization problem (LABEL:admm) is in the form of (3) (regarding $\mathbf{x}_{n}$’s and $\mathbf{v}_{n,i}$’s as the variable $\mathbf{x}$ in (3) and $\mathbf{w}_{n}$’s as the variable $\mathbf{z}$ in (3)). Thus, we can apply ADMM to problem (LABEL:admm). Introducing Lagrange multiplier $\mathbf{y}_{n}\in\mathbb{R}^{M},\mathbf{z}_{n,i}\in\mathbb{R}^{M}$, we can form the augmented Lagrangian of (LABEL:admm) as follows:

[TABLE]

In the following, for ease of notation, we use $\mathbf{x}$ to represent all the $\{\mathbf{x}_{n}\}$ and similarly for $\mathbf{v,w,y,z}$. We apply the ADMM updates (5), (6) and (7) to problem (LABEL:admm) as follows:

[TABLE]

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 $\mathbf{x}$ and $\mathbf{v}$

The update of $\mathbf{x}$ and $\mathbf{v}$ in (10) can be decomposed across nodes. For each node $n$, the subproblem is:

[TABLE]

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

[TABLE]

Note that the objective function in (LABEL:xv1) is a convex quadratic function. Hence, the necessary and sufficient condition for optimality of problem (LABEL:xv1) is that the gradient of the objective function vanishes. The gradient of the objective function, which is denoted as $J_{n}^{k}(T)$, with respect to $\mathbf{x}_{n}$ and $\mathbf{v}_{n,i}$ can be computed as follows:

[TABLE]

Letting the gradients with respect to $\mathbf{x}_{n}$ and $\mathbf{v}_{n,i}$ be zero, we rewrite the update in (LABEL:xv1) as:

[TABLE]

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

Lemma 1**.**

For arbitrary matrices $\mathbf{A}\in\mathbb{R}^{m\times m},\mathbf{B}\in\mathbb{R}^{m\times n},\mathbf{C}\in\mathbb{R}^{n\times m},\mathbf{D}\in\mathbb{R}^{n\times n}$ such that all the matrix inversions at the R.H.S. of (24) exist, we have:

[TABLE]

Define a new matrix:

[TABLE]

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

[TABLE]

III-B2 Updating $\mathbf{w}$

We note that the update for $\mathbf{w}$ in (11) can be decomposed not only across nodes but also across each entry of the vector $\mathbf{w}_{n}$. For each node $n$, the $l$-th entry of $\mathbf{w}_{n}$ can be updated as follows:

[TABLE]

where the soft-threshold function $\mathcal{S}$ is defined for $a\in\mathbb{R},\kappa>0$ as follows:

[TABLE]

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

Lemma 2**.**

For any $\lambda>0,\rho>0,v\in\mathbb{R}$, we have:

[TABLE]

Once we extend the definition of $\mathcal{S}$ to vectors in a entrywise way, we can write the update for $\mathbf{w}_{n}$ compactly as:

[TABLE]

III-B3 Online Algorithm with Varying $T$

So far, the derived ADMM algorithm is only suitable for one particular time shot $T$. Since it takes iterations for ADMM to converge to the optimal point, for each time $T$, we ought to run multiple rounds of ADMM iterations $k=1,...,K$ for some sufficiently large $K$. After the ADMM has converged for this particular time $T$, we update the data related quantities ($\mathbf{R}_{n}(T)$ and $\mathbf{p}_{n}(T)$) 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 $K=1$, 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 $k$ with $T-1$ in the previously derived updates (26), (27), (33) and get updates that are suitable for varying time $T$:

[TABLE]

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

[TABLE]

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

[TABLE]

And $\mathbf{F}_{n}(T)$ is computed according to (25).

Remark 1**.**

The computation of $\mathbf{F}_{n}(T)$ in (25) necessitates inversion of an $M\times M$ matrix, which incurs a computational complexity of $\mathcal{O}(M^{3})$ unless special structure is present. For the special case of $\lambda=1$ (which is suitable for time-invariant weight vectors), this burden can be alleviated as follows. According to (25), (39) and the condition that $\lambda=1$, we have:

[TABLE]

However, in the general case where $\lambda<1$, the matrix inversion incurred by the computation of $\mathbf{F}_{n}(T)$ is inevitable, which is the most computationally intensive part of the proposed ADMM algorithm.

III-B4 Simplification of the ADMM Updates

So far, the ADMM updates involve primal variables $\{\mathbf{v}_{n,i}\}$ and dual variables $\{\mathbf{z}\}_{n,i}$. For each node $n$, $\{\mathbf{v}_{n,i}\}$ and $\{\mathbf{z}\}_{n,i}$ include $2|\Omega_{n}|$ $M$-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 (34)-(40) so that the number of vectors at each node is independent of its degree. To this end, we first define the following auxiliary variables:

[TABLE]

Thus, the update for $\mathbf{x}$ in (34) can be rewritten as:

[TABLE]

Using (35) yields the update for $\mathbf{\underline{v}}_{n}(T)$ and $\mathbf{\overline{v}}_{n}(T)$:

[TABLE]

[TABLE]

The update for $\mathbf{w}_{n}(T)$ can be rewritten as:

[TABLE]

Similarly, from (38), we can spell out the updates for $\mathbf{\underline{z}}_{n}(T)$ and $\mathbf{\overline{z}}_{n}(T)$:

[TABLE]

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 $M\times M$ matrix $\mathbf{R}_{n}(T)$ and six $M$ dimensional vectors $\mathbf{p}_{n}(T),\mathbf{w}_{n}(T),\mathbf{\overline{w}}_{n}(T),\mathbf{y}_{n}(T),\mathbf{\underline{z}}_{n}(T),\mathbf{\overline{z}}_{n}(T)$. 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, $\mathbf{w}_{n}(T)$, and the optimal point of problem (2), which we denote as $\mathbf{w}_{n}^{*}(T)$, converges to zero. We make the following assumptions.

Assumption 1**.**

The true weight vector $\mathbf{\widetilde{w}}_{n}$ is time-invariant, i.e., the linear regression data model is $d_{n}(t)=\mathbf{u}_{n}(t)^{\mathsf{T}}\mathbf{\widetilde{w}}_{n}+e_{n}(t)$.

Assumption 2**.**

For each node $n$, the input process $\{\mathbf{u}_{n}(t)\}_{t=1,2,...}$ is independent across time with time-invariant correlation matrix $\mathbf{R}_{n}=\mathbb{E}\left[\mathbf{u}_{n}(t)\mathbf{u}_{n}(t)^{\mathsf{T}}\right]$.

Assumption 3**.**

For each node $n$, the noise process $\{e_{n}(t)\}_{t=1,2,...}$ has zero mean, i.e., $\mathbb{E}[e_{n}(t)]=0$ and is independent across time and independent from the input process $\{\mathbf{u}_{n}(t)\}$.

Note that all of these assumptions are standard when analyzing the performance of adaptive algorithms in the literature [1]. From the definition of $\mathbf{R}_{n}(T)$ and $\mathbf{p}_{n}(T)$, 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 $n\rightarrow\infty$, $\mathbf{R}_{n}(T)$ converges to $\lim_{T\rightarrow\infty}\mathbb{E}[\mathbf{R}_{n}(T)]=\frac{\mathbf{R}_{n}}{1-\lambda}$. Similarly, $\mathbf{p}_{n}(T)$ converges to $\lim_{T\rightarrow\infty}\mathbb{E}[\mathbf{p}_{n}(T)]=\frac{\mathbf{R}_{n}\mathbf{\widetilde{w}}_{n}}{1-\lambda}$. When $T\rightarrow\infty$, the optimization problem at time $T$, i.e., problem (2), is to minimize (w.r.t. $\{\mathbf{w}_{n}\}$):

[TABLE]

Note that the R.H.S. of (61) does not depend on $T$. 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 (61) can be transformed into the form of (3) as we do in Subsection III-B. So, the output of Algorithm 1, $\mathbf{w}_{n}(t)$, converges to the minimum point of the R.H.S. of (61). Due to (61) and the definition of $\mathbf{w}_{n}^{*}(T)$, we know that $\mathbf{w}_{n}^{*}(T)$ also converges to the minimum point of the R.H.S. of (61). Hence, the difference between the output of Algorithm 1, i.e., $\mathbf{w}_{n}(t)$, and the optimal point of (2), i.e., $\mathbf{w}_{n}^{*}(T)$, converges to zero.

IV The Decentralized Online Subgradient Method

The implementation of the proposed Algorithm 1 necessitates an inversion of an $M\times M$ 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 $T$, i.e., problem (2). Denote the objective function of (2) as $H_{T}(\mathbf{w})$. We derive the subdifferential (the set of subgradients [27]) of $H_{T}$ at $\mathbf{w}$ to be:

[TABLE]

where the sign (set) function is defined as:

[TABLE]

The extension of the $\operatorname{sgn}$ function to vectors is entrywise. The subgradient method is to simply use the iteration $\mathbf{w}(T)=\mathbf{w}(T-1)-\alpha\mathbf{g}$, where $\mathbf{g}\in H_{T}(\mathbf{w}(T-1))$ is any subgradient of $H_{T}$ at $\mathbf{w}(T-1)$ and $\alpha>0$ is the step size [27]. This naturally leads to the following decentralized online update:

[TABLE]

where $\operatorname{sgn}(0)$ is any number within the interval $[-1,1]$111There is a standard abuse of notation for the $\operatorname{sgn}$ function: in (65) and (66), $\operatorname{sgn}(0)$ is defined to be the interval $[-1,1]$ while in (IV-A), $\operatorname{sgn}(0)$ is defined to be any arbitrary number within $[-1,1]$. In the following, the latter definition will be used.. By introducing an auxiliary variable $\mathbf{\overline{w}}_{n}(T)$, we propose the decentralized online subgradient method for (2) in Algorithm 2. We observe that Algorithm 2 is completely decentralized as every node only communicates with its neighbors. It is also online since each node only needs to store and update one $M\times M$ matrix and three $M$ dimensional vectors. More importantly, Algorithm 2 is free of any matrix inversion, which is a major burden of Algorithm 1.

IV-B Convergence Analysis of Algorithm 2

In this subsection, we analyze the convergence behavior of Algorithm 2. We still make the same assumptions as in the analysis of the proposed ADMM algorithm in Subsection III-C, i.e., Assumptions 1-3. As explained in Subsection III-C, for large $T$, we can approximate $\mathbf{R}_{n}(T)$ by $\frac{\mathbf{R}_{n}}{1-\lambda}$ and $\mathbf{p}_{n}(T)$ by $\frac{\mathbf{R}_{n}\mathbf{\widetilde{w}}_{n}}{1-\lambda}$. Define the error vector of node $n$ at time $T$ to be:

[TABLE]

Thus, substituting (72) into the definition of $\mathbf{\overline{w}}_{n}(T-1)$ yields:

[TABLE]

Hence, using (70), (72) and the approximation of $\mathbf{R}_{n}(T),\mathbf{p}_{n}(T)$, we can derive a recursive equation for the error vector:

[TABLE]

Taking expectations yields:

[TABLE]

where $\mathbf{r}_{n}(T)$ is defined as:

[TABLE]

Denote the adjacency matrix of the network as $\mathbf{A}\in\mathbb{R}^{M\times M}$, i.e., $A(i,j)=1$ if $i,j$ are connected with an edge; otherwise $A(i,j)=0$. Define $\mathbf{C}=\mathbf{A}\otimes\mathbf{I}_{M\times M}\in\mathbb{R}^{MN\times MN}$ to be the block adjacency matrix, where $\otimes$ means Kronecker product. Define the matrix $\mathbf{B}\in\mathbb{R}^{MN\times MN}$ as follows:

[TABLE]

Furthermore, stack $\mathbf{f}_{n}(T)$ and $\mathbf{r}_{n}(T)$ into long vectors respectively:

[TABLE]

Hence, (77) can be written in a more compact form:

[TABLE]

Define $\boldsymbol{\Phi}=\frac{1}{\alpha}(\mathbf{I}-\mathbf{B})$. Since $\boldsymbol{\Phi}-4\beta\mathbf{C}$ is symmetric, we can perform eigendecomposition for it, i.e., $\boldsymbol{\Phi}-4\beta\mathbf{C}=\mathbf{Q\Lambda Q}^{\mathsf{T}}$ for some orthogonal matrix $\mathbf{Q}\in\mathbb{R}^{MN\times MN}$ and diagonal matrix $\boldsymbol{\Lambda}=\operatorname{diag}(\delta_{1},...,\delta_{MN})$, where $\delta_{i}$’s are the eigenvalues of $\boldsymbol{\Phi}-4\beta\mathbf{C}$. Therefore, $\mathbf{B}+4\alpha\beta\mathbf{C}=\mathbf{I}-\alpha(\mathbf{\Phi}-4\beta\mathbf{C})=\mathbf{Q\Sigma Q}^{\mathsf{T}}$, where we define $\mathbf{\Sigma}=\operatorname{diag}(1-\alpha\delta_{1},...,1-\alpha\delta_{MN})$. Hence, (91) can be rewritten as:

[TABLE]

or

[TABLE]

where we define $\boldsymbol{\epsilon}(T)=\mathbf{Q}^{\mathsf{T}}\mathbb{E}[\mathbf{f}(T)]$ and $\mathbf{\widetilde{r}}(T)=\mathbf{Q}^{\mathsf{T}}\mathbf{r}(T)$. We first want to bound $\mathbf{\widetilde{r}}(T)$. To this end, from (78) we derive:

[TABLE]

Therefore,

[TABLE]

Moreover, recursive applications of (93) yields:

[TABLE]

where the superscript $T$ of $\mathbf{\Sigma}$ means power instead of transposition. So,

[TABLE]

where $\|\mathbf{\Sigma}\|_{2}$ means the spectral norm of $\mathbf{\Sigma}$, i.e., the maximum singular value of $\mathbf{\Sigma}$. Since $\mathbf{\Sigma}$ is a diagonal matrix, it is easy to calculate its spectral norm as $\|\mathbf{\Sigma}\|_{2}=\max_{i=1,...,NM}|1-\alpha\delta_{i}|$. Now, we make two additional assumptions in the following.

Assumption 4**.**

$\frac{2}{1-\lambda}\min_{n}\lambda_{\min}(\mathbf{R}_{n})-4\beta(\max_{n}|\Omega_{n}|-\min_{n}|\Omega_{n}|)>0$.

Assumption 5**.**

$\alpha<\frac{2}{\max_{i=1,...,NM}\delta_{i}}$.

We note the following Weyl’s theorem [28].

Lemma 3**.**

(Weyl) Given two symmetric matrices $\mathbf{A,B}\in\mathbb{R}^{n\times n}$. Arrange the eigenvalues in increasing order, i.e., $\lambda_{1}(\mathbf{X})\leq...\leq\lambda_{n}(\mathbf{X})$ for $\mathbf{X}\in\{\mathbf{A,B,A+B}\}$. Then for any $k\in\{1,...,n\}$, we have:

[TABLE]

Making use of Weyl’s theorem, we obtain:

[TABLE]

Recall that $\mathbf{\Phi}=\operatorname{diag}(4\beta|\Omega_{1}|\mathbf{I}+\frac{2\mathbf{R}_{1}}{1-\lambda},...,4\beta|\Omega_{N}|\mathbf{I}+\frac{2\mathbf{R}_{N}}{1-\lambda})$. We have $\lambda_{\min}(\Phi)=\min_{n}\left\{4\beta|\Omega_{n}|+\frac{2}{1-\lambda}\lambda_{\min}(\mathbf{R}_{n})\right\}$. Further noting that $\mathbf{C}=\mathbf{A}\otimes\mathbf{I}$ and defining $\|\mathbf{A}\|_{1}=\max_{1\leq j\leq N}\sum_{i=1}^{N}|A_{ij}|$ to be the maximum-column-sum norm of matrix $\mathbf{A}$, we have

[TABLE]

where $\rho(\mathbf{A})=\max_{i=1,...,N}|\lambda_{i}(\mathbf{A})|$ is the spectral radius of $\mathbf{A}$. So, continuing from (100) we have:

[TABLE]

where the last step follows from Assumption 4. Since $\delta_{i}$’s are eigenvalues of $\mathbf{\Phi}-4\beta\mathbf{C}$, we thus have $\delta_{i}>0,i=1,...,MN$. This together with Assumption 5 implies that $\max_{i=1,...,NM}|1-\alpha\delta_{i}|<1$, i.e., $\|\mathbf{\Sigma}\|_{2}<1$. Hence, as $T\rightarrow\infty$, the first time of the R.H.S. of (98) converges to zero. Furthermore, since (96) holds for any $T$, the second term of the R.H.S. of (98) is always bounded above by $\frac{\sqrt{N}\alpha}{1-\|\mathbf{\Sigma}\|_{2}}\left(8\beta\max_{n=1,...,N}|\Omega_{n}|\max_{n=1,...,N}\|\mathbf{\widetilde{w}}_{n}\|_{2}+\gamma\sqrt{M}\right)$. So, in all, from (98) and the fact that $\|\mathbb{E}[\mathbf{f}(T)]\|_{2}=\|\mathbf{Q}^{\mathsf{T}}\mathbb{E}[\mathbf{f}(T)]\|_{2}=\|\boldsymbol{\epsilon}(T)\|_{2}$, we have:

[TABLE]

To summarize, we have proved the following theorem.

Theorem 1**.**

For Algorithm 2, under Assumptions 1-5, we have:

[TABLE]

where $\|\mathbf{\Sigma}\|_{2}<1$ is a constant related to the network structure and the error vector $\mathbf{f}(T)$ is defined as:

[TABLE]

Two remarks regarding the result and assumptions of the theorem are in order, respectively.

Remark 2**.**

From the R.H.S. of (110), it seems that the smaller the values of $\beta,\gamma$, the better the performance of Algorithm 2. However, this is not true. The reason is that the R.H.S. of (110) is only an upper bound on the error and when deriving this upper bound, the regularization terms are treated as biases and we want to bound their influences (similar remarks can be made to the analysis of most regularized adaptive algorithms, e.g., [11, 10]). However, in practice, as we have argued, the facts that (i) neighbors have similar weights; and (2) weights are sparse; are important prior knowledge, which can enhance the performance of the algorithm. To embody this prior knowledge, we need to assign appropriate positive values to $\beta$ and $\gamma$. Another way to decrease the R.H.S. of (110) is to decrease $\alpha$. However, this comes with a cost: smaller $\alpha$ leads to slower convergence speed. So, in practice, a positive value for $\alpha$ cannot be too small.

Remark 3**.**

Assumptions 1-3 are standard hypotheses to analyze the convergence behaviors of adaptive algorithms. Algorithm 5, i.e., small enough step size $\alpha$, is also common in gradient-descent-type algorithms such as LMS and subgradient method. The only nonconventional assumption is Assumption 4. In particular, Assumption 4 holds when $\lambda$ is sufficiently close to 1, which is the case (e.g., $\lambda=0.995,0.999$) for most applications of RLS [1, 10, 11]. This is not surprising because in most applications the weight vectors vary fairly slowly (if the weight vectors change too drastically across time, then virtually no adaptive algorithms can track them well as the past data become useless and mere one shot data at current time are not enough for estimating the high dimensional weight vectors) and thus $\lambda$ is chosen very close to 1. Another interesting observation is that for networks with uniform degree, i.e., $\max_{n}|\Omega_{n}|=\min_{n}|\Omega_{n}|$, Assumption 4 is always satisfied.

V Numerical Evaluation

In this section, numerical simulations are conducted to verify the effectiveness of the proposed decentralized online ADMM algorithm (or ADMM algorithm in short), Algorithm 1, and the proposed decentralized online subgradient method (or subgradient method in short), Algorithm 2. The performance of the global offline optimization of problem (2) (global optimizor henceforth) is shown as a benchmark. For the sake of comparison, the performance of the distributed single task sparse RLS algorithm in [11] (DSPARLS henceforth) is also presented to highlight the impact of multitask.

We consider a network with $N=20$ nodes and 40 random edges so that the average node degree is 4. The network topology is illustrated in Fig. 1. The dimension of the input data is $M=20$. Each entry of the input data sequence $\{\mathbf{u}_{n}(t)\}$ is generated according to the uniform distribution over the interval $[0,1]$ independently. The noise sequence $\{e_{n}(t)\}$ is generated according to the uniform distribution on $[0,N_{0}]$ independently, where $N_{0}$ is a constant controlling the noisy level of the observations. To achieve sparsity, we let 18 entries (whose positions are randomly selected) of the true weight vectors $\widetilde{\mathbf{w}}_{n}(t)$ be zero. The two remaining entries $\widetilde{\mathbf{w}}_{n}^{\text{part}}(0)\in\mathbb{R}^{2}$ of the initial weight vectors are generated in a way that enforces similarity between neighbors. Specifically, we first generate $N$ i.i.d. two dimensional random vectors $\{\boldsymbol{\phi}_{n}\}_{n=1,...,N}$ uniformly distributed on $[0,1]^{2}$. Then, we solve the following optimization problem to obtain $\widetilde{\mathbf{w}}_{n}^{\text{part}}(0)$:

[TABLE]

which promotes similarity between neighbors and can be easily solved as the objective function is a convex quadratic function. To capture the slowly time-variant trait of the weight vectors, the increment from $\widetilde{\mathbf{w}}_{n}(t)$ to $\widetilde{\mathbf{w}}_{n}(t+1)$, i.e., $\widetilde{\mathbf{w}}_{n}(t+1)-\widetilde{\mathbf{w}}_{n}(t)$ is generated by uniform distribution on $[-0.5N_{1},0.5N_{1}]$ independently across time and nodes, where $N_{1}$ is a constant controlling the varying rate of the weight vectors.

Now, we choose the regularization parameters and forgetting factor as $\beta=\gamma=1,\gamma=0.995$. We consider two scenarios in terms of the noise level $N_{1}$ and the varying rate of weight vectors $N_{2}$. In Scenario 1, $N_{1}=0.1,N_{2}=0.02$ while in Scenario 2, $N_{1}=0.3,N_{2}=0.05$. The latter scenario has noisier observations and weight vectors which vary faster. Thus, the weight vectors of Scenario 2 are more difficult to track than those of Scenario 1. For the proposed ADMM algorithm, the proposed subgradient method, the DSPARLS algorithm in [11] and the global optimizor, we plot the relative errors (defined to be $\|\mathbf{w}(t)-\widetilde{\mathbf{w}}(t)\|_{2}/\|\widetilde{\mathbf{w}}(t)\|_{2}$, where $\mathbf{w}(t)$ and $\widetilde{\mathbf{w}}(t)$ are concatenations of $\mathbf{w}_{n}(t)$ and $\widetilde{\mathbf{w}}_{n}(t)$ of all nodes, respectively) as functions of time indices, i.e., the learning curves, under Scenario 1 (Fig. 2-(a)) and Scenario 2 (Fig. 2-(b)), respectively. Each learning curve is the average of 300 independent trials. Several interesting observations can be made from Fig. 2. First, the relative errors of both the proposed ADMM algorithm and the proposed subgradient method can converge to that of the global optimizor, i.e., the performance benchmark, as the observation data accumulate. On the contrary, the relative error of DSPARLS does not converge to that of the global optimizor. This highlights the effectiveness of the two proposed algorithms when tracking multitask weight vectors, which cannot be tracked well by existing method (DSPARLS in this case) for the single task situation. Second, comparisons between the learning curves of the proposed two algorithms indicate that the proposed ADMM algorithm needs much fewer observations, or equivalently much less time (about 100 time units), to track the weight vectors accurately than the proposed subgradient method does (about 600 time units). This is not surprising as dual domain methods generally converge faster than primal domain methods in the literature of optimization theory [17, 29]. However, the advantage of the proposed ADMM algorithm in convergence speed comes at the cost of higher computational overhead per time unit than the proposed subgradient method. This accuracy-complexity tradeoff makes the proposed two algorithms appropriate for different applications depending on the computational capability of devices and needed tracking accuracy. Third, as one expects, Scenario 1 has better tracking performance than Scenario 2: the ultimate relative error of the proposed algorithms in Scenario 1 is about 0.067 while that of Scenario 2 is about 0.17. So, higher noise level and faster varying speed of the weight vectors do result in lower tracking accuracy.

Next, we investigate the tracking performance of each individual node. To this end, in Fig. 3, we show the relative tracking errors of each node at time $200$ and $500$ under Scenario 1 and Scenario 2. Several remarks are in order. First, we note that in all four cases of Fig. 3, the red curve (the proposed ADMM algorithm) and the blue curve (the global optimizor) coincide precisely for every node. This further confirms the previous observation from Fig. 2 that the performance of the proposed ADMM can converge to that of the global optimizor quickly. Second, the proposed subgradient method, though performs poorly at time 200, has relative errors close to those of the global optimizor at time 500. This suggests that the proposed subgradient method eventually has performance close to the benchmark (the global optimizor). But this good performance necessitates longer time (or equivalently more data) compared to the proposed ADMM algorithm. Third, the performance of the single task learning algorithm DSPARLS never converge to that of the global optimizor. In particular, from Fig. 3-(b)(d), the performance of DSPARLS is worst at node 5 and node 17. Recall the network topology in Fig. 1 and we see that these two nodes are loosely connected to other nodes, i.e., their degree is low. Thus, the weight vectors at these two nodes can potentially deviate far from the weight vectors at the rest of the nodes and thus violate the single task assumption of DSPARLS the most among all nodes. This partially explains the poor performance of DSPARLS at nodes 5 and 17.

Previous experiments indicate that the proposed ADMM algorithm possesses faster and more accurate tracking performance than the proposed subgradient method. Next, we conduct a more thorough performance comparison between the proposed two algorithms for different regularization parameters $\beta,\gamma$ and different forgetting factors $\lambda$. We note that the global optimizor usually converges well before time $1000$ and we denote the steady relative error of the global optimizor at time $1000$ as $\check{e}$. We say a simulation trial of an algorithm (either the proposed ADMM algorithm or the proposed subgradient method) is successful if, before time 1000, there exists a time window (i.e., interval) of length 20 over which the average relative error of the algorithm is lower than $1.1\check{e}$. The basic parameter setup is $\gamma=\beta=1,\lambda=0.995$. In each of the subfigures in Fig. 4-(a)(c)(e), we vary one parameter while keep the remaining two parameters the same as the basic setup. For each parameter setup, we conduct 100 independent trials and plot the number of successful trials in Fig. 4-(a)(c)(e). We observe that (i) the proposed ADMM algorithm is always successful, i.e., it can always converge to the steady performance of the global optimizor; (ii) the proposed subgradient method is successful in most trials as long as the forgetting factor $\lambda$ is sufficiently close to 1, which is the case in most applications (e.g., $\lambda=0.995$) as the weight vectors are varying very slowly and a large $\lambda$ is needed for tracking them. Note that the large $\lambda$ is reminiscent of the Assumption 4 in the performance analysis of the proposed subgradient method. See Remark 3 for a justification of large $\lambda$. Moreover, we investigate the average time needed to reach success (defined to be the middle point of the first time window over which the average relative error is less than $1.1\check{e}$) among successful trials. The results are shown in Fig. 4-(b)(d)(f). We remark that the proposed ADMM mostly needs no more than 150 time units to be successful, i.e., be close to the steady performance of the global optimizor, while it takes the proposed subgradient method a much longer time (around 600 time units) to be successful. This further confirms our previous assertion that the proposed ADMM algorithm possesses faster tracking performance than the proposed subgradient method.

VI Conclusion

In this paper, we study the decentralized sparse multitask RLS problem. We first propose a decentralized online ADMM algorithm for the formulated RLS problem. We simplify the algorithm so that each iteration of the ADMM algorithm consists of simple closed form computations and each node only needs to store and update one $M\times M$ matrix and six $M$ dimensional vectors. We show that the gap between the iterations of the proposed ADMM algorithm and the optimal points of the formulated RLS problem will converge to zero. Moreover, to further reduce the computational complexity, we propose a decentralized online subgradient method. We show that the iterations of the proposed subgradient method will track the true weight vectors with errors upper bounded by some constant related to the network topology and algorithm parameters. Compared with the ADMM algorithm, the subgradient method enjoys lower computational complexity at the expense of slower convergence speed. The proposed algorithms are corroborated by numerical experiments, which highlight the effectiveness of the proposed algorithms and the accuracy-complexity tradeoff between the proposed two algorithms.