On Distributed Estimation in Hierarchical Power Constrained Wireless Sensor Networks

TL;DR

This paper investigates optimal power allocation and signal fusion strategies in hierarchical wireless sensor networks to minimize estimation error under power constraints, providing theoretical bounds and practical schemes.

Contribution

It introduces a comprehensive framework for optimizing training, transmission, and fusion weights in hierarchical sensor networks, including analysis of special cases and complexity tradeoffs.

Findings

Derived the MSE and Bayesian CRB for the estimation problem.

Proposed optimal power and fusion weight schemes to minimize MSE.

Quantified performance gains and complexity tradeoffs of different schemes.

Abstract

We consider distributed estimation of a random source in a hierarchical power constrained wireless sensor network. Sensors within each cluster send their measurements to a cluster head (CH). CHs optimally fuse the received signals and transmit to the fusion center (FC) over orthogonal fading channels. To enable channel estimation at the FC, CHs send pilots, prior to data transmission. We derive the mean square error (MSE) corresponding to the linear minimum mean square error (LMMSE) estimator of the source at the FC, and derive the Bayesian Cramer-Rao bound (CRB). Our goal is to find (i) the optimal training power, (ii) the optimal power that sensors in a cluster spend to transmit their amplified measurements to their CH, and (iii) the optimal weight vector employed by each CH for its linear signal fusion, such that the MSE is minimized, subject to a network power constraint. To untangle the performance gain that optimizing each set of these variables provide, we also analyze three special cases of the original problem, where in each special case, only two sets of variables are optimized across clusters. We define three factors that allow us to quantify the effectiveness of each power allocation scheme in achieving an MSE-power tradeoff that is close to that of the Bayesian CRB. Combining the information gained from the factors and Bayesian CRB with our computational complexity analysis provides the system designer with quantitative complexity-versus-MSE improvement tradeoffs offered by different power allocation schemes.

Tables

Table 1. TABLE I : Notations and their corresponding definitions.

Notation	Vector and Matrix Definitions
${𝒙}_l,{𝒕}_l$	${𝒙}_l={[x_{l,1},\dots ,x_{l,K_l}]}^T$, ${𝒕}_l={[t_{l,1},\dots ,t_{l,K_l}]}^T$
${𝒏}_l,{𝒒}_l$	${𝒏}_l={[n_{l,1},\dots ,n_{l,K_l}]}^T$, ${𝒒}_l={[q_{l,1},\dots ,q_{l,K_l}]}^T$
$\sqrt{{𝑨}_l}$	$\sqrt{{𝑨}_l}=\text{diag}(\sqrt{{\alpha }_{l,1}},\dots ,\sqrt{{\alpha }_{l,K_l}})$
$𝒙,𝒕$	$𝒙={[𝒙_1^{}{}_{}{}^T,\dots ,𝒙_L^{}{}_{}{}^T]}^T$, $𝒕={[𝒕_1^{}{}_{}{}^T,\dots ,𝒕_L^{}{}_{}{}^T]}^T$
$𝒚,𝒛$	$𝒚={[y_1,\dots ,y_L]}^T$, $𝒛={[z_1,\dots ,z_L]}^T$
$𝒏,𝒒$	$𝒏={[𝒏_1^{}{}_{}{}^T,\dots ,𝒏_L^{}{}_{}{}^T]}^T$, $𝒒={[𝒒_1^{}{}_{}{}^T,\dots ,𝒒_L^{}{}_{}{}^T]}^T$
$𝒗,𝑯$	$𝒗={[v_1,\dots ,v_L]}^T$, $𝑯=\text{diag}\left([h_1,\dots ,h_L]\right)$
$𝑴,𝑾$	$𝑴=\text{diag}(\sqrt{{𝑨}_1},\dots ,\sqrt{{𝑨}_L})$, $𝑾=\text{diag}(𝒘_1^{}{}_{}{}^T,\dots ,𝒘_L^{}{}_{}{}^T)$
${𝚺}_n,{𝚺}_{n_l}$	${𝚺}_n=\text{diag}({𝚺}_{n_1},\dots ,{𝚺}_{n_L})$, ${𝚺}_{n_l}$ is arbitrary
${𝚺}_q,{𝚺}_{q_l}$	${𝚺}_q=\text{diag}({𝚺}_{q_1},\dots ,{𝚺}_{q_L})$, ${𝚺}_{q_l}=\text{diag}({\sigma }_{q_{l,1}}^2,\dots ,{\sigma }_{q_{l,K_l}}^2)$
${𝚺}_v$	${𝚺}_v=\text{diag}\left([2{\sigma }_{v_1}^2,\dots ,2{\sigma }_{v_L}^2]\right)$
$\widehat{𝑯},\widetilde{𝑯}$	$\widehat{𝑯}=\text{diag}\left([{\hat{h}}_1,\dots ,{\hat{h}}_L]\right),\widetilde{𝑯}=\text{diag}\left([{\tilde{h}}_1,\dots ,{\tilde{h}}_L]\right)$
$𝚪,𝚺$	$𝚪=\text{diag}\left([{\zeta }_1,\dots ,{\zeta }_L]\right)$, $𝚺=\text{diag}({𝚺}_1,\dots ,{𝚺}_L)$
$𝝁,{𝚲}_1$	$𝝁={[𝝁_1^{}{}_{}{}^H,\dots ,𝝁_L^{}{}_{}{}^H]}^H$, ${𝚲}_1=\text{diag}({𝚲}_{1_1},\dots ,{𝚲}_{1_L})$
${𝑫}_l$	${𝑫}_l=\text{diag}\left([\sqrt{d_{l,1}},\dots ,\sqrt{d_{l,K_l}}]\right)$

Equations

$$x_{l,k}=\theta +n_{l,k},l=1,...,L,k=1,...,K_l,\text{vspace}-.1cm$$

$$t_{l,k}=\sqrt{{\alpha }_{l,k}}{x}_{l,k}+q_{l,k},l=1,...,L,k=1,...,K_l,\text{vspace}-.1cm$$

$${\mathbf{x}}_l=\theta 1_l+{\mathbf{n}}_l,{\mathbf{t}}_l=\sqrt{{\mathbf{A}}_l}{\mathbf{x}}_l+{\mathbf{q}}_l,l=1,...,L.\text{vspace}-.1cm$$

$$z_l=h_l{y}_l+v_l,l=1,...,L,\text{vspace}-.1cm$$

$${\hat{z}}_l=h_l\sqrt{{\psi }_l}+{\nu }_l,l=1,...,L,\text{vspace}-.1cm$$

$$P_{l,k}={\alpha }_{l,k}\mathbb{E}\{x_{l,k}^2\}={\alpha }_{l,k}({\sigma }_{\theta }^2+{\sigma }_{n_{l,k}}^2),k=1,...,K_l.\text{vspace}-.1cm$$

$${\mathcal{P}}_l=\mathbb{E}\{y_l^2\}={{\mathbf{w}}_l}^T\underset{={\mathbf{R}}_{t_l}}{\underbrace{\mathbb{E}\{{\mathbf{t}}_l{{\mathbf{t}}_l}^T\}}}{\mathbf{w}}_l.\text{vspace}-.1cm$$

$${\mathbf{R}}_{t_l}=P_l{\mathbf{\Omega }}_l+{\mathbf{\Sigma }}_{q_l},\text{vspace}-.1cm$$

$${\mathbf{\Omega }}_l={\mathbf{\Delta }}_l+{\sigma }_{\theta }^2{\mathbf{\Pi }}_l,{\mathbf{\Delta }}_l={\mathbf{D}}_l{\mathbf{\Sigma }}_{n_l}{\mathbf{D}}_l,{\mathbf{\Pi }}_l={\mathbf{\rho }}_l{{\mathbf{\rho }}_l}^T,{\mathbf{\rho }}_l={\mathbf{D}}_l{1}_l.\text{vspace}-.1cm$$

$${\mathcal{P}}_l={{\mathbf{w}}_l}^T(P_l{\mathbf{\Omega }}_l+{\mathbf{\Sigma }}_{q_l}){\mathbf{w}}_l.\text{vspace}-.1cm$$

$$P_{trn}+\sum _{l=1}^L{P}_l+{\mathcal{P}}_l\le P_{tot}.\text{vspace}-.1cm$$

$$P_{trn}+\sum _{l=1}^L{{\mathbf{w}}_l}^T{\mathbf{\Sigma }}_{q_l}{\mathbf{w}}_l+P_l(1+{{\mathbf{w}}_l}^T{\mathbf{\Omega }}_l{\mathbf{w}}_l)\le P_{tot}.\text{vspace}-.1cm$$

$${\min }_{P_{trn},\{P_l,{\mathbf{w}}_l{\}}_{l=1}^L}$$

$$P_{trn}+\sum _{l=1}^L{{\mathbf{w}}_l}^T{\mathbf{\Sigma }}_{q_l}{\mathbf{w}}_l+P_l(1+{{\mathbf{w}}_l}^T{\mathbf{\Omega }}_l{\mathbf{w}}_l)\le P_{tot},$$

$$P_{trn}\in {\mathbb{R}}^{+},P_l\in {\mathbb{R}}^{+},{\mathbf{w}}_l\in {\mathbb{R}}^{K_l},\forall l.\text{vspace}-.1cm$$

$$\mathbf{x}=\theta 1+\mathbf{n},\mathbf{t}$$

$$\mathbf{z}$$

$${\hat{h}}_l=\frac{{\sigma }_{h_l}^2\sqrt{{\psi }_l}{\hat{z}}_l}{{\sigma }_{v_l}^2+{\psi }_l{\sigma }_{h_l}^2},{\zeta }_l^2=\frac{2{\sigma }_{h_l}^2{\sigma }_{v_l}^2}{{\sigma }_{v_l}^2+{\psi }_l{\sigma }_{h_l}^2}.\text{vspace}-.1cm$$

$$\mathbf{z}=\underset{={\mathbf{z}}_1}{\underbrace{[\hat{\mathbf{H}}\mathbf{W}\mathbf{M}1]\theta }}+\underset{={\mathbf{z}}_2}{\underbrace{(\tilde{\mathbf{H}}\mathbf{W}\mathbf{M}1)\theta }}+\underset{={\mathbf{z}}_3}{\underbrace{(\hat{\mathbf{H}}+\tilde{\mathbf{H}})\mathbf{W}(\mathbf{q}+\mathbf{M}\mathbf{n})+\mathbf{v}}}.$$

$$\hat{\theta }$$

$$D$$

$${\mathbf{C}}_{\mathbf{z}}$$

$$+\hat{\mathbf{H}}\mathbf{W}({\mathbf{\Sigma }}_q+\mathbf{M}{\mathbf{\Sigma }}_n\mathbf{M}){\mathbf{W}}^T{\hat{\mathbf{H}}}^H$$

$$+\mathbf{\Gamma }\mathbf{W}({\mathbf{\Sigma }}_q+\mathbf{M}{\mathbf{\Sigma }}_n\mathbf{M}){\mathbf{W}}^T\mathbf{\Gamma }+{\mathbf{\Sigma }}_v,\text{vspace}-.1cm$$

$${\mathbf{\Lambda }}_{1_l}={\sigma }_{\theta }^2{\zeta }_l^2{P}_l{\mathbf{\Pi }}_l+({|{\hat{h}}_l|}^2+{\zeta }_l^2)({\mathbf{\Sigma }}_{q_l}+P_l{\mathbf{\Delta }}_l),$$

$${\mathbf{\Lambda }}_2=|\mathbf{\mu }|{|\mathbf{\mu }|}^T,{\mathbf{\mu }}_l=\sqrt{P_l}{\hat{h}}_l{\mathbf{\rho }}_l,\forall l.\text{vspace}-.1cm$$

$${\mathbf{C}}_{\mathbf{z}}=\mathbf{W}({\mathbf{\Lambda }}_1+{\sigma }_{\theta }^2{\mathbf{\Lambda }}_2){\mathbf{W}}^T+{\mathbf{\Sigma }}_v.\text{vspace}-.1cm$$

$$\mathbb{E}\{\theta \mathbf{z}\}=\mathbb{E}\{\theta {\mathbf{z}}_1\}+\mathbb{E}\{\theta {\mathbf{z}}_2\}\stackrel{(a)}{=}{\sigma }_{\theta }^2\hat{\mathbf{H}}\mathbf{W}\mathbf{M}1={\sigma }_{\theta }^2\mathbf{W}\mathbf{\mu },$$

$$\hat{\theta }$$

$$D$$

$$D_1=({\sigma }_{\theta }^{-2}+\sum _{l=1}^L\frac{P_l{|h_l|}^2{{\mathbf{w}}_l}^T{\mathbf{\Pi }}_l{\mathbf{w}}_l}{{\sigma }_{v_l}^2+{|h_l|}^2{{\mathbf{w}}_l}^T({\mathbf{\Sigma }}_{q_l}+P_l{\mathbf{\Delta }}_l){\mathbf{w}}_l}{)}^{-1}.\text{vspace}-.1cm$$

$$D_2=({\sigma }_{\theta }^{-2}+\sum _{l=1}^L\frac{{|h_l|}^2{{\mathbf{w}}_l}^T{\mathbf{\Sigma }}_l{\mathbf{w}}_l}{{\sigma }_{v_l}^2+{|h_l|}^2{{\mathbf{w}}_l}^T{\mathbf{\Sigma }}_{n_l}{\mathbf{w}}_l}{)}^{-1}.\text{vspace}-.1cm$$

$$D_3=({\sigma }_{\theta }^{-2}+\sum _{l=1}^L{1_l}^T{\mathbf{\Sigma }}_{n_l}^{-1}{1}_l{)}^{-1}.\text{vspace}-.1cm$$

Full text

I Introduction

The plethora of wireless sensor network (WSN) applications, with stringent power constraints, raises challenging technical problems for system-level engineers, one of which is distributed estimation (DES) in a power constrained WSN [1, 2, 3, 4, 5]. In this work, we address DES of a random signal $\theta$ in a WSN, where sensors are deployed in a large field and make noisy measurements of $\theta$. Due to limited communication range, however, the battery-powered sensors cannot directly communicate with the fusion center (FC). Hence, the field is divided into $L$ geographically disjoint zones (clusters) and hierarchically into three tiers: sensors, cluster-heads (CHs) one per cluster, and the FC [6, 7, 8, 9]. The implicit assumption is that the communication ranges of CHs are larger, and their energy and computational resources are higher (compared with sensors). After local signal processing, CHs transmit signals received from sensors over orthogonal fading channels to the FC, whose task is to find an estimate of $\theta$, based on the received signals from CHs [10, 11].

There is a rich body of literature on DES and distributed detection in a power constrained WSN, where the researchers study and optimize an estimation-theoretic-based or a detection-theoretic-based performance metric, subject to power constraints. Examples in the context of distributed detection are [12, 13, 14, 15, 16]. An alternative direction is to study the outlier contamination of the data in WSNs by outlier detection methods such as [17, 18, 19] caused by imperfect sensors and power deteriorations [20, 21]. We focus on power optimization and to conserve space, we elaborate only the most related ones to our current work in the following. The authors in [6, 7, 8] studied DES in a three-layered hierarchical power constrained WSN, assuming that the FC forms the linear minimum mean square error (LMMSE) estimate of random $\theta$, and the objective is to minimize the MSE of this estimator. The authors in [22] considered DES in a WSN, where sensors transmit to the FC over orthogonal fading channels and the FC finds the LMMSE estimate of $\theta$. The authors studied how partial channel state information (CSI) at the sensors affects the MSE performance and the optimal power allocation among the sensors. The authors in [23] considered DES in a hierarchical WSN, where the CHs amplify and forward their received signals over orthogonal Nakagami fading channels to the FC. Assuming the FC finds the LMMSE estimate of $\theta$, the authors studied how partial CSI at the CHs impacts the outage probability of the MSE. None of the works in [22, 23] consider the cost of channel estimation at the FC. To enable channel estimation at the FC, each CH needs to transmit a training (pilot) symbol, prior to data symbol. In a hierarchical WSN, where there is a cap on the network transmit power, the cost of channel estimation cannot be overlooked. Note that training symbol transmission consumes the power that could have been used otherwise for data symbol transmission. Hence, training and data transmit power should be optimized judiciously, such that the estimation accuracy of $\theta$ at the FC is maximized.

Assuming the FC employs the LMMSE estimator of $\theta$, we address this problem, by formulating and solving a new optimization problem that allows us to analyze the effect of channel estimation on the MSE performance and transmit power allocation. The optimization problem is novel, since considering training transmit power introduces a new dimension to the network performance analysis and power allocation optimization. In this regard, the most relevant works are [24, 25], where the authors considered channel estimation for DES in a WSN with one FC only. Our work is different from [24, 25], since in the hierarchical WSN, our problem formulation considers power distribution among different clusters for sensor-CH data transmission as well as power allocation among different CHs for CH-FC data and training transmissions. Moreover, we obtain the optimal linear fusion rules at CHs as the by-product of solving the network power allocation problem111We note that there is a rich body of literature on clustering algorithms and energy efficient routing protocols [26]. Similar to [6, 7, 8, 9], we assume that clusters and their CHs are given. Given this network structure, our goal is designing (sub-)optimal distributed signal processing such that the MSE distortion at the FC is minimized, under a network power constraint..

Contribution: We derive the MSE corresponding to the LMMSE estimator of $\theta$ at the FC, denoted as $D$, and establish lower bounds on $D$, including the Bayesian Cramér-Rao bound (CRB). We then formulate a new constrained optimization problem that minimizes $D$, subject to network transmit power constraint $P_{tot}$, where the optimization variables are: i) training power for CHl, ii) total power that sensors in cluster $l$ spend to transmit their amplified measurements to CHl (which we refer to as intra-cluster power), iii) power that CHl spends to send its fused signal to the FC. We demonstrate the superior performance of our proposed power allocation scheme with respect to the following spacial case schemes: scheme (i) allots a fixed percentage of $P_{tot}$ for training power and distributes this power equally among CHs, however, it optimally allocates intra-cluster power among clusters, and optimally allocates power among CHs for data transmission, scheme (ii) optimally allocates power among CHs for training, equally allocates intra-cluster power among clusters, and optimally allocates power among CHs for data transmission, scheme (iii) optimally allocates power among CHs for training, optimally allocates intra-cluster power among clusters, and equally allocates power among CHs for data transmission. We analytically and numerically compare the power allocation scheme obtained from solving the original problem with the special case schemes, and show their effectiveness in providing an MSE-power tradeoff that is close to that of the Bayesian CRB. Our numerical results demonstrate that power allocations among CHs for training and CH-FC data transmission are always beneficial for low-region of $P_{tot}$, and power allocation among clusters for sensor-CH data transmission is beneficial for low-region to moderate-region of $P_{tot}$.

**Organization: The rest of the paper is organized as follows. Section II describes our system model and power constraints and states the problem we aim to solve (i.e., the constrained minimization of MSE $D$ at the FC, with respect to three sets of optimization variables). Section III characterizes $D$ and its lower bounds. We also derive the Bayesian CRB. In Section IV we solve our proposed constrained MSE minimization problem. We also briefly discuss the constrained minimization of MSE lower bounds. In Section V we solve three special cases of the original problem, where in each special case, only two (of three) sets of variables are optimized across clusters. This analysis allows us to entangle the performance gain that optimizing each set of these variables provide. Section VI compares the computational complexity of the proposed algorithms for solving the original problem as well as its three special cases. In Section VII we discuss the convergence analysis of our proposed algorithms. Section VIII presents our numerical and simulation results. Section IX concludes the work and outlines our future research directions. **

Notations: Matrices are denoted by bold uppercase letters, vectors by bold lowercase letters, and scalars by normal letters. $\mathbb{E}$ denotes the mathematical expectation operator, $[.]^{T}$ represents the matrix-vector transpose operation, and $|\cal A|$ is the cardinality of set $\cal A$. The real and imaginary parts of a complex random variable $x$ are represented by $x_{r}\!=\!{\cal{R}}e\{x\}$ and $x_{i}\!=\!{\cal{I}}m\{x\}$. The probability distribution function (pdf) of $x$, denoted as $f(x)$, is defined as the joint pdf of $x_{r}$ and $x_{i}$, i.e., we have $f(x)\!=\!f(x_{r},x_{i})$ [27].

II System Model and Problem Formulation

II-A System Model Description

We consider a DES problem in a hierarchical power constrained WSN (see Fig. 1), consisting of $K$ spatially-distributed sensors deployed in $L$ disjoint clusters, $L$ cluster heads (CHs), and a FC. Each sensor makes a noisy measurement of an unknown random variable $\theta$, that we wish to estimate at the FC. Cluster $l$ includes $K_{l}$ sensors and its associated CH, denoted as CHl, and we have $\sum_{l=1}^{L}K_{l}=K$. We assume $\theta$ is zero-mean with variance $\sigma^{2}_{\theta}$. Let $x_{l,k}$ denote the measurement of sensor $k$ in cluster $l$. We have:

[TABLE]

where $n_{l,k}\!$ denotes zero-mean additive measurement noise with variance $\sigma_{n_{l,k}}^{2}\!$. We assume that $n_{l,k}$’s are correlated across sensors, due to their proximity within cluster $l$. Sensors within a cluster amplify and forward their measurements to their respective CH over orthogonal AWGN channels222The AWGN channel model is equivalent to the channel model with a static and known channel gain. Given the channel gain, CHl can equalize it, which is equivalent to scaling the communication noise variance $\sigma_{q_{l,k}}^{2}$. The AWGN channel model for communication channels within a cluster is reasonable, since sensors are closely located and typically there are direct line of sight transmissions between sensors and their CH [28, 9]. On the other hand, we model the communication channels between CHs and the FC as randomly-varying fading channels that require channel estimation. The reason is that the transmission distances between CHs and the FC are large and hence communication becomes subject to multipath fading effect., such that the received signal at CHl from sensor $k$ within cluster $l$ is:

[TABLE]

where ${\alpha}_{l,k}\geq 0$ is an amplifying factor (to be determined) used by sensor $k$, and $q_{l,k}\sim{\cal N}(0,\sigma_{q_{l,k}}^{2})$ is the additive communication channel noise. We assume that $q_{l,k}$’s are uncorrelated across the sensors. For a compact representation, we define the column vectors $\boldsymbol{x}_{l}$ and $\boldsymbol{t}_{l}$ in Table I corresponding to cluster $l$ and rewrite (1) and (2) as:

[TABLE]

where $\boldsymbol{1}_{l}$ is a column vector of $K_{l}$ ones, column vectors $\boldsymbol{n}_{l}$, $\boldsymbol{q}_{l}$, with covarinace matrices $\boldsymbol{\Sigma}_{n_{l}},\boldsymbol{\Sigma}_{q_{l}}$, matrix $\sqrt{\boldsymbol{A}_{l}}$ defined in Table I.

We assume $\boldsymbol{n}_{l}$, $\boldsymbol{q}_{l}$, $\theta$ are uncorrelated, i.e., \mathbb{E}\{\boldsymbol{n}_{l}{\boldsymbol{q}_{l}}^{T}\}\!=\!\boldsymbol{0},~{}{\color[rgb]{0.2,0.3,0.8}\mathbb{E}\{\boldsymbol{n}_{l}\theta\}\!=\!\boldsymbol{0},\ \mathbb{E}\{\boldsymbol{q}_{l}\theta\}\!=\!\boldsymbol{0}},~{}\forall l, and the noise vectors across different clusters are mutually uncorrelated, i.e., $\mathbb{E}\{\boldsymbol{n}_{i}{\boldsymbol{n}_{j}}^{T}\}\!=\!\boldsymbol{0}$ and $\mathbb{E}\{\boldsymbol{q}_{i}{\boldsymbol{q}_{j}}^{T}\}\!=\!\boldsymbol{0},\ \forall i\!\neq\!j$.

Each CH linearly fuses the signals received from the sensors within its cluster. Let $y_{l}\!=\!{\boldsymbol{w}_{l}}^{T}\boldsymbol{t}_{l}$, where $y_{l}$ is the scalar fused signal at CHl and $\boldsymbol{w}_{l}$ is the linear weight vector employed by CHl for linear fusion333When the pdf of $\theta$ is unknown, it is reasonable to assume that CHl applies a linear fusion rule $\boldsymbol{w}_{l}$ and we seek the best $\boldsymbol{w}_{l}$. In Section IV-A1 we show that $\boldsymbol{w}_{l}^{opt}$ is equal to the linear operator corresponding to the LMMSE estimation of $\theta$ based on $\boldsymbol{t}_{l}$, multiplied by an optimized scalar $\chi_{l}$. When $\theta\sim N(0,\sigma_{\theta}^{2})$ the MMSE and LMMSE estimates of $\theta$ based on $\boldsymbol{t}_{l}$ coincide. (to be optimized). CHs transmit these fused signals to the FC over orthogonal Rayleigh fading channels, such that the received signal at the FC from CHl is:

[TABLE]

where $h_{l}\sim\mathcal{CN}\left(0,2\sigma_{h_{l}}^{2}\right)$ is fading channel coefficient corresponding to the link between CHl and the FC and $v_{l}\sim\mathcal{CN}\left(0,2\sigma_{v_{l}}^{2}\right)$ is the additive communication channel noise. We assume $v_{l}$ is uncorrelated with $\theta,\boldsymbol{n}_{l},\boldsymbol{q}_{l}$, $\forall l$.

To enable estimating $h_{l}$ at the FC, CHl transmits a pilot symbol [14] with power $\psi_{l}$ to the FC, prior to sending its signal $y_{l}$. Without loss of generality, we assume training symbols are all ones. Assuming $h_{l}$ does not change during transmission of $y_{l}$ and the training symbol444 We assume time-division-duplex transmission and channel reciprocity. We also assume that the channel coherence time is larger than the overall duration of pilot transmission, channel estimation, power optimization, information feedback, and data transmission., the received signal at the FC from CHl corresponding to the training symbol is:

[TABLE]

where $\nu_{l}$ in (5) is independent of $v_{l}$ in (4) and is identically distributed. The FC adopts the following two-stage strategy to process the received signals from the CHs and reconstruct $\theta$: stage 1) the FC uses the received signals $\{\hat{z}_{l}\}_{l=1}^{L}$ corresponding to training symbol transmissions to estimate $\{h_{l}\}_{l=1}^{L}$ and obtain the channel estimates $\{\hat{h}_{l}\}_{l=1}^{L}$, stage 2) the FC uses these channel estimates and the received signals $\{z_{l}\}_{l=1}^{L}$ corresponding to $\{y_{l}\}_{l=1}^{L}$ transmissions and find the LMMSE estimate of $\theta$, denoted as $\hat{\theta}$. Finding the LMMSE estimator has a lower computational complexity, compared with the optimal MMSE estimator, and it requires only the knowledge of first and second order statistics. Let $D\!=\!\mathbb{E}\{(\theta-\hat{\theta})^{2}\}$ denote the MSE corresponding to the LMMSE estimator $\hat{\theta}$. Our main objective is to study power allocation among different clusters, subject to a network transmit power constraint (including power for training and data transmissions), such that $D$ is minimized. Section II-B provides a formal description of our constrained optimization problem, including the power constraints and the set of our optimization variables555Comparing orthogonal channel model and multiple-access channel (MAC) model adopted in [25, 7, 13], the former consumes more time or bandwidth for transmission, however, it does not require symbol-level synchronization for compensating complex channel phase at transmitter. We note that the complexity of the sequence of operation in our work (pilot transmission, channel estimation, power optimization, information feedback, and data transmission) is comparable with that of those works that rely upon perfect CSI, since implementing power allocation solutions obtained based on perfect CSI [11, 6, 8, 7, 13] requires pilot transmission and channel estimation, prior to data transmission..

II-B Power Constraints

We describe our power constraints. Let $P_{l,k}$ denote the average power that sensor $k$ consumes to send its amplified measurement to CHl and $P_{l}\!=\!\sum_{k=1}^{K_{l}}P_{l,k}$ be the total power that sensors in cluster $l$ spend to send their amplified measurements to CHl. From (2) we have:

[TABLE]

For tractability, similar to [29, 9] we assume $P_{l}$ is equally divided between sensors within cluster $l$, i.e., $P_{l,k}\!=\!P_{l}/K_{l}$. Under this assumption from (6) we obtain ${\alpha}_{l,k}=P_{l}d_{l,k}$ where $d_{l,k}\!=\!\frac{1}{K_{l}(\sigma^{2}_{\theta}+\sigma_{n_{l,k}}^{2})}$, or equivalently in matrix form, we find $\sqrt{\boldsymbol{A}_{l}}\!=\!\sqrt{P_{l}}\boldsymbol{D}_{l}$, where $\boldsymbol{D}_{l}$ is given in Table I. Let ${\mathcal{P}}_{l}$ represent the average power that CHl spends to send its fused signal $y_{l}$ to the FC. We have:

[TABLE]

Applying (3) and noting that $\boldsymbol{x}_{l},\boldsymbol{q}_{l}$ in (3) are zero mean and uncorrelated, it is easy to verify that:

[TABLE]

where

[TABLE]

Combining (7) and (8) we obtain:

[TABLE]

Let $P_{trn}\!=\!\sum_{l=1}^{L}\psi_{l}$ be the total power that CHs spend to transmit their pilot symbols to the FC for channel estimation. We assume there is a constraint on the network transmit power, such that:

[TABLE]

Substituting ${\mathcal{P}}_{l}$ in (9) into the constraint in (10) we reach:

[TABLE]

II-C Problem Statement

Under the network power constraint in (10), our goal is to find the optimal $P_{trn},\{P_{l},{\cal P}_{l}\}_{l=1}^{L}$ such that $D$ is minimized. The constraint in (11) shows that finding the optimal $\{P_{l},{\cal P}_{l}\}_{l=1}^{L}$ in our problem is equivalent to finding the optimal $\{P_{l},\boldsymbol{w}_{l}\}_{l=1}^{L}$, since given $\{P_{l},\boldsymbol{w}_{l}\}_{l=1}^{L}$ one can find $\{{\cal P}_{l}\}_{l=1}^{L}$ using (9). Therefore, our goal is to find the optimal total training power $P_{trn}$, the optimal total power that sensors in cluster $l$ spend to transmit their measurements to their CH $P_{l}$, and the optimal $\boldsymbol{w}_{l}$ employed by CHl for its linear fusion, such that $D$ is minimized. In other words, we are interested in solving the following constrained optimization problem:

[TABLE]

We note that $\boldsymbol{\Sigma}_{q_{l}}$ and $\boldsymbol{\Omega}_{l}$ in the network transmit power constraint do not depend on our optimization variables.

