Study of Robust Distributed Diffusion RLS Algorithms with Side Information for Adaptive Networks

TL;DR

This paper introduces robust diffusion RLS algorithms for adaptive networks that effectively handle impulsive noise by incorporating side information and reducing computational complexity, with proven convergence and superior performance.

Contribution

The paper presents novel diffusion RLS algorithms with side information and reduced complexity, enhancing robustness against impulsive noise in adaptive networks.

Findings

Algorithms outperform existing methods in impulsive noise scenarios.

Proposed methods demonstrate mean-square convergence.

Reduced complexity makes the algorithms practical for real-world applications.

Abstract

This work develops robust diffusion recursive least squares algorithms to mitigate the performance degradation often experienced in networks of agents in the presence of impulsive noise. The first algorithm minimizes an exponentially weighted least-squares cost function subject to a time-dependent constraint on the squared norm of the intermediate update at each node. A recursive strategy for computing the constraint is proposed using side information from the neighboring nodes to further improve the robustness. We also analyze the mean-square convergence behavior of the proposed algorithm. The second proposed algorithm is a modification of the first one based on the dichotomous coordinate descent iterations. It has a performance similar to that of the former, however its complexity is significantly lower especially when input regressors of agents have a shift structure and it is well suited to practical implementation. Simulations show the superiority of the proposed algorithms over previously reported techniques in various impulsive noise scenarios.

Tables

Table 1. TABLE I: Proposed R- d d \rm d RLS Algorithm Allied with the NC Method

Parameters: , $\lambda$, $\delta$ and $E_c$ (R-dRLS part);

$\varrho$, $\tau$ and $t_{\text{th}}$ (NC part)

Initialization: ${𝒘}_{k,0}=\mathrm{𝟎}$, ${𝑷}_{k,0}={\delta }^{-1}{𝑰}_M$ and ${\xi }_k\left(0\right)=E_c\frac{{\sigma }_{d,k}^2}{M{\sigma }_{u,k}^2}$ (R-dRLS part);

${\Theta }_{\text{old},k}={\Theta }_{\text{new},k}=0$, ${\sigma }_{e,k}^2=0$, $V_t=\varrho M$

and $V_d=0.75V_t$ (NC part)

for iteration $i=1,2,3,\dots$

for node $k=1,2,3,\dots ,N$

[R-dRLS part:]

$e_k\left(i\right)=d_k\left(i\right)-{𝒖}_{k,i}^T{𝒘}_{k,i-1}$

${𝑷}_{k,i}=\frac{1}{\lambda }\left({𝑷}_{k,i-1}-\frac{{𝑷}_{k,i-1}{𝒖}_{k,i}{𝒖}_{k,i}^T{𝑷}_{k,i-1}}{\lambda +{𝒖}_{k,i}^T{𝑷}_{k,i-1}{𝒖}_{k,i}}\right)$

${𝒈}_{k,i}={𝑷}_{k,i}{𝒖}_{k,i}$

${𝝍}_{k,i}={𝒘}_{k,i-1}+\min [\frac{\sqrt{{\xi }_k\left(i-1\right)}}{{\parallel {𝒈}_{k,i}\parallel }_2\left|e_k\left(i\right)\right|},\mathrm{\hspace{0.33em}1}]{𝒈}_{k,i}{e}_k\left(i\right)$

${𝒘}_{k,i}=\sum _{m\in {𝒩}_k}{c}_{m,k}{𝝍}_{m,i}$

[NC part:]

Step 1: to compute ${\Delta }_k\left(i\right)$

if $i=n{V}_t,n=0,1,2,\dots$

${𝒂}_{k,i}^T=ℛ\left([\frac{e_k^2\left(i\right)}{{\Vert {𝒖}_{k,i}\Vert }_2^2},\frac{e_k^2\left(i-1\right)}{{\Vert {𝒖}_{k,i-1}\Vert }_2^2},\text{…},\frac{e_k^2\left(i-V_t+1\right)}{{\Vert {𝒖}_{k,i-V_t+1}\Vert }_2^2}]\right)$

${\sigma }_{e,k}^2\leftarrow \tau {\sigma }_{e,k}^2+\left(1-\tau \right){𝒂}_{k,i}^T𝒆$

${\Theta }_{\text{new},k}=\frac{1}{V_t-V_d}\sum _{m\in {𝒩}_k}{c}_{m,k}{\sigma }_{e,m}^2$

${\Delta }_k\left(i\right)=\frac{{\Theta }_{\text{new},k}-{\Theta }_{\text{old},k}}{{\xi }_k\left(i-1\right)}$

end

Step 2: to reset ${\xi }_k\left(i\right)$

if ${\Delta }_k\left(i\right)>t_{\text{th}}$

${\zeta }_k\left(i\right)={\xi }_k\left(0\right)$, ${𝑷}_{k,i}={𝑷}_{k,0}$

elseif   ${\Theta }_{\text{new},k}>{\Theta }_{\text{old},k}$

${\zeta }_k\left(i\right)={\xi }_k\left(i-1\right)+\left({\Theta }_{\text{new},k}-{\Theta }_{\text{old},k}\right)$

else

${\zeta }_k\left(i\right)=\beta {\xi }_k\left(i-1\right)+\left(1-\beta \right)\min [{\parallel {𝒈}_{k,i}\parallel }_2^2{e}_k^2\left(i\right),{\xi }_k\left(i-1\right)]$

end

${\Theta }_{\text{old},k}={\Theta }_{\text{new},k}$

${\xi }_k\left(i\right)=\sum _{m\in {𝒩}_k}{c}_{m,k}{\zeta }_m\left(i\right)$

end

end

Table 2. TABLE II: DCD Method for Solving ( 51 ).

Parameters:$H,N_u,M_b$,

Initialization: $\Delta {\widehat{𝒘}}_{k,i}=\mathrm{𝟎},{𝒓}_{k,i}={𝒃}_{k,i},y=1,\mu =H/2$

for $j=1,\dots ,N_u$

$l=\arg \underset{j=1,\dots ,M}{\max }\left\{\right|{\left[{𝒓}_{k,i}\right]}_j|\}$

while $\left|{\left[{𝒓}_{k,i}\right]}_l\right|\le (\mu /2)\left[{𝚽}_{k,i}\right]{}_{l,l}$ and $y\le M_b$

$y=y+1$, $\mu =\mu /2$

end

if $y>M_b$

break

else

${\left[\Delta {\widehat{𝒘}}_{k,i}\right]}_l\leftarrow {\left[\Delta {\widehat{𝒘}}_{k,i}\right]}_l+\mu \text{sign}\left({\left[{𝒓}_{k,i}\right]}_l\right)$

${𝒓}_{k,i}\leftarrow {𝒓}_{k,i}-\mu \text{sign}\left({\left[{𝒓}_{k,i}\right]}_l\right){\left[{𝚽}_{k,i}\right]}_{:,l}$

end

end

Table 3. TABLE III: Proposed DCD-R- d d \rm d RLS Algorithm

Parameters: , $\lambda$, $\delta$ and $E_c$

Initialization: ${𝒘}_{k,0}=\mathrm{𝟎}$, ${𝚽}_{k,0}=\delta {𝑰}_M$ and ${\xi }_k\left(0\right)=E_c\frac{{\sigma }_{d,k}^2}{M{\sigma }_{u,k}^2}$

for each node $k$:

$e_k\left(i\right)=d_k\left(i\right)-{𝒖}_{k,i}^T{𝒘}_{k,i-1}$

${𝚽}_{k,i}=\lambda {𝚽}_{k,i-1}+{𝒖}_{k,i}{𝒖}_{k,i}^T$

${𝒃}_{k,i}=\lambda {𝒓}_{k,i-1}+e_k\left(i\right){𝒖}_{k,i}$

Using DCD to solve ${𝚽}_{k,i}\Delta {𝒘}_{k,i}={𝒃}_{k,i}$, yielding

$\Delta {\widehat{𝒘}}_{k,i}=\Delta {\widehat{𝒘}}_{k,i}^{\left(1\right)}$ and ${𝒓}_{k,i}={𝒓}_{k,i}^{\left(1\right)}$

if ${\parallel \Delta {\widehat{𝒘}}_{k,i}\parallel }_2^2>{\xi }_k\left(i-1\right)$

${𝒃}_{k,i}=\lambda {𝒓}_{k,i-1}+\frac{\sqrt{{\xi }_k\left(i-1\right)}}{{\parallel \Delta {\widehat{𝒘}}_{k,i}\parallel }_2}{e}_k\left(i\right){𝒖}_{k,i}$

Using DCD to solve ${𝚽}_{k,i}\Delta {𝒘}_{k,i}={𝒃}_{k,i}$, yielding

$\Delta {\widehat{𝒘}}_{k,i}=\frac{\sqrt{{\xi }_k\left(i-1\right)}}{{\parallel \Delta {\widehat{𝒘}}_{k,i}^{\left(2\right)}\parallel }_2}\Delta {\widehat{𝒘}}_{k,i}^{\left(2\right)}$ and ${𝒓}_{k,i}={𝒓}_{k,i}^{\left(2\right)}$

end

${𝝍}_{k,i}={𝒘}_{k,i-1}+\Delta {\widehat{𝒘}}_{k,i}$

${𝒘}_{k,i}=\sum _{m\in {𝒩}_k}{c}_{m,k}{𝝍}_{m,i}$

${\zeta }_k\left(i\right)=\beta {\xi }_k\left(i-1\right)+\left(1-\beta \right){\parallel \Delta {\widehat{𝒘}}_{k,i}\parallel }_2^2$

${\xi }_k\left(i\right)=\sum _{m\in {𝒩}_k}{c}_{m,k}{\zeta }_m\left(i\right)$

Table 4. TABLE IV: Computational complexity of algorithms for node k per time instant.

