Second-Order Nonlinearity Estimated and Compensated Diffusion LMS Algorithm: Theoretical Upper Bound, Cramer-Rao Lower bound, and Convergence Analysis

The problem of distributed estimation of an unknown vector from linear measurements is a well-known subject in signal processing community which has numerous applications in wireless sensor network (WSN), channel estimation, spectrum estimation, massive MIMO communication, and target tracking problems [1], [2]. The distributed estimation algorithms benefit from the inter-collaboration of sensor nodes. The cooperation strategies for distributed estimation are incremental, consensus, and diffusion approaches [1]. Among them, the diffusion strategy is more versatile due to its simplicity, scalability, and low storage demands.

Many distributed diffusion algorithms are proposed in the literature, e.g., diffusion LMS [3], [4], diffusion LMP [5], [6], diffusion Affine Projection Algorithm (APA) [7], [8], diffusion CMPN [9], and diffusion correntropy [10]-[13], to name a few. Among diffusion algorithms, the Diffusion Least Mean Square (DLMS) algorithm is the basic algorithm which uses mean square error (MSE) as its cost function. There are also numerous variants of the DLMS algorithm, which aim to either reduce the communication load [14], [15], make the algorithm robust against impulsive noise [16], [17], make the algorithm secure with respect to adversaries [18], [19], [20], or in sparse setting [21]. Unfortunately, the performance of the aforementioned algorithms in a WSN deteriorates when the sensors have some nonlinearity effect due e.g. to their power amplifiers. This is because they are designed for the linear measurement model. The main objective of this paper is to make the DLMS algorithm robust against nonlinearities.

In the literature of distributed estimation, there are some works that consider a nonlinear model for the measurements [22]-[28]. In the pioneering work of [22], two distributed algorithms are suggested for estimation in a nonlinear observation model. Moreover, a diffusion based kernel least mean squares (KLMS) is presented in a nonlinear measurement setup [23]. In addition, a distributed estimation algorithm with nonlinear sensors with one bit measurements are proposed in [24]. Besides, [25] suggests two algorithms for estimating the parameters of nonlinear Hammerstein systems with missing data. A distributed nonlinear parameter estimation algorithm is further developed in unbalanced multi-agent networks [26]. Nonlinear model is partially used in [27], in which a method for distributed solution of robust estimation problems is proposed with equality constraints based on the augmented Lagrangian method. [28] discusses both linear and nonlinear models for secure distributed estimation in the presence of attackers in the network.

In this paper, we deal with second-order non-linear model for sensors. It allows to model a linear system which shows some small degree of nonlinearities. The challenges of the second-order nonlinear model is in adaptiveness in which the nonlinear coefficient may change over time. In the proposed solution, the adaptiveness are taken into account. Thanks to the second-order nonlinear model, we can investigate the nonlinearity effect on the performance of a DLMS algorithm. Hence, an upper bound for the error is calculated in the paper. Also, an improved version of DLMS algorithm is suggested which incorporates nonlinearity estimation and compensation units. Further, the Cramer-Rao bound (CRB) is calculated for the distributed estimation problem in the presence of second-order nonlinearity. Simulation results show the benefit of the proposed method especially when there are small nonlinear coefficients.

Consider a WSN with $N$ sensors (nodes) collecting a scalar measurement $d_{k,i}$, where $1\leq k\leq N$ is the node index and $1\leq i\leq I$ is the time instant. Each sensor contains its own $L\times 1$ input regression vector $\mathbf{u}_{k,i}$. The model of measurements is linear, i.e., $d_{k,i}={\mathbf{u}^{T}_{k,i}}\mathbf{\omega}_{o}+v_{k,i}$, with the unknown $L\times 1$ vector $\mathbf{\omega}_{o}$, where $v_{k,i}$ denotes the measurement noise. It is assumed that the sensor equipment has a power amplifier with a second-order nonlinear model. So, if the nonlinear function is $f(\cdot)$ then nonlinear measurements are

