Robust Euclidean embedding via EDM optimization

This paper aims to propose an efficient numerical method for the most challenging problem known as the robust Euclidean embedding (REE) in the family of multi-dimensional scaling (MDS). The problem is notoriously known to be nonsmooth, nonconvex and its objective is non-Lipschitzian. We first explain that the semidefinite programming (SDP) relaxations and Euclidean distance matrix (EDM) approach, popular for other types of problems in the MDS family, failed to provide a viable method for this problem. We then propose a penalized REE (PREE), which can be economically majorized. We show that the majorized problem is convex provided that the penalty parameter is above certain threshold. Moreover, it has a closed-form solution, resulting in an efficient algorithm dubbed as PREEEDM (for Penalized REE via EDM optimization). We prove among others that PREEEDM converges to a stationary point of PREE, which is also an approximate critical point of REE. Finally, the efficiency of PREEEDM is compared with several state-of-the-art methods including SDP and EDM solvers on a large number of test problems from sensor network localization and molecular conformation.


Introduction
This paper aims to propose an efficient numerical method for the most challenging problem in the Multi-Dimensional Scaling (MDS) family, which has found many applications in social and engineering sciences [6,10]. The problem is known as the Robust Euclidean Embedding, a term borrowed from [8]. In the following, we first describe the problem and its three variants. We then explain our approach and main contribution. We will postpone the relevant literature review to the next section in order to shorten the introduction.

Problem description
The problem can be described as follows. Suppose we are given some dissimilarity measurements (e.g., noisy distances), collectively denoted as δ i j , for some pairs (i, j) among n items. The problem is to find a set of n points x i ∈ r , i = 1, . . . , n such that where x is the Euclidean norm (i.e., 2 norm) in r and E is the set of the pairs (i, j), whose dissimilarities δ i j > 0 are known (E can be thought of the edge set if we treat δ i j as a weighted edge distance between vertex i and vertex j, resulting in a weighted graph.) Throughout, we use ":=" or "=:" to mean "define". The space r is called an embedding space and it is most interesting when r is small (e.g., r = 2, 3 for data visualization). One may also try to find a set of embedding points such that: A great deal of effort has been made to seek the best approximation from (1) or (2). The most robust criterion to quantify the best approximation is the Robust Euclidean Embedding (REE) defined by min X f (d,1) where W i j > 0 if δ i j > 0 and W i j ≥ 0 otherwise (W i j can be treated as a weight for the importance of δ i j ), and X := [x 1 , . . . , x n ] with each x i being a column vector. In [1,8], Problem (3) was referred to as a robust variant of MDS and is denoted as rMDS.
We will reserve rMDS for the Robust MDS problem: The reference rMDS for the problem (4) is more appropriate because it involved the squared distances D i j , which were used by the classical MDS [22,29,43,49,53]. The preceding two problems are robust because of the robustness of the 1 norm used to quantify the errors [31,Sect. IV]. When the least squares criterion is used to (1), we have the popular model known as the Kruskal's stress [30] minimization: Similarly, when the least-squares criterion was applied to (2), we get the so-called squared stress [6]: In many applications such as molecular conformation [21], lower and upper bounds data on the distances can also be collected: where 0 ≤ L i j ≤ U i j . In applications such as nonlinear dimensionality reduction [47] and sensor network localization [44,54], upper bounds U i j can be computed by the shortest path distances and L i j are simply set to be zero. According to [8,Sect. 5.1], all of those problems are NP-hard. However, some problems are computationally more "difficult" to solve than the others. The most challenging one, which is also the main focus of this paper, is the problem (3) with/without the constraint (7). The difficulty comes from the nonsmooth term of 1 norm and the distance terms d i j used. All other problems either involve the squared distances D i j or the squared 2 norm, which make them "easier" to approximate. We will explain the reasons in the literature review part.
In contrast to all other three problems, there lacks efficient methods for the REE problem (3). One of the earliest computational papers that discuss this problem is Heiser [23], which is followed up by [28], where the Huber smoothing function was used to approximate the 1 norm near zero with a majorization technique. It was emphasized in [28] that "the function is not differentiable at its minimum and is hard to majorize, leading to a degeneracy that makes the problem numerically unstable". Another important method is the PlaceCenter (PC for short) algorithm studied in [1]. We will compare with it in the numerical part. The difficulty in solving (3) is also well illustrated by a sophisticated Semi-definite Programming (SDP) approach in [35,Sect. IV] (see the literature review part). We now describe our approach proposed in this paper.

Our approach and main contributions
Our approach heavily makes use of the concept of Euclidean Distance Matrix (EDM). We need some notation. Let S n denote the space of all n × n symmetric matrices, endowed with the standard inner product. The induced norm is the Frobenius norm, denoted by A for A ∈ S n . The (i, j)th element of A ∈ S n is often written as A i j . Let S n + be the cone of positive semidefinite matrices in S n and we write A 0 for A ∈ S n + . A matrix D ∈ S n is called an EDM if there exists a set of points x i ∈ r , i = 1, 2 . . . , n such that the (i, j)th element of D is given by D i j := x i − x j 2 , i, j = 1, . . . , n. The smallest dimension r is called the embedding dimension of D and r = rank(J D J ), where J := I − 1 n 11 T is known as the centring matrix with I being the identity matrix in S n and 1 being the vector of all ones in n . We use D n to denote the set of all Euclidean distance matrices of size n × n.
The reformulation well captures four difficulties in solving the REE problem (3).
(i) The objective function f (D) is not convex. The term | D i j − δ i j | is convex when δ 2 i j > D i j and concave otherwise. (ii) The objective function is nonsmooth. It is not differentiable at certain points due to the 1 norm and the square root operation involved. (iii) The objective function is not Lipschizian. The Lipschitz constant goes to infinity as D i j goes to zero. The implication is that the subdifferential of the objective function [42,Def. 8.3] may be unbounded. This would create a huge obstacle in establishing any convergence results of iterative algorithms for (11).
(iv) The rank constraint is not convex and is hard to approximate. This is a common issue for any optimization problem with a rank constraint.
We note that no matter what reformulations one may use for (3), those four difficulties would appear in different forms and won't go away. We also note that all other three problems, when reformulated in terms of EDM, have a convex objective function. This distinctive feature alone makes the problem (11) the most challenging one to solve.
Existing numerical experiments have shown that the MDS embedding (9) and (10) works well as long as D is close to a true EDM. A typical example is when the data sits on a lower-dimensional manifold [47]. Motivated by this, we are going to generate an approximate EDM instead of a true EDM in our algorithm. It follows from (8) that (also see [32, Theorem A]): where K n + is known to be the conditionally positive semidefinite cone: and 1 ⊥ is the subspace in n orthogonal to 1. The diagonal constraint in (12) can be integrated to the set B with the choice L ii = U ii = 0 for i = 1, . . . , n. We combine K n + with the rank constraint into the set K n + (r ): We call it the conditionally positive semidefinite cone with the rank-r cut. Consequently, the constraints in (11) become −D ∈ K n + (r ) and D ∈ B. Next, we quantify the feasibility of −D belonging to K n + (r ) as follows. Let Π B K n + (r ) (A) be the set of all nearest points in K n + (r ) from a give matrix A ∈ S n . That is Since g(A) is just the half of the squared distance from (−A) to K n + (r ), it does not depend on which element Π K n + (r ) (A) is being used. It is easy to see that Hence, the problem (11) is equivalent to This is a classical constrained optimization problem with an equality constraint and a simple box constraint. Therefore, the quadratic penalty method [34,Chp. 17] can be applied to get the following problem: where ρ > 0 is the penalty parameter. We refer to this problem as the penalized REE problem (PREE). The quadratic penalty method is used often in practice [34, P. 497]. In fact, it is particularly suitable to (11) because it overcomes all four difficulties discussed above. We will need two more important tools to help us efficiently solve the penalty problem (16). One is the majorization technique that has recently become very popular in engineering sciences [46] (also see [6,Chp. 8] for its extensive use in MDS). Suppose we have the current iterate D k . We construct a majorization function g m (D, and The majorization is constructed in such a way that it is easier to solve the majorized problem: It can be seen that Hence, the algorithm generates a sequence {D k } that is nonincreasing in f ρ (D). Since f ρ (D) is bounded below by 0, the functional sequence { f ρ (D k )} converges. However, we are more concerned where the iterate sequence {D k } converges. The second concern is how the subproblem (18) has to be solved. This brings out the second technique, which is to solve the following one-dimensional problem: for given δ > 0 and 0 ≤ a ≤ b. We will show that the solution of this problem will lead to a close-form solution of (18). Since our method is for the Penalized REE by EDM optimization, we call it PREEEDM. The major contribution of this paper is to make the outlined solution procedure water-tight. In particular, we will investigate the relationship between the PREE problem (16) and the original problem (11) in terms of the -optimality (Proposition 1). We will also show that the majorization function g m (·, ·) can be economically constructed (Sect. 3.2). Moreover, the majorized function f k ρ (D) is guaranteed to be convex provided that the penalty parameter is above certain threshold and the subdifferentials at the generated sequences are bounded (Proposition 4). Furthermore, each majorization subproblem has a closed form solution (Theorem 1). We are also able to prove that any accumulation of the generated sequence by PREEEDM is a stationary point of (16), which is also an approximate stationary point of (11) (Theorem 2). Built upon its solid convergence results and simple implementation, PREEEDM is demonstrated to be comparable to six state-of-the-art software packages in terms of solution quality and outperform them in terms of the computational time for a large number of tested problems from sensor network localizations and molecular conformations.

Organization of the paper
In the next section, we give a selective literature review mainly on the Semi-Definite Programming (SDP) and EDM approaches. In Sect. 3. we introduce some necessary background and prove a key technical result (Lemma 1) that is crucial to the convexity of the majorization subproblem. We study the relationship between the penalized REE (16) and the original REE in Sect. 4, where the majorized subproblem is shown to have a closed-form solution. In Sect. 5, we provide a complete set of convergence results for the proposed PREEEDM algorithm. Numerical experiments are included in Sect. 6. The paper concludes in Sect. 7. All proofs except that of Theorem 2 can be found in "Appendix".

Literature review
One can find a thorough review on all of the four problems in [17] by France and Carroll, mainly from the perspective of applications. One can also find valuable discussion on some of those problems in [2]. So the starting point of our review is that those problems have their own reasons to be studied and we are more concerned how they can be efficiently solved.
Most of existing algorithms can be put in three groups. The first group consists of alternating coordinates descent methods, whose main variables are x i , i = 1, . . . , n. A famous representative in this group is the method of SMACOF for the stress minimization (5) [13,14]. The key idea is to alternatively minimize the function f (d,2) with respect to each x i , while keeping other points x j ( j = i) unchanged, and each minimization problem is relatively easier to solve by employing the technique of majorization. SMACOF has been widely used and the interested reader can refer to [6] for more references and to [54] for some critical comments on SMACOF when it is applied to the sensor network localization problem. The second and third group consist respectively the methods of SDP and EDM optimization. We will give a more detailed review on the two groups because of their close relevance to our proposed method in this paper. The main purpose of our review is to show that there lacks efficient numerical methods for the REE problem (3).

On SDP approach
We note that each of the four objective functions either involves the Euclidean distance d i j or its squared D i j = d 2 i j . A crucial observation is that constraints on them often have SDP relaxations. For example, it is easy to see where Y := X T X 0. Hence, the squared distance d 2 i j is a linear function of the positive semidefinite matrix Y . Consequently, the EDM cone D n can be described through linear transformations of positive semidefinite matrices. One can further relax the constraint Y = X T X to Y X T X . By the Schur-complement, one has By dropping the rank constraint, the robust MDS problem (4) can be relaxed to a SDP, which was initiated by Biswas and Ye [15]. For the Euclidean distance d i j , we introduce a new variable T i j = d i j . One may relax this constraint to T i j ≤ d i j , which has a SDP representation: Combination of (20), (21) and (22) leads to a large number of SDP relaxations. Typical examples, for the robust MDS problem (4), are the SDP relaxation method [5] and the edge-based SDP relaxation method [38,50] and [27], which leads to a comprehensive Matlab package SFSDP. For the squared stress (6), one may refer to [16,25]. For the stress problem (5), a typical SDP relaxation can be found in [35,Problem (8)]. However, unlike the problems (4), (5) and (6), the REE problem (3) does not have a straightforward SDP relaxation. We use an attempt made in [35] to illustrate this point below.
First, it is noted that problem (3) can be written in terms of EDM: The term | D i j − δ i j | is convex if δ i j > D i j and is concave otherwise. A major obstacle is how to efficiently deal with the concavity in the objective.
Secondly, by dropping the rank constraint and through certain linear approximation to the concave term, a SDP problem is proposed for (3) (see [35,Eq. (20)]): where the quantities a i j and b i j can be computed from δ i j . We note that each quadratic constraint in (23) is equivalent to a positive semidefinite constraint on S 2 + and D ∈ D n is a semidefinite constraint on S n + by (8). Therefore, the total number of the semidefinite constraints is |E| + 1. Finally, the optimal solution of (23) is then refined through a second-stage algorithm (see [35, Sect. IV(B)]). Both stages of the algorithmic scheme above would need sophisticated implementation skills and its numerical efficiency and solution quality are yet to be confirmed due to the twice relaxation schemes (dropping the rank constraint and linearization). The lack of efficient algorithms for (3) motivated our research in this paper.

On EDM approach
A distinguishing feature from the SDP approach is that this approach treats EDM D as the main variable, without having to rely on its SDP representation. This approach works because of the characterization (12) and that the orthogonal projection onto K n + has a closed-form formula [19,20]. Several methods are based on this formula. The basic model for this approach is the so-called nearest EDM problem: which is a convex relaxation of (6) with the special choice W i j ≡ 1. Here the elements of the matrix Δ (2) are given by Δ (2) i j := δ 2 i j . The relaxation is obtained by dropping the rank constraint rank(J D J ) ≤ r . Since the constraints of (24) are the intersection of a subspace and a convex cone, the method of alternation projection was proposed in [19,20] with applications to molecule conformation [21]. Newton's method for (24) was developed in [39]. Extensions of Newton's method for the model (24) with more constraints including general weights W i j , the rank constraint rank(J D J ) ≤ r or the box constraints (7) can be found in [3,11,40]. A recent application of the model (24) with a regularization term to Statistics is [55], where the problem is solved by an SDP, similar to that proposed by Toh [48].
There are two common features in this class of methods. One is that they require the objective function to be convex, which is true for the problems (4), (5) and (6) when formulated in EDM. The second feature is that the nonconvexity is only caused by the rank constraint. However, as already seen in Sect. 1.2, the REE problem (3) in terms of EDM has a nonconvex objective coupled with the distance d i j (not squared distances) being used. This has caused various difficulties in applying the existing EDM-based methods mentioned above to (3). A latest research [56] by the authors has tried to extend the EDM approach to the stress minimization problem (5) along a similar line as outlined in Sect. 1.2. Once again, we emphasize that the key difference between the problem (3) and (5) is about nonconvex objective vs convex objective and non-differentiability vs differentiability. Hence, the problem (3) is significantly more difficult to solve than (5). Nevertheless, we will show that it can be efficiently solved by the proposed EDM optimization.

Background and technical lemmas
In this part, we introduce the necessary background about subgradient and positive roots of a special depressed cubic equation. In particular, we will prove a technical result about a composite function between the absolute value and the square root functions. This result (Lemma 1) is in the style of Taylor-expansion for differentiable functions.

Subgradients of functions
An important function appearing in our EDM reformulation (11) of the REE problem (3) is φ δ (·) : + → + defined for a given constant δ > 0 by where + is the set of all nonnegative numbers. We will need to compute its subgradient in the sense of Rockafellar and Wets [42].
Definition 1 [42,Def. 8.3] Consider a function f : n → ∪ {−∞, +∞} and a pointx with f (x) finite. For a vector v ∈ n , one says that where the little 'o' term is a short-hand notation for the one-sided limit condition: We call ∂ f (x) the subdifferential of f atx. For a given number x ∈ , we define its sign by Apparently, φ δ (x) is continuous for x > 0 and its subdifferential at x > 0 is given by directly applying Def. 1 (note δ > 0) We note that the subdifferential of φ δ (x) at x = 0 is more complicated to describe. Fortunately, we won't need it in our analysis. We state our key lemma below.

Lemma 1
Let δ > 0 be given. It holds Compared with the definition of classical convex function, the inequality above has an extra negative quadratic term, which is not sufficient for φ δ (x) to be convex on (0, +∞) (actually, φ δ (x) is quasi-convex.) Fortunately, this property allows us to make the function φ δ + μϕ convex provided ϕ is a strongly convex and μ is chosen properly. And this is one of the key ideas developed in this paper.

Construction of the majorization function
A major building block in our algorithm is the majorization function g m (D, D k ) at a given point D k for the function g(A) defined in (14). We construct it below. Suppose A ∈ S n has the following eigenvalue-eigenvector decomposition: where λ 1 ≥ λ 2 ≥ · · · ≥ λ n are the eigenvalues of A in non-increasing order, and p i , i = 1, . . . , n are the corresponding orthonormal eigenvectors. We define a PCA-style matrix truncated at r : Recall the definition of Π B K n + (r ) (A) in (13). We let Π K n + (r ) (A) be an element of Π B K n + (r ) (A) and note that the function g(A) in (14) does not depend on the choice of As seen from the known results below, one particular element Π K n + (r ) (A) can be computed through PCA + r (A).

Lemma 2 For a given matrix A ∈ S n and an integer r ≤ n. The following results hold.
(i) [40,Eq. (22), Proposition 3.3] One particular Π K n + (r ) (A) can be computed through is well defined and is convex. Moreover, where ∂h(A) is the subdifferential of h(·) at A. (iv) [56, Lemma 2.2] Let g(A) be defined in (14). We have for any A ∈ S n Since h(·) is convex and Π K n This, with Lemma 2(iii), implies It is straightforward to check that the function g m (·, ·) in (31) satisfies the majorization properties (17).

Positive roots of depressed cubic equations
In our algorithm, we will encounter the positive root of a depressed cubic equation [7,Chp. 7], which arises from the optimality condition of the following problem where ν > 0 and t ∈ are given. A positive stationary point x must satisfy the optimality condition Let y := √ x. The optimality condition above becomes This is in the classical form of the so-called depressed cubic equation [7,Chp. 7]. Its roots (complex or real) and their computational formulae have a long history with fascinating and entertaining stories. A comprehensive revisit of this subject can be found in Xing [51] and a successful application of the depressed cubic equation to the compressed sensing can be found in [36,52]. The following lemma says that, under certain conditions, the Eq. (33) has two distinctive positive roots and its proof is a specialization of [9, Lemma 2.1(iii)] when p = 1/2 therein.

Penalized REE model and its' majorization subproblem
With the preparation above, we are ready to address our penalized REE problem (16) and its majorization subproblem (18). We first address the relationship between (16) and its original problem (11). We then show how the subproblem (18) is solved.

-optimal solution and critical point
The classical results on penalty methods in [34] on the differentiable case (i.e., all functions involved are differentiable) are not applicable here. Recently, the penalty approach was studied by Gao in her PhD thesis [18] in the context of semidefinite programming, which motivated our investigation below. The main result is that (16) provides an -optimal solution and -critical point for the original problem when the penalty parameter is above certain threshold. We first introduce the concept of critical point of (11) [i.e., (15)], which is a first order optimality condition. It is associated with the Lagrange function of (15) defined by where β ∈ is the Lagrangian multiplier and I B (D) is the indicator function taking value 0 if D ∈ B and +∞ otherwise. We say D ∈ S n is a critical point of (15) if it satisfies 0 ∈ ∂L(D,β) for someβ ∈ , namely, there exists Γ ∈ ∂ f (D) and where we used the following facts: [40,Proposition 3.4] and Conv(Ω) denotes the convex hull of a given set Ω) , and the sub-differential of the indicator function of a convex set Ω is its normal cone, namely, is enough for our analysis and hence it is used in our definition of approximate critical point.

Definition 2 (i) ( -Optimal solution)
Suppose D * is an optimal solution of (11). For a given error tolerance > 0, a point D is called an -optimal solution of (11) if it satisfies Obviously, if = 0, -optimal solution D would be an optimal solution of (11). We will show that the optimal solution of (16) is -optimal as well as an -critical point provided that ρ is large enough. Let D * ρ be an optimal solution of the penalized REE (16) and D r be any feasible solution of the original problem (11). If the lower bound matrix L ≡ 0, then we can simply choose D r = 0. Define We have the following result.

Proposition 1
For any ρ ≥ ρ , D * ρ must be an -optimal solution as well as ancritical point of (11).

Solving the subproblem
Having constructed the majorization function in (31), we now focus on how to solve the majorization subproblem (18), which is equivalent to the solution of the following problem. Given the current iterate Z ∈ B, the majorization subproblem aims to compute an improved iterate, denoted by Z + , by solving where the matrix Z K : . This subproblem has a perfect separability property that allows it to be computed elementwise: For the ease of our description, we denote the subproblem solution process by Here, PREEEDM stands for the Penalized REE by EDM optimization. We will show how (40) can be computed. Let us consider a simplified one-dimensional optimization problem, whose solution will eventually give rise to PREEEDM.
The acronym dcroot stands for the root of depressed cubic equation, which will eventually give rise to the solution formula of (41). It suffices to consider the case that matters to us: Before solving the above problem, we define some notation for convenience and Obviously, q(x) has a representation of two pieces: may not necessarily so. We will show that both pieces have a closed-form formula for their respective minimum.

Proposition 2
Consider the optimization problem: Then (43) has a unique solution x * − given by Now we characterize the optimal solution of q + (x) over B + .

Proposition 3
Assume that β < 4δ 3 and consider the optimization problem: Define Then q + (x) is strictly convex over the interval [δ 2 , ∞) and Putting together Propositions 2 and 3 gives rise to the optimal solution of (41). The optimal solution is either x * − or x * + , whichever gives a lower functional value of q(x). This is the first result of our major theorem below. We note that both Propositions 2 and 3 make use of the convexity of q − (x) and q + (x) on the respective interval [a, δ 2 ] and [δ 2 , b]. In fact, we can establish a stronger result that when the two pieces join together, the resulting function q(x) is still convex on the whole interval [a, b]. This result is very important to our convergence analysis in the next section and is the second result of the theorem below. A key tool for the proof is Lemma 1.
Then, the following hold.
(i) The optimal solution of the problem (41) is given by (iii) Let γ ω,β be defined in (42), as a function of ω. Suppose C > 0 is an arbitrarily given constant. Then there exists a constant κ > 0 such that However, its location could be within the interval [a, depending on the magnitudes of the parameters (ω, β and δ) involved. The dependence is illustrated in Fig. 1. We also note that the function q(x) may not be convex if the condition β < 4δ 3 is violated. It now follows from Theorem 1 that the optimal solution Z + i j in (39) can be computed by: Consequently, (40) is well defined and its elements can be computed by (47). q(x)

Algorithm PREEEDM and its convergence
With the preparations above, we are ready to state our algorithm. Let D k ∈ B be the current iterate. We update it by solving the majorization subproblem of the type (38) with Z replaced by D k : which can be computed by In more detail, we have , and the elements of D k+1 are computed as follows: Our algorithm PREEEDM is formally stated as follows.
A major obstacle in analysing the convergence for the penalized EDM model (16) is the non-differentiability of the objective function. We need the following two reasonable assumptions: Assumption 1 can be easily satisfied (e.g., setting the upper bound to be n 2 max{δ 2 i j }). Assumption 2 means that if δ i j = 0 (e.g., value missing), the corresponding weight W i j should be 0. This is a common practice in applications. If δ i j > 0, then we require δ 2 i j to be between L i j and U i j . We further define a quantity that bounds our penalty parameter ρ from below: Our first result in this section is about the boundedness of the subdifferential of f (·) along the generated sequence {D k }. (i) There exists a constant c 1 > 0 such that such that the first-order optimality condition for (49) is Theorem 1(i) ensures that D k i j > 0 for all k = 1, . . . , . Hence, we can apply Lemma 1 to each function φ δ i j (·) with x = D k+1 i j and y = D k i j . This yields for any ζ k+1 Multiplying W i j on both sides and adding those inequalities over (i, j), we get where Γ k+1 We note that the inequality (53) holds for any Γ k+1 ∈ ∂ f (D k+1 ).
Theorem 2 Let ρ > ρ o and {D k } be the sequence generated by Algorithm 1. Suppose Assumptions 1 and 2 hold.
That is, D is a critical point of the problem (16). Moreover, for a given > 0, if Proof (i) We are going to use the following facts that are stated on D k+1 and D k . The first fact is the identity: The second fact is due to the convexity of h(D) (see Lemma 2(ii)): The last fact is that there exists Γ k+1 ∈ ∂ f (D k+1 ) such that (52). Those facts yield the following chain of inequalities: This proves that the sequence {F ρ (D k )} is non-increasing and it is also bounded below by 0. Taking the limits on both sides yields D k+1 − D k → 0.
(ii) Suppose D is the limit of a subsequence {D k }, = 1, . . . ,. Since we have established in (i) that (D k +1 − D k ) → 0, the sequence {D k +1 } also converges to D. Furthermore, there exist a sequence of Γ k +1 ∈ ∂ f (D k +1 ) such that (52) holds. Proposition 4(ii) ensures that there exists a constant c 2 > 0 such that Γ k +1 ≤ c 2 for all k . Hence, there exists a subsequence of {k } (we still denote the subsequence by {k } for simplicity) such that Γ k +1 converges to some Γ ∈ ∂ f ( D). Now taking the limits on both sides of (52) on {k }, we reach the desired inequality (54). We now prove D is an -critical point of (11). Since we already have (54), we only need to show g( D) ≤ . It follows from D 0 ∈ K n (because of (i)) Taking the limit on the right-hand side yields where we used f ( D) ≥ 0. Therefore, thanks to ρ > ρ , it has (iii) We note that we have proved in (i) that (D k+1 − D k ) → 0. The convergence of the whole sequence to D follows from [26,Proposition 7].

Numerical experiments
In this part, we will conduct extensive numerical experiments of our algorithm PREEEDM by using MATLAB (R2014a) on a desktop of 8GB of memory and Inter(R) Core(TM) i5-4570 3.2Ghz CPU, against 6 leading solvers on the problems of sensor network localizations (SNL) in 2 (r = 2) and Molecular Conformation (MC) in 3 (r = 3). This section is split into the following parts. The MATLAB package is available at DOI: https://doi.org/10.5281/zenodo.3343047. Our implementation of PREEEDM was described in Sect. 6.1. We will give a brief explanation how the six benchmark methods were selected in Sect. 6.2. Descriptions of how the test data of SNL and MC were collected and generated, and extensive numerical comparisons are reported in Sect. 6.3.

Implementation
The PREEEDM Algorithm 1 is easy to implement. We first address the issue of its stopping criterion that is to be used in Step 3 of Algorithm 1. We monitor two quantities. One is on how close of the current iterate D k is to be Euclidean (belonging to K n + (r )). This can be computed by using (28) as follows.
where λ 1 ≥ λ 2 ≥ · · · ≥ λ n are the eigenvalues of (−J D k J ). The smaller Kprog k is, the closer D k is to K n + (r ). The benefit of using Kprog over g(D) is that the former is independent of any scaling of D.
The other quantity is to measure the progress in the functional values f ρ (·) by the current iterate D k . In theory (see Theorem 2), we should require ρ > ρ o , which is defined as (51) and is potentially large. As with the most penalty methods [34,Chp. 17], starting with a very large penalty parameter may degrade the performance of the method (e.g., causing ill-conditioning). We adopt a dynamic updating rule for ρ. In particular, we choose ρ 0 = κ max δ i j n 3/2 and update it as where and Ftol = ln(κ) × 10 −4 and Ktol = 10 −2 with κ being the number of non-zero elements of Δ. We terminate PREEEDM when  (27)). For the problem data input, Δ, L and U will be described below. For the initial point, we follow the popular choice used in [44,47] √ D 0 := Δ, where Δ is the matrix obtained by the shortest path distances among Δ. If Δ has no missing values, then Δ = Δ.

Benchmark methods
We select six representative state-of-the-art methods for comparison. They are ADMMSNL [37], ARAP [54], EVEDM (short for EepVecEDM) [12], PC [1], PPAS (short for PPA Semismooth) [24] and SFSDP [27]. Those methods have been shown to be capable of returning satisfactory localization/embedding in many applications. We will compare our method PREEEDM with ADMMSNL, ARAP, EVEDM, PC and SFSDP for Sensor Network Localization (SNL, r = 2) problems and with EVEDM, PC, PPAS and SFSDP for Molecular Conformation (MC, r = 3) problems since the current implementations of ADMMSNL, ARAP do not support the embedding for r ≥ 3.
We note that ADMMSNL is motivated by [45] and aims to enhance the package diskRelax of [45] for the SNL problems (r = 2). Both methods are based on the stress minimization (5). As we mentioned before, SMACOF [13,14] has been a very popular method for (5). However, we will not compare it with other methods here since its performance demonstrated in [54,56] was not very satisfactory (e.g., when comparing with ARAP) for either SNL or MC problems. To our best knowledge, PC is the only viable method, whose code is also publicly available for the model (3). We select SFSDP and PPAS because of their high reputation in the field of SDP and quadratic SDP in returning quality localizations and conformations. We note that SFSDP is for the model (4) and the methods PPAS and EVEDM are proposed for the model (6). It is worth mentioning that the MADMM package in [29] is capable of solving the Robust MDS (4) as well as other nonsmooth optimization problems. However, MADMM does not contain the implementation of its listed example: Robust MDS. So we were not able to compare it with ours here. We also implemented the subgradient method of Cayton and Dasgupta [8] for their robust Euclidean embedding. Numerical experiments showed that its performance was similar to PC on our tested problems. It works well when a large number of the dissimilarities in Δ are available and it often performs poorly otherwise. Hence, we omitted it from our reported results.
In our tests, we used all of their default parameters except one or two in order to achieve the best results. In particular, for PC, we terminate it when | f (D k−1 ) − f (D k )| < 10 −4 × f (D k ) and set its initial point to be the embedding by cMDS on Δ. For SFSDP which is a high-level MATLAB implementation of the SDP approach initiated in [50], we set pars.SDPsolver = "sedumi" because it returns the best overall performance, and pars.objSW = 1 when m > r + 1 and = 3 when m = 0. We also note that the parameter pars.minDegree controls the degree of a graph and thus enhances the strength of the SDP relaxation. Numerical experiments have shown that the larger it is, the more accurate solutions might be generated by SFSDP. However, the computational time shoots up dramatically when it increases even for small n.
Our extensive experiments suggest that its default value (pars.minDegree = r +2) is a balanced choice between solution quality and time of computation for large n. Hence we choose to use its default setting in our test. For ARAP, in order to speed up the termination, we let tol = 0.05 and IterNum = 20 to compute its local neighbour patches. Numerical performance demonstrated that ARAP could yield satisfactory embedding, but would take very long time for some examples with large n.

Numerical comparison
To assess the embedding quality, we adopt a widely used measure RMSD (Root of the Mean Squared Deviation) defined by , where x i 's are the true positions of the sensors or atoms in our test problems and x i 's are their corresponding estimates. The x i 's were obtained by applying the classical MDS (cMDS) method to the final output of the distance matrix, followed by aligning them to the existing anchors through the well-known Procrustes procedure (see [54], [6,Chp. 20] or [41, Proposition 4.1] for more details). Furthermore, upon obtaining x i 's, a heuristic gradient method can be applied to improve their accuracy and it is called the refinement step in [5]. We report rRMSD to highlight its contribution. As we will see, all tested methods benefit from this step, but with varying degrees.
The quality of the general performance of each method can be better appreciated through visualizing their key indicators: RMSD, rRMSD, rTime (time for the refinement step) and the CPU Time (in s) which is the total time including rTime. Hereafter, for all examples, we test 20 randomly generated instances for each case (n, m, R, nf) in SNL or each case (n, R, nf) in MC, and record the average results.

Comparison on SNL
SNL has been widely used to test the viability of many existing methods for the stress minimization. In such a problem, we typically have m anchors (e.g., sensors with known locations) and the rest sensors need to be located. We will test two types of SNL problems. One has a regular topological layout (Examples 1 and 2 below). The other has an irregular layout (Example 3).
where R is known as the radio range, i j 's are independent standard normal random variables, and nf is the noise factor (e.g., nf = 0.1 was used and it corresponds to 10% noise level). In literature (e.g., [5]), this type of perturbation in δ i j is known to be multiplicative and follows the unit-ball rule in defining N x and N a (see [3,Sect. 3.1] for more detail). The corresponding weight matrix W and the lower and upper bound matrices L and U are given as in the table below. Here, M is a large positive quantity. For example, M := n max i j Δ i j is the upper bound of the longest shortest path if the network is viewed as a graph.
Example 2 (Square Network with m random anchors) This example also tested in [5] is similar to Example 1 but with randomly generated anchors. The generation of n points follows the uniform distribution over the square region [−0.5, 0.5] 2 . Then the first m points are chosen to be anchors and the last (n − m) points to be sensors. The rest of the data generation is same as in Example 1.

Example 3 (EDM word network)
This problem has a non-regular topology and was first used in [3] to challenge existing methods. In this example, n points are randomly generated in a region whose shape is similar to the letters "E", "D" and "M". The ground truth network is depicted in Fig. 2. We choose the first m points to be the anchors. The rest of the data generation is same as in Example 1.
(a) Effect of the radio range R It is easy to see that the radio range R decides the amount of missing dissimilarities among all elements of Δ. The smaller R is, the  Fig. 3. It can be seen that ARAP and PREEEDM were joint winners in terms of both RMSD and rRMSD. However, the time used by ARAP was the longest. When R became bigger than 0.6, ADMMSNL, SFSDP and EVEDM produced similar rRMSD as ARAP and PREEEDM, while the time consumed by ADMMSNL was significantly larger than that by SFSDP, EVEDM and PREEEDM. By contrast, PC only worked well when R ≥ 1. Next we test a number of instances with larger size n ∈ {300, 500, 1000, 2000}. For Example 1, the average results were recorded in Table 1. When R = √ 2 under which no dissimilarities were missing because Example 1 was generated in a unit region, PC, ARAP and PREEEDM produced the better RMSD ( almost in the order of 10 −3 ). But with the refinement step, all methods led to similar rRMSD. This meant SFSDP and EVEDM benefited a lot from the refinement step. For the computational speed, PREEEDM outperformed others, followed by PC, EVEDM and SFSDP. By contrast, ARAP consumed too much time even for n = 500. When R = 0.2, the picture was significantly different since there were large amounts of unavailable dissimilarities in Δ. Basically, ADMMSNL, PC and SFSDP failed to localize even with the refinement due to undesirable RMSD and rRMSD (both in the order of 10 −1 ). Clearly, ARAP and PREEEDM produced the best RMSD and rRMSD, and EVEDM got comparable rRMSD but inaccurate RMSD. In terms of the computational speed, EVEDM and PREEEDM were very fast, consuming about 30 s to solve problems with n = 2000 nodes. By contrast, ARAP still was the slowest, followed by ADMMSNL and PC. Now we test those methods for the irregular network in Example 3. The average results were recorded in Table 2. We note that this example was generated in the region [0, 1] × [0, 0.5] as presented in Fig. 2. It implies that no dissimilarities in Δ were missing when R = √ 1.25 while a large number of dissimilarities in Δ were missing when R = 0.1. When R = √ 1.25, it can be clearly seen that SFSDP and EVEDM failed to localize before the refinement step due to their large RMSD (in the order of 10 −1 ), whilst the rest four methods succeeded. However, they all achieved a  Finally we test a number of instances of Example 2 with choosing n ∈ {300, 500, 1000, 2000} and m ∈ {10, 50}. The average results were recorded in Table 3. When m = 10, ADMMSNL and PC produced undesirable RMSD and rRMSD (both in the order of 10 −1 ). SFSDP benefited greatly from the refinement because it generated relatively inaccurate RMSD. By contrast the rest three methods enjoyed the successful localization except for EVEDM under the case n = 300. With regard to the computational speed, EVEDM and PREEEDM were the fastest, followed by SFSDP, PC, ADMMSNL and ARAP. When m = 50, more information was known, the results were better than before, especially for the methods ADMMSNL and PC. But PC still heavily relied on the refinement step to get the satisfactory localization. The rest five methods produced a satisfactory localization with varying degree of accuracy. It is encouraging to see that PREEEDM produced the most accurate rRMSD for all cases. The comparison of the computational speed is similar to the case of m = 10. We repeated the test for Example 3 and the average results were recorded in Table 4, where we observed a similar performance of the six methods as for Example 2. We omit the details.  Fig. 6, in terms of RMSD it can be seen that ARAP got the smallest ones, whilst EVEDM and PC obtained the worst ones. The line of ADMMSNL dropped down from 0.1 ≤ nf ≤ 0.3 and then ascended. By contrast the line of PREEEDM reached the peak at nf = 0.3 but declined afterwards and gradually approached to RMSD of ARAP. However, after the refinement step, ARAP, SFSDP and PREEEDM all derived a similar rRMSD while the other three methods produced undesirable ones. Apparently, EVEDM was indeed the fastest (yet with the worst rRMSD), followed by PC, SFSDP and PREEEDM. Again, ARAP and ADMMSNL were quite slow. Next, we test Example 2 with a moderate size (for the visualization purpose in Fig. 7) n = 200, m = 4 and R = 0.3 and varying nf ∈ {0.1, 0.3, 0.5}. The actual embedding by each method was shown in Fig. 7, where the four anchors were plotted in green square and x i in pink points were jointed to its ground truth location (blue circle). It can be clearly seen that ARAP and PREEEDM were quite robust to the noise factor since their localization matched the ground truth well. EVEDM failed to locate when nf = 0.5. By contrast, SFSDP generated worse results when nf got bigger, and ADMMSNL and PC failed to localize for all cases. Finally, we test Example 1 with larger sizes n ∈ {300, 500, 1000, 2000} and fixed m = 4, R = 0.3. The average results were recorded in Table 5. When nf = 0.1, ADMMSNL and PC failed to render accurate embedding. Compared with ARAP, EVEDM and PREEEDM, SFSDP generated lager RMSD and rRMSD. Again, EVEDM and PREEEDM ran faster than ARAP. When nf = 0.7, the results were different. ARAP and PREEEDM were still able to produce high-quality RMSD and rRMSD. However, the former took extremely long time (16617 vs. 83 s). By contrast, ADMMSNL and PC again failed to reconstruct the network. Furthermore, EVEDM got large RMSD but comparable rRMSD when n ≤ 1000, but it failed when n = 2000.

Comparison on MC
MC has long been an important application of EDM optimization [2,21,33]. We will test two types of MCs respectively from an artificial data set and a real data set in Protein Data Bank (PDB) [4]. For the former, we adopt the rule of generating data from [2,33]. For the latter, we used the real data of 12 molecules derived from 12 structures of proteins from PDB. They are 1GM2, 304D, 1PBM, 2MSJ, 1AU6, 1LFB, 104D, 1PHT, 1POA, 1AX8, 1RGS, 2CLJ. They provide a good set of test problems in terms of the size n, which ranges from a few hundreds to a few thousands (the smallest n = 166 for 1GM and the largest n = 4189 for 2CLJ). The distance information was obtained in a realistic way as done in [24]. Example 4 (Artificial data) As described in [2,33], the artificial molecule has n = s 3 atoms (x 1 , . . . , x n ) located in the three-dimensional lattice for some integer s ≥ 1, i.e., x i = (i 1 , i 2 , i 3 ) T . We define N x for the index set on which δ i j are available as: where p(x i ) := 1 + (1, s, s 2 ) T x i = 1 + i 1 + si 2 + s 2 i 3 and R is a given constant (e.g., R = s 2 ). The corresponding dissimilarity matrix Δ, weight matrix W and the lower and upper bound matrices L and U are given as in the table below. Here the generation of δ i j is the same as Example 1.
Example 5 (Real PDB data) Each molecule comprises n atoms {x 1 , . . . x n } in 3 and its distance information is collected as follows. If the Euclidean distance between two of the atoms is less than R, the distance is chosen; otherwise no distance information about this pair is known. For example, R = 6Å (1Å = 10 −8 cm) is nearly the maximal distance that the nuclear magnetic resonance (NMR) experiment can measure between two atoms. For realistic molecular conformation problems, not all the distances below R are known from NMR experiments, so one may obtain c% (e.g., c = 50%) of all the distances below R. Denote N x the set formed by indices of those measured distances. Moreover, the distances in N x can not be exactly measured. Instead, only lower bounds i j and upper bounds u i j are provided, that is for where i j , ε i j ∼ N (0, nf 2 × π/2) are independent normal random variables. In our test, we set the noise factor nf = 0.1 and the parameters W , Δ, L, U ∈ S n are given as in the table below, where M > 0 is the upper bound (e.g., M := n max i j Δ i j ).
As we mentioned before, the current implementations of ADMMSNL, ARAP do not support the embedding for r ≥ 3 and thus are removed in the following comparison, where the method PPAS will be added. The main reason for adding PPAS is that it is particularly suitable and credible for the MC problems [24,25].
(d) Test on Example 4. To see the performance of each method on this problem, we first test it with fixing s = 6 (n = 6 3 ), nf = 0.1 but varying R ∈ {36, 38, . . . , 48}. We note that the percentage of available dissimilarities increased from 32.47 to 39.87% with R increasing from 36 to 48, making the problem become 'easier' for conformation. The Average results were recorded in Fig. 8. Clearly, PREEEDM and PPAS outperformed the other three methods in terms of RMSD and rRMSD. The former generated the best RMSD when R ≥ 42 while the latter got the best RMSD when R ≤ 42, but they both obtained similar rRMSD. As for the computational speed, PREEEDM ran much faster than PPAS. By contrast, the other three methods failed to produce accurate embeddings due to the worse RMSD and rRMSD. Notice that the refinement would not always make the final results better. For instance, rRMSD yielded by SFSDP was bigger than RMSD for each s.
We then test the example with fixing s = 6 (n = 6 3 ), R = s 2 and varying nf ∈ {0.1, 0.2, . . . , 0.5}. As illustrated in Fig. 9, in terms of RMSD and rRMSD, it can be clearly seen that PREEEDM and PPAS were the joint winners. In particular, our method rendered the best RMSD when nf ≥ 0.2 and also ran much faster than PPAS. Obviously, the other three methods again failed to obtain desirable RMSD and rRMSD irrelevant of the time they used. Finally, for larger size problems with n = s 3 and s ∈ {7, 8, . . . , 13}, the average results were presented in Fig. 10, where we omitted the results by PPAS for s > 10 because it took too much time to terminate. It is worth mentioning that the percentage of the available dissimilarities over all elements of Δ decreases from 26.78% to 14.83% when s increasing from 7 to 13, making the problems more and more challenging. Clearly PC, SFSDP and EVEDM failed to locate all atoms in 3 . PPAS rendered the most accurate RMSD when s ≤ 10 whilst PREEEDM achieved the most accurate RMSD when s > 10 and the most accurate rRMSD for all cases. Equally important for PREEEDM is that it spent less than 50 s for all tested cases, while PPAS took much more time to terminate (e.g., consuming over 2000 s when s ≥ 10).
(e) Test on Example 5 For the 12 collected real data, we fixed R = 6, c = 50% and nf = 0.1. The generated embeddings by the five methods for the three molecules 1GM2, 1AU6 and 1LFB were shown in Fig. 11, where the true and estimated positions of the atoms were plotted by blue circles and pink stars respectively. Each pink star was linked to its corresponding blue circle by a pink line. For these three data, PREEEDM and PPAS almost conformed the shape of the original data. Clearly, the other three methods failed to conform. The complete numerical results for the 12 problems were reported in Table 6. It can be clearly seen that PREEEDM and PPAS performed significantly better in terms of the RMSD and rRMSD than the other methods. What is more impressive is that PREEEDM only used a small fraction of the time by PPAS, which in general took relatively long time to terminate. For example, PREEEDM only used 22.64 s for 2CLJ, which is a very large data set with n = 4189. In contrast, we had to omit the result of PPAS for this instance (as well as to omit for other tested instances, and the missed results were indicated as "-" in Table 6) because it took too long to terminate.

Robustness of PREEEDM
The excellent performance of PREEEDM reported above was actually due to its robustness to noise. Previous examples all had Gaussian noise. We now demonstrate below that PREEEDM works much better than the other methods when the noise is from a heavy-tailed distribution, for instance, t-distribution with a small degree of freedom.
We also take this opportunity to test SQREDM solver of our own [56], which also made use of penalty, majorization and minimization techniques, yet for the least squares problem (5). We will see that PREEEDM outperforms SQREDM for both types of noise (Gaussian and t distributions).
To shorten the presentation, we restrict our numerical tests on two representative examples: Example 1 with n = 100, R = 0.3 and Example 4 with s = 5, R = s 2 . For each example, we generate 20 instances under two types of noise from standard normal distribution and Student-t distribution with the degree of freedom being 1. We alter nf from {0.1, 0.2, . . . , 0.9} and from {0.01, 0.02, . . . , 0.09} for the Gaussian and the Student-t noises respectively. Average RMSD were recorded in Fig. 12. We have the following observations.
(i) PREEEDM is competitive under Gaussian noise. For Example 1, Fig. 12a showed ARAP yielded the best RMSD followed by PREEEDM and SQREDM. For Example 4, Fig. 12c showed that PREEEDM rendered the smallest RMSD for most cases, followed by PPAS and SQREDM (note that the current implementation of ARAP is only for r = 2 and hence is not applicable to this example).
In particular, when the nf is over 50%, PPAS and PREEEDM closely follow each other. The behaviour of those methods under Gaussian noise is expected as the least-squares formulation is equivalent to the maximum-likelihood criterion. On the one hand, least squares favour large distances. On the other hand, under Gaussian (thin-tailed), the number of large distance errors is relatively small and hence would not cause significant distortion in locating the unknown sensors. (ii) PREEEDM performs the best under heavy-tailed noise (from Student t 1 distribution). For Example 1, both PREEEDM and SQREDM behaved much better than the other methods, see Fig. 12b. For Example 4, PREEEDM stood out as the best method when nf is bigger than 0.02 and is much better than SQREDM, see Fig. 12d. The test data now has more numbers of large distance errors than under the Gaussian distribution and the absolute value criterion alleviates the tendency   We conclude that PREEEDM based on the model (3) is robust to the noise in terms of these two examples.

Conclusion
The purpose of this paper is to develop an efficient method for one of the most challenging distance embedding problems in a low-dimensional space, which have been widely studied under the framework of multi-dimensional scaling. The problem employs 1 norm to quantify the embedding errors. Hence, the resulting model (3) appears to be robust to outliers and is referred to as the robust Euclidean embedding (REE) model. To the best knowledge of the authors, the only viable method, whose matlab code is also publicly available for REE is the PlaceCenter (PC) algorithm proposed in [1]. Our extensive numerical results on the SNL and MC test problems convincingly demonstrated that the proposed PREEEDM method outperform PC in terms of both the embedding quality and the CPU time. Moreover, PREEEDM is also comparable to several state-of-the-art methods for other embedding models in terms of the embedding quality, but is far more efficient in terms of the CPU time. The advantage becomes even more superior as the size of the problem gets bigger.
The novelty of the proposed PREEEDM lies with its creative use of the Euclidean distance matrix and a computationally efficient majorization technique to derive its subproblem, which has a closed-form solution closely related to the positive root of the classical depressed cubic equation. Furthermore, a great deal of effort has been devoted to its convergence analysis, which well justifies the numerical performance of PREEEDM. We feel that PREEEDM will become a very competitive embedding method in the field of SNL and MC and expect its wide use in other visualization problems.
For y ≥ δ 2 , we have the following chain of inequalities Hence, we proved the conclusion for this case.

Proof of Theorem 1
Proof (i) is a direct consequence of Propositions 2 and 3. We now prove (ii). For any x, y > 0 and any ξ x ∈ ∂q(x), it follows that where the first inequality above used Lemma 1 and the last inequality used the fact 4δ 3 > β > 0. Swapping the role of x and y above yields q(x) − q(y) > ξ y (x − y) ∀ x, y > 0, ξ y ∈ ∂q(y).
This together with Theorem 12.17 of [42] proves that q(x) is strictly convex over [a, b]. The rest in (ii) is just the first order optimality condition of the convex optimization problem (41) because we just proved the convexity of q(x) over [a, b]. Finally, we prove (iii). It follows from (42) that and from Propositions 2 and 3 that x * − ≥ min{δ 2 , 1, κ 0 } and x * + ≥ δ 2 .
We finish our proof.

Proof of Proposition 4
Proof (i) Let us pick a pair (i, j) such that W i j > 0, which implies δ i j > 0 (Assumption 2). It follows from (50) that where Z k−1 K := −Π K n + (r ) (−D k−1 ). Since B is bounded (Assumption 1) and D k ∈ B, the sequence {D k } is bounded. Lemma 2 implies We let ∂ i j f (D) denote the subdifferential of f with respect to its (i, j)th element D i j . We consider two cases. Case 1: W i j = 0. This implies that f (D) is a constant function (≡ 0) of D i j and hence f (D) is continuously differentiable with respect to D i j . Consequently, ∂ i j f (D k ) = {0}.
Case 2: W i j > 0, which implies δ i j > 0 (Assumption 2). It follows from (i) that there exists c 1 > 0 such that D k i j ≥ c 1 for all k = 1, 2, . . . , . The Eqs. (64) and (25) yield which implies that for any ξ k i j ∈ ∂ i j f (D k ) there exists ζ k i j ∈ sign((D k i j ) 1/2 − δ i j ) such that |ξ k i j | = W i j |ζ k i j | 2 D k i j ≤ W i j / 4c 1 .
In other words, ∂ i j f (D k ) is bounded by W i j / √ c 1 , which is independent of the index k. It follows directly from the definition of subdifferential [42,Chp. 8.3] that in the sense that for any Γ k ∈ ∂ f (D k ), there exist ξ k i j ∈ ∂ i j f (D k ) such that Γ k i j = ξ k i j , i, j = 1, . . . , n.
(iii) Since ρ > ρ o , for each pair (i, j) we have β i j := W i j /ρ < 4δ 3 i j . It then follows from Theorem 1(ii) that each separable function f k i j (D i j ) is convex and hence the function f k ρ (D) is convex over D ∈ B. Consequently, subproblem (48) is convex. The first-order necessary and sufficient optimality condition is just (52).