extfieldsofstudy:
Computer Science
Mathematics

s2fieldsofstudy:
Computer Science
Mathematics

text:
Quickest Change Point Detection and Identification Across a Generic Sensor Array

In this paper, we consider the problem of quickest change point detection and identification over a linear array of $N$ sensors, where the change pattern could first reach any of these sensors, and then propagate to the other sensors. Our goal is not only to detect the presence of such a change as quickly as possible, but also to identify which sensor that the change pattern first reaches. We jointly design two decision rules: a stopping rule, which determines when we should stop sampling and claim a change occurred, and a terminal decision rule, which decides which sensor that the change pattern reaches first, with the objective to strike a balance among the detection delay, the false alarm probability, and the false identification probability. We show that this problem can be converted to a Markov optimal stopping time problem, from which some technical tools could be borrowed. Furthermore, to avoid the high implementation complexity issue of the optimal rules, we develop a scheme with a much simpler structure and certain performance guarantee.

I. Introduction
The standard quickest change detection problem is set to detect some unknown time point at which certain signal probability distribution changes over a sequence of observations. Recently, with the development of wireless sensor networks, multiple sensors can be deployed to execute the quickest change detection, and the sensors can send quantized or unquantized observations or certain local decisions to a control center, who then makes a final decision [1], [2]. Most of the existing work is based on an assumption that the statistical properties of observations at all sensors change simultaneously. However, in certain scenarios, this assumption may not hold well. For instance, when multiple sensors are used to detect the occurrence of the chemical leakage, the sensors that are closer to the leakage source usually observe the change earlier than those far away from the source. In such cases, two interesting problems arise: one is to detect the change as soon as possible; the other is to identify which sensor is the closest to the source, such that we could have a first-order inference over the leakage source location.
As far as we know, currently there are few work studying the case of change occurring non-simultaneously. In the related work, the authors in [3] proposed a scheme that each sensor makes a local decision with the computing burden at the local sensors, where they did not consider the identification problem. In [4], the authors modeled the change propagation process as a Markov process to derive the optimal stopping rule and assumed that the change pattern always first reaches a predetermined sensor, such that the identification problem is ignored. In [5], the identification problem for the special case of two sensors was studied, where the sufficient statistic is proven as a Markov process and a joint optimal stopping rule and terminal decision rule are proposed.
In this paper, we study the joint change point detection and identification problem over a linear array of N sensors, where the change first occurs near an unknown sensor, then propagates to sensors further away. We assume that all sensors send their observations to a control center. With the sequential observation signals, the control center first operates a stopping rule to decide when to alarm that the change has occurred; then the control center deploys a terminal decision rule to determine which sensor that the change pattern reaches first. In our setup, three performance metrics are of interest: i) detection delay, which is the time interval between the moment that the change occurs and the moment that an alarm is raised; ii) false alarm probability, which is the probability that an alarm is raised before the actual change occurs; and iii) false identification probability, which is the probability that the control center does not correctly identify the sensor that the change pattern first reaches. We apply the Markov optimal stopping time theory to design the optimal decision rules to minimize a weighted sum of the above three metrics. Furthermore, we derive a scheme with a much simpler structure and certain performance guarantee.
The rest of this paper is organized as follows. In Section II, we introduce the system model. In Section III, we derive the optimal decision rules. In Section IV, we propose a scheme approximate to the optimal decision rules with a much lower complexity. In Section V, we present some numerical results, with conclusions in Section VI.

II. System Model
We consider a scenario with N sensors constructing a linear array to monitor the environment, as shown in Fig. 1. At an unknown time point, a change occurs at an unknown location and propagates, where we use change point time Γ i to denote the time that the change pattern reaches sensor i. We further use S to denote the index of the sensor that the change pattern first reaches. We focus on the Bayesian setup and use p i to denote the prior probability of {S = i}, simply with p i = 1/N. Conditioned on the event that the change pattern first reaches sensor i, Γ i is assumed to bear a geometric distribution [4], [5] with parameter ρ, 0 < ρ ≤ 1 i.e., where k denotes the discretized time and takes integer values. We consider the practical factors in the environment, such as the wind or the blockers, which will affect the propagation speed of the change. For instance, see in Fig. 1, if the direction of the wind is from the left to the right side in the monitored scenario, then the propagation of the air pollution will be much faster at the right side of sensor S = i than that of the left side. And at the same side, the propagation follows the deterministic order shown as We further assume that after the change patten reaches the first sensor i, for the right side of sensor i, it will propagate from one sensor to another sensor following the geometric propagation models as while for the left side of sensor i, the propagation follows where ρ 1 and ρ 2 are used to model possibly different propagation speed along each direction, e.g., ρ 1 > ρ 2 means the propagation speed is higher at the right side that that of the left side.
Taking above assumption, for S = i, we define all possible i × (N − i + 1) + 1 events at time k as follows: where T 1,k ∼ T N−i+1,k denote the events that after the change pattern first reaching sensor i, it propagates across the sensors sequentially at the right side of sensor i. The number of events equals to the number of sensors at the right side plus 1, i.e., ,k denote the events that after the change pattern reaching sensor i and i − 1, it propagates across the sensors sequentially at the right side of sensor i, and the number of events is also N − i + 1, which is the same for the case that after the change pattern reaching sensor i, i − 1, i − 2 and so on. Since there are i−1 sensors at the left side of sensor i and the event of no change pattern reaches any sensor is T 0,k , the total number of possible events is i × (N − i + 1) + 1.
At each time k, we assume that the observations z k = [z 1,k , ..., z N,k ] from all sensors are available at a control center. For each sensor i, conditioned on Γ i , z i,k is Identically and Independently Distributed (IID) according to f 0 before Γ i , and IID according to f 1 after Γ i , i.e., The observation sequence We use P to denote the probability measure that specifies the prior distribution of {S = i}, the distribution of change point time, and the distribution of {z k ; k = 1, 2, ...}. We also use E to denote the expectation under the probability measure P. Specifically, we use P i to denote the probability measure when S = i.
With above setups, the control center needs to detect the earliest change point time Γ = min{Γ 1 , ..., Γ N } as soon as it occurs. A stopping time τ will be decided for when to stop sampling and alarm that a change has occurred, where a false alarm may happen if τ < Γ. We target to minimize the averaged detection delay E{(τ − Γ) + } with keeping the false alarm probability P[τ < Γ] small. In addition, we also require the control center to identify which sensor the change pattern reaches first. We adopt δ k to denote the F k -measurable terminal decision rule used by the control center to make the identification, andŜ to denote the index of the sensor identified, i.e.,Ŝ = δ k (F k ). A false identification occurs if S S , such that we also want to keep P[Ŝ S ] small. We use δ = {δ 1 , δ 2 , ...} to denote the sequence of terminal decision rules. Summarizing above, our goal is to design a stopping time τ and a terminal decision rule δ that minimize the aggregated risk function defined as where c 1 and c 2 are appropriate constants that balance the three costs.

III. Optimal Rules
In this section, we present the optimal stopping and terminal decision rules. To proceed, we define the following posterior probabilities at time k: We also define the following matrices and vectors: For each i with the number of events less than M, the extra elements in π k are set as 0.
With the posterior probabilities defined above, we denote We first have the following theorem regarding the optimal terminal decision rule δ.
Theorem 1: For any stopping time τ, the optimal terminal decision rule isŜ and we have The proof follows from Proposition 4.1 of [6]. Theorem 1 implies that the optimal terminal decision rule is simply to choose the sensor that has the largest posterior probability. A similar situation also arises in the multiple hypothesis testing problem considered in [7].
Using above optimal terminal decision rule, we can further express the optimization objective in (7) as a function of the posterior probabilities defined in (8) and (9), as shown below.
Lemma 1: For any stopping time τ, (7) can be written as Proof: Based on the Bayesian's rule, we have
Furthermore, we have the following lemma regarding q k .

Lemma 2:
There is a time-invariant function g such that q k = g(q k−1 , z k ).