Algorithms	Multiplications	Additions	Divisions	Square-root
dLMS	$n_{k}M+2M+1$	$n_{k}M+M$	-	-
dRLS	$n_{k}M+4M^2+3M$	$n_{k}M+3M^2$	$M$	-
DCD-dRLS
without shift structure in input	$n_{k}M+2M^2+3M$	$n_{k}M+M^2+2M+C_{dcd}^{+}$	-	-
DCD-dRLS
with shift structure in input	$n_{k}M+5M$	$n_{k}M+3M+C_{dcd}^{+}$	-	-
R-dRLS	$n_k\left(M+1\right)+4M^2+4M+5$	$n_k\left(M+1\right)+3M^2+M+1$	$M+1$	1
DCD-R-dRLS
without shift structure in input	$n_k\left(M+1\right)+2M^2+4M+3\kappa M+2$	$n_k\left(M+1\right)+M^2+3M+\kappa \left(2M-1+C_{dcd}^{+}\right)+C_{dcd}^{+}$	$2\kappa$	$2\kappa$
DCD-R-dRLS
with shift structure in input	$n_k\left(M+1\right)+6M+3\kappa M+2$	$n_k\left(M+1\right)+4M+\kappa \left(2M-1+C_{dcd}^{+}\right)+C_{dcd}^{+}$	$2\kappa$	$2\kappa$

Equations

$$d_k(i)={\mathbf{u}}_{k,i}^T{\mathbf{w}}^o+v_k(i),$$

$$\begin{array}[]{rcl}\begin{aligned} &\bm{w}_{i}=\arg\min\limits_{\bm{w}}\\ &\left\{\lambda^{i+1}\delta\lVert\bm{w}\rVert_{2}^{2}+\sum\limits_{j=0}^{i}\lambda^{i-j}\sum\limits_{k=1}^{N}\left(d_{k}(j)-\bm{u}_{k,i}^{T}\bm{w}\right)^{2}\right\},\end{aligned}\end{array}$$

$$\begin{array}[]{rcl}\begin{aligned} J_{k}(\bm{\psi}_{k,i})=&\lVert\bm{\psi}_{k,i}-\bm{w}_{k,i-1}\rVert_{\bm{B}_{k,i}}^{2}\\ &+[d_{k}(i)-\bm{u}_{k,i}^{T}\bm{\psi}_{k,i}]^{2},\end{aligned}\end{array}$$

$$\begin{array}[]{rcl}\begin{aligned} \bm{\Phi}_{k,i}\triangleq&\lambda^{i+1}\delta\bm{I}_{M}+\sum\limits_{j=0}^{i}\lambda^{i-j}\bm{u}_{k,j}\bm{u}_{k,j}^{T}\\ =&\lambda\bm{\Phi}_{k,i-1}+\bm{u}_{k,i}\bm{u}_{k,i}^{T}\end{aligned}\end{array}$$

$$\begin{array}[]{rcl}\begin{aligned} \bm{\psi}_{k,i}=\bm{w}_{k,i-1}+\bm{P}_{k,i}\bm{u}_{k,i}e_{k}(i),\end{aligned}\end{array}$$

$$\begin{array}[]{rcl}\begin{aligned} e_{k}(i)=d_{k}(i)-\bm{u}_{k,i}^{T}\bm{w}_{k,i-1}\end{aligned}\end{array}$$

$$\begin{array}[]{rcl}\begin{aligned} \bm{P}_{k,i}&\triangleq\bm{\Phi}_{k,i}^{-1}\\ &=\frac{1}{\lambda}\left(\bm{P}_{k,i-1}-\frac{\bm{P}_{k,i-1}\bm{u}_{k,i}\bm{u}_{k,i}^{T}\bm{P}_{k,i-1}}{\lambda+\bm{u}_{k,i}^{T}\bm{P}_{k,i-1}\bm{u}_{k,i}}\right),\end{aligned}\end{array}$$

$$\begin{array}[]{rcl}\begin{aligned} \bm{w}_{k,i}=\sum\limits_{m\in\mathcal{N}_{k}}c_{m,k}\bm{\psi}_{m,i},\end{aligned}\end{array}$$

$$\sum _{m\in {\mathcal{N}}_k}{c}_{m,k}=1\text{, and }{c}_{m,k}=0\text{ if }m\notin {\mathcal{N}}_k.$$

$$\begin{array}[]{rcl}\begin{aligned} \lVert\bm{\psi}_{k,i}-\bm{w}_{k,i-1}\rVert_{2}^{2}\leq\xi_{k}(i-1),\end{aligned}\end{array}$$

$$\begin{array}[]{rcl}\begin{aligned} \lVert\bm{g}_{k,i}\rVert_{2}\lvert e_{k}(i)\rvert\leq\sqrt{\xi_{k}(i-1)},\end{aligned}\end{array}$$

$$\begin{array}[]{rcl}\begin{aligned} \bm{\psi}_{k,i}=\bm{w}_{k,i-1}+\sqrt{\xi_{k}(i-1)}\frac{\bm{g}_{k,i}}{\lVert\bm{g}_{k,i}\rVert_{2}}\text{sign}(e_{k}(i)),\\ \end{aligned}\end{array}$$

$$\begin{array}[]{rcl}\begin{aligned} \bm{\psi}_{k,i}=\bm{w}_{k,i-1}+\min\left[\frac{\sqrt{\xi_{k}(i-1)}}{\lVert\bm{g}_{k,i}\rVert_{2}\lvert e_{k}(i)\rvert},\;1\right]\bm{g}_{k,i}e_{k}(i).\\ \end{aligned}\end{array}$$

$$\begin{array}[]{rcl}\begin{aligned} \zeta_{k}(i)=&\beta\xi_{k}(i-1)+(1-\beta)\left\|\bm{\psi}_{k,i}-\bm{w}_{k,i-1}\right\|_{2}^{2}\\ =\beta&\xi_{k}(i-1)+(1-\beta)\min[\lVert\bm{g}_{k,i}\rVert_{2}^{2}e_{k}^{2}(i),\xi_{k}(i-1)],\\ \xi_{k}(i)=&\sum\limits_{m\in\mathcal{N}_{k}}c_{m,k}\zeta_{m}(i),\end{aligned}\end{array}$$

$${\sigma }_{e,k}^2\leftarrow \tau {\sigma }_{e,k}^2+(1-\tau ){\mathbf{a}}_{k,i}^T\mathbf{e},$$

$$\begin{array}[]{rcl}\begin{aligned} \widetilde{\bm{w}}_{k,i}&\triangleq\bm{w}^{o}-\bm{w}_{k,i},\\ \widetilde{\bm{\psi}}_{k,i}&\triangleq\bm{w}^{o}-\bm{\psi}_{k,i}.\end{aligned}\end{array}$$

$$\begin{array}[]{rcl}\begin{aligned} \widetilde{\bm{\psi}}_{k,i}=\widetilde{\bm{w}}_{k,i-1}-\sqrt{\frac{\zeta_{k}(i)-\beta\xi_{k}(i-1)}{1-\beta}}\frac{\bm{g}_{k,i}}{\lVert\bm{g}_{k,i}\rVert_{2}}\text{sign}(e_{k}(i)),\end{aligned}\end{array}$$

$$\begin{array}[]{rcl}\begin{aligned} \widetilde{\bm{w}}_{k,i}=\sum\limits_{m\in\mathcal{N}_{k}}c_{m,k}\widetilde{\bm{\psi}}_{m,i}.\end{aligned}\end{array}$$

$$\begin{array}[]{rcl}\begin{aligned} E&\left\{\lVert\widetilde{\bm{\psi}}_{k,i}\rVert_{2}^{2}\right\}=E\left\{\lVert\widetilde{\bm{w}}_{k,i-1}\rVert_{2}^{2}\right\}\\ &-2E\left\{\sqrt{\frac{\zeta_{k}(i)-\beta\xi_{k}(i-1)}{1-\beta}}\frac{\widetilde{\bm{w}}_{k,i-1}^{T}\bm{g}_{k,i}}{\lVert\bm{g}_{k,i}\rVert_{2}}\text{sign}(e_{k}(i))\right\}\\ &+\frac{E\left\{\zeta_{k}(i)\right\}-\beta E\left\{\xi_{k}(i-1)\right\}}{1-\beta}.\end{aligned}\end{array}$$

$$\begin{array}[]{rcl}\begin{aligned} E\left\{\lVert\widetilde{\bm{w}}_{k,i}\rVert_{2}^{2}\right\}\leq\sum\limits_{m\in\mathcal{N}_{k}}c_{m,k}E\left\{\lVert\widetilde{\bm{\psi}}_{m,i}\rVert_{2}^{2}\right\}.\end{aligned}\end{array}$$

$$\begin{array}[]{rcl}\begin{aligned} &E\left\{\sqrt{\frac{\zeta_{k}(i)-\beta\xi_{k}(i-1)}{1-\beta}}\frac{\widetilde{\bm{w}}_{k,i-1}^{T}\bm{g}_{k,i}}{\lVert\bm{g}_{k,i}\rVert_{2}}\text{sign}(e_{k}(i))\right\}\approx\\ &\sqrt{\frac{E\left\{\zeta_{k}(i)\right\}-\beta E\left\{\xi_{k}(i-1)\right\}}{1-\beta}}E\left\{\frac{\widetilde{\bm{w}}_{k,i-1}^{T}\bm{g}_{k,i}}{\lVert\bm{g}_{k,i}\rVert_{2}}\text{sign}(e_{k}(i))\right\}.\end{aligned}\end{array}$$

$$\begin{array}[]{rcl}\begin{aligned} E\left\{\lVert\widetilde{\bm{\psi}}_{k,i}\rVert_{2}^{2}\right\}=&E\left\{\lVert\widetilde{\bm{w}}_{k,i-1}\rVert_{2}^{2}\right\}-\\ &2\sqrt{\frac{E\left\{\zeta_{k}(i)\right\}-\beta E\left\{\xi_{k}(i-1)\right\}}{1-\beta}}\times\\ &\underbrace{E\left\{\frac{\widetilde{\bm{w}}_{k,i-1}^{T}\bm{g}_{k,i}}{\lVert\bm{g}_{k,i}\rVert_{2}}\text{sign}(e_{k}(i))\right\}}\limits_{(a)}+\\ &\frac{E\left\{\zeta_{k}(i)\right\}-\beta E\left\{\xi_{k}(i-1)\right\}}{1-\beta}.\end{aligned}\end{array}$$

$$\begin{array}[]{rcl}\begin{aligned} \lim_{i\rightarrow\infty}E\left\{\bm{\Phi}_{k,i}\right\}=\frac{\bm{R}_{k}}{1-\lambda}\triangleq\bar{\bm{\Phi}}_{k}.\end{aligned}\end{array}$$