where $b_{k}$ is the second-order nonlinearity coefficient of $k$th sensor and $\theta_{k,i}$ is the measurement noise which is assumed to be zero-mean Gaussian with variance $\sigma^{2}_{\theta,k}$. The constant term in the above model (1) is omitted since the system is an approximately linear system with a second-order nonlinear term.

The main objective of the distributed estimation problem in the WSN is to estimate the unknown vector $\mathbf{\omega}_{o}$ using nonlinear measurements $\tilde{d}_{k,i}$ and regression vectors $\mathbf{u}_{k,i}$ of sensors. The other objective is to estimate the nonlinearity of the sensors in the network.

Diffusion algorithms are usually suitable solutions for distributed estimation problems, of which DLMS is the most basic. The two steps of the DLMS algorithm is the adaptation and combination steps. It can be implemented in two ways: Adapt Then Combine (ATC) and Combine Then Adapt (CTA). The ATC version of DLMS is as follows [1]:

where $\cal N_{\mathrm{k}}$, $\tilde{\mathbf{\varphi}}_{k,i}$, and $\tilde{\tilde{\mathbf{\varphi}}}_{l,i}=f(\tilde{\mathbf{\varphi}}_{l,i})$ denote the neighborhood set of the $k$’th sensor, the intermediate estimation of $k$’th sensor in the presence of nonlinearity at time index $i$, and the received intermediate estimation of sensor $l$, respectively. Further, $a_{lk}$ and $c_{lk}$ are the combination coefficients from node $l$ to node $k$ in the adaptation and combination steps, respectively. The local cost function of node $k$ in the DLMS algorithm in presence of nonlinearity is defined as

where expectation operator $\mathrm{E}\{\}$ can be neglected and the point estimation can be replaced for expectation.

Diffusion algorithms are usually suitable solutions for distributed estimation problems, of which DLMS is the most basic. The two steps of the DLMS algorithm is the adaptation and combination steps. It can be implemented in two ways: Adapt Then Combine (ATC) and Combine Then Adapt (CTA). The ATC version of DLMS is as follows [1]:

where $\cal N_{\mathrm{k}}$, $\tilde{\mathbf{\varphi}}_{k,i}$, and $\tilde{\tilde{\mathbf{\varphi}}}_{l,i}=f(\tilde{\mathbf{\varphi}}_{l,i})$ denote the neighborhood set of the $k$’th sensor, the intermediate estimation of $k$’th sensor in the presence of nonlinearity at time index $i$, and the received intermediate estimation of sensor $l$, respectively. Further, $a_{lk}$ and $c_{lk}$ are the combination coefficients from node $l$ to node $k$ in the adaptation and combination steps, respectively. The local cost function of node $k$ in the DLMS algorithm in presence of nonlinearity is defined as

where expectation operator $\mathrm{E}\{\}$ can be neglected and the point estimation can be replaced for expectation.

In this section, we investigate the effect of second-order nonlinearity on the performance of the DLMS algorithm. We derive an upper bound for the $l^{2}$-norm of the error term, which is the difference between the estimated vector after combination step in the presence of nonlinearity and without the nonlinearity. To that end, we write the formula of the intermediate estimation in (4) in the following form

where the noise term $\theta_{l}$ is neglected in comparison to nonlinear term $b_{l}d^{2}_{l,i}$ and $\Delta\mathbf{\varphi}_{k,i}\triangleq\mu_{k}\sum_{l\in\cal N_{\mathrm{k}}}c_{lk}(b_{l}\tilde{d}^{2}_{l,i})\mathbf{u}_{l,i}$ is the error vector of the intermediate estimation. In fact, we assume that the noise level is much lower than the nonlinearity term. Then, the intermediate estimations of $\tilde{\mathbf{\varphi}}_{l,i}$ is received by the node $k$ as $\tilde{\tilde{\mathbf{\varphi}}}_{l,i}=f(\tilde{\mathbf{\varphi}}_{k=l,i})+\eta_{l,i}$. So, we have

