Beamforming Design towards Mutual Information Maximization for ...
Beamforming Design towards Mutual Information Maximization for Centralized Wireless Sensor Network
arXiv:1412.3448v2 [cs.IT] 21 Jun 2015
†
Yang Liu† , Jing Li† , Xuanxuan Lu† and Chau Yuen‡ Electrical and Computer Engineering Department, Lehigh University, Bethlehem, PA 18015, USA ‡ Singapore University of Technology and Design, 20 Dover Drive, 138682, Singapore Email: {[email protected], [email protected], [email protected], [email protected]}
Abstract—This paper focuses on joint beamforming design towards maximizing mutual information (MI) in a coherent multiple access channel (MAC) wireless sensor network (WSN) with nodes equipped with multiple antennae. We consider the scenario where the same source signal is observed by all sensors, with each sensor having independent noisy observation and individual power constraint. All the sensors transmit their observations to a preassigned node, called fusion center (FC), to perform further processing and data fusion. To attack this nonconvex hard problem, we adopt the weighted minimum mean square error(WMMSE) method to transform the original problem by introducing intermediate variables and then consult to the block coordinate ascent (BCA) methodology to develop iterative algorithms for solutions. Specifically, we design a 3-block BCA algorithm with each of its three subproblems having closed form solutions or being efficiently solvable by standard numerical solvers. As will be shown, this proposed 3-BCA algorithm exhibits a satisfactory convergence property. Besides that, a multiple block BCA algorithm is also developed, which has closed form solution to each subproblem (possibably up to a simple bisection search). This multiple block BCA algorithm cherishes low complexity, does not depend on numerical solver and can even give out fully analytical update in special circumstance. Extensive numerical results are presented to test our algorithms.
I. I NTRODUCTION Recently the wireless sensor network (WSN) has attracted great attentions due to its wide applications in practice [1]– [11]. A typical wireless sensor network has multiple sensors which are spacially distributed and wirelessly connected. Sensors in the same neighborhood monitor the same physical event or measure some common environmental parameters and transmit their (usually contaminated) observations to a preassigned fusion center (FC) to perform further data processing and fusion. The goal of data transmission and fusion in wireless sensor networks can be achieved more effectively by leveraging multiple antenna and linear beamforming techniques. It is always an interesting and meaningful problem to collaboratively design the beamformers so that the wireless sensor network can reliably transmit and recover the observed signals. Magnitudes of studies have been performed on the beamforming design problem in wireless sensor networks and solutions are provided from various perspectives. For example, the papers [1]–[4] aim at designing effective beamformers for
signal compression. [1] and [2] consider the perfect channel case, i.e. there exists no noise or fading in transmission from the sensors to the FC. Although the perfect channel assumption in [1] and [2] excludes power constraints and greatly simplifies the problem, it is too restrictive for wireless settings. More practical models with noisy channels are considered in [3]– [11] [3] considers the problem of transmission and fusion of scalar source signal for noisy multiple access channels (MAC), where all sensors share one total transmission power. In practice total power constraint could still be stringent since the sensors are usually sufficiently distributed within a large area and therefore power sharing is hard to realize. [4] studies noisy and fading channels and separate power constraint for each sensor, under the assumptions that all channel matrices are square and nonsingular. [5] considers the special case of scalar source signal, where separate power constraints and noisy channels are assumed. Compared to the forementoined literature, the most generic model for centralized wireless sensor network is first introduced in [6]. The model proposed by [6] considers fading and noisy channels, separate power constraint for each sensor and both orthogonal and coherent MAC. Besides the above, no additional assumptions are imposed on the dimension/rank of beamformers or channel matrices, i.e. beamformers can be compressive, redundancy-added or rate-1 and channel matrices can be slim, flat or square (singular or nonsingular). Due to the difficulty of the problem, [6] provides solutions to several important special cases subsumed in the generic model for coherent MAC, including the scalar source signal case, the noiseless sensor-FC channel case and the no-intersymbolinterference (no-ISI) noisy channel case. Following the exact generic model in [6], [7] develops an iterative block coordinate descent (BCD) method that is applicable to any general case for coherent MAC. Recently various strong-convergenceguaranteed BCD-based algorithms have been proposed in [8], which can solve the most generic model in [6] for coherent MAC and subsume the algorithm by [7] as a specialized realization. All of the above mentioned papers [1]–[8] adopt mean square error (MSE) as performance metric. Besides the MSE criterion, signal to noise ratio (SNR) is another crucial and commonly used metric for scalar signal
recovery. For the coherent MAC wireless sensor networks proposed in [6], joint beamforming design towards maximizing SNR is reported in [10] and [11]. Recently joint beamforming design to maximize mutual information (MI) for orthogonal MAC is considered in [9]. It worth noting that the beamforming design problems in MIMO multi-sensor decision-fusion system are closely related with those in other multi-agent communication networks, e.g. MIMO multi-relay and multiuser communication systems. Plenty of exciting results exist in literature, see, for example, [12]–[14] and the reference therein. The contribution of this paper is follows: 1) In this paper we research the joint beamforming design in coherent MAC wireless sensor network towards MI maximization. Just as the MSE and SNR metric, MI is also a very meaningful design criterion, which is commonly adopted in communication theory to evaluate the average information transmission rate of a system. In wireless sensor network, MI represents the average information of the source signal which can be extracted at the fusion center from the sensors’ observations for each use of channel. Compared to the great deal of existing literature focusing on MSE, however, not many results have been reported on MI optimization in the WSN context due to its difficult nature. One recent inspiring paper [9] provides the beamforming solution to maximize MI in the orthogonal MAC wireless sensor network. As will be seen, the original MI optimization problem in coherent MAC wireless sensor networks is also a highly nonconvex hard problem and efficient solutions are meaningful and desirable. 2) Inspired by the seminal idea of weighted minimum mean square error (WMMSE) method in [15] and [16], we introduce a weight matrix and a virtual FC receiver as intermediate variables (the original MI maximization problem does not assume the presence of linear filter at FC, since, according to the data-processing inequality in [20], MI will never increase whatever processing procedure is performed at the receiver) and develop block coordinate ascent (BCA) algorithms to efficiently solve the original problem. Here we decompose the MI problem into three subproblems—one subproblem to update the virtual FC receiver, one subproblem to update the weight matrix and the third one to jointly optimize the entire beamformers of all sensors. The two former subproblems have closed form solutions and the third one can be proved to be a standard second order cone programming (SOCP) problem. The convergence analysis shows that the limit points of our solutions satisfy Karush-Kuhn-Tucker (KKT) conditions of the original MI maximization problem. 3) Besides the above 3 block BCA algorithm, we also come out a multiple block BCA algorithm, which has closed form solutions (possibly up to a simple bisection search) for each subproblem and is consequently highly efficient for implementation and not reliant on numerical solvers. Moreover we show that, in special circumstance, fully analytical update is even possible for the multiple BCA algorithm. Complexity of this algorithm is examined and extensive numerical results show that this multiple block BCA algorithm exhibits quite
good convergence performance. The rest of the paper is organized as follows: Section II introduces the system model of the coherent MAC wireless sensor network and formulates the joint beamforming problem towards maximizing mutual information. In section III we propose two BCA based algorithms to solve our original problem, with the convexity, closed form solutions and convergence being discussed in full details. Section IV provides numerical experiment results. Section V concludes the article. Notations: In the sequel, we use bold lowercase letters to denote complex vectors and bold capital letters to denote complex matrices. 0, Om×n , and Im are used to denote zero vectors, zero matrices of dimension m × n, and identity matrices of order m respectively. AT , A∗ , AH , and A† are used to denote the transpose, the conjugate, the conjugate transpose(Hermitian transpose), and the Moore-Penrose pseudoinverse respectively of an arbitrary complex matrix A. Tr{·} denotes the trace operation of a square matrix. | · | denotes the modulus of a complex scalar, and k · k2 denotes the l2 -norm of a complex vector. vec(·) means vectorization operation of a matrix, which is performed by packing the columns of a matrix into a long one column. ⊗ denotes the Kronecker product. Diag{A1 , · · · , An } denotes the block diagonal matrix with its i-th diagonal block being the square complex matrix Ai , n n represent the cones of positive and H++ i ∈ {1, · · · , n}. H+ semidefinite and postive definite matrices of dimension n respectively. Here 0 and ≻ 0 denote that a square complex n n matrix belongs to H+ and H++ respectively. Re{x} means taking the real part of a complex value x. II. S YSTEM M ODEL Here we consider the centralized wireless sensor network as illustrated in Fig.1.
Fig. 1: Multi-Sensor System Model This system has L sensors and one fusion center. We assume that all sensors and the FC are equipped with multi-antenna.
Denote the number of the antennae of the i-th sensor as Ni , i ∈ {1, · · · , L}, and that of the FC as M . The source signal s is a complex vector of dimension K, i.e. s ∈ CK×1 . Each sensor utilizes a linear beamformer(transmiter/precoder) Fi ∈ CNi ×K to transmit its observed data. In the beamforming problem to optimize MSE or SNR, a linear receiver(postcoder) is usually employed on the side of fusion center. In fact the presence of linear receiver at the FC leads to joint optimization of transmitters and receivers and can greatly improve the performance metric in terms of MSE or SNR. However linear receiver is not considered for MI maximization problem. According to information nonincreasing principle, any kinds of processing at the receiver will not increase the mutual information between the source and receiver. Thus, without loss of optimality, no filter is necessary at FC.
fusion center is presented as:
Here we adopt the assumption that the source signal follows zero mean circularly symmetric complex Gaussian distribution, i.e. s ∼ CN(0, Σs ) with Σs being positive definite. The meaning of the Gaussian signaling assumption has many folds as follows: generally MI lacks of analytical expression and the Gaussian signaling is one of the few exceptions having elegant closed form, which makes our problem analyzable. Moreover, for given source covariance and Gaussian vector channel, Gaussian source maximizes the mutual information [20], i.e. achieves the channel capacity. So in practice sensors can perform signaling transformation to approximate the transformed signals’ distribution to Gaussian distribution to improve transmission efficiency [18], [19]. At least, Gaussian source can provide an upper bound of the reduced uncertainty of the source at the fusion center. At the same time by central limit theory, the Gaussian signaling can serve as a good approximation for a large number of observations which follow independent and identical distribution.
It should be pointed out that the whiteness assumption of the Gaussian noise n0 at the receiver does not undermine generality of the model. Indeed if n0 ∼ CN(0, Σ0 ) has coloured 1 −1 e i , Σ− 2 Hi and covariance Σ0 , by redefining e r , Σ0 2 r, H 0 −1 e 0 , Σ0 2 n0 , the received signal can be equivalently written n as
The channel status {Hi }L i=1 are assumed to be known at the receiver, which can be achieved by standard channel estimation technique via pilots. We denote Hi ∈ CM×Ni as the channel coefficients from the i-th sensor to the fusion center. Due to interference from surroundings or thermal noise from the sensor device, the observed signals at the sensors are typically contaminated. We assume that the corruptions are additive zero mean circularly symmetric complex Gaussian noise, i.e. ni ∼ CN(0, Σi ), i ∈ {1, · · · , L} with Σi ∈ CK×K being covariance matrix. Since the sensors are spatially distributed, it is reasonable to assume that the noise ni at different sensors are mutually uncorrelated. Here we consider the coherent multiple access channels (MAC), which means the data from different sensors are superimposed at the fusion center. Here we assume that the transmissions is the network are timesynchronous, i.e. the FC receives data from different sensors in the same time slot, which can be realized via GPS system. The collected data at the fusion center is corrupted by additive Gaussian noise n0 . Without loss of generality, n0 is assumed to be white and zero-mean, i.e. n0 ∈ CM×1 ∼ CN(0, σ02 IM ).
Based on the system model above, the transmitted signal at the i-th sensor is Fi (s + ni ), and the received signal at the
L X
Hi Fi s + ni + n0
(1)
L X Hi Fi ni + n0 , Hi Fi s +
(2)
r=
i=1
=
L X i=1
i=1
|
{z n
}
where the compound noise vector n is still Gaussian, i.e. n ∼ CN(0, Σn ) with its covariance matrix Σn as Σn = σ02 IM +
L X
H Hi Fi Σi FH i Hi .
(3)
i=1
e r=
L X i=1
e i Fi s + ni + n e0, H
(4)
e 0 ∼ CN(0, IM ), which coincides with the model in (1). with n The mutual information between the source signal and the received signal at FC can be given in equation (5) on the top of next page [17]. In practice, each sensor has independent transmission power according to its own battery condition. The average transmitted 2
s+n power for the i-th sensor is E F = Tr Fi Σs + i i 2 , which must respect its power constraint Pi . Thus Σi FH i the beamforming problem of the multiple sensor system can be formulated as the following optimization problem: (6a) (P0) : max .MI {Fi }L i , H s.t.Tr Fi Σs + Σi Fi ≤ Pi , i ∈ {1, · · · , L}. (6b)
The above optimization problem is nonconvex, which can be easily seen by examining the convexity of the special case where {Fi }L i=1 are all scalars. Efficient solutions to (P0) are desirable. Since the above problem can hardly be solved in one shot, we propose iterative algorithms which fall in the framework of block coordinate descent/ascent (BCD/A) algorithms [22], also known as alternative minimization/maximization algorithm(AMA) [21] or Gauss-Seidel (GS) algorithms somewhere else [22] [23]. III. A LGORITHM D ESIGN In this section, we focus on solutions to the problem (P0). Note that directly utilizing BCA method to partition the beamformers into groups does not help to simplify our problem. Even if only one separate beaformer is considered, the objective is still hard. Inspired by the weighted mean
MI
Fi
L i=1
L L L −1 H X X X H Hi Fi Σi FH H σ02 I + Hi Fi Hi Fi Σs = log det IM + i i i=1
square error(WMMSE) method proposed by the seminal papers [15] and [16], we introduce auxiliary variables to convert the objective into a BCA-friendly form and then decompose the problem into solvable subproblems. Interestingly, although mutual information is independent of processing techniques at the receiver, our solution actually introduces a virtual linear filter at the fusion center to achieve our goal. Firstly we introduce two useful lemmas which pave the way for transforming the original hard problem (P0). Lemma 1 ( [15], [16]). For any positive definite matrix E ∈ n H++ , the following fact holds true − log det(E) = maxn. log det W − Tr{WE} + n (7)
Notice that W is positive definite, it can be cancelled and thus the equation (11) has been obtained. By substituting (11) into (9), (12) can be proved. Comment III.1. For the special case W = I, the result in lemma 2 is the well known Wiener filter. Here lemma 2 actually slightly generalizes this well known result. As we have shown above, when the mean square error is weighted by a matrix W, the Wiener filter maintains its optimality as long as the weighted parameter W is positive definite. Now by introducing the notation
W∈H++
e , H
with the optimal solution W⋆ given as W⋆ = E−1 .
(8) Lemma 2. Define a matrix function E G of variable G as H E G , I − GH H Σs I − GH H + GH Σn G, (9)
(5)
i=1
i=1
L X
Hi Fi ,
(16)
i=1
and the notations in equation (3), we can transform our objective function MI {Fi }L i=1 as the following:
e H Σ−1 e sH (17) = log det I + HΣ MI {Fi }L M n i=1 with Σs and Σn being positive definite matrices. Then for any −1 e H −1 e (18) positive definite matrix W, the following optimization problem = log det H Σn H+Σs Σs −1 −1 e H Σ−1 H+Σ e min .Tr WE(G) (10) = −log det H +log det Σs (19) n s G n −1 −1 can be solved by the optimal solution e H Σ−1 H+Σ e = max . log det W −Tr W H n s K W∈H++ −1 o HΣs . (11) G⋆ = HΣs HH + Σn (20) +K +log det Σs At the same time, E(G⋆ ) is given as n n He H He −1 . (12) = max . log det W −Tr W I−G H Σs I−G H E(G⋆ ) = HH Σ−1 H + Σ−1 n
s
Proof: The problem in (10) is a convex problem. To see this, notice that the objective function in (10) is a quadratic function of G with its quadratic terms being given as Tr WGH HΣs HH G + Tr WGH Σn G . (13)
By the identities Tr{AB} = Tr{BA} and Tr{ABCD} = vecT (DT ) CT ⊗ A vec(B), the first term of the above quadratic terms can be rewritten as h i Tr WGH HΣs HH G =vecH(G) W∗⊗ HΣs HH vec(G). (14)
W∈HK , ++
G
oo +K +log det Σs , +GH Σn G
(21)
where the last two steps follow lemma 1 and 2 respectively. Thus the optimization problem (P0) maximizing MI has been transformed into an equivalent problem (P1) in (22), which is shown on the top of next page. As a straightforward consequence of the above two lemmas, we have obtained the optimal solutions to the following two subproblems of (P1). ⋆ When {Fi }L i=1 and G are given, the optimal W is given as W⋆ = arg max . MI W {Fi }L ,G (23) i=1 K
H Notice i s H are both positive semi-definite, h that W and HΣ so W∗⊗ HΣs HH is positive semi-definite [24] and thus the first quadratic term is a convex function of G. Similarly the second quadratic term in (13) can also be proved to be W∈H++ convex function of G. Thus (10) is non-constrained convex −1 L L X X H H problem of G. By setting the derivative with respective to G H H Hi Fi +G Σn G . Hi Fi Σs I−G = I−G to zero [25], we obtain i=1 i=1 h i ∂Tr WE(G) H W= O. (15) G − HΣ = HΣ H +Σ s s n ⋆ When {Fi }L ∂G∗ i=1 and W are given, the optimal G is given
(P1) max . W∈HK , ++
{Fi }L i=1 ,G
h L L i X X H H H H +G Σ G +log det(Σs )+K, H F I−G Σ H F MI {Fi}L ,W,G = log det W −Tr W I−G n i i s i i i=1 (22a)
≤ Pi , s.t. Tr Fi Σs+Σn FH i
i ∈ {1, · · · , L}.
as
G⋆ = arg max .MI G {Fi }L , W i=1
(24)
G
=
L h X i=1
i=1
i=1
L L i−1 X X H Hi Fi Σs Hi Fi +Σn Hi Fi Σs , i=1
i=1
with Σn being given in equation (3). Now we focus on the subproblem of optimizing {Fi }L i=1 with W and G given. Towards this end, we have two options—we can either jointly optimize {Fi }L i=1 in one shot, or we can further consult to BCA methodology again to partition the entire variables {Fi }L i=1 into L blocks, {F1 }, · · · , {FL } and attack L smaller problems one by one in a cyclic manner. For both of these two options, solutions, hopefully in a closed form, are desirable and complexity are concerned. In the following, we discuss these two alternatives in details. A. Jointly Optimizing {Fi }L i=1
The subproblem of (P1) maximizing MI {Fi }L i=1 W, G with W and G given is rewritten as follows h L L X X H Hi Fi Hi Fi Σs I−GH (P2) min . Tr W I−GH {Fi }L i=1
i=1
i H +G Σn G ,
i=1
≤ Pi , i ∈ {1, · · · , L}. s.t. Tr Fi Σs+Σi FH i
(25a) (25b)
The following theorem identifies the convexity of (P2). Theorem 1. The problem (P2 ) is convex. Proof: To begin with, we first look at the function f X : Cm×n 7→ R given as follows: f X , Tr Σ1 XΣ2 XH (26)
with constant matrices Σ1 and Σ2 being positive semidefinite and having appropriate dimensions. By the identity Tr{ABCD} = vecT (DT ) CT ⊗ A vec(B), f X can be equivalently written as (27) f X = vecH X) Σ∗2 ⊗ Σ1 vec X .
Since Σ1 and Σ2 are positive semi-definite, [Σ∗1 ⊗ Σ2 ] is positive semi-definite [24]. Thus f X is actually a convex function with respect to X. PL For a further step, we replace X = i=1 Hi Fi . Since PL H F is an affine (linear actually) transformation of i i i=1
(22b)
variables {Fi }L i=1 , and affine operations preserve convexity by [26], the following function L L n X X H o (28) H F Σ H F = Tr Σ f {Fi }L i i 2 i i 1 i=1 i=1
i=1
is a convex function with respect to variables {Fi }L i=1 jointly. To identify the convexity of the objective in (25a), it suffices to prove the nonlinear terms of {Fi }L i=1 are convex, which are given as Tr
n
GWGH
L X i=1
+
L X i=1
L X H o Hi Fi Hi Fi Σs
Tr
i=1
n
o H H . (29) HH i GWG Hi Fi Σi Fi
Based on the discussion at the beginning of this proof, each of above terms is convex and thus the objective is convex. Similarly the convexity of each power constraint function can also be recognized. Thus the problem (P2) is convex. After identifying the convexity of problem (P2), we reformulate it into a standard quadratic constrained quadratic problem(QCQP) problem. To this end, we introduce the following notations (30a) fi , vec Fi ; g , vec G ; (30b) ∗ H H (30c) Aij , Σs ⊗ Hi GWG Hj ; ∗ (30d) Bi , WΣs ⊗ Hi ; ∗ H H (30e) Ci , Σi ⊗ Hi GWG Hi ; T T f , f1 , · · · , fLT ; (30f) L A , Aij i,j=1 ; (30g) (30h) B , B1 , · · · , BL ; (30i) C , Diag C1 , · · · , CL ; o n ∗ Di,Diag OK(Pi−1Nj), Σs+Σi ⊗INi,OK(PL Nj) ; (30j) j=i+1 j=1 (30k) c , Tr WΣs + σ02 Tr GWGH . Based on the above notations, problem (P2) can be equivalently written as the following QCQP problem, (31a) (P3) : min f H A+C f −2Re gH Bf +c, f
s.t. f H Di f ≤ Pi ,
i ∈ {1, · · · , L}.
(31b)
By theorem 1, (P3) is convex, thus (A + C) is positive 1 semidefinite, which implies that its square root (A+C) 2 exists. Therefore the above problem can be further rewritten in a standard SOCP form as follows: (P3SOCP ) : min . t,
(32a)
f ,t,s H
s.t. s − 2Re{g Bf } + c ≤ t;
(A+C) 12 f
≤ s+1 ; s−1
2
2 1 2
Pi +1
Di2 f , i ∈ {1, · · · , L}.
Pi −1 ≤
2 2
(32b) (32c)
Now the problem of updating one separate beamformer becomes (P3i ) min . fiH Aii+Ci fi − 2Re gH Bi −qH i fi , (33a) fi
s.t. fiH Ei fi ≤ Pi
with the definitions of qi and Ei as follows X ∗ qi , Aij fj ; Ei , Σs+Σi ⊗INi .
(33b)
(34)
j6=i
We introduce the following notations λi,1 (32d) 1 1 − − H .. 2 Ei 2 Aii+Ci Ei 2 = Ui Ui ; (35a) . The above problem can be solved by standard numerical tools λi,KNi like CVX [27]. 1 − 2 pi = UH BH (35b) The method discussed above is summarized in algorithm 1. i Ei i g − qi , Algorithm 1: 3-Block BCA Algorithm to solve (P0) 1
2 3
4
5
6
(0)
Initialization: randomly generate feasible {Fi }L i=1 ; obtain G(0) by (24); obtain W(0) by (23); repeat with G(j−1) and W(j−1) being fixed, solve (P3) in (j) (32), obtain {Fi }L i=1 ; (j) L with {Fi }i=1 and W(j−1) being fixed, obtain G(j) by (24); (j) (j) being fixed, obtain W(j) with {Fi }L i=1 and G by (23); until increase of MI is sufficiently small or predefined number of iterations is reached;
For the proposed 3-BCA algorithm, we have the following conclusion on its convergence Theorem 2. Assume that the covariance matrix Σs ≻ 0. Algorithm 1 generates increasing MI sequence. Its solution sequence has limit points, and each limit point of the solution sequence is a KKT point of the original problem (P0 ). Proof: Refer to appendix A. B. Cyclic (L+1)-BCA Algorithm Although the above proposed 3-block BCA algorithm guarantees a satisfactory convergence, the subproblem (P3) relies on standard numerical solvers, e.g. interior point method [27], to obtain solutions. Closed form solutions to (P3) is unknown. According to the complexity analysis performed in next subsection III-C, when the number of sensors and/or antenna number of each sensor grows, the problem (P3) can be very large size and consequently highly computation demanding. So effective algorithms with lower complexity are desirable. In this subsection, we consult to BCA methodology again to further partition the variables {Fi }L i=1 into L singleton sets: {F1 },· · · ,{FL }. This results in a cyclic (L+1)-BCA algorithm, where only one separate beamformer Fi is optimized at each time and different beamformers are updated in an round robin manner.
i with the eigenvalues {λi,j }KN j=1 arranged in an decreasing order, i.e. λi,1 ≥ · · · ≥ λi,KNi . We denote the k-th element of pi as pi,k and assume that ri = rank Aii+Ci . Then the solution to problem (P3i ) is given by the following theorem.
Theorem 3. Under the assumption that Σs ≻ 0 or Σi ≻ 0, i ∈ {1, · · · , L}, the optimal solution of problem (P3 i ) is given as follows: CASE(I)—if either of the following two conditions holds: i) ∃k ∈ {ri + 1, · · · , KNi } such that |pi,k | 6= 0; PKNi Pri |pi,k |2 > Pi . or ii) k=ri+1 |pi,k | = 0 and k=1 λ2i,k The optimal solution to (P3 i ) is given by −1 H (36) Bi g − qi , fi⋆ = Aii+Ci+µ⋆i Ei
with the positive value µ⋆i being the unique solution to the following equation: ri X |pi,k |2 f (µi ) = = Pi . (37) (λi,k + µi )2 k=1
P i Pri |pi,k |2 CASE(II)—otherwise,i.e. KN k=ri+1 |pi,k | = 0, k=1 λ2i,k ≤ Pi , The optimal solution to (P3 i ) is given by 1 − 1 † − 1 −1 − fi⋆ = Ei 2 Ei 2 Aii+Ci Ei 2 Ei 2 BH i g − qi . (38)
Proof: For limit of space, please refer to Theorem 3 in [8] for detailed proof. In CASE(I) of theorem 3, equation (37) generally has no closed form solution. Notice that f (µi ) is a one-dimension strictly decreasing function of µi . So the determination of µ⋆i can be efficiently performed by a bisection search. Thus a finite interval containing µ⋆i is necessary from which the bisection search can start. The following lemma provides us bounds for µ⋆i . Lemma 3. A lower bound lbdi and upper bound ubdi for positive µ⋆i in (36) in theorem 3 can be given as follows: i) For subcase i) of CASE(I) + kpi k2 kpi k2 lbdi = √ − λi,1 , ubdi = √ ; (39) Pi Pi
ii) For subcase ii) of CASE(I) + kpi k2 − λi,1 , lbdi = √ Pi
with unitary matrix Ui , ui,1 , ui,2 , · · · , ui,KNi having its i columns {ui,j }KN j=1 satisfying the following properties
kpi k2 ubdi = √ − λi,ri , HH H i g Pi , and uH ui,1 = i,j Hi g = 0, for j = 2, · · · , KNi . (44) kHH (40) i gk2 It can be readily checked that the parameter pi in theorem 3 is given as:
where [x]+ , max(x, 0).
Proof: For limit of space, please refer to Lemma 1 in [8] ei )kHH pi,1 = σs2 (1 − gH q i gk2 ; pi,j = 0, j = {2, · · · , KNi }. (45) for detailed proof. By theorem 3 and lemma 3, we have obtained a nearly At this time, the function f (µi ) in (37) reduces to an elegant closed form solution to the problem (P3i ). Here we claim the form 2
2 above solution nearly closed form since it involves a bisection
ei HH σs4 1 − gH q i g 2 search. (46) f (µi ) = 2 . It is worth noting that the fully closed form solution to g µi + (σs2 + σi2 )gH Hi HH i problem (P3i ) does exist in the special but important case Based on the above observations, it can be concluded that of scalar source signal, i.e. K = 1. As will be seen in the the subcase i) of CASE(I) in theorem 3 will never occur. The complexity analysis in subsection III-C, the increase of source two cases for positive and zero µ⋆i can be specified as follows: signal dimension can extensively enlarge the beamforming ⋆ problem size and therefore its complexity. So in practice, when CASE(I)— µi > 0 2 2 ei > (σs2 +σi2 )2 P¯i kHH This is equivalent to σs4 1−gH q the wireless sensor network has adequate bandwidth, it is prei gk2 ferred to transmit the sensed data component by component to and optimal solutions are determined by
1 decrease the processing and beamforming design complexity ⋆ H H 2 ¯− 2
ei Hi g 2 −(σs2 +σi2 ) HH (47a) i g 2, at the fusion center. Another tempting reason to do so is that µi = σs Pi 1−g q −1 ⋆ when the source signal is scalar, µi can be obtained in an f ⋆ = σ 2 (1−gH q H ei ) µ⋆i I+(σs2 +σi2 )HH HH i g. (47b) s i gg Hi explicit way and therefore the bisection search is not needed. i Moreover, at this time, the eigenvalue decomposition (recall CASE(II)— µ⋆i = 0 2 (35a)) will not be involved in solving (P3i ). This conclusion ei ≤ (σs2 + This holds if and only if σs4 1 − gH q 2 ⋆ reads as the following corollary. σi2 )2 P¯i kHH i gk2 and the optimal fi is given by Corollary 1. For scalar transmission (K = 1 ), fully analytic solution to subproblem (P3 i ) can be obtained without evoking bisection search or eigenvalue decomposition. Proof: In the special case where the signal source is scalar, the variables and parameters in the subproblem optimizing one separate beamformer are specialized as follows W → w; Fi → fi ; G → g; Σs → σs2 ; Σi → σi2 .
(41)
ei , j6=i Hj fj , ignoring the terms independent By defining q of fi and omitting the constant positive factor w in the objective, the problem (P3i ) is rewritten as follows H (P3i ) : min . (σs2 + σi2 )fiH HH i gg Hi fi fi H eH (42a) − 2σs2 Re (1 − q i g)g Hi fi P i s.t. kfi k2 ≤ 2 , P¯i . (42b) σs + σi2 Solving the problem (P3i ) just follows the outline of theorem 3. Here, the key point leading to a closed form solution H is the fact that the quadratic matrix HH i gg Hi has rank-1, i.e. ri = 1 in theorem 3. Thus we obtain
H
0 O(KNi −1)×(KNi −1)
ei )HH σs2 (1−gH q i g . 2 H 2 (σs +σi )g Hi HH i g
(48)
Thus we have seen that for scalar transmission case, fully closed form solution to (P3i ) can be obtained without bisection search or eigenvalue decomposition. The cyclic (L+1)-BCA algorithm is summarized in algorithm 2. C. Complexity
P
(σs2 + σi2 )HH ggH Hi 2i (σs + σi2 )gH Hi HH i g = Ui 0
fi⋆ =
(43) UH i ,
In this subsection, we discuss the complexity of the proposed algorithms. The two subproblems optimizing G and W have closed form solutions in (24) and (23), their complexities come from matrix inversion and are given as O K 3 . For the 3-block BCA algorithm, the SOCP problem (P3SOCP ) in (32) is solved by jointly optimizing all beamformers. The complexity of solving an SOCP is [28] kX kX SOC SOC 1 2 3 2 2 nSOC,i , (49) nSOC,i+ O kSOC mSOC+mSOC i=1
i=1
where kSOC is the number of second order cone constraints, mSOC is the dimension of optimization problem and nSOC,i denotes the dimension of the i-th second order cone constraint. problem in (32), PL kSOC = L + 1, mSOC = PL For the , n = K K N SOC,1 i=1 Ni +1 for the first second i=1 i order cone constraint in (32c) and nSOC,i+1 = KNi + 1 for
Algorithm 2: Cyclic (L + 1)-BCA Algorithm to Solve (P0) Optimizing MI 1
2 3 4
5 6 7
8 9
(0)
Initialization: randomly generate feasible {Fi }L i=1 ; obtain G(0) by (24); obtain W(0) by (23); repeat for i = 1; i Pi then else if k=rii +1 |pi,k |2 = 0, k=1 λ2 i,k
10 11
12 13 14 15 16 17 18 19
determine bounds lbd i and ubdi by (40); bisection search on lbdi , ubdi to determine µ⋆i satisfying (37); obtain Fi by (36); else obtain Fi by (38); end update G by (24) ; update W by (23) ; end until increase of MI is sufficiently small or predefined number of iterations is reached;
the i-th power constraint in (32d), i ∈ {1, · · · , L}. Substituting these into (49), the complexity of solving (P3) is √ parameters P 3 O LK 3 ( L , this is also the complexity for each N ) i=1 i loop of 3-block BCA algorithm. For the cyclic (L + 1)-BCA algorithm, the problem (P3i ) optimizing one separate sensor’s beamformer has its major complexity coming from eigenvalue decomposition, 3 3 which i ). Thus the complexity for each loop is P is O(K N L 3 3 O i=1 K Ni . Clearly by fully decomposing the original problem and researching the solution structure of the subproblems, the (L + 1)-block BCA algorithm effectively lowers the computation complexity.
IV. N UMERICAL R ESULTS In this section, numerical results are presented to verify the algorithms proposed in the previous section. In our following experiments, we test the case where the source signal and all observation noise are colored. Specifically, we set the covariance matrices of the source signal and observation noise as Σs = σs2 Σ0 ,
Σi = σi2 Σ0 ,
i ∈ {1, · · · , L},
(50)
where the K × K Toeplitz matrix Σ0 is defined as .. 2 . ρK−1 1 ρ ρ .. .. ρ . . 1 ρ .. 2 Σ0 = ρ2 . . ρ ρ 1 . . . . .. .. .. .. ρ .. 2 K−1 . ρ ρ 1 ρ
(51)
The parameter ρ in the above equation is used to adjust the correlation level between different components of the signal or noise. In our test, ρ is set as ρ = 0.5. Here we define the observation signal to noise ratio at the i-th sensor as SNRi , σs2 σ2 and the channel signal to noise ratio as SNR , σs2 . σi2 0 In figure 2 and 3 we test the performance of the 3-block BCA and cyclic (L + 1)-BCA algorithms for multiple dimension source signal. Here two cases are tested—heterogeneous network and homogeneous network, in figures 2 and 3 respectively. In the heterogeneous network, the transmission power, observation noise level and numbers of antennae of each sensor are different. We set up a wireless sensor network with three sensors, i.e. L = 3. The dimension of the source signal and the number of antenna of the FC are chosen as 3 and 4 respectively, i.e. K = 3 and M = 4. We randomly set the antenna number for each sensor as N1 = 3, N2 = 4, and N3 = 5 respectively, the transmission power constraint for each sensor as P1 = 2, P2 = 2 and P3 = 3 respectively and the observation signal to noise ratio for each sensor as SNR1 = 8dB, SNR2 = 9dB and SNR3 = 10dB respectively. Comparatively, in homogeneous sensor network each sensor has the same transmission power, observation noise level and number of antenna. In this test case we assume that K = 4, M = 4, each sensor has Ni = 5 antennae and transmission power Pi = 2, with observation noise level SNRi = 9dB. In our test, to take into account the impact of the channel parameters, for the above system set-up and any specific channel SNR we randomly generate 500 channel realizations. For each channel realization, two proposed algorithms are run, both of which start from one common random feasible solution. The progress of MI with respect to outer-loop iteration numbers are recorded. For one given iteration number, the average MI performance over all 500 channel realizations is presented in figure 2 and 3. For the implementation of 3-BCA algorithm, SDPT3 solver of CVX is chosen. The blue solid curves represent the average MI performance obtained by 3-block BCA algorithm with different numbers of iterations and the red dotted ones represent those obtained by cyclic (L+1) BCA algorithm. The black dotted curve represents the average MI obtained by random full-power-transmission solutions, which are actually the average MI performance for feasible solutions which make all power constraints active. From figures 2 and 3, we see that the optimized beamformers obtained by the proposed algorithms present significant MI improvement compared to nonoptimized beamformers. Usually 40 to 50 iteration loops are sufficient to make the two algorithms
3−Block BCA and Cyclic (L+1) BCA with Different Numbers of Iterations−Heterogeneous Network 10
3−Block BCA and Cyclic (L+1) BCA with Different Numbers of Iterations−Scalar Source 3
2.8 Mutal Information(Nats per use of channel)
Mutal Information(Nats per use of channel)
9
8 9.1 7
9 8.9 8.8
6
3−blk BCA, 3 Iters 3−blk BCA, 7 Iters 3−blk BCA, 20 Iters 3−blk BCA 40 Iters Cyc. (L+1)−BCA, 3 Iters Cyc. (L+1)−BCA, 7 Iters Cyc. (L+1)−BCA, 20 Iters Cyc. (L+1)−BCA, 40 Iters Random Initials
8.7 8.6 5 8.5 4 4
0
2
4.5
5
4
5.5
6
6
8
10
12
14
16
2.6
2.95
2.4
2.9 2.85
2.2
2.8 2
2.75 2.7
1.8 −6
18
0 −4
−2
0
2
SNR(σ2/σ2) s
0
s
0
Fig. 4: Scalar Source Signal Case: 3-Block BCA Algorithm and Cyclic (L + 1)-Block BCA Algorithm with Different Numbers of Iterations.
3−Block BCA and Cyclic (L+1) BCA with Different Numbers of Iterations−Homogeneous Network 15
3 Block BCA and Cyclic (L+1) BCA Algorithms Starting from Different Initial Points−Test Case 1 9.5
14
9
13
Mutual Information(Nats per use of channel)
Mutal Information(Nats per use of channel)
6
SNR(σ2/σ2)
Fig. 2: Heterogenous Test Case: 3-Block BCA Algorithm and Cyclic (L+1)-Block BCA Algorithm with Different Numbers of Iterations.
12 11
3−blk BCA, 3 Iters 3−blk BCA, 7 Iters 3−blk BCA, 30 Iters 3−blk BCA 50 Iters Cyc. (L+1)−BCA, 3 Iters Cyc. (L+1)−BCA, 7 Iters Cyc. (L+1)−BCA, 30 Iters Cyc. (L+1)−BCA, 50 Iters Random Initials
14.2 14
10
13.8
9
13.6 8 13.4 7
8.5
9.3 9.2
8
9.1 9
7.5
8.9 8.8
7
8.7 8.6
6.5
8.5 8.4
6
8.3
13.2
2
5.5
6
4
0
2
4
7 6
8 8
10 2
6
8
10
12
14
16
18
20
3−blk BCA, Starting from One Specific Initial Point Cyc. (L+1)−BCA, Starting from One Specific Initial Point
13 6
5
4
3−blk BCA, 3 Iters 3−blk BCA, 10 Iters 3−blk BCA, 30 Iters 3−blk BCA 50 Iters Cyc. (L+1)−BCA, 3 Iters Cyc. (L+1)−BCA, 10 Iters Cyc. (L+1)−BCA, 30 Iters Cyc. (L+1)−BCA, 50 Iters Random Initials 2 4 6 8 10 12
12
14
16
18
2
SNR(σs /σ0)
5
0
10
20
30
40 Iteration Index
50
60
70
80
Fig. 3: Homogeneous Test Case:3-Block BCA Algorithm and Cyclic (L+1)-Block BCA Algorithm with Different Numbers of Iterations
Fig. 5: Heterogenous Test Case: Optimizing Generalized MI by 3-Block BCA Algorithm and Cyclic (L + 1)-Block BCA Algorithm with Different Initial Points
converge and the two algorithms finally converge to almost identical MI performance. In figure 4 we test the special case of scalar source signal (K = 1), where (L+1)-block BCA algorithm has fully closed form solution, which is summarized in corollary 1. In this experiment, we have the system setup as follows M = 4, N1 = 3, N2 = 4, N3 = 5, P1 = 1, P2 = 2, P3 = 3 SNR1 = 7dB, SNR2 = 8dB and SNR3 = 9dB. Similar results as in the multiple dimension source signal case have been obtained. In figure 5 and 6, we check the impact of the random initials to the proposed algorithms. We use the same system setup as those in figure 2 and 3 respectively. Here the channel param-
eters are randomly chosen and fixed. 10 feasible solutions, each of which makes all the power constraints active(satisfied with equality) are randomly generated. For each random initial point, we invoke the 3-block BCA and cyclic (L + 1)block BCA algorithms to optimize the beamformers. The MI progress of the proposed algorithms with different initials are illustrated in 5 and 6. It can be seen that both the 3-block BCA and cyclic (L+1)-block BCA are rather insensitive to selection of initial points. The two different algorithms with different initial points finally converge to almost identical value. Last we test the complexity of the proposed algorithms. In Table I the average MATLAB running time for each out-loop is
3 Block BCA and Cyclic (L+1) BCA Algorithms Starting from Different Initial Points−Homogeneous Network 14
Mutual Information(Nats per use of channel)
13 12 11
13
10 12.5 9 12
8 7
11.5
6 11 2
5 4
4
6
8
10
12
14
16
18
20
3−blk BCA, Starting from One Specific Initial Point Cyc. (L+1)−BCA, Starting from One Specific Initial Point 0
10
20
30
40 Iteration Index
50
60
70
80
Fig. 6: Homogeneous Test Case: Optimizing Generalized MI by 3-Block BCA Algorithm and Cyclic (L + 1)-Block BCA Algorithm with Different Initial Points TABLE I: MATLAB Running Time Per Outer-Loop (in Sec.) PP P L Algorithms L=5 L = 10 L = 20 Dim. PP P K=1 3-BCA M = Ni = 2 Cyc. (L+1) K=4 3-BCA M = Ni = 4 Cyc. (L+1) K=8 3-BCA M = Ni = 8 Cyc. (L+1) Notes: SDPT3 solver of CVX is
0.2748 0.0125 0.5079 0.0319 10.41 0.0488 chosen to
0.4018 0.0446 2.168 0.0999 90.51 0.1310 implement
0.7627 0.1792 12.93 0.3761 729.2 0.4747 3-BCA.
being fixed, the objective value obtained by solving the current sub-problem cannot be smaller than previous one. Thus the entire MI sequence keeps increasing. Under the positive definiteness assumption of Σs , Σs+ Σi ≻ 0. Thus ∀i ∈ {1, · · · , L} we have (52) kFi k2F λmin Σs+Σi ≤ Tr{Fi Σs+Σi FH i } ≤ Pi ,
where λmin (·) denotes the minimum eigenvalue of a Hermitian matrix. Since λmin Σs+Σi > 0, kFi k2F is finite for all i. Thus the variable {Fi }L i=1 is bounded. By Bolzano-Weierstrass the(kj ) L }i=1 orem, there exits a subsequence {kj }∞ j=1 such that {Fi converges. Since G and W are updated by continuous func- (kj ) L tions of {Fi }L }i=1 , W(kj ) , G(kj ) i=1 in (24) and (23), {Fi converges. Thus the existence of limit points in solution sequence has been proved. ¯ i }L , W, ¯ G ¯ is any limit point We assume that {F i=1 (k) (k) , G(k) . Then there exists a subseof {Fi }L i=1 , W j→∞ (k ) (kj ) quence {kj } such that {Fi j }L , G(kj ) −→ i=1 , W (k) L L ¯ ¯ ¯ {Fi }i=1 , W, G . Since {Fi }i=1 is bounded, by possibably restricting can assume that to a subsequence, we L (k +1) bi . F {Fi j }L i=1 converges to a limit i=1 (k +1)
Since for each j, {Fi j }L i=1 are feasible, i.e. (k +1) H (k +1) Tr{Fi j Σs+Σi Fi j )} ≤ Pi , i ∈ {1, · · · , L}. (53)
By taking j → ∞ in the above inequalities, we obtain H b i )} ≤ Pi , i ∈ {1, · · · , L}. b i Σs+Σi F (54) Tr{F L bi are feasible. So F i=1 L For any feasible Fi i=1 , we have presented. Here we consider the homogeneous sensor network. (k ) (k ) (k ) (k ) Different values of K, L and Ni are tested, which result in (k +1) L j j j ≤ MI {Fi j }L ,G j . (55) i=1 W different sizes of problem. SDPT3 solver of the CVX is used MI {Fi }i=1 W ,G to implement 3-BCA algorithm. As shown in the table, the Noticing that MI function is continuous and taking j → ∞ in cyclic (L+1) BCA algorithm is highly efficient, since each of the above, we obtain its update step can be performed in an almost analytical way. W, ¯ G ¯ ≤ MI {F b i }L W, ¯ G ¯ , MI {Fi }L (56) i=1 i=1 V. C ONCLUSION In this paper, we consider the linear beamforming design for any feasible {Fi }L i=1 . (k) L problem for a coherent MAC wireless sensor network to maxNotice the {Fi }i=1 generated by algorithm 1 are feasible, ¯ i }L are imize the mutual information. As we have seen, the original by continuity of power constraint functions, {F i=1 problem is nonconvex and difficult. To solve this problem, feasible. Thus we have we adopt the weighted minimum mean square error method ¯ i }L W, ¯ G ¯ ≤ MI {F b i }L W, ¯ G ¯ . MI {F (57) i=1 i=1 and block coordinate ascent method to decompose the original difficult problem into subproblems and examine the solution At the same time, since the MI sequence is increasing and to each subproblem, especially their closed form solution. ¯ i }L , W, ¯ G ¯ is a limit point of the solution sequence, {F i=1 The complexity and convergence of proposed algorithms are also discussed in details. Extensive numerical results are ¯ i }L W, ¯ G ¯ ≥ MI {F(k) }L W, ¯ G ¯ , MI {F (58) i=1 i=1 i presented to verify and compare the behaviors of the proposed for any integer k. Substitute k with kj in (58), take limit algorithms. ¯ i }L j → ∞ and combine it with (57), we have shown that {F i=1 A PPENDIX is actually an optimal solution to the problem (P2) with A. Proof of Theorem 2 ¯ and G. ¯ So {F ¯ i }L satisfy KKT conditions parameters W i=1 ¯ and G, ¯ which are listed in (59) Proof: Since for each sub-problem, we solve an optimiza- of (P2) with parameters W tion problem with respect to a subset of variables with others shown on the top of next page.
L X H ¯ ¯ ¯ ¯H ¯ ¯ ¯ ¯ ¯ Σs +HH F H −HH G W I− G i i i GWG Hi Fi Σi +λi Fi Σs+Σi = O, i i=1
n Ho ¯ ¯ i Σs+Σi F = 0, − P λi Tr F i i n Ho ¯ i ≤ Pi , ¯ i Σs+Σi F Tr F λi ≥ 0,
To simplify the following exposition, we introduce the following two notations: ¯ , H
L X
¯ i; Hi F
(60a)
i=1
¯ n , σ2 I + Σ 0
L X
¯ H HH . ¯ i Σi F Hi F i i
(60b)
i=1
According to the update step in algorithm L 1, the limit points ¯ and G ¯ have the relations with F ¯i W as follows. i=1 ¯ s, ¯ H+Σ ¯ n −1 HΣ ¯ = HΣ ¯ sH (61a) G H ¯ −1 ¯ −1 ¯ ¯ W =H Σ H+Σ . (61b) n
s
Utilizing (61) we can prove two identities in (62) and (63) respectively on the top of next page. Substituting equations (62) and (63) into (59a), we can rewrite the first order KKT conditions associated with only ¯ i }L as in equation (64) shown in the next page. {F i=1 To check the conditions of the original problem (P0), we need to determine the derivative of its Lagrangian function, or equivalently the derivative of MI with respect to {Fi }. By defining H,
L X
Hi Fi ,
(65)
i=1
the derivative of MI is calculated in (66) (shown in next page) with C1 (dFi ) and C2 (dFi ) being uninteresting terms involved dFi only and independent of d(F∗i ). By comparing the equations (64) with the derivative in (66b), it is easily to recognize that (64) is actually the first order KKT condition of problem (P0) optimizing MI. Together with equations (59b), (59c) and (59d), the KKT conditions of ¯ i }L . original problem have been proved to be satisfied by {F i=1 Thus the proof is complete. R EFERENCES [1] Y. Zhu, E. Song, J. Zhou, and Z. You, “Optimal dimensionality reduction of sensor data in multisensor estimation fusion,” IEEE Trans. Signal Process., vol. 53, no. 5, pp. 1631-1639, May 2005. [2] J. Fang and H. Li, “Optimal/near-optimal dimensionality reduction for distributed estimation in homogeneous and certain inhomogeneous scenarios,” IEEE Trans. Signal Process., vol. 58, no. 8, pp. 4339-4353, 2010.
i ∈ {1, · · · , L}; i ∈ {1, · · · , L}; i ∈ {1, · · · , L}; i ∈ {1, · · · , L}.
(59a) (59b) (59c) (59d)
[3] J. Fang and H. Li, “Power constrained distributed estimation with cluster-based sensor collaboration,” IEEE Trans. Wireless Commun., vol. 8, no. 7, pp. 3822-3832, July 2009. [4] I. D. Schizas, G. B. Giannakis, and Z.-Q. Luo, “Distributed estimation using reduced-dimensionality sensor observations,” IEEE Trans. Signal Process., vol. 55, no. 8, pp. 4284-4299, Aug. 2007. [5] Y. Liu, J. Li, X. Lu and C. Yuen, “Optimal linear precoding and postcoding for MIMO multi-sensor noisy observation problem,”, IEEE International Conference on Communications(ICC), Sydney , Jun. 2014. [6] J. Xiao, S. Cui, Z. Luo, and A. J. Goldsmith, “Linear coherent decentralized estimation,” IEEE Trans. Signal Process., vol. 56, no. 2, pp. 757-770, Feb. 2008. [7] A. S. Behbahani, A. M. Eltawil, H. Jafarkhani, “Linear decentralized estimation of correlated data for power-constrained wireless sensor networks,” IEEE Trans. Signal Process., vol. 60, no. 11, pp. 6003-6016, Nov. 2012. [8] Y. Liu, J. Li, X. Lu and C. Yuen, “Joint transceiver design for wireless sensor networks through block coordinate descent optimization,” available online: http://arxiv.org/abs/1409.7122. [9] J. Fang, H. Li, Z. Chen and Y. Gong, “Joint precoder design for distributed transmission of correlated sources in sensor networks,”, IEEE Trans. Wireless Commun., vol. 12, no. 6, pp. 2918-2929, June 2013. [10] Y. Liu, J. Li, and X. Lu, “Transceiver design for clustered wireless sensor networks — Towards SNR Maximization,” available online: http://arxiv.org/abs/1504.05311. [11] Y. Liu, J. Li, and X. Lu “Multi-terminal joint transceiver design for MIMO systems with contaminated source and individual power constraint,” IEEE Intl Symposium on Information Theory(ISIT), Honolulu, Jun. 2014. [12] S. Serbetli and A. Yener, “Transceiver optimization for multiuser MIMO systems,” IEEE Trans. Signal Process., vol. 52, pp. 214-226, Jan. 2004. [13] M. R. A. Khandaker and Y. Rong, “Joint transceiver optimization for multiuser MIMO relay communication systems”, IEEE Trans. Signal Processing, vol. 60, pp. 5977-5986, Nov. 2012. [14] Y. Rong, X. Tang, and Y. Hua, “A unified framework for optimizing linear non-regenerative multicarrier MIMO Relay Communication Systems”, IEEE Trans. Signal Process., vol. 57, pp. 4837-4851, Dec. 2009. [15] S. S. Christensen, R. Argawal, E. de Carvalho, and J. M. Cioffi, “Weighted sum-rate maximization using weighted MMSE for MIMO-BC beamforming design,” IEEE Trans. Wireless Commun., vol. 7, no. 12, pp. 1-7, Dec. 2008. [16] Q. Shi, M. Razaviyayn, Z. Luo, and C. He, “An iteratively weighted MMSE approcah to distributed sum-utility maximization for a MIMO iterfering broadcast channel,” IEEE Trans. Signal Process., vol. 59, no. 9, pp. 4331-4340, Dec. 2011. [17] I. E. Telatar, “Capacity of multi-antenna Gaussian channels,” Eur. Trans. Telecommun., vol. 10, no. 6, pp. 585-595, Nov.-Dec. 1999.
−1 ¯ H+Σ ¯ n −1 H ¯ ¯ sH ¯ H HΣ ¯ HΣ ¯ −1 H+Σ ¯ ¯ H+Σ ¯ n −1 HΣ ¯ s H ¯ sH ¯W ¯ I−G ¯ HH ¯ = HΣ I−Σs H G n s h i −1 ¯ H HΣ ¯ ¯ HΣ ¯ −1 H+Σ ¯ H+Σ ¯n − H ¯ sH ¯ HΣ ¯ −1 HΣ ¯ s H ¯ H+Σ ¯ n −1 HΣ ¯ ¯ sH ¯ H+Σ ¯ n −1 H+ ¯ sH ¯ ¯ H+Σ ¯ n −1 H ¯ sH H Σ HΣ = HΣ s n s n | {z } ¯ H+Σ ¯n ¯ sH = HΣ
−1
=O
¯ H
(62)
−1 ¯ H+Σ ¯ n −1 ¯ sH ¯ H HΣ ¯ HΣ ¯ −1 H+Σ ¯ ¯ s H Σs H HΣ n s i h ¯ H+Σ ¯ n −1 ¯ H HΣ ¯ sH ¯ H+HΣ ¯ sH ¯ HΣ ¯ −1 HΣ ¯ sH ¯ H+Σ ¯ n −1 HΣ ¯ sH ¯ sH = HΣ n i h ¯ H+Σ ¯ n −1 ¯ sH ¯ H+Σ ¯ n HΣ ¯ sH ¯ HΣ ¯ −1 HΣ ¯ sH ¯ H+Σ ¯ n −1 HΣ ¯ sH = HΣ n ¯ HΣ ¯ −1 ¯ sH ¯ H+Σ ¯ n −1 HΣ ¯ sH = HΣ n
¯ H+Σ ¯n ¯W ¯G ¯ H = HΣ ¯ sH G
HH i
h
−1
L L L L X X h i−1 X X H ¯ i Σs I ¯i H ¯ i Σs ¯H ¯ i Σi F Hi F Hi F Hi F σ02 I+ Hi F i Hi + i=1 L X
−
i=1
¯i Hi F
n
σ02 I +
L X
i=1
i=1
i=1
H
H −1
¯HH ¯ i Σi F Hi F i i
i=1
i
¯ i Σs+Σi = O, ¯ i Σi − λi F Hi F
i ∈ {1, · · · , L};
o d HΣs HH Σ−1 n io n −1 h H −1 H −1 HΣ d(H )Σ +HΣ H d(Σ ) + C1 (dFi ) = Tr I+HΣs HH Σ−1 s s n n n h n i o H −1 H HΣs I − HH Σ−1 = Tr HH + C2 (dFi ), n Hi Fi Σi d(Fi ) i Σn+HΣs H
d(MI) = Tr
⇒
(63)
I+HΣs HH Σ−1 n
(64)
−1
(66a)
L L L L h X X h H i−1X X ∂MI H H 2 H H F H F Σ H F + H F Σ F H σ I+ = H i i Σs I− i i s i i i i i i i 0 i ∂F∗i i=1 i=1 i=1 i=1 L X i=1
L i H 2 X H −1 Hi Fi Σi , i ∈ {1, · · · , L}. (66b) σ0 I+ Hi Fi Σi FH Hi Fi i Hi
[18] S.-Y. Chung, “On the construction of some capacityapproaching coding schemes,” Ph.D. dissertation, Dept. EECS, Massachusetts Institute of Technology, 2000. [19] Y. Hu, J. Garcia-Frias, and M. Lamarca, “Analog joint sourcechannel coding using non-linear curves and MMSE decoding,” IEEE Trans. Commun., vol. 59, no. 11, pp. 3016-3026, Nov. 2011. [20] T. M. Cover and J. A. Thomas Elements of Information Theory, 2nd ed., John Wiley & Sons, 2006. [21] D. Bertsekas and J. Tsitsiklis, Parallel and Distributed Computation, Englewood Cliffs, NJ: Prentice-Hall, 1989. [22] D. Bertsekas, Nonlinear Programming, 2nd ed., Belmont, MA:Athena Scientific, 1999. [23] L. Grippo, and M. Sciandrone, “On the convergence of the block nonlinear Gauss-Seidel method under convex constraints,” Operations Research Letters, vol. 26, pp. 127-136, 2000. [24] T. P. Minka. Old and New Matrix Algebra Useful for Statistics,
i=1
Notes, December 2000. [25] A. Hjorungnes, Complex-Valued Matrix Derivatives: With Applications in Signal Processing and Communications. Cambridge University Press, 2011. [26] S. Boyd, and L. Vandenberghe, Convex Optimization. New York: Cambridge University Press, 2004. [27] M. Grant and S. Boyd, “CVX: Matlab software for disciplined convex programming (web page and software),” [Online]. Available: http://cvxr.com/cvx Apr. 2010 [28] I. Polik and T. Terlaky, “Interior Point Methods for Nonlinear Optimization,” in Nonlinear Optimization, 1st edition, G. Di Pillo, F. Schoen, editors. Springer, 2010.
arXiv:1412.3448v2 [cs.IT] 21 Jun 2015
†
Yang Liu† , Jing Li† , Xuanxuan Lu† and Chau Yuen‡ Electrical and Computer Engineering Department, Lehigh University, Bethlehem, PA 18015, USA ‡ Singapore University of Technology and Design, 20 Dover Drive, 138682, Singapore Email: {[email protected], [email protected], [email protected], [email protected]}
Abstract—This paper focuses on joint beamforming design towards maximizing mutual information (MI) in a coherent multiple access channel (MAC) wireless sensor network (WSN) with nodes equipped with multiple antennae. We consider the scenario where the same source signal is observed by all sensors, with each sensor having independent noisy observation and individual power constraint. All the sensors transmit their observations to a preassigned node, called fusion center (FC), to perform further processing and data fusion. To attack this nonconvex hard problem, we adopt the weighted minimum mean square error(WMMSE) method to transform the original problem by introducing intermediate variables and then consult to the block coordinate ascent (BCA) methodology to develop iterative algorithms for solutions. Specifically, we design a 3-block BCA algorithm with each of its three subproblems having closed form solutions or being efficiently solvable by standard numerical solvers. As will be shown, this proposed 3-BCA algorithm exhibits a satisfactory convergence property. Besides that, a multiple block BCA algorithm is also developed, which has closed form solution to each subproblem (possibably up to a simple bisection search). This multiple block BCA algorithm cherishes low complexity, does not depend on numerical solver and can even give out fully analytical update in special circumstance. Extensive numerical results are presented to test our algorithms.
I. I NTRODUCTION Recently the wireless sensor network (WSN) has attracted great attentions due to its wide applications in practice [1]– [11]. A typical wireless sensor network has multiple sensors which are spacially distributed and wirelessly connected. Sensors in the same neighborhood monitor the same physical event or measure some common environmental parameters and transmit their (usually contaminated) observations to a preassigned fusion center (FC) to perform further data processing and fusion. The goal of data transmission and fusion in wireless sensor networks can be achieved more effectively by leveraging multiple antenna and linear beamforming techniques. It is always an interesting and meaningful problem to collaboratively design the beamformers so that the wireless sensor network can reliably transmit and recover the observed signals. Magnitudes of studies have been performed on the beamforming design problem in wireless sensor networks and solutions are provided from various perspectives. For example, the papers [1]–[4] aim at designing effective beamformers for
signal compression. [1] and [2] consider the perfect channel case, i.e. there exists no noise or fading in transmission from the sensors to the FC. Although the perfect channel assumption in [1] and [2] excludes power constraints and greatly simplifies the problem, it is too restrictive for wireless settings. More practical models with noisy channels are considered in [3]– [11] [3] considers the problem of transmission and fusion of scalar source signal for noisy multiple access channels (MAC), where all sensors share one total transmission power. In practice total power constraint could still be stringent since the sensors are usually sufficiently distributed within a large area and therefore power sharing is hard to realize. [4] studies noisy and fading channels and separate power constraint for each sensor, under the assumptions that all channel matrices are square and nonsingular. [5] considers the special case of scalar source signal, where separate power constraints and noisy channels are assumed. Compared to the forementoined literature, the most generic model for centralized wireless sensor network is first introduced in [6]. The model proposed by [6] considers fading and noisy channels, separate power constraint for each sensor and both orthogonal and coherent MAC. Besides the above, no additional assumptions are imposed on the dimension/rank of beamformers or channel matrices, i.e. beamformers can be compressive, redundancy-added or rate-1 and channel matrices can be slim, flat or square (singular or nonsingular). Due to the difficulty of the problem, [6] provides solutions to several important special cases subsumed in the generic model for coherent MAC, including the scalar source signal case, the noiseless sensor-FC channel case and the no-intersymbolinterference (no-ISI) noisy channel case. Following the exact generic model in [6], [7] develops an iterative block coordinate descent (BCD) method that is applicable to any general case for coherent MAC. Recently various strong-convergenceguaranteed BCD-based algorithms have been proposed in [8], which can solve the most generic model in [6] for coherent MAC and subsume the algorithm by [7] as a specialized realization. All of the above mentioned papers [1]–[8] adopt mean square error (MSE) as performance metric. Besides the MSE criterion, signal to noise ratio (SNR) is another crucial and commonly used metric for scalar signal
recovery. For the coherent MAC wireless sensor networks proposed in [6], joint beamforming design towards maximizing SNR is reported in [10] and [11]. Recently joint beamforming design to maximize mutual information (MI) for orthogonal MAC is considered in [9]. It worth noting that the beamforming design problems in MIMO multi-sensor decision-fusion system are closely related with those in other multi-agent communication networks, e.g. MIMO multi-relay and multiuser communication systems. Plenty of exciting results exist in literature, see, for example, [12]–[14] and the reference therein. The contribution of this paper is follows: 1) In this paper we research the joint beamforming design in coherent MAC wireless sensor network towards MI maximization. Just as the MSE and SNR metric, MI is also a very meaningful design criterion, which is commonly adopted in communication theory to evaluate the average information transmission rate of a system. In wireless sensor network, MI represents the average information of the source signal which can be extracted at the fusion center from the sensors’ observations for each use of channel. Compared to the great deal of existing literature focusing on MSE, however, not many results have been reported on MI optimization in the WSN context due to its difficult nature. One recent inspiring paper [9] provides the beamforming solution to maximize MI in the orthogonal MAC wireless sensor network. As will be seen, the original MI optimization problem in coherent MAC wireless sensor networks is also a highly nonconvex hard problem and efficient solutions are meaningful and desirable. 2) Inspired by the seminal idea of weighted minimum mean square error (WMMSE) method in [15] and [16], we introduce a weight matrix and a virtual FC receiver as intermediate variables (the original MI maximization problem does not assume the presence of linear filter at FC, since, according to the data-processing inequality in [20], MI will never increase whatever processing procedure is performed at the receiver) and develop block coordinate ascent (BCA) algorithms to efficiently solve the original problem. Here we decompose the MI problem into three subproblems—one subproblem to update the virtual FC receiver, one subproblem to update the weight matrix and the third one to jointly optimize the entire beamformers of all sensors. The two former subproblems have closed form solutions and the third one can be proved to be a standard second order cone programming (SOCP) problem. The convergence analysis shows that the limit points of our solutions satisfy Karush-Kuhn-Tucker (KKT) conditions of the original MI maximization problem. 3) Besides the above 3 block BCA algorithm, we also come out a multiple block BCA algorithm, which has closed form solutions (possibly up to a simple bisection search) for each subproblem and is consequently highly efficient for implementation and not reliant on numerical solvers. Moreover we show that, in special circumstance, fully analytical update is even possible for the multiple BCA algorithm. Complexity of this algorithm is examined and extensive numerical results show that this multiple block BCA algorithm exhibits quite
good convergence performance. The rest of the paper is organized as follows: Section II introduces the system model of the coherent MAC wireless sensor network and formulates the joint beamforming problem towards maximizing mutual information. In section III we propose two BCA based algorithms to solve our original problem, with the convexity, closed form solutions and convergence being discussed in full details. Section IV provides numerical experiment results. Section V concludes the article. Notations: In the sequel, we use bold lowercase letters to denote complex vectors and bold capital letters to denote complex matrices. 0, Om×n , and Im are used to denote zero vectors, zero matrices of dimension m × n, and identity matrices of order m respectively. AT , A∗ , AH , and A† are used to denote the transpose, the conjugate, the conjugate transpose(Hermitian transpose), and the Moore-Penrose pseudoinverse respectively of an arbitrary complex matrix A. Tr{·} denotes the trace operation of a square matrix. | · | denotes the modulus of a complex scalar, and k · k2 denotes the l2 -norm of a complex vector. vec(·) means vectorization operation of a matrix, which is performed by packing the columns of a matrix into a long one column. ⊗ denotes the Kronecker product. Diag{A1 , · · · , An } denotes the block diagonal matrix with its i-th diagonal block being the square complex matrix Ai , n n represent the cones of positive and H++ i ∈ {1, · · · , n}. H+ semidefinite and postive definite matrices of dimension n respectively. Here 0 and ≻ 0 denote that a square complex n n matrix belongs to H+ and H++ respectively. Re{x} means taking the real part of a complex value x. II. S YSTEM M ODEL Here we consider the centralized wireless sensor network as illustrated in Fig.1.
Fig. 1: Multi-Sensor System Model This system has L sensors and one fusion center. We assume that all sensors and the FC are equipped with multi-antenna.
Denote the number of the antennae of the i-th sensor as Ni , i ∈ {1, · · · , L}, and that of the FC as M . The source signal s is a complex vector of dimension K, i.e. s ∈ CK×1 . Each sensor utilizes a linear beamformer(transmiter/precoder) Fi ∈ CNi ×K to transmit its observed data. In the beamforming problem to optimize MSE or SNR, a linear receiver(postcoder) is usually employed on the side of fusion center. In fact the presence of linear receiver at the FC leads to joint optimization of transmitters and receivers and can greatly improve the performance metric in terms of MSE or SNR. However linear receiver is not considered for MI maximization problem. According to information nonincreasing principle, any kinds of processing at the receiver will not increase the mutual information between the source and receiver. Thus, without loss of optimality, no filter is necessary at FC.
fusion center is presented as:
Here we adopt the assumption that the source signal follows zero mean circularly symmetric complex Gaussian distribution, i.e. s ∼ CN(0, Σs ) with Σs being positive definite. The meaning of the Gaussian signaling assumption has many folds as follows: generally MI lacks of analytical expression and the Gaussian signaling is one of the few exceptions having elegant closed form, which makes our problem analyzable. Moreover, for given source covariance and Gaussian vector channel, Gaussian source maximizes the mutual information [20], i.e. achieves the channel capacity. So in practice sensors can perform signaling transformation to approximate the transformed signals’ distribution to Gaussian distribution to improve transmission efficiency [18], [19]. At least, Gaussian source can provide an upper bound of the reduced uncertainty of the source at the fusion center. At the same time by central limit theory, the Gaussian signaling can serve as a good approximation for a large number of observations which follow independent and identical distribution.
It should be pointed out that the whiteness assumption of the Gaussian noise n0 at the receiver does not undermine generality of the model. Indeed if n0 ∼ CN(0, Σ0 ) has coloured 1 −1 e i , Σ− 2 Hi and covariance Σ0 , by redefining e r , Σ0 2 r, H 0 −1 e 0 , Σ0 2 n0 , the received signal can be equivalently written n as
The channel status {Hi }L i=1 are assumed to be known at the receiver, which can be achieved by standard channel estimation technique via pilots. We denote Hi ∈ CM×Ni as the channel coefficients from the i-th sensor to the fusion center. Due to interference from surroundings or thermal noise from the sensor device, the observed signals at the sensors are typically contaminated. We assume that the corruptions are additive zero mean circularly symmetric complex Gaussian noise, i.e. ni ∼ CN(0, Σi ), i ∈ {1, · · · , L} with Σi ∈ CK×K being covariance matrix. Since the sensors are spatially distributed, it is reasonable to assume that the noise ni at different sensors are mutually uncorrelated. Here we consider the coherent multiple access channels (MAC), which means the data from different sensors are superimposed at the fusion center. Here we assume that the transmissions is the network are timesynchronous, i.e. the FC receives data from different sensors in the same time slot, which can be realized via GPS system. The collected data at the fusion center is corrupted by additive Gaussian noise n0 . Without loss of generality, n0 is assumed to be white and zero-mean, i.e. n0 ∈ CM×1 ∼ CN(0, σ02 IM ).
Based on the system model above, the transmitted signal at the i-th sensor is Fi (s + ni ), and the received signal at the
L X
Hi Fi s + ni + n0
(1)
L X Hi Fi ni + n0 , Hi Fi s +
(2)
r=
i=1
=
L X i=1
i=1
|
{z n
}
where the compound noise vector n is still Gaussian, i.e. n ∼ CN(0, Σn ) with its covariance matrix Σn as Σn = σ02 IM +
L X
H Hi Fi Σi FH i Hi .
(3)
i=1
e r=
L X i=1
e i Fi s + ni + n e0, H
(4)
e 0 ∼ CN(0, IM ), which coincides with the model in (1). with n The mutual information between the source signal and the received signal at FC can be given in equation (5) on the top of next page [17]. In practice, each sensor has independent transmission power according to its own battery condition. The average transmitted 2
s+n power for the i-th sensor is E F = Tr Fi Σs + i i 2 , which must respect its power constraint Pi . Thus Σi FH i the beamforming problem of the multiple sensor system can be formulated as the following optimization problem: (6a) (P0) : max .MI {Fi }L i , H s.t.Tr Fi Σs + Σi Fi ≤ Pi , i ∈ {1, · · · , L}. (6b)
The above optimization problem is nonconvex, which can be easily seen by examining the convexity of the special case where {Fi }L i=1 are all scalars. Efficient solutions to (P0) are desirable. Since the above problem can hardly be solved in one shot, we propose iterative algorithms which fall in the framework of block coordinate descent/ascent (BCD/A) algorithms [22], also known as alternative minimization/maximization algorithm(AMA) [21] or Gauss-Seidel (GS) algorithms somewhere else [22] [23]. III. A LGORITHM D ESIGN In this section, we focus on solutions to the problem (P0). Note that directly utilizing BCA method to partition the beamformers into groups does not help to simplify our problem. Even if only one separate beaformer is considered, the objective is still hard. Inspired by the weighted mean
MI
Fi
L i=1
L L L −1 H X X X H Hi Fi Σi FH H σ02 I + Hi Fi Hi Fi Σs = log det IM + i i i=1
square error(WMMSE) method proposed by the seminal papers [15] and [16], we introduce auxiliary variables to convert the objective into a BCA-friendly form and then decompose the problem into solvable subproblems. Interestingly, although mutual information is independent of processing techniques at the receiver, our solution actually introduces a virtual linear filter at the fusion center to achieve our goal. Firstly we introduce two useful lemmas which pave the way for transforming the original hard problem (P0). Lemma 1 ( [15], [16]). For any positive definite matrix E ∈ n H++ , the following fact holds true − log det(E) = maxn. log det W − Tr{WE} + n (7)
Notice that W is positive definite, it can be cancelled and thus the equation (11) has been obtained. By substituting (11) into (9), (12) can be proved. Comment III.1. For the special case W = I, the result in lemma 2 is the well known Wiener filter. Here lemma 2 actually slightly generalizes this well known result. As we have shown above, when the mean square error is weighted by a matrix W, the Wiener filter maintains its optimality as long as the weighted parameter W is positive definite. Now by introducing the notation
W∈H++
e , H
with the optimal solution W⋆ given as W⋆ = E−1 .
(8) Lemma 2. Define a matrix function E G of variable G as H E G , I − GH H Σs I − GH H + GH Σn G, (9)
(5)
i=1
i=1
L X
Hi Fi ,
(16)
i=1
and the notations in equation (3), we can transform our objective function MI {Fi }L i=1 as the following:
e H Σ−1 e sH (17) = log det I + HΣ MI {Fi }L M n i=1 with Σs and Σn being positive definite matrices. Then for any −1 e H −1 e (18) positive definite matrix W, the following optimization problem = log det H Σn H+Σs Σs −1 −1 e H Σ−1 H+Σ e min .Tr WE(G) (10) = −log det H +log det Σs (19) n s G n −1 −1 can be solved by the optimal solution e H Σ−1 H+Σ e = max . log det W −Tr W H n s K W∈H++ −1 o HΣs . (11) G⋆ = HΣs HH + Σn (20) +K +log det Σs At the same time, E(G⋆ ) is given as n n He H He −1 . (12) = max . log det W −Tr W I−G H Σs I−G H E(G⋆ ) = HH Σ−1 H + Σ−1 n
s
Proof: The problem in (10) is a convex problem. To see this, notice that the objective function in (10) is a quadratic function of G with its quadratic terms being given as Tr WGH HΣs HH G + Tr WGH Σn G . (13)
By the identities Tr{AB} = Tr{BA} and Tr{ABCD} = vecT (DT ) CT ⊗ A vec(B), the first term of the above quadratic terms can be rewritten as h i Tr WGH HΣs HH G =vecH(G) W∗⊗ HΣs HH vec(G). (14)
W∈HK , ++
G
oo +K +log det Σs , +GH Σn G
(21)
where the last two steps follow lemma 1 and 2 respectively. Thus the optimization problem (P0) maximizing MI has been transformed into an equivalent problem (P1) in (22), which is shown on the top of next page. As a straightforward consequence of the above two lemmas, we have obtained the optimal solutions to the following two subproblems of (P1). ⋆ When {Fi }L i=1 and G are given, the optimal W is given as W⋆ = arg max . MI W {Fi }L ,G (23) i=1 K
H Notice i s H are both positive semi-definite, h that W and HΣ so W∗⊗ HΣs HH is positive semi-definite [24] and thus the first quadratic term is a convex function of G. Similarly the second quadratic term in (13) can also be proved to be W∈H++ convex function of G. Thus (10) is non-constrained convex −1 L L X X H H problem of G. By setting the derivative with respective to G H H Hi Fi +G Σn G . Hi Fi Σs I−G = I−G to zero [25], we obtain i=1 i=1 h i ∂Tr WE(G) H W= O. (15) G − HΣ = HΣ H +Σ s s n ⋆ When {Fi }L ∂G∗ i=1 and W are given, the optimal G is given
(P1) max . W∈HK , ++
{Fi }L i=1 ,G
h L L i X X H H H H +G Σ G +log det(Σs )+K, H F I−G Σ H F MI {Fi}L ,W,G = log det W −Tr W I−G n i i s i i i=1 (22a)
≤ Pi , s.t. Tr Fi Σs+Σn FH i
i ∈ {1, · · · , L}.
as
G⋆ = arg max .MI G {Fi }L , W i=1
(24)
G
=
L h X i=1
i=1
i=1
L L i−1 X X H Hi Fi Σs Hi Fi +Σn Hi Fi Σs , i=1
i=1
with Σn being given in equation (3). Now we focus on the subproblem of optimizing {Fi }L i=1 with W and G given. Towards this end, we have two options—we can either jointly optimize {Fi }L i=1 in one shot, or we can further consult to BCA methodology again to partition the entire variables {Fi }L i=1 into L blocks, {F1 }, · · · , {FL } and attack L smaller problems one by one in a cyclic manner. For both of these two options, solutions, hopefully in a closed form, are desirable and complexity are concerned. In the following, we discuss these two alternatives in details. A. Jointly Optimizing {Fi }L i=1
The subproblem of (P1) maximizing MI {Fi }L i=1 W, G with W and G given is rewritten as follows h L L X X H Hi Fi Hi Fi Σs I−GH (P2) min . Tr W I−GH {Fi }L i=1
i=1
i H +G Σn G ,
i=1
≤ Pi , i ∈ {1, · · · , L}. s.t. Tr Fi Σs+Σi FH i
(25a) (25b)
The following theorem identifies the convexity of (P2). Theorem 1. The problem (P2 ) is convex. Proof: To begin with, we first look at the function f X : Cm×n 7→ R given as follows: f X , Tr Σ1 XΣ2 XH (26)
with constant matrices Σ1 and Σ2 being positive semidefinite and having appropriate dimensions. By the identity Tr{ABCD} = vecT (DT ) CT ⊗ A vec(B), f X can be equivalently written as (27) f X = vecH X) Σ∗2 ⊗ Σ1 vec X .
Since Σ1 and Σ2 are positive semi-definite, [Σ∗1 ⊗ Σ2 ] is positive semi-definite [24]. Thus f X is actually a convex function with respect to X. PL For a further step, we replace X = i=1 Hi Fi . Since PL H F is an affine (linear actually) transformation of i i i=1
(22b)
variables {Fi }L i=1 , and affine operations preserve convexity by [26], the following function L L n X X H o (28) H F Σ H F = Tr Σ f {Fi }L i i 2 i i 1 i=1 i=1
i=1
is a convex function with respect to variables {Fi }L i=1 jointly. To identify the convexity of the objective in (25a), it suffices to prove the nonlinear terms of {Fi }L i=1 are convex, which are given as Tr
n
GWGH
L X i=1
+
L X i=1
L X H o Hi Fi Hi Fi Σs
Tr
i=1
n
o H H . (29) HH i GWG Hi Fi Σi Fi
Based on the discussion at the beginning of this proof, each of above terms is convex and thus the objective is convex. Similarly the convexity of each power constraint function can also be recognized. Thus the problem (P2) is convex. After identifying the convexity of problem (P2), we reformulate it into a standard quadratic constrained quadratic problem(QCQP) problem. To this end, we introduce the following notations (30a) fi , vec Fi ; g , vec G ; (30b) ∗ H H (30c) Aij , Σs ⊗ Hi GWG Hj ; ∗ (30d) Bi , WΣs ⊗ Hi ; ∗ H H (30e) Ci , Σi ⊗ Hi GWG Hi ; T T f , f1 , · · · , fLT ; (30f) L A , Aij i,j=1 ; (30g) (30h) B , B1 , · · · , BL ; (30i) C , Diag C1 , · · · , CL ; o n ∗ Di,Diag OK(Pi−1Nj), Σs+Σi ⊗INi,OK(PL Nj) ; (30j) j=i+1 j=1 (30k) c , Tr WΣs + σ02 Tr GWGH . Based on the above notations, problem (P2) can be equivalently written as the following QCQP problem, (31a) (P3) : min f H A+C f −2Re gH Bf +c, f
s.t. f H Di f ≤ Pi ,
i ∈ {1, · · · , L}.
(31b)
By theorem 1, (P3) is convex, thus (A + C) is positive 1 semidefinite, which implies that its square root (A+C) 2 exists. Therefore the above problem can be further rewritten in a standard SOCP form as follows: (P3SOCP ) : min . t,
(32a)
f ,t,s H
s.t. s − 2Re{g Bf } + c ≤ t;
(A+C) 12 f
≤ s+1 ; s−1
2
2 1 2
Pi +1
Di2 f , i ∈ {1, · · · , L}.
Pi −1 ≤
2 2
(32b) (32c)
Now the problem of updating one separate beamformer becomes (P3i ) min . fiH Aii+Ci fi − 2Re gH Bi −qH i fi , (33a) fi
s.t. fiH Ei fi ≤ Pi
with the definitions of qi and Ei as follows X ∗ qi , Aij fj ; Ei , Σs+Σi ⊗INi .
(33b)
(34)
j6=i
We introduce the following notations λi,1 (32d) 1 1 − − H .. 2 Ei 2 Aii+Ci Ei 2 = Ui Ui ; (35a) . The above problem can be solved by standard numerical tools λi,KNi like CVX [27]. 1 − 2 pi = UH BH (35b) The method discussed above is summarized in algorithm 1. i Ei i g − qi , Algorithm 1: 3-Block BCA Algorithm to solve (P0) 1
2 3
4
5
6
(0)
Initialization: randomly generate feasible {Fi }L i=1 ; obtain G(0) by (24); obtain W(0) by (23); repeat with G(j−1) and W(j−1) being fixed, solve (P3) in (j) (32), obtain {Fi }L i=1 ; (j) L with {Fi }i=1 and W(j−1) being fixed, obtain G(j) by (24); (j) (j) being fixed, obtain W(j) with {Fi }L i=1 and G by (23); until increase of MI is sufficiently small or predefined number of iterations is reached;
For the proposed 3-BCA algorithm, we have the following conclusion on its convergence Theorem 2. Assume that the covariance matrix Σs ≻ 0. Algorithm 1 generates increasing MI sequence. Its solution sequence has limit points, and each limit point of the solution sequence is a KKT point of the original problem (P0 ). Proof: Refer to appendix A. B. Cyclic (L+1)-BCA Algorithm Although the above proposed 3-block BCA algorithm guarantees a satisfactory convergence, the subproblem (P3) relies on standard numerical solvers, e.g. interior point method [27], to obtain solutions. Closed form solutions to (P3) is unknown. According to the complexity analysis performed in next subsection III-C, when the number of sensors and/or antenna number of each sensor grows, the problem (P3) can be very large size and consequently highly computation demanding. So effective algorithms with lower complexity are desirable. In this subsection, we consult to BCA methodology again to further partition the variables {Fi }L i=1 into L singleton sets: {F1 },· · · ,{FL }. This results in a cyclic (L+1)-BCA algorithm, where only one separate beamformer Fi is optimized at each time and different beamformers are updated in an round robin manner.
i with the eigenvalues {λi,j }KN j=1 arranged in an decreasing order, i.e. λi,1 ≥ · · · ≥ λi,KNi . We denote the k-th element of pi as pi,k and assume that ri = rank Aii+Ci . Then the solution to problem (P3i ) is given by the following theorem.
Theorem 3. Under the assumption that Σs ≻ 0 or Σi ≻ 0, i ∈ {1, · · · , L}, the optimal solution of problem (P3 i ) is given as follows: CASE(I)—if either of the following two conditions holds: i) ∃k ∈ {ri + 1, · · · , KNi } such that |pi,k | 6= 0; PKNi Pri |pi,k |2 > Pi . or ii) k=ri+1 |pi,k | = 0 and k=1 λ2i,k The optimal solution to (P3 i ) is given by −1 H (36) Bi g − qi , fi⋆ = Aii+Ci+µ⋆i Ei
with the positive value µ⋆i being the unique solution to the following equation: ri X |pi,k |2 f (µi ) = = Pi . (37) (λi,k + µi )2 k=1
P i Pri |pi,k |2 CASE(II)—otherwise,i.e. KN k=ri+1 |pi,k | = 0, k=1 λ2i,k ≤ Pi , The optimal solution to (P3 i ) is given by 1 − 1 † − 1 −1 − fi⋆ = Ei 2 Ei 2 Aii+Ci Ei 2 Ei 2 BH i g − qi . (38)
Proof: For limit of space, please refer to Theorem 3 in [8] for detailed proof. In CASE(I) of theorem 3, equation (37) generally has no closed form solution. Notice that f (µi ) is a one-dimension strictly decreasing function of µi . So the determination of µ⋆i can be efficiently performed by a bisection search. Thus a finite interval containing µ⋆i is necessary from which the bisection search can start. The following lemma provides us bounds for µ⋆i . Lemma 3. A lower bound lbdi and upper bound ubdi for positive µ⋆i in (36) in theorem 3 can be given as follows: i) For subcase i) of CASE(I) + kpi k2 kpi k2 lbdi = √ − λi,1 , ubdi = √ ; (39) Pi Pi
ii) For subcase ii) of CASE(I) + kpi k2 − λi,1 , lbdi = √ Pi
with unitary matrix Ui , ui,1 , ui,2 , · · · , ui,KNi having its i columns {ui,j }KN j=1 satisfying the following properties
kpi k2 ubdi = √ − λi,ri , HH H i g Pi , and uH ui,1 = i,j Hi g = 0, for j = 2, · · · , KNi . (44) kHH (40) i gk2 It can be readily checked that the parameter pi in theorem 3 is given as:
where [x]+ , max(x, 0).
Proof: For limit of space, please refer to Lemma 1 in [8] ei )kHH pi,1 = σs2 (1 − gH q i gk2 ; pi,j = 0, j = {2, · · · , KNi }. (45) for detailed proof. By theorem 3 and lemma 3, we have obtained a nearly At this time, the function f (µi ) in (37) reduces to an elegant closed form solution to the problem (P3i ). Here we claim the form 2
2 above solution nearly closed form since it involves a bisection
ei HH σs4 1 − gH q i g 2 search. (46) f (µi ) = 2 . It is worth noting that the fully closed form solution to g µi + (σs2 + σi2 )gH Hi HH i problem (P3i ) does exist in the special but important case Based on the above observations, it can be concluded that of scalar source signal, i.e. K = 1. As will be seen in the the subcase i) of CASE(I) in theorem 3 will never occur. The complexity analysis in subsection III-C, the increase of source two cases for positive and zero µ⋆i can be specified as follows: signal dimension can extensively enlarge the beamforming ⋆ problem size and therefore its complexity. So in practice, when CASE(I)— µi > 0 2 2 ei > (σs2 +σi2 )2 P¯i kHH This is equivalent to σs4 1−gH q the wireless sensor network has adequate bandwidth, it is prei gk2 ferred to transmit the sensed data component by component to and optimal solutions are determined by
1 decrease the processing and beamforming design complexity ⋆ H H 2 ¯− 2
ei Hi g 2 −(σs2 +σi2 ) HH (47a) i g 2, at the fusion center. Another tempting reason to do so is that µi = σs Pi 1−g q −1 ⋆ when the source signal is scalar, µi can be obtained in an f ⋆ = σ 2 (1−gH q H ei ) µ⋆i I+(σs2 +σi2 )HH HH i g. (47b) s i gg Hi explicit way and therefore the bisection search is not needed. i Moreover, at this time, the eigenvalue decomposition (recall CASE(II)— µ⋆i = 0 2 (35a)) will not be involved in solving (P3i ). This conclusion ei ≤ (σs2 + This holds if and only if σs4 1 − gH q 2 ⋆ reads as the following corollary. σi2 )2 P¯i kHH i gk2 and the optimal fi is given by Corollary 1. For scalar transmission (K = 1 ), fully analytic solution to subproblem (P3 i ) can be obtained without evoking bisection search or eigenvalue decomposition. Proof: In the special case where the signal source is scalar, the variables and parameters in the subproblem optimizing one separate beamformer are specialized as follows W → w; Fi → fi ; G → g; Σs → σs2 ; Σi → σi2 .
(41)
ei , j6=i Hj fj , ignoring the terms independent By defining q of fi and omitting the constant positive factor w in the objective, the problem (P3i ) is rewritten as follows H (P3i ) : min . (σs2 + σi2 )fiH HH i gg Hi fi fi H eH (42a) − 2σs2 Re (1 − q i g)g Hi fi P i s.t. kfi k2 ≤ 2 , P¯i . (42b) σs + σi2 Solving the problem (P3i ) just follows the outline of theorem 3. Here, the key point leading to a closed form solution H is the fact that the quadratic matrix HH i gg Hi has rank-1, i.e. ri = 1 in theorem 3. Thus we obtain
H
0 O(KNi −1)×(KNi −1)
ei )HH σs2 (1−gH q i g . 2 H 2 (σs +σi )g Hi HH i g
(48)
Thus we have seen that for scalar transmission case, fully closed form solution to (P3i ) can be obtained without bisection search or eigenvalue decomposition. The cyclic (L+1)-BCA algorithm is summarized in algorithm 2. C. Complexity
P
(σs2 + σi2 )HH ggH Hi 2i (σs + σi2 )gH Hi HH i g = Ui 0
fi⋆ =
(43) UH i ,
In this subsection, we discuss the complexity of the proposed algorithms. The two subproblems optimizing G and W have closed form solutions in (24) and (23), their complexities come from matrix inversion and are given as O K 3 . For the 3-block BCA algorithm, the SOCP problem (P3SOCP ) in (32) is solved by jointly optimizing all beamformers. The complexity of solving an SOCP is [28] kX kX SOC SOC 1 2 3 2 2 nSOC,i , (49) nSOC,i+ O kSOC mSOC+mSOC i=1
i=1
where kSOC is the number of second order cone constraints, mSOC is the dimension of optimization problem and nSOC,i denotes the dimension of the i-th second order cone constraint. problem in (32), PL kSOC = L + 1, mSOC = PL For the , n = K K N SOC,1 i=1 Ni +1 for the first second i=1 i order cone constraint in (32c) and nSOC,i+1 = KNi + 1 for
Algorithm 2: Cyclic (L + 1)-BCA Algorithm to Solve (P0) Optimizing MI 1
2 3 4
5 6 7
8 9
(0)
Initialization: randomly generate feasible {Fi }L i=1 ; obtain G(0) by (24); obtain W(0) by (23); repeat for i = 1; i Pi then else if k=rii +1 |pi,k |2 = 0, k=1 λ2 i,k
10 11
12 13 14 15 16 17 18 19
determine bounds lbd i and ubdi by (40); bisection search on lbdi , ubdi to determine µ⋆i satisfying (37); obtain Fi by (36); else obtain Fi by (38); end update G by (24) ; update W by (23) ; end until increase of MI is sufficiently small or predefined number of iterations is reached;
the i-th power constraint in (32d), i ∈ {1, · · · , L}. Substituting these into (49), the complexity of solving (P3) is √ parameters P 3 O LK 3 ( L , this is also the complexity for each N ) i=1 i loop of 3-block BCA algorithm. For the cyclic (L + 1)-BCA algorithm, the problem (P3i ) optimizing one separate sensor’s beamformer has its major complexity coming from eigenvalue decomposition, 3 3 which i ). Thus the complexity for each loop is P is O(K N L 3 3 O i=1 K Ni . Clearly by fully decomposing the original problem and researching the solution structure of the subproblems, the (L + 1)-block BCA algorithm effectively lowers the computation complexity.
IV. N UMERICAL R ESULTS In this section, numerical results are presented to verify the algorithms proposed in the previous section. In our following experiments, we test the case where the source signal and all observation noise are colored. Specifically, we set the covariance matrices of the source signal and observation noise as Σs = σs2 Σ0 ,
Σi = σi2 Σ0 ,
i ∈ {1, · · · , L},
(50)
where the K × K Toeplitz matrix Σ0 is defined as .. 2 . ρK−1 1 ρ ρ .. .. ρ . . 1 ρ .. 2 Σ0 = ρ2 . . ρ ρ 1 . . . . .. .. .. .. ρ .. 2 K−1 . ρ ρ 1 ρ
(51)
The parameter ρ in the above equation is used to adjust the correlation level between different components of the signal or noise. In our test, ρ is set as ρ = 0.5. Here we define the observation signal to noise ratio at the i-th sensor as SNRi , σs2 σ2 and the channel signal to noise ratio as SNR , σs2 . σi2 0 In figure 2 and 3 we test the performance of the 3-block BCA and cyclic (L + 1)-BCA algorithms for multiple dimension source signal. Here two cases are tested—heterogeneous network and homogeneous network, in figures 2 and 3 respectively. In the heterogeneous network, the transmission power, observation noise level and numbers of antennae of each sensor are different. We set up a wireless sensor network with three sensors, i.e. L = 3. The dimension of the source signal and the number of antenna of the FC are chosen as 3 and 4 respectively, i.e. K = 3 and M = 4. We randomly set the antenna number for each sensor as N1 = 3, N2 = 4, and N3 = 5 respectively, the transmission power constraint for each sensor as P1 = 2, P2 = 2 and P3 = 3 respectively and the observation signal to noise ratio for each sensor as SNR1 = 8dB, SNR2 = 9dB and SNR3 = 10dB respectively. Comparatively, in homogeneous sensor network each sensor has the same transmission power, observation noise level and number of antenna. In this test case we assume that K = 4, M = 4, each sensor has Ni = 5 antennae and transmission power Pi = 2, with observation noise level SNRi = 9dB. In our test, to take into account the impact of the channel parameters, for the above system set-up and any specific channel SNR we randomly generate 500 channel realizations. For each channel realization, two proposed algorithms are run, both of which start from one common random feasible solution. The progress of MI with respect to outer-loop iteration numbers are recorded. For one given iteration number, the average MI performance over all 500 channel realizations is presented in figure 2 and 3. For the implementation of 3-BCA algorithm, SDPT3 solver of CVX is chosen. The blue solid curves represent the average MI performance obtained by 3-block BCA algorithm with different numbers of iterations and the red dotted ones represent those obtained by cyclic (L+1) BCA algorithm. The black dotted curve represents the average MI obtained by random full-power-transmission solutions, which are actually the average MI performance for feasible solutions which make all power constraints active. From figures 2 and 3, we see that the optimized beamformers obtained by the proposed algorithms present significant MI improvement compared to nonoptimized beamformers. Usually 40 to 50 iteration loops are sufficient to make the two algorithms
3−Block BCA and Cyclic (L+1) BCA with Different Numbers of Iterations−Heterogeneous Network 10
3−Block BCA and Cyclic (L+1) BCA with Different Numbers of Iterations−Scalar Source 3
2.8 Mutal Information(Nats per use of channel)
Mutal Information(Nats per use of channel)
9
8 9.1 7
9 8.9 8.8
6
3−blk BCA, 3 Iters 3−blk BCA, 7 Iters 3−blk BCA, 20 Iters 3−blk BCA 40 Iters Cyc. (L+1)−BCA, 3 Iters Cyc. (L+1)−BCA, 7 Iters Cyc. (L+1)−BCA, 20 Iters Cyc. (L+1)−BCA, 40 Iters Random Initials
8.7 8.6 5 8.5 4 4
0
2
4.5
5
4
5.5
6
6
8
10
12
14
16
2.6
2.95
2.4
2.9 2.85
2.2
2.8 2
2.75 2.7
1.8 −6
18
0 −4
−2
0
2
SNR(σ2/σ2) s
0
s
0
Fig. 4: Scalar Source Signal Case: 3-Block BCA Algorithm and Cyclic (L + 1)-Block BCA Algorithm with Different Numbers of Iterations.
3−Block BCA and Cyclic (L+1) BCA with Different Numbers of Iterations−Homogeneous Network 15
3 Block BCA and Cyclic (L+1) BCA Algorithms Starting from Different Initial Points−Test Case 1 9.5
14
9
13
Mutual Information(Nats per use of channel)
Mutal Information(Nats per use of channel)
6
SNR(σ2/σ2)
Fig. 2: Heterogenous Test Case: 3-Block BCA Algorithm and Cyclic (L+1)-Block BCA Algorithm with Different Numbers of Iterations.
12 11
3−blk BCA, 3 Iters 3−blk BCA, 7 Iters 3−blk BCA, 30 Iters 3−blk BCA 50 Iters Cyc. (L+1)−BCA, 3 Iters Cyc. (L+1)−BCA, 7 Iters Cyc. (L+1)−BCA, 30 Iters Cyc. (L+1)−BCA, 50 Iters Random Initials
14.2 14
10
13.8
9
13.6 8 13.4 7
8.5
9.3 9.2
8
9.1 9
7.5
8.9 8.8
7
8.7 8.6
6.5
8.5 8.4
6
8.3
13.2
2
5.5
6
4
0
2
4
7 6
8 8
10 2
6
8
10
12
14
16
18
20
3−blk BCA, Starting from One Specific Initial Point Cyc. (L+1)−BCA, Starting from One Specific Initial Point
13 6
5
4
3−blk BCA, 3 Iters 3−blk BCA, 10 Iters 3−blk BCA, 30 Iters 3−blk BCA 50 Iters Cyc. (L+1)−BCA, 3 Iters Cyc. (L+1)−BCA, 10 Iters Cyc. (L+1)−BCA, 30 Iters Cyc. (L+1)−BCA, 50 Iters Random Initials 2 4 6 8 10 12
12
14
16
18
2
SNR(σs /σ0)
5
0
10
20
30
40 Iteration Index
50
60
70
80
Fig. 3: Homogeneous Test Case:3-Block BCA Algorithm and Cyclic (L+1)-Block BCA Algorithm with Different Numbers of Iterations
Fig. 5: Heterogenous Test Case: Optimizing Generalized MI by 3-Block BCA Algorithm and Cyclic (L + 1)-Block BCA Algorithm with Different Initial Points
converge and the two algorithms finally converge to almost identical MI performance. In figure 4 we test the special case of scalar source signal (K = 1), where (L+1)-block BCA algorithm has fully closed form solution, which is summarized in corollary 1. In this experiment, we have the system setup as follows M = 4, N1 = 3, N2 = 4, N3 = 5, P1 = 1, P2 = 2, P3 = 3 SNR1 = 7dB, SNR2 = 8dB and SNR3 = 9dB. Similar results as in the multiple dimension source signal case have been obtained. In figure 5 and 6, we check the impact of the random initials to the proposed algorithms. We use the same system setup as those in figure 2 and 3 respectively. Here the channel param-
eters are randomly chosen and fixed. 10 feasible solutions, each of which makes all the power constraints active(satisfied with equality) are randomly generated. For each random initial point, we invoke the 3-block BCA and cyclic (L + 1)block BCA algorithms to optimize the beamformers. The MI progress of the proposed algorithms with different initials are illustrated in 5 and 6. It can be seen that both the 3-block BCA and cyclic (L+1)-block BCA are rather insensitive to selection of initial points. The two different algorithms with different initial points finally converge to almost identical value. Last we test the complexity of the proposed algorithms. In Table I the average MATLAB running time for each out-loop is
3 Block BCA and Cyclic (L+1) BCA Algorithms Starting from Different Initial Points−Homogeneous Network 14
Mutual Information(Nats per use of channel)
13 12 11
13
10 12.5 9 12
8 7
11.5
6 11 2
5 4
4
6
8
10
12
14
16
18
20
3−blk BCA, Starting from One Specific Initial Point Cyc. (L+1)−BCA, Starting from One Specific Initial Point 0
10
20
30
40 Iteration Index
50
60
70
80
Fig. 6: Homogeneous Test Case: Optimizing Generalized MI by 3-Block BCA Algorithm and Cyclic (L + 1)-Block BCA Algorithm with Different Initial Points TABLE I: MATLAB Running Time Per Outer-Loop (in Sec.) PP P L Algorithms L=5 L = 10 L = 20 Dim. PP P K=1 3-BCA M = Ni = 2 Cyc. (L+1) K=4 3-BCA M = Ni = 4 Cyc. (L+1) K=8 3-BCA M = Ni = 8 Cyc. (L+1) Notes: SDPT3 solver of CVX is
0.2748 0.0125 0.5079 0.0319 10.41 0.0488 chosen to
0.4018 0.0446 2.168 0.0999 90.51 0.1310 implement
0.7627 0.1792 12.93 0.3761 729.2 0.4747 3-BCA.
being fixed, the objective value obtained by solving the current sub-problem cannot be smaller than previous one. Thus the entire MI sequence keeps increasing. Under the positive definiteness assumption of Σs , Σs+ Σi ≻ 0. Thus ∀i ∈ {1, · · · , L} we have (52) kFi k2F λmin Σs+Σi ≤ Tr{Fi Σs+Σi FH i } ≤ Pi ,
where λmin (·) denotes the minimum eigenvalue of a Hermitian matrix. Since λmin Σs+Σi > 0, kFi k2F is finite for all i. Thus the variable {Fi }L i=1 is bounded. By Bolzano-Weierstrass the(kj ) L }i=1 orem, there exits a subsequence {kj }∞ j=1 such that {Fi converges. Since G and W are updated by continuous func- (kj ) L tions of {Fi }L }i=1 , W(kj ) , G(kj ) i=1 in (24) and (23), {Fi converges. Thus the existence of limit points in solution sequence has been proved. ¯ i }L , W, ¯ G ¯ is any limit point We assume that {F i=1 (k) (k) , G(k) . Then there exists a subseof {Fi }L i=1 , W j→∞ (k ) (kj ) quence {kj } such that {Fi j }L , G(kj ) −→ i=1 , W (k) L L ¯ ¯ ¯ {Fi }i=1 , W, G . Since {Fi }i=1 is bounded, by possibably restricting can assume that to a subsequence, we L (k +1) bi . F {Fi j }L i=1 converges to a limit i=1 (k +1)
Since for each j, {Fi j }L i=1 are feasible, i.e. (k +1) H (k +1) Tr{Fi j Σs+Σi Fi j )} ≤ Pi , i ∈ {1, · · · , L}. (53)
By taking j → ∞ in the above inequalities, we obtain H b i )} ≤ Pi , i ∈ {1, · · · , L}. b i Σs+Σi F (54) Tr{F L bi are feasible. So F i=1 L For any feasible Fi i=1 , we have presented. Here we consider the homogeneous sensor network. (k ) (k ) (k ) (k ) Different values of K, L and Ni are tested, which result in (k +1) L j j j ≤ MI {Fi j }L ,G j . (55) i=1 W different sizes of problem. SDPT3 solver of the CVX is used MI {Fi }i=1 W ,G to implement 3-BCA algorithm. As shown in the table, the Noticing that MI function is continuous and taking j → ∞ in cyclic (L+1) BCA algorithm is highly efficient, since each of the above, we obtain its update step can be performed in an almost analytical way. W, ¯ G ¯ ≤ MI {F b i }L W, ¯ G ¯ , MI {Fi }L (56) i=1 i=1 V. C ONCLUSION In this paper, we consider the linear beamforming design for any feasible {Fi }L i=1 . (k) L problem for a coherent MAC wireless sensor network to maxNotice the {Fi }i=1 generated by algorithm 1 are feasible, ¯ i }L are imize the mutual information. As we have seen, the original by continuity of power constraint functions, {F i=1 problem is nonconvex and difficult. To solve this problem, feasible. Thus we have we adopt the weighted minimum mean square error method ¯ i }L W, ¯ G ¯ ≤ MI {F b i }L W, ¯ G ¯ . MI {F (57) i=1 i=1 and block coordinate ascent method to decompose the original difficult problem into subproblems and examine the solution At the same time, since the MI sequence is increasing and to each subproblem, especially their closed form solution. ¯ i }L , W, ¯ G ¯ is a limit point of the solution sequence, {F i=1 The complexity and convergence of proposed algorithms are also discussed in details. Extensive numerical results are ¯ i }L W, ¯ G ¯ ≥ MI {F(k) }L W, ¯ G ¯ , MI {F (58) i=1 i=1 i presented to verify and compare the behaviors of the proposed for any integer k. Substitute k with kj in (58), take limit algorithms. ¯ i }L j → ∞ and combine it with (57), we have shown that {F i=1 A PPENDIX is actually an optimal solution to the problem (P2) with A. Proof of Theorem 2 ¯ and G. ¯ So {F ¯ i }L satisfy KKT conditions parameters W i=1 ¯ and G, ¯ which are listed in (59) Proof: Since for each sub-problem, we solve an optimiza- of (P2) with parameters W tion problem with respect to a subset of variables with others shown on the top of next page.
L X H ¯ ¯ ¯ ¯H ¯ ¯ ¯ ¯ ¯ Σs +HH F H −HH G W I− G i i i GWG Hi Fi Σi +λi Fi Σs+Σi = O, i i=1
n Ho ¯ ¯ i Σs+Σi F = 0, − P λi Tr F i i n Ho ¯ i ≤ Pi , ¯ i Σs+Σi F Tr F λi ≥ 0,
To simplify the following exposition, we introduce the following two notations: ¯ , H
L X
¯ i; Hi F
(60a)
i=1
¯ n , σ2 I + Σ 0
L X
¯ H HH . ¯ i Σi F Hi F i i
(60b)
i=1
According to the update step in algorithm L 1, the limit points ¯ and G ¯ have the relations with F ¯i W as follows. i=1 ¯ s, ¯ H+Σ ¯ n −1 HΣ ¯ = HΣ ¯ sH (61a) G H ¯ −1 ¯ −1 ¯ ¯ W =H Σ H+Σ . (61b) n
s
Utilizing (61) we can prove two identities in (62) and (63) respectively on the top of next page. Substituting equations (62) and (63) into (59a), we can rewrite the first order KKT conditions associated with only ¯ i }L as in equation (64) shown in the next page. {F i=1 To check the conditions of the original problem (P0), we need to determine the derivative of its Lagrangian function, or equivalently the derivative of MI with respect to {Fi }. By defining H,
L X
Hi Fi ,
(65)
i=1
the derivative of MI is calculated in (66) (shown in next page) with C1 (dFi ) and C2 (dFi ) being uninteresting terms involved dFi only and independent of d(F∗i ). By comparing the equations (64) with the derivative in (66b), it is easily to recognize that (64) is actually the first order KKT condition of problem (P0) optimizing MI. Together with equations (59b), (59c) and (59d), the KKT conditions of ¯ i }L . original problem have been proved to be satisfied by {F i=1 Thus the proof is complete. R EFERENCES [1] Y. Zhu, E. Song, J. Zhou, and Z. You, “Optimal dimensionality reduction of sensor data in multisensor estimation fusion,” IEEE Trans. Signal Process., vol. 53, no. 5, pp. 1631-1639, May 2005. [2] J. Fang and H. Li, “Optimal/near-optimal dimensionality reduction for distributed estimation in homogeneous and certain inhomogeneous scenarios,” IEEE Trans. Signal Process., vol. 58, no. 8, pp. 4339-4353, 2010.
i ∈ {1, · · · , L}; i ∈ {1, · · · , L}; i ∈ {1, · · · , L}; i ∈ {1, · · · , L}.
(59a) (59b) (59c) (59d)
[3] J. Fang and H. Li, “Power constrained distributed estimation with cluster-based sensor collaboration,” IEEE Trans. Wireless Commun., vol. 8, no. 7, pp. 3822-3832, July 2009. [4] I. D. Schizas, G. B. Giannakis, and Z.-Q. Luo, “Distributed estimation using reduced-dimensionality sensor observations,” IEEE Trans. Signal Process., vol. 55, no. 8, pp. 4284-4299, Aug. 2007. [5] Y. Liu, J. Li, X. Lu and C. Yuen, “Optimal linear precoding and postcoding for MIMO multi-sensor noisy observation problem,”, IEEE International Conference on Communications(ICC), Sydney , Jun. 2014. [6] J. Xiao, S. Cui, Z. Luo, and A. J. Goldsmith, “Linear coherent decentralized estimation,” IEEE Trans. Signal Process., vol. 56, no. 2, pp. 757-770, Feb. 2008. [7] A. S. Behbahani, A. M. Eltawil, H. Jafarkhani, “Linear decentralized estimation of correlated data for power-constrained wireless sensor networks,” IEEE Trans. Signal Process., vol. 60, no. 11, pp. 6003-6016, Nov. 2012. [8] Y. Liu, J. Li, X. Lu and C. Yuen, “Joint transceiver design for wireless sensor networks through block coordinate descent optimization,” available online: http://arxiv.org/abs/1409.7122. [9] J. Fang, H. Li, Z. Chen and Y. Gong, “Joint precoder design for distributed transmission of correlated sources in sensor networks,”, IEEE Trans. Wireless Commun., vol. 12, no. 6, pp. 2918-2929, June 2013. [10] Y. Liu, J. Li, and X. Lu, “Transceiver design for clustered wireless sensor networks — Towards SNR Maximization,” available online: http://arxiv.org/abs/1504.05311. [11] Y. Liu, J. Li, and X. Lu “Multi-terminal joint transceiver design for MIMO systems with contaminated source and individual power constraint,” IEEE Intl Symposium on Information Theory(ISIT), Honolulu, Jun. 2014. [12] S. Serbetli and A. Yener, “Transceiver optimization for multiuser MIMO systems,” IEEE Trans. Signal Process., vol. 52, pp. 214-226, Jan. 2004. [13] M. R. A. Khandaker and Y. Rong, “Joint transceiver optimization for multiuser MIMO relay communication systems”, IEEE Trans. Signal Processing, vol. 60, pp. 5977-5986, Nov. 2012. [14] Y. Rong, X. Tang, and Y. Hua, “A unified framework for optimizing linear non-regenerative multicarrier MIMO Relay Communication Systems”, IEEE Trans. Signal Process., vol. 57, pp. 4837-4851, Dec. 2009. [15] S. S. Christensen, R. Argawal, E. de Carvalho, and J. M. Cioffi, “Weighted sum-rate maximization using weighted MMSE for MIMO-BC beamforming design,” IEEE Trans. Wireless Commun., vol. 7, no. 12, pp. 1-7, Dec. 2008. [16] Q. Shi, M. Razaviyayn, Z. Luo, and C. He, “An iteratively weighted MMSE approcah to distributed sum-utility maximization for a MIMO iterfering broadcast channel,” IEEE Trans. Signal Process., vol. 59, no. 9, pp. 4331-4340, Dec. 2011. [17] I. E. Telatar, “Capacity of multi-antenna Gaussian channels,” Eur. Trans. Telecommun., vol. 10, no. 6, pp. 585-595, Nov.-Dec. 1999.
−1 ¯ H+Σ ¯ n −1 H ¯ ¯ sH ¯ H HΣ ¯ HΣ ¯ −1 H+Σ ¯ ¯ H+Σ ¯ n −1 HΣ ¯ s H ¯ sH ¯W ¯ I−G ¯ HH ¯ = HΣ I−Σs H G n s h i −1 ¯ H HΣ ¯ ¯ HΣ ¯ −1 H+Σ ¯ H+Σ ¯n − H ¯ sH ¯ HΣ ¯ −1 HΣ ¯ s H ¯ H+Σ ¯ n −1 HΣ ¯ ¯ sH ¯ H+Σ ¯ n −1 H+ ¯ sH ¯ ¯ H+Σ ¯ n −1 H ¯ sH H Σ HΣ = HΣ s n s n | {z } ¯ H+Σ ¯n ¯ sH = HΣ
−1
=O
¯ H
(62)
−1 ¯ H+Σ ¯ n −1 ¯ sH ¯ H HΣ ¯ HΣ ¯ −1 H+Σ ¯ ¯ s H Σs H HΣ n s i h ¯ H+Σ ¯ n −1 ¯ H HΣ ¯ sH ¯ H+HΣ ¯ sH ¯ HΣ ¯ −1 HΣ ¯ sH ¯ H+Σ ¯ n −1 HΣ ¯ sH ¯ sH = HΣ n i h ¯ H+Σ ¯ n −1 ¯ sH ¯ H+Σ ¯ n HΣ ¯ sH ¯ HΣ ¯ −1 HΣ ¯ sH ¯ H+Σ ¯ n −1 HΣ ¯ sH = HΣ n ¯ HΣ ¯ −1 ¯ sH ¯ H+Σ ¯ n −1 HΣ ¯ sH = HΣ n
¯ H+Σ ¯n ¯W ¯G ¯ H = HΣ ¯ sH G
HH i
h
−1
L L L L X X h i−1 X X H ¯ i Σs I ¯i H ¯ i Σs ¯H ¯ i Σi F Hi F Hi F Hi F σ02 I+ Hi F i Hi + i=1 L X
−
i=1
¯i Hi F
n
σ02 I +
L X
i=1
i=1
i=1
H
H −1
¯HH ¯ i Σi F Hi F i i
i=1
i
¯ i Σs+Σi = O, ¯ i Σi − λi F Hi F
i ∈ {1, · · · , L};
o d HΣs HH Σ−1 n io n −1 h H −1 H −1 HΣ d(H )Σ +HΣ H d(Σ ) + C1 (dFi ) = Tr I+HΣs HH Σ−1 s s n n n h n i o H −1 H HΣs I − HH Σ−1 = Tr HH + C2 (dFi ), n Hi Fi Σi d(Fi ) i Σn+HΣs H
d(MI) = Tr
⇒
(63)
I+HΣs HH Σ−1 n
(64)
−1
(66a)
L L L L h X X h H i−1X X ∂MI H H 2 H H F H F Σ H F + H F Σ F H σ I+ = H i i Σs I− i i s i i i i i i i 0 i ∂F∗i i=1 i=1 i=1 i=1 L X i=1
L i H 2 X H −1 Hi Fi Σi , i ∈ {1, · · · , L}. (66b) σ0 I+ Hi Fi Σi FH Hi Fi i Hi
[18] S.-Y. Chung, “On the construction of some capacityapproaching coding schemes,” Ph.D. dissertation, Dept. EECS, Massachusetts Institute of Technology, 2000. [19] Y. Hu, J. Garcia-Frias, and M. Lamarca, “Analog joint sourcechannel coding using non-linear curves and MMSE decoding,” IEEE Trans. Commun., vol. 59, no. 11, pp. 3016-3026, Nov. 2011. [20] T. M. Cover and J. A. Thomas Elements of Information Theory, 2nd ed., John Wiley & Sons, 2006. [21] D. Bertsekas and J. Tsitsiklis, Parallel and Distributed Computation, Englewood Cliffs, NJ: Prentice-Hall, 1989. [22] D. Bertsekas, Nonlinear Programming, 2nd ed., Belmont, MA:Athena Scientific, 1999. [23] L. Grippo, and M. Sciandrone, “On the convergence of the block nonlinear Gauss-Seidel method under convex constraints,” Operations Research Letters, vol. 26, pp. 127-136, 2000. [24] T. P. Minka. Old and New Matrix Algebra Useful for Statistics,
i=1
Notes, December 2000. [25] A. Hjorungnes, Complex-Valued Matrix Derivatives: With Applications in Signal Processing and Communications. Cambridge University Press, 2011. [26] S. Boyd, and L. Vandenberghe, Convex Optimization. New York: Cambridge University Press, 2004. [27] M. Grant and S. Boyd, “CVX: Matlab software for disciplined convex programming (web page and software),” [Online]. Available: http://cvxr.com/cvx Apr. 2010 [28] I. Polik and T. Terlaky, “Interior Point Methods for Nonlinear Optimization,” in Nonlinear Optimization, 1st edition, G. Di Pillo, F. Schoen, editors. Springer, 2010.
