Distance geometry in linearizable norms

Publication GSI2017
OAI : oai:www.see.asso.fr:17410:22560
contenu protégé  Document accessible sous conditions - vous devez vous connecter ou vous enregistrer pour accéder à ou acquérir ce document.
- Accès libre pour les ayants-droit


Distance Geometry puts the concept of distance at its center. The basic problem in distance geometry could be described as drawing an edge-weighted undirected graph in RK for some given K such that the positions for adjacent vertices have distance which is equal to the corresponding edge weight. There appears to be a lack of exact methods in this eld using any other norm but ℓ2. In this paper we move some rst steps using the ℓ1 and ℓ norms: we discuss worst-case complexity, propose mixed-integer linear programming formulations, and sketch a few heuristic ideas.

Distance geometry in linearizable norms


application/pdf Distance geometry in linearizable norms Claudia D'Ambrosio, Leo Liberti
Détails de l'article
contenu protégé  Document accessible sous conditions - vous devez vous connecter ou vous enregistrer pour accéder à ou acquérir ce document.
- Accès libre pour les ayants-droit

Distance geometry in linearizable norms


Voir la vidéo


295.14 Ko


Creative Commons Aucune (Tous droits réservés)


Sponsors Platine


Sponsors Bronze


Sponsors scientifique





<resource  xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
                xsi:schemaLocation="http://datacite.org/schema/kernel-4 http://schema.datacite.org/meta/kernel-4/metadata.xsd">
        <identifier identifierType="DOI">10.23723/17410/22560</identifier><creators><creator><creatorName>Leo Liberti</creatorName></creator><creator><creatorName>Claudia D'Ambrosio</creatorName></creator></creators><titles>
            <title>Distance geometry in linearizable norms</title></titles>
        <resourceType resourceTypeGeneral="Text">Text</resourceType><subjects><subject>distance geometry</subject><subject>norms</subject><subject>mathematical programming</subject></subjects><dates>
	    <date dateType="Created">Wed 7 Mar 2018</date>
	    <date dateType="Updated">Wed 7 Mar 2018</date>
            <date dateType="Submitted">Fri 20 Apr 2018</date>
	    <alternateIdentifier alternateIdentifierType="bitstream">2916b3754046423534e0a79d4be579b03935951b</alternateIdentifier>
            <description descriptionType="Abstract">Distance Geometry puts the concept of distance at its center. The basic problem in distance geometry could be described as drawing an edge-weighted undirected graph in RK for some given K such that the positions for adjacent vertices have distance which is equal to the corresponding edge weight. There appears to be a lack of exact methods in this eld using any other norm but ℓ2. In this paper we move some rst steps using the ℓ1 and ℓ∞ norms: we discuss worst-case complexity, propose mixed-integer linear programming formulations, and sketch a few heuristic ideas.

Distance geometry in linearizable norms Claudia D’Ambrosio1 and Leo Liberti1? LIX CNRS (UMR7161), École Polytechnique, 91128 Palaiseau, France {dambrosio,liberti}@lix.polytechnique.fr Abstract. Distance Geometry puts the concept of distance at its center. The basic problem in distance geometry could be described as drawing an edge-weighted undirected graph in RK for some given K such that the positions for adjacent vertices have distance which is equal to the corresponding edge weight. There appears to be a lack of exact methods in this field using any other norm but `2. In this paper we move some first steps using the `1 and `∞ norms: we discuss worst-case complexity, propose mixed-integer linear programming formulations, and sketch a few heuristic ideas. Keywords: distance geometry, norms, mathematical programming 1 Introduction We discuss the following basic problem in Distance Geometry (DG) Distance Geometry Problem (DGP). Given an integer K > 0 and a simple, edge-weighted, undirected graph G = (V, E, d), where d : E → R+, determine whether there exists a realization function x : V → RK such that: ∀{i, j} ∈ E kxi − xjk = dij, (1) where k · k is either the `1 or the `∞ norm. We assume all along that, without loss of generality, G is connected, otherwise it suffices to realize the disconnected components independently. Most existing work concerning the DGP focuses on the `2 (or Euclidean) norm, the only exception being [4] on `∞ norm DGPs. In this paper, we move some steps forward in the direction of the `1 and `∞ norms, which we call linearizable norms, since their unit spheres are polyhedral. DG in the `2 norm recently received a lot of attention [14, 7] due to its widespread use in engineering and science applications, such as, for example, finding the structure of proteins from Nuclear Magnetic Resonance (NMR) inter- atomic distances [20] and many others. It was shown in [18] that the DGP (with the Euclidean norm and K = 1) is NP-hard, by reduction from Partition. The limited attention to DG in other norms stems from a scarcity of ap- plications. Yet, recently, we were made aware of applications for both of the ? Partly supported by the ANR “Bip:Bip” project under contract ANR-10-BINF-0003. linearizable norms. The DGP with the `1 norm arises in the positioning of mo- bile sensors in urban areas [2]. The DGP with the `∞ or with the `1 can be used in order to fill a hypercube with a pre-fixed number of “well-distributed” points, which is relevant in the design of experiments [5, 19]. The rest of this paper is organized as follows. We recall some notions relating to linearizable norms in Sect. 2. We propose two new formulations for linearizable norms in Sect. 3, and prove that the DGP in linearizable norms is NP-complete for any K. Lastly, we sketch some new ideas for solving the linearizable norm DGP using heuristics in Sect. 4. 2 Known results for `1 and `∞ norms Complexity: Since the `1, `2, `∞ norms coincide for K = 1, the reduction from the Partition problem given in [18] shows that the DGP is NP-complete for these three norms in K = 1. Since a realization in one dimension can be em- bedded isometrically in any number of dimensions, this shows by inclusion that the DGP is NP-hard for these three norms. We show in Thm. 1 that the `1 and `∞ norm variants are also NP-complete for K > 1. This strikes a remarkable difference with the Euclidean norm DGP, for which the status of membership in NP is currently unknown [1]. Isometric embeddings: For the `∞ norm, the isometric embedding problem can be solved in a very elegant way [10]: any finite metric (X, d) with X = {x1, . . . , xn} can be embedded in Rn using the `∞ norm by means of the Fréchet embedding: ∀i, j ≤ n T(xi) = (d(xi, x1), . . . , d(xi, xn)). The proof is short and to the point: for any i, j ≤ n we have: kT(xi) − T(xj)k∞ = max k≤n |d(xi, xk) − d(xj, xk)| ≤ max k≤n d(xi, xj) = d(xi, xj) by the triangle inequality on the given metric d. Moreover, the maximum over k of |d(xi, xk) − d(xj, xk)| is obviously achieved when k ∈ {i, j}, in which case |d(xi, xk) − d(xj, xk)| = d(xi, xj). So we have kT(xi) − T(xj)k∞ = d(xi, xj) as claimed. We remark that (X, d) need not be given explicitly: the distance matrix is enough. We also remark that a Fréchet embedding can be constructed for any given square symmetric matrix A: if the Fréchet embedding of A is infeasible w.r.t. A, it means that A is not a valid distance matrix. For the `1 norm there no such general result is known. It is known that `2 metric spaces consisting n points can be embedded in a vector space of O(n) dimensions in `1 norm almost isometrically [16, §2.5] (the “almost” refers to a multiplicative distortion measure of the form: (1 − ε)kxk2 ≤ kT(x)k1 ≤ (1 + ε)kxk2 for some ε ∈ (0, 1)). 3 MILP formulations for linearizable norms The DGP is a nonlinear feasibility problem. As such, it can be modelled by means of Mathematical Programming (MP), which is a formal language for describing optimization problems. When the norm is linearizable, it is possible to replace nonlinear functions by piecewise linear forms, which are routinely modelled in MP using binary variables and linear forms. This yields MPs of the Mixed-Integer Linear Programming (MILP) class. A convenient feature of MILP is that solution technology is very advanced — MILP solvers are currently at the forefront for their generality and empirical efficiency. To the best of our knowledge, no MILP formulations have ever been proposed for the DGP in linearizable norms. We give here two MILP formulations for `1 and `∞ norms. We first re-write Eq. (1) as follows: min x P {i,j}∈E | kxi − xjkp − dij |, for p ∈ {1, ∞}. Obviously, even if the unconstrained optimization problem above always has a feasible solution, it has global optimal value zero if and only if the global optimum is a solution of Eq. (1). The MILP formulations below can be solved using any off-the-shelf MILP solver, such as, e.g., CPLEX [13]. The `1 norm. For p = 1 we write: min x X {i,j}∈E X k≤K |xik − xjk| − dij . (2) The MILP reformulation we propose is the following: min x,s,t,z P {i,j}∈E sij ∀{i, j} ∈ E −sij ≤ P k≤K (t+ ijk + t− ijk) − dij ≤ sij ∀k ≤ K, {i, j} ∈ E t+ ijk − t− ijk = xik − xjk ∀k ≤ K, {i, j} ∈ E t+ ijk ≤ dijzijk ∀k ≤ K, {i, j} ∈ E t− ijk ≤ dij(1 − zijk) ∀k ≤ K P i∈V xik = 0 ∀{i, j} ∈ E sij ∈ [0, dij] ∀k ≤ K, {i, j} ∈ E t+ ijk, t− ijk ∈ [0, dij] ∀k ≤ K, {i, j} ∈ E zijk ∈ {0, 1}.                                    (3) Additional variables sij ≥ 0, for each {i, j} ∈ E, are non-negative decision variables sij ≥ 0 that represent the outermost absolute value of (2) thanks to the well known fact that min |f| is equivalent to min ˆ f≥0 ˆ f subject to − ˆ f ≤ f ≤ ˆ f. Moreover, additional slack and surplus variables t+ , t− were introduced to reformulate the innermost absolute value terms. In order to do this, they are subject to complementarity constraints t+ ijk t− ijk = 0, ∀k ≤ K, {i, j} ∈ E, that in (3) were linearized by adding binary variables z and the two sets of standard “big-M” constraints that links the t and the z variables. Note that, as each realization can be translated at will, constraint P i∈V xi = 0 can be safely added. It means that we can, without loss of generality, impose that the barycenter is zero and it is not necessary but useful in practice, see [6] for detailed empirical results. Again for practical efficiency, we let U = P {i,j}∈E dij and use it to bound x, so that for each i ∈ V and k ≤ K, we have xik ∈ [−U, U]. Proposition 1. Eq. (3) is a valid formulation for the DGP using the `1 norm. We omit the proof for lack of space. The `∞ norm. For p = ∞ we write: min x P {i,j}∈E max k≤K |xik − xjk| − dij . The MILP reformulation that we propose is the following: min x,w,z,s,t P {i,j}∈E sij ∀{i, j} ∈ E, k ∈ K t+ ijk + t− ijk − dij ≤ sij ∀{i, j} ∈ E, k ∈ K t+ ijk + t− ijk + sij ≥ dijwijk ∀{i, j} ∈ E P k≤K wijk ≥ 1 ∀k ≤ K, {i, j} ∈ E t+ ijk − t− ijk = xik − xjk ∀k ≤ K, {i, j} ∈ E t+ ijk ≤ dijzijk ∀k ≤ K, {i, j} ∈ E t− ijk ≤ dij(1 − zijk) ∀k ≤ K P i∈V xik = 0 ∀{i, j} ∈ E sij ∈ [0, dij] ∀k ≤ K, {i, j} ∈ E t+ ijk, t− ijk ∈ [0, dij] ∀k ≤ K, {i, j} ∈ E wijk, zijk ∈ {0, 1}                                              (4) where the outermost and the innermost absolute values are reformulated as before and the barycenter constraint is added. At this point similarities with the `1 norm stop. Note that the first three set of constraints model the linearization of constraints −sij ≤ max k≤K (t+ ijk +t− ijk)−dij ≤ sij, ∀{i, j} ∈ E. The center and right-hand-side can be reformulated as the first set of constraints. The left-hand-side and the center is more complicated as it corresponds to a non-convex constraint. For linearizing it, we need to add auxiliary binary variables w and the second and third sets of constraints that express the fact that at least one out of K components satisfies the constraint (that component being the maximum, of course). Proposition 2. Eq. (4) is a valid formulation for the DGP using the `∞ norm. We omit the proof for lack of space. The DGP with linearizable norms is in NP. Our MILP formulations make it easy to prove the following. Theorem 1. The DGP in linearizable norms is NP-complete. Proof. It was already remarked in Sect. 2 that the DGP in linearizable norms is NP-hard by reduction from Partition to the DGP in K = 1, where the three norm `1, `2, `∞ coincide. Now by Prop. 1 and 2 we know that certificates to the DGP in linearizable norms are rational for rational input, since the MP formu- lations in (3) and (4) can be solved by a Branch-and-Bound (BB) algorithm, each node of which involves the solution of a Linear Program (LP), which is known to be in NP. Simple BB implementations will continue branching until the incumbent, found by solving the LP at some node, is proven optimal. So the certificate (solution) provided by the BB is polynomially sized, as claimed. u t 3.1 Computational results Both formulations (Eq. (3) and (4)) belong to the MILP class, and can be solved by several existing MILP solvers. We employ CPLEX 12.6.2 [13] on a MacBook Pro mid-2015 running Darwin 15.5.0 on a (virtual) quad-core i7 at 3.1GHz with 16 GB RAM. We generated a set of (feasible) random DGP instances in both norms as follows: for each cardinality value n = |V | of the vertex set of G ranging over {10, 15, 20, 25, 30, 35, 40} we sampled n points in the plane, bounded by a box. We then added a Hamiltonian cycle to the edge set, in order to guarantee con- nectedness. Lastly, we add the remaining edges with an Erdős-Rényi generation process [8] having probability s ranging over {0.1, 0.2, 0.3, 0.5, 0.8}. This yields 35 instances per norm p ∈ {1, ∞}. We deployed CPLEX on these 70 instances with a time limit of 10 minutes enforced on the “wall clock”, meaning the actual elapsed time. We used this measure instead of the user time since CPLEX exploits all four processor cores, which means that the system time taken for parallel execution tasks is essential. Since the running time is limited, we do not always find feasible solutions. We evaluate the error by employing two well known measures: the scaled Mean Distance Error (MDE) and the scaled Largest Distance Error (LDE). MDE(x) = 1 |E| X {i,j}∈E |kxi − xjkp − dij| dij , LDE(x) = max {i,j}∈E |kxi − xjkp − dij| dij . Intuitively, the scaled MDE gives an idea of the percentual average discrepancy of the given realization from the given data. The scaled LDE gives an idea of the percentual average worst error over all edges. Obviously, the LDE is generally higher than the MDE. The detailed results are not reported for lack of space. For each instance we computed scaled MDE and LDE scores, and compared the wall clock CPU time taken by CPLEX running on Eq. (3) and Eq. (4). From these results, it appears clear that DGP instances for linearizable norms can be solved in practice up to at least n = 40, though the densest instances are more difficult with the `∞ norm. A good `1 norm solution (with MDE and LDE smaller than 10−4 ) has been obtained using CPLEX on Eq. (3) for a 50 vertex instance with s = 0.8 in just over 304s of wall clock time, but no solution was obtained within 10 minutes of wall clock time for a 60 vertex instance and s = 0.8. This seems to indicate the need for heuristic methods (discussed below) to address larger instance sizes. 4 Heuristic ideas In this section we provide some ideas to design heuristic methods for the DGP problem with `∞ and `1 norm and report computational results. 4.1 Solving the DMCP in the `∞ norm A problem related to DGP, the Distance Matrix Completion Problem (DMCP), asks whether a partially defined matrix can be completed to a (full) distance matrix in a specified norm. Although a solution of the DMCP is a distance matrix, finding the missing distance values usually entails finding a realization that is consistent with the given values: this realization is then used to compute the missing entries in the given partial matrix. In this sense, realizations provide certificates for both DGP and DMCP. In the DMCP, however, differently from the DGP, the dimensionality K of the embedding space is not given — realizations into RK for any K > 0 will provide a feasibility certificate for DMCP instances. In the following we propose a heuristic for DMCP that will be the starting point for a heuristic for DGP. Our first proposal concerns the exploitation of the Fréchet embedding to “approximately solve” the DMCP (rather than the DGP) in the `∞ norm. Our algorithm works as follows: 1. let A0 be the n×n weighted adjacency matrix of G, where off-diagonal zeros are to be considered as “unspecified entries” 2. complete A0 to a full symmetric matrix A using the Floyd-Warshall all- shortest-paths algorithm [17] 3. output the realization x : V → Rn given by the Fréchet embedding xi = Ai for each i ∈ V . By “approximately solve” we mean that the output realization x is generally not going to be a valid certificate for the given DMCP instance, but that its MDE and LDE error measures are hopefully going to be low enough. The worst- case complexity of this heuristic is dominated by the Floyd-Warshall algorithm, which is generally O(n3 ). Experimentally, we found that this heuristic is useless for sparse graphs (where the reconstruction of A0 has the highest chances of being wrong), but is both fast and successful for dense graphs, notably those with s = 0.8, which represented the hardest instances in the tests. This heuristic took 2.42s to find all realizations for 17 random graphs obtained with s = 0.8 and n varying in {5 + 5` | 1 ≤ ` ≤ 17}. The cumulative scaled MDE over all instances is 0.019, with cumulative scaled LDE at 9.61 mostly due to one bad outlier. Overall, it found MDE and LDE smaller than 10−4 for 11 over 17 instances. In conclusion, it is both fast and effective. 4.2 Solving the DGP in the `∞ norm The second idea turns the above heuristic into a method for solving DGP in `∞ norm for a given dimensionality K: it consists in selecting the K columns from the realization x ∈ Rn×n obtained in Sect. 4.1 which best match the given distances. For this purpose, we solve the following MILP (based on Eq. (4)): min w,y,s P {i,j}∈E sij ∀{i, j} ∈ E, k ≤ n |xik − xjk|yk − sij ≤ dij ∀{i, j} ∈ E, k ≤ n |xik − xjk|yk + sij ≥ dijwijk ∀{i, j} ∈ E P k≤K wijk ≥ 1 P k≤n yk = K ∀{i, j} ∈ E sij ∈ [0, dij] ∀k ≤ K, {i, j} ∈ E wijk, yk ∈ {0, 1}.                            (5) In Eq. (5), note that x are no longer decision variables, but constants (output of the Fréchet heuristic for the DMCP). The new decision variables y ∈ {0, 1}n decide whether coordinate k should be in the realization in RK or not. Since solving MILP takes exponential time, and we would like the heuristic to be fast, we set a 120s time limit on CPLEX. For instances up to n = 70 with s fixed at 0.8, the average scaled MDE error is non-negligible but still acceptable (around 0.18) while the average scaled LDE is disappointingly close to 0.9. However, this heuristic allows us to solve instances which the MP formulations of Sect. 3 cannot solve. 4.3 Solving the DGP in the `1 norm At first, we tested a Variable Neighbourhood Search (VNS) [12] type algorithm with neighbourhoods defined by centers given by infeasible realization x0 , and radii given by the maximum number of coordinates of x0 that are allowed to change during a BB search. This constraint is enforced by a Local Branching (LB) [9] mechanism added to Eq. (3). The centers are computed using the prob- abilistic constructive heuristic proposed in [15], which extends the concept of Fréchet embeddings to sets. Unfortunately, this idea yielded very poor results, both quality-wise and in terms of CPU time. We only report it to prevent other researchers from pursuing this same direction. The second idea we tested is based on alternating solutions between two continuous reformulations of Eq. (2), denoted A (Eq. (2) with | · | replaced by p (·)2 + ε) and B ((3) with nonlinear complementarity constraints): we ran- domly sample an infeasible realization as a starting point for A, then use A’s solution as a starting point to B. We repeat this loop, updating the “best solu- tion found so far”, until a CPU time-based termination condition (set to 600s) becomes true. We used SNOPT [11] to solve A and IPOPT [3] to solve B. The results, obtained on 85 instances with n ∈ {10, 15, 20, . . . , 90} of any density s ∈ {0.1, 0.2, 0.3, 0.5, 0.8}, count 13 failures (generally instances with high den- sity where SNOPT had convergence issues) and 8 feasible realizations (scaled MDE/LDE scores < 10−4 ). If we consider only the results for the instances ended with no failures, overall, this heuristic displayed average scaled MDE and LDE of 0.13 and 1.16. References 1. N. Beeker, S. Gaubert, C. Glusa, and L. Liberti. Is the distance geometry problem in NP? In A. Mucherino, C. Lavor, L. Liberti, and N. Maculan, editors, Distance Geometry: Theory, Methods, and Applications. Springer, New York, 2013. 2. W.-Y. Chiu and B.-S. Chen. Mobile positioning problem in Manhattan-like urban areas: Uniqueness of solution, optimal deployment of BSs, and fuzzy implementa- tion. IEEE Transactions on Signal Processing, 57(12):4918–4929, 2009. 3. COIN-OR. Introduction to IPOPT: A tutorial for downloading, installing, and using IPOPT, 2006. 4. Gordon M. Crippen. An alternative approach to distance geometry using l dis- tances. Discrete Applied Mathematics, 197:20 – 26, 2015. Distance Geometry and Applications. 5. C. D’Ambrosio, G. Nannicini, and G. Sartor. MILP models for the selection of a small set of well-distributed points. working paper. 6. C. D’Ambrosio, K. Vu, C. Lavor, L. Liberti, and N. Maculan. New error measures and methods for realizing protein graphs from distance data. Technical Report submit-1604830, arXiv, 2016. 7. I. Dokmanić, R. Parhizkar, J. Ranieri, and M. Vetterli. Euclidean distance matrices: Essential theory, algorithms and applications. IEEE Signal Processing Magazine, 1053-5888:12–30, Nov. 2015. 8. P. Erdős and A. Rényi. On random graphs i. Publicationes Mathematicae (Debre- cen), 6:290–297, 1959. 9. M. Fischetti and A. Lodi. Local branching. Mathematical Programming, 98:23–37, 2005. 10. M. Fréchet. Sur quelques points du calcul fonctionnel. Rendiconti del Circolo Matematico di Palermo, 22:1–74, 1906. 11. P.E. Gill. User’s guide for SNOPT version 7.2. Systems Optimization Laboratory, Stanford University, California, 2006. 12. P. Hansen and N. Mladenović. Variable neighbourhood search: Principles and applications. European Journal of Operations Research, 130:449–467, 2001. 13. IBM. ILOG CPLEX 12.6 User’s Manual. IBM, 2014. 14. L. Liberti, C. Lavor, N. Maculan, and A. Mucherino. Euclidean distance geometry and applications. SIAM Review, 56(1):3–69, 2014. 15. J. Matousek. On the distortion required for embedding finite metric spaces into normed spaces. Israel Journal of Mathematics, 93:333–344, 1996. 16. J. Matoušek. Lecture notes on metric embeddings. Technical report, ETH Zürich, 2013. 17. K. Mehlhorn and P. Sanders. Algorithms and Data Structures. Springer, Berlin, 2008. 18. J. Saxe. Embeddability of weighted graphs in k-space is strongly NP-hard. Pro- ceedings of 17th Allerton Conference in Communications, Control and Computing, pages 480–489, 1979. 19. K. Vu, C. D’Ambrosio, Y. Hamadi, and L. Liberti. Surrogate-based methods for black-box optimization. International Transactions in Operational Research, to appear. 20. K. Wüthrich. Protein structure determination in solution by nuclear magnetic resonance spectroscopy. Science, 243:45–50, 1989.