Proof: We have
in which Since the geometric distribution based model owns the memoryless property, the transition probabilities of the events are shown in Fig. 2. According to these transition probabilities, we have Therefore, π k can be computed by q k−1 and z k .
For another element p k in q k , we have where which can be calculated by using (18) and (19). Hence, p k can also be computed by q k−1 and z k . This lemma implies that the posterior probabilities q k can be recursively computed from q k−1 and z k . Combined with Lemma 1, we know that q k is a sufficient statistic for the problem of minimizing (14). Thus, the problem at the hand is a Markov stopping time problem.
Therefore, we could borrow results from the optimal stopping time theory to design the optimal decision rules for our problem. We first consider a finite time horizon case, in which one has to make a decision before a deadline T , i.e., τ ≤ T . It is easy to check that the cost-to-go functions are where Applying the optimal stopping time theory [6], we have the following theorem for the optimal decision rules.
Theorem 2: The optimal stopping time is obtained as with the optimal terminal decision rule is given in (12). In the infinite time horizon case when T → ∞, we have since we have J T k (q) > 0, J T +1 k (q) ≤ J T k (q), and the fact that all strategies allowed with deadline T are also allowed with deadline T + 1. Since the observations are memoryless and conditionally IID, J ∞ k (q) is the same for all k; we then use J(q) to denote J ∞ k (q). Thus, A(q) is derived as in which the interchange of lim and is allowed due to the dominated convergence theorem. Therefore, when the deadline is infinite, the optimal stopping rule becomes with the optimal terminal decision rule is given in (12).

IV. Approximation to The Optimal Stopping Rule
When N is large, the optimal stopping rule does not have a simple structure, which makes the implementation highly costly. In this section, we propose a much simpler rule which approximates to the optimal stopping rule.
The proof follows from page 477 of [8], by using Fatou's lemma. We can use Lemma 3 to derive the following approximation of the optimal stopping rule.
Theorem 3: In the asymptotic case of the rare change occurring with ρ → 0, one approximation of the optimal stopping rule has the following simple structure And we use the optimal terminal decision rule specified in (12). Proof: The proof follows closely with the proof of Theorem 2 of [4]. First, we have where i k = arg max For the first part of (29), after interchanging the integral and sum, by using (17), (19), (20), and (21), we have For the second part of (29), according to Lemma 3, Plugging the above two results (30) and (31) into (29), we have In the sequel, we assume that A T T −1 (q T −1 ) equals to the right side of (32).
We define the following transformation as Then, J T T −1 can be rewritten as We define Φ k and Ψ k as Then straightly we see that Φ T = 0, Ψ T −1 = 0, and For the next steps, we follow the proof of Theorem 2 of [4], which is skipped here. And additionally we use Lemma 3. Finally, it can be derived that And the test structure reduces to stopping when Therefore, we have the structure of the stopping rule as stated in Theorem 3.
Regarding to Theorem 3, we have several notes as follows. 1) From Lemma 3 and Theorem 2, we see that τ app is a lower bound of the optimal stopping time, i.e. τ app ≤ τ opt , in the case of ρ → 0. The supermartingale property shown in Lemma   close, which indicates that τ app would be close to τ opt .
2) From (15) and (33), we have the testing statistic This structure conforms to the well-known Shiryaev's procedure [9], which is the optimal stopping rule for single sensor with IID observations and Bayesian setting.

V. Numerical Simulation
Given that it is hard to efficiently compute the solution structure in (23), we compute the approximate optimal stopping rule in (28) and simulate its performance. We assign 5 nodes constructing a linear sensor array and assume that f 0 ∼ N(0, 1) and f 1 ∼ N(1, 1). The change point time is generated according to the geometric distribution with ρ = 0.01, ρ 1 = 0.1 and ρ 2 = 0.05, respectively. According to (15), the false alarm probability with τ app is Thus we have c 1 = ρα 1−α , where α is the maximum allowance for the false alarm probability, which could determine the required select c 1 value.
In Fig. 3, we illustrate the relationships among the false alarm probability, the false identification probability, and the averaged detection delay. We see that as the averaged detection delay increases, the false alarm probability decreases. When the averaged detection delay becomes large, the false identification probability does not decrease much and a probability floor appears, which is due to the fact that only the samples between the time when the change pattern reaches the first sensor and the time when it reaches the second sensor can be used to effectively distinguish the sensor that the change pattern first reaches. Since this part of the samples is limited, which will not increase with the detection delay, a false identification probability floor exists. In Fig. 4, we draw the posterior probability p i k over time, where we assume that the change pattern first reaches node 3, and then propagates to node 4. We see that as time goes, p 3 k gradually becomes larger than the others, which indicates that node 3 should be identified. In Fig. 5, we show the relation between max{p 1 k , ..., p N k } and E{max{p 1 k+1 , ..., p N k+1 }|F k } in (31). Since (31) is the key in deriving the our simplified rule, the fact that these two curves are close suggests that the performance of our low-complexity rule might be close to that of the optimal stopping rule in (23) and (26).

VI. Conclusions
We have studied the quickest change point detection problem and the closest-node identification problem over a sensor


text:
Linear Precoding and Equalization for Network MIMO with Partial Cooperation

A cellular multiple-input multiple-output (MIMO) downlink system is studied in which each base station (BS) transmits to some of the users, so that each user receives its intended signal from a subset of the BSs. This scenario is referred to as network MIMO with partial cooperation, since only a subset of the BSs are able to coordinate their transmission towards any user. The focus of this paper is on the optimization of linear beamforming strategies at the BSs and at the users for network MIMO with partial cooperation. Individual power constraints at the BSs are enforced, along with constraints on the number of streams per user. It is first shown that the system is equivalent to a MIMO interference channel with generalized linear constraints (MIMO-IFC-GC). The problems of maximizing the sum-rate(SR) and minimizing the weighted sum mean square error (WSMSE) of the data estimates are non-convex, and suboptimal solutions with reasonable complexity need to be devised. Based on this, suboptimal techniques that aim at maximizing the sum-rate for the MIMO-IFC-GC are reviewed from recent literature and extended to the MIMO-IFC-GC where necessary. Novel designs that aim at minimizing the WSMSE are then proposed. Extensive numerical simulations are provided to compare the performance of the considered schemes for realistic cellular systems.

I. INTRODUCTION
Interference is known to be major obstacle for realizing the spectral efficiency increase promised by multiple-antenna techniques in wireless systems. Indeed, the multiple-input multiple-output (MIMO) capacity gains are severely degraded and limited in cellular environments due to the deleterious effect of interference [1], [2]. Therefore, network-level interference management appears to be of fundamental importance to overcome this limitation and harness the gains of MIMO technology. Confirming this point, multi-cell cooperation, also known as network MIMO, has been shown to significantly improve the system performance [3]. Network MIMO involves cooperative transmission by multiple base stations (BSs) to each user. Depending on the level of multi-cell cooperation, network MIMO reduces to a number of scenarios, ranging from a MIMO broadcast channel (BC) [4] in case of full cooperation among all BSs, to a MIMO interference channel [5], [6] in case no cooperation at the BSs is allowed. In general, network MIMO allows cooperation only among a cluster of BSs for transmission to a certain user [7], [8] (see also references in [3]).
In this paper, we consider a MIMO interference channel with partial cooperation at the BSs. It is noted that all BSs cooperating for transmission to a certain user have to be informed about the message (i.e., the bit string) intended for the user. This can be realized using the backhaul links among the BSs and the central switching unit. We focus on the sumrate maximization (SRM) and on the minimization of weighted sum-MSE (WSMSE) under per-BS power constraints and constraints on the number of streams per user. Moreover, although non-linear processing techniques such as vector precoding [9], [10] may generally be useful, we focus on more practical linear processing techniques. Both the SRM and WSMSE minimization (WSMMSE) problems are non-convex [11], and thus suboptimal design strategies of reasonable complexity are called for.
The contributions of this paper are as follows: (i) It is first shown in Sec. II-A that network MIMO with partial BS cooperation, that is, with partial message knowledge, is equivalent to a MIMO interference channel in which each transmitter knows the message of only one user under generalized linear constraints, which we refer to as MIMO-IFC-GC; (ii) We review the available suboptimal techniques that have been proposed for the SRM problem [12]- [14] and extend them to the MIMO-IFC-GC scenario where necessary in Sec. V. Since these techniques are generally unable to enforce constraints on the number of streams, we also review and generalize techniques that are based on the idea of interference alignment [5] and are able to impose such constraints; (iii) Then, we propose two novel suboptimal solutions for the WSMMSE problem in Sec. VI under arbitrary constraints on the number of streams. It is noted that the WSMMSE problem without such constraints would be trivial, as it would result in zero MMSE and no stream transmitted. The proposed solutions are based on a novel insight into the single-user MMSE problem with multiple linear constraints, which is discussed in Sec. IV; (iv) Finally, extensive numerical simulations are provided in Sec. VII to compare performance of the proposed schemes in realistic cellular systems.
Linear MMSE precoding and equalization techniques proposed in this paper were discussed briefly in [15]. The detailed analysis and discussion (including the proofs to the lemmas) are included in this paper. Additionally, we have also reviewed and extended available solutions to the SRM problem. Furthermore, we have included discussions of the complexity and overhead of the proposed techniques and previously available (and/or extended) solutions. Notation: We denote the positive definite matrices as A 0.
[·] + denotes max(·, 0). Capital bold letters represent matrices and small bold letters represent vectors. We denote the transpose operator with (·) T and conjugate transpose (Hermitian) with (·) H . A − 1 2 represents the inverse square of positive definite matrix A.

II. SYSTEM MODEL AND PRELIMINARIES
Consider the MIMO downlink system illustrated in Fig.  1 with M base stations (BSs) forming a set M, and K users forming a set K. Each BS is equipped with n t transmit antennas and each mobile user employs n r receive antennas. The mth BS is provided with the messages of its assigned users set K m ⊆ K. In other words, the kth user receives its message from a subset of M k BSs M k ⊆ M. Notice that, if K m contains one user for each transmitter m and M k = 1, then the model at hand reduces to a standard MIMO interference channel. Moreover, when all transmitters cooperate in transmitting to all the users, i.e., K m = K for all m ∈ M or equivalently M k = M , then we have a MIMO broadcast channel (BC).
We now detail the signal model for the channel at hand, which is referred to as MIMO interference channel with partial message sharing. Define as u k = [u k,1 . . . u k,d k ] T ∈ C d k the d k × 1 complex vector representing the d k ≤ min(M k n t , n r ) independent information streams intended for user k. We assume that u k ∼ CN (0, I). The data streams u k are known to all the BS in the set M k . In particular, if m ∈ M k , the mth BS precodes vector u k via a matrix B k,m ∈ C nt×d k , so that the signalx m ∈ C nt sent by the mth BS can be expressed as Imposing a per-BS power constraint, the following constraint must be then satisfied where P m is the power constraint of the mth BS.
The signal received at the kth user can be written as where H k,m ∈ C nr×nt is the channel matrix between the mth BS and kth user and n k is additive complex Gaussian noisẽ n k ∼ CN (0, I) 1 . We assume ideal channel state information at all nodes. In (3b), we have distinguished between the first term, which represents useful signal, the second term, which accounts for interference, and the noise.

A. Equivalence with MIMO-IFC-GC
We now show that the MIMO interference channel with partial message sharing and per-BS power constraints described above is equivalent to a specific MIMO interference channel with individual message knowledge and generalized linear constraints, which we refer to as MIMO-IFC-GC.
Definition 1: (MIMO-IFC-GC) The MIMO-IFC-GC consists of K transmitters and K receiver, where the kth transmitter has m t,k antennas and the kth receiver has m r,k antennas. The received signal at the kth receiver is where n k ∼ CN (0, I), the inputs are x k ∈ C m t,k and the channel matrix between the lth transmitter and the kth receiver is H k,l ∈ C m r,k ×m t,k . The data vector intended for user k is The precoding matrix for user k is defined as B k ∈ C m t,k ×d k so that x k = B k u k . The inputs x k have to satisfy M generalized linear constraints for given weight matrices Φ k,m ∈ C m t,k ×m t,k and m = 1, . . . , M. The weight matrices are such that matrices M m=1 Φ k,m are positive definite for all k = 1, . . . , K.
We remark that the positive definiteness of matrices M m=1 Φ k,m guarantees that the system is not allowed to transmit infinite power in any direction [16].
Lemma 1: Let (l) k be the lth BS in subset M k of BSs that know user k's message. The MIMO interference channel with partial message sharing (and per-transmitter power constraints) is equivalent to a MIMO-IFC-GC. This equivalent MIMO-IFC-GC is defined with m t,k = M k n t , m r,k = n r , channel matrices beamforming matrices and weight matrices Φ k,m being all zero except that their lth n t × n t submatrix on the main diagonal is I nt , if m = (l) k (If k / ∈ K m then Φ k,m = 0). We emphasize that the definition of MIMO-IFC-GC and this equivalence rely on the assumption of linear processing at the transmitters.
Proof: : The proof follows by inspection. Notice that matrices M m=1 Φ k,m are positive definite by construction. Given the generality of the MIMO-IFC-GC, which includes the scenario of interest of MIMO interference channel with partial message sharing as per the Lemma above, in the following we focus on the MIMO-IFC-GC as defined above and return to the cellular application in Sec. VII. It is noted that a model that subsumes the MIMO-IFC-GC has been studied in [16], as discussed below.

B. Linear Receivers and Mean Square Error
In this paper, we focus on the performance of the MIMO-IFC-GC under linear processing at the receivers. Therefore, the kth receiver estimates the intended vector u k using the receive processing (or equalization) matrix A k ∈ C d k ×m r,k asû The most common performance measures, such as weighted sum-rate or bit error rate, can be derived from the estimation error covariance matrix for each user k, which is referred to as Mean Square Error (MSE)-matrix (see [17] for a review). The name comes from the fact that that the jth term on the main diagonal of E k is the MSE on the estimation of the kth user's jth data stream u k,j . Using the definition of MIMO-IFC-GC, it is easy to see that the MSE-matrix can be written as a function of the equalization matrix A k and all the transmit matrices {B k } K k=1 as where Ω k is the covariance matrix that accounts for noise and interference at user k

III. PROBLEM DEFINITION AND PRELIMINARIES
In this paper, we consider the optimization of the sum of some specific functions f k (E k ) of the MSE-matrices E k of all users k = 1, . . . , K for the MIMO-IFC-GC. Specifically, we address the following constrained optimization problem where the optimization is over all transmit beamforming matrices B k and equalization matrices A k . Specifically, we focus on the weighted sum-MSE functions (WSMSE) with given diagonal weight matrices W k ∈ C d k ×d k where the main diagonal of W k is given by [w k,1 , ..., w k,d k ] with nonnegative weights w kj ≥ 0. With cost function (14), we refer to the problem (13) as the weighted sum-MSE minimization (WSMMSE) problem. Of more direct interest for communications systems is the maximization of the sum-rate. This is obtained from (13) by selecting the sum-rate (SR) functions With cost function (15), problem (13) is referred to as the sumrate maximization (SRM) problem. In fact, from informationtheoretic considerations, it can be seen that (15) is the maximum achievable rate (in bits per channel use) for the kth user where the signals of the other users are treated as noise (see, e.g., [17]). Remark 1: Consider an iterative algorithm where at each iteration a WSMMSE problem is solved with the weight matrices W k assumed to be non-diagonal and selected based on the previous MSE-matrix E k . This algorithm can approximate the solution of (13) for any general cost function f k (E k ). This was first pointed out in [18] for the weighted SRM problem in a MIMO BC, then in [19] for the single-antenna interference channel and a general utility function, and has been generalized to a MIMO (broadcast) interference channel in [20] with conventional power constraints. It is not difficult to see that this result extends also to the MIMO-IFC-GC, which is not subsumed in the model of [20] due to the generalized linear constraints. We explicitly state this conclusion below.
Lemma 2 [20]: For strictly concave utility functions f k (·) for all k, the global optimal solution of problem (13) and the solution of where g k (·) is the inverse function of the ∇f k (·), are the same.
Consequently, in order to find an approximate solution of (13), at each step matrices W k for k = 1. . . . , K are updated by solving (16) with respect to W k only (i.e., we keep (A k , B k ), ∀k unchanged in this step). Then, using the obtained matrices W k , for k = 1, . . . , K, the problem (16) reduces to a WSMMSE problem with respect to matrices A k and B k for k = 1, . . . , K (i.e., matrices W k are kept fixed). This results in the iterative algorithm, that is discussed in Remark 1 and that leads to a suboptimal solution of (13). In the special case of the SRM problem, we have The optimization problem (17) can be solved in an iterative fashion, where at each iteration the weights are selected as W ⋆ k = E −1 k and then the WSMMSE problem is solved with respect to matrices (A k , B k ) for k = 1, . . . , K.

IV. THE SINGLE-USER CASE (K = 1)
The WSMMSE and SRM problems are non-convex and thus global optimization is generally prohibitive. In this section, we address the case of a single user (K = 1). In particular, the SRM problem with K = 1 is non-convex if one includes constraints on the number of streams d 1 , but is otherwise convex and in this special case can be solved efficiently [17]. The global optimal solution for the single-user problem with multiple linear power constraint (and a rank constraint) is still unknown [21]. The WSMMSE problem is trivial without rank constraint, as explained above, and is non-convex. Here we first review a key result in [17][22] that shows with K = 1 and a single constraint (M = 1) the solution of the WSMMSE problem can be, however, found efficiently. We then discuss that with multiple constraints (M > 1), this is not the case, and a solution of the WSMMSE problem even with K = 1 must be found through some complex global optimization strategies. One such technique was recently proposed in [21] based on a sophisticated gradient approach. At the end of this section we then propose a computationally and conceptually simpler solution based on a novel result (Lemma 5), that our numerical result have shown to have excellent performance. This will be then leveraged in Sec. VI-B to propose a novel solution for the general multiuser case.
To elaborate, consider a scenario where the noise-plusinterference matrix Ω k (12) is fixed and given (i.e., not subject to optimization). Now, we solve the WSMMSE problem with K = 1 for specified weight matrices W and Φ m . For the rest of this section, we drop the index k = 1 from all quantities for simplicity of notation. We proceed by solving the problem at hand, first with respect to A for fixed B, and then with respect to B without loss of optimality. The first optimization, over A, is easily seen to be a convex problem (without constraints) whose solution is given by the minimum MSE equalization Plugging (18) in the MSE matrix (11). we obtain We now need to optimize over B the following problem Consider first the single-constraint problem, i.e., M = 1. The global optimal solution for single-user WSMMSE problem with M = 1 is given in [22][21] and reported below. Recall that, according to Definition 1, matrix Φ 1 is positive definite.
Lemma 3 [22]: The optimal solution of the WSMMSE problem with K = 1 and a single trace constraint (M = 1) is given by where U ∈ C mt×d is the matrix of eigenvectors of matrix Φ corresponding to its largest eigenvalues with µ ≥ 0 being the "waterfilling" level chosen so as to satisfy the single power constraint tr Φ 1 BB H = P 1 . Proof: Introducing the "effective" precoding matrixB = Φ 1/2 1 B and "effective" channel matrixH = HΦ − 1 2 1 , the problem is equivalent to the one discussed in [22,Theorem 1].
In the case of multiple constraints the approach used in Lemma 3 cannot be leveraged. Here we propose a simple, but effective, approach, which is based on the following considerations summarized in the following two lemmas.
Lemma 4: The precoding matrix (21)-(22) for a given fixed µ > 0 minimizes the Lagrangian function whereB is the effective precoding matrix defined above. Proof: We first note that (23) is the Lagrangian function of the single-user single-constraint problem solved in Lemma 2. Then, we prove (23) by contradiction. Assume that the minimum of the Lagrangian function is attained at where the corresponding E is not diagonal. Then, one can always find a unitary matrix Q ∈ C d×d such that the matrix . The function tr {WE} is Schur concave, and therefore the matrix B * does not decrease the function tr {WE} with respect toB, whilē BB H =B * B * H . This implies that the minimum of tr {WE} is reached when the MSE matrix is diagonalized. Therefore, we can set without loss of generalityB = UΣ where U is defined as in Lemma 3 and Σ is diagonal with non-negative elements on the main diagonal. Substituting this form ofB into the Lagrangian function, we obtain a convex problem in the diagonal elements of Σ, whose solution can be easily shown to be given by (22) for the given µ. This concludes the proof.
Lemma 5: Let p ⋆ be the optimal value of the single-user WSMMSE problem with multiple constraints (K = 1, M ≥ 1). We have where is the Lagrangian function of the single-user WSMMSE problem at hand and λ = (λ 1 , . . . , λ M ). Moreover, if there exists an optimal solutionB achieving p ⋆ that, together with a strictly positive Lagrange multiplierλ > 0, satisfies the conditions then (24) holds with equality. Proof: The proof is given in the appendix. Lemma 5 suggests that to solve the single-user multipleconstraint problem, under some technical conditions, one can minimize instead the dual problem on the right-hand side of (24). Lemma 3 showed that this is always possible with a single constraint. The conditions in Lemma 5 hold in most cases where the power constraints for the optimal solution are satisfied with equality. While this may not be always the case, in practice, e.g., if the power constraints represent per-BS power constraints as in the original formulation of Sec. II, this condition can be shown to hold [23].
Inspired by Lemma 5, here we propose an iterative approach to the solution of the WSMMSE problem with K = 1 that is based on solving the dual problem max λ≥0 min B L(B; λ). Specifically, in order to maximize inf B L(B; λ) over λ 0, in the proposed algorithm, the auxiliary variables λ is updated at the jth iteration via a subgradient update given by [16] so as to attempt to satisfy the power constraints. Having fixed the vector λ (j) , problem min B L(B, λ) reduces to minimizing (23) m Φ m and µ = 1. This can be done using Lemma 3, so that from (21)- (22), at the jth iteration, B (j) is obtained as The SRM problem for a number of users K > 1 is nonconvex even when removing the constraints on the number of streams per user. The general problem in fact remains non-convex and is NP-hard [24]. Therefore, since finding the global optimal has prohibitive complexity, one needs to resort to suboptimal solutions with reasonable complexity. In this section, we review several suboptimal solutions to the SRM problem that have been proposed in the literature. Since some of these techniques were originally proposed for a scenario that does not subsume the considered MIMO-IFC-GC, we also propose the necessary modifications required for application to the MIMO-IFC-GC. Note that these techniques perform an optimization over the transmit covariance matrices by relaxing the rank constraint due to the number of users per streams (see discussion below). Therefore, we also review and modify when necessary a different class of algorithms that solve problems related to SRM but are able to enforce constraints on the number of transmitting streams per user. The WSMMSE problem does not seem to have been addressed previously for the MIMO-IFC-GC and will be studied in the next section.

A. Soft Interference Nulling
A solution to the SRM problem for the MIMO-IFC-GC was proposed in [12]. In this technique the optimization is over all transmit covariance matrices The constraints on the number of streams would impose a rank constraint on Σ k as rank(Σ k ) = d k . Here, and in all the following reviewed techniques below, unless stated otherwise, such rank constraints are relaxed by assuming that the number of transmitting data streams is equal to the transmitting antennas to that user, i.e. d k = m t,k . From (15) and (18), we can rewrite the (negative) sum-rate as where Ω k is defined in (12). Notice that it is often convenient to work with the covariance matrices instead of the beamforming matrices B k , since this change of variables may render the optimization problem convex as, for instance, when minimizing the first term only in (29). It can then be seen that the SRM problem is, however, non-convex due to the presence of the − log |Ω k | term, which is indeed a concave function of the matrices Σ k . An approximate solution is then be found in [12] via an iterative scheme, whereby at each (j + 1)th iteration, given the previous solution Σ (j) k the non-convex term − log |Ω k | is approximated using a first-order Taylor expansion as where Ω is convex, a solution can be found efficiently. Following the original reference [12], we refer to this scheme as "soft interference nulling". We refer to [12] for further details about the algorithm.

B. SDP Relaxation
A related approach is taken in [13] for the SRM problem 2 for a MIMO-IFC with regular per-transmitter, rather than generalized, power constraints. Similarly to the previous technique, the optimization is over the transmit covariance matrices and under the relaxed rank constraints. In particular, the authors first approximate the problem by using the approach in [18]. Then, an iterative solution is proposed by linearizing a non-convex term similar to soft interference nulling as reviewed above. It turns out that such linearized problem can be solved using semi-definite programming (SDP). Specifically, denoting with Ω for the (j + 1)th iteration is given by and Y k is an auxiliary optimization variable, defined using the Schur complement as Y k = to convert the original optimization problem to an SDP problem [13]. The derivation requires minor modifications with respect to [13] and is 2 More generally, the reference studies the weighted SRM problem. therefore not detailed. The scheme is referred to as "SDP relaxation" in the following. We refer to [13] for further details about the algorithm.

C. Polite Waterfilling
Reference [16] studied the (weighted) SRM problem for a general model that includes the MIMO-IFC-GC. We review the approach here for completeness. Two algorithms are proposed, whose basic idea is to search iteratively for a solution of the KKT conditions [11] for the (weighted) SRM problem. Notice that, since the problem is non-convex, being a solution of the KKT conditions is only necessary (as proved in [16]) but not sufficient to guarantee global optimality. It is shown in [16] that any solution Σ k , k = 1, . . . , K, of the KKT conditions must have a specific structure that is referred to as "polite waterfilling", which is reviewed below for the SRM problem.
Lemma 6 [16]: For a given set of Lagrange multipliers λ = (µλ 1 , ..., µλ M ), where µ > 0 and λ i ≥ 0 for i = 1, ..., M , associated with the M power constraints in (13), define the covariance matriceŝ An optimal solution Σ k , k = 1, ..., K, of the SRM problem must have the "polite waterfilling" form where the columns of V k are the right singular vectors of the "pre-and post-whitened channel matrix" Ω with (12) for k = 1, . . . , K, and P k is a diagonal matrix with diagonal elements p k,i . The powers p k,i must satisfy where √ γ k,i is the ith singular value of the whitened matrix which implied by the constraints of the original problem (13). Moreover, parameters λ i ≥ 0 are to be chosen so as to satisfy each individual constraint in (13).
In order to obtain a solution Σ k , k = 1, . . . , K, according to polite waterfilling form as described in Lemma 6, [16] proposes to use the interpretation ofΩ k in (34) as the interference plus noise covariance matrix andΣ k in (35) as the transmit covariance matrix both at the "dual" system 3 .
Based on this observation, the algorithm proposed in [16] works as follows. At each jth iteration, first one calculates the covariance matrices Σ (j) k in the original system using the polite waterfilling solution of Lemma 6; then one calculates the matricesΣ (j) k using again polite waterfilling in the dual system as explained above. Finally, at the end of each jth iteration, one updates the Lagrange multipliers as thus forcing the solution to satisfy the constraints of the SRM problem (13). For details on the algorithm, we refer to [16]. Remark 2: Other notable algorithms designed to solve the SRM problem for the special case of a MIMO-BC with generalized constraints are [25], [26]. As explained in [16], these schemes are not easily generalized to the scenario at hand where the cost function is not convex. As such, they will not be further studied here.

D. Leakage Minimization
While the techniques discussed above do not enforce constraints on the number of stream per users, here we extend a technique previously proposed in [27] that aims at aligning interference through minimizing the interference leakage and is able to enforce the desired rank constraints. It is known that this approach is solves the SRM problem for high signal-tonoise-ratio (SNR). In this algorithm, it is assumed that the power budget is divided equally between the data streams and the precoding matrix of user k from BS m is given as B k,m = Pm Kmd kB k,m whereB k,m is a n t × d k matrix of orthonormal columns (i.e.B H k,mB k,m = I). The equalization matrices are also assumed to have orthonormal columns (i.e. A H k A k = I). Hence, there is no inter-stream interference for each user. Total interference leakage at user k is given by where To minimize the interference leakage, the equalization matrix A k for user k can be obtained as represents a matrix with columns given by the eigenvectors corresponding to the d k smallest eigenvalues of A. Now, for fixed matrices A k , the cost function (40) can be rewritten as 4 Minimizing over the matrices B k leads to choosingB k,m = v d k (Q k,m ). The algorithm iterates until convergence. We refer to this scheme as "min leakage" in the following.

E. Max-SINR
Another algorithm called "max-SINR" has been proposed in [27] which is based on the maximization of SINR, rather than directly the sum-rate. This algorithm is also able to enforce rank constraints. The max-SINR algorithm assumes equal power allocated to the data streams and attempts at maximizing the SINR for each stream by selecting the receive filters. Then, it exchanges the role of transmitter and receiver to obtain the transmit precoding matrices which maximizes the max-SINR. This iterates until convergence. A modification of this algorithm is given in [28] by maximizing the ratio of the average signal power to the interference plus noise power (SINR-like) term. However, these techniques are only given for standard MIMO interference channels and not for MIMO-IFC-GC.

VI. MSE MINIMIZATION
In this section, we propose two suboptimal techniques to solve the WSMMSE problem. We recall that with the WSMMSE problem enforcing the constraint on d k is necessary in order to avoid trivial solutions. Performance comparison among all the considered schemes will be provided in Sec. VII for a multi-cell system with network MIMO.

A. MMSE Interference Alignment
A technique referred to as MMSE interference alignment (MMSE-IA) was presented in [19] for an interference channel with per-transmitter power constraints and where each receiver is endowed with a single antenna. Here we extend the approach to to the MIMO-IFC-GC.
The idea is to approximate the solution of the WSMMSE problem by optimizing the precoding matrices B k followed by the equalization matrices A k and iterating the procedure. Specifically, initialize B k arbitrarily. Then, at each iteration j: (i) For each user k, evaluate the equalization matrices using the MMSE solution (18), obtaining A k , ∀k, this problem is convex in B k and can be solved by enforcing the KKT conditions. Therefore, matrices B where µ m are Lagrangian multipliers satisfying and the power constraints Once obtained the matrices B (j) k using the results in Lemma 7, the iterative procedure continues with the (j +1)th iteration. We refer to this scheme as extended MMSE-IA, or eMMSE-IA.
Remark 3: The algorithm proposed above reduces to the one introduced in [19] in the special case of per-transmitter power constraints and single-antenna receivers. It is noted that in such case, problem (42) can be solved in a distributed fashion, so that each transmitter k can calculate its matrix (more precisely vector, given the single antenna at the receivers) independently from the other transmitters. In the MIMO-IFC-GC, the power constraints couple the solutions of the different users and thus make a distributed approach infeasible.

B. Diagonalized MMSE
We now propose an iterative optimization strategy inspired by the single-user algorithm that we put forth in Sec. IV. At the (j + 1)th iteration, given the matrices obtained at the previous iteration, we proceed as follows. The weighted sum-MSE (14) with the definition of MSE-matrices (11) is a convex function in each A k and B k when (B j , A j ), ∀j = k are fixed. Nevertheless, it is not jointly convex in terms of both (A k , B k ). Inspired by Lemma 5 for the corresponding single-user problem, we propose a (suboptimal) solution based on the solution of the dual problem for calculation of (A k , B k ). To this end, we first obtain A k as (18). Then, we simplify the Lagrangian function with respect to B k by removing the terms independent of B k . Specifically, by defining Υ k = l =k H H l,k A l W l A H l H l,k , we have that the Lagrangian function at hand is given by This Lagrangian function for user k is the same as the Lagrangian function (25) of single-user WSMMSE problem when Φ(λ) is replaced with F k (λ) = Υ k + λ m Φ k,m . Matrix F k (λ) is non-singular and therefore, using the same argument as in the proof of Lemma 5, for a given Lagrange multipliers λ and given other users' transmission strategies (A l , B l ), ∀l = k, the optimal transmit precoding matrix can be obtained as corresponding to its largest eigenvalues γ k,1 ≥ . . . ≥ γ k,d k and Σ k is diagonal matrices with the elements √ p k,i given by with λ 0 being the Lagrangian multipliers satisfy the power constraints. Since this scheme diagonalizes the MSE matrices defined in (9), it is referred to as diagonalized MMSE (DMMSE).
To summarize, the proposed algorithm at each iteration j (i) evaluates the transmit precoding matrices B to satisfy the power constraints. Remark 4: In this paper, we assume perfect knowledge of channel state information (CSI). Therefore, each transmitter and receiver has sufficient information to calculate the resulting precoders and equalizers by running the proposed algorithms. Under this assumption, which is common to other reviewed works such as [12] [13], no exchange of precoder and equalizer vectors is required between the transmitters and receivers. Nevertheless, in practice, the CSI may only be available locally, in the sense that transmitter k knows channel matrices H l,k , for all l = 1, . . . , K, whereas receiver k is aware of channel matrices H k,l , for all l = 1, . . . , K. The proposed DMMSE and the reviewed PWF [14][16] algorithms require, beside the local CSI, that the transmitter k has available also the interference plus noise covariance matrix, Ω k , and the current equalization matrices A l for all l = 1, . . . , K in order to update the precoder for user k. Hence, to enable DMMSE and PWF with local CSI, exchange of the equalizer matrices is needed between the nodes. Similarly, the proposed eMMSEIA, and min leakage and Max-SINR algorithms [27], require the transmitters to know the equalizing matrices A l for l = 1, . . . , K at each iteration, in addition to the local CSI. Moreover, each receiver must know the current precoders B l for all l = 1, . . . , K. Therefore, the overhead for the proposed eMMSEIA and the min leakage and Max-SINR algorithms involves the exchange of equalizer and precoder matrices between the transmitters and receivers. However, these latter algorithms can also be adapted using the bidirectional optimization process proposed in [29]. This process involves bi-directional training followed by data transmission. In the forward direction, the training sequences are sent using the current precoders. Then, at the user receivers the equalizers are updated to minimize the least square error cost function. In the backward training phase, the current equalizers are used to send the training sequences and the precoders are updated accordingly. Finally, the SIN [12] and SDP relaxation [13] techniques are applied in a centralized fashion (rather than by updating the transmitter and receiver for each user at each iteration), and they require centralized full knowledge of all channel matrices.
Remark 5: Reference [13] addresses the SRM problem for a MIMO-IFC with regular per-transmitter, rather than generalized, power constraints. The problem is addressed by solving an SDP problem at each iteration. Moreover, the optimization is over the transmit covariance matrices and under the relaxed rank constraint. This enforces a constraint on the number of transmitted streams per user. References [14]- [16] study the (weighted) SRM problem by decomposing the multiuser problem into single-user problems for each user. Each singleuser problem is a standard single-user SRM problem with an additional interference power constraint. The approach used in [14]- [16] assumes that the number of transmitted streams is equal to n r . In this paper, we address WSMMSE problem and allow for an arbitrary number of streams (d k ≤ n r ).
Remark 6: Our algorithms consists of an inner loop, which solves the WSMMSE problem, and an outer loop, which is the subgradient algorithm to update λ. The subgradient algorithm in the outer loop is convergent (with a proper selection of the step sizes [30]) due to the fact that the dual function inf B L(B; λ) is a concave function with respect to λ [11]. The inner loops of the proposed algorithms in this paper (i.e. eMMSEIA and DMMSE) are convergent since the objective function decreases at each iteration. A discussion of the convergence for a special case of the eMMSEIA algorithm can be found in [19]. However, the original problem is non-convex and our solutions are only local minima. Nevertheless, the DMMSE algorithm is shown to converge to a local minimum with better performance compared to the previously known schemes in Sec. VII.

VII. NUMERICAL RESULTS
We consider a hexagonal cellular system where each BS is equipped with n t transmit antennas and each user has n r receive antennas. The users are located uniformly at random. Two tiers of surrounding cells are considered as interference for each cluster. We consider the worst-case scenario for the inter-cluster interference, which will be the condition that interfering BSs transmit at the full allowed power [7], [8], [31], [32]. We define the cooperation factor κ as a number of BSs cooperating on transmission to each user. The κ BSs are assigned to each user so that the corresponding channel norms (or, alternatively, the corresponding received SNRs) are the largest.
The propagation channel between each BS's transmit antennas and mobile user's receive antenna is characterized by path loss, shadowing and Rayleigh fading. The path loss component is proportional to d −β km , where d km denotes distance from base station m to mobile user k and β = 3.8 is the path loss exponent. The channel from the transmit antenna t of the base station b at the receive antenna r of the kth user is given by [7] H (r,t) where α radius. γ 0 is the interference-free SNR at the cell boundary. We consider one user randomly located per cell for the numerical results.
When sectorization is employed, the transmit antennas are equally divided among the sectors of a cell. Each transmit antenna has a parabolic beam pattern as a function of the direction of the user from the broadside direction of the antenna (For more details refer to [7], [33]). The antenna gain is a function of the direction of the user k from the broadside direction of the tth transmit antenna of the bth base station denoted by Θ (t) k,b ∈ [−π, π]; Θ 3dB is the half-power angle and A s is the sidelobe gain. The antenna gain is given as [33] A Θ For the 3,6-sector cells A s = 20, 23 dB and Θ 3dB = 70π 180 , 35π 180 , respectively [7], [33], [34]. When there is no sectorization we set A = 1.
We first compare different algorithms (for the solution of the SRM problem) without enforcing rank constraints on SIN, PWF, SDP relaxation and setting d k = min(m t,k , m r,k ) = n r for the eMMSEIA and DMMSE algorithms. To solve the SRM problem, the weight matrices in the eMMSEIA and DMMSE algorithms are updated at each iteration as W k = E −1 k using the current MSE-matrix E k . Fig. 2 compares the per-cell sumrate of the algorithms discussed in this paper for a cluster with M = 3 cells and a cooperation factor κ = 2. The results show that our proposed DMMSE algorithm outperforms other techniques, while the polite water-filling algorithm (PWF) [14], [16] has a similar performance. Our proposed eMMSEIA scheme converges to a poorer local optimum value compared to these two schemes. The soft interference nulling (SIN) [12] and SDP relaxation [13] algorithms, which use the approximation of the non-convex terms in the objective function, perform worse in this example.
In Fig. 3, we evaluate the effect of partial cooperation for the DMMSE, eMMSEIA, and PWF algorithms in a cluster of size M = 5 where each BS is equipped with n t = 4 transmit antennas, each user employs n r = 2 receive antennas, and 2 users are dropped randomly in each cell. Recall that the cooperation factor κ represents the number of BSs cooperating in transmission to each user. It can be seen that as κ increases the performance improves with diminishing returns as κ grows large. Moreover, the relative performance of the algorithms confirms the considerations above.
In Fig. 4, we compare again the performance of the schemes considered in Fig. 3 but with a stricter requirement on the number of streams, namely d k = 1. It can be seen that the proposed DMMSE tends to perform better than PWF, which was not designed to handle rank constraints. We have adopted the PWF algorithm to support d k < min(m t,k , m r,k ) by using a thin SVD ofΩ  Recall that M represents the set of BSs whose transmission is coordinated, but only κ BSs cooperate for transmission to a given user. As an example, for a cluster size of M = 7 a cooperation factor of κ = 4 performs almost as well as the full cooperation scenario with κ = 7. Moreover, the performance gains with respect to the non-cooperative case κ = 1 are evident. We also show the performance with a cluster containing a single cell, i.e., M = 1. This highlights the performance gains attained even in the absence of message sharing among the BSs (i.e., κ = 1) due to the coordination of the BSs within the cluster.
Finally, the effect of sectorization is studied in Fig. 6 where n t = 6 transmit antennas at each BS are divided equally into S = 1, 3, 6 sectors. Each cell contains 6 users, each equipped with n r = 2 receive antennas. The users are randomly located at the distance of 2 3 d 0 from its BS. For a given channel realization the DMMSE algorithm is used to obtain the percell sum rate. The cumulative distribution functions (CDFs) of per-cell sum rates are computed using large number of channel realizations. The gains of sectorization and cooperation are compared. For example, the system with coordination of 7 cells and κ = 3 cooperation factor and without sectorization performs better than the sectorized system with S = 6 and without any coordination between the BSs.

VIII. CONCLUSIONS
In this paper, we have studied a MIMO interference channel with partial cooperation at the BSs and per-BS power constraints. We have shown that the channel at hand is equivalent to a MIMO interference channel under generalized linear constraints (MIMO-IFC-GC). Focusing on linear transmission strategies, we have reviewed some of the available techniques for the maximization of the sum-rate and extended them to the MIMO-IFC-GC when necessary. Moreover, we have proposed two novel strategies for minimization of the weighted mean square error on the data estimates. Specifically, we have proposed an extension of the recently introduced MMSE interference alignment strategy and a novel strategy termed diagonalized MSE-matrix (DMMSE). Our proposed strategies support transmission of any arbitrary number of data streams per user. Extensive numerical results show that the DMMSE outperforms most previously proposed techniques and performs just as well as the best known strategy. Moreover, our results bring insight into the advantages of partial cooperation and sectorization and the impact of the size of the cooperating cluster of BSs and sectorization.
We conclude with a brief discussion on the complexity of the algorithms. Due to the difficulty of complete complexity analysis, especially in terms of speed of convergence, we present a discussion based on our simulation experiments. The PWF algorithm converges in almost the same number of iterations as the DMMSE algorithm. The complexity per iteration of PWF and DMMSE is also almost the same as K O(κn t n 2 r ) + O(n 3 r ) (required for the thin SVD operation). However, the PWF algorithm contains additional operations (matrix inversion and SVD) to obtain the precoding matrices from the calculated transmit covariance matrices. 5 Also, the PWF algorithm includes a water-filling algorithm within its inner loop, which is not required in the DMMSE algorithm. The eMMSEIA algorithm has lower complexity per iteration (i.e. KO(n 3 r )) than the PWF and DMMSE algorithms, since its complexity is due to a matrix inversion per iteration per user. However, eMMSEIA converges in a larger number of iterations than DMMSE and PWF. The complexity per iteration for the SDP relaxation is higher than for the SIN algorithm (this is because of the extra auxiliary positive semidefinite matrix variable, Y, introduced in the SDP relaxation algorithm). The SIN algorithm also converges in a smaller number of iterations than the SDP relaxation algorithm. 5 This can be performed together with finding the MMSE receive matrices. APPENDIX PROOF OF LEMMA 5 The inequality (24) follows from weak Lagrangian duality. We now prove the second part of the statement. Recognizing now that tr {WE} with (19) is a Schur-concave function of the diagonal elements of (19) 6 , it can be argued that the minimum is attained when E is diagonalized as we did for Lemma 4. Defining R = H H Ω −1 H, we can conclude that B H RB must be also diagonal in this search domain. Now assume that an optimal solution of the single-user WSMMSE problem is denoted asB. Without loss of generality we can assume that this solution diagonalizes the MSE matrices. The necessity of the KKT conditions can be proved as in [16] and in special cases such as the MIMO interference channel with partial message sharing of Sec. II, it also follows from linear independence constraint qualification conditions [30].
Hence, there exists a Lagrange multiplier vectorλ which together withB satisfies the KKT conditions of the WSMMSE problem (20) [18] [30]. As it is stated in the Lemma, we consider the case thatλ m are also strictly positive (i.e.λ m > 0 for all m). Simplifying the KKT condition (26), we have 7 Left-multiplying (52) byB H gives us SinceB H RB and correspondinglyẼ are diagonal matrices, B H mλ m Φ m B must also be diagonal. For simplicity, we introduce Φ(λ) = M m=1λ m Φ m . Sinceλ m > 0 for every m, therefore Φ(λ) is a non-singular matrix. This can be easily verified due to the structure of Φ m . Hence, we can writeB H Φ(λ)B =∆ where∆ ∈ C d×d is a diagonal matrix. Therefore, we can write whereŨ ∈ C mt×d consists of orthonormal columns (i.e. U HŨ ) andΣ ∈ C d×d is a diagonal matrix with the diagonal terms of √p i . Hence, we can writẽ Replacing the structure ofB given in (55), we can writẽ Thus, we can conclude from the equation above thatŨ must contain the eigenvectors of Φ(λ) − 1 2 RΦ(λ) − 1 2 .
6 A Schur-concave function f (x) of vector x = (x 1 , ..., x d ) is such that  7 We use differentiation rule ∇ X tr AX H B = BA and ∇ X tr Y −1 = −Y −1 (∇ X Y) Y −1 . For the complex gradient operator each matrix and its conjugate transpose are treated as independent variables [35]. Now, plugging (55) into (26) and left-multiply it with Φ − 1 2 , we getΓΣ whereΓ(λ) = diag[γ 1 (λ) · · · γ d (λ)] is a diagonal matrix with the diagonal terms of the d largest eigenvalues of Φ(λ) − 1 2 RΦ(λ) − 1 2 . Since all the matrices are diagonal, (57) reduces to the scalar equations: Solving these equations gives us the optimalp i given bỹ Thus, for the given Lagrange multiplierλ which together with B, satisfying the KKT conditions of (20),B must satisfy (55) and (59). If all power constraints are satisfied with equality by this solution, then (55) and (59) The solution of this problem is given in Lemma 3 as where U consists of d eigenvectors of Φ(λ) − 1 2 RΦ(λ) − 1 2 corresponding to its largest eigenvalues and Σ is a diagonal matrix with the diagonal elements of √ p i , which is given by If we set µ = 1 and comparing (59) and (62), we can conclude thatp i = p i , ∀i which together with comparison of (61) and (55) we can conclude that B(λ) =B and the µ = 1 is the optimal Lagrange multiplier of the single-constraint WSMMSE problem (60). Following Lemma 4, this precoding matrix is also a result of minimization of the Lagrangian function (23) thus concluding the proof.


extfieldsofstudy:
Mathematics
Computer Science

text:
Output sensitive algorithm for covering many points

A set of points and a positive integer $m$ are given and our goal is to cover the maximum number of these point with $m$ disks. We devise the first output sensitive algorithm for this problem. We introduce a parameter $\rho$ as the maximum number of points that one disk can cover. In this paper first we solve the problem for $m=2$ in $O({n\rho} + {\rho ^3}\log \rho ))$ time. The previous algorithm for this problem runs in $O({n^3}\log n)$ time. Our algorithm outperforms the previous algorithm because $\rho$ is much smaller than $n$ in many cases. Then we extend the algorithm for any value of $m$ and we solve the problem in $O(m{n\rho} + {(m\rho )^{2m - 1}}\log m\rho )$ time. The previous algorithm for this problem runs in $O({n^{2m - 1}}\log n)$ time. Our algorithm runs faster than the previous algorithm because $m\rho$ is smaller than $n$ in many cases. Our technique to obtain an output sensitive algorithm is to use a greedy algorithm to confine the areas that we should search to obtain the result. Our technique in this paper may be applicable in other set covering problems that deploy a greedy algorithm, to obtain faster solutions.

Introduction
In the classic covering problem, a set of points are given, and the goal is to place the minimum number of unit disks to cover all the input points. In the maximum coverage problem, the number of disks to be used, m, is given and the goal is to cover the maximum number of points with m disks. We call this problem M ostP oints(P, m). This problem was introduced by Drezner [1] for m = 1 and it was solved in O(n 2 log n) time. To solve the problem, he replaced every point with a disk centered at that point and he obtained the maximum depth in the arrangement of disks. Later, Chazelle and Lee [2] solved the problem in O(n 2 ) time. This problem belongs to the 3-SUM hard complexity class [3], which means that, the O(n 2 ) running time algorithm is the best one for this problem. However, there are some algorithms which approximately solve this problem in less than O(n 2 ) time. These algorithms approximate the disk radius and the number of covered points. Figueiredo and Fonesca [5] gave an algorithm to cover the maximum number of points that a unit disk in d can cover, with a (1 + ε)-radius disk in O(n/ε d−1 ) time. In [4] Aronov and Har-peled gave a (1 − ε)-approximation algorithm in the number of covered points which runs in O(nε −2 log n) time. M ostP oints(P, m), for m > 1, is NP-Hard unless we consider m as a constant. A trivial greedy algorithm is a (1 − 1 e )-approximation algorithm for it [6]. The greedy algorithm first finds a disk which covers the maximum number of points in O(n 2 ) time. To pick the next disks, it removes the points located in the first disk, and finds the disk which covers the maximum total weight of points as the second disk. It repeats this process until m disks are picked up. This yields a (1 − 1 e )-approximation algorithm which runs in O(mn 2 ) time. The first (1 − ε)approximation algorithm for this problem was given in [7] which runs in O(n log n + nε −6m+6 log( 1 ε )) time. We presented a polynomial time approximation schema for this problem in our previous work [8], The trivial method to obtain the optimal result of the M ostP oints(P, m) is to consider all subsets of size m of n 2 possible disks and finds the subset which covers the most number of points. This method takes O(n 2m ) time. In [7] de Berg et. al gave an algorithm for this problem. They first solved the problem for m = 2 in O(n 3 log n) time. To solve the problem for m > 2, they fixed every subset of m − 2 disks, and they found the best 2 disks after removing the points contained in the m − 2 disks. This takes O(n 2m−1 log n) time.
In section 2, we present an output sensitive algorithm to obtain the optimal result of M ostP oints(P, 2) which runs in O(nρ + ρ 3 log ρ) time where ρ is the maximum number of points that one disk can cover. In Section 3 we extends the algorithm for M ostP oints(P, m) . In section 4 we compare the implementation results of our algorithm and the algorithm of [7] and in section 5 we conclude the paper. Our innovation in this algorithm is to use the greedy algorithm to find the regions in which the m resultant disks may reside, and searching in those regions only. In the rest of the paper when we use disks we mean unit disks.

Output sensitive algorithm for MostPoints(P,2)
In this section we describe our algorithm for M ostP oints(P, 2). Let g 1 = M ostP oints(P, 1) be the disk that covers the maximum number points from the point set P . Let P oints(g 1 ) refers to the points of P which are located inside g 1 . Define g 2 = M ostP oints(P − P oints(g 1 ), 1) as the disk that covers the maximum number of points after removing the points located in g 1 . g 1 and g 2 are the result of the greedy algorithm for M ostP oints(P, 2). An example is given in Figure 1. In this example, the optimal solution covers 18 points whereas the greedy algorithm covers 15 points. As illustrated in this example, the two disks of the optimal solution have common points with g 1 . Intuitively, either the disks of the optimal solution have common points with g 1 or, the greedy algorithm obtains the optimal solution. Lemma 1 proves this claim. By the best two disks we mean the optimal solution of M ostP oints(P, 2). Let D be a set of disks. The function Cover(D) denotes the number of points covered by the disks in D.
Lemma 1: Let g 1 be one of the disks that covers the maximum number of points and g 2 be one of the disks that covers the maximum number of points after removing the points located in g 1 . Let o 1 and o 2 be the two disks that together cover the maximum number of points among any combination of two disks: Either both o 1 and o 2 have common points with g 1 , or g 1 and g 2 cover the maximum number of points.
Proof: Suppose that at least one of o 1 and o 2 , say o 2 , does not have any common point with g 1 . Then we have Cover({g 1 , g 2 }) ≥ Cover({g 1 , o 2 }), because g 2 covers the maximum number of points not In this example, the greedy algorithm returns g1 and g2 which cover 15 points together, whereas the disks in the optimal solution, the two disks specified with the dash lines, cover 18 points together. covered by g 1 . We also have Cover({g 1 , o 2 }) ≥ Cover({o 1 , o 2 }) because g 1 covers at least the same weight of points as o 1 (since it is a disk that covers the maximum weight of points) and g 1 does not have any point in common with o 2 . Thus Cover({g 1 , g 2 }) ≥ Cover({o 1 , o 2 }), which implies that the total number of the points covered by g 1 and g 2 is maximal. 2 Based on Lemma 1, we should look for the best two disks in the disks which have common points with g 1 . These disks are located around g 1 in a circle having the same center as g 1 and radius 3, where the radius of g 1 is 1. We call this region N g 1 . Figure 2 shows N g 1 . Lemma 2 bounds the number of points in N g 1 .
Lemma 2: Let ρ be the maximum number of points that can be covered by a disk. The number of points in N g 1 is O(ρ).
Proof: We can pack N g 1 with at most 21 unit disks, which is shown in Figure 3. Every one of these 21 disks covers equal or less than ρ points. So the number of points in N g 1 is at most 21ρ points, which is O(ρ) points. 2 We can also devise an output sensitive algorithm for M ostP oints(P, 1) that solves the problem in O(nρ) time. The algorithm is as follows. Cover the plane using four shifted grids G 1 , ..., G 4 . Each grid covers the whole plane, but the grids are shifted relative to each other, 2 units from each side. The cells of the grids have size 4x4. For any unit disk d there is a grid G i such that d is fully contained in a single cell of G i . Now, for each non-empty grid cell C of every grid, compute the optimal disk on the set of points lying inside C. C contains at most cρ points, for a constant value of c and the time spent for cell C is O(|C| 2 ) where |C| denotes the number of points in the cell. The total number of points over all cells C is 4n, since each point is contained in exactly one cell per grid. To obtain the running time of the algorithm, we count the cells of the grids C 1 , ..., C l . The total running time of the algorithm is: Using Lemma 1 and Lemma 2, the best two disks can be g 1 and g 2 , or the two disks around g 1 that cover the maximum number of points. To find the best two disks in P , first we find g 1 and g 2 in O(nρ) time. Then, we determine the points in N g 1 , and we run the algorithm of [7] for the points located in N g 1 . This algorithm obtains the best two disks in N g 1 , in O(ρ 3 log ρ) time. Finally, we compare this result, the best two disks in N g 1 , with g 1 and g 2 and consider the one that covers more weight as the optimal solution for M ostP oints(P, 2). Since ρ can be much smaller than n in practice, our algorithm which runs in O(nρ + ρ 3 log ρ) time, is an improvement over the algorithm of [7] that requires O(n 3 log n) time.

The algorithm for MostPoints(P,m)
The algorithm can be extended to cover the maximum number of points with m disks. From now on, by the best j disks, we mean a set of j disks that cover the maximum total weight of points. The proposed algorithm is an iterative algorithm. At the first iteration, the algorithm finds two disks, which cover the maximum number of points as explained before. To find the best 3 disks, the algorithm uses the following fact, which will be proven in Lemma 3: the best 3 disks are among the disks that have some common points with the best 2 disks; Otherwise, the best 2 disks plus one disk that covers the maximum number of points, after removing the points located in the best two disks, are the best 3 disks. The algorithm repeats this process until the best m disks are found. Now we explain the algorithm formally. Remember that the function M ostP oints(Q, j) returns a set of j disks that cover the maximum number of points in the point set Q. Let D be a set of unit disks. Define P oints(D) = {p|p ∈ (P ∩ D)} as the subset of P that is located in the disks of D. Define N eighbor(D) as the points in the circles with the same centers as the disks in D and radius 3. Let OP T i denotes the best i disks. The outline of the algorithm is as follows: The exact output sensitive algorithm for the most points covering problem Inputs: A set P of n points in the plane, and a positive integer m. Outputs: A set of m disks which cover the maximum number of points from P. Definitions: The function MostPoints(Q,j), Points(D), and Neighbor(D) as defined in the text.

end
In Lemma 3 we prove the correctness of the above algorithm. Before proving Lemma 3 we have the following definitions. Proof: We prove the correctness of the algorithm by induction on the number of iterations. Assume that at the end of iteration i we have correctly computed OP T i . Consider OP T i+1 , which consists of i + 1 disks. These disks may have common points with P oints(OP T i ) or not. We consider both cases. In case 1 we consider that all i + 1 disks of OP T i+1 have common points with P oints(OP T i ), and in case 2 we consider that some of these i disks do not have any common point with P oints(OP T i ).
Case1: All i + 1 disks of OP T i+1 have common points with P oints(OP T i ): All disks that have common points with P oints(OP T i ), have some points in N eighbor(OP T i ). The algorithm, in line 8, obtains the best i + 1 disks in the points of N eighbor(OP T i ). So in this case the algorithm exactly computes OP T i+1 .
Case 2: Some of the disks in OP T i+1 do not have any common point with P oints(OP T i ): Call one of these disks F . We can partition OP T i+1 to {OP T i+1 − F } and F . Based on the Definition 1, we have: As OP T i is a set of i disks which cover the maximum number of points, we have: Furthermore, g i+1 is the disk that covers the maximum number of points after removing the points in P oints(OP T i ). Therefore, g i+1 covers the maximum number of points among the disks that have not common points with OP T i . Thus, we have: From (2) and (3) we have: From (1) and (4) we have: Thus in case 2, OP T i and g i+1 cover the maximum number of points. The algorithm compares the results of case 1 with the result of case 2, and considers the best one as OP T i+1 . So, the algorithm correctly computes OP T i+1 . The base case, i = 1, is correct and can be satisfied using the algorithm devised for M ostP oints(P, 1). 2 In the following lemma, we obtain the running time of the algorithm. Lemma 4: The exact output sensitive algorithm for the most points covering problem runs in O(mnρ+ (mρ) 2m−1 log mρ) time, where ρ is the maximum number of points that one disk can cover in the points set P .

Proof:
The algorithm runs in m iterations. At the start of the iteration i, we have found the best i − 1 disks, and we are getting to find the best i disks. To find the best i disks, we obtain g i in line 7, which takes O(nρ) time. We also look for the i disks covering the maximum number of points in N eighbor(OP T i−1 ), in line 8. N eighbor(OP T i−1 ) consists of i − 1 disks, and based on Lemma 2, the neighborhood of each disk can be packed with at most 21 disks. So the maximum number of points in N eighbor(OP T i−1 ) is 21ρ(i−1), which is O(iρ) points. We search for the i disks that cover the maximum number of points in O(iρ) points. The algorithm of [7] is applied for these points, which takes O((iρ) 2i−1 log(iρ)) time. The algorithm runs in m iterations so, the running time is bounded by O(nρ) + m i=2 (O(nρ) + (iρ) 2i−1 log(mρ)) which is O(mnρ + (mρ) 2m−1 log(mρ)). 2

Implementation results
We have implemented our algorithm and the algorithm of [7] for M ostP oints(P, 2). In each algorithm, we observed the number of pairs of disks that each algorithms consider to obtain the optimal solution. We distributed different number of points randomly in squares with different side length and we run both algorithms on the points. Table 1 shows the result of the comparison. The results show that our algorithm computes significantly less number of pairs than the algorithm of [7] in most cases. When ρ n becomes close to 1, the number of pair of disks that are considered in both algorithms become close to each other.

Extending the algorithm for other shapes and weighted point sets
This algorithm can be applied for other types of maximum coverage problem such as covering point sets with other shapes and maximum coverage problem in sets. The algorithm uses the parameter ρ which is Tab. 1: Number of pairs of disks processed in our algorithm an algorithm of [7] Num of points square side ρ Num of pairs in algorithm of [7] Num of pairs in our algorithm  1000  200  19  7007  201  1000  100  44  22878  1131  1000  50  125 45414  4259  500  50  105 26041  6402  500  100  34  8605  252 the cardinality of the largest subset. Lemma 1 and 3 can be applied for general set systems. So we should search in the subsets that have common members with the largest subsets. To bound the running time of the algorithm there should be a constant number of subsets that have common points with the largest subsets.

Conclusion and future works
In this paper we presented an exact output sensitive algorithm for covering many points problem. In the algorithm, we confined the region in which the resultant disks may reside using the greedy algorithm, and we searched in that regions only. This led to an output sensitive algorithm for this problem which runs in O(mnρ + (mρ) 2m−1 log mρ) time, where ρ is the maximum number of points that one disk can cover. Our algorithm can be improved by improving the running time of the algorithm in [7]. If better algorithms for covering most points with two disks are suggested, they can improve the running time of our algorithm and the algorithm of [7]. Furthermore, considering ρ as a parameter, and using our technique, it may be possible to improve the analysis of other covering algorithms.