A Semideﬁnite Relaxation Method for Elliptical Location

: Wireless location is a supporting technology in many application scenarios of wireless communication systems. Recently, an increasing number of studies have been conducted on range-based elliptical location in a variety of backgrounds. Speciﬁcally, the design and implementation of position estimators are of great signiﬁcance. The difﬁculties arising from implementing a maximum likelihood estimator for elliptical location come from the nonconvexity of the negative log-likelihood functions. The need for computational efﬁciency further enhances the difﬁculties. Traditional algorithms suffer from the problems of high computational cost and low initialization justiﬁability. On the other hand, existing closed-form solutions are sensitive to the measurement noise levels. We recognize that the root of these drawbacks lies in an oversimpliﬁed linear approximation of the nonconvex model, and accordingly design a maximum likelihood estimator through semideﬁnite relaxation for elliptical location. We relax the elliptical location problems to semideﬁnite programs, which can be solved efﬁciently with interior-point methods. Additionally, we theoretically analyze the complexity of the proposed algorithm. Finally, we design and carry out a series of simulation experiments, showing that the proposed algorithm outperforms several widely used closed-form solutions at a wide range of noise levels. Extensive results under extreme noise conditions verify the deployability of the algorithm.


Introduction
During the deployment of powerful location-based applications such as enhanced 911 services, asset management, and workflow automation, wireless location is a classical problem and has attracted intensive research recently [1][2][3][4][5][6]. Technically, wireless location can be based on a variety of different kinds of measurements, such as range [7], angle [8], energy [9], visible light [10,11], or fingerprinting [12][13][14]. The common range-based measurements include time of arrivals (TOAs) [15,16] and time difference of arrivals (TDOAs) [17,18]. TDOAs (or TOAs) are combined to produce intersecting hyperbolic (or circular) lines from which the target location is determined. In addition to these two measurement schemes, there is a less well-known range-based approach called elliptical location [19][20][21][22][23][24][25][26][27]. In particular, in a wireless location system with multiple transmitters and multiple receivers, the sum of each pair of transmitter-target range and target-receiver range defines an ellipse, and the target resides on the intersection of all these ellipses.
Depending on whether the transmitters and receivers are synchronized in time, we have the synchronous elliptical location (SEL) [20,23,24] and the asynchronous elliptical location 1.
We present the negative log-likelihood (NLL) function of the elliptical location estimation problem with multiple transmitters and multiple receivers in a concise form, and then derive the Cramér-Rao bound (CRB) of the problem. This is the theoretical basis for the feasibility analysis of the problem.

2.
We identify the NLL function minimization as a quadratic optimization problem with multiple quadratic equality/inequality constraints and design its semidefinite relaxation. 3.
We find the dual problem of the primal SDP problem and give the pseudo-code of the interior-point primal-dual path-following algorithm. With these results, engineers can design customized algorithms themselves, instead of relying on codes provided by mathematicians that are not optimized for a specific problem. 4.
We theoretically analyze the polynomial time complexity of the algorithm and experimentally confirm that the proposed estimator attains the CRB approximately. 5.
We cautiously verify the deployability of the algorithm by setting the measurement noise levels far higher than what can be achieved by real systems. Simulation experiments with extreme parameter settings are alternative to the costly field tests.
The rest of this paper is organized as follows. Section 2 provides the location scenario and its observation model. Section 3 performs the CRB analysis for the elliptical location problem. Section 4 presents the SDP-based location algorithm and its complexity analysis. Section 5 contains the simulation results and compares our algorithm with state-of-the-art ones, and Section 6 draws conclusions.
We shall use bold uppercase letter to denote a matrix and bold lowercase letter to denote a vector. The notation M 0 means that M is a positive semidefinite matrix, and v ≥ 0 means that each element of v is non-negative. 1 p×q is a p × q ones matrix, and I p is an identity matrix of size p. The symbol ⊗ is the Kronecker product operator. v is the Euclidean norm of v. cond(M) is the condition number of M. diag(v) is the square diagonal matrix with the elements of vector v on the main diagonal. B = blkdiag(M 1 , . . . , M N ) is the block diagonal matrix created by aligning the matrices M 1 , . . . , M N along the diagonal of B. If X and Y are symmetric, then X, Y is a Frobenius inner product of X and Y. We use ρ x,y as a shorthand for (x − y)/ x − y .

Problem Formulation and Data Model
We are interested in a d-dimensional (d = 2 or d = 3) elliptical location scenario where the transmitters and receivers are stationary. Let M be the number of transmitters and N be the number of receivers. The position vectors of transmitters and receivers are t i and s j , for i = 1, 2, . . . , M and j = 1, 2, . . . , N, which are available in practice. Let c be the known signal propagation speed. The unknown position vector of the target of concern is denoted by u o . A typical location scenario for a two-dimensional case is illustrated in Figure 1. The equations of the ellipses in Figure 1 are shown in Equation (1). The target is located at the intersection of these ellipses.
The transmitter at t i radiates a signal and the receiver at s j observes the signal from direct propagation (from t i to s j ) and from indirect reflection of the target at u o (from t i to u o and then from u o to s j ). The ideal differential TOA measurement between the indirect path and direct path is Obviously, the graph of u o is an ellipse with foci t i and s j for each pair of i and j.
In an actual elliptical location system, we have the measurements where n i,j is the measurement noise. For notation simplicity, we collect τ i,j to form a measurement vector. For the transmitter at t i (i = 1, 2, . . . , M), all the related measurements can be collected in a column vector τ i = [τ i,1 , τ i,2 , . . . , τ i,N ] T . Further, the measurements related to all the transmitters can be denoted as τ = [τ T 1 , τ T 2 , · · · , τ T M ] T . It is customary to assume that where τ o i is the mean vector of τ i and Q τ i is a known covariance matrix of τ i [28]. Furthermore,

Cramér-Rao Bound
We now give the CRB for the estimation problem proposed in Section 2. The CRB can serve as the benchmark for our proposed estimator in Section 4. The CRB gives a lower bound for the mean square error (MSE) matrix of an unbiased estimator of unknown parameter vector u o [61]. The square root of the trace of CRB matrix has the physical meaning that it is the minimum possible root-mean-square error (RMSE) for positioning the target unbiasedly. Evaluation of CRB is very useful because it not only gives us a performance benchmark, but also provides us insight into the system design.
Because the observation model (4) is a Gaussian model, it is easy to get the CRB for the estimation of u o as follows [61].

Design of the Estimator
Maximum likelihood estimation is overwhelmingly the most popular approach to obtaining practical estimators. Our design procedure of MLE based on semidefinite relaxation (SDR) technique is divided into three steps. We first establish a quadratic optimization problem with multiple quadratic equality/inequality constraints based on the maximum likelihood criterion. In the second step we rewrite the quadratic optimization problem using Frobenius inner product. Finally, we drop the nonconvex constraints in the problem.

Formulation of MLE as a Quadratic Optimization Problem
Under the Gaussian noise assumption in (4) and with the matrix reformulation of (1), the NLL function of the estimation problem is where In (7), u is the optimization variable and ξ 1 does not depend on u. Here we write The NLL function is nonconvex with respect to the optimization variable u.
In the MLE sense, by introducing auxiliary optimization variable x, the estimator can be considered as the optimal solution of the constrained optimization problem as below.
subject to The objective function (8a) can be expanded as Then the original problem (8) can be rewritten as follows.
subject to It is ready to recognize the problem (9) as a quadratic optimization problem with multiple quadratic equality constraints [38,62]. While the objective function (9a) is convex, the constraint functions (9b) are nonconvex. In short, the MLE involves a nonconvex optimization problem. It is well known that computing the MLE via a grid search results in computationally prohibitive solutions and Gauss-Newton implementation of the MLE requires carefully chosen initial guesses [61]. For such a computationally difficult problem, our strategy is to approximate the problem (9) to a convex one using SDR technique. The related SDP can be solved in a polynomial complexity using interior-point methods without causing the risk of local convergence [63,64].

Some Equivalent Transformation
In order to introduce semidefinite constraints, some necessary equality transformations are made in this subsection. Rewriting the objective function (9a) in matrix form and reformulating the constraints (9b), the problem (9) is recast as subject to where Note that the parameter matrices F and G i are positive semidefinite. Then we introduce two new matrix variables U and X, reformulating the problem (10) as where We emphasize that by now the problem (12) is still equivalent to the original MLE problem (8), but not a convex one.

Semidefinite Relaxation
In this subsection, we relax the problem (12) to an SDP problem. The semidefinite programs can be solved almost as easily as linear programs with interior-point methods [64] and several advanced SDP solvers are readily available [65]. By the Schur complement condition for positive semidefiniteness [66], In view of (14) and (15), if the constraints (12b) and (12c) in (12) are relaxed to drop the (nonconvex) rank-one constraints, then we obtain the following (convex) SDR problem.
where Z = blkdiag( X, U, diag(x)) is the primal variable. Now note that the problem (16) is an approximation of the original problem (8), but a more tractable one. We call the minimizer (if unique) of the SDP problem (16) as an SDR solution to the original elliptical location problem.

Semidefinite Programming Solver
In order to design an efficient numerical optimization algorithm, we first find the dual problem of the primal problem (16). By the conic duality theorem [41], the dual of problem (16) is as follows. maximize S,y y 1 + y 2 , subject to −c 2 diag(β) α/2 where y = [y 1 , y 2 , α T , β T ] T and S = blkdiag(S 1 , S 2 , diag(s)) are the dual variables. The lengths of vectors α and β are both (M + N).
From the perspective of optimization theory, the rest of the work is to reduce the duality gap of primal-dual problem pairs (16) and (17) to find the optimal solution. We will solve the primal-dual problem pairs with an interior-point primal-dual path-following method.
To develop efficient interior-point algorithms, the key is to find suitable barrier functions. With an appropriate choice of barrier functions, the interior-point algorithm admits a polynomial time implementation. We choose the canonical barrier functions for our problems (16) and (17), and then the parameter of barrier is n = 2(M + N + 1) + d [41]. The barrier functions are added to the objective functions, forming primal and dual barrier problems. Given a barrier weight parameter µ > 0, the central path associated with these barrier problems is defined as the set of points C SDP (µ) = {(Z, S, y) strictly feasible in (16) and (17) Now all we need for obtaining good primal and dual approximate solutions is to trace fast the central path (18). In practice, path following (i.e., computing primal-dual search directions) can be implemented in a moderate (few tens) number of iterations, each iteration reducing to assembling and solving a Newton system of linear equations resulted from linearization of (19) below.
the equality constraints in the primal problem (16) the equality constraints in the dual problem (17) In summary, our algorithm is shown in Algorithm 1, where (Z 0 , S 0 , y 0 ) are the initial values of primal and dual variables, and > 0 is a tolerance parameter for balancing the accuracy and speed. One advantage of using SDR technique is that the convergence of the algorithm is guaranteed, regardless of the initial values. Compute primal-dual search directions ∆y, ∆S and ∆Z using Newton's method on (19) 6: Z ← Z + ∆Z, S ← S + ∆S, y ← y + ∆y, µ ← Z, S /n 7: end while 8: end procedure Using (13), the minimizer of (16) (i.e., the SDR solution) contains an estimate of target location u o . Strictly speaking, if rank( X) = rank( U) = 1 holds, we do not need to further refine the SDR solution. Otherwise, some post-processing techniques are required to improve the accuracy of the estimate [45]. Fortunately, our experiments in Section 5 show that our algorithms require little post-processing.

Complexity Analysis
Now we give a complexity analysis to show that our algorithm complexity is controllable as the numbers of transmitters and receivers increase. Since our estimator is implemented using the interior-point method, its complexity analysis can be based on some established theoretical results [40,41]. The worst-case Newton complexity of an interior-point primal-dual path-following method is O(1)n 1 2 ln 1 . We should also know how heavy a Newton step is, i.e., what is its arithmetic cost of each iteration. We count the numbers of semidefinite matrix blocks and equality constraints in the SDP problem (16). In total, there is one semidefinite matrix block of size (M + N + 1), one semidefinite matrix block of size (d + 1), and (M + N) non-negative scalar variables. Furthermore, there are 2(M + N + 1) equality constraints on these matrix/scalar variables. When (M + N + 1) d, the flop count in each Newton step is approximately 6(M + N + 1) 4 .
In conclusion, the total complexity of our algorithm is on the order of (M + N + 1) 4.5 ln 1 .

Results and Discussion
We conducted a series of Monte Carlo simulations to evaluate the proposed SDR estimator comprehensively. Generally, when using convex relaxation techniques for parameter estimation, three aspects of the estimator are of concern: What is the statistical performance of the estimator, What is the degree of relaxation of the related optimization model, and whether the relaxed solution needs post-processing.

1.
To ascertain the statistical performance of the proposed SDR estimator, we compared it with several closed-form methods recently proposed, i.e., Yang's TS-WLS method [21] and Amiri's CWLS method [20]. For a given estimatorû of the unknown parameter u o , its performance is measured by the empirical bias and RMSE, which are defined as follows.
where L is the number of Monte Carlo simulations and u is the -th random realization ofû. In our following experiments, we set L = 10 4 and compared the biases and RMSEs of all these estimators at different noise levels.

2.
To measure the tightness of the relaxed optimization model in (16), we recorded the condition numbers of the optimal solutions of X and U in (17) at different noise levels. The motivation is that if the condition number of a matrix is greater than 10 4 , then it can be approximated as a rank-1 matrix [44].

3.
To verify whether the SDR solution needs post-processing, we checked the performance improvement gifted by post-processing using local search under different noise conditions. We used the SDR solution as the initial point for the commercial optimization solver MATLAB function lsqnonlin, and the algorithm for lsqnonlin was set as 'trust-region-reflective'. Our estimator with post-processing is simply referred to as SDR-LS.

Localization Scenarios
Our experimental design was as follows. We used M = 3 transmitters and N = 5 receivers to determine the position u o of a stationary object. Their positions are given in Table 1. The location geometry was the same as the one in [22] and depicted in Figure 2.
The signal propagation speed was c = 1500 m/s. The covariance matrix of the observation vector τ i related to the i-th transmitter was Q τ i = σ 2 τ R for i = 1, 2, . . . , M, where σ τ was a given positive constant, and R = 1 2 (1 N×N + I N ) [34].

Visualization of Cramér-Rao Bound and Selection of Testing Points
We visualized the CRB given by (5), which sets a lower bound for the MSE of our estimation problem proposed in Section 2. We selected two representative positions as test points according to the square root of the trace of CRB matrix. Based on these considerations, we plotted the contours of the square root of the trace of CRB matrix for σ τ = 0.2 as in Figure 3. The positioning difficulty of the target at different locations indicated by the contours in Figure 3 is consistent with our intuitive experience. According to the contours, we selected easier point [0, −1000] m (with the minimum square root of trace of CRB matrix) and harder point [−5000, 5000] m (with the maximum square root of trace of CRB matrix) as the two representative test points for further in-depth comparative study in next subsection.

Results on the Two Representative Test Points
In order to examine the performance of various algorithms at different noise levels, we increased measurement noise level σ τ from 0.01 s to 0.1 s with a step size of 0.01 s in the following simulation experiments.  It can be seen that when the noise level increases, there is a sharp increase in the bias of CWLS estimator, while the other three algorithms control the bias very well over a wide range of noise levels. Furthermore, a closer look at Figure 4 reveals that the bias of TS-WLS estimator is still moderately larger than that of SDR estimator and SDR-LS estimator. We believe that these findings are due to the following two facts. First, the square terms of observation errors are ignored in the models of CWLS estimator and TS-WLS estimator. In addition, CWLS estimator involves solving high-order polynomial equations with numerical instability.

Results on Easier Test Point
Due to their excellent biases, we limited our attention to SDR estimator and SDR-LS estimator. Figure 5 shows that the RMSEs of both SDR estimator and SDR-LS estimator agree well with CRBs at all noise levels. This implies that our SDR estimator has achieved the theoretically optimal RMSE performance at this test point and it does not need to use local search for post-processing.
To further confirm that the SDR model is tight enough, we showed the maxima and minima of the condition numbers of optimal X and U in the Monte Carlo experiments at different noise levels in Figure 6. As can be seen, the condition numbers are much larger than the critical value of 10 4 proposed by [44]. As a practical guideline for deploying our SDR estimator in actual systems, we drew the scatter plots of its Monte Carlo random realizations at even higher noise levels in Figure 7. It can be seen that, no matter how high the noise level is, the average value of a huge large number of estimation results is highly consistent with the actual coordinates of the target. As shown in Figure 8, the bias of SDR-LS estimator is the smallest, and that of CWLS estimator is the largest. What is worth noting in Figure 8 is that at all noise levels, SDR estimator's bias is slightly larger than that of SDR-LS estimator. Furthermore, a similar situation was observed when we examined the RMSEs in Figure 9. It must be pointed out that SDR-LS estimator's RMSE agrees well with CRB here. These results suggest that SDR estimator provides a reliable initial solution for local search, and local search based on this initial solution achieves theoretically optimal performance. Since SDR estimator is near optimal, we suggest that the practitioners decide whether post-processing is required based on time requirements and accuracy requirements in practice. As in Figure 6 for the easier test point, Figure 10 confirms from the perspective of condition numbers that the SDR model is also tight enough at this harder test point. As with the easier test point, to ensure the reliability of the actual deployment, we drew the scatter plots of SDR estimator's Monte Carlo random realizations at even higher noise levels in Figure 11. Not surprisingly, increasing the standard deviation of the measurement noise did not cause significant deviations in the position estimation results.

Conclusions
This work develops an SDP-based estimator (i.e., SDR estimator) for elliptical location problems with multiple transmitters and multiple receivers, implementing advanced optimization techniques in critical signal processing in wireless communication systems. The performance degradation of the SDR estimator is almost negligible compared with the statistical performance of the optimal estimator. Its performance is more pronounced at higher noise levels, ensuring its deployability. As with the other closed-form algorithms from the literature, the SDR solution has neither the difficulties in initialization, nor the problems of divergence and local convergence. In particular, the SDR solution is computationally efficient owing to the polynomial complexity of the interior-point method.
Based on the modelling flexibility of SDP, our future efforts are to extend the algorithm in this paper to scenarios (1) where the locations of transmitters or receivers are uncertain and (2) where signal propagation paths involve non-line-of-sight interference.
Author Contributions: Conceptualization, X.W. and L.Y.; methodology, X.W.; software, X.W. and Y.D.; validation, L.Y.; writing-original draft preparation, X.W. and Y.D.; writing-review and editing, X.W. and L.Y.; funding acquisition, X.W. and L.Y. All authors have read and agreed to the published version of the manuscript.

Acknowledgments:
The authors would like to thank the anonymous reviewers and the editor for their careful reviews and constructive suggestions to help us improve the quality of this paper.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: