Consensus-based sparse signal reconstruction algorithm for wireless sensor networks

Abstract

This article presents a distributed Bayesian reconstruction algorithm for wireless sensor networks to reconstruct the sparse signals based on variational sparse Bayesian learning and consensus filter. The proposed approach is able to address wireless sensor network applications for a fusion-center-free scenario. In the proposed approach, each node calculates the local information quantities using local measurement matrix and measurements. A consensus filter is then used to diffuse the local information quantities to other nodes and approximate the global information at each node. Then, the signals are reconstructed by variational approximation with resultant global information. Simulation results demonstrate that the proposed distributed approach converges to their centralized counterpart and has good recovery performance.

Introduction

The ongoing concerns about environment and global warming pose more challenges in meeting the increasing demand for deployment of wireless communication networks.^1,2 The green communication (GC) has become a new trend in the design and operation of wireless communication networks. In the future communication system, the Internet of Things (IoT) will play an important role. As a key enabling technology of IoT, wireless sensor networks (WSNs) help IoT to flourish by fusing sensing and wireless communication.^3,4 Generally, a WSN consists of a large number of sensor nodes with low processing, limited power, low storage capacity, and unreliable communication over short-range radio links.^5,6 Based on the WSNs, the potential applications of IoT in industrial automation, habitat monitoring, and smart cities are numerous and diverse.^7–10 However, the energy demand for IoT will increase dramatically in the near future considering the widespread interest and adoption of various organizations, which will lead to higher carbon footprint and other environmental issues.

The recently developed compressive sensing (CS) theory^11,12 is a new sampling paradigm that can achieve acquisition of information contained in a large-scale data using much fewer samples than that are required by Nyquist sampling theorem. By exploiting sparsity, which is inherent characteristic of many natural signals, CS enables the signal to be stored in few samples and subsequently be recovered accurately. Moreover, the CS has been extensively applied in WSNs since the signals in many applications exhibit sparsity. Sparse Bayesian learning (SBL) was introduced in Tipping ¹³ and has become a popular method for sparse signal recovery in CS.^14,15 In SBL, the sparse signal recovery problem is formulated from a Bayesian perspective, while the sparsity information is exploited by assuming a hierarchical sparse prior to the signal of interest. To recover the values of model parameters, an inference algorithm is derived based on well-known Type-II Maximum Likelihood which assumes that the hyperprior are uninformative. ¹⁶ In contrast to uninformative, a fully Bayesian treatment was introduced in Bishop and Tipping ¹⁷ with the variational rendition of sparse Bayesian learning (VSBL), where both the model parameters and hyperparameters can be estimated for the distribution. ¹⁶

On the other hand, due to the high fault tolerance and scalability, distributed processing is becoming increasingly popular in WSNs applications. Different from the centralized approach that rely on a fusion center (FC), distributed processing requires no central coordinator and only single-hop communications among neighbors that aim to achieve consensus on local estimates. However, most sparse signal recovery algorithms operate in a centralized manner. Recently, distributed processing for CS applications has received considerable attention.^18–21 Moreover, the VSBL algorithm is a centralized method and may not be used for sparse signal reconstruction in a distributed WSN directly.

Average-consensus algorithms have lately investigated as a family of low-complexity iterative distributed algorithms, where sensors in a group communicate with each other to reach a consensus. ²² In more detail, each sensor receives information from others and adjusts their own information state with the goal to reach an agreement in a scalable and fault-tolerant manner. ²³ Consensus was initially elaborated in Tsitsiklis et al. ²⁴ and has received a considerable attention in many subjects due to its wide range of applications such as load balancing in parallel calculation, ²⁵ coordination of autonomous agents, distributed control, ²⁶ and data fusion.^27–29

In this article, we develop a distributed sparse signal reconstruction algorithm using probabilistic graphical models in the Bayesian framework. First, three global information quantities are particularly designed for distributed Sparse Bayesian inference by centralized update equations. Then, several average-consensus iterations are needed to reach a consensus on global information quantities in each local variational Bayesian (VB) step. In comparison with the centralized VSBL algorithms, the proposed algorithm allows each sensor to parallelly reconstruct sparse signal with local information and moderate inter-node communication.

The rest of this article is structured as follows: the fundamental of compressive sampling is provided in section “Background,” and introduces the sparse signal recovery using SBL. The system model is described in section “Problem statement and system model.” The centralized variational Bayesian inference for the system model is developed in section “Variational approximation.” Then, the proposed distributed sparse signal reconstruction is presented in section “Distributed variational SBL algorithm.” Numerical results are provided in section “Simulations,” followed by conclusions in the final section.

Notation

Throughout this article, we use b, B , and b for scalars, matrices, and column vectors, respectively. The superscripts $(·{)}^T$ and $(·{)}^{-1}$ denote the transpose and the inverse of a matrix, respectively. $E_{p(x)}(·)$ denotes the expectation with respect to $p(x)$ . ${\mathcal{U}}_{\mathrm{int}}[a,b]$ and $\mathcal{N}(\mathbf{\mu },\mathbf{\sum })$ denote the integer uniform distribution in the interval $[a,b]$ and multivariate Gaussian distribution with mean $\mathbf{\mu }$ and variance $\mathbf{\sum }$ , respectively. ${\mathit{I}}_N$ and $\mathrm{tr}$ denote the $N\times N$ identity matrix and trace of matrix, respectively. $||·|{|}_0$ , $||·|{|}_p$ , and $||·|{|}_2$ denote ${\ell }_0\mathrm{-}\mathrm{norm}$ , ${\ell }_p\mathrm{-}\mathrm{norm}$ , and ${\ell }_2\mathrm{-}\mathrm{norm}$ , respectively.