III Characterizing $D$ and its Lower Bounds

III-A Characterization of $D$ in terms of Channel Estimates

We characterize the objective function $D$ in (12), in terms of our optimization variables. Before delving into the derivations of $D$, we introduce the following notations. Considering our signal model in Section II, we define column vectors $\boldsymbol{x}$, $\boldsymbol{t}$, $\boldsymbol{y}$, $\boldsymbol{z}$ in Table I, which are obtained from stacking the signals corresponding to all clusters. We have:

[TABLE]

where $\boldsymbol{1}$ is a column vector of $K$ ones, column vectors $\boldsymbol{n}$, $\boldsymbol{q}$, $\boldsymbol{v}$, and matrices $\boldsymbol{M}$, $\boldsymbol{W}$, $\boldsymbol{H}$ are defined in Table I. The noise vectors $\boldsymbol{n}$, $\boldsymbol{q}$, $\boldsymbol{v}$ are zero-mean with covariance matrices $\boldsymbol{\Sigma}_{n}$, $\boldsymbol{\Sigma}_{q}$, $\boldsymbol{\Sigma}_{v}$, respectively, given in Table I. We model the fading coefficient as $h_{l}\!=\!\hat{h}_{l}+\tilde{h}_{l}$, where $\hat{h}_{l}$ is the MMSE channel estimate and $\tilde{h}_{l}$ is the corresponding zero-mean estimation error with the variance $\zeta_{l}^{2}$. The expressions for $\hat{h}_{l}$ and $\zeta_{l}^{2}$ in terms of training power $\psi_{l}$ are [30]:

[TABLE]

We define matrices $\hat{\boldsymbol{H}},\tilde{\boldsymbol{H}}$ in Table I and thus we have $\boldsymbol{H}\!=\!\hat{\boldsymbol{H}}\!+\!\tilde{\boldsymbol{H}}$. Substituting this channel model into (13b), we can rewrite the received signal $\boldsymbol{z}$ as the following:

[TABLE]

We proceed with characterizing $D$ in terms of the channel estimates. From optimal linear estimation theory, we have:

[TABLE]

where $\hat{\theta}$ and $D$ depend on the channel estimates $\{\hat{h}_{l}\}_{l=1}^{L}$. In the following, we find $\mathbb{E}\{\boldsymbol{z}{\boldsymbol{z}}^{H}\}$ and $\mathbb{E}\{\theta\boldsymbol{z}\}$ in (III-A) by examining the statistics of channel estimation error. By the orthogonality principle of LMMSE estimation [31], $\tilde{h}_{l}$ is orthogonal to $\hat{h}_{l}$, that is $\mathbb{E}\{\tilde{h}_{l}\hat{h}_{l}\}\!=\!0,\ \forall l$, and therefore, $\mathbb{E}\{\boldsymbol{z}_{1}{\boldsymbol{z}_{2}}^{H}\}\!=\!\boldsymbol{0}$. Using the fact that $\theta$, $\boldsymbol{n}$, $\boldsymbol{q}$, $\boldsymbol{v}$ are mutually uncorrelated, we have $\mathbb{E}\{\boldsymbol{z}_{1}{\boldsymbol{z}_{3}}^{H}\}\!=\!\boldsymbol{0},\ \mathbb{E}\{\boldsymbol{z}_{2}{\boldsymbol{z}_{3}}^{H}\}\!=\!\boldsymbol{0}$. Combined these with the fact that $\mathbb{E}\{\boldsymbol{z}\}\!=\!\boldsymbol{0}$, the covariance matrix $\boldsymbol{C}_{\boldsymbol{z}}\!=\!\mathbb{E}\{\boldsymbol{z}{\boldsymbol{z}}^{H}\}$ given $\hat{\boldsymbol{H}}$ can be expressed as:

[TABLE]

where $\boldsymbol{\Gamma}$ and $\boldsymbol{\Sigma}$ are defined in Table I and $\boldsymbol{\Sigma}_{l}$ is a $K_{l}\times K_{l}$ matrix of all ones. We define $\boldsymbol{\Lambda}_{1_{l}}$ and $\boldsymbol{\Lambda}_{2}$ as bellow:

[TABLE]

where $\boldsymbol{\mu}$ is defined in Table I. It is straightforward to simplify (III-A) and write it as the following:

[TABLE]

where $\boldsymbol{\Lambda}_{1}$ is defined in Table I. To find $\mathbb{E}\{\theta\boldsymbol{z}\}$ we consider (15) and we realize that $\mathbb{E}\{\theta\boldsymbol{z}_{3}\}\!=\!\boldsymbol{0}$. Therefore:

[TABLE]

where ($a$) in (20) is obtained from the fact that $\mathbb{E}\{\tilde{\boldsymbol{H}}\}\!=\!\boldsymbol{0}$. Based on (19), (20), the LMMSE estimator $\hat{\theta}$ and its corresponding MSE in (III-A) can be written as:

[TABLE]

in which $\boldsymbol{\mu}$ and $\boldsymbol{C}_{\boldsymbol{z}}$ depend on the channel estimates. Substituting (19) in (III-A) and applying the Woodbury identity666For matrices $\boldsymbol{A},\boldsymbol{B},\boldsymbol{C}$, $\boldsymbol{D}$ the Woodbury identity states that ${(\boldsymbol{A}+\boldsymbol{B}\boldsymbol{C}\boldsymbol{D})}^{-1}\!=\!{\boldsymbol{A}}^{-1}-{\boldsymbol{A}}^{-1}\boldsymbol{B}{({\boldsymbol{C}}^{-1}+\boldsymbol{D}{\boldsymbol{A}}^{-1}\boldsymbol{B})}^{-1}\boldsymbol{D}{\boldsymbol{A}}^{-1}$ [32]. yields:

[TABLE]

Examining $D$ in (22) we notice that $\boldsymbol{\Pi}_{l}$ does not depend on our optimization variables. However, $\boldsymbol{\Lambda}_{1_{l}}$ depends on $P_{l}$ and $\psi_{l}$ (through the channel estimate ${|\hat{h}_{l}|}^{2}$ and the channel estimation error variance $\zeta_{l}^{2}$). Clearly, $D$ depends on $\boldsymbol{w}_{l}$.

III-B Three Lower Bounds on $D$

We provide three lower bounds on $D$, denoted as $D_{1},D_{2},D_{3}$. To obtain $D_{1}$ we consider the scenario when $\{h_{l}\}_{l=1}^{L}$ are available at the FC (perfect CSI). This implies $\hat{h}_{l}\!=\!h_{l}$ and $\zeta_{l}^{2}\!=\!0,\forall l$, in (22), and the MSE becomes:

[TABLE]

To obtain $D_{2}$ we consider the scenario when in addition to perfect CSI, sensors’ noisy measurement vector $\boldsymbol{x}_{l}$ is available at CHl (i.e., error-free channels between sensors and their CHs). Therefore, $\boldsymbol{A}_{l}\!=\!\boldsymbol{I}_{l}$, where $\boldsymbol{I}_{l}$ denotes the identity matrix, and $\boldsymbol{\Sigma}_{q_{l}}\!=\!\boldsymbol{0},\forall l$. In this scenario (23) simplifies to:

[TABLE]

To obtain $D_{3}$ we consider the scenario when $\boldsymbol{x}_{l}$ is available at CHl and $y_{l}$ is available at the FC. This is equivalent to having all measurements $\{\boldsymbol{x}_{l}\}_{l=1}^{L}$ available at the FC (i.e., error-free channels between sensors and their CHs, and between CHs and the FC). Therefore, the MSE becomes:

[TABLE]

Clearly, we have $D_{3}<D_{2}<D_{1}<D$.

III-C Bayesian CRB

Let $G$ denote the Bayesian Fisher information corresponding to estimating $\theta$, given $\boldsymbol{z}$ and the vector of channel estimates $\hat{\boldsymbol{h}}\!=\![\hat{h}_{1},...,\hat{h}_{L}]$ at the FC. The inverse of $G$ is the Bayesian CRB and it sets an estimation-theoretic lower bound on the MSE of any Bayesian estimation of $\theta$, given $\boldsymbol{z},\hat{\boldsymbol{h}}$ [33, 34, 35]. Using the definition in [33, 34, 35] in our problem $G\!=\!\mathbb{E}\{{\!(\frac{\partial\ln f(\boldsymbol{z},\hat{\boldsymbol{h}},\theta)}{\partial\theta})\!}^{2}\}$, where $f(\boldsymbol{z},\hat{\boldsymbol{h}},\theta)$ denotes the joint pdf of $\boldsymbol{z},\hat{\boldsymbol{h}},\theta$ and and the expectation is taken over $f(\boldsymbol{z},\hat{\boldsymbol{h}},\theta)$.

Lemma 1**.**

The Bayesian Fisher information corresponding to estimating $\theta$, given $\boldsymbol{z},\hat{\boldsymbol{h}}$ is:

[TABLE]

where $G_{1}(\theta)=-\frac{\partial^{2}\ln f(\theta)}{\partial\theta^{2}}$ and $G_{2}(\theta)$ is given below. For $\theta\sim N(0,\sigma_{\theta}^{2})$ we have $\mathbb{E}\{G_{1}(\theta)\}=\sigma_{\theta}^{-2}$. Both expectations in (26) are taken over $f(\theta)$, which represents the pdf of $\theta$.

[TABLE]

where $f(z_{l}|\hat{h}_{l},\theta)$ and its derivative with respect to $\theta$ are:

[TABLE]

[TABLE]

and the parameters $a_{1},a_{2},a_{3},c_{m,n,p}(\theta),s_{m,n,p,b}(\theta),\bar{\phi}$ are:

[TABLE]

Proof.

See Appendix -A. ∎

IV Solving the Constrained Minimization of $D$

We consider the constrained optimization problem in (12), where $D$ is provided in (22). We define:

[TABLE]

Using the two definitions in (31) we can replace the problem in (12) with its equivalent, problem (P1), that has a simpler presentation. In particular, we can write $D^{-1}\!=\!\sigma^{-2}_{\theta}+\sum_{l=1}^{L}{\cal J}_{l}(P_{trn},P_{l},\boldsymbol{w}_{l})$. Hence, problem (P1) becomes:

[TABLE]

It is easy to show that the solution of (P1) holds with active constraint $P_{trn}+\sum_{l=1}^{L}C_{l}(P_{l},\boldsymbol{w}_{l})=P_{tot}$. We further note that due to the cap on the network transmit power, only a subset of the clusters may become active at each observation period. We refer to this active subset as ${\cal A}\!=\!\{l:\!P_{l}\!>\!0,\ l\!=\!1,\dots,L\}$, where $|{\cal A}|\leq L$. Regarding the objective function ${\cal J}_{l}$ in (P1) we note that it depends on $\hat{h}_{l}$ (through ${|\hat{h}_{l}|}^{2}$ in the numerator and $\boldsymbol{\Lambda}_{1_{l}}$ in the denominator of (31)). Regarding the optimization variables in (P1) we notice that, since pilot transmission proceeds data transmission, $P_{trn}$ cannot depend on the channel estimates $\{\hat{h}_{l}\}_{l=1}^{L}$ and can only depend on the statistical information of communication channels and the observation model. Examining (P1), we note however, that solving it for $P_{trn}$ provides an answer that depends on $\hat{h}_{l}$ (which is unrealizable). On the other hand, the variables $P_{l},\boldsymbol{w}_{l}$ should be chosen according to the available CSI $\hat{h}_{l}$. Based on these observations, we propose to consider two problems (PA) and (PB) stemming from (P1). problem (PA) finds the optimal $\{P_{l},\boldsymbol{w}_{l}\}_{l=1}^{L}$ that minimizes $D$, given $P_{trn}$. Let $\sigma\in(0,1)$ such that $P_{trn}\!=\!(1-\sigma)P_{tot}$. Given $P_{trn}$ (and thus $\sigma$), we define ${\cal F}_{l}(P_{l},\boldsymbol{w}_{l})\!=\!{\cal J}_{l}(P_{trn},P_{l},\boldsymbol{w}_{l})$. Problem (PA) becomes:

[TABLE]

Section IV-A is devoted to solving (PA). Problem (PB) finds the optimal $P_{trn}$ that, instead of minimizing $D$, it minimizes a modified objective function $\mathbb{E}\{D\}$, where an average is taken over the channel estimates. In Section IV-B we address (PB) and find $P_{trn}$ as well as training power distribution $\{\psi_{l}\}_{l=1}^{L}$ among the CHs such that $\sum_{l=1}^{L}\psi_{l}\!=\!P_{trn}$.

IV-A Finding Optimal $\{P_{l},\boldsymbol{w}_{l}\}_{l=1}^{L}$ Given Total Training Power

We start with (PA). By taking the second derivative of $\sum_{l=1}^{L}{\cal F}_{l}(P_{l},\boldsymbol{w}_{l})$ w.r.t $\{P_{l},\boldsymbol{w}_{l}\}$, it is straightforward to show that (PA) is not jointly concave over the optimization variables. Alternatively, we propose a solution approach that converges to a stationary point of (PA). Problem (PA) contains the constraint $\sum_{l=1}^{L}C_{l}(P_{l},\boldsymbol{w}_{l})\leq\sigma P_{tot}$, which is referred to as coupling or complicating constraint in the literature [36]. By introducing additional auxiliary variables $\{{\cal V}_{l}\}_{l=1}^{L}$, problem (PA) becomes:

[TABLE]

Note that the auxiliary variable ${\cal V}_{l}$ represents the total amount of power allocated to cluster $l$ (for sensors within cluster $l$ to transmit their observations to CHl and for CHl to transmit $y_{l}$ to the FC). According to the primal decomposition [36], problem (P2) can be decomposed as the following:

[TABLE]

[TABLE]

where ${\cal F}^{opt}_{l}$ denotes the maximum of ${\cal F}_{l}(P_{l},\boldsymbol{w}_{l})$, which depends on ${\cal V}_{l}$. The solution can be reached by iteratively solving sub-problems (SP2-1) and (SP2-2). In the following, we provide the detailed solutions for (SP2-1) and (SP2-2).

IV-A1 Solving Optimization Problem (SP2-1)

We start with a brief overview of this section. Let ${\boldsymbol{w}}^{opt}_{l},P^{opt}_{l}$ denote the solution of (SP2-1). We will show how to compute ${\boldsymbol{w}}^{opt}_{l}$ in terms of $P_{l}$ using (42) and how to compute $P^{opt}_{l}$ in terms of $\boldsymbol{w}_{l}$ using (47). Having two equations (42), (47), we substitute $\boldsymbol{w}_{l}$ from (42) into (47) to reach (IV-A1), which is a function of $P^{opt}_{l}$ only. Employing a numerical line search method we obtain $P^{opt}_{l}$ from (IV-A1). Having $P^{opt}_{l}$, we find ${\boldsymbol{w}}^{opt}_{l}$ using (42). The detailed explanations follow.

Examining ${\cal F}_{l}(P_{l},\boldsymbol{w}_{l})$ and $C_{l}(P_{l},\boldsymbol{w}_{l})$ expressions given in (31), it is evident that scaling up equally $P_{l},\boldsymbol{w}_{l}$ increases both ${\cal F}_{l}(P_{l},\boldsymbol{w}_{l})$ and $C_{l}(P_{l},\boldsymbol{w}_{l})$. Therefore, (SP2-1) is equivalent to its converse formulation, where $C_{l}(P_{l},\boldsymbol{w}_{l})$ is minimized subject to a constraint on ${\cal F}_{l}(P_{l},\boldsymbol{w}_{l})$:

[TABLE]

Let $C^{opt}_{l}$ be the minimum of $C_{l}(P_{l},\boldsymbol{w}_{l})$, which depends on ${\cal U}_{l}$. To solve (CSP2-1) we simplify its constraint by substituting $\boldsymbol{\Lambda}_{1_{l}}$ from (18) into ${\cal F}_{l}(P_{l},\boldsymbol{w}_{l})$ in (31). Let $\boldsymbol{B}_{l}\!=\!\sigma^{2}_{\theta}\zeta_{l}^{2}\boldsymbol{\Pi}_{l}+({|\hat{h}_{l}|}^{2}+\zeta_{l}^{2})\boldsymbol{\Delta}_{l}$. The constraint in (CSP2-1) becomes:

[TABLE]

Consider (CSP2-1) where its constraint is now replaced with the inequality in (32). To solve (CSP2-1) we use the Lagrange multiplier method. Let ${\mathcal{L}}(\gamma,\eta,P_{l},\boldsymbol{w}_{l})$ be the Lagrangian for this problem and $\gamma$ and $\eta$ be the lagrange multipliers for the constraint in (32) and the constraint $P_{l}\!\geq\!0$, respectively. Equation (33) shows ${\mathcal{L}}(\gamma,\eta,P_{l},\boldsymbol{w}_{l})$. The corresponding Karush-Kuhn-Tucker (KKT) optimality conditions are [37, pp. 243-244]:

[TABLE]

where $\boldsymbol{R}_{t_{l}}$, defined in (8), depends on $P_{l}$. Similar to the solution of (P1), one can show that the solutions of (SP2-1) and (CSP2-1) must satisfy the equality constraints $C_{l}(P_{l},\boldsymbol{w}_{l})\!=\!{\cal V}_{l}$ and ${\cal F}_{l}(\!P_{l},\boldsymbol{w}_{l}\!)\!=\!{\cal U}_{l}$ (or equivalently (34b)), respectively. Thus we find:

[TABLE]

Combining (35a) and (35b) we reach:

[TABLE]

From (34a) and (36) we find the lagrange multiplier $\gamma$:

[TABLE]

$\bullet$ Computing ${\boldsymbol{w}}^{opt}_{l}$ given $P_{l}$: Substituting (37) into (34a) and conducting some mathematical manipulations result in:

[TABLE]

where $\boldsymbol{\mu}_{l}$ is defined in (18). Since $\boldsymbol{R}_{t_{l}}\!\succ\!\boldsymbol{0},\boldsymbol{\Sigma}_{q_{l}}\!\succ\!\boldsymbol{0},\boldsymbol{B}_{l}\!\succ\!\boldsymbol{0}$, the matrix $\boldsymbol{\cal B}_{1}$ is positive definite and full rank and hence invertible. Multiplying both sides of (38) with ${\boldsymbol{\cal B}_{1}}^{-1}$, we find:

[TABLE]

Also, multiplying both sides of (38) with $\frac{1}{{\cal U}_{l}}{\boldsymbol{R}_{t_{l}}}^{-1}$ we reach:

[TABLE]

Inspecting (39) and (40), and aiming at finding vector $\boldsymbol{w}_{l}$, we realize that (39) and (40) are ordinary eigenvalue problems. Since the solutions to (SP2-1) and (CSP2-1) satisfy the equality constraints $C_{l}(P_{l},\boldsymbol{w}_{l})\!=\!{\cal V}_{l}$ and ${\cal F}_{l}(\!P_{l},\boldsymbol{w}_{l}\!)\!=\!{\cal U}_{l}$, respectively, from (39) and (40) we find:

[TABLE]

Let ${\boldsymbol{s}}^{opt}_{l}$ be the eigenvector corresponding to $\lambda_{max}({\boldsymbol{\cal B}_{1}}^{-1}|\boldsymbol{\mu}_{l}|{|\boldsymbol{\mu}_{l}|}^{T})$. We note that ${\cal F}^{opt}_{l}$ is achieved when $\boldsymbol{w}_{l}$ is an appropriately scaled version of ${\boldsymbol{s}}^{opt}_{l}$, i.e., ${\boldsymbol{w}}^{opt}_{l}\!=\!r_{l}{\boldsymbol{s}}^{opt}_{l}$, where scalar $r_{l}$ is such that (35a) is satisfied. Also recall $\boldsymbol{\Pi}_{l}\!=\!\boldsymbol{\rho}_{l}{\boldsymbol{\rho}_{l}}^{T}$ is rank-1. Thus ${\cal F}^{opt}_{l}$ is the only non-zero eigenvalue of ${\boldsymbol{\cal B}_{1}}^{-1}|\boldsymbol{\mu}_{l}|{|\boldsymbol{\mu}_{l}|}^{T}$ and ${\boldsymbol{s}}^{opt}_{l}$ is the corresponding eigenvector. Proposition 1 gives expressions for ${\boldsymbol{w}}^{opt}_{l}$ and ${\cal F}^{opt}_{l}$ in terms of $P_{l}$.

Proposition 1**.**

Considering problem (SP2-1), the optimal fusion vector ${\boldsymbol{w}}^{opt}_{l}$ and the maximum value of the objective function ${\cal F}^{opt}_{l}$ in terms of $P_{l}$ are:

[TABLE]

where $\tau_{l}\!=\!{\boldsymbol{\rho}_{l}}^{T}{\boldsymbol{R}_{t_{l}}}^{-1}\boldsymbol{\rho}_{l},\ \beta_{l}\!=\!\frac{\sigma_{v_{l}}^{2}}{{|\hat{h}_{l}|}^{2}(1-\sigma^{2}_{\theta}P_{l}\tau_{l})+\zeta_{l}^{2}}$.

**

Proof.

See Appendix -B. ∎

For our system model $\boldsymbol{R}_{t_{l}\theta}\!=\!\mathbb{E}\{\theta{\boldsymbol{t}_{l}}\}\!=\!\sigma^{2}_{\theta}\sqrt{P_{l}}\boldsymbol{\rho}_{l}$. Hence, we can rewrite ${\boldsymbol{w}}^{opt}_{l}$ in (42) as:

[TABLE]

Since $\boldsymbol{R}_{t_{l}}^{-1}\boldsymbol{R}_{t_{l}\theta}$ is the linear operator corresponding to the LMMSE estimator, (43) implies that the optimal linear fusion rule at CHl is equal to the linear operator corresponding to the LMMSE estimation of $\theta$ based on $\boldsymbol{t}_{l}$, multiplied by the amplification factor $\chi_{l}$.

$\bullet$ Computing $P^{opt}_{l}$ given $\boldsymbol{w}_{l}$: Note that (34d) results in $\eta\!=\!0$ for active clusters with $P_{l}>0$. Letting $\eta\!=\!0$ in (34c) and solving for $\gamma$ we find:

[TABLE]

Equating (44) with (37) and solving for ${\cal U}_{l}$ we get:

[TABLE]

On the other hand, solving (35b) for ${\cal U}_{l}$ results in:

[TABLE]

Combining (45) and (46), we obtain $P^{opt}_{l}$ in terms of ${\boldsymbol{w}}_{l}$ as the following:

[TABLE]

At this point, we have obtained two equations: (42) provides ${\boldsymbol{w}}^{opt}_{l}$ in terms of $P_{l}$, and (47) provides $P^{opt}_{l}$ in terms of $\boldsymbol{w}_{l}$. Substituting ${\boldsymbol{w}}^{opt}_{l}$ from (42) in (47) yields in:

[TABLE]

Note that $\tau_{l},\!\boldsymbol{R}_{t_{l}}$ in (IV-A1) depend on $P^{opt}_{l}\!\!$, and thus, a closed-form solution for $P^{opt}_{l}$ remains elusive. One can employ a line search method (e.g., the Golden section method [38, p. 216]) to solve (IV-A1) in the interval $(0,{\cal V}_{l})$. Having $P^{opt}_{l}\!$ we find ${\boldsymbol{w}}^{opt}_{l}\!$ using (42).

IV-A2 Solving Optimization Problem (SP2-2)

By substituting ${\cal F}^{opt}_{l}$ from (42) in the objective function, problem (SP2-2) becomes:

[TABLE]

The maximization problem in (IV-A2) is concave and its solution can be found via solving the KKT conditions. In particular, we find (see Appendix -C for derivations):

[TABLE]

Note that the first term of the right side of the equality in (50a) is ${\cal P}_{l}$ introduced in Section II-B. Given $\lambda,|\hat{h}_{l}|,\sigma_{v_{l}}$ in (50a) and the easy-to-prove fact that $\tau_{l}\!+\!P_{l}\frac{\partial\tau_{l}}{\partial P_{l}}\!>\!0$, it is straightforward to show that $\frac{\partial{\cal P}_{l}}{\partial P_{l}}\!>\!0$ for active clusters, i.e., increasing $P_{l}$ increases ${\cal P}_{l}$. Having the solutions to problems (SP2-1) and (SP2-2), Algorithm 1 summarizes our proposed solution to problem (PA). Essentially, this algorithm iteratively solves (SP2-1) and (SP2-2) in a block-coordinate ascent manner until the convergence is reached. In Section VII, we argue that the algorithm output converges to a stationary point of (PA).

IV-B Finding Optimal Total Training Power and its Distribution Among CHs

In this section, we focus on (PB) and find $P_{trn}$ as well as training power distribution $\{\psi_{l}\}_{l=1}^{L}$ among the CHs such that $\sum_{l=1}^{L}\psi_{l}\!=\!P_{trn}$. As we mentioned earlier, to find $P_{trn}$ we consider a modified objective function, i.e., instead of $\sum_{l=1}^{L}{\cal J}_{l}$ in (P1) we consider $\sum_{l=1}^{L}\mathbb{E}\{{\cal J}_{l}\}$, where the expectation is taken over the channel estimates ${|\hat{h}_{l}|}^{2}$. Since solving this problem analytically is still intractable, we use the Jensen’s inequality for concave functions [37, pp. 77-78], to establish a lower bound on $\mathbb{E}\{{\cal J}_{l}(P_{trn},P_{l},\boldsymbol{w}_{l})\}$:

[TABLE]

where ${\cal G}_{l}(P_{trn},P_{l},\boldsymbol{w}_{l})$ is obtained from ${\cal J}_{l}(P_{trn},P_{l},\boldsymbol{w}_{l})$, after replacing ${|\hat{h}_{l}|}^{2}$ with $\mathbb{E}\{{|\hat{h}_{l}|}^{2}\}$. To find $\mathbb{E}\{{|\hat{h}_{l}|}^{2}\}$ needed for ${\cal G}_{l}(P_{trn},P_{l},\boldsymbol{w}_{l})$ we revisit the error corresponding to the LMMSE channel estimation in (14). Note that $\hat{h}_{l}$ is a zero-mean complex Gaussian. Let $2\sigma_{\hat{h}_{l}}^{2}$ denote the variance of $\hat{h}_{l}$. For the model $h_{l}\!=\!\hat{h}_{l}\!+\!\tilde{h}_{l}$, we invoke the orthogonality principle from the linear estimation theory [30], that states $var({\hat{h}_{l}})\!=\!var({h_{l}})-var({\tilde{h}_{l}})\!=\!2\sigma_{h_{l}}^{2}-\zeta_{l}^{2}$, where $\zeta_{l}^{2}$ in (14) depends on $\psi_{l}$. Since $\hat{h}_{l}$ is zero-mean, we have $\mathbb{E}\{{|\hat{h}_{l}|}^{2}\}\!=\!var({\hat{h}_{l}})$. Thus, ${\cal G}_{l}(P_{trn},P_{l},\boldsymbol{w}_{l})\!=\!\frac{(2\sigma_{h_{l}}^{2}-\zeta_{l}^{2})P_{l}{\boldsymbol{w}_{l}}^{T}\boldsymbol{\Pi}_{l}\boldsymbol{w}_{l}}{\sigma_{v_{l}}^{2}+{\boldsymbol{w}_{l}}^{T}\boldsymbol{\Lambda}_{1_{l}}\boldsymbol{w}_{l}}$, where $\boldsymbol{\Lambda}_{1_{l}}\!=\!\sigma^{2}_{\theta}\zeta_{l}^{2}P_{l}\boldsymbol{\Pi}_{l}+2\sigma_{h_{l}}^{2}(\boldsymbol{\Sigma}_{q_{l}}+P_{l}\boldsymbol{\Delta}_{l})$. Notice that ${\cal G}_{l}$ depends on the optimization variable $P_{trn}$ through $\zeta_{l}^{2}$ in the numerator and $\boldsymbol{\Lambda}_{1_{l}}$ in the denominator. We reconsider (P1) in which ${\cal J}_{l}$ is now replaced with ${\cal G}_{l}$:

[TABLE]

Examining (P${}_{B^{\prime}}$), we realize that solving it for $P_{trn}$ provides an answer that depends on ${P}_{l},{\boldsymbol{w}}_{l}$ (which is undesirable). To circumvent this problem we propose a method to find $P_{trn}$ based on the following observation. We observe that, although (P${}_{B^{\prime}}$) is a non-concave maximization problem, given $\{{P}_{l},{\boldsymbol{w}}_{l}\}_{l=1}^{L}$ and letting $\sigma=1-\frac{P_{trn}}{P_{tot}}$, the problem (P${}_{B^{\prime}}$) under these conditions becomes strictly concave with respect to the variable $\sigma$ over the interval $(0,1)$, and hence, the objective function has a unique global maximum in this interval. Let ${\sigma}^{opt}$ denote the solution to this problem, which can be efficiently found using numerical line search methods (e.g., the Golden section777Let $x^{opt}$ denote the maximum value that a concave function $f(x)$ attains over a search interval $x\!\in\!(x_{b},x_{e})$. This numerical method finds $x^{opt}$ via successively narrowing the range of the search interval. Let $i$ be the iteration index, $({\cal I}_{b}^{(i)},{\cal I}_{e}^{(i)})$ be the starting and ending points of the search interval at iteration $i$, $\alpha_{b}^{(i)}\!=\!0.382({\cal I}_{e}^{(i)}-{\cal I}_{b}^{(i)})+{\cal I}_{b}^{(i)}$ and $\alpha_{e}^{(i)}\!=\!0.618({\cal I}_{e}^{(i)}-{\cal I}_{b}^{(i)})+{\cal I}_{b}^{(i)}$ be the evaluating points. Also let $f_{b}^{(i)},f_{e}^{(i)}$, denote the values of the function $f(x)$ when it is evaluated at the evaluating points $\alpha_{b}^{(i)},\alpha_{e}^{(i)}$, respectively. For initialization, we let $i\!=\!0,{\cal I}_{b}^{(0)}\!=\!x_{b},{\cal I}_{e}^{(0)}\!=\!x_{e}$. At iteration $i$, we compute $f_{b}^{(i)}$ and $f_{e}^{(i)}$ and then update the search interval to find $x^{opt}$ as the following: if $f_{b}^{(i)}\!>\!f_{e}^{(i)}$, then ${\cal I}_{b}^{(i+1)}\!=\!{\cal I}_{b}^{(i)},{\cal I}_{e}^{(i+1)}\!=\!\alpha_{e}^{(i)}$, if $f_{b}^{(i)}\!=\!f_{e}^{(i)}$, then ${\cal I}_{b}^{(i+1)}\!=\!\alpha_{b}^{(i)},{\cal I}_{e}^{(i+1)}\!=\!\alpha_{e}^{(i)}$, and if $f_{b}^{(i)}\!<\!f_{e}^{(i)}$, then ${\cal I}_{b}^{(i+1)}\!=\!\alpha_{b}^{(i)},{\cal I}_{e}^{(i+1)}\!=\!{\cal I}_{e}^{(i)}$. As the stopping criterion, we check whether the length of the search interval exceeds a pre-determined threshold $\epsilon$. If the stopping criterion is met at iteration $j$, the algorithm returns the optimal solution $x^{opt}\!=\!{\cal I}_{b}^{(j)}$. Otherwise, we update the search interval and continue the iterations until the stopping criterion is met. method [38, p. 216]). Since this problem is concave over $(0,1)$, the convergence of Golden section method to ${\sigma}^{opt}$ is guaranteed.

