A Block Coordinate Descent Method for Sensor Network Localization

The problem of sensor network localization (SNL) can be formulated as a semidefinite programming problem with a rank constraint. We propose a new method for solving such SNL problems. We factorize a semidefinite matrix with the rank constraint into a product of two matrices via the Burer--Monteiro factorization. Then, we add the difference of the two matrices, with a penalty parameter, to the objective function, thereby reformulating SNL as an unconstrained multiconvex optimization problem, to which we apply the block coordinate descent method. In this paper, we also provide theoretical analyses of the proposed method and show that each subproblem that is solved sequentially by the block coordinate descent method can also be solved analytically, with the sequence generated by our proposed algorithm converging to a stationary point of the objective function. We also give a range of the penalty parameter for which the two matrices used in the factorization agree at any accumulation point. Numerical experiments confirm that the proposed method does inherit the rank constraint and that it estimates sensor positions faster than other methods without sacrificing the estimation accuracy, especially when the measured distances contain errors.


Introduction
Sensor network localization (SNL) is the problem of estimating the unknown positions of m sensors from the known positions of n anchors and the measured distances (which may contain measurement errors) between sensor-sensor or sensor-anchor pairs. In this problem, let E ss and E sa be the sets of sensorsensor and sensor-anchor pairs, respectively, with known measured distances, and let d ij and d ik be the measured distances for each ij := {i, j} ∈ E ss and ik := {i, k} ∈ E sa , respectively. Let the anchor coordinates be a k = (a k1 , . . . , a kd ) ⊤ ∈ R d (k = m + 1, . . . , m + n). Then, the SNL problem is formulated as the following system of equations with variables x i ∈ R d (i = 1, . . . , m): When a matrix variable Z is introduced, finding the x 1 , . . . , x m satisfying system (1) is known to be equivalent to finding a solution of the following semidefinite programming (SDP) problem with a rank constraint [4,21]: Here, for each ij ∈ E ss and ik ∈ E sa , where e 1 , . . . , e m is the canonical basis of R m . See Subsection 1.2 for the definitions of the other symbols. When d ij and d ik contain measurement errors, problem (2) generally does not have a solution. Therefore, to account for the case in which d ij and d ik contain errors, researchers have often considered problem (3) defined below [4,8,9,10,16,19]. In this problem, the distance constraints are incorporated in the objective function in the form of quadratic errors, and a penalty is imposed for violating those constraints. In this paper, we also seek to estimate sensor positions in SNL by solving this problem: The SNL problem is generally known to be NP-hard [1], and the formulations of the SNL problem as optimization problems (2) and (3) are also nonconvex. This nonconvexity is due to the rank constraint that appears in problems (2) and (3). Therefore, many previous SNL studies removed the rank constraint and relaxed the problem into an SDP problem to estimate approximate sensor positions [2,3,4,5,14,22]. Among these methods, sparse full SDP (SFSDP) as proposed by Kim et al. [13] is an especially representative one. However, a solution of the SDP relaxation problem is not always a solution of the original problem. So and Ye [18] referred to problem (1), which has a unique solution and does not have sensor positions satisfying all of the given distances in a higher-dimensional space, as "uniquely localizable." They proved that problem (1) is uniquely localizable if and only if the maximum rank of the solutions of the SDP relaxation problem (2) is d. Therefore, when the SDP problem is solved by the interior-point method, which is a representative method for solving general SDP problems, if the problem is uniquely localizable, then the exact sensor positions can be determined by solving the SDP relaxation problem. Otherwise, an optimal solution of the SDP relaxation problem corresponds to a configuration of the sensors in a higher-dimensional space because of the max-rank property of the interior-point method [11], and it thus might give poorly estimated sensor positions in d-dimensional space.
On the other hand, the exact sensor positions in d-dimensional space can be estimated only if problem (1) has a unique solution (Wan et al. [21] referred to a problem satisfying this condition as "locatable"). Therefore, methods have recently emerged for estimating the sensor positions by solving problem (2) or (3) itself. Wan et al. [21] proposed a method that obtains a solution of problem (2) by solving SDP problems multiple times. This approach is based on the fact that the rank of a matrix being less than or equal to d is equivalent to its (d+1)th and subsequent eigenvalues all being zero. Numerical experiments showed that their method's estimation accuracy was better than that of an SDP relaxation-based method. However, their method took more time than the latter method, because it requires solving the SDP problem via an SDP solver at each iteration. Wan et al. [20] also proposed a method that transforms the SDP problem with the rank constraint into an SDP problem with a complementarity constraint and alternately performs minimization with regard to the two semidefinite matrices that appear in the problem. Numerical experiments also confirmed that this method was more accurate than SDP relaxation methods such as SFSDP. However, as in [21], it took too long to estimate the sensor positions.
Another method for solving general SDP problems is the Burer-Monteiro factorization, in which a semidefinite matrix Z is factorized into the form V V ⊤ and a nonconvex optimization problem is solved after the factorization [7]. If the number of columns of V is chosen as r in this factorization, then it introduces the constraint that the rank must be less than or equal to r. Therefore, this method is suitable for obtaining low-rank solutions of SDP problems. In a series of studies, Chang and colleagues [8,9,10] attempted to estimate sensor positions by using the Burer-Monteiro factorization. First, Chang and Xue [9] proposed a method that applies the limited-memory Broyden-Fletcher-Goldfarb-Shanno method to the problem after the Burer-Monteiro factorization. However, they set the number of columns of V used in the factorization to that of the semidefinite matrix before the factorization, so their method does not consider the rank constraint. Second, Chang et al. [10] used the same method as in [9] to estimate sensor positions in three-dimensional space, but unlike in [9], they set the number of columns of V to three, so we can say that this method takes the rank constraint into account. They compared it with SFSDP through numerical experiments and reported that the sensor positions could be estimated more quickly and with the same level of accuracy as SFSDP. However, those experiments involved only small problems with up to 200 sensors. Finally, Chang and Liu [8] proposed a method that they called NLP-FD, which solves the optimization problem obtained from the Burer-Monteiro factorization by the curvilinear search algorithm [23]. Their numerical experiments showed the superiority of NLP-FD over SFSDP when a problem is large in scale and the measured distances include errors.
In this paper, we propose a new method for SNL that accounts for the rank constraint. First, we factorize Z in problem (3) into a product of two matrices through the Burer-Monteiro factorization: This factorization is equivalent under the constraint U − V = O. Therefore, problem (3) can be transformed into an unconstrained multiconvex optimization problem by adding the difference between the two matrices, with a penalty parameter γ, to the objective function. Then, the block coordinate descent method can be applied to the new objective function, and optimization can be performed sequentially for each column of U and V . We formalize this procedure as Algorithm 1.
We also analyze the proposed method theoretically. First, we show that each subproblem in Algorithm 1 is an unconstrained convex quadratic optimization problem and can be solved analytically (Theorem 2). Second, we show that any accumulation point of the sequence generated by Algorithm 1 is a stationary point of the objective function (Theorem 3). Third, we give a range of γ for which the two matrices U and V used in the factorization coincide at any accumulation point (Theorem 4). Finally, we explain the relationship between the objective function in the reformulated problem and the augmented Lagrangian. Numerical experiments confirm that the proposed method does inherit the rank constraint; furthermore, the results demonstrate not only that our method estimates sensor positions faster than SFSDP and NLP-FD without sacrificing estimation accuracy, especially when the measured distances include errors, but also that our method does not run out of memory even for large-scale SNL problems.
The rest of this paper is organized as follows. In Section 2, we present the proposed method and analyze it theoretically. In Section 3, we compare it with the other methods to confirm its effectiveness. Finally, in Section 4, we present our conclusions and suggest possible future work.

Notation
-N denotes the set of natural numbers without zero. R p denotes the set of p-dimensional real vectors, and 0 p denotes the zero vector of R p . When the size is clear from the context, we omit the size subscript at the lower right. R p×q denotes the set of p × q real matrices. Let I p be the identity matrix of R p×p and O be the zero matrix of an appropriate size. S p + denotes the set of p × p symmetric positive semidefinite matrices.
-For x ∈ R p , x 2 denotes the 2-norm of x. respectively, that are connected directly to sensor i.

Proposed method
Problem (3) is an SDP problem with a rank constraint and is difficult to solve directly. In this subsection, we propose a new method that transforms problem (3) into an unconstrained multiconvex optimization problem and solves the latter problem sequentially to estimate sensor positions. First, a matrix Z satisfies the three constraints in problem (3) if and only if it can be factorized into the product of two matrices as follows: Thus, problem (3) is equivalent to To make problem (4) easier to solve, we remove the constraint U − V = O and add a quadratic penalty term γ/2 U − V F with a penalty parameter γ (> 0) to the objective function; as a result, the objective function takes larger values as the constraint U − V = O is more strongly violated. In other words, we adopt the following unconstrained optimization problem: In the proposed algorithm, we let U = (u 1 , . . . , u m ) and V = (v 1 , . . . , v m ) and then perform minimization with regard to u 1 , . . . , u m , v 1 , . . . , v m , i.e., each column of U and V sequentially. Specifically, the procedure is as listed in Algorithm 1.

Algorithm 1 Proposed algorithm for problem (5)
Require: an initial point The proposed method has the following advantages over other methods: (i) The SDP problem with the rank constrained (3) is equivalent to problem (4), from which we obtained problem (5) by incorporating U −V = O in the objective function as a quadratic penalty term with a penalty parameter γ. As we will see from Theorem 4, if γ is larger than a real-valued threshold, then the U (p) and V (p) generated by Algorithm 1 coincide with each other at any accumulation point, thereby satisfying the constraint of problem (4). Therefore, the proposed method inherits the rank constraint in problem (3) and retains the potential capability to estimate sensor positions accurately for problems that are not uniquely localizable. This advantage will be verified in Subsection 3.1. (ii) As we will see from Theorem 2, each subproblem appearing inside a for statement in Algorithm 1 is an unconstrained convex quadratic optimization problem. The solution of each subproblem can be obtained analytically, because the solution process can be reduced to solving a system of linear equations with an invertible coefficient matrix of size d.
Because d is at most three in real situations, the system can be solved rapidly and without running out of memory, regardless of the number of sensors m. Moreover, the subproblems only need to be solved 2m times (i.e., a number proportional to m) for each outer loop. Therefore, especially in the case of large-scale SNL problems, we expect faster estimates of the sensor positions as compared with other methods. This advantage will be verified in Subsection 3.2.

Analyses of the proposed method
In this subsection, we present theoretical analyses of problem (5) and Algorithm 1. First, we impose an assumption about the problem that we are examining.
Assumption 1 All sensors are connected to an anchor either directly or indirectly.
The same assumption was also made in [8,18,22] and is very natural when estimating sensor positions: if a sensor is not connected to any anchors, then its absolute position cannot be determined uniquely. First, we prove that the optimal solution of each subproblem in Algorithm 1 can be obtained uniquely as an analytical solution.
Proof If we focus on only u i in F (U, V ; γ) in particular, then we can represent From equation (6), we can see that and A ui is positive definite and invertible in particular, we can obtain u * i . We can also obtain v * i from the same calculation. ⊓ ⊔ Next, we show that the sequence generated by Algorithm 1 converges to a stationary point of the objective function F . Theorem 3 Fix the penalty parameter γ in problem (5) arbitrarily. Let N be the set of stationary points of F . Then, the sequence {(U (p) , V (p) )} ∞ p=1 generated by Algorithm 1 satisfies In particular, any accumulation point (U * , V * ) of the generated sequence is a stationary point of F .
The consequence of Theorem 3 is based on the result of [25]. By Corollary 2.4 in [25], if all three of the following conditions are satisfied, then equation (7) holds when the stationary point of F in the definition of N is replaced by the Nash equilibrium of F (see Definition 2 below).
Condition (a): F is continuous, bounded below, and has a Nash equilibrium. Condition (b): The objective function of each subproblem is strongly convex. 1 Condition (c): The sequence generated by Algorithm 1 is bounded.
In the following, we prove Theorem 3 by showing the equivalence between the Nash equilibrium and the stationary point of F and then verifying that all three conditions are satisfied. We can see from equation (6) that the function F (U, V ; γ) is convex on R d if we focus only on each column of U and V . Such a function F with this property is called multiconvex [25]. For simplicity, let X := R n1 × · · · × R ns (where n 1 , . . . , n s ∈ N), 2 and when x ∈ X is represented as Definition 1 A function g : X → R is called multiconvex on X (with respect to the block division x = (x 1 , . . . , x s ) ∈ X ) if for all i = 1, . . . , s and all x j ∈ R nj (j = 1, . . . , i − 1, i + 1, . . . , m), the function One of the concepts of minimality for a multiconvex function is the Nash equilibrium [25], which appears in Condition (a).
Definition 2 For a function g : X → R, (x * 1 , . . . , x * s ) ∈ X is called a Nash equilibrium of g (with respect to the block division as in Definition 1) if Gorski et al. [12] proved the equivalence between the stationary point and the Nash equilibrium in the case of s = 2. 3 Herein, we extend the equivalence to the case of arbitrary s. Lemma 1 Let g : X → R be once differentiable and multiconvex. Then, Proof We begin by proving the "if" part. If we assume that x * is a stationary point of g, then because is convex on R ni for all i = 1, . . . , m, holds for all x i ∈ R ni . Because ∇ xi g i (x * i ) = 0 follows from the assumption of x * being a stationary point of g, we can say from inequality (8) Because i is arbitrary, we can conclude that x * is a Nash equilibrium of g.
Next, we prove the "only if" part. If we assume that x * is a Nash equilibrium of g, then for each i = 1, . . . , s, g i (x i ) attains its minimum value at Thus, x * is a stationary point of g.
⊓ ⊔ Lemma 2 Suppose that Assumption 1 holds. Then, for all α, the level set Because similar (but not the same) results were already pointed out in [10,18], we omit the proof of Lemma 2 because of the page limit. Note that if Assumption 1 does not hold, then S F (α) is always not bounded.
Proof The continuity and below-boundedness of F are evident from its definition. When combined with Lemma 2, the global optimal solution of problem (5) is guaranteed to exist, from which the existence of a Nash equilibrium can be proved. Therefore, we can say that Condition (a) holds. In addition, from equation (6), the objective function of each subproblem is γ-strongly convex, and thus Condition (b) holds. Finally, Condition (c) is also satisfied from Corollary 1. Therefore, the conditions of Corollary 2.4 in [25] are all satisfied, from which we can show equation (7) by using the equivalence between the stationary point and the Nash equilibrium that was shown in Lemma 1.
Because F is of class C 1 , N is closed, from which we can easily show the last part of Theorem 3.

⊓ ⊔
Finally, we show that the U (p) and V (p) generated by Algorithm 1 coincide with each other at any accumulation point if γ is larger than a real-valued threshold, which does not generally hold in the quadratic penalty method.
Proof For each p ∈ N, ij ∈ E ss , and ik ∈ E sa , let Then, because the initial point satisfies U (0) = V (0) and the value of the objective function F decreases monotonically by Algorithm 1, we can conclude that for all p ∈ N, By taking a subsequence, without loss of generality, we can assume that {(U (p) , V (p) )} ∞ p=1 itself converges to (U * , V * ). Then, it follows from inequality (10) that By Theorem 3, (U * , V * ) is a stationary point of F , from which we have Thus, let U l , V l ∈ R m be the respective lth-column vectors of U ⊤ , V ⊤ for each l = 1, . . . , d; then, ∇ U l F (U * , V * ; γ) = ∇ V l F (U * , V * ; γ) = 0 holds. For each ik ∈ E sa and l = 1, . . . , d, let b l ik be an m-dimensional vector such that its ith component is −a kl and all other components are zeros. Furthermore, for each ij ∈ E ss and ik ∈ E sa , letĀ ij andĀ ik respectively bē A ij := (A ij ) (d+1:d+m,d+1:d+m) ,Ā ik := (A ik ) (d+1:d+m,d+1:d+m) .
Using these symbols, we can represent F (U, V ; γ) as = 0 for all l = 1, . . . , d, we obtain For convenience, letĀ Next, we seek to prove the following inequality: In fact, if inequality (13) can be shown, then because γ satisfies inequality (9), 2γI m −Ā is a positive definite matrix, and thus, U * l = V * l from equality (12). Because of the arbitrariness of l, we can eventually conclude that U * = V * . Therefore, we need only prove inequality (13).
It follows from the Gershgorin circle theorem that For each i = 1, . . . , m, let v(i) be the optimal value of the following optimization problem: 4 where the right side of inequality (14) does not exceed max 1≤i≤m v(i) because of inequality (11). We can easily check that v(i) is equal to the optimal value of the following optimization problem for each i = 1, . . . , m: Using the method of Lagrange multipliers, we can see that the optimal value of problem (15) which implies inequality (13). ⊓ ⊔

Relationship to the augmented Lagrangian
In this paper, we adopt the quadratic-penalty-based method, in which the equality constraint U − V = O is incorporated in the objective function as a quadratic penalty term and the resulting new objective function is minimized.
On the other hand, there are also methods such as the augmented Lagrangian method [17] and the alternating direction method of multipliers [6] that minimize the augmented Lagrangian, which contains not only the quadratic penalty term but also the Lagrange multiplier term. It is known that the augmented Lagrangian method is more efficient than the quadratic penalty method. For example, while the quadratic penalty method requires the penalty parameter to diverge to positive infinity, the augmented Lagrangian method does not require it to diverge, and the sequence obtained by the augmented Lagrangian method converges faster than that obtained by the quadratic penalty method [17,Example 17.4]. Hence, we explain that problem (5) can be regarded as a minimization problem of the augmented Lagrangian with an exact Lagrangian multiplier. Λ = O is the exact Lagrange multiplier of problem (4). In fact, for all local optimum solutions (U * , V * ) of problem (4), because problem (4) satisfies the linear independence constraint qualification, there exists a Lagrange multiplier Λ * ∈ R d×m satisfying the Karush-Kuhn-Tucker condition. In other words, if we let the Lagrangian for problem (4) be hold. We get U * = V * from equation (16c) and denote both of them as W * . Be- Using this equation and equations (16a) and (16b), we obtain Λ * = O. Therefore, the augmented Lagrangian with the Lagrange multiplier Λ = O and penalty parameter γ is which is the definition of F (U, V ; γ) itself. Hence, problem (5) can be regarded as a minimization problem of the augmented Lagrangian with the exact Lagrange multiplier Λ = O for problem (4).

Numerical experiments
In this section, we use numerical simulation to verify the advantages (i) and (ii) described in Subsection 2.1 for the proposed method. We begin by confirming that our method does inherit the rank constraint; to confirm this, we compare it with an SDP relaxation-based method for a problem that is locatable but not uniquely localizable. Next, to confirm the effectiveness of the proposed method, we compare its estimation time and estimation accuracy with those of other methods by using artificial data under various conditions. All experiments were conducted on a computer with the macOS Catalina operating system, an Intel Core i5-8279U 2.40 GHz CPU, and 16 GB of memory. All the algorithms were implemented using MATLAB (R2020a). The parameter ǫ in Algorithm 1 was set to 10 −5 throughout the experiments.

Comparison with SFSDP for a problem that is not uniquely localizable
In this subsection, we demonstrate that the proposed method has the capability to estimate sensor positions accurately for a problem that is locatable but not uniquely localizable. Specifically, we examine a problem from [18]: = (0.6, 0.7) ⊤ . For this problem, Algorithm 1 was executed after fixing the penalty parameter γ as 2f (U (0) , V (0) ) max 1≤i≤m 4|E ss | + |E sa |/2 according to Theorem 4. We examined the two cases of whether the initial points u in Algorithm 1 are in the interior or the exterior of the convex hull of the three anchors. Figure 1 shows the sensor positions estimated by SFSDP, the SDP relaxation-based method described in Section 1, and the proposed method. Note that when we estimated the sensor positions with the proposed method, the randomness of the initial points was varied 10 times. The results were similar to those of Figure 1b in all cases in which the initial points were in the interior of the convex hull of the anchors. On the other hand, the results were similar to those of either Figure 1c or Figure 1d in all cases in which the initial points were in the exterior of the convex hull. Accordingly, only these three cases are included in Figure 1.
When we used SFSDP, the sensor positions were not estimated correctly (Figure 1a). For the proposed method, the estimation accuracy depended on the initial points. When the initial points were in the interior of the convex hull of the anchors, the sensor positions were estimated accurately (Figure 1b). On the other hand, when the initial points were in the exterior of the convex hull, the sensor positions were estimated accurately in some cases (Figure 1c) but not in others (Figure 1d).
Of course, because problem (5) examined in this paper is a nonconvex optimization problem, whether the sensor positions can be estimated accurately depends on the initial points. However, as shown in Figure 1b and Figure 1c, the proposed method still has the capability to estimate sensor positions accurately even for problems that are not uniquely localizable, although the example here is quite simple. On the other hand, when we use SDP relaxation-based methods, if a given problem is not uniquely localizable, there is no capability for accurate estimation because of the max-rank property of the interior-point method, as described in Section 1. Therefore, we can say that the proposed method does inherit the rank constraint.

Comparison of estimation time and accuracy
In this subsection, we quantitatively compare the estimation time and the estimation accuracy of the proposed method with those of existing methods for sensors located in two-or three-dimensional space. The compared methods are SFSDP, which was also used in Subsection 3.1, and NLP-FD, which takes the rank constraint into account as our proposed method does. Although we introduced the methods proposed by Wan et al. [20,21] in Section 1, which also account for the rank constraint, we do not compare them here because of their extremely low scalability. In these experiments, m = 1000, 3000, 5000, and 20000 sensors and n = 0.1m anchors were placed randomly in [0, 1] d . E ss and E sa were defined as where the x true i (i = 1, ..., m) are the sensors' true positions. In other words, we considered a model in which the distance between two sensors or between a sensor and an anchor is observed if and only if it is less than a radio range ρ (> 0). We set ρ to 0.1 and 10/m in the case of d = 2 and to 0.25 and 3 15/m in the case of d = 3. The measured distances d ij (ij ∈ E ss ) and d ik (ik ∈ E sa ) were given by where ǫ ij , ǫ ik were selected independently from the standard normal distribution, and σ is a noise factor determining the influence of the error. σ was set to 0, 0.1, and 0.2. As an indicator to measure the estimation accuracy, we used the root-mean-square distance (RMSD), which has been used in many other papers on SNL [8,13,20,22] and is defined as wherex i is the estimated position of sensor i. For each set of (m, n, ρ, σ), five different problems of varying randomness were created, and the final results were the averages of five measurements of the estimation time (CPU time) and the estimation accuracy (RMSD). The initial point (U (0) , V (0) ) ∈ R d×m × R d×m in Algorithm 1 was decided similarly to the method in [8]. That is, for each sensor i (i = 1, . . . , m), if it was connected directly to an anchor, then u The penalty parameter γ was updated dynamically according to Theorem 4 by the following procedure.
Step 4. Let W := (U (p) + V (p) )/2 and Restart Algorithm 1 with W as the initial point and γ as the penalty parameter.
The method of updating γ described above consists of two components. First, in Steps 1-3, γ is updated so that the value of f , which represents the squared error of the squared distances, decreases rather than the penalty term U − V F . However, if we keep reducing the value of f rather than the penalty term, then the penalty term does not decrease much, and U and V may end up taking very different values from each other. Therefore, in Step 4, we try to reduce the difference between U and V by fixing γ according to Theorem 4. First, the results for d = 2 are given in Table 1, 6 wherein the proposed method is referred to as "BCD." We can see that when the measured distances included no errors (σ = 0), the estimation time of the proposed method was the lowest in most cases; furthermore, even when the proposed method was not the fastest, its estimation time was almost the same as that of the fastest method. In terms of the estimation accuracy, SFSDP estimated the sensor positions with the best accuracy of all the methods, by an order of magnitude. However, comparing the proposed method and NLP-FD shows that there was no appreciable difference between their estimation accuracies. When the measured distances included errors (σ = 0.1 and 0.2), the proposed method estimated the sensor positions the most rapidly of all the methods, by an order of magnitude in all cases, and the estimation accuracy was about the same as those of the other two methods. Next, the results for d = 3 are given in Table 2. 7 The results for the threedimensional scenario were similar to those for the two-dimensional scenario; that is, when the measured distances did not include errors, the estimation time of the proposed method was the lowest in each case. In terms of the estimation accuracy, SFSDP estimated the sensor positions with the highest accuracy of all the methods, by an order of magnitude; however, comparing the proposed method and NLP-FD again shows that there was no appreciable difference between their estimation accuracies. When the measured distances included errors, the estimation time of the proposed method was the lowest in all cases except (m, n, ρ, σ) = (1000, 100, 0.1, 0.2), and even in that case, its estimation time was also almost the same as that of the fastest method (SFSDP). The estimation accuracy of the proposed method was also comparable to those of the other two methods. In addition, for m = 20000, SFSDP could not estimate the sensor positions because of insufficient memory, while the proposed method could estimate the positions without running out of memory.
Overall, these results for the two-and three-dimensional cases show that the proposed method has practical advantages over the other methods: it can estimate sensor positions faster than those methods can without sacrificing the estimation accuracy, especially when measurement errors are included, and it does not run out of memory even for large-scale SNL problems.
It is interesting to consider why the proposed method and NLP-FD, both of which account for the rank constraint, could not estimate the sensor positions with as much accuracy as SFSDP when there were no measurement errors. The reason was probably because a formulation that accounts for the rank constraint is a nonconvex optimization problem and thus might converge to a stationary point that is not a global optimal solution. In contrast, if a problem is uniquely localizable, then SFSDP can estimate accurate sensor positions because the convergence to the global optimum of a relaxation problem is guaranteed. On the other hand, when measurement errors are included, even if a global minimum solution of the objective function f is obtained and the optimal value is zero, it does not mean that the true sensor positions are estimated, but rather that the positions are estimated incorrectly. In other words, even if the objective function is strictly minimized, it does not necessarily mean that a good estimate of the sensor positions is obtained; thus, when measurement errors are included, estimation accuracy comparable to that of SFSDP can be obtained even with methods that account for the rank constraint and may cause the generated sequence to fall into a stationary point that is not a global optimal solution.

Conclusion
In this paper, we proposed a new method that transforms the formulation of problem (3), which appears in SNL, into an unconstrained multiconvex optimization problem (5), to which the block coordinate descent method is applied. We also presented theoretical analyses of the proposed method. First, we showed that each subproblem that appears in Algorithm 1 can be solved analytically. In addition, we showed that any accumulation point (U * , V * ) of the sequence {(U (p) , V (p) )} ∞ p=1 generated by the proposed algorithm is a stationary point of the objective function of problem (5), and we gave a range of γ such that (U * , V * ) satisfies U * = V * . We also pointed out the relationship between the objective function of problem (5) and the augmented Lagrangian. Numerical experiments showed that our method does inherit the rank constraint and that it can estimate sensor positions faster than other methods without sacrificing the estimation accuracy, especially when the measured distances contain errors, and without running out of memory.
The present study suggests three directions for future work. First, Algorithm 1 uses a cycle rule in which the 2m subproblems are solved in the order of u 1 , . . . , u m , v 1 , . . . , v m . However, in the general coordinate descent method, there are other update rules such as a random rule and a greedy rule [15,24]. In SNL, the strategy of updating from variables corresponding to sensors that are connected directly to anchors is also expected to improve the estimation accuracy and time. Therefore, there is still room to consider how the order of solving the 2m subproblems affects the estimation time and accuracy. Second, we performed the minimization sequentially with respect to each column of U and V for computational efficiency, but updating some columns of U and V together is also possible, and the manner of block division in applying the block coordinate descent method should be examined further. Finally, the proposed method could be extended to general quadratic SDP problems with a rank constraint.

Acknowledgment
We thank Dr. Xiaokai Chang for providing the NLP-FD MATLAB code. This is a preprint of an article published in Optimization Letters. The final authenticated version is available online at: https://doi.org/10.1007/s11590-021-01762-9.

Funding
This work was supported by JSPS KAKENHI Grant Number JP20H02385.