Background

In this following, we will briefly review the principle of CS and SBL,^13,30 which is a centralized sparse signal recovery algorithm.

Let $\mathbf{s}\in {\mathbb{R}}^N$ be an original signal, consider the following noisy measurement model

$$\mathit{y}=\mathbf{\Phi }\mathit{s}+\mathit{v}$$

(1)

where $\mathit{y}\in {\mathbb{R}}^M$ is the measurement vector, $\mathbf{\Phi }\in {\mathbb{R}}^{M\times N}$ denotes the measurement matrix, and $\mathit{v}\in R^M$ denotes the measurement noise. For the CS signal reconstruction, the s is assumed to have a sparse representation $\mathit{x}\in {\mathbb{R}}^N$ on some basis $\mathbf{\Psi }\in {\mathbb{R}}^{N\times N}$ , that is, $\mathit{s}=\mathbf{\Psi }\mathit{x}$ . According to the characteristics of signal, the basis has a predefined structure, for example, wavelet basis and Fourier basis. While $||\mathit{x}|{|}_0=s<<N$ , the signal can be considered as sparse over the basis. Therefore, equation (1) can be modified as follows

$$\mathit{y}=\mathit{Hx}+\mathit{v}$$

(2)

where $\mathit{H}=\mathbf{\Phi }\mathbf{\Psi }\in {\mathbb{R}}^{M\times N}$ is also referred as equivalent measurement matrix.

In order to recover x from the noisy measurement y , following optimization problem can be used

$$\underset{\mathit{x}}{min}||\mathit{x}|{|}_0,\mathrm{s}.\mathrm{t}.||\mathit{Hx}-\mathit{y}||\le ϵ$$

(3)

where $ϵ>0$ is an estimate of the measurement noise level. Since the above optimization problem (3) is non-deterministic polynomial (NP)-hard and cannot be solved effectively, the conventional CS approaches generally resort to solve the following optimization problem

$$\underset{\mathit{x}}{min}\frac{1}{2}||\mathit{Hx}-\mathit{y}|{|}_2+\lambda ||\mathit{x}|{|}_p$$

(4)

where $0<p\le 1$ and $\lambda >0$ are regularized parameters. When $p=1$ , the problem in equation (4) becomes a convex problem, and the solution of equation (3) can be obtained with overwhelming probability. For $p<1$ case, equation (4) is non-convex but a closer approximation to sparsity than $p=1$ case and shows superior performance.^31,32

From a Bayesian perspective, the CS problem can be also formulated by SBL whose close relationship to non-convex ${\ell }_p\mathrm{-}\mathrm{norm}$ minimization problem is revealed in Wipf and colleagues.^33,34 In SBL framework, a zero-mean Gaussian prior distribution is considered

$$p(\mathit{x}|\mathbf{\Gamma })=\underset{i=1}{\overset{N}{\Pi }}\mathcal{N}(x_i|0,{\gamma }_i)$$

(5)

where $\mathbf{\Gamma }\in {\mathbb{R}}^{N\times N}$ is a diagonal matrix that is composed of N hyperparameters ${\gamma }_i(i=1,\dots ,N)$ controlling the prior variance of each component. In Tipping, ¹³ and Wipf and colleagues,^30,34 the rationale of using this prior has been elaborated. With uniform hyperpriors $p({\gamma }_i)$ and $p({\sigma }^2)$ , one can infer these hyperparameters by maximizing

$$\begin{array}{c}\multicolumn{1}{c}{\log p(\mathbf{\Gamma },{\sigma }^2|\mathit{y})\propto \log p(\mathit{y}|\mathbf{\Gamma },{\sigma }^2)}\\ \multicolumn{1}{c}{=\log \int p(\mathit{y}|\mathit{x},{\sigma }^2)p(\mathit{x}|\mathbf{\Gamma })d\mathit{x}}\end{array}$$

(6)

As per the analysis in Tipping, ¹³ solving equation (6) is equivalent to minimizing the following cost function

$$\mathcal{L}=\log |\mathbf{\sum }|+{\mathit{y}}^T{\mathbf{\sum }}^{-1}\mathit{y}$$

(7)

where $\mathbf{\sum }={\sigma }^2{\mathit{I}}_M+\mathit{H}\mathbf{\Gamma }{\mathit{H}}^T$ . In Tipping, ¹³ the expectation maximization (EM) algorithm is employed to solve equation (17). Given these hyperparameters, x can be obtained by maximizing the posterior distribution

$$\begin{array}{cc}\multicolumn{1}{c}{\mathit{x}}& =\arg \underset{\mathit{x}}{max}p(\mathit{x}|\mathit{y},\mathbf{\Gamma },{\sigma }^2)\hfill \\ \multicolumn{1}{c}{}& =\arg \underset{\mathit{x}}{max}p(\mathit{y}|\mathit{x},{\sigma }^2)p(\mathit{x}|\mathbf{\Gamma })\hfill \\ \multicolumn{1}{c}{}& =\mathbf{\Gamma }{\mathit{H}}^T{\mathbf{\sum }}^{-1}\mathit{y}\hfill \end{array}$$