$$\begin{array}[]{rcl}\begin{aligned} &E\left\{\frac{\widetilde{\bm{w}}_{k,i-1}^{T}\bm{g}_{k,i}}{\lVert\bm{g}_{k,i}\rVert_{2}}\text{sign}(e_{k}(i))\right\}\\ &\;\;\;\;\;\;\;\;\;\;\approx E\left\{\frac{\widetilde{\bm{w}}_{k,i-1}^{T}\bar{\bm{\Phi}}_{k}^{-1}\bm{u}_{k,i}}{\sqrt{\bm{u}_{k,i}^{T}\bar{\bm{\Phi}}_{k}^{-2}\bm{u}_{k,i}}}\text{sign}(e_{k}(i))\right\}\\ &\;\;\;\;\;\;\;\;\;\;=E\left\{\frac{\widetilde{\bm{w}}_{k,i-1}^{T}\bm{R}_{k}^{-1}\bm{u}_{k,i}}{\sqrt{\bm{u}_{k,i}^{T}\bm{R}_{k}^{-2}\bm{u}_{k,i}}}\text{sign}(e_{k}(i))\right\}.\end{aligned}\end{array}$$

$$\begin{array}[]{rcl}\begin{aligned} &E\left\{\frac{\widetilde{\bm{w}}_{k,i-1}^{T}\bm{R}_{k}^{-1}\bm{u}_{k,i}}{\sqrt{\bm{u}_{k,i}^{T}\bm{R}_{k}^{-2}\bm{u}_{k,i}}}\text{sign}(e_{k}(i))\right\}\approx\\ &\;\;\;\;\;\;\;\;\;\;\chi_{k}E\left\{e_{a,k}^{\bm{R}_{k}^{-1}}(i)\text{sign}\left(e_{a,k}(i)+v_{k}(i)\right)\right\},\end{aligned}\end{array}$$

$$\begin{array}[]{rcl}\begin{aligned} e_{a,k}(i)&\triangleq\widetilde{\bm{w}}_{k,i-1}^{T}\bm{u}_{k,i},\\ e_{a,k}^{\bm{R}_{k}^{-1}}(i)&\triangleq\widetilde{\bm{w}}_{k,i-1}^{T}\bm{R}_{k}^{-1}\bm{u}_{k,i},\\ \chi_{k}&=E\left\{\frac{1}{\sqrt{\bm{u}_{k,i}^{T}\bm{R}_{k}^{-2}\bm{u}_{k,i}}}\right\}.\end{aligned}\end{array}$$

$$\begin{array}[]{rcl}\begin{aligned} E\{\text{sign}(e)u\}=p_{r}E\{\text{sign}(e_{1})u\}+(1-p_{r})E\{\text{sign}(e_{2})u\}.\end{aligned}\end{array}$$

$$\begin{array}[]{rcl}\begin{aligned} E\left\{\left.e_{a,k}^{\bm{R}_{k}^{-1}}(i)\text{sign}\left(e_{a,k}(i)+v_{k}(i)\right)\right|\widetilde{\bm{w}}_{k,i-1}\right\}\\ =\varpi_{k,i}E\left\{\left.e_{a,k}^{\bm{R}_{k}^{-1}}(i)e_{a,k}(i)\right|\widetilde{\bm{w}}_{k,i-1}\right\},\end{aligned}\end{array}$$

$$\begin{array}[]{rcl}\begin{aligned} \varpi_{k,i}=\sqrt{\frac{2}{\pi}}\left\{\frac{p_{r,k}}{\sqrt{E\{e_{a,k}^{2}(i)\}+(\hbar+1)\sigma_{\theta}^{2}\}}}\right.&\\ +\left.\frac{1-p_{r,k}}{\sqrt{E\{e_{a,k}^{2}(i)\}+\sigma_{\theta}^{2}\}}}\right\}\neq 0&,\end{aligned}\end{array}$$

$$\begin{array}[]{rcl}\begin{aligned} E&\left\{e_{a,k}^{\bm{R}_{k}^{-1}}(i)\text{sign}\left(e_{a,k}(i)+v_{k}(i)\right)\right\}\\ &=E\left\{E\left\{\left.e_{a,k}^{\bm{R}_{k}^{-1}}(i)\text{sign}\left(e_{a,k}(i)+v_{k}(i)\right)\right|\widetilde{\bm{w}}_{k,i-1}\right\}\right\}\\ &=\varpi_{k,i}E\left\{e_{a,k}^{\bm{R}_{k}^{-1}}(i)e_{a,k}(i)\right\}\\ &\stackrel{{\scriptstyle(a)}}{{=}}\varpi_{k,i}E\left\{\lVert\widetilde{\bm{w}}_{k,i-1}\rVert_{2}^{2}\right\},\end{aligned}\end{array}$$

Full text

I Introduction

Over the past decade, distributed parameter estimation over wireless sensor networks with multiple nodes (agents) has attracted much attention. It only relies on the local data exchange between interconnected nodes, and therefore removes the requirement of a powerful central processor and, as such, reduces communications bandwidth of the traditional centralized estimation whilst retaining similar estimation performance [2, 3]. Distributed estimation has been applied to target localization [4], clustering [5], frequency estimation [6] and spectrum estimation in Cognitive radio (CR) [7, 8].

I-A Prior and Related Work

According to the cooperation strategies between interconnected nodes, existing algorithms can be categorized as incremental [9], consensus [10], and diffusion [11, 12, 13, 14] types. Among these, the diffusion strategy is popular, because it does not require a Hamiltonian cycle path as in the incremental type, thereby it is more robust to nodes/links failures; it is stable and shows a faster convergence rate and a lower mean-square error than that of the consensus approach. Several diffusion algorithms were proposed, e.g., diffusion least mean square (dLMS) algorithm [11] and its variable step size variants [15, 16].

In practice, the measurements can be corrupted by non-Gaussian noise with impulsiveness. Impulsive noise has small occurrence probability but much higher amplitude than the nominal measurements. It may occur due to atmospheric phenomena, or man-made due to either electric machinery in the operation environment [17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43]. Other examples are keyboard clicking or pen dropping in teleconference [44], double-talk in echo cancellation [45], biological noise [46] or ice cracking [47] in various underwater signals, out-of-band spectral leakage in CR [48], etc. In such scenarios, the conventional algorithms like the dLMS designed for Gaussian noise would undergo a significant performance deterioration. To this end, many robust distributed algorithms have been proposed. Some algorithms are based on the instantaneous gradient-descent method to minimize different robust criteria, for instance, the diffusion error nonlinearity (dEN) [49], diffusion least mean p-th power (dLMP) [50], diffusion sign error LMS (dSE-LMS) [51], and diffusion maximum correntropy [52] algorithms. Moreover, because of the insensitivity of correntropy to impulsive noise, the maximum total correntropy diffusion algorithm was proposed in [53] for the case of large outliers in communication links. Nevertheless, their main limitation is slow convergence especially when the nodes’ input signals are colored (highly correlated). As shown in [49], the dEN algorithm converges slower than the dSE-LMS algorithm. In [54], by resorting to the adaptive projected subgradient method, a robust diffusion algorithm was developed which projects the output errors onto halfspaces defined by Huber’s error function at each node, thereby speeding up the convergence. However, the setting of the parameters controlling the algorithm’s robustness requires prior knowledge of the noise distribution which is often unavailable.

It is well-known that due to the exponentially weighted least squares (EWLS) criterion, the diffusion recursive LS (dRLS) algorithm provides fast convergence even for colored signals  [55, 56]. By means of the alternating direction method of multipliers to solve the EWLS problem, Mateos et al. proposed another type of distributed RLS algorithm [57]. Following this algorithm, to reduce computation and communication costs, its variants were presented by censoring observations with small innovations [58]. Likewise, these algorithms might experience convergence issues in impulsive noise environments, because impulsive noise samples are directly involved in the adaptation through output errors of nodes. For the single-agent case, many works have proposed RLS algorithms robust against impulsive noise, e.g., [59, 60]. However, distributed RLS-based techniques that are robust to impulsive noise have not been well investigated. The study in [61] develops the diffusion recursive least p-th power (dRLP) algorithm, while its robustness relies on the value of $p$ as the dLMP does.

Analogous to the RLS, distributed RLS requires high computational complexity. Apart from this, it may also suffer from numerical instability due to accumulation of round-off errors in finite-precision implementations [62]. Aiming to address these problems, an efficient alternative method is the dichotomous coordinate-descent (DCD) that solves a system of normal equations associated with the RLS-type algorithms [63, 62, 64, 65]. In particular, the DCD method only involves shift and addition operations, thus the DCD-based RLS algorithms reduce the computational cost and improve the numerical stability in contrast with the original RLS counterparts, whilst preserving comparable estimation performance. For this reason, reference [66] also explored the use of the DCD in distributed networks, and developed the DCD-dRLS algorithm. It is worth mentioning that, however, the development of the DCD-based algorithms in impulsive noise environments has not been studied in single nor multi -agent scenarios.

I-B Contributions

The focus of this paper is on developing robust distributed RLS algorithms for scenarios with impulsive noise. Specifically, our contributions are listed as follows:

1. A robust dRLS (R-dRLS) algorithm is developed by extending the framework of [59] to multi-agent scenarios with a diffusion distributed strategy. To ensure that the proposed R-dRLS algorithm has good convergence performance after an abrupt change in the set of parameters to be estimated, we also propose a diffusion-based non-stationary control (NC) method.

2. Theoretical insights into the mean square steady-state and evolution behaviors of the R-dRLS algorithm in impulsive noise environment are presented.

3. We employ the DCD method for developing recursions used in the adaptation step of the R-dRLS algorithm, resulting in the DCD-R-dRLS algorithm with similar learning performance. Remarkably, the DCD-R-dRLS algorithm brings a reduction in computational complexity; especially for shift structured input regressors, it reduces the order of complexity from $\mathcal{O}(M^{2})$ to $\mathcal{O}(M)$, where $M$ is the length of the estimated vector.

4. Simulation examples are presented to demonstrate the performance of the proposed algorithms in impusive noise scenarios described by either Bernoulli-Gaussian (BG) or $\alpha$-stable processes.

In comparison to the preliminary results [1] related to this work, the current version is further developed due to the main contributions 2) and 3). We slightly improve the NC method by a smoothing operation as shown in (15). Moreover, the effectiveness of the proposed algorithms are also verified in an application to distributed spectrum estimation.

This paper is organized as follows. In Section II, the estimation problem is described. The R-dRLS algorithm is derived in Section III. Analyses of its mean square behavior are presented in Section IV. In Section V, we review the DCD algorithm and propose the DCD-R-dRLS algorithm. In Section VI, extensive simulations are presented to verify the proposed algorithms. Finally, conclusions are given in Section VII.

Notation: Throughout the paper, all vectors are column vectors. We use the parenthesis on $i$ to denote matrices and vectors, and the subscript on $i$ to denote scalars. The superscript $(\cdot)^{T}$ denotes the transpose, $\lVert\cdot\rVert_{2}$ denotes the $l_{2}$-norm of a vector, and $E\{\cdot\}$ denotes the expectation of random variables. We use $\text{col}\{\cdot\cdot\cdot\}$ to denote an enlarged column vector structured by stacking its columns on top of each other, $\text{diag}\{\cdot\cdot\cdot\}$ to yield a diagonal matrix with its arguments, and $\text{Tr}\{\cdot\}$ to denote the trace of a matrix. $\bm{I}_{M}$ is the identity matrix of size $M\times M$, $\otimes$ is the Kronecker product, and $\bm{1}$ is the column vector of length $M$ with all entries being one. For symmetric matrices $\bm{X}$ and $\bm{Y}$, the notation $\bm{X}\geq\bm{Y}$ stands for $\bm{X}-\bm{Y}\geq 0$, meaning that the matrix difference $\bm{X}-\bm{Y}$ is positive semi-definite.

II Problem Formulation

Let us consider a diffusion network with $N$ nodes located at different positions in space, as shown in Fig. 1, where each node communicates only with its neighboring nodes by a link (single-hop communication). All nodes connected directly to node $k$ (including itself) are referred to as its neighborhood, denoted as $\mathcal{N}_{k}$. At every time instant $i\geq 0$, every node has access to an $M\times 1$ input regressor vector $\bm{u}_{k,i}$ and an output measurement $d_{k}(i)$, which are related as:

[TABLE]

where $\bm{w}^{o}$ is a parameter vector of size $M\times 1$ to be estimated, and $v_{k}(i)$ is the additive noise at node k. The additive noises $v_{k}(i)$ and $v_{l}(j)$ are spatially and temporally independent for $k\neq l$ and $i\neq j$. Moreover, any $\bm{u}_{k,i}$ is independent of any $v_{l}(j)$. The model (1) is used in many applications [3, 67]. The objective of the in-network processing is to estimate $\bm{w}^{o}$, using the available data $\{\bm{u}_{k,i},d_{k}(i)\}$ collected at nodes.

For this purpose, the global EWLS estimation problem is described as [55]:

[TABLE]

where $\delta>0$ is a regularization constant, and $\lambda$ ($0<\lambda\leq 1$) is the forgetting factor. The dRLS algorithm solves (2) in a diffusion-based distributed manner [55]. As already mentioned in the introduction, the noise $v_{k}(i)$ may be non-Gaussian with impulsiveness so that the algorithms derived from (2), e.g., the dRLS algorithm, would exhibit poor convergence and even diverge. In general, when studying robust adaptive algorithms, both contaminated-Gaussian (CG) [51, 68] and $\alpha$-stable [69, 70] random processes are often used for modeling impulsive noise.

III Proposed R-$\rm d$RLS Algorithm

In this section, we derive the R-dRLS algorithm and propose a control method for endowing it with tracking capability. The diffusion strategy has two alternatives: the adapt-then-combine (ATC) and the combine-then-adapt (CTA). However, we focus only on the ATC policy, which performs first the adaptation step and then the combination step. This is based on the fact that the extension to CTA is straightforward by reversing the order of the adaptation and combination steps [2, 5]. In what follows, we neglect the notation ATC for brevity.

III-A dRLS Algorithm

To conveniently develop the R-dRLS algorithm, we re-derive here the dRLS algorithm from the following method instead of directly solving (2).

In the adaptation step, every node $k$, at time instant $i$, finds an intermediate estimate $\bm{\psi}_{k,i}$ of $\bm{w}^{o}$ by minimizing the individual local cost function:

[TABLE]

with $\bm{B}_{k,i}=\bm{\Phi}_{k,i}-\bm{u}_{k,i}\bm{u}_{k,i}^{T}$, where

[TABLE]

is the time-averaged correlation matrix for the input vector at node $k$, and $\bm{w}_{k,i-1}$ is an estimate of $\bm{w}^{o}$ at node $k$ at time instant $i-1$. Notice that the quadratic form $\lVert\bm{x}\rVert_{\bm{B}_{k,i}}^{2}\triangleq\bm{x}^{T}\bm{B}_{k,i}\bm{x}$ in (3) defines the Riemmanian distance between vectors $\bm{\psi}_{k,i}$ and $\bm{w}_{k,i-1}$, where $\bm{B}_{k,i}$ is a Riemannian metric tensor characterizing that the distance properties are not uniform along the $M$-dimensional space [71, 72].

Setting the derivative of $J_{k}(\bm{\psi}_{k,i})$ with respect to $\bm{\psi}_{k,i}$ to zero, we obtain

[TABLE]

where

[TABLE]

stands for the output error at node $k$, and

[TABLE]

with $\bm{P}_{k,i}$ initialized as $\left.\bm{P}_{k,0}=\delta^{-1}\bm{I}_{M}\right.$. The recursion (7) is the result of applying the matrix inversion lemma [67].

At the combination step, the intermediate estimates $\bm{\psi}_{m,i}$, $\left.m\in\mathcal{N}_{k}\right.$ from the neigborhood of node $k$ are linearly weighted, yielding a combined estimate $\bm{w}_{k,i}$ [3]:

[TABLE]

where the combination coefficients $\{c_{m,k}\}$ are non-negative, and satisfy:

[TABLE]

Note that $c_{m,k}$ is a weight that node $k$ assigns to the intermediate estimate $\bm{\psi}_{m,i}$ received from its neighbor node $m$. If one assumes $\bm{w}_{k,i-1}=\bm{\Phi}_{k,i-1}^{-1}\bm{z}_{k,i-1}$ in (5), where $\bm{z}_{k,i}=\lambda\bm{z}_{k,i-1}+\bm{u}_{k,i}d_{k}(i)$, (5) is a standard RLS update for node $k$. In summary, (5)-(8) formulate the dRLS algorithm. It is noteworthy that the term $\bm{P}_{k,i}\bm{u}_{k,i}$ in (5) provides the decorrelating ability for colored inputs, thus speeding up the convergence.

Remark 1: In general, $\{c_{m,k}\}$ in (9) are determined by one of many static rules (e.g., the Metropolis rule [73] that we adopt in this paper) which keeps them constant during the estimation. Considering that nodes may be working under different signal-to-noise ratios (SNRs), several adaptive rules have been proposed to optimize the algorithm behavior [73, 74, 75]. However, these adaptive rules are severely polluted when impulsive noise samples appear, since the output errors at nodes directly participate in the adaptation of $c_{m,k}$. Designing a robust adaptive rule is an alternative, but it is not the focus of this paper. In another approach, based on the detection of impulsive noise, Ahn et al. proposed a robust variable weighting coefficients dLMS (RVWC-dLMS) algorithm which sets the weighting coefficients to zero at nodes disturbed by impulsive noise [76]. Likewise, the RVWC scheme can be extended to dRLS in a straightforward way, resulting in the RVWC-dRLS algorithm with robustness in impulsive noises111In the literature, the RVWC scheme was presented for more general diffusion strategies (namely, also exchanging information among nodes in adaptation step). However, here we do not consider this case for a fair comparison. Besides, such general diffusion strategies require higher computational complexity and communication load [5]., as can be seen in the simulations later on.

III-B R-dRLS Algorithm

An impulsive noise sample at time instant $i$ might lead the dRLS algorithm to diverge via $e_{k}(i)$ in (5) due to its large amplitude and the propagation of its effect. This degradation effect can last for many iterations. To endow the algorithm with robustness in impulsive noise scenarios, we propose to minimize (3) under the following constraint:

[TABLE]

where $\xi_{k}(i-1)$ is a positive bound. A similar constraint appeared in an adaptive filter for a single agent scenario [59], but when generalizing to the distributed version with multiple agents, the constraint could be imposed on the adaptation at all the nodes. This constraint represents that the energy (squared norm) of the update at every node $k$ from $\bm{w}_{k,i-1}$ to $\bm{\psi}_{k,i}$ always does not exceed the amount $\xi_{k}(i-1)$ regardless of the type of noise (possibly, impulsive noise), thereby guaranteeing the robustness of the algorithm. In doing so, if (5) satisfies (10), i.e.,

[TABLE]

where $\left.\bm{g}_{k,i}\triangleq\bm{P}_{k,i}\bm{u}_{k,i}\right.$ represents the Kalman gain vector, then (5) is a solution of the above constrained minimization problem. Conversely, if (10) is not satisfied (usually in the case of appearance of impulsive noise), i.e., $\left.\lVert\bm{g}_{k,i}\rVert_{2}\lvert e_{k}(i)\rvert>\sqrt{\xi_{k}(i-1)}\right.$, we propose to replace the update (5) by its normalized form to satisfy the equality in (10), which is described by

[TABLE]

where $\text{sign}(\cdot)$ is the sign function. Thus, combining (5), (11) and (12), we obtain the adaptation step for each node $k$ as:

[TABLE]

Evidently, the crucial problem is how to properly choose the bound $\xi_{k}(i)$ as it controls the robustness of the algorithm against impulsive noise and influences its dynamic behavior. To be more specific, we wish $\xi_{k}(i)$ to have larger values at the earlier adaptation stage to provide a fast initial convergence, while for enforcing good robustness against impulsive noise, its values cannot be too large. In addition, we also wish to obtain a small estimation error at steady-state, so $\xi_{k}(i)$ should be reduced to a small value. Based on these requirements, we consider the equality in (10) to propose a useful recursive method for adjusting $\xi_{k}(i)$, as described by

[TABLE]

where $\beta$ is a memory factor with $0<\beta<1$. At every node $k$, $\xi_{k}(i)$ can be initialized by $\xi_{k}(0)=E_{c}\sigma_{d,k}^{2}/(M\sigma_{u,k}^{2})$, where $E_{c}$ is a positive integer, and $\sigma_{d,k}^{2}$ and $\sigma_{u,k}^{2}$ are powers of the output measurement $d_{k}(i)$ and the input regressor $\bm{u}_{k,i}$, respectively. As one can see in (14), every node $k$ not only uses its own adaptive rule to update $\xi_{k}(i)$, but also exploits the side information $\zeta_{m}(i)$ transmitted from its neighboring nodes by the diffusion cooperation. In doing so, the proposed R-dRLS algorithm is more effective at computing consistent estimates at all nodes, which will be observed in Section VI-A. Table I details the proposed R-dRLS algorithm together with the NC method.

Remark 2: As can be seen from (13), the operation mode of the proposed R-dRLS algorithm in the adaptation step can be as follows. At time instant $i$, if $\lVert\bm{g}_{k,i}\rVert_{2}^{2}e_{k}^{2}(i)\leq\xi_{k}(i-1)$, the classical RLS update is performed. If not, the squared norm of the RLS increment is limited to the amount $\xi_{k}(i-1)$ as in (12) for guaranteeing the robustness in impulsive noise. At the early iterations, the values of $\xi_{k}(i-1)$ can be high compared to $\lVert\bm{g}_{k,i}\rVert_{2}^{2}e_{k}^{2}(i)$ so that the algorithm will behave as the dRLS algorithm, providing a fast initial convergence. On the other hand, whenever an impulsive noise sample appears, due to its significant magnitude, the R-dRLS algorithm will work as an dRLS update multiplied by a very small scaling factor $\frac{\sqrt{\xi_{k}(i-1)}}{\lVert\bm{g}_{k,i}\rVert_{2}|e_{k}(i)|}$. It has been shown in [77, 78] that in the adaptation update term, the multiplication of a small scaling factor can reduce the negative influence of impulsive noise on the estimation. Thus, this also indirectly implies that the R-dRLS algorithm has robustness against impulsive noise. Moreover, the robustness is further maintained over the iterations, due to the decreasing property of $\xi_{k}(i)$ given by (14). In addition to this, the diminishing $\xi_{k}(i)$ also leads to a reduction in the steady-state error of the algorithm. To sum up, the R-dRLS algorithm can be considered as an improved dRLS algorithm with a variable ’step-size’ scheme which has an automatic switch between 1 and $\frac{\sqrt{\xi_{k}(i-1)}}{\lVert\bm{g}_{k,i}\rVert_{2}|e_{k}(i)|}$, as can be observed in (13).

III-C NC Method

As a consequence of the diminishing sequence $\{\xi_{k}(i)\}$, the R-dRLS algorithm has poor ability of tracking (i.e., re-convergence of the algorithm) after $\bm{w}^{o}$ undergoes an abrupt change. In order to overcome this problem, inspired by the idea in [45] for the single-agent scenario, we propose here a diffusion-based NC method, as summarized in Table I. The NC method is implemented in the following two steps.

Step 1: A variable $\varDelta_{k}(i)$ at node $k$ is computed once for every $V_{t}$ iterations, to judge whether the unknown vector changed or not. In this step, $\bm{a}_{k,i}^{T}=\mathcal{R}\left(\left[\frac{e_{k}^{2}(i)}{\|\bm{u}_{k,i}\|_{2}^{2}},\frac{e_{k}^{2}(i-1)}{\|\bm{u}_{k,i-1}\|_{2}^{2}},\text{...},\frac{e_{k}^{2}(i-V_{t}+1)}{\|\bm{u}_{k,i-V_{t}+1}\|_{2}^{2}}\right]\right)$ with $\mathcal{R}(\cdot)$ denoting the ascending arrangement for its arguments. With $\bm{e}=[1,...,1,0,...,0]^{T}$ being a vector whose first $V_{t}-V_{d}$ elements set to one, where $V_{d}$ is a positive integer with $V_{d}<V_{t}$, the product $\bm{a}_{k,i}^{T}\bm{e}$ can remove the effect of outliers (e.g., impulsive noise samples) when computing $\varDelta_{k}(i)$. We use a smooth filtering of $\bm{a}_{k,i}^{T}\bm{e}$ to avoid large fluctuations in computing $\varTheta_{\text{new},k}$ (see Table I), as follows:

[TABLE]

where $\tau$, $0<\tau<1$, is a memory factor. Note that, every node $k$ to compute $\varTheta_{\text{new},k}$ also combines the information from its neighboring nodes based on a diffused cooperation; $\varTheta_{\text{old},k}$ stores the value of $\varTheta_{\text{new},k}$ at the last time instant.

From Step 1, one can see that using a larger $V_{t}$, the algorithm has lower steady-state error but a higher delay in tracking. Moreover, for a large occurrence probability of impulsive noise, the value of $V_{d}$ should be increased to better discard the impulsive noise samples in the computation of $\varDelta_{k}(i)$. From our extensive simulations, we found out that for both $V_{t}$ and $V_{d}$, good choices are $V_{t}=\varrho M$ with $1\leq\varrho\leq 3$ and $V_{d}=0.75V_{t}$ [45].

Step 2: If $\varDelta_{k}(i)>t_{\text{th}}$, where $t_{\text{th}}$ is a predefined threshold, it is decided that a change of $\bm{w}^{o}$ has occurred. Then, we reset $\xi_{k}(i)$ to its initial value $\xi_{k}(0)$ so that the R-dRLS algorithm can track this change rapidly. Meanwhile, $\bm{P}_{k,i}$ should also be re-initialized with $\bm{P}_{k,0}$.

It is worth noting that in this scheme the parameters $\tau,\;\varrho$, and $t_{\text{th}}$ are not affected by each other so that their choices are simplified.

IV Mean Square Performance Analyses

IV-A Steady-state Behavior

In this section, we discuss the steady-state behavior of the R-dRLS algorithm in impulsive noise. Assuming that the vector $\bm{w}^{o}$ is invariant, then we define the estimate deviation and intermediate estimate deviation vectors respectively as:

[TABLE]

Using these definitions and (14), it is easy to rearrange (13) and (8), respectively, as:

[TABLE]

and

[TABLE]

Equating the squared $l_{2}$-norm of both sides of (17) and then taking the expectation, we obtain

[TABLE]

Likewise treating (18) and applying Jensen’s inequality [79, p.77], we obtain

[TABLE]

Typically, $\beta$ is close to 1 so that the variances of $\zeta_{k}(i)$ and $\xi_{k}(i)$ given in (14) would be small enough. Accordingly, it can be assumed that

[TABLE]

Then, with this approximation, (19) is changed to

[TABLE]

To deal with the (a) term in (22), some assumptions are helpful.

Assumption 1: The input regressors $\bm{u}_{k,i}$ are zero-mean with covariance matrices $\bm{R}_{k}=E\{\bm{u}_{k,i}\bm{u}_{k,i}^{T}\}$ and spatially independent.

Assumption 2: The regressors $\{\bm{u}_{k,i}\}$ are independent of the estimates $\{\bm{w}_{m,j}\}$ for $j\leq i$ and all $k,m$, referred to as the independence assumption, which is known as useful in the analysis of adaptive algorithms [67] and distributed estimation algorithms [3, 80].

Assumption 3: There is an iteration number $i_{0}$ such that for all $i>i_{0}$, the time-averaged matrix $\bm{\Phi}_{k,i}$ at every node $k$ can be replaced by its expected value $E\left\{\bm{\Phi}_{k,i}\right\}$. This is an ergodicity assumption since $0\ll\lambda<1$, and from (4) we have

[TABLE]

Correspondingly, we can also replace the random matrix $\bm{\Phi}_{k,i}^{-1}$ by $\bar{\bm{\Phi}}_{k}^{-1}\triangleq E\left\{\bm{\Phi}_{k,i}^{-1}\right\}$ for a sufficiently large number of iterations $i$. Note that such replacements are commonly used in the performance analysis of RLS-type algorithms, see [55, 67, 81, 56] and the references therein.

Applying assumption 3, we are able to represent the term (a) in (22) as:

[TABLE]

In the light of assumption 1, if the dimension of $\bm{w}^{o}$ is large, i.e., $M\gg 1$, the fluctuation of the denominator term in (24) from one iteration to the next can be assumed to be small. So, we could make the following approximation (which is also verified in Appendix A):

[TABLE]

where

[TABLE]

Considering the presence of impulsive noise, we need the following assumptions to continue the analysis.

Assumption 4: At every node $k$, the additive noise $v_{k}(i)$ is drawn from a CG random process, $v_{k}(i)=\theta_{k}(i)+\eta_{k}(i)$, where $\theta_{k}(i)$ is the background noise assumed to be zero-mean white Gaussian with variance $\sigma_{\theta,k}^{2}$. The impulsive part $\eta_{k}(i)$ is described as $\eta_{k}(i)=b_{k}(i)g_{k}(i)$, where $b_{k}(i)$ is drawn from a Bernoulli random process with the probability $\left.P[b_{k}(i)=1]=p_{r,k}\right.$, and $g_{k}(i)$ is drawn from a white Gaussian random process with zero-mean and variance $\sigma_{g,k}^{2}=\hbar\sigma_{\theta,k}^{2}$, $\hbar\gg 1$. Usually, $p_{r,k}$ is also called the appearance probability of an impulsive noise sample.

Then, the mean and variance of $v_{k}(i)$ are zero and $\sigma_{v,k}^{2}=p_{r,k}(\hbar+1)\sigma_{\theta,k}^{2}+(1-p_{r,k})\sigma_{\theta,k}^{2}$, respectively. Note that, only when $p_{r,k}=0$ or 1, $v_{k}(i)$ is Gaussian; otherwise, $v_{k}(i)$ is non-Gaussian. Also, $v_{k}(i)$ conditioned on $b_{k}(i)$ is Gaussian [68]. Although the $\alpha$-stable process is more appropriate for modeling impulsive noise in practice [44, 69, 70], one would not consider it in the algorithms’ analysis because its probability density function has no explicit form. Accordingly, the above assumption was used frequently for performance analysis of adaptive algorithms in impulsive noise environments, providing mathematical tractability [51, 45, 82, 68].

Furthermore, as pointed out in [83], when $M\gg 1$, then by using the central limit theorem, it can be assumed that $e_{a,k}(i)$ and $e_{a,k}^{\bm{\Sigma}}(i)$ are zero mean Gaussian variables for any constant matrix $\bm{\Sigma}$. Then, we can employ the following Lemma.

Lemma: Let $e_{a}$ and $u$ be jointly Gaussian zero-mean random variables. Let $e=e_{a}+v$, where $v$ is a zero-mean CG random variable with variance $\sigma_{v}^{2}=p_{r}(\hbar+1)\sigma_{\theta}^{2}+(1-p_{r})\sigma_{\theta}^{2}$, and $v$ is independent of $e_{a}$ and $u$. If $e_{1}=e_{a}+\omega_{1}$ and $e_{2}=e_{a}+\omega_{2}$, where $\omega_{1}$ and $\omega_{2}$ are zero-mean Gaussian random variables with variances $\sigma_{\omega_{1}}^{2}=(\hbar+1)\sigma_{\theta}^{2}$ and $\sigma_{\omega_{2}}^{2}=\sigma_{\theta}^{2}$, and are independent of $u$ and $e_{a}$, then

[TABLE]

Such a Lemma has been commonly used in the past for analyzing the sign-based algorithms [82, 51]. Based on Price’s theorem in [84], Lemma and assumption 2, we can establish the following equation

[TABLE]

where

[TABLE]

and the notation $E\{s|q\}$ accounts for the expectation of $s$ conditioned on $q$. Subsequently, the right-hand term in equality (25) becomes

[TABLE]

where the equality (a) is the result of using (26) under assumption 2.

Substituting (25) and (29) into (22), it is rearranged as

[TABLE]

Next, we introduce the following network global vectors:

[TABLE]

and the network global matrices

[TABLE]

Also, we define the matrix $\bm{C}$ to collect the combination coefficients, i.e., $[\bm{C}]_{m,k}=c_{m,k}$. Following (31) and (32), we can formulate (20) and (30) for all nodes as follows:

[TABLE]

Taking the $\infty$-norm for both sides of (33) leads to

[TABLE]

where the equality (a) uses the fact that $\|\bm{C}^{T}\|_{\infty}=1$ in that the summation of each column of $\bm{C}$ is 1. Since $\bm{\Lambda}_{i}$ and $\bm{\Omega}_{i}$ are diagonal matrices with positive entries, (34) can be equivalently expressed as [12]:

[TABLE]

for $k=1,...,N$. When the algorithm reaches the steady-state, i.e., $E\left\{\lVert\widetilde{\bm{w}}_{k,i}\rVert_{2}^{2}\right\}=E\left\{\lVert\widetilde{\bm{w}}_{k,i-1}\rVert_{2}^{2}\right\}$ as $i\rightarrow\infty$, from (35) we will get:

[TABLE]

In view of the result that $E\{\zeta_{k}(i)\}$ and $E\{\xi_{k}(i)\}$ converge approximately to 0 as $i\rightarrow\infty$ (see Appendix B) as well as $\chi_{k}\neq 0$ and $\varpi_{k,i}\neq 0$, thus, from (LABEL:036) we are able to deduce that

[TABLE]

As a result, (37) illustrates that based on given assumptions, the R-dRLS algorithm can converge to the true parameter vector in the mean-square sense after enough iterations even in impulsive noise environments.

IV-B Analysis of Evolution Behavior

The result (37) is qualitative so that it does not predict the steady-state performance of the algorithm, due mainly to the use of the upper bound relation (20). In this subsection, we will establish a recursive model to describe the evolution behavior of the algorithm in impulsive noise. We start by defining the following network vectors collected from all nodes:

[TABLE]

where

[TABLE]

for nodes $k=1,...,N$. By these defined vectors, we can associate (17) with (18) at all the nodes:

[TABLE]

where $\mathcal{C}=\bm{C}\otimes I_{M}$. Post-multiplying (40) by its transpose and taking the expectation, we have

[TABLE]

where $\bm{W}_{i}\triangleq E\{\widetilde{\bm{w}}_{i}\widetilde{\bm{w}}_{i}^{T}\}$ denotes the covariance matrix of the deviation vector $\widetilde{\bm{w}}_{i}$, and its $k$-th diagonal block of size $M\times M$, i.e., $\bm{W}_{k,i}\triangleq E\{\widetilde{\bm{w}}_{k,i}\widetilde{\bm{w}}_{k,i}^{T}\}$, represents the covariance matrix of the deviation vector $\widetilde{\bm{w}}_{k,i}$ at node $k$.

To evaluate terms I-III in (41), in addition to the spatially independence in assumption 1, we also require the input regressors $\bm{u}_{k,i}$ to be statistically independent in time, which is also often used in analysis of distributed estimation algorithms [3, 2]. Therefore, performing similar manipulations as in Section IV-A on the expectations under assumptions 2-4, Lemma and Price’s theorem, we can compute these three terms. Specifically, the term I in (41) becomes

[TABLE]

where we rewrite $\varpi_{k,i}$ contained in $\bm{\Lambda}_{i}$ as

[TABLE]

The term II in (41) is the transpose of (43). For any $k$ and $m$ belonging to the set $\{1,...,N\}$, we define the $(m,k)$-th $M\times M$ matrix $E\{\bm{\varXi}_{i}\bm{\varXi}_{i}^{T}\}$ as follows:

[TABLE]

When $k=m$, (44) represents the $k$-th diagonal block of $E\{\bm{\varXi}_{i}\bm{\varXi}_{i}^{T}\}$, which is described as

[TABLE]

where $\Omega_{k,i}$ is the $k$-th element of $\bm{\Omega}_{i}$. When $k\neq m$, the off-diagonal blocks will be simplified as

[TABLE]

From (45) and (LABEL:046), we obtain the term III in (41):

[TABLE]

where

[TABLE]

By substituting (42) and (47) into (41), we obtain the recursive expression for $\bm{W}_{i}$:

[TABLE]

The mean square deviation (MSD) at node $k$ is defined as $\text{MSD}_{k}(i)\triangleq\text{Tr}\{\bm{W}_{k,i}\}$, and the network MSD over all the nodes is defined as $\text{MSD}_{\text{net}}(i)=\frac{1}{N}\sum_{k=1}^{N}\text{MSD}_{k}(i)=\text{Tr}\{\bm{W}_{i}\}/N$ [2]. Equation (49) models the MSD evolution behavior of the algorithm. It needs to be mentioned that to implement the model (49), $E\{\xi_{k}\}$ and $E\{\zeta_{k}\}$ defined in $\bm{\Omega}_{i}$ still need to be evaluated further. However, as shown in (14), $\xi_{k}$ and $\zeta_{k}$ between interconnected nodes are affected by each other and there is a comparison operation, so it is difficult to provide an evolution expression for them. In this paper, we suggest that $E\{\xi_{k}\}$ and $E\{\zeta_{k}\}$ are obtained by the ensemble average using simulations. Consequently, although (49) is a semi-analytic result, it can also be used to evaluate the convergence of the proposed algorithm.

V DCD-Based Algorithms

In this section, we review the DCD-dRLS algorithm from [66], and then develop a robust DCD-dRLS algorithm.

V-A The Original DCD-dRLS Algorithm

Since the dRLS algorithm involves the matrix operation of size $M\times M$ in the computations of $\bm{g}_{k,i}$ and (7) at every node, it requires a computational complexity that scales as a quadratic function of $M$ in terms of additions and multiplications per iteration $i$. To reduce the complexity, the adaptation step of the DCD-dRLS algorithm is described as [66]:

[TABLE]

where the increment $\Delta\bm{w}_{k,i}$ is obtained by solving the normal equation:

[TABLE]

[TABLE]

$\bm{r}_{k,i-1}$ defines the residual vector at node $k$ at time instant $i-1$:

[TABLE]

For reducing the complexity of computing $\Delta\hat{\bm{w}}_{k,i}$ and $\bm{r}_{k,i}$, the DCD method presented in Table II is used; see [63, 62, 64] for details. In Table II, $[\bm{r}_{k,i}]_{l}$ is the $l$-th entry of a vector $\bm{r}_{k,i}$, and $[\bm{\Phi}_{k,i}]_{l,l}$ and $[\bm{\Phi}_{k,i}]_{:,l}$ are the $(l,l)$-th entry and the $l$-th column of $\bm{\Phi}_{k,i}$, respectively.

The accuracy and complexity of the DCD method are dependent on three parameters: $H$, $M_{b}$, and $N_{u}$. In general, $H$ is chosen as a power-of-two number; $M_{b}$ is the number of bits being enough for a fixed-point representation of $\Delta\hat{\bm{w}}_{k,i}$ within an amplitude range $[-H,H]$; and $N_{u}$ defines a maximum number of elements in $\Delta\hat{\bm{w}}_{k,i}$ that can be updated at a time instant. The DCD method only requires $2N_{u}M+M_{b}$ additions at most at each time instant with no multiplication [62]. Also, a larger $N_{u}$ makes the solution $\Delta\hat{\bm{w}}_{k,i}$ closer to the optimal solution $\Delta\bm{w}_{k,i}$ in (51), but increases the number of additions. It follows that if $N_{u}<M$, the DCD-based algorithm implements a selective partial update [85].

Similar to the dRLS algorithm, however, the DCD-dRLS algorithm will also encounter the performance deterioration when impulsive noise happens.

V-B Proposed DCD-R-dRLS Algorithm

To achieve robustness against impulsive noise, we present here the DCD-R-dRLS algorithm.

Step 1: At every node $k$, we firstly use the DCD method to solve the normal equation (51) with (4) and (52), yielding a solution $\Delta\hat{\bm{w}}_{k,i}^{(1)}$ and residual vector $\bm{r}_{k,i}^{(1)}$. In the presence of impulsive noise, we also impose a constraint similar to that in (10):

[TABLE]

Step 2: If $\lVert\Delta\hat{\bm{w}}_{k,i}^{(1)}\rVert_{2}^{2}\leq\xi_{k}(i-1)$, we set $\Delta\hat{\bm{w}}_{k,i}=\Delta\hat{\bm{w}}_{k,i}^{(1)}$ and $\bm{r}_{k,i}=\bm{r}_{k,i}^{(1)}$ and then perform the update (50). Otherwise, we need to recalculate $\bm{b}_{k,i}$ in (51) as:

[TABLE]

Subsequently, based on the DCD method, we obtain the solution $\Delta\hat{\bm{w}}_{k,i}^{(2)}$ and the residual vector $\bm{r}_{k,i}^{(2)}$ from the normal equation (51) under (4) and (55), thereby performing the update (50) with the increment

[TABLE]

and $\bm{r}_{k,i}=\bm{r}_{k,i}^{(2)}$.

Step 3: The combination step (12) is performed.

Step 4: The bound parameter $\xi_{k}(i)$ in the DCD-R-dRLS algorithm is updated according to

[TABLE]

Table III summarizes the DCD-R-dRLS algorithm.

Remark 3: An impulsive noise sample appearing at time instant $i$ would yield a mismatch solution $\Delta\hat{\bm{w}}_{k,i}^{(1)}$ so that $\left.\lVert\Delta\hat{\bm{w}}_{k,i}^{(1)}\rVert_{2}^{2}>\xi_{k}(i-1)\right.$. In this case, the scaling factor $\frac{\sqrt{\xi_{k}(i-1)}}{\lVert\Delta\hat{\bm{w}}_{k,i}^{(1)}\rVert_{2}}$ in (55) is small enough to eliminate impulsive noise hidden in $e_{k}(i)$. A similar scaling factor $\frac{\sqrt{\xi_{k}(i-1)}}{\lVert\Delta\hat{\bm{w}}_{k,i}^{(2)}\rVert_{2}}$ in (56) is to make the increment satisfy the constraint (54). Consequently, the DCD-R-dRLS algorithm improves the robustness to impulsive noise relative to the DCD-dRLS algorithm. Moreover, the decreasing sequence $\{\xi_{k}(i)\}$ shown in (57) further guarantees the robustness. It is worth noting that due to $\rVert\bm{g}_{k,i}\lVert_{2}\lvert e_{k}(i)\rvert\approx\lVert\Delta\hat{\bm{w}}_{k,i}\rVert_{2}$, the DCD-R-dRLS algorithm is a DCD-based variant of the R-dRLS algorithm. Unlike the R-dRLS algorithm, based on the NC method we re-initialize $\bm{\Phi}_{k,i}$ with $\bm{\Phi}_{k,0}$ to endow the DCD-R-dRLS algorithm with the tracking capability when $\bm{w}^{o}$ suddenly changes.

Remark 4: Let $C_{dcd}^{+}$ denote the only required number of additions for the DCD algorithm, with $C_{dcd}^{+}\leq 2N_{u}M+M_{b}$. In Table IV, we provide the computational complexity of the existing dLMS, dRLS, DCD-dRLS, and both proposed R-dRLS and DCD-R-dRLS algorithms at node k per time instant $i$, where $n_{k}$ denotes the cardinality of $\mathcal{N}_{k}$. For shift structured input regressor at node $k$ [54, 9], i.e., $\bm{u}_{k,i}=[u_{k}(i),u_{k}(i-1),...,u_{k}(i-M+1)]^{T}$, where $u_{k}(i)$ is an input sample at time instant $i$, implementing $\bm{\Phi}_{k,i}$ in (4) is very simplified. In this situation, by copying the upper-left $(M-1)\times(M-1)$ block of $\bm{\Phi}_{k,i-1}$ leads to the lower-right $(M-1)\times(M-1)$ block of $\bm{\Phi}_{k,i}$. The remaining part of $\bm{\Phi}_{k,i}$ that needs to be updated is the first row and first column. Owing to symmetry of $\bm{\Phi}_{k,i}$, only calculating the first column is sufficient, which is formulated as:

[TABLE]

Note that, in the DCD-R-dRLS algorithm, $\left.\kappa=1\right.$ represents the case $\left.\lVert\Delta\hat{\bm{w}}_{k,i}\rVert_{2}^{2}>\xi_{k}(i-1)\right.$ at time instant $i$ (which leads to the maximum complexity), otherwise $\kappa=0$. The comparisons required in the R-dRLS and DCD-R-dRLS algorithms are counted as additions.

Consider an example with $n_{k}=10$, $M_{b}=32$ and $\kappa=1$, Fig. 2 depicts the number of operations of some diffusion algorithms in terms of multiplications and additions at node $k$ at each time instant versus $M$. It is clear that the computational complexity of the dLMS algorithm, with the order of $\mathcal{O}(M)$, is much lower than that of the dRLS algorithm. As expected, since $N_{u}<M$, compared with the standard dRLS and R-dRLS algorithms, their DCD versions obtain about $50\%$ reduction in both multiplications and additions for the case of general input regressors. However, for shift structured input regressors, the computational cost is drastically reduced from the order $\mathcal{O}(M^{2})$ to $\mathcal{O}(M)$, which is more pronounced in scenarios with large $M$. Moreover, the multiplications required in the DCD-based algorithms are not dependent of $N_{u}$. On the other hand, in contrast with the existing dRLS and DCD-dRLS algorithms, the additional complexities in the proposed R-dRLS and DCD-R-dRLS algorithms resulted from the computations of the scaling factor and the bound parameter are small. In addition to the complexity, for both proposed algorithms, each node $k$ increases communication cost of $n_{k}-1$ numbers for transmitting $\zeta_{k}$ to its neighbors.

Remark 5: From the DCD-R-dRLS algorithm, we can directly obtain its special form for a single-agent scenario, referred it to as the DCD-R-RLS algorithm. In other words, the DCD-R-RLS algorithm is the DCD implementation of the algorithm presented in [59].

VI Simulation Results

Simulation examples are presented for a diffusion network with $N=20$ nodes on distributed parameter estimation and distributed spectrum estimation. The network topology adopted for all simulations is shown in Fig. 3(a), unless otherwise specified. Herein, we do not consider the measurement sharing in the adaptation step for all diffusion algorithms. The Metropolis rule [73] used for computing the combination coefficients $\{c_{m,k}\}$ in combination step is expressed as:

[TABLE]

VI-A Distributed Parameter Estimation

The vector $\bm{w}^{o}$ to be estimated has a length of $M=16$ and a unit norm; it is generated randomly from a zero-mean uniform distribution. The input regressor $\bm{u}_{k,i}$ has a shift structure, where $u_{k}(i)$ is colored and generated by a second-order autoregressive system:

[TABLE]

where $\epsilon_{k}(i)$ is a zero-mean white Gaussian process with variance $\sigma_{\epsilon,k}^{2}$. The background noise $\theta_{k}(i)$ is zero-mean white Gaussian noise with variance $\sigma_{\theta,k}^{2}$. Variances $\sigma_{\epsilon,k}^{2}$ and $\sigma_{\theta,k}^{2}$ are shown in Fig. 3(b) and (c), respectively, for all the nodes. We employ the network MSD to assess the performance of algorithms. All results are the average over 200 independent trials.

Example 1: Except for the background noise $\theta_{k}(i)$, a cluster of impulses with length $200$ is also added to corrupt $d_{k}(i)$ at iteration $i=5001$ 222Such a scenario is similar to double-talk in echo cancellation.. The cluster is drawn from a zero-mean white Gaussian process, but with a large variance $1000\sigma_{y,k}^{2}$ to generate impulsive samples, where $\sigma_{y,k}^{2}$ denotes the power of $y_{k}(i)=\bm{u}_{k,i}^{T}\bm{w}^{o}$. Fig. 4 compares the performance of the proposed R-dRLS algorithm with that of the dRLS and both LTVFF-dRLS and LCTVFF-dRLS algorithms presented in [81]. The parameters of the algorithms are set to make a comparable convergence rate. The regularization constant for all RLS-type algorithms is chosen as $\delta=0.01$. It is clear to see, for a small forgetting factor $\lambda=0.98$, the conventional dRLS algorithm converges faster but has a higher estimation error; conversely, by increasing the forgetting factor, it has a lower estimation error but its convergence rate becomes slower. In particular, using a large forgetting factor $\lambda=0.998$, the dRLS will need more time to converge again after a cluster of impulses enforces the algorithm to diverge. Due to the use of variable forgetting factor schemes, both LTVFF-dRLS and LCTVFF-dRLS algorithms solve this performance trade-off to a certain extent. As stated in Remark 2, the R-dRLS algorithm also overcomes this performance trade-off since it employs a variable ’step-size’ factor in the adaptation step. Besides, unlike the dRLS, LTVFF-dRLS and LCTVFF-dRLS algorithms, even though a cluster of impulses does not happen until the algorithms reach the steady-state, the R-dRLS algorithm also does not undergo divergence. This is because the R-dRLS algorithm can judge by (11) whether impulses occur or not and perform corresponding updates.

Example 2: The additive noise $v_{k}(i)$ is a CG process given in assumption 4. At every node $k$, we set $p_{r,k}$ as a random number in the range of $[0.001,0.05]$ and $\sigma_{g,k}^{2}=1000\sigma_{y,k}^{2}$. For a fair comparison of RLS-type algorithms, we choose the same forgetting factor $\lambda$=0.985 and regularization constant $\delta$=0.01, except $\delta$=0.5 in the dRLP and RVWC-dRLS algorithms.

Fig. 5 checks the validity of the semi-analytic result (49), where we plot $E\{\xi_{1}(i)\}$ at node 1 (having similar results at other nodes). To take into account the assumption on input regressors $\bm{u}_{k,i}$ in analysis, here its entries are generated from a white Gaussian process $\epsilon_{k}(i)$. To compute (49), we use the same impulsive noise parameters: $p_{r,k}=0.01$ or 0.05, and $\sigma_{g,k}^{2}=10000\sigma_{\theta,k}^{2}$ at all the nodes. As one can see, the theoretical results have good fit with the simulated results. Moreover, $E\{\xi_{1}(i)\}$ obtained by the ensemble average of simulations is a decreasing function of the iteration $i$, which further supports the theory in Appendix B.

Fig. 6 investigates the effect of the NC method on the R-dRLS algorithm. It can be seen that the R-dRLS algorithm will not re-converge after $\bm{w}^{o}$ changes to $-\bm{w}^{o}$ at iteration $i=2501$. In this scenario, all algorithms have a large sharp phase transition of MSD due to the mismatch between $-\bm{w}^{o}$ and its estimate $\bm{w}_{k,i}$ at that moment. The NC method can endow the R-dRLS algorithm with good tracking capability for such a change of $\bm{w}^{o}$. Benefited from the smoothing operation (15), the NC ($\left.\tau=0.96\right.$) only slightly degrades the steady-state performance of the R-dRLS algorithm compared with the non-smooth version in [1] (i.e., $\left.\tau=0\right.$).

In Fig. 7, we compare the performance of the dRLS, dSE-LMS, dLMP, RVWC-dRLS, and dRLP algorithms with that of the proposed R-dRLS with NC algorithm. Note that, the R-dRLS (no cooperation) is that each node performs a standalone adaptive algorithm presented in [59]. As expected, the dRLS algorithm has a poor performance in the presence of impulsive noise, while other algorithms are robust. Among these robust algorithms, the convergence of dSE-LMS and dLMP algorithms is slow. Thanks to the decorrelation property of dRLS, the RVWC-dRLS, dRLP, and R-dRLS with NC algorithms obtain fast convergence. In particular, the proposed R-dRLS with NC algorithm has also a large reduction in the steady-state MSD. This is due mainly to the fact that its updated energy described by (10) and (14) diminishes with iterations.

Example 3: The additive noise $v_{k}(i)$ here is generated by the $\alpha$-stable process, also called the $\alpha$-stable noise. Its characteristic function is given by $\varphi(t)=\exp(-\gamma\lvert t\lvert^{\alpha})$ [72, 69], where the characteristic exponent $\alpha\in(0,2]$ describes the impulsiveness of the noise (smaller $\alpha$ leads to more impulsive noise samples) and $\gamma>0$ represents the dispersion level of the noise. In particular, when $\alpha=2$, it reduces to the Gaussian noise. It is rare to find $\alpha$-stable noise with $\alpha<1$ in practice [72, 69]. In this example, thus we set $\alpha=1.2$ and $\gamma=2/15$. The learning performance of the algorithms is shown in Fig. 8. Fig. 9 shows the node-wise steady-state MSD of the robust algorithms (i.e., excluding the dRLS), by averaging MSD values from iteration 2 400 to 2500. As can be seen from Figs. 8 and 9, the proposed R-dRLS algorithm with NC outperforms the known robust diffusion algorithms in terms of convergence rate, steady-state accuracy and tracking capability. As shown in Fig. 7 to Fig. 9, due to the cooperation of interconnected nodes, the R-dRLS algorithm improves the estimation performance compared with its non-cooperative counterpart.

We also perform the simulations for the network in Fig. 10 with less connections among nodes. Fig. 11 shows the node-wise steady-state MSD of those algorithms in Fig. 9. By comparing these two figures, it is seen that the proposed R-dRLS algorithm is more likely to reach the same estimates at all nodes.

Example 4: Comparison of DCD-algorithms. Figs. 12 and 13 compare the DCD-R-dRLS algorithm using different $N_{u}$ values with its standard version in CG-noise and $\alpha$-noise scenarios333Here the curves of both the R-dRLS and DCD-dRLS algorithms are omitted due to their divergence performance in impulsive noise.. The DCD parameters are $H=4$ and $M_{b}=16$. It is seen that, the proposed DCD-R-dRLS algorithm is also robust to impulsive noises, and approaches the R-dRLS performance with increase in $N_{u}$. In this example, the DCD-R-dRLS algorithm with $N_{u}=4\;(<M)$ has a good approximation to the R-dRLS algorithm, while the complexity of the former is significantly lower than that of the latter. Moreover, many simulations have been carried out in different impulsive noise scenarios by prolonging the iteration $i$ to a larger number than the one in Fig. 5, e.g., $5\times 10^{5}$, using MATLAB R2013A on a Intel(R) Core(TM) i5-4590 CPU @ 3.30 GHz processor. We did not observe any numerical instability during the simulations for both proposed R-dRLS and DCD-dRLS algorithms.

VI-B Application: Distributed Spectrum Estimation

We have also tested the proposed algorithms’ performance in an application of distributed spectrum estimation in CR, in which the objective is to estimate the spectrum of a transmitted signal source $s$ in the network with $N$ nodes [7, 8, 81]. We use $\bm{\phi}_{s}(f)=\sum_{m=1}^{M}q_{m}(f)w_{m}^{o}=\bm{q}^{T}(f)\bm{w}^{o}$ to denote the power spectral density (PSD) of the signal $s$ at frequency $f$, where $\bm{q}(f)=[q_{1}(f),...,q_{M}(f)]^{T}$ is a vector consisting of basis functions evaluated at normalized frequency $f$, and $\bm{w}^{o}=[w_{1}^{o},...,w_{M}^{o}]^{T}$ stands for the power that transmits the signal $s$ over each of $M$ basis functions and needs to be estimated. Such basis expansion can accurately model the spectrum of the signal $s$ for large enough $M$. Considering $\mathcal{H}_{k}(f,i)$ is the transfer function of the channel between the station emitting the signal $s$ and receiver node $k$ at time instant $i$, the PSD of the received signal at node $k$ can be expressed as

[TABLE]

where $\bm{q}_{k,i}(f)=|\mathcal{H}_{k}(f,i)|\bm{q}(f)$, and $\sigma_{r,k}^{2}$ is the received noise power at node $k$.

At time instant $i$, each node $k$ observes the received PSD expressed in (58) over $N_{c}$ frequency samples $f_{\iota}=f_{\min}:(f_{\max}-f_{\min})/N_{c}:f_{\max}$ for $\iota=1,...,N_{c}$; accordingly, the output measurements of node obey the following relation:

[TABLE]

where $v_{k}^{\iota}(i)$ denotes the observation noise at frequency $f_{\iota}$. The noise power $\sigma_{r,k}^{2}$ can be estimated with high accuracy before the spectrum estimation, using, for example, an energy estimator over an idle band, and then subtracted from (59) [7, 8, 81]. Then, by collecting the output measurements over $N_{c}$ frequencies, we obtain a data model at every node $k$ for distributed spectrum estimation:

[TABLE]

where $\bm{Q}_{k,i}=[\bm{q}_{k,i}(f_{1}),...,\bm{q}_{k,i}(f_{N_{c}})]^{T}$, $\bm{d}_{k}(i)=[d_{k}^{1}(i),...,d_{k}^{N_{c}}(i)]^{T}$, and $\bm{v}_{k}(i)=[v_{k}^{1}(i),...,v_{k}^{N_{c}}(i)]^{T}$.

Based on this model, we estimate the unknown spectrum $\bm{w}^{o}$ of the signal $s$ using different diffusion algorithms over the network given in Fig. 3(a). In the simulation [8, 81], we use $M=50$ nonoverlapping rectangular basis functions444Other basis functions are also possible, e.g., raised cosines, or Gaussian bells [7]. with amplitude equal to one to model the PSD of the signal $s$. The nodes scan $N_{c}=100$ frequencies over the normalized frequency axis between 0 and 1. We assume that $\bm{w}^{o}$ has only 8 non-zero elements, meaning that the unknown spectrum is transmitted over 8 basis functions, and the power transmitted over each basis function is set to 0.7. The observation noise $v_{k}^{\iota}(i)$ is an $\alpha$-stable process as in the previous Example 3 [48]. In Fig. 14, we compare the network MSD performance of different algorithms considered for the distributed spectrum estimation. As depicted, the dRLS algorithm can not identify the spectrum coefficients $\bm{w}^{o}$ due to its divergence in an $\alpha$-stable noise environment. In comparison with the dSE-LMS, dLMP, RVWC-dRLS and dRLP algorithms, the proposed R-dRLS and DCD-R-dRLS (with $N_{u}=4$) algorithms still obtain better estimation performance. We also notice from this figure that the DCD-R-dRLS algorithm with lower computational complexity approaches the R-dRLS performance. In Fig. 15, we also select the robust dRLS-type algorithms to show their performance in terms of PSD at node $1$. From the results, the proposed R-dRLS and DCD-R-dRLS algorithms have lower side lobes in the PSD curves than those of the other two algorithms, thus fitting much better the true spectrum.

VII Conclusion

In this paper, we have derived a new dRLS algorithm which is robust in impulsive noise, based on the minimization of a local RLS cost function with a time-dependent constraint on the squared norm of the intermediate estimate update. Following the diffusion strategy, the constraint is dynamically adjusted with the help of side information from the neighboring nodes. We also analyze the convergence of the proposed algorithm in the mean square sense under impulsive noise. Then, its DCD version was developed to reduce the computational complexity. Moreover, to adapt the proposed algorithms to an abrupt change of the unknown parameter vector, a non-stationary control approach has also been designed. Simulation results have verified that the proposed algorithms perform better than the known algorithms in impulsive noise scenarios.

Appendix A Verification of (25)

From Fig. 16, one can see that the left side of (25) has a good agreement with the right side of that555 Similar results at other nodes have not been shown here because of the page limitation.. This reveals that the simplification from (24) to (25) is reasonable.

Appendix B Convergence of $E\{\xi_{k}(i)\}$ to 0

It is evident from (14) that $\xi_{k}(i)$ as a function of $i$ is non-increasing in adaptation process, with positive values. So, the limit of $E\{\xi_{k}(i)\}$ at $i\rightarrow\infty$ is existent. Applying the expectation operator to (14), we obtain

[TABLE]

[TABLE]

Again using the assumption that the variance of $\xi_{k}(i)$ is small enough since $\beta$ closes to 1, we are able to make the approximation,

[TABLE]

where $m_{k}\doteq\lVert\bm{g}_{k,i}\rVert_{2}^{2}e_{k}^{2}(i)$ means that both $m_{k}$ and $\lVert\bm{g}_{k,i}\rVert_{2}^{2}e_{k}^{2}(i)$ have the same distribution, $P_{k,i}[\cdot]$ denotes the probability of event in the argument, and $F_{k,i}(m_{k})$ denotes the distribution function of $m_{k}$ at time instant $i$.

Let us define the network global vectors as follows:

[TABLE]

Therefore, according to (LABEL:A40) and (B.4), we reformulate (B.1) and (B.2) for all nodes as:

[TABLE]

where

[TABLE]

and

[TABLE]

Taking the $\infty$-norm of both sides of (B.5) and recalling $\lVert\bm{C}^{T}\rVert_{\infty}=1$, it is found the following inequality:

[TABLE]

Based on the diagonal definition in (B.6), we deduce an equivalent form from (B.8), i.e., for $k=1,...,N$,

[TABLE]

It is supposed that there is a limit for $F_{k,i}(m_{k})$ when $i\rightarrow\infty$, (B.9) further reduces to

[TABLE]

In (B.10), the relation $E\{\xi_{k}(i)\}=E\{\xi_{k}(i-1)\}$ as $i\rightarrow\infty$ is also used. Herein, we consider the equality case in (B.10), i.e.,

[TABLE]

It is shown in Appendix A in [45], for a similar equation (B.11), its solution is $E\{\xi_{k}(\infty)\}]=0$. Since (B.11) is an upper bound of (B.10), we can conclude that $E\{\xi_{k}(i)\}$ given by (14) would also converge to zero. Moreover, as its intermediate quantity, $E\{\zeta_{k}(i)\}$ also converges to zero.