where $\mathbf{x}^{2}\triangleq\mathbf{x}\odot\mathbf{x}$ is the element-wise square of a vector in which $\odot$ is the hadamard operator. Then, the output of combination unit will be

where $\omega_{k,i}=\sum_{l\in\cal N_{\mathrm{k}}}a_{lk}\mathbf{\varphi}_{l,i}$ is the new true estimation without nonlinearity. To upper bound the $l^{2}$-norm of the error vector of $\tilde{\omega}_{k,i}-\omega_{k,i}$, we write

The upper bound for the error term is expressed in the following theorem.

where

where $b_{l,\mathrm{max}}$ is the upper bound of $|b_{l}|$.

In this section, the mean convergence analysis of DLMS in presence of second-order nonlinearity is performed. In the first case, we assume the nonlinearity both in measurements and links. In the second case, we consider the nonlinearity just in measurements.

DLMS algorithm performance is degraded by nonlinearity. To improve the performance of the DLMS algorithm in the presence of nonlinearity, we propose to estimate the nonlinear coefficients and then compensate their effects. This process is adaptive and online as it shows its benefits in comparison to a pre-calibration process. We call the proposed algorithm second-order nonlinearity estimated and compensated DLMS (SONEC-DLMS) algorithm. It consists of five different steps. The details of steps are:

1. 1.

   Nonlinear coefficient estimation: If we define the coefficient vector of node $k$ as $\mathbf{b}_{k}=[b_{l}]_{l\in\cal N_{\mathrm{k}}}$, we estimate this vector at node $k$ by suggesting the following cost function as:

   $$J_{k}(\mathbf{\omega})=\sum_{l\in\cal N_{\mathrm{k}}}a_{lk}\mathrm{E}\{||\tilde{d}_{l,i}-b_{l}d^{2}_{l,i}-\mathbf{u}^{T}_{k,i}\mathbf{\omega}||^{2}_{2}\}=$$

   $$\bar{J}_{k}(\mathbf{\omega},\mathbf{b}_{k},\mathbf{d}_{k}),$$

   (34)

   where $\mathbf{d}_{k}=[d_{l,i}]_{l\in\cal N_{\mathrm{k}}}$ is the linear ground truth measurement vector of node $k$ with incorporating the measurement noise and is estimated in the second step of the algorithm, and $\mathbf{b}_{k}=[b_{l,k}]_{l\in\cal N_{\mathrm{k}}}$ is the nonlinear coefficients estimated by node $k$. To find $\mathbf{b}_{k}$, we use a steepest descent of $\bar{J}_{k}(\mathbf{\omega},\mathbf{b}_{k},\hat{\mathbf{d}_{k}})$ by assuming $\hat{\mathbf{d}_{k}}$ is known. So, we have

   $$\hat{\mathbf{b}}_{k,i}=\hat{\mathbf{b}}_{k,i-1}-\mu_{b}\nabla_{\mathbf{b}_{k}}\bar{J}=\hat{\mathbf{b}}_{k,i-1}+\mu_{b}{\bf{c}}_{k}\tilde{\mathbf{e}}_{k,i}\odot\mathbf{d}^{2}_{k,i},$$

   (35)

   where ${\bf{c}}_{k}=[c_{lk}]_{l\in\cal N_{\mathrm{k}}}$, $\tilde{\mathbf{e}}_{k,i}=[\tilde{e}_{l,i}]_{l\in\cal N_{\mathrm{k}}}$, and $\mathbf{d}^{2}=\mathbf{d}\odot\mathbf{d}$.

2. 2.

   True measurement estimation: In this step, since we estimate $d_{l,i}$ from $\tilde{d}_{l,i}$, i.e., the second-order nonlinear equation, we can equivalently call this step compensation of nonlinear measurements. We use the second-order equation $b_{l}d^{2}_{l,i}+d_{l,i}-\tilde{d}_{l,i}=0$. So, neglecting the other incorrect solution, the correct solution is

   $$\hat{d}_{l,i}=\frac{-1+\sqrt{1+4\hat{b}_{l}\tilde{d}_{l,i}}}{2\hat{b}_{l}}.$$

   (36)

3. 3.

   Adaptation step: This step is equivalent to a classical adaptation step in DLMS algorithm, i.e.,

   $$\mathbf{\varphi}_{k,i}=\mathbf{\omega}_{k,i-1}-\mu\nabla_{\mathbf{\omega}}\bar{J}(\mathbf{\omega}_{k,i-1})=$$

   $$\mathbf{\omega}_{k,i-1}+\mu\sum_{l\in\cal N_{\mathrm{k}}}c_{lk}\tilde{e}_{l,i}\mathbf{u}_{l,i},$$

   (37)

   where $\tilde{e}_{l,i}=\tilde{d}_{l,i}-\hat{b}_{l}\hat{d}^{2}_{l,i}-\mathbf{u}^{T}_{l,i}\mathbf{\omega}_{k,i-1}$.

4. 4.

   Compensation of nonlinearity in the intermediate estimation: After exchanging the intermediate estimations of $\mathbf{\varphi}_{l,i}$, the received intermediate estimations, i.e., $\tilde{\mathbf{\varphi}}_{l,i}=f(\mathbf{\varphi}_{l,i})$, should be compensated by the relation of second order nonlinear equation. Similar to step 2, we write

   $$\hat{\varphi}_{l,i,j}=\frac{-1+\sqrt{1+4\hat{b}_{l}\tilde{\varphi}_{l,i,j}}}{2\hat{b}_{l}}.$$

   (38)

5. 5.

   Combination step: In this step, we update the final estimation of node $k$ as $\mathbf{\omega}_{k,i}=\sum_{l\in\cal N_{\mathrm{k}}}a_{lk}\hat{\mathbf{\varphi}}_{l,i}$.

Hence, the SONEC-DLMS algorithm is a five-step diffusion algorithm having three extra steps compared to DLMS algorithm. These three steps are, one step for estimation of nonlinear coefficients and, two other steps for compensation of nonlinearity of measurements and intermediate estimations. We call this algorithm fully-distributed SONEC-DLMS algorithm. The computational complexity of the proposed SONEC-DLMS in terms of number of addition, multiplication, and nonlinear operator of square root, in comparison to DLMS algorithm is depicted in Table 1. It is seen that the computational complexity of the proposed algorithm is approximately three times that of the DLMS algorithm.

To further improve the performance of the SONEC-DLMS algorithm, we propose a semi-distributed SONEC-DLMS algorithm, in which the nonlinear coefficients are estimated in a centralized manner by a fusion center by using training data. This version of the proposed algorithm, which can be considered as a combination of distributed and centralized algorithms, contains only four steps and the nonlinearity estimation step is done separately in a centralized manner. As we will see in the simulation results, the performance of this semi-distributed algorithm is close to that of DLMS without nonlinearity.

In this section, the Cramer-Rao bound for estimating the nonlinear coefficient vector $\mathbf{b}=[b_{l}]$ and unknown vector $\mathbf{\omega}=\mathbf{\omega}_{o}$ is derived. Let us define an $(L+N)\times 1$ parameter vector $\mathbf{\theta}=[\mathbf{\omega}^{T}\mathbf{b}^{T}]^{T}$. We intend to find the CRB of estimated $\mathbf{\theta}$ based on the measurements $\tilde{d}_{l,i}$. All the measurements are defined in a matrix $\mathbf{X}=[\mathbf{x}_{1},\mathbf{x}_{2},...,\mathbf{x}_{N}]$ where $\mathbf{x}_{l}=[\tilde{d}_{l,1},\tilde{d}_{l,2},...,\tilde{d}_{l,I}]^{T}$ is the total observations of node $l$. From (1), we have $\tilde{d}_{l,i}\sim\mathrm{N}(\mathbf{u}^{T}_{l,i}\mathbf{\omega}+b_{l}(\mathbf{u}^{T}_{l,i}\mathbf{\omega})^{2},\sigma^{2}_{\theta,l})$, where $\mathrm{N}(a,b)$ represents the Gaussian distribution with mean $a$ and variance $b$. The Fisher Information Matrix (FIM) of measurements $\mathbf{X}$ for estimating $\mathbf{\theta}$ is calculated as $\mathrm{FIM}_{\mathbf{\theta}}=\left(\begin{array}[]{cc}\mathbf{F}_{\omega}&\mathbf{F}^{T}_{b\omega}\\ \mathbf{F}_{b\omega}&\mathbf{F}_{b}\\ \end{array}\right)$, where $\mathbf{F}_{\omega,i,j}=-\mathrm{E}\Big{\{}\frac{\partial^{2}\ln p(\mathbf{X};\mathbf{\theta})}{\partial\omega_{i}\omega_{j}}\Big{\}}$, $\mathbf{F}_{b,i,j}=-\mathrm{E}\Big{\{}\frac{\partial^{2}\ln p(\mathbf{X};\mathbf{\theta})}{\partial b_{i}b_{j}}\Big{\}}$, and $\mathbf{F}_{b\omega,i,j}=-\mathrm{E}\Big{\{}\frac{\partial^{2}\ln p(\mathbf{X};\mathbf{\theta})}{\partial b_{i}\omega_{j}}\Big{\}}$.

To further proceed and to calculate the likelihood $p(\mathbf{X};\mathbf{\theta})=p(\mathbf{x}_{1},\mathbf{x}_{2},...,\mathbf{x}_{N};\mathbf{\theta})$, we write the observation vector $\mathbf{x}_{k}$ as

where $\mathbf{U}^{T}_{k}=[\mathbf{u}_{k,1}|\mathbf{u}_{k,2}|...|\mathbf{u}_{k,I}]$ and $\mathbf{z}_{k}=[\theta_{k,1},\theta_{k,2},...,\theta_{k,I}]^{T}$. Assuming the independence of measurement noises in different nodes $\mathbf{z}_{k}$, we have $\mathbf{x}_{k}\sim\mathrm{N}(\mathbf{U}_{k}\mathbf{\omega}+b_{k}(\mathbf{U}_{k}\mathbf{\omega})^{2},\mathrm{diag}(\sigma^{2}_{\theta,k}))$. Hence, the log-likelihood can be written as

So, the partial derivative is $\frac{\partial\ln p(\mathbf{X};\mathbf{\theta})}{\partial\omega_{i}\partial\omega_{j}}=\sum_{k=1}^{N}\frac{-1}{2\sigma^{2}_{\theta,k}}\frac{\partial}{\partial\omega_{i}\omega_{j}}||\mathbf{x}_{k}-\mathbf{U}_{k}\mathbf{\omega}-b_{k}(\mathbf{U}_{k}\mathbf{\omega})^{2}||^{2}$. To further proceed, we define $\mathbf{r}_{k}=\mathbf{U}_{k}\mathbf{\omega}+b_{k}(\mathbf{U}_{k}\mathbf{\omega})^{2}$.

Let $A=\frac{\partial^{2}}{\partial\omega_{i}\omega_{j}}(\mathbf{r}^{T}_{k}\mathbf{r}_{k})=\frac{\partial^{2}}{\partial\omega_{i}\omega_{j}}\sum_{p=1}^{I}r^{2}_{k,p}$, then some calculations lead to

Using (41) and taking expectations and considering that $\mathrm{E}(\mathbf{x}_{k})=\mathbf{r}_{k}$, we have