40

I have a set of points (with unknown coordinates) and the distance matrix. I need to find the coordinates of these points in order to plot them and show the solution of my algorithm.

I can set one of these points in the coordinate (0,0) to simplify, and find the others. Can anyone tell me if it's possible to find the coordinates of the other points, and if yes, how?

Thanks in advance!

EDIT Forgot to say that I need the coordinates on x-y only

amWhy
  • 209,954
BrunoB
  • 403
  • 5
    You still have a rotational degree of freedom allowing you to assume that the second point is on the positive $x$-axis. Then there is still mirror symmetry so that you can assume that the third point (unless on the $x$-axis as well) can have positive $y$-coordinate. From that point on the system is rigid. Use the law of cosines heavily to get the angles. Add one point at a time. Work with triangles otherwise except use a fourth point to check whether the new angle is cw or ccw from an earlier reference direction. The system is heavily overdetermined, so you can use extra inputs for checks. – Jyrki Lahtonen Jun 09 '12 at 17:40
  • @JyrkiLahtonen this is a fantastic comment--thank you! – duhaime Apr 07 '18 at 02:13

2 Answers2

54

Doing this with angles, as Jyrki suggested, is cumbersome and difficult to generalize to different dimensions. Here is an answer that's essentially a generalization of WimC's, which also fixes an error in his answer. In the end, I show why this works, since the proof is simple and nice.

The algorithm

Given a distance matrix $D_{ij}$, define $$M_{ij} = \frac {D^2_{1j}+D^2_{i1}-D^2_{ij}} 2 \,.$$ One thing that is good to know in case the dimensionality of the data that generated the distance matrix is not known is that the smallest (Euclidean) dimension in which the points can be embedded is given by the rank $k$ of the matrix $M$. No embedding is possible if $M$ is not positive semi-definite.

The coordinates of the points can now be obtained by eigenvalue decomposition: if we write $M = USU^T$, then the matrix $X = U \sqrt S$ (you can take the square root element by element) gives the positions of the points (each row corresponding to one point). Note that, if the data points can be embedded in $k$-dimensional space, only $k$ columns of $X$ will be non-zero (corresponding to $k$ non-zero eigenvalues of $M$).

Why does this work?

If $D$ comes from distances between points, then there are $\mathbf x_i \in \mathbb R^m$ such that $$D_{ij}^2 = (\mathbf x_i - \mathbf x_j)^2 = \mathbf x_i^2 + \mathbf x_j^2 - 2\mathbf x_i \cdot \mathbf x_j \,.$$ Then the matrix $M$ defined above takes on a particularly nice form: $$M_{ij} = (\mathbf x_i - \mathbf x_1) \cdot (\mathbf x_j - \mathbf x_1) \equiv \sum_{a=1}^m \tilde x_{ia} \tilde x_{ja}\,,$$ where the elements $\tilde x_{ia} = x_{ia} - x_{1a}$ can be assembled into an $n \times m$ matrix $\tilde X$. In matrix form, $$M = \tilde X \tilde X^T \,.$$ Such a matrix is called a Gram matrix. Since the original vectors were given in $m$ dimensions, the rank of $M$ is at most $m$ (assuming $m \le n$).

The points we get by the eigenvalue decomposition described above need not exactly match the points that were put into the calculation of the distance matrix. However, they can be obtained from them by a rotation and a translation. This can be proved for example by doing a singular value decomposition of $\tilde X$, and showing that if $\tilde X \tilde X^T = X X^T$ (where $X$ can be obtained from the eigenvalue decomposition, as above, $X = U\sqrt S$), then $X$ must be the same as $\tilde X$ up to an orthogonal transformation.

Legendre17
  • 714
  • 2
    I hate to necro such an old thread, but can't find any other way to communicate. Legendre17, I'll be publishing a paper quite soon making use of this algorithm, and I'd like to give you due credit. You can find my contact info in my profile. – DaveTheScientist Nov 20 '15 at 22:02
  • 1
    Hi, @DaveTheScientist. This method is in no way invented by me. I don't remember exactly where I got it from when I wrote the answer, but some Googling led to this paper https://www.researchgate.net/publication/24061253_Multidimensional_Scaling_I._Theory_and_Method. I think they cite an even older paper for the method. – Legendre17 Nov 22 '15 at 22:50
  • 1
    Is this algorithm numerically stable? For example, suppose there is a small error bar associated with the distances $d_{ij}$. If I know that the points are embedded in 3-dimensions, when I compute $X$ with this method, probably more than 3 columns of $X$ will be distinct from zero, due to the error. Will these extra columns be very small, if the distance error bars are small? – a06e Jul 19 '16 at 13:20
  • 1
    Found this: Crippen, G. . (1978). Note rapid calculation of coordinates from distance matrices. Journal of Computational Physics, 26(3), 449–452. http://doi.org/10.1016/0021-9991(78)90081-5 – a06e Jul 19 '16 at 14:20
  • 2
    I haven't thought about this a whole lot, but my guess is that adding noise to the distances would essentially just push the vanishing eigenvalues of the $M$ matrix slightly away from zero. If you put a threshold on them, you should be able to get rid of the noise. Note that if you add noise to the distances, you will likely spoil the positive-definiteness of the $M$ matrix, so the threshold would have to be placed on the absolute value of the eigenvalues. This also means that you need to do the thresholding before taking the square root $\sqrt S$ and calculating $X$, to avoid complex entries. – Legendre17 Jul 28 '16 at 21:17
  • 1
    @Legendre17 Good answer, thanks! Also, you remember correct, they cite another paper - Young and Householder, Discussion of a set of points in terms of their mutual distances. Paychometrika, 1938, 3, 19-22. – Itay Feb 27 '17 at 23:25
  • So I'm imagining a coordinate system where you have multiple floors and obviously the closest point in the upper floor to the downstairs exit would be the elevator no matter how close other points are in terms of euclidean space. Is it possible to take this into account? – pwned Jul 29 '18 at 18:33
  • 1
    @pwned This calculation refers specifically to Euclidean distances. I don't know how to generalize it to cases in which you have floors and elevators. – Legendre17 Jul 30 '18 at 19:46
  • @Legendre17 I think it goes down to graph theory, thanks. – pwned Jul 31 '18 at 05:29
  • When you say then X is up to $\hat{X}$ an orthogonal transformation. Is there a way to calculate this orthogonal transformation matrix? – FantasticAI Nov 20 '22 at 03:39
  • 1
    No, the precise orientation and origin of the coordinate system is lost when you calculate distances. Since all you have access to is the matrix of distances, there is no way to recover the origin and rotation information. – Legendre17 Nov 21 '22 at 14:33
9

I assume this is in $\mathbb{R}^2$. You can compute possible points $v_j \in \mathbb{R}^2$ as follows. First compute the Gram matrix $M$ from the distance matrix $D$:

$$M_{i,j} = \frac{(D_{1,i})^2 + (D_{1, j})^2 - (D_{i, j})^2}{2}$$

for $i, j \in \{1,\dotsc, n\}$. (Note that the first row and column of $M$ consist of zeroes.) This is a positive semi-definite matrix of rank at most two. I'll assume the generic case where the rank is two (not all points on a single line).

Find an orthonormal basis $\{b_1, b_2\}$ for the column space of $M$ (e.g. by applying Gram-Schmidt). Let $m_j$ denote the $j$-th column of $M$. Then take $v_j = (\langle b_1, m_j \rangle, \langle b_2, m_j \rangle) \in \mathbb{R}^2$ where $\langle \cdot, \cdot \rangle$ is the euclidean inner product on $\mathbb{R}^n$. The points $v_j$ will satisfy your distance matrix.

WimC
  • 32,192
  • 2
  • 48
  • 88