Reconstructing a point set from a Euclidean Distance Matrix

Recovering point locations from Euclidean Distance Matrices has practical applications in ultrasound tomography and other fields that must determine the position of recording devices from acoustic signals. This post explains the solution to a concrete problem, drawing on the explanations in Linear Algebra and Learning from Data and a paper on Euclidean distance geometry.

I want to answer the question: How to reconstruct the locations of the original vectors given the Euclidean distance matrix, ?

We assume really is a Euclidean distance matrix, so that a position matrix exists, where has rows for points in -dimensional space. The triangle inequality is necessary for that but not sufficient; the test that settles it is whether the Gram matrix built below comes out positive semidefinite, which is also what tells you the dimension . After identifying the point set, you can align it with known locations using Procrustes analysis, also known as the Orthogonal Procrustes problem.

Theory

This procedure, known as classical multidimensional scaling (MDS), comes from this paper. The approach extends to real-world problems like noisy data.

Consider a collection of points in a -dimensional Euclidean space. The squared distance between points and is and is calculated using:

Expanding this norm yields:

The matrix equation for the distance matrix is calculated using:

is the column vector of all ones and is a column vector of the diagonal entries of .

The operator is defined which is equivalent to that operates directly on the Gram matrix . This now gives the equation below:

Let the first point be the origin, then the first column of contains the squared norms of the point vectors,

It is now possible to construct the term and its transpose in the equation to calculate .

The Gram matrix can now be found from like so:

The final stage is to identify the point set using Eigenvalue Decomposition (EVD), for example:

Remember that

Example problem

This problem is Q5 from Problem Set from the Linear Algebra and Learning from Data textbook. Given a Euclidean distance matrix, , find the locations of the points:

Let’s start with some useful imports:

import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import distance_matrix
from scipy.spatial.distance import cdist

And now let’s define the matrix :

D = np.array([[0, 9, 25], [9, 0, 16], [25, 16, 0]])

Calculate the Gram matrix. Any point can serve as the origin; the code below uses the second, because it makes diagonal and the arithmetic easy to follow. On why the result is positive semidefinite whenever is a genuine distance matrix, see J. C. Gower, “Euclidean Distance Geometry,” Math. Sci., vol. 7, pp. 1–14, 1982.

G = -0.5 * (D - np.outer(np.ones(3), D[1, :]) - np.outer(D[:, 1], np.ones(3)))
G

array([[9., -0., -0.],
       [-0., -0., -0.],
       [-0., -0., 16.]])

Now use np.linalg.svd to solve .

Q, Lambda, _ = np.linalg.svd(G)
print(Q)
print(Lambda)

[[0. 1. 0.]
 [0. 0. 1.]
 [1. 0. 0.]]
[16. 9. -0.]

Return the original point set using . The multiplication order matters here: np.sqrt(Lambda) * Q.T scales the columns of and quietly hands back a point set whose distance matrix is a permutation of , which passes a casual glance because the entries are all still there.

np.diag(np.sqrt(Lambda)) @ Q.T

array([[0., 0., 4.],
       [3., 0., 0.],
       [0., 0., 0.]])

The columns are the points. The third row is zero because the third eigenvalue is, so the three points lie in a plane: , and . Their squared distances are 9, 25 and 16, which is .