(8)

In Wipf and colleagues,^30,34 authors provide the theoretical justification for applying SBL to sparse signal recovery and demonstrate superior performance compared with other algorithms.

Problem statement and system model

A network with K nodes modeled by an undirected graph $\mathcal{G}=(\mathcal{V},\mathcal{E},\mathcal{A})$ of the order K is considered in this study. In the considered system, the nodes in $\mathcal{V}=\{1,2,\dots ,K\}$ represent the sensors, the edge $(k,l)$ in the set $\mathcal{E}\subseteq \mathcal{E}\times \mathcal{E}$ models that sensor l can transmit the information to sensor k, and the adjacency matrix $\mathcal{A}=[a_{\mathrm{kl}}{]}_{K\times K}$ with non-negative adjacency element $a_{\mathrm{ij}}$ . An edge of $\mathcal{G}$ is positive, that is, $a_{\mathrm{kl}}>0\iff (k,l)\in \mathcal{E}$ . Node j is called a neighbor of node i if $(k,l)\in \mathcal{E}$ and $l\ne k$ . The neighbor set of node k is denoted by ${\mathcal{N}}_k$ . Each node is able to process the data that are stored locally and exchange messages with its neighbors. An example graph is illustrated in Figure 1.

Figure 1.

An example of network structure.

Assuming that each node observes a linear combinations of unknown sparse signal $\mathit{x}\in {\mathbb{R}}^N$ , the resultant matrix will be sparse. Then, the measurement corrupted by some noise at node k is as follows

$${\mathit{y}}_k={\mathit{H}}_k\mathit{x}+{\mathit{e}}_k$$

(9)

where ${\mathit{y}}_k\in {\mathbb{R}}^{m_k}$ is the local measurement of node k, $m_k$ is the number of simultaneous measurements made at node k, ${\mathit{H}}_k\in {\mathbb{R}}^{m_k\times N}$ is the local measurement matrix of node k, $m_k$ is the number of simultaneous measurements made at node k, and ${\mathit{e}}_k\in {\mathbb{R}}^{m_k}$ is the measurement noise at node k. The measurements at all K nodes are stacked as follows. Let M be the total number of measurements at all nodes, that is, $M=\sum _{k=1}^K{m}_k$ . The global measurement vector, $\mathit{Y}\in {\mathbb{R}}^M$ , the global observation matrix, $\mathit{H}\in {\mathbb{R}}^{M\times N}$ , and the global observation noise vector, $\mathit{e}\in {\mathbb{R}}^M$ , are defined as follows

$$\mathit{Y}=\left[\begin{array}{c}{\mathit{y}}_1\hfill \\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{{\mathit{y}}_K}\end{array}\right],\mathit{H}=\left[\begin{array}{c}{\mathit{H}}_1\hfill \\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{{\mathit{H}}_K}\end{array}\right],\mathit{e}=\left[\begin{array}{c}{\mathit{e}}_1\hfill \\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{{\mathit{e}}_K}\end{array}\right]$$

(10)

Then the global observation model is given by

$$\mathit{Y}=\mathit{Hx}+\mathit{e}$$

(11)

where $\mathit{e}~\mathcal{N}(0,{\beta }^{-1}{\mathit{I}}_M)$ . The elements of H are derived from a Gaussian distribution with zero mean and variance $1/M$ . Note that this construction of H satisfies the restricted isometry property (RIP) in the design of CS schemes. From the measurement model (11) and noise statistics, the measurements likelihood function is set by

$$\begin{array}{cc}\multicolumn{1}{c}{p(\mathit{Y}|\mathit{x},\beta )}& =\underset{k=1}{\overset{K}{\Pi }}p({\mathit{y}}_k|\mathit{x},\beta )\hfill \\ \multicolumn{1}{c}{}& =\underset{k=1}{\overset{K}{\Pi }}\frac{{\beta }^{m_k/2}}{{(2\pi )}^{m_k/2}}\exp [-\frac{\beta }{2}||{\mathit{y}}_k-{\mathit{H}}_k\mathit{x}|{|}^2]\hfill \end{array}$$

(12)

Moreover, an appropriate conjugate prior with respect to equation (12) is further attached on the parameter $\beta$ so as to complementing the likelihood. The prior for the noise precision $\beta$ is selected to be a Gamma distribution with parameters $\varrho$ and $\delta$ , that is

$$p(\beta |\varrho ,\delta )=\mathcal{G}(\beta |\varrho ,\delta )=\frac{{\delta }^{\varrho }{\beta }^{\varrho -1}\exp [-\delta \beta ]}{\mathbf{\Gamma }(\varrho )}$$

(13)

To reflect our knowledge about sparsity of x , a hierarchically heavy-tailed prior is selected for the signal x next. There are two levels of hierarchy set for the model. For the first level, a Gaussian prior is attached to x , that is

$$\begin{array}{cc}\multicolumn{1}{c}{p(\mathit{x}|\mathbf{\alpha })}& =\mathcal{N}(\mathit{x}|0,{\mathit{A}}^{-1})=\underset{i=1}{\overset{N}{\Pi }}p(x_i|{\alpha }_i)\hfill \\ \multicolumn{1}{c}{}& =\underset{i=1}{\overset{N}{\Pi }}{(2\pi )}^{-1/2}{\alpha }_i^{1/2}\exp [-\frac{1}{2}{x}_i^2{\alpha }_i]\hfill \end{array}$$

(14)

where $\mathbf{\alpha }=[{\alpha }_1,{\alpha }_2,\dots ,{\alpha }_N{]}^T$ is a vector of precision parameters of $x_i\mathrm{s}$ , $\mathit{A}=\mathrm{diag}(\mathbf{\alpha })$ and the $x_i\mathrm{s}$ have been assumed to be a priori independent. Next, as the second level, the precision parameters ${\alpha }_i\mathrm{s}$ also depend on a prior distribution which is assumed to follow a Gamma distribution, expressed as

$$p({\alpha }_i|a_i,b_i)=\mathcal{G}({\alpha }_i|a_i,b_i)=\frac{b_i^{a_i}{\alpha }_i^{a_i-1}\exp [-b_i{\alpha }_i]}{\mathbf{\Gamma }(a_i)}$$

(15)

So far, the system model is developed. For the Bayesian model described above, we illustrate the directed acyclic graph (DAG) in Figure 2, where $\mathit{a}=[a_1,\dots ,a_N{]}^T$ and $\mathit{b}=[b_1,\dots ,b_N{]}^T$ . Now, the centralized variational approximation method for the developed hierarchical Bayesian model is presented in the next section.

Figure 2.

DAG of the hierarchical Bayesian model.

Variational approximation

According to Bayes’ theorem

$$p(\mathbf{\phi }|\mathit{Y},\varphi )=\frac{p(\mathit{Y}|\mathbf{\phi })p(\mathbf{\phi }|\varphi )}{p(\mathit{Y}|\varphi )}$$

(16)

where $\mathbf{\phi }=[x_1,\dots ,x_n,\beta ,{\alpha }_1,\dots ,{\alpha }_m{]}^T$ are the unknown parameters and hidden variables of the model, and $\varphi =[a_1,\dots ,a_n,b_1,\dots ,b_n,\delta ,\varrho {]}^T$ are the hyperparameters of the imposed prior. In order to estimate the value of the hyperparameters, that is, $\varphi$ , the following log-likelihood is maximized

$$\log p(\mathit{Y}|\varphi )=F(q(\mathbf{\phi }))+\mathrm{KL}(q(\mathbf{\phi })||p(\mathbf{\phi }|\mathit{Y},\varphi ))$$

(17)

where F is the free energy given by the following expression

$$F(q(\mathbf{\phi }),\varphi )=\int q(\mathbf{\phi })\log \left(\frac{p(\mathit{Y}|\mathbf{\phi }\mathbf{)}\mathbf{p}\mathbf{(}\mathbf{\phi }\mathbf{|}\mathbf{\varphi }\mathbf{)}}{q(\mathbf{\phi })}\right)d\mathbf{\phi }$$

(18)

and

$$\mathrm{KL}(\mathit{q}(\mathbf{\phi })||p(\mathbf{\phi }|\mathit{Y},\varphi ))=\int q(\mathbf{\phi })\log \left(\frac{q(\mathbf{\phi })}{p(\mathbf{\phi }|\mathit{Y},\varphi )}\right)$$

(19)

where KL is the Kullback–Leibler divergence between the true posterior $p(\mathbf{\phi }|\mathit{Y},\varphi )$ and the variational distribution $q(\mathbf{\phi })$ . It is noted here that the KL divergence is greater than or equal to zero and is minimized when $\mathit{q}(\mathbf{\phi })=p(\mathit{x}|\mathit{Y},\varphi )$ . Hence, $F(q(\mathbf{\phi }),\varphi )$ also can be regarded as evidence lower bound. The minimization of KL divergence is equivalent to maximizing the evidence lower bound. From an optimization point of view, the model parameters of $q(\mathbf{\phi })$ is well-suited so that the lower bound can be minimized. However, due to the complexity of model, computing the posterior of interest is troublesome. Thus, we resort to a simpler variational free form $q(\mathit{x},\beta ,\mathbf{\alpha })$ to approximate the posterior in equation (16).

Based on the mean-field theory from statistical physics, $q(x,\beta ,\mathbf{\alpha })$ can be fully factorized into a family of q-distribution with respect to the parameters as follows

$$\begin{array}{cc}\multicolumn{1}{c}{q(\mathit{x},\beta ,\mathbf{\alpha })}& =q(\mathit{x})q(\beta )q(\mathbf{\alpha })\hfill \\ \multicolumn{1}{c}{}& =\underset{i=1}{\overset{N}{\Pi }}q(x_i)q(\beta )\underset{i=1}{\overset{N}{\Pi }}q({\alpha }_i)\hfill \end{array}$$

(20)

that is, all model parameters are assumed to be a posteriori independent. This fully factorized form of distribution $q(\mathit{x},\beta ,\mathbf{\alpha })$ turns out to be computationally tractable. In fact, if ${\phi }_i$ denote the ith component of the vector $\mathbf{\phi }=[x_1,\dots ,x_N,\beta ,{\alpha }_1,\dots ,{\alpha }_N{]}^T$ containing the parameters of the Bayesian hierarchical model, then ${\mathbf{\phi }}_{-i}$ refers to all parameters after removing the ith component. Maximizing the free energy in equation (17) is realized by computing the functional derivatives with respect to each of the $q(·)$ distributions while fixing the other distributions and setting $\partial F(q)/\partial q(·)=0$ . Furthermore, the computation of $\partial F(q)/\partial q(·)=0$ can be expressed as follows ³⁵

$$\log q({\phi }_i)\propto E_{q({\phi }_{-i})}[\log (p(\mathit{Y}|\mathbf{\phi })p(\mathbf{\phi }|\varphi ))]$$

(21)

where $E_{q({\mathbf{\phi }}_{\mathbf{-}\mathbf{i}})}$ denotes the expectation with respect to $\underset{j\ne i}{\Pi }q({\phi }_j)$ . This cannot be solved directly, as all the factors $q({\phi }_i)$ are interdependent on each other. Then, the coupled $q({\phi }_i)$ involves an optimization process in turn, where the factors are initialized appropriately, and each one of them is updated iteratively. Moreover, the lower bound is known to be gradually increased until convergence in this process.

The log joint distribution of Y and $\mathbf{\phi }$ is expressed as follows

$$\log (p(Y|\phi )p(\phi |\varphi ))={\displaystyle \sum _{k=1}^K\log }p(y_k|x,\beta )+{\displaystyle \sum _{i=1}^{N}p}(x_i|{\alpha }_i)+{\displaystyle \sum _{k=1}^K{\displaystyle \sum _{i=1}^{m_k}p}}({\alpha }_i|a_i,b_i)+p(\beta |\delta ,\varrho )$$

(22)

Due to the conjugacy properties of the chosen distributions, as mentioned previously, a general solution (21) can be derived analytically as follows

$$\begin{array}{cc}\multicolumn{1}{c}{q(\mathit{x})}& \propto \exp (E_{q(\mathbf{\alpha }\mathbf{)}\mathbf{q}\mathbf{(}\mathbf{\beta }\mathbf{)}}[\log (p(\mathit{Y}|\mathbf{\phi })p(\mathbf{\phi }|\varphi ))])\hfill \\ \multicolumn{1}{c}{}& \propto \exp (E_{q(\beta )}[\sum _{k=1}^K\log p({\mathit{y}}_k|\mathit{x},\beta )]+E_{q(\mathbf{\alpha }\mathbf{)}}[\sum _{i=1}^{N}p(x_i|{\alpha }_i)])\hfill \\ \multicolumn{1}{c}{}& \propto \mathcal{N}(\mathit{x}|\mathit{m},\mathbf{\sum })\hfill \end{array}$$

(23)

$$\begin{array}{cc}\multicolumn{1}{c}{q(\beta )}& \propto \exp (E_{q(\mathit{x})q(\mathbf{\alpha }\mathbf{)}}[\log (p(\mathit{Y}|\mathbf{\phi })p(\mathbf{\phi }|\varphi ))])\hfill \\ \multicolumn{1}{c}{}& \propto \exp (E_{q(\mathit{x})}[\sum _{k=1}^K\log p({\mathit{y}}_k|\mathit{x},\beta )]+E_{q(\mathbf{\alpha }\mathbf{)}}[\sum _{i=1}^{N}p(x_i|{\alpha }_i)])\hfill \\ \multicolumn{1}{c}{}& \propto \mathcal{G}(\beta |\tilde{\varrho },\tilde{\delta })\hfill \end{array}$$

(24)

$$\begin{array}{cc}\multicolumn{1}{c}{q(\alpha )}& \propto \exp (E_{q(\mathit{x})q(\mathbf{\beta }\mathbf{)}}[\log (p(\mathit{Y}|\mathbf{\phi })p(\mathbf{\phi }|\varphi ))])\hfill \\ \multicolumn{1}{c}{}& \propto \exp \left(E_{q(\mathit{x})q(\beta )}[\sum _{k=1}^K\log p({\mathit{y}}_k|\mathit{x},\beta )]\right)\hfill \\ \multicolumn{1}{c}{}& \propto \mathcal{G}(\mathbf{\alpha }|\tilde{\mathit{a}},\tilde{\mathit{b}})=\underset{i=1}{\overset{N}{\Pi }}\mathcal{G}({\alpha }_i|\widetilde{a_i},\widetilde{b_i})\hfill \end{array}$$

(25)

where

$$\mathbf{\sum }={(\mathrm{diag}(E[{\alpha }_i])+E[\beta ]\sum _{k=1}^K{\mathit{H}}_k^T{\mathit{H}}_k)}^{-1}$$

(26)

$$\mathit{m}=E[\beta ]\mathbf{\sum }\sum _{k=1}^K{\mathit{H}}_k^T{\mathit{y}}_k$$

(27)

$${\tilde{a}}_i=a+\frac{1}{2}$$

(28)

$${\tilde{b}}_i=b+\frac{E[x_i^2]}{2}$$

(29)

$$\tilde{\varrho }=\varrho +\frac{M}{2}$$

(30)

$$\begin{array}{c}\multicolumn{1}{c}{\tilde{\delta }=\delta +\frac{1}{2}\sum _{k=1}^K{\mathit{y}}_k^T{\mathit{y}}_k-E[\mathit{x}{]}^T\sum _{k=1}^K{\mathit{H}}_k^T{\mathit{y}}_k}\\ \multicolumn{1}{c}{+\frac{1}{2}\mathrm{tr}(\sum _{k=1}^K{\mathit{H}}_k^T{\mathit{H}}_k·(\mathbf{\sum }+{\mathit{m}}^T\mathit{m}))}\end{array}$$

(31)

The required moments can be easily evaluated using the following results

$$\begin{array}{c}E[\mathit{x}]=\mathit{m}\hfill \\ E[x_i^2]={\mathbf{\sum }}_{\mathrm{ii}}+m_i^2\hfill \\ E[{\alpha }_i]=\frac{{\tilde{a}}_i}{\widetilde{b_i}}\hfill \\ E[\beta ]=\frac{\tilde{\varrho }}{\tilde{\delta }}\hfill \end{array}$$

(32)

From the above, it is noted that the aforementioned formulas can be used to compute the parameters of model in a centralized manner when all the measurements can be gathered in a FC. However, in the distributed scenario, there is no FC in the network and the formulas derived above cannot be implemented directly. In order to develop the distributed algorithm for VSBL, the formulas so far discussed will be reformulated in the following section such that the VSBL can be used for distributed sparse learning in a WSN.

Distributed variational SBL algorithm

Since the measurements are divided into K different nodes, the following global quantities can be defined by inspecting equations (26)–(31)

$$\begin{array}{c}\mathit{A}=\sum _{k=1}^K{\mathit{A}}_k\hfill \\ \mathit{B}=\sum _{k=1}^K{\mathit{B}}_k\hfill \\ C=\sum _{k=1}^K{C}_k\hfill \end{array}$$

(33)

where ${\mathit{A}}_k$ , ${\mathit{B}}_k$ , and $C_k$ are defined as local quantities and are defined as follows

$$\begin{array}{c}{\mathit{A}}_k={\mathit{H}}_k^T{\mathit{H}}_k\hfill \\ {\mathit{B}}_k={\mathit{H}}_k^T{\mathit{y}}_k\hfill \\ C_k={\mathit{y}}_k^T{\mathit{y}}_k\hfill \end{array}$$

(34)

Hence, equations (26)–(31) can be reformulated as follows

$$\mathbf{\sum }=(\mathrm{diag}(E[{\alpha }_i])+E[\beta ]\mathit{A}{)}^{-1}$$

(35)

$$\mathit{m}=E[\beta ]\mathbf{\sum }\mathit{B}$$

(36)

$${\tilde{a}}_i=a+\frac{1}{2}$$

(37)

$${\tilde{b}}_i=b+\frac{E[x_i^2]}{2}$$

(38)

$$\tilde{\varrho }=\varrho +\frac{M}{2}$$

(39)

$$\tilde{\delta }=\delta +\frac{1}{2}C-E[\mathit{x}{]}^T\mathit{B}+\frac{1}{2}\mathrm{tr}(\mathbf{A}·(\mathbf{\sum }+{\mathit{m}}^T\mathit{m}))$$

(40)

It should be noted that all the required parameters can be computed based on the variational approximation method and using the global quantities. In the sense of distributed calculation, each sensor communicates with its neighbors and operates accordingly. In the distributed variational sparse Bayesian learning (DVSBL) algorithm, the global quantities A , B , and C cannot be calculated locally. However, one can compute them by averaging the local quantities from all nodes in equation (34). This is because the global quantities in equation (34) can be redetermined as follows

$$\begin{array}{c}\mathit{A}=\frac{1}{K}\sum _{k=1}^K{\mathit{H}}_k^T{\mathit{H}}_k\hfill \\ \mathit{B}=\frac{1}{K}\sum _{k=1}^K{\mathit{H}}_k^T{\mathit{y}}_k\hfill \\ C=\frac{1}{K}\sum _{k=1}^K{\mathit{y}}_k^T{\mathit{y}}_k\hfill \end{array}$$

(41)

It is easy to see that above redefinition has no impact on the parameter approximations considered in equations (26)–(31). After inspecting the average formulas mentioned in equation (41), the idea of an average-consensus filter suggested in Kingston and Beard ³⁶ can be employed to approximate the global information quantities. In particular, the local information quantities of each sensor are interchanged with their neighbors, then each sensor’s global information quantities changes depending on the local information quantities input from others using consensus filter. Hence, the DVSBL algorithm can be developed by employing such average-consensus filter.

According to Kingston and Beard, ³⁶ a consensus filter can be formulated as follows in a continuous compact form

$$\stackrel{·}{\xi }=(\mathit{I}-ϵ\mathit{L})\xi$$

(42)

where L is the Laplacian matrix of graph and is defined as follows

$${\ell }_{\mathrm{ij}}=\{\begin{array}{cc}\sum _{k=1,k\ne i}^K{a}_{i,k},\hfill & j=1\hfill \\ \multicolumn{1}{c}{-a_{i,j},}& j\ne i\hfill \end{array}$$

(43)

The discrete-time form of consensus filter suggested in Kingston and Beard ³⁶ is as follows

$${\xi }_k(\tau +1)={\xi }_k(\tau )+\sum _{k=1}^K{\beta }_{\mathrm{kl}}(\tau )({\xi }_l(\tau )-{\xi }_k(\tau ))$$

(44)

where $\tau$ denotes the iteration step of discrete-time consensus filter and ${\beta }_{\mathrm{kl}}(\tau )$ is the linear weight on ${\xi }_l$ at node k. Thus, the filtering algorithm can be carried out in a distributed manner if the averages given by equation (41) are obtained at every node. Here, ${\mathit{A}}_k$ , ${\mathit{B}}_k$ , and $C_k$ are treated as input states of consensus filter mentioned in equation (44), the outputs ${\hat{\mathit{A}}}_k$ , ${\hat{\mathit{B}}}_k$ , and ${\hat{C}}_k$ can asymptotically track the values of A , B , and C, respectively. Then, the hyperparameters at each node are approximated using its estimated global quantities. Thus

$${\mathbf{\sum }}_k^t={(\mathrm{diag}(E^t[{\alpha }_i])+E^t[\beta ]K{\hat{\mathit{A}}}_k)}^{-1}$$

(45)

$${\mathit{m}}_k^t=E^t[\beta ]{\mathbf{\sum }}_k^{t}K{\hat{\mathit{B}}}_k$$

(46)

$${\tilde{a}}_i^t=a+\frac{1}{2}$$

(47)

$${\tilde{b}}_i^t=b+\frac{E[x_i^2]}{2}$$

(48)

$${\tilde{\varrho }}^t=\varrho +\frac{M}{2}$$

(49)

$${\tilde{\delta }}^t=\delta +\frac{1}{2}K{\widehat{C}}_k-E{[x]}^{T}K{\widehat{B}}_k+\frac{1}{2}\mathrm{tr}\left(K{\widehat{A}}_k·\left({\Sigma }_k^t+m_k^{tT}{m}_k^t\right)\right)$$

(50)

It can be noted from the above equations that each node should communicate with its neighbors several times before implementing variational approximation. However, the exchange of messages among the nodes iteratively is inevitable to consume time and energy and it is hard to know the number of iterations that are necessary to achieve consensus in different application scenarios, which are the main limitations of the proposed distributed VSBL. Moreover, another concerned problem is the design of weight matrix $\mathbf{\Xi }=[{\beta }_{\mathrm{kl}}(\tau ){]}_{K\times K}$ that leads to a fast convergence rate of average consensus. Usually, the weight matrix $\mathbf{\Xi }$ must be subject to the constraints of algebraic connectivity and graph topology

$$\{\begin{array}{c}\sum _{k=1}^K{\beta }_{\mathrm{kl}}(\tau )=1\hfill \\ \sum _{l=1}^K{\beta }_{\mathrm{kl}}(\tau )=1\hfill \\ {\beta }_{\mathrm{kl}}(\tau )>0\hfill \end{array}$$

(51)

Moreover, it has been shown that the second smallest eigenvalue of $\mathbf{\Xi }$ determines the convergence speed of the algorithm in the case of fixed network topology. ³⁷ Similar results have been obtained for time-varying network topology ³⁸ and many algorithms have been proposed to accelerate the convergence in Aysal et al. ³⁹ and Sardellitti et al. ⁴⁰ In this article, metropolis weight is used in the simulations ⁴¹

$${\beta }_{\mathrm{kl}}(\tau )=\{\begin{array}{c}\frac{1}{1+\max \{d_k(\tau ),d_l(\tau )\}},\mathrm{if}\{k,l\}\in \mathcal{E}(\tau )\hfill \\ \multicolumn{1}{c}{1-\sum _{l\in {\mathcal{N}}_k(\tau )}{\beta }_{\mathrm{kl}}(\tau ),\mathrm{if}k=l}\\ \multicolumn{1}{c}{0,\mathrm{otherwise}}\end{array}$$

(52)

The performance of the proposed algorithm is validated using simulations, and results are presented in the next section. Moreover, the intact algorithm is given in Algorithm 1 (Figure 3).

Algorithm 1. Distributed variational Bayesian sparse learning (DVSBL) algorithm
1: Initialization: Input measurement vectors $\{{\mathit{y}}_k\}$ , measurement matrices $\{{\mathit{H}}_k\}$ , and weight matrix $\mathbf{\Xi }$ ; initialize ${\mathbf{\sum }}_k$ by identity matrix, ${\mathit{m}}_k$ by $0$ , and $a,b,\varrho ,\delta$ with ${10}^{-6}$ .
2: for $t=1,\dots ,L$ do
3:  Compute ${\mathit{A}}_k={\mathit{H}}_k^T{\mathit{H}}_k$ , ${\mathit{B}}_k={\mathit{H}}_k^T{\mathit{y}}_k$ , $C_k={\mathit{y}}_k^T{\mathit{y}}_k$ and send $\{{\mathit{A}}_k,{\mathit{B}}_k,C_k\}$ to node $l\in {\mathcal{N}}_k$ .
4:  Receive $\{{\mathit{A}}_l,{\mathbf{B}}_l,C_l\}$ from node $l\in N_k$ and fuse the local information quantities using the consensus filter (equation (44));
5:  Compute equations (45)–(50) at each node;
6:  If stop criterion is met, return ${\mathit{x}}_k(k=1,\dots ,K)$ ; otherwise go to step 5;
7: end for

Figure 3.

Algorithm architecture.

Simulations

First, a sensor network with six nodes is used to verify the performance of the proposed algorithm for distributed WSNs (Figure 4). Without loss of generality, the considered six-node network is represented by an undirected graph $\mathcal{G}=(\mathcal{V},\mathcal{E},\mathcal{A})$ with a set of nodes $\mathcal{A}=(1,2,3,4,5,6)$ , a set of edges $\mathcal{E}=\{(1,1),(1,2),(1,3),(2,3),(2,4),(2,5),(3,3),(3,5),(3,6),(4,4),(4,5),(5,5),(5,6),(6,6)\}$ , and an adjacency matrix

$$\mathcal{A}=\left[\begin{array}{cccccc}1\hfill & 1\hfill & 1\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{1}& 1\hfill & 1\hfill & 1\hfill & 1\hfill & 0\hfill \\ \multicolumn{1}{c}{1}& 1\hfill & 1\hfill & 0\hfill & 1\hfill & 1\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 0\hfill & 1\hfill & 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 1\hfill & 1\hfill & 1\hfill & 1\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill & 0\hfill & 1\hfill & 1\hfill \end{array}\right]$$

(53)

Figure 4.

Random topology.

In the considered example, the signal $\mathit{x}\in {\mathbb{R}}^{256}$ is assumed to be sparse itself. There are altogether 10 non-zero elements $\mathit{x}\{i\}\ne 0$ in the sparse signal, where i is the index of support and $\mathit{x}\{i\}$ denotes the value of support. Here, the index and value of support are unknown and sampled over $i~{\mathcal{U}}_{\mathrm{int}}[1,256]$ and $\mathit{x}\{i\}~\mathcal{N}(0,5^2)$ , respectively. The measurement matrix ${\mathit{H}}_k\in {\mathbb{R}}^{12\times 256}$ of sensor node k is constructed by the entries sampled from $\mathcal{N}(0,(1/72))$ . Then, the local measurement ${\mathit{y}}_k\in {\mathbb{R}}^{12}$ is obtained using equation (1) and the global measurement matrix $\mathit{H}=\mathrm{col}\{{\mathit{H}}_1,\dots ,{\mathit{H}}_6\}\in {\mathbb{R}}^{72\times 256}$ satisfies the RIP with overwhelming probability. The measurement noise ${\mathit{e}}_k~\mathcal{N}(0,{10}^{-4}\times {\mathit{I}}_{12})$ . The normalized mean square error (MSE) and averaged relative error are employed to evaluate the performance, which are defined as follows

$$\mathrm{Normalized\; MSE}=\frac{||\mathit{x}-\hat{\mathit{x}}|{|}_2}{||\mathit{x}|{|}_2}$$

(54)

$$\mathrm{Average\; relative\; error}=\frac{\sum _{k=1}^K||{\hat{\mathit{x}}}_k-\mathit{x}|{|}_2^2}{K·||\mathit{x}|{|}_2^2}$$

(55)

The convergence rate of the proposed algorithm is shown in Figures 5 and 6. It can be seen from the figures that both the local information quantity $C_k$ and the normalized MSE reach the consensus as the iteration progresses. It can also be observed that the normalized MSE has the same convergence tendency with the local information quantity. The results confirm that the consensus iterations enhance the observability by information, which is in agreement with the analysis.

Figure 5.

Convergence for local information quantity $C_k$ .

Figure 6.

Normalized MSE versus iteration.

The estimated x from sensor nodes using DVSBL algorithm is presented in Figure 7. It can be observed that the estimated values of actual sparse signals at all nodes are satisfactory.

Figure 7.

Estimation of x .

In order to prove the scalability of the proposed distributed algorithm, we consider two networks of different sizes, one of which is an L-connected Harary graph formed by 24 nodes, and the other is an L-connected Harary graph formed by 72 nodes. The L-connection denotes the number of neighbors for each node. In the considered example, the L is set as 3. The error performance is shown in Figures 8 and 9. It can be noted that irrespective of the size, both the networks have the comparable reconstruction performance compared to the centralized VSBL.

Figure 8.

Comparison of Normalized MSE (24 nodes).

Figure 9.

Comparison of Normalized MSE (72 nodes).

The average relative error of 72 nodes is presented in Figure 10. It can be seen that the relative error after two iterations is negligible. This demonstrates that the proposed algorithm is salable.

Figure 10.

ARE versus iteration.

Finally, the real temperature signals obtained from Intel Berkeley Research laboratory are considered to test the effectiveness of the proposed algorithm. The considered temperature signals are represented using discrete cosine transform (DCT) as shown in Figure 11. In this context, a three-connected Harary graph with 24 nodes is employed. The effectiveness of our distributed algorithm is demonstrated in Figure 12 where both the original signal and the reconstructed signal of one node are provided. It can be observed from Figure 12 that both the proposed distributed algorithm and the centralized VSBL are successful in reconstructing the temperature signals.

Figure 11.

The coefficients of DCT.

Figure 12.

Reconstructed results of temperature signal.

In a word, all the results presented so far have demonstrated the effectiveness of the distributed sparse signal reconstruction algorithm presented in this article.

Conclusion

In this article, a new distributed variational SBL algorithm for sparse signal reconstruction in a WSNs is presented. By combining VSBL and average-consensus algorithm, the DVSBL is consistent with centralized VSBL, where all data are available at FC. Experimental results demonstrate the superior recovery performance and convergence properties of the proposed distributed algorithm.