Based on the above observation, we propose the method described in Algorithm 2 to solve (P${}_{B^{\prime}}$) and find $P_{trn}^{opt}$. The proposed method is basically Golden section method, where in each iteration we apply Algorithm 1 to find $\{{P}_{l},{\boldsymbol{w}}_{l}\}_{l=1}^{L}$, only for the purpose of successively narrowing the search interval for $\sigma$. The output of Algorithm 2 converges to ${\sigma}^{opt}$ and thus $P_{trn}^{opt}=(1-{\sigma}^{opt})P_{tot}$.

Given $P_{trn}^{opt}$, we find $\{\psi_{l}\}_{l=1}^{L}$, via minimizing the MSE of the LMMSE channel estimates for all clusters:

[TABLE]

The above is a convex minimization problem. Solving the associated KKT conditions, we obtain:

[TABLE]

The solution in (52) is based on the assumption that all CHs participate in pilot transmission and $P_{trn}^{opt}$ satisfies the inequality $P_{trn}^{opt}\!\geq\!\frac{\sigma_{v_{L}}}{\sigma_{h_{L}}^{2}}\sum_{l=1}^{L}\sigma_{v_{l}}\!-\!\sum_{l=1}^{L}\frac{\sigma_{v_{l}}^{2}}{\sigma_{h_{l}}^{2}}\!=\!\Upsilon$. However, when $P_{trn}^{opt}\!<\!\Upsilon$, the solutions in (52) imply that $\psi_{l}\!=\!0$ for some clusters. In this case, we propose to choose $\psi_{l}\!=\!a\frac{\sigma_{h_{l}}^{2}}{\sigma_{v_{l}}}$, in which $a$ is a common factor. Imposing the constraint $\sum_{l=1}^{L}\psi_{l}\!=\!P_{trn}^{opt}$ results in:

[TABLE]

Fig. 2 shows a block diagram that summarizes our approach to solve the original constrained optimization problem (P1). Overall, the sequence of algorithm implementations and network operation follow. The FC implements Algorithm 2 to obtain $P_{trn}^{opt}$, and consequently to find $\{\psi_{l}\}_{l=1}^{L}$ given in (52). The FC feeds back this information to CHs (all the obtained $\{{P}_{l},{\boldsymbol{w}}_{l}\}_{l=1}^{L}$ values during the execution of Algorithm 2 are discarded at this point). CHs send their pilot symbols to the FC and the FC estimates the channels $\{\hat{h}_{l}\}_{l=1}^{L}$. Now, given $P_{trn}^{opt},\{\hat{h}_{l}\}_{l=1}^{L}$, the FC implements Algorithm 1, finds $\{P_{l}^{opt},\boldsymbol{w}_{l}^{opt}\}_{l=1}^{L}$, feeds back888Similar to [24, 25] we assume that the FC energy resource is much larger than those of the sensors/CHs. Therefore, the overhead required for feeding back the necessary information from the FC to the sensors/CHs is neglected. this new information to CHs, and feeds back $P_{l,k}=\frac{P_{l}^{opt}}{K_{l}}$ to sensors. Sensors send their amplified measurements to their CHs. CHs send their fused signals to the FC. Finally, the FC estimates $\theta$.

IV-C Minimizing Lower Bounds on MSE $D$

This section discusses constrained minimization of the lower bounds $D_{1},D_{2}$ we derived in Section III-B. The lower bound $D_{1}$ depends on $\{P_{l},\boldsymbol{w}_{l}\}_{l=1}^{L}$ and hence its constrained minimization becomes:

[TABLE]

This is similar to (P2), with the difference that $P_{trn}\!=\!0$, and hence in (50a) and (50b) expressions we let $\zeta_{l}^{2}\!=\!0,{|\hat{h}_{l}|}^{2}\!=\!{|h_{l}|}^{2},\sigma\!=\!1$. Algorithm 1 can be followed to find the solution to (P3), using $\boldsymbol{w}_{l}^{opt}$ in (42). The lower bound $D_{2}$ depends on $\{\boldsymbol{w}_{l}\}_{l=1}^{L}$ and hence its constrained minimization becomes:

[TABLE]

This is similar to (P2), with the differences that $P_{trn}\!=\!0$ and $P_{l}\!=\!0,\forall l$. Following similar steps we took in Section IV-A to solve (P2), we find that (50a) and (50b) become:

[TABLE]

The optimal weight vector $\boldsymbol{w}_{l}^{opt}$ corresponding to the solution of (P4) is computed as ${\boldsymbol{w}}^{opt}_{l}\!=\!\sqrt{\frac{{\cal V}_{l}}{\tau_{l}^{{}^{\prime}}}}{(\sigma^{2}_{\theta}\boldsymbol{\Sigma}_{l}\!+\!\boldsymbol{\Sigma}_{n_{l}})}^{-1}\boldsymbol{1}_{l}$.

V Solving the Special Cases of the Original Problem

The original problem (P1) aims at constrained minimization of $D$, with respect to three sets of optimization variables: $P_{trn}$ total training power, $P_{l}$ power allocated to sensors in cluster $l$ to send their measurements to CHl, and ${\cal P}_{l}$ power allocated to CHl to transmit its signal to the FC. To untangle the performance gain that optimizing each set of these optimization variables provides, we consider the following three special cases of (P1). In problem (P1-SC1) assuming $P_{trn}$ is given and $\psi_{l}\!=\!P_{trn}/L$, we optimize $\{P_{l},{\cal P}_{l}\}_{l=1}^{L}$. In problem (P1-SC2) assuming $P_{l}\!=\!P,\forall l$, we optimize $P_{trn},P,\{{\mathcal{P}}_{l}\}_{l=1}^{L}$. In problem (P1-SC3) assuming ${\mathcal{P}}_{l}\!=\!{\mathcal{P}},\forall l$, we optimize $P_{trn},{\mathcal{P}},\{P_{l}\}_{l=1}^{L}$. Note that problem (P1-SC1) is the same as problem (PA) addressed in Section IV-A. In the following we address problems (P1-SC2) and (P1-SC3).

V-A Solving Special Case (P1-SC2): When Intra-Cluster Powers of all Clusters are Equal

Problem (P1-SC2) becomes:

[TABLE]

where $\boldsymbol{\Lambda}_{1_{l}}\!=\!\sigma^{2}_{\theta}\zeta_{l}^{2}P\boldsymbol{\Pi}_{l}+({|\hat{h}_{l}|}^{2}+\zeta_{l}^{2})(\boldsymbol{\Sigma}_{q_{l}}+P\boldsymbol{\Delta}_{l})$. To address (P1-SC2) we consider the following two sub-problems: ($a$) finding $P^{*},\{\boldsymbol{w}_{l}^{*}\}_{l=1}^{L}$ given $P_{trn}$, ($b$) finding $P_{trn}^{*}$ as well as $\{\psi_{l}^{*}\}_{l=1}^{L}$ such that $\sum_{l=1}^{L}\psi_{l}^{*}\!=\!P_{trn}^{*}$. Sub-problem ($a$) is a special case of (PA) in which, for finding $P^{*}$, we use Golden section method, and sub-problem ($b$) is similar to (PB). Recall that $P_{trn}\!=\!(1-\sigma)P_{tot}$ and thus $\sum_{l=1}^{L}(P\!+\!{\cal P}_{l})\!=\!\sigma P_{tot}$. We let $\sigma_{c}\in(0,1)$ such that $P\!=\!(1-\sigma_{c})\sigma P_{tot}$. It is simple to show that sub-problems ($a$) and ($b$) are both concave and hence $P^{*}$ and $P_{trn}^{*}$ are unique. Next, we summarize our proposed solutions for solving sub-problems ($a$) and ($b$) in Algorithms 3-a and 3-b, respectively.

Description of Algorithm 3-a: Let $P^{*}\!=\!(1-\sigma_{c}^{*})\sigma P_{tot}$ denote the optimal $P$. We apply Golden section method to find $\sigma_{c}^{*}\in(0,1)$ and thus $P^{*}$ that maximizes the objective function in (P1-SC2), denoted as ${\cal F}(\sigma_{c})$. At iteration $i$, for each evaluating point we first compute the optimal ${\cal V}_{l}^{(i)}$, denoted as $\{\bar{{\cal V}}_{l}^{(i)}\}_{l=1}^{L}$ using (50a), and substitute $\bar{{\cal V}}_{l}^{(i)}$ into (42) to obtain $\{\bar{\boldsymbol{w}}_{l}^{(i)}\}_{l=1}^{L}$. Next we compute ${\cal F}_{b}^{(i)}$ and ${\cal F}_{e}^{(i)}$. The stopping criterion is similar to Algorithm 2. Algorithm 3-a returns the optimal $\sigma_{c}^{*},\{\boldsymbol{w}_{l}^{*}\}_{l=1}^{L}$.

Description of Algorithm 3-b: We address sub-problem $(b)$ similar to problem (PB) in Section IV-B. More specifically, we consider problem (P${}_{B^{\prime}}$), where $P_{l}$ is substituted by $P$, and apply a modified version of Algorithm 2 to solve it. In particular, at iteration $i$ of Algorithm 2, we use Algorithm 3-a to obtain the optimal variables $\bar{P}^{(i)},\{\bar{\boldsymbol{w}}_{l}^{(i)}\}_{l=1}^{L}$, and then compute ${\cal R}_{b}^{(i)}$ and ${\cal R}_{e}^{(i)}$. The rest is similar to Algorithm 2. Algorithm 3-b returns the optimal $P_{trn}^{*},\{\psi_{l}^{*}\}_{l=1}^{L}$.

V-B Solving Special Case (P1-SC3): When Powers of all CHs for Their Data Transmission to the FC are Equal

To incorporate the constraint ${\mathcal{P}}_{l}\!=\!{\mathcal{P}}$ in the cost function of problem (P1-SC3), from Section IV-A1 we recall that $\boldsymbol{w}_{l}^{opt}\!=\!\chi_{l}({\boldsymbol{R}_{t_{l}}}^{-1}\sigma^{2}_{\theta}\sqrt{P_{l}}\boldsymbol{\rho}_{l})$. Therefore from ${\mathcal{P}}_{l}\!=\!{\boldsymbol{w}_{l}}^{T}\boldsymbol{R}_{t_{l}}\boldsymbol{w}_{l}$ in (9) and ${\mathcal{P}}_{l}\!=\!{\mathcal{P}}$, we conclude $\chi_{l}^{2}\!=\!{\mathcal{P}}/{\sigma^{4}_{\theta}P_{l}\tau_{l}}$. Substituting for $\boldsymbol{w}_{l}$ in (P1), problem (P1-SC3) becomes:

[TABLE]

To address (P1-SC3) we consider the following two sub-problems: ($a$) finding ${\cal P}^{*},\{P_{l}^{*}\}_{l=1}^{L}$ given $P_{trn}$, ($b$) finding $P_{trn}^{*}$ as well as $\{\psi_{l}^{*}\}_{l=1}^{L}$ such that $\sum_{l=1}^{L}\psi_{l}^{*}\!=\!P_{trn}^{*}$. Sub-problem ($a$) is a special case of (PA) in which, for finding ${\cal P}^{*}$, we use Golden section method, and sub-problem ($b$) is similar to (PB). We let $\sigma_{d}\in(0,1)$ such that ${\cal P}\!=\!(1-\sigma_{d})\sigma P_{tot}$. It is easy to show that finding ${\cal P}^{*},P_{trn}^{*}$ in sub-problems ($a$) and ($b$), respectively, are concave problems, and hence ${\cal P}^{*}$ and $P_{trn}^{*}$ are unique. In Appendix -D, we prove that finding $\{P_{l}^{*}\}_{l=1}^{L}$ in sub-problem ($a$) is jointly concave over $P_{l}$’s and therefore its solution is unique. In the absence of a closed form expression we use gradient-ascent algorithm to find the solution. Algorithms 4-a and 4-b summarize how we solve sub-problems ($a$) and ($b$), respectively.

Description of Algorithm 4-a: Let ${\cal P}^{*}\!=\!(1-\sigma_{d}^{*})\sigma P_{tot}$ denote the optimal ${\cal P}$. We apply Golden section method to find $\sigma_{d}^{*}\in(0,1)$ and thus ${\cal P}^{*}$ that maximizes the objective function in (P1-SC3), denoted as ${\cal F}(\sigma_{d})$. At iteration $i$, for each evaluating point we compute the optimal $P_{l}^{(i)}$, denoted as $\{\bar{P}_{l}^{(i)}\}_{l=1}^{L}$ using gradient-ascent algorithm, and substitute them in (P1-SC3) to compute ${\cal F}_{b}^{(i)}$ and ${\cal F}_{e}^{(i)}$. The stopping criterion is similar to Algorithm 2. Algorithm 4-a returns the optimal $\sigma_{d}^{*},\{P_{l}^{*}\}_{l=1}^{L}$.

Description of Algorithm 4-b: We address sub-problem $(b)$ similar to problem (PB) in Section IV-B. Specifically, we consider problem (P${}_{B^{\prime}}$), where ${\cal P}_{l}$ is substituted by ${\cal P}$ and apply a modified version of Algorithm 2 to solve it. In particular, at iteration $i$ of Algorithm 2, we use Algorithm 4-a to obtain the optimal variables $\bar{\cal P}^{(i)},\{\bar{P}_{l}^{(i)}\}_{l=1}^{L}$, and then compute ${\cal R}_{b}^{(i)}$ and ${\cal R}_{e}^{(i)}$. The rest is similar to Algorithm 2. Algorithm 4-b returns the optimal $P_{trn}^{*},\{\psi_{l}^{*}\}_{l=1}^{L}$.

VI Complexity of Algorithms

We discuss the computational complexity of Golden section method as well as Algorithms 1, 2, 3-a, 3-b, 4-a, 4-b, which allows us to compare the computational complexity of solving (P1) versus those of (P1-SC1), (P1-SC2), (P1-SC3).

$\bullet$ Golden section method: This method includes a one-dimensional search to find the optimal point. If no matrix inversion is required, its complexity order for convergence to an $\epsilon$-accurate solution is $\bar{\epsilon}$, where $\bar{\epsilon}\!=\!\log(1/\epsilon)$ [38, p. 217]. We use this method for solving (IV-A1). In each iteration, to compute the left side of (IV-A1) we employ the matrix inversion algorithm in [39] to calculate ${\boldsymbol{R}_{t_{l}}}^{-1}$ with complexity order of ${\cal O}(K_{l}^{2.37})$. Therefore, the overall complexity order of finding $P_{l}^{opt}\!\in\!(0,{\cal V}_{l})$ becomes ${\cal O}(\bar{\epsilon}K_{l}^{2.37})$.

$\bullet$ Algorithm 1 for solving (PA): We switch between solving (SP2-1) and (SP2-2) until the stopping criteria is met. In each iteration, we need to (i) find $\{P_{l}\}_{l=1}^{L}$ using Golden section method, with the overall complexity order of ${\cal O}(\bar{\epsilon}\bar{K})$, where $\bar{K}\!=\!\sum_{l=1}^{L}K_{l}^{2.37}$, and (ii) calculate $\{{\cal V}_{l}\}_{l=1}^{L}$ using (50), which needs $\tau_{l},\beta_{l}$ that are found in (i) and hence, the complexity order of finding $\{{\cal V}_{l}\}_{l=1}^{L}$ is ${\cal O}(L)$. The overall complexity order of Algorithm 1 becomes ${\cal O}(\bar{\epsilon}(L\!+\!\bar{\epsilon}\bar{K}))$.

$\bullet$ Algorithm 2 for solving (P${}_{B^{\prime}}$): In each iteration, for each evaluating point we use Algorithm 1 to obtain $\{P_{l},\boldsymbol{w}_{l}\}_{l=1}^{L}$. Therefore, the overall complexity order of Algorithm 2 becomes ${\cal O}({\bar{\epsilon}}^{2}(L\!+\!\bar{\epsilon}\bar{K}))$.

$\bullet$ Algorithm 3-a for solving sub-problem ($a$) of (P1-SC2): In each iteration, for each evaluating point computing $\tau_{l}$ in (50), (42) involves the matrix inversion ${\boldsymbol{R}_{t_{l}}}^{-1}$, and thus, the complexity order of finding $\{{\cal V}_{l}\}_{l=1}^{L}$ and then $\{\boldsymbol{w}_{l}\}_{l=1}^{L}$ is ${\cal O}(\bar{K})$. Therefore, the overall complexity order of Algorithm 3-a is ${\cal O}(\bar{\epsilon}\bar{K})$.

$\bullet$ Algorithm 3-b for solving sub-problem ($b$) of (P1-SC2): In each iteration, for each evaluating point we use Algorithm 3-a to obtain $P,\{\boldsymbol{w}_{l}\}_{l=1}^{L}$. Therefore, the overall complexity order of Algorithm 3-b is ${\cal O}({\bar{\epsilon}}^{2}\bar{K})$.

$\bullet$ Algorithm 4-a for solving sub-problem ($a$) of (P1-SC3): Note that the complexity order of the gradient-ascent algorithm to maximize a general non-smooth convex function $f(x)$ and converge to an $\epsilon$-accurate solution is ${\cal O}(1/\epsilon)$, if no matrix inversion is required for finding $f(x)$ and its gradient $\triangledown f(x)$ [38, p. 232]. In each iteration of Algorithm 4-a, for each evaluating point, since computing the objective function in (P1-SC3) and its derivative with respect to $P_{l}$ involves the matrix inversion ${\boldsymbol{R}_{t_{l}}}^{-1}$, the complexity order of finding $\{P_{l}\}_{l=1}^{L}$ using the gradient-ascent algorithm is ${\cal O}(\bar{K}/\epsilon)$. Therefore, the overall complexity order of Algorithm 4-a becomes ${\cal O}({\bar{\epsilon}}\bar{K}/\epsilon)$.

$\bullet$ Algorithm 4-b for solving sub-problem ($b$) of (P1-SC3): In each iteration, for each evaluating point we use Algorithm 4-a to obtain ${\cal P},\{P_{l}\}_{l=1}^{L}$. Therefore, the overall complexity order of Algorithm 4-b is ${\cal O}({\bar{\epsilon}}^{2}\bar{K}/\epsilon)$.

To solve (P1) we need to solve (PA), (P${}_{B^{\prime}}$). Therefore, the complexity order of solving (P1) is $e_{0}\!=\!{\cal O}(\bar{\epsilon}(1\!+\!\bar{\epsilon})(L\!+\!\bar{\epsilon}\bar{K}))$. To solve (P1-SC1) we need to solve (PA). Therefore, the complexity order of solving (P1-SC1) is $e_{1}\!=\!{\cal O}(\bar{\epsilon}(L\!+\!\bar{\epsilon}\bar{K}))$. To solve (P1-SC2) we need to solve sub-problems ($a$) and ($b$) of (P1-SC2). Therefore, the complexity order of solving (P1-SC2) is $e_{2}\!=\!{\cal O}(\bar{\epsilon}(1\!+\!\bar{\epsilon})\bar{K})$. To solve (P1-SC3) we need to solve sub-problems ($a$) and ($b$) of (P1-SC3). Therefore, the complexity order of solving (P1-SC3) is $e_{3}\!=\!{\cal O}({\bar{\epsilon}(1\!+\!\bar{\epsilon})}\bar{K}/\epsilon)$ It is clear that $e_{1}\!<\!e_{2}\!<\!e_{0}\!<\!e_{3}$.

VII Convergence Analysis

We discuss the convergence analysis of Algorithms 1 and 2 which solve problems (PA) and (P${}_{B^{\prime}}$), respectively.

$\bullet$ Convergence of Algorithm 1: Problems (PA) and (P2) are equivalent. In (P2), the cost function is non-concave and the constraint is a closed convex set w.r.t. the optimization variables $\{{\cal V}_{l},P_{l},\boldsymbol{w}_{l}\}_{l=1}^{L}$. Algorithm 1 is indeed a block-coordinate ascent type algorithm with two blocks. The first block solves (SP2-1) for all clusters to obtain $\{P_{l},\boldsymbol{w}_{l}\}_{l=1}^{L}$. (SP2-1) is a non-concave maximization problem for which we have a numerical solution for $P_{l}$ using Golden Section method and a closed-form solution for $\boldsymbol{w}_{l}$. Since (SP2-1) is a non-concave maximization problem, we cannot claim that our proposed solution for $P_{l},\boldsymbol{w}_{l}$ is globally optimal and unique. The second block solves (SP2-2) to obtain $\{{\cal V}_{l}\}_{l=1}^{L}$. (SP2-2) is a concave maximization problem for which we have a closed-form solution for ${\cal V}_{l}$. Since (SP2-2) is a concave maximization problem, its solution is globally optimal and unique.

The authors in [40] proved that in a block-coordinate descent algorithm with only two blocks, which solves the unconstrained minimization problem

[TABLE]

given we have the global minimizer $\boldsymbol{x}_{1}^{(k+1)}\!\!=\!\underset{\boldsymbol{x}_{1}}{\text{argmin}}\ f(\boldsymbol{x}_{1},\boldsymbol{x}_{2}^{(k)}),\forall k$, the algorithm converges to a stationary point if we can find a point $\boldsymbol{x}_{2}^{(k+1)}$ such that $f(\boldsymbol{x}_{1}^{(k+1)}\!,\boldsymbol{x}_{2}^{(k+1)})\!\leq\!f(\boldsymbol{x}_{1}^{(k+1)}\!,\boldsymbol{x}_{2}^{(k)})$ and $\nabla_{2}f(\boldsymbol{x}_{1}^{(k+1)}\!,\boldsymbol{x}_{2}^{(k+1)})\!=\!0,\forall k$. The authors also proved the convergence when $f$ in (54) is minimized subject to a convex constraint set (see Corollary 1 and Section 4 in [41]). Equipped with this result from [40, 41], we return to our own problem. Let $f\!=\!-\sum_{l=1}^{L}{\cal F}_{l},\ \boldsymbol{x}_{1}\!=\!\{{\cal V}_{l}\}_{l=1}^{L},\ \boldsymbol{x}_{2}\!=\!\{P_{l},\boldsymbol{w}_{l}\}_{l=1}^{L}$ in problem (P2). When solving problem (P2) using the block-coordinate method with two blocks, we note that (SP2-2) has a globally optimal solution and thus $\boldsymbol{x}_{1}^{(k+1)}\!\!=\!\underset{\boldsymbol{x}_{1}}{\text{argmin}}\ f(\boldsymbol{x}_{1},\boldsymbol{x}_{2}^{(k)})$ is completely known. Also, our proposed solution for (SP2-1) satisfies the condition $\nabla_{2}\ f(\boldsymbol{x}_{1}^{(k+1)},\boldsymbol{x}_{2}^{(k+1)})\!=\!0$ (because it is the solution of KKT conditions for (SP2-1)). Furthermore, our extensive simulations indicate that the condition $f(\boldsymbol{x}_{1}^{(k+1)},\boldsymbol{x}_{2}^{(k+1)})\leq f(\boldsymbol{x}_{1}^{(k+1)},\boldsymbol{x}_{2}^{(k)})$ is always satisfied $\forall k$. Hence, we conclude that the output of the block-coordinate ascent method between (SP2-1) and (SP2-2) converges to a stationary point.

Regarding the convergence speed of the block-coordinate descent method, few works have obtained a convergence rate under special conditions on $f$ in (54). However, for the general case of non-convex $f$, even under convex constraints, no convergence rate is established in the literature. Our extensive simulations indicate that the average number of iterations needed for Algorithm 1 to converge to an $\epsilon$-accurate solution for $\{P_{l}^{opt},\boldsymbol{w}_{l}^{opt}\}_{l=1}^{L}$ is 30.

$\bullet$ Convergence of Algorithm 2: In this algorithm, we employ Golden section method to obtain $P_{trn}^{opt}$, where in each iteration we apply Algorithm 1 to find $\{{P}_{l},{\boldsymbol{w}}_{l}\}_{l=1}^{L}$, only for the purpose of successively narrowing the search interval of Golden section method. Consider solving the following non-convex minimization problem under convex constraints:

[TABLE]

where Golden section method is used to obtain $x^{opt}$. If $\boldsymbol{y}^{(k)}=\underset{\boldsymbol{y}}{\text{argmin}}\ f(x^{(k)},\boldsymbol{y})$ given $x^{(k)},\ k=0,1,...$ is known instantly, Golden section method converges linearly, and the rate of convergence is approximately 0.62 [38, p. 217]. Equipped with this result from [38], we return to our own problem. Let $f\!=\!-\sum_{l=1}^{L}{\cal G}_{l},\ x=P_{trn},\ \boldsymbol{y}\!=\!\{P_{l},\boldsymbol{w}_{l}\}_{l=1}^{L}$ in problem (P${}_{B^{\prime}}$). We obtain $x^{opt}$ using Golden section method, where in each iteration Algorithm 1 is applied to obtain $\boldsymbol{y}^{opt}$. Note that in Section IV-B, we proved that $f$ is strictly convex w.r.t. $x$ and thus, the convergence of Algorithm 2 to $x^{opt}$ is guaranteed. Since Algorithm 1 is an iterative algorithm with an unknown convergence rate, the exact convergence rate of Algorithm 2 is unknown. We only know that convergence rate of Algorithm 2 is less than 0.62. Our extensive simulations indicate that the average number of iterations needed for Algorithm 2 to converge to an $\epsilon$-accurate solution for $x^{opt}$ is 15.

VIII Numerical and Simulation Results

In this section, we corroborate our analytical results with numerical simulations, compare the effectiveness of different proposed power optimization schemes in acheiveing an MSE distorion-power tradeoff which is close to the Bayesian CRB, and investigate how the allocated power across clusters vary as signal-to-noise ratio (SNR) changes.

VIII-A Comparing $D$ and its Lower Bounds

Suppose $\theta$ is zero-mean with $\sigma^{2}_{\theta}\!=\!1$ and $L\!=\!10$ clusters. To enforce the heterogeneity in the network, we randomly choose $\sigma_{h_{l}},\sigma_{v_{l}},\sigma_{n_{l,k}},\sigma_{q_{l,k}}\!\in\!(0,1)$, and $K_{l}\!\in\!\{1,2,...,10\},l=1,...,L,k=1,...,K_{l}$. To capture the effect of randomness in flat fading channel coefficients and communication noise, the numerical results are computed based on $10^{6}$ Monte-Carlo trials, where in each trial, one realization of $|h_{l}|,\nu_{l},\forall l$ are generated. We also assume $\epsilon\!=\!10^{-3}$. In Section III-B we derived three lower bounds on $D$, of which we optimized $D_{1},D_{2}$ in problems (P3), (P4), respectively. Fig. 3 plots optimized $D$, optimized $D_{1}$, optimized $D_{2}$ versus $P_{tot}$. Note that $D_{3}\!=\!0.0043$ is constant. Clearly, $D_{3}\!<\!D_{2}\!<\!D_{1}\!<\!D\!<\!\sigma^{2}_{\theta}$. Also, $D_{2},D_{1},D$ decrease as $P_{tot}$ increases.

VIII-B Comparing Different Power Allocation Schemes

We compare the effectiveness of power optimization schemes, obtained from solving (P1) and its special cases (P1-SC1), (P1-SC2), (P1-SC3), in decreasing the MSE of the LMMSE estimator. We also compare the optimized MSE with the Bayesian CRB $G^{-1}$ derived in Section III-C. Let $D_{t},D_{c},D_{d}$ denote the MSE corresponding to the optimal solutions of (P1-SC1), (P1-SC2), (P1-SC3), respectively. We know $D_{3}\!<\!G^{-1}\!<\!D\!<\!D_{t},D_{c},D_{d}\!<\!\sigma_{\theta}^{2}$. To quantify the efficacy of different power allocation (w.r.t three sets of optimization variables $P_{trn}$, $P_{l}$’s, ${\mathcal{P}}_{l}$’s) in closing the MSE performance gap $\sigma_{\theta}^{2}-G^{-1}$, we define three factors as the following:

[TABLE]

where $0\!\leq g_{t},g_{c},g_{d}\leq\!1$. A larger factor $g$ means that the particular power allocation is more effective in reducing the MSE performance gap (closing the MSE performance gap). Fig. 4 and Fig. 5 plot $g_{t},g_{c},g_{d}$ versus $P_{tot}$ for two sets of noise variances (in Fig.  5 $\sigma_{h_{l}},\sigma_{q_{l,k}}$ are chosen from a smaller interval $(0,0.5)$). For $g_{t}$ we plot three curves corresponding to $P_{trn}\!=\!5\%,25\%,60\%P_{tot}$. Fig. 4 shows $g_{c}\!>\!g_{t}(P_{trn}\!=\!5\%P_{tot})\!>\!g_{t}(P_{trn}\!=\!60\%P_{tot})\!>\!g_{d}\!>\!g_{t}(P_{trn}\!=\!25\%P_{tot})$. Whereas Fig. 5 shows $g_{t}(P_{trn}\!=\!5\%P_{tot})\!>\!g_{c}\!>\!g_{t}(P_{trn}\!=\!25\%P_{tot})\!>\!g_{d}\!>\!g_{t}(P_{trn}\!=\!60\%P_{tot})$. Evidently, a more accurate channel estimation does not necessarily lead into a smaller $D_{t}$. Two takeaway messages are: (1) $g_{t},g_{c},g_{d}\!>\!0$, i.e., the solution obtained from solving (P1) always leads into an MSE improvement, (2) the actual values of $g_{t},g_{c},g_{d}$ depend on the system parameters and $P_{tot}$. Note that in Fig. 4 at $P_{tot}\!=\!0\text{dB}$, $g_{t}\!=\!0.15$ (for $P_{trn}\!=\!25\%P_{tot}$), $g_{d}\!=\!0.17,g_{c}\!=\!0.48$, meaning that power allocation among CHs for training and ${\cal P}_{l}$, and among clusters for obtaining $P_{l}$ reduce the MSE performance gap to $15\%,17\%,48\%$, respectively. Combining the information given by $g_{t},g_{c},g_{d},G$ with the computational complexity analysis in Section VI provides the system designer with quantitative complexity-versus-MSE improvement tradeoffs offered by different power optimization schemes.

VIII-C Behavior of Power Allocation Across Clusters

We study the effect of heterogeneous clusters on the behavior of our proposed power allocation scheme to solve (P1) as $P_{tot}$ increases. Consider a network consisting $L\!=\!3$ clusters with $K_{l}\!=\!6,\sigma_{n_{l,k}}\!=\!\sigma_{n_{l}},\sigma_{q_{l,k}}\!=\!\sigma_{q_{l}},\forall l,k$. We define ${\gamma}_{l}^{o}\!=\!\frac{\sigma^{2}_{\theta}}{\sigma_{n_{l}}^{2}}$ as observation SNR of sensors within cluster $l$, ${\gamma}_{l}^{c}\!=\!\frac{1}{\sigma_{q_{l}}^{2}}$ as channel-to-noise ratio (CNR) corresponding to sensors-CHl links, and ${\gamma}_{l}^{d}\!=\!\frac{\sigma^{2}_{h_{l}}}{\sigma_{v_{l}}^{2}}$ as CNR corresponding to CHl-FC link. Let $\psi_{l}\ (\text{dB})\!=\!10\text{log}_{10}(\psi_{l}),P_{l}\ (\text{dB})\!=\!10\text{log}_{10}(P_{l}),{\cal P}_{l}\ (\text{dB})\!=\!10\text{log}_{10}({\cal P}_{l}),{\cal V}_{l}\ (\text{dB})\!=\!10\text{log}_{10}({\cal V}_{l})$, where ${\cal V}_{l}\!=\!P_{l}\!+\!{\cal P}_{l}$ represents the allocated power to cluster $l$, excluding its training power $\psi_{l}$. In the following we consider three scenarios: (i) when observation SNR ${\gamma}_{l}^{o}$ and CNR ${\gamma}_{l}^{c}$ are equal and CNR ${\gamma}_{l}^{d}$ are different across clusters, (ii) when observation SNR ${\gamma}_{l}^{o}$ and CNR ${\gamma}_{l}^{d}$ are equal and CNR ${\gamma}_{l}^{c}$ are different across clusters, (iii) when CNRs ${\gamma}_{l}^{c}$ and ${\gamma}_{l}^{d}$ are equal and observation SNR ${\gamma}_{l}^{o}$ are different across clusters.

Figs. 6(a), 6(b), 6(c), 6(d), respectively, depict $\psi_{l}\!$ (dB),${\cal V}_{l}\!$ (dB),$P_{l}\!$ (dB),${\cal P}_{l}\!\ (\text{dB}),\forall l$, versus $P_{tot}$ for ${\gamma}_{l}^{o}\!=\!5\!\ \text{dB},{\gamma}_{l}^{c}\!=\!5\!\ \text{dB},\forall l$ and ${\gamma}_{1}^{d}\!=\!14\!\ \text{dB},{\gamma}_{2}^{d}\!=\!8\!\ \text{dB},{\gamma}_{3}^{d}\!=\!2\!\ \text{dB}$. Regarding Fig. 6 we make the following observations: 1) all powers increase as $P_{tot}$ increases, 2) when $P_{tot}$ is small, only cluster 1 is active, and as $P_{tot}$ increases, clusters 2 and 3 become active in a sequential order, 3) in all regions of $P_{tot}$, a cluster with a larger ${\gamma}_{l}^{d}$ is allotted a larger $\psi_{l}$ (water filling), 4) in low-region to moderate-region of $P_{tot}$, a cluster with a larger ${\gamma}_{l}^{d}$ is allocated a larger ${\cal V}_{l}$ (water filling), and in high-region of $P_{tot}$, ${\cal V}_{l}$ of all clusters converge (uniform power allocation), 5) in all regions of $P_{tot}$, a cluster with a larger ${\gamma}_{l}^{d}$ is assigned a larger $P_{l}$ (water filling), 6) in low-region of $P_{tot}$, a cluster with a larger ${\gamma}_{l}^{d}$ is allocated a larger ${\cal P}_{l}$ (water filling), and in high-region of $P_{tot}$, a cluster with a larger ${\gamma}_{l}^{d}$ is allotted a smaller ${\cal P}_{l}$ (inverse of water filling). The behavior of $P_{l}$ and ${\cal P}_{l}$ in high-region of $P_{tot}$ can be explained by examining the behavior of ${\cal V}_{l}$. Note that, although CNRs ${\gamma}_{1}^{d},{\gamma}_{2}^{d},{\gamma}_{3}^{d}$ are different, the differences are compensated as $P_{tot}$ increases and ${\cal V}_{l}$ of all clusters converge. This fact implies the behaviors of $P_{l}$ and ${\cal P}_{l}$ in high-region of $P_{tot}$ are opposite, i.e., water filling and inverse of water filling power allocation for $P_{l}$ and ${\cal P}_{l}$, respectively.

Figs. 7(a), 7(b), 7(c), 7(d), respectively, depict $\psi_{l}$ (dB), ${\cal V}_{l}$ (dB), $P_{l}$ (dB), ${\cal P}_{l}\ (\text{dB}),\forall l$, versus $P_{tot}$ for ${\gamma}_{l}^{o}\!=\!5\ \text{dB},{\gamma}_{l}^{d}\!=\!5\ \text{dB},\forall l$ and ${\gamma}_{1}^{c}\!=\!14\ \text{dB},{\gamma}_{2}^{c}\!=\!8\ \text{dB},{\gamma}_{3}^{c}\!=\!2\ \text{dB}$. The following observations can be made for Fig. 7: comments 1) and 2) for Fig. 6 also hold for Fig. 7, 3) in all regions of $P_{tot}$, $\psi_{l}$ of all clusters are equal (uniform power allocation) since ${\gamma}_{l}^{d}$’s are equal, 4) behavior of ${\cal V}_{l}$ in Fig. 7(b) is the same as that of Fig. 6(b), 5) in low-region of $P_{tot}$, a cluster with a larger ${\gamma}_{l}^{c}$ is allocated a larger $P_{l}$ (water filling), and in high-region of $P_{tot}$, a cluster with a larger ${\gamma}_{l}^{c}$ is allocated a smaller $P_{l}$ (inverese of water filling), 6) in all regions of $P_{tot}$, a cluster with a larger ${\gamma}_{l}^{c}$ is allocated a larger ${\cal P}_{l}$ (water filling). Note that, although CNRs ${\gamma}_{1}^{c},{\gamma}_{2}^{c},{\gamma}_{3}^{c}$ are different, the differences are compensated as $P_{tot}$ increases and ${\cal V}_{l}$ of all clusters converge. This fact implies the behaviors of $P_{l}$ and ${\cal P}_{l}$ in high-region of $P_{tot}$ are opposite, i.e., inverse of water filling and water filling power allocation for $P_{l}$ and ${\cal P}_{l}$, respectively.

Figs. 8(a), 8(b), 8(c), 8(d), respectively, depict $\psi_{l}$ (dB), ${\cal V}_{l}$ (dB), $P_{l}$ (dB), ${\cal P}_{l}\ (\text{dB}),\forall l$, versus $P_{tot}$ for ${\gamma}_{l}^{c}\!=\!5\ \text{dB},{\gamma}_{l}^{d}\!=\!5\ \text{dB},\forall l$ and ${\gamma}_{1}^{o}\!=\!14\ \text{dB},{\gamma}_{2}^{o}\!=\!8\ \text{dB},{\gamma}_{3}^{o}\!=\!2\ \text{dB}$. The following observations can be made for Fig. 8: comments 1) and 2) for Figs. 6 and 7 also hold for Fig. 8, 3) in all regions of $P_{tot}$, $\psi_{l}$ of all clusters are equal (uniform power allocation) since ${\gamma}_{l}^{d}$’s are equal, 4) in all regions of $P_{tot}$ a cluster with a larger ${\gamma}_{l}^{o}$ is allocated a larger ${\cal V}_{l}$, a larger $P_{l}$, and a larger ${\cal P}_{l}$ (water filling). The behaviors of ${\cal V}_{l}$, $P_{l}$, ${\cal P}_{l}$ in high-region of $P_{tot}$ are different from the two previous scenarios (CNRs across clusters were different), in which ${\cal V}_{l}$ of all clusters converge as $P_{tot}$ increases. Here the difference in observation SNR across clusters cannot be compensated as $P_{tot}$ increases. Hence, ${\cal V}_{l}$ of clusters are different, such that a cluster with a larger (smaller) ${\gamma}_{l}^{o}$ is allocated a larger (smaller) ${\cal V}_{l}$.

One may wonder given our proposed power allocation scheme, how the powers allocated to a CH and a sensor would be different. To answer this question, we let $P_{{\mbox{\tiny CH}}_{l}}={\cal P}_{l}+\psi_{l}$ denote the sum of power that CHl consumes for transmitting its fused signal $y_{l}$ as well as its training symbol to the FC. Fig. 9 plots $P_{{\mbox{\tiny CH}}_{l}}$ and $P_{l,k}$ versus $P_{tot}$, using the same setup parameters of Fig. 8. We observe that for all clusters $P_{{\mbox{\tiny CH}}_{l}}\!>>\!P_{l,k},k=1,...,K_{l}$, i.e., the power allocated to each sensor is much smaller than the power allocated to each CH.

IX Conclusions

We studied distributed estimation of a random source in a hierarchical power constrained WSN, where CHs linearly fuse the received signals from sensors within their clusters, and transmit over orthogonal fading channels to the FC. Prior to data transmission, CHs send pilot symbols to the FC to enable channel estimation at the FC. We derived the MSE $D$ corresponding to the LMMSE estimator of the source at the FC, and established lower bounds on $D$, including the Bayesian CRB. We addressed constrained minimization of $D$ under the constraint on $P_{tot}$, where the optimization variables are: i) training power $P_{trn}$ and $\{\!\psi_{l}\!\}_{l=1}^{L}$, ii) sensor-CH data transmission powers $\{\!P_{l}\!\}_{l=1}^{L}$, iii) CH-FC data transmission powers $\{\!{\cal P}_{l}\!\}_{l=1}^{L}$. We demonstrated the superior performance of our proposed power allocation scheme, comparing with schemes obtained from solving special case problems where subsets of these variables are optimized. Our simulations revealed that

1. when CNR corresponding to CHl-FC link varies across clusters, $\psi_{l},P_{l}$ allocation follow water filling fashion in all regions of $P_{tot}$, ${\cal P}_{l}$ follows (inverse of) water filling fashion in (high-region) low-region of $P_{tot}$,
2. when CNR corresponding to sensors-CHl links varies across clusters, $P_{l}$ allocation follows (inverse of) water filling fashion in (high-region) low-region of $P_{tot}$, ${\cal P}_{l}$ allocation follows water filling fashion in all regions of $P_{tot}$,
3. when observation SNR varies across clusters, both $P_{l},{\cal P}_{l}$ allocation follow water filling fashion in all regions of $P_{tot}$, and they diverge from uniform power allocation scheme as $P_{tot}$ increases. Leveraging on this work, we discuss three future research directions as follows. First direction is considering a coherent multiple access channel model (instead of orthogonal channels) for intra-cluster communication, where sensors within a cluster transmit their amplified measurements to their CH simultaneously. Second direction is exploring distributed estimation of a random vector source with correlated components. Similar to our work, all sensors can make noisy measurements of a common vector source, or sensors of different clusters can make partial observations of the vector source. Third direction is studying a system where the FC is equipped with multiple antennas (MIMO system model).

-A Derivation of Bayesian CRB

Using the Bayes’ rule $f(\boldsymbol{z},\hat{\boldsymbol{h}},\theta)\!=\!f(\boldsymbol{z},\hat{\boldsymbol{h}}|\theta)f(\boldsymbol{\theta})$, we can decompose $G$ into two terms [33]:

[TABLE]

in which the outer expectations are taken over the pdf of $\theta$, denoted as $f(\theta)$. Note that $\mathbb{E}\{G_{1}(\theta)\}$ depends on $f(\theta)$ [35]. For instance, if $\theta$ is Gaussian with variance $\sigma^{2}_{\theta}$, we obtain $\mathbb{E}\{G_{1}(\theta)\}\!=\!\sigma^{-2}_{\theta}$. Since $\hat{\boldsymbol{h}}$ and $\theta$ are independent, the Bayes’ rule says $f(\boldsymbol{z},\hat{\boldsymbol{h}}|\theta)\!=\!f(\boldsymbol{z}|\hat{\boldsymbol{h}},\theta)f(\hat{\boldsymbol{h}})$, and we can rewrite G_{2}(\theta)\!=\!-\mathbb{E}\{\mathbb{E}\{\frac{\partial^{2}\ln f(\boldsymbol{z}|\hat{\boldsymbol{h}},\theta)}{\partial\theta^{2}}\big{|}\hat{\boldsymbol{h}}\}\}, where the outer and inner expectations are taken over the pdfs $f(\hat{\boldsymbol{h}})$ and $f(\boldsymbol{z}|\hat{\boldsymbol{h}},\theta)$, respectively. We note that $G_{2}(\theta)$ depends on the parameters of the observation model at the sensors as well as the physical layer parameters corresponding to sensors-CHs and CHs-FC links. One can show that $z_{l}$’s conditioned on $\hat{\boldsymbol{h}},\theta$ are independent, i.e., $f(\boldsymbol{z}|\hat{\boldsymbol{h}},\theta)=\prod_{l=1}^{L}f(z_{l}|\hat{h}_{l},\theta)$. Moreover, since channel estimation is performed independently for each cluster, we have $f(\hat{\boldsymbol{h}})\!=\!\prod_{l=1}^{L}f(\hat{h}_{l})$. Hence $G_{2}(\theta)$ becomes:

[TABLE]

Using the following two facts:

[TABLE]

[TABLE]

we find that $G_{2}(\theta)$ reduces to (27). Examining (27) we realize that we need to find two terms in order to fully characterize $G_{2}(\theta)$: the conditional pdf $f(z_{l}|\hat{h}_{l},\theta)$, and its first derivative with respect to $\theta$, $\partial f(z_{l}|\hat{h}_{l},\theta)/\partial\theta$. In the following, we derive these two terms. Using (15) we can write the received signal at the FC from CHl as:

[TABLE]

in which $u_{1_{l}},u_{2_{l}},v_{l}$ are mutually independent conditioned on $\hat{h}_{l},\theta$. Let $\bar{z}_{l}\!=\!u_{1_{l}}u_{2_{l}}$. Hence, $z_{l}\!=\!\bar{z}_{l}+v_{l}$. Next, we find the conditional pdf of $\bar{z}_{l}$, conditioned on $\hat{h}_{l},\theta$. Considering (5), we note that $h_{l}$, $\nu_{l}$ are zero-mean independent complex Gaussian, and hence from (14) we find that $\hat{h}_{l}$ is also a zero-mean complex Gaussian. Since $h_{l}\!=\!\hat{h}_{l}\!+\!\tilde{h}_{l}$, we have $\tilde{h}_{l}\!\sim\!\mathcal{CN}\left(0,\zeta_{l}^{2}\right)$. Also, $u_{1_{l}}\!\sim\!\mathcal{CN}\left(\hat{h}_{l},\zeta_{l}^{2}\right)$ and $u_{2_{l}}\!\sim\!\mathcal{N}\left(\bar{\mu}_{l},\bar{\sigma}_{l}^{2}\right)$ in (58), where $\bar{\mu}_{l}\!=\!\theta{\boldsymbol{w}_{l}}^{T}\sqrt{\boldsymbol{A}_{l}}\boldsymbol{1}_{l}$, $\bar{\sigma}_{l}^{2}\!=\!{\boldsymbol{w}_{l}}^{T}(\sqrt{\boldsymbol{A}_{l}}\boldsymbol{\Sigma}_{n_{l}}\sqrt{\boldsymbol{A}_{l}}+\boldsymbol{\Sigma}_{q_{l}})\boldsymbol{w}_{l}$. To find the conditional pdf of $\bar{z}_{l}$ we use the following lemma from [42].

Lemma 2**.**

If $X\!\sim\!\mathcal{CN}\left(\mu_{x}e^{j\phi_{x}},\sigma_{x}^{2}\right)$ and $Y\!\sim\!\mathcal{CN}\left(\mu_{y}e^{j\phi_{y}},\sigma_{y}^{2}\right)$ are independent complex Gaussian random variables, the pdf of $Z\!=\!XY$ (which is equal to the joint pdf of its real and imaginary parts) is:

[TABLE]

where $k_{x}\!=\!\mu_{x}/\sigma_{x},k_{y}\!=\!\mu_{y}/\sigma_{y},|Z|\!=\!\sqrt{z_{r}^{2}+z_{i}^{2}},\angle Z\!=\!\arctan(z_{i}/z_{r})$, and $K_{r}(x)$ is the modified Bessel function of the second kind with order $r$ and argument $x$. **

Therefore, we can write the conditional joint pdf $f(\bar{z}_{l_{r}},\bar{z}_{l_{i}}|\hat{h}_{l},\theta)$ using (59). Recall $v_{l}\!\sim\!\mathcal{CN}\left(0,2\sigma_{v_{l}}^{2}\right)$. Hence $f(v_{l})\!=\!f(v_{l_{r}},v_{l_{i}})\!=\!\frac{1}{(2\pi\sigma_{v_{l}}^{2})}\exp{(-\frac{v_{l_{r}}^{2}+v_{l_{i}}^{2}}{2\sigma_{v_{l}}^{2}})}$. Since $\bar{z}_{l}$ and $v_{l}$ are independent, the conditional joint pdf $f(z_{l_{r}},z_{l_{i}}|\hat{h}_{l},\theta)$ is computed as $f(z_{l_{r}},z_{l_{i}}|\hat{h}_{l},\theta)\!=\!f(\bar{z}_{l_{r}},\bar{z}_{l_{i}}|\hat{h}_{l},\theta)*f(v_{l_{r}},v_{l_{i}})$, in which $*$ is the operator for two-dimensional convolution. Substituting for $f(\bar{z}_{l_{r}},\bar{z}_{l_{i}}|\hat{h}_{l},\theta),f(v_{l_{r}},v_{l_{i}})$ from above and defining $b=|b|e^{j\angle b}$, after some mathematical manipulations, we reach $f(z_{l}|\hat{h}_{l},\theta)$ and $\frac{\partial f(\!z_{l}|\hat{h}_{l},\theta)}{\partial\theta}$ in (1) and (1), respectively, whose parameters are defined in (30).

Substituting (1) and (1) in (27), we compute $G_{2}(\theta)$.

-B Proof of Proposition 1: Finding ${\boldsymbol{w}}^{opt}_{l},{\cal F}^{opt}_{l}$ in terms of $P_{l}$

According to (39), the only non-zero eigenvalue of $\boldsymbol{\cal B}_{1}$ and its corresponding eigenvector are:

[TABLE]

Define $\delta_{l}\!=\!{\cal V}_{l}-P_{l},\ \xi_{l}\!=\!\frac{\sigma^{2}_{\theta}}{{|\hat{h}_{l}|}^{2}}(\frac{\sigma_{v_{l}}^{2}}{\delta_{l}}+\zeta_{l}^{2}),\ \boldsymbol{\Sigma}_{\mu_{l}}\!=\!|\boldsymbol{\mu}_{l}|{|\boldsymbol{\mu}_{l}|}^{T},\ \phi_{l}\!=\!{|\hat{h}_{l}|}^{2}+\zeta_{l}^{2}+\frac{\sigma_{v_{l}}^{2}}{\delta_{l}},\ \boldsymbol{\Sigma}_{P_{l}}\!=\!\boldsymbol{\Sigma}_{q_{l}}\!+\!P_{l}\boldsymbol{\Delta}_{l},\ \boldsymbol{\Sigma}_{\phi_{l}}\!=\!\phi_{l}\boldsymbol{\Sigma}_{P_{l}}$. By substituting $\boldsymbol{\Delta}_{l},\boldsymbol{\Pi}_{l}$ into $\boldsymbol{\Omega}_{l}$ and $\boldsymbol{B}_{l}$, and $\boldsymbol{\Omega}_{l}$ into $\boldsymbol{R}_{t_{l}}$, $\boldsymbol{\cal B}_{1}$ in (38) becomes $\boldsymbol{\cal B}_{1}\!=\!\boldsymbol{\Sigma}_{\phi_{l}}+\xi_{l}\boldsymbol{\Sigma}_{\mu_{l}}$. Using the Binomial inversion Lemma [32] we compute ${\boldsymbol{s}}^{opt}_{l}$ in (60):

[TABLE]

From (61), we obtain ${\boldsymbol{w}}^{opt}_{l}$:

[TABLE]

where $\tau_{l}$ is defined in Proposition 1. To obtain $(a)$ in (62), we use the fact that ${|\boldsymbol{\mu}_{l}|}^{T}\boldsymbol{\Sigma}_{\phi_{l}}^{-1}\boldsymbol{R}_{t_{l}}\boldsymbol{\Sigma}_{\phi_{l}}^{-1}|\boldsymbol{\mu}_{l}|\!=\!\frac{\epsilon_{l}}{\phi_{l}^{2}}(1+\frac{\sigma^{2}_{\theta}}{{|\hat{h}_{l}|}^{2}}\epsilon_{l})$, where $\epsilon_{l}\!=\!{|\boldsymbol{\mu}_{l}|}^{T}\boldsymbol{\Sigma}_{P_{l}}^{-1}|\boldsymbol{\mu}_{l}|$. To obtain $(b)$ in (62), we use ${\boldsymbol{R}_{t_{l}}}^{-1}\boldsymbol{\rho}_{l}\!=\!\frac{\boldsymbol{\Sigma}_{P_{l}}^{-1}\boldsymbol{\rho}_{l}}{1+\sigma^{2}_{\theta}P_{l}\boldsymbol{\rho}_{l}^{T}\boldsymbol{\Sigma}_{P_{l}}^{-1}\boldsymbol{\rho}_{l}}$, which is established using the Binomial inversion lemma. We have ${\cal F}^{opt}_{l}\!=\!{|\boldsymbol{\mu}_{l}|}^{T}{\boldsymbol{s}}^{opt}_{l}$. Substituting ${\boldsymbol{s}}^{opt}_{l}$ from (61) in (60) and using the fact that $1\!-\!\sigma^{2}_{\theta}P_{l}\tau_{l}\!=\!\frac{1}{1+\sigma^{2}_{\theta}P_{l}\boldsymbol{\rho}_{l}^{T}\boldsymbol{\Sigma}_{P_{l}}^{-1}\boldsymbol{\rho}_{l}}$ we reach:

[TABLE]

-C Solution of the Problem in (IV-A2)

Define $\delta_{l}\!=\!{\cal V}_{l}-P_{l}$ and let $T$ denote the objective function in (IV-A2). We have $\frac{\partial T}{\partial\delta_{l}}\!=\!\frac{{|\hat{h}_{l}|}^{2}\beta_{l}^{2}P_{l}\tau_{l}}{\sigma_{v_{l}}^{2}{(\beta_{l}+\delta_{l})}^{2}}\!>\!0$, implying that the solution to (IV-A2) must satisfy the equality constraint $\sum_{l=1}^{L}\delta_{l}+P_{l}\!=\!\sigma P_{tot}$. Also, $\frac{\partial^{2}T}{\partial\delta_{i}\partial\delta_{j}}\!=\!0,\forall i\!\neq\!j$, and $\frac{\partial^{2}T}{\partial\delta_{l}^{2}}\!=\!\frac{-2{|\hat{h}_{l}|}^{2}\beta_{l}^{2}P_{l}\tau_{l}}{\sigma_{v_{l}}^{2}(\beta_{l}+\delta_{l})^{3}}\!<\!0,\forall l$. Thus the Hessian of $T$ with respect to $\delta_{l}$’s is diagonal and negative definite, proving that $T$ is jointly concave over $\delta_{l}$’s. Since the constraint is linear in $\delta_{l}$, the problem in (IV-A2) is concave. The Lagrangian function $\cal L$ associated with (IV-A2) is:

[TABLE]

where $\lambda,\eta_{l}$’s are the Lagrange multipliers. The KKT optimality conditions are:

[TABLE]

The condition (64c) implies $\eta_{l}\!=\!0$ for active clusters with $\delta_{l}\!>\!0$. From (64a) we infer:

[TABLE]

in which $\left[x\right]^{+}\!=\!max\{x,0\}$. Having $\delta_{l}^{opt}$, we find ${\cal V}_{l}^{opt}\!=\!\delta_{l}^{opt}\!+\!P_{l}$ given in (50a). Substituting (65) in the active constraint condition $\sum_{l=1}^{L}\!\delta_{l}\!+\!P_{l}\!=\!\sigma P_{tot}$, the Lagrange multiplier $\lambda$ becomes equal to the expression given in (50b), in which $\cal A$ is the set of active clusters. To uniquely determine $\cal A$, we carry out the following procedure. Let $L_{\cal A}\!=\!|{\cal A}|$ where $L_{\cal A}\!\leq\!L$. Suppose the clusters are indexed in the descending order of $\frac{{|\hat{h}_{1}|}^{2}P_{1}\tau_{1}}{\sigma^{2}_{v_{1}}}\!\geq\!\frac{{|\hat{h}_{2}|}^{2}P_{2}\tau_{2}}{\sigma^{2}_{v_{2}}}\!\geq\!...\!\geq\!\frac{{|\hat{h}_{L}|}^{2}P_{L}\tau_{L}}{\sigma^{2}_{v_{L}}}$. Choosing an $L_{\cal A}$ value we find $\lambda$ and compute $\delta_{l}^{opt}\!=\!\beta_{l}(\frac{|\hat{h}_{l}|}{\sigma_{v_{l}}}\sqrt{\frac{P_{l}\tau_{l}}{\lambda}}\!-\!1),\forall l$. If $\delta_{l}^{opt}\!>\!0,\ l\!=\!1,...,L_{\cal A}$ and $\delta_{l}^{opt}\!\leq\!0,\ l\!=\!L_{\cal A}\!+\!1,...,L$, then we have identified the set of active clusters $\cal A$ with their corresponding $P_{l},l\!\in\!{\cal A}$. Otherwise, we repeat this process for another $L_{\cal A}$ value. It is proved that the solution always exists and is unique[43].

-D Proof of Concavity of sub-problem $\!(\!a\!)$ of (P1-SC3) over $P_{l}\!$’s

We rewrite the cost function of sub-problem $(a)$, denoted as ${\cal F}$, as:

[TABLE]

where $b_{l}\!=\!\frac{1}{{|\hat{h}_{l}|}^{2}}(\frac{\sigma_{v_{l}}^{2}}{{\mathcal{P}}}\!+\!\zeta_{l}^{2}),s_{l}\!=\!\frac{1+b_{l}}{\sigma^{2}_{\theta}b_{l}},m_{l}\!=\!{\boldsymbol{\rho}_{l}}^{T}\boldsymbol{\Sigma}_{P_{l}}^{-1}\boldsymbol{\rho}_{l},\boldsymbol{\Sigma}_{P_{l}}\!=\!\boldsymbol{\Sigma}_{q_{l}}\!+\!P_{l}\boldsymbol{\Delta}_{l}$. We have $b_{l},s_{l},m_{l}\!>\!0$ and $\boldsymbol{\Sigma}_{P_{l}}\!\succ\!\boldsymbol{0}$. Also, $\frac{\partial m_{l}}{\partial P_{l}}\!=\!-{\boldsymbol{\rho}_{l}}^{T}\boldsymbol{\Sigma}_{P_{l}}^{-1}\boldsymbol{\Delta}_{l}\boldsymbol{\Sigma}_{P_{l}}^{-1}\boldsymbol{\rho}_{l}\!<\!0,\frac{\partial^{2}m_{l}}{\partial P_{l}^{2}}\!=\!2{\boldsymbol{\rho}_{l}}^{T}\boldsymbol{\Sigma}_{P_{l}}^{-1}\boldsymbol{\Delta}_{l}\boldsymbol{\Sigma}_{P_{l}}^{-1}\boldsymbol{\Delta}_{l}\boldsymbol{\Sigma}_{P_{l}}^{-1}\boldsymbol{\rho}_{l}\!>\!0$. One can obtain $\frac{\partial{\cal F}_{l}}{\partial P_{l}}\!=\!\frac{s_{l}(m_{l}+P_{l}\frac{\partial m_{l}}{\partial P_{l}})}{b_{l}(s_{l}+P_{l}m_{l})^{2}}$ and prove that $m_{l}\!+\!P_{l}\frac{\partial m_{l}}{\partial P_{l}}\!>\!0$ which infers $\frac{\partial{\cal F}_{l}}{\partial P_{l}}\!>\!0$, i.e.,

[TABLE]

${\cal F}$ in (66)

is an increasing function of $P_{l}$, and thus, the solution of sub-problem $(a)$ of (P1-SC3) must satisfy the equality constraint $P_{trn}\!+\!\sum_{l=1}^{L}\{P_{l}\!+\!{\mathcal{P}}\}\!=\!P_{tot}$. Furthermore

[TABLE]

The denominator of the right-hand side is positive. The numerator of the right-hand side can be simplified as $\text{num}\!=\!I_{1}\!+\!I_{2}\!+\!I_{3}$, where

[TABLE]

One can prove that $I_{1}\!<\!0,I_{2}\!<\!0,I_{3}\!=\!0$. Hence, $\text{num}\!<\!0$ and $\frac{\partial^{2}{\cal F}_{l}}{\partial P_{l}^{2}}\!<\!0$. The following sequences of inequalities are easy to verify:

[TABLE]

where ($a$) comes by the fact that ${\boldsymbol{\rho}_{l}}^{T}\boldsymbol{\Sigma}_{P_{l}}^{-1}\boldsymbol{\Delta}_{l}\boldsymbol{\rho}_{l}$ is scalar. The Hessian of ${\cal F}$ with respect to $P_{l}$’s is diagonal and negative definite, which proves that ${\cal F}$ is jointly concave over $P_{l}$’s. Moreover, the constraint is linear in $P_{l}$, and therefore finding $P_{l}$’s in sub-problem $(a)$ of (P1-SC3) is jointly concave over $P_{l}$’s and has a